Small Box FFT学习笔记(四):FFT

本篇主要介绍程序中底层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指的是离散傅里叶变换而非密度泛函理论。

发表回复

您的电子邮箱地址不会被公开。 必填项已用 * 标注