本篇主要介绍程序中底层FFT调用,主要参考的是FFTP.f90文件(串行版本是FFTS.f90)中的Fourier_new模块。
程序中实际是调用FFTW库,由于在Fortran 03以后的版本中可以使用C语言API,这也是FFTW鼓励的做法(参考FFTW官方文档),为此需要在开头加入一行:
use, intrinsic :: iso_c_binding
程序中将FFT相关的信息封装成了一个FFT_CONFIG数据类型,其定义如下
type FFT_CONFIG
character(len=4) id ! 标记,例如DGBX、DSBX
integer fftAlgo ! FFT计算模式
integer(c_intptr_t) totalDimX, totalDimX, totalDimZ ! 三个方向的阶数
integer(c_intptr_t) localDimZ ! z方向的局域阶数
integer(c_intptr_t) localDimZOff
integer iCountFFT ! FFT计数器
integer offset
real*8 normConst ! 如果定义了FFT_NORM_EXACT,则类型为integer*8
type(c_ptr) fftw3GlobFWD,fftw3GlobBWD
type(c_ptr) cdata
real(c_double), pointer, dimension( *, *, * ) in ! 指向实空间数组的指针
complex(c_double_complex), pointer, dimension( *, *, * ) out ! 指向倒空间数组的指针
end type FFT_CONFIG
可以看到其中使用了诸如c_intptr_t和c_double这样的数据类型。基于这一数据类型,定义了一个统一的FFT_NEW接口用于不同数据类型的正向和逆向傅里叶变换:
interface FFT_NEW
module procedure ForwardFFT_4D
module procedure BackFFT_4D
module procedure ForwardFFT_3D
module procedure BackFFT_3D
end interface
这里通过输入参数的类型来选择使用的具体子程序,例如当第二个参数是实三维数组而第三个参数是复三维数组时(第一个参数是FFT_CONFIG类型)则调用ForwardFFT_3D。由于使用的是FFTW库,在实际计算前需要建立一个FFTW_PLAN,模块里定义了一个PlanFFT_NEW子程序用于初始化并创建 FFTW_PLAN ,其定义为
subroutine PlanFFT_NEW(
integer MPI_COMM,
type(FFT_CONFIG) config,
integer mode,
integer dimX,dimY,dimZ
)
这里mode是选择创建FFTW_PLAN的方式,这一子程序中关键的调用为:
totalLocalSize = fftw_mpi_local_size_3d( config%totalDimZ, config%totalDimY, config%totalDimX/2+1, MPI_COMM, config%localDimZ, config%localDimZOff )
config%cdata = fftw_alloc_complex( totalLocalSize )
call c_f_pointer( config%cdata, config%in, [2*(config%totalDimX/2+1), config%totalDimY, config%localDimZ] )
call c_f_pointer( config%cdata, config%out, [config%totalDimX/2+1, config%totalDimY, config%localDimZ] )
config%fftw3GlobFWD = fftw_mpi_plan_dft_r2c_3d( config%totalDimZ, config%totalDimY, config%totalDimX, config%in, config%out, MPI_COMM, planDepth)
config%fftw3GlobBWD = fftw_mpi_plan_dft_c2r_3d( config%totalDimZ, config%totalDimY, config%totalDimX, config%out, config%in, MPI_COMM, planDepth )
这里的planDepth是由输入参数mode决定。最后,我们以ForwardFFT_3D为例来说明如何调用FFTW库中的进行计算,这一子程序的定义为:
subroutine ForwardFFT_3D(
type(FFT_CONFIG) config,
real*8, dimension( *, *, * ) array,
complex*8, dimension( *, *, * ) transform
)
这一子程序首先将输入的数组array传入config%in,并进行如下调用:
call fftw_mpi_execute_dft_r2c( config%fftw3GlobFWD, config%in, config%out )
再将config%out传到输出数组transform,就完成了array到transform的(正向)FFT计算。注意这里所有中的DFT指的是离散傅里叶变换而非密度泛函理论。