FFTW 3 分段错误

Question

在将 FFTW 与 fortran 一起使用时，我遇到了一个奇怪的分段错误，我无法弄清楚原因。在线搜索也没有运气。这是我的代码：

 integer, parameter :: Nx=128, Ny=Nx, Nz=Nx
 integer, parameter :: NORM=Nx*Ny*Nz,Ntop=MAX(Nx,Ny,Nz)/2, LX1=Nx/2+1
 double precision, parameter :: Etot=0.5
 double complex, dimension(LX1,Ny,Nz) :: fbx, fby, fbz
 double precision, dimension(Nx,Ny,Nz) :: bx, by, bz
 real
 ! ...
 ! Assign values to fbx, fby, fbz 
 ! ...
 call dfftw_plan_dft_c2r_3d(plan,Nx,Ny,Nz,fbx,bx,FFTW_ESTIMATE)
 call dfftw_execute_dft_c2r(plan, fbx, bx)
 call dfftw_execute_dft_c2r(plan, fby, by)
 call dfftw_execute_dft_c2r(plan, fbz, bz)
 call dfftw_destroy_plan(plan)
 ! et = sum(abs(fbx)**2+abs(fby)**2+abs(fbz)**2)
 ! bx=bx*sqrt(Etot/et); by=by*sqrt(Etot/et); bz=bz*sqrt(Etot/et)

 ! Write bx, by, bz to files

代码以这种方式运行良好。当我取消注释计算 et 和重新缩放 bx 行时，问题就出现了。当我取消注释这些行时，我得到分段错误，说“无效的内存引用”。第一次执行dfftw_execute_dft_c2r()时出现seg错误

我已经在两台不同的机器上尝试了 fftw 3.2 和 fftw 3.4 以及 ifort 以及 gfortran（这些机器的所有可能组合），但没有任何成功。

快把我逼疯了！！我需要帮助！为什么会这样？

谢谢！

Answer 1

您正在使用应用于不同阵列的相同计划。 FFTW 仅在它们在内存中对齐的情况下才有效 - 参见例如final point here，以及文档中的相关部分，说明为什么这很难保证。

我在文档或其他地方看到的两个建议是

1. 创建临时数组，为这些临时数组构建计划，并在调用 execute_dft 之前/之后将数据复制到临时数组中或从临时数组中复制出来。我过去曾用它来避免类似的错误，但这并不是最好的性能。但是，这是一个简单的更改，可以轻松检查这是否是问题的原因。
2. 创建多个计划。

此外，您似乎正在使用legacy Fortran interface。如果可能的话，最好使用"Modern" interface，但我过去的经验是一些编译器不能很好地使用它。