利用MKL库FFT求导

wt_tree

在利用mkl提供的fft模块计算周期性函数的导数时进行测试

1.  利用傅里叶变换的微分性质：\mathscr{F}(f‘(x))=i\omega \mathscr{F}(f(x))
[latex]$$f'(x)=\mathscr{F}^{-1}\{\mathscr{F}[f'(x)]\}=\mathscr{F}^{−1}\{i\omega F[f(x)]\}$$[/latex]

2. 调用MKL库中的FFT模块，具体调用过程可以参考oneAPI官方手册，代码如下


        

[Fortran] syntaxhighlighter_viewsource syntaxhighlighter_copycode

module fft

        use MKL_DFTI

        implicit none

        real(8), parameter :: pi = dacos(-1.0D0)

        

        contains

            subroutine fft1d(n, x_in, x_out)

                implicit none

                integer, intent(in) :: n

                complex(8), intent(in) :: x_in(n)

                complex(8), intent(out) :: x_out(n)

                type(DFTI_DESCRIPTOR), pointer :: My_Desc_Handle

                Integer :: Status

                Status = DftiCreateDescriptor( My_Desc_Handle, DFTI_DOUBLE, DFTI_COMPLEX, 1, n )

                Status = DftiSetValue( My_Desc_Handle, DFTI_PLACEMENT, DFTI_NOT_INPLACE)

                Status = DftiCommitDescriptor( My_Desc_Handle )

                Status = DftiComputeForward( My_Desc_Handle, x_in, x_out )

                Status = DftiFreeDescriptor(My_Desc_Handle)

        

            end subroutine

        

            subroutine ifft1d(n, x_in, x_out)

                implicit none

                integer, intent(in) :: n

                complex(8), intent(in) :: x_in(n)

                complex(8), intent(out) :: x_out(n)

          

                type(DFTI_DESCRIPTOR), pointer :: My_Desc_Handle

                Integer :: Status

                Status = DftiCreateDescriptor( My_Desc_Handle, DFTI_DOUBLE,DFTI_COMPLEX, 1, n )

                Status = DftiSetValue( My_Desc_Handle, DFTI_PLACEMENT, DFTI_NOT_INPLACE)

                Status = DftiSetValue( My_Desc_Handle, DFTI_BACKWARD_SCALE, 1.0D0/N)

                Status = DftiCommitDescriptor( My_Desc_Handle )

                Status = DftiComputeBackward( My_Desc_Handle, x_in, x_out )

                Status = DftiFreeDescriptor(My_Desc_Handle)

                

            end subroutine

            

            subroutine fft_1st_deriv_1d(N, x_in, deriv_x, L)

                        implicit none

                        integer, intent(in) :: N

                        complex(8), intent(in) :: x_in(N)

                        complex(8), intent(out) :: deriv_x(N)

                        real(8), intent(in) :: L

                complex(8) :: fft_x(N), dfft_x(N)

                integer :: k

                real(8) :: kx(N)

                call fft1d(N, x_in, fft_x)

                do k = 1, N

                    if(mod(N,2) .eq. 0) then

                        write(*,*) N

                        if( k .eq. N/2+1) then

                            kx(k) = 2*pi/L*0

                            dfft_x(k) = cmplx(0.0D0, 1.0D0,kind=8)*kx(k)*fft_x(k)

                        else if(k .le. N/2) then

                            kx(k) = 2*pi/L*(k-1)

                            dfft_x(k) = cmplx(0.0D0, 1.0D0,kind=8)*kx(k)*fft_x(k)

                        else

                            kx(k) = 2*pi/L*(k-1-N)

                            dfft_x(k) = cmplx(0.0D0, 1.0D0,kind=8)*kx(k)*fft_x(k)

                        end if

                    else

                        if(k .le. N/2) then

                            dfft_x(k) = 2*pi*dcmplx(0.0D0, 1.0D0)/L*(k-1)*fft_x(k)

                        else

                            dfft_x(k) = 2*pi*dcmplx(0.0D0, 1.0D0)/L*(k-N-1)*fft_x(k)

                        end if

                    end if

                end do

        

                call ifft1d(N, dfft_x, deriv_x)

                open(unit = 16, file="dfft.txt")

        

            end subroutine

        end module fft

        

        program  test

            use fft

            implicit none

            integer, parameter :: n =128

            complex(8) :: x_out(n), z(n), deriv_x(n), x_in(n)

            ! real(8), parameter :: pi = 4.0D0*datan(1.0D0)

            real(8) :: l =2*pi

            integer(8) :: i

            real(8) :: dt

            write(*,*) pi

            do i = 1, n

                ! x_in(i) = exp(dcmplx(0.0D0,2.0D0*pi)*real(i)/n)

                x_in(i) = cmplx(cos(2*pi*(i-1)/(n-1)),kind = 8)

            end do

            call fft1d(n,x_in,x_out)

            call ifft1d(n,x_out,z)    

            call fft_1st_deriv_1d(n,x_in, deriv_x, l)

        end program  test

在用此办法求导过程中，改变输入函数x_in(i)时发现：

[Fortran] syntaxhighlighter_viewsource syntaxhighlighter_copycode

 x_in(i) = cmplx(sin(2*pi*(i-1)/(n-1)),kind = 8)

计算得到的导数两端会有振荡，如下图所示：

()

sinx的导数，红色为FFT计算得到的结果，黑色为cosx

但当输入函数改变为：

[Fortran] syntaxhighlighter_viewsource syntaxhighlighter_copycode

 x_in(i) = cmplx(sin(2*pi*(i-1)/n),kind = 8)

振荡消失，结果如下图所示：

()

输入函数为cosx没有遇到这一问题。结果均收敛。但是结果对比matlab的计算和python中numpy模块的fft精度上略有不足。所以这个是什么原因导致的？

li913

cmplx函数有一个bug，反常识，就是不管参数是单精度还是双进度，cmplx(a,b) 的结果都是单精度，cmplx(a,b,8) 才是双精度。

楚香饭

li913 发表于 2022-11-24 22:31
cmplx函数有一个bug，反常识，就是不管参数是单精度还是双进度，cmplx(a,b) 的结果都是单精度，cmplx(a,b,8 ...

还真是这样，妈呀，吓得我赶快把以前的project全局搜索了一遍，发现了好多处cmplx，赶快都改成了dcmplx

wt_tree

li913 发表于 2022-11-24 22:31
cmplx函数有一个bug，反常识，就是不管参数是单精度还是双进度，cmplx(a,b) 的结果都是单精度，cmplx(a,b,8 ...

貌似不是这个问题，我调用的时候都是用的cmplx(a,b,kind=8)，而且求导后精度问题也很麻烦。