c++ - 与 fftw3 相比,错误的 2D CuFFT 逆变换
问题描述
我正在尝试产生一些 FFT 数学,特别是它进行两个 2D 正向变换,将它们相乘,然后进行逆变换。在逆变换之前一切都很好。我已经通过 fftw3 做到了,但是在 CuFFT 中出现了问题。大多数值是相似的,但有些是错误的,这对未来的数学很重要。这段代码有什么问题?
std::vector<complex> conv2dCUDA(complex *ui_anomaly, double *ds2,
complex *u0, int anx, int any, double factor) {
cufftComplex *b1, *b2;
int size = 2 * anx * 2 * any;
int memsize = size * sizeof(cufftComplex);
b1 = (cufftComplex *)calloc(size, sizeof(cufftComplex));
b2 = (cufftComplex *)calloc(size, sizeof(cufftComplex));
// filling the matrixes
cufftHandle plan;
cufftComplex *ui, *g;
checkCudaErrors(cudaMalloc((void**)&ui, memsize));
checkCudaErrors(cudaMalloc((void**)&g, memsize));
checkCudaErrors(cufftPlan2d(&plan, 2 * anx, 2 * any, CUFFT_C2C));
checkCudaErrors(cudaMemcpy(ui, (cufftComplex *)&b1[0], memsize, cudaMemcpyHostToDevice));
checkCudaErrors(cudaMemcpy(g, (cufftComplex *)&b2[0], memsize, cudaMemcpyHostToDevice));
checkCudaErrors(cufftExecC2C(plan, ui, ui, CUFFT_FORWARD));
checkCudaErrors(cufftExecC2C(plan, g, g, CUFFT_FORWARD));
int blockSize = 16;
dim3 dimGrid(int(2 * any / blockSize) + 1, int(2 * anx / blockSize) + 1);
dim3 dimBlock(blockSize, blockSize);
ComplexMulAndScale<<<dimGrid, dimBlock>>>(ui, g, size, 1.0f);
getLastCudaError("Kernel execution went wrong");
checkCudaErrors(cudaMemcpy(b1, ui, memsize, cudaMemcpyDeviceToHost));
std::cout << "After mult Cuda" << std::endl;
for (auto i = 0; i < 2 * any; i++) {
for (auto j = 0; j < 2 * anx; j++) {
std::cout << b1[i * 2 * anx + j].x << " ";
}
std::cout << std::endl;
}
checkCudaErrors(cufftExecC2C(plan, ui, ui, CUFFT_INVERSE));
cuComplex *inversed;
inversed = (cuComplex*)malloc(memsize);
checkCudaErrors(cudaMemcpy(inversed, ui, memsize, cudaMemcpyDeviceToHost));
std::vector<complex> res(anx * any);
for (auto i = 0; i < any; i++) {
for (auto j = 0; j < anx; j++) {
res[i * anx + j] = complex(inversed[i * anx * 2 + j].x * factor, inversed[i * anx * 2 + j].y * factor);
}
}
std::cout << "CUDA" << std::endl;
for (auto i = 0; i < 2 * any; i++) {
for (auto j = 0; j < 2 * anx; j++) {
std::cout << inversed[i * 2 * anx + j].x << " ";
}
std::cout << std::endl;
}
checkCudaErrors(cudaFree(ui));
checkCudaErrors(cudaFree(g));
checkCudaErrors(cufftDestroy(plan));
free(b1);
free(b2);
free(inversed);
return res;
}
std::vector<complex> conv2d(complex *ui_anomaly, double *ds2, complex *u0, int anx, int any, double factor) {
std::vector<complex> b1(anx * 2 * 2 * any, complex(0., 0.)), b2(anx * 2 * 2 * any, complex(0., 0.));
// filling matrixes
// forward fft 1 in-place
fftw_plan p;
p = fftw_plan_dft_2d(2 * any, 2 * anx, (fftw_complex *) (&b1[0]), (fftw_complex *) (&b1[0]),
FFTW_FORWARD, FFTW_ESTIMATE);
fftw_execute(p);
fftw_destroy_plan(p);
// forward fft 2 in-place
p = fftw_plan_dft_2d(2 * any, 2 * anx, (fftw_complex *) (&b2[0]), (fftw_complex *) (&b2[0]),
FFTW_FORWARD, FFTW_ESTIMATE);
fftw_execute(p);
fftw_destroy_plan(p);
std::vector<complex> out(2 * anx * 2 * any, complex(0., 0.));
for (auto i = 0; i < 2 * any * 2 * anx; i++) {
out[i] = b1[i] * b2[i];
}
std::cout << "After mult fftw" << std::endl;
for (auto i = 0; i < 2 * any; i++) {
for (auto j = 0; j < 2 * anx; j++) {
std::cout << out[i * 2 * anx + j].real() << " ";
}
std::cout << std::endl;
}
// inverse fft in-place
p = fftw_plan_dft_2d(2 * (int) any, 2 * (int) anx, (fftw_complex *) (&out[0]), (fftw_complex *) (&out[0]),FFTW_BACKWARD, FFTW_ESTIMATE);
fftw_execute(p);
fftw_destroy_plan(p);
std::vector<complex> res(anx * any);
for (auto i = 0; i < any; i++) {
for (auto j = 0; j < anx; j++) {
res[i * anx + j] = out[i * anx * 2 + j] * factor;
}
}
std::cout << "FFTW" << std::endl;
for (auto i = 0; i < 2 * any; i++) {
for (auto j = 0; j < 2 * anx; j++) {
std::cout << out[i * 2 * anx + j].real() << " ";
}
std::cout << std::endl;
}
return res;
}
所以,这是我的代码。输出应该在两个函数中
After mult fftw
8.34304e-08 -5.99259e-07 -4.7876e-07 5.30254e-07 9.55877e-07 4.28985e-07
-1.56375e-07 1.19699e-07 2.39276e-07 -1.68662e-08 -7.56988e-08 -3.69897e-07
-2.66505e-07 -2.33361e-07 -5.21763e-07 -5.29126e-07 1.8915e-07 1.68158e-07
-9.01859e-07 -2.37453e-07 -3.50661e-08 -4.11154e-07 4.14802e-07 -7.9879e-07
2.09404e-07 6.52034e-08 1.8915e-07 4.97805e-07 3.32612e-07 -2.33361e-07
-1.95738e-07 -3.69897e-07 -1.63577e-07 1.07737e-07 2.39276e-07 2.50198e-07
FFTW
-1.57349e-06 -7.5964e-06 -1.57349e-06 1.68876e-06 5.82335e-22 1.68876e-06
2.37158e-06 6.35275e-22 2.37158e-06 -1.18579e-06 1.05879e-22 -1.18579e-06
-1.57349e-06 -7.5964e-06 -1.57349e-06 1.68876e-06 1.97573e-22 1.68876e-06
3.14928e-06 2.37158e-06 3.14928e-06 -4.94164e-07 5.82335e-22 -4.94164e-07
2.11758e-22 -8.47033e-22 -1.05879e-22 5.29396e-22 1.41851e-23 1.05879e-22
3.14928e-06 2.37158e-06 3.14928e-06 -4.94164e-07 1.05879e-22 -4.94164e-07
After mult Cuda
8.34303e-08 -5.99259e-07 -4.78761e-07 5.30254e-07 9.55877e-07 4.28985e-07
-1.56375e-07 1.19699e-07 2.39276e-07 -1.68662e-08 -7.56988e-08 -3.69897e-07
-2.66505e-07 -2.33361e-07 -5.21763e-07 -5.29126e-07 1.8915e-07 1.68158e-07
-9.01859e-07 -2.37453e-07 -3.50661e-08 -4.11154e-07 4.14802e-07 -7.9879e-07
2.09404e-07 6.52034e-08 1.8915e-07 4.97805e-07 3.32612e-07 -2.33361e-07
-1.95738e-07 -3.69897e-07 -1.63577e-07 1.07737e-07 2.39276e-07 2.50198e-07
CUDA
-1.57349e-06 -7.5964e-06 -1.57349e-06 1.68876e-06 3.33981e-13 1.68876e-06
2.37158e-06 2.84217e-13 2.37158e-06 -1.18579e-06 1.10294e-13 -1.18579e-06
-1.57349e-06 -7.5964e-06 -1.57349e-06 1.68876e-06 -9.03043e-14 1.68876e-06
3.14928e-06 2.37158e-06 3.14928e-06 -4.94164e-07 4.62975e-13 -4.94164e-07
-3.2685e-13 -1.03562e-13 -3.59548e-13 -2.13163e-13 4.27658e-15 -2.43358e-14
3.14928e-06 2.37158e-06 3.14928e-06 -4.94164e-07 3.49288e-13 -4.94164e-07
可以看出,正向 fft 和乘法都是正确的,但是在 cuda smth 的逆 fft 的情况下出错了。
PS抱歉代码风格不佳
解决方案
由于使用了 FFTW,卷积后的信号具有很多约 1e-6 的数字和少数约 1e-22 的数字。它很可能应该是零,但不是因为这些零是使用双精度浮点数计算的。双精度数字大约精确到 15 位,因此可能会出现接近 1e(-6-15)=1e-21 的错误。
当使用 cufft 时,这些应该为零的数字大约是 1e-13,就好像使用单精度浮点数执行了计算一样。
就是这样:类型cuComplex
和cufftComplex
是单精度复数,而fftw_complex
是双精度复数。虽然complex
可能默认为双精度,但可以明确指定为double complex
.
要获得 1e-22 附近的数字,请尝试类型cufftDoubleComplex
和cuDoubleComplex
。一次执行多次乘法的blockSize
引入可能需要减少到 8。然而,虽然很可能可以获得 1e-22 级的数字,但这些数字也很可能与 FFTW 的数字不同。事实上,由于算法可能不同,因此可能执行了不同的操作,并且精度使得结果中大约 1e-22 的任何值与 0 没有显着差异。
然而,更改双精度数可能会增加计算时间并明显增加内存占用。如果卷积结果的六位数精度对您的应用程序来说足够好,那么坚持单精度复 DFT 可能是正确的方法。
推荐阅读
- node.js - 将图像上传到 Azure Blob 存储
- ubuntu - ubuntu apt-get 未满足的 libsmbios2v5 和 libsmbios2 依赖项
- json - 如何使用 jq 提取此 JSON 片段中 href 的值?
- ruby-on-rails - Rails 活动模型序列化程序如何处理嵌套资源?
- mysql - 点燃和参考表
- html - 滤镜亮度过渡 CSS3
- user-interface - Inkscape的主题定制
- javascript - 使用 Javascript 数组创建 HTML 表
- sonarqube - Sonarqube 扫描仪错误 DirectoryNotEmptyException
- php - OTP 硬件符合 OATH TOTP 和 php 库