使用pyculib进行3D FFT的结果是错误的

Question

我使用pyculib在Anaconda 3.5中对矩阵执行3D FFT。我刚跟着the example code发布在网站上。但我发现了一些有趣的东西并且不明白为什么。

仅当使用numpy.arange创建矩阵时，才能使用pyculib对矩阵执行3D FFT。

这是代码：

from pyculib.fft.binding import Plan, CUFFT_C2C
import numpy as np
from numba import cuda
data = np.random.rand(26, 256, 256).astype(np.complex64)
orig = data.copy()
d_data = cuda.to_device(data)
fftplan = Plan.three(CUFFT_C2C, *data.shape)
fftplan.forward(d_data, d_data)
fftplan.inverse(d_data, d_data)
d_data.copy_to_host(data)
result = data / n
np.allclose(orig, result.real)

最后，结果证明是假的。而orig和结果之间的差异不是一个小数字，不可忽略不计。我尝试了一些其他数据集（不是随机数），并得到一些错误的结果。

另外，我测试没有逆FFT：

from pyculib.fft.binding import Plan, CUFFT_C2C
import numpy as np
from numba import cuda
from scipy.fftpack import fftn,ifftn
data = np.random.rand(26,256,256).astype(np.complex64)
orig = data.copy()
orig_fft = fftn(orig)
d_data = cuda.to_device(data)
fftplan = Plan.three(CUFFT_C2C, *data.shape)
fftplan.forward(d_data, d_data)
d_data.copy_to_host(data)
np.allclose(orig_fft, data)

结果也是错误的。

网站上的测试代码，他们使用numpy.arange来创建矩阵。我试过了：

n = 26*256*256
data = np.arange(n, dtype=np.complex64).reshape(26,256,256)

并且该矩阵的FFT结果是正确的。

任何人都可以帮忙指出原因吗？

Answer 1

我不使用CUDA，但我认为你的问题本质上是数字的。不同之处在于您使用的两个数据集。 random.rand的动态范围为0-1，范围为0-26 * 256 * 256。 FFT尝试按值/点数范围的顺序解析空间频率分量。对于范围，这变为单位，并且FFT在数值上是准确的。对于rand，这是1/26 * 256 * 256~5.8e-7。

在不使用CUDA的情况下在numpy阵列上运行FFT / IFFT会显示出类似的差异。