FFTW3计划在同名同尺寸数组的循环迭代中是否需变更？

问题背景 #

使用FEAP8.4（Fortran）结合FFTW3库，对同名同尺寸数组执行傅里叶变换，主程序通过时间循环调用如下子routine：

subroutine XYZ()
 implicit none

 integer*8 plan_r2c_p
 integer*8 N
 PARAMETER(N=8)
 real*8 psi(N), psiHat(N/2+1)

 call dfftw_plan_dft_r2c_1d(plan_r2c_p,N,psi,psiHat,FFTW_ESTIMATE)
 write(*,*) 'plan_r2c_p=',plan_r2c_p
end subroutine XYZ

运行发现每次循环迭代中plan_r2c_p的数值会变化（量级相同），存在以下疑问：

1. 同尺寸数组在同一机器上的FFTW3计划应保持一致，该理解是否正确？
2. 若理解错误，计划数值变化的原因是什么？
3. 是否可以通过复用计划而非每次创建来节省计算资源？

解答 #

1. 初始理解错误：FFTW的plan变量本质是指向内部计算策略结构体的内存句柄/指针，其数值（内存地址）是否一致不代表计划的计算逻辑是否相同。只要创建计划时的参数（数组尺寸、数据类型、变换方向、优化flags等）完全一致，同一机器上生成的计划底层计算策略是一致的，和plan变量的数值无关。

2. 计划数值变化的原因：每次调用dfftw_plan_dft_r2c_1d时，FFTW都会在堆内存中分配新的计划结构体，返回的plan_r2c_p是这个新结构体的内存地址，因此每次循环的数值自然不同。但这些不同地址的计划，只要参数一致，执行变换的逻辑和性能是完全相同的。

3. 完全可以复用计划，且这是FFTW性能优化的核心手段：创建计划的过程会涉及算法搜索、性能评估（即使是FFTW_ESTIMATE也会做基础的策略选择），重复创建会浪费大量计算资源，复用计划能显著提升循环内的执行效率。

复用计划的实现方法 #

需要将计划变量提升为全局变量（通过模块），或在主程序中提前创建并作为参数传递给子routine，避免在循环内重复创建。示例如下：

方法1：使用模块定义全局计划

! 定义存储计划和数组的模块
module fftw_setup
    implicit none
    integer*8, parameter :: N=8
    integer*8 plan_r2c_p
    real*8 psi(N), psiHat(N/2+1)
end module fftw_setup

! 修改后的子routine，直接复用已创建的计划
subroutine XYZ()
    use fftw_setup
    implicit none
    ! 直接执行变换，无需重复创建计划
    call dfftw_execute_dft_r2c(plan_r2c_p, psi, psiHat)
end subroutine XYZ

! 主程序：仅初始化一次计划，循环调用子routine
program main
    use fftw_setup
    implicit none
    integer i

    ! 仅在程序启动时创建一次计划
    call dfftw_plan_dft_r2c_1d(plan_r2c_p, N, psi, psiHat, FFTW_ESTIMATE)

    ! 时间循环
    do i = 1, 100
        ! 此处更新psi数组的数值
        call XYZ()
        ! 处理变换后的psiHat
    end do

    ! 程序结束前销毁计划，避免内存泄漏
    call dfftw_destroy_plan(plan_r2c_p)
end program main

方法2：主程序创建计划并传递给子routine

如果数组是主程序的局部变量，可将计划作为参数传入子routine：

subroutine XYZ(plan_r2c_p, psi, psiHat, N)
    implicit none
    integer*8, intent(in) :: plan_r2c_p, N
    real*8, intent(in) :: psi(N)
    real*8, intent(out) :: psiHat(N/2+1)

    call dfftw_execute_dft_r2c(plan_r2c_p, psi, psiHat)
end subroutine XYZ

program main
    implicit none
    integer*8, parameter :: N=8
    integer*8 plan_r2c_p
    real*8 psi(N), psiHat(N/2+1)
    integer i

    call dfftw_plan_dft_r2c_1d(plan_r2c_p, N, psi, psiHat, FFTW_ESTIMATE)

    do i = 1, 100
        ! 更新psi
        call XYZ(plan_r2c_p, psi, psiHat, N)
        ! 处理psiHat
    end do

    call dfftw_destroy_plan(plan_r2c_p)
end program main

内容的提问来源于stack exchange，提问作者R2197