I'm working over transforming some code from using FFTW library to CUFFT (CPU computing to GPU computing). I need to transform a matrix of forces, make some math on it and transform it back. Operation in FFTW looks like it:
fftw_real u0[DIM * 2*(DIM/2+1)], v0[DIM * 2*(DIM/2+1)];
static rfftwnd_plan plan_rc, plan_cr;
void init_FFT(int n) {
plan_rc = rfftw2d_create_plan(n, n, FFTW_REAL_TO_COMPLEX, FFTW_IN_PLACE);
plan_cr = rfftw2d_create_plan(n, n, FFTW_COMPLEX_TO_REAL, FFTW_IN_PLACE);
}
#define FFT(s,u)\
if(s==1) rfftwnd_one_real_to_complex(plan_rc,(fftw_real *)u,(fftw_complex*)u);\
else rfftwnd_one_complex_to_real(plan_cr,(fftw_complex *)u,(fftw_real *)u)
and finally:
FFT(1,u0);
FFT(1,v0);
//math
...
//and transforming back
FFT(-1,u0);
FFT(-1,v0);
After moving to CUFFT:
#define OURARRAYSIZE (DIM * 2*(DIM/2+1))
#define DIM 16
cufftHandle planR2C;
cufftHandle planC2R;
cufftReal forcesX[OURARRAYSIZE];
cufftReal forcesY[OURARRAYSIZE];
cufftReal *dev_forcesX;
cufftReal *dev_forcesY;
Init:
cufftPlan2d(&planR2C, DIM, DIM, CUFFT_R2C);
cufftPlan2d(&planC2R, DIM, DIM, CUFFT_C2R);
cufftSetCompatibilityMode(planR2C, CUFFT_COMPATIBILITY_FFTW_ALL);
cufftSetCompatibilityMode(planC2R, CUFFT_COMPATIBILITY_FFTW_ALL);
cudaMalloc( (void**)&dev_forcesX, OURARRAYSIZE*sizeof(cufftReal) );
cudaMalloc( (void**)&dev_forcesY, OURARRAYSIZE*sizeof(cufftReal) );
And finally:
cufftExecR2C(planR2C, (cufftReal*) dev_forcesX, (cufftComplex*)dev_forcesX);
cufftExecR2C(planR2C, (cufftReal*) dev_forcesY, (cufftComplex*)dev_forcesY);
cudaMemcpy( forcesX, dev_forcesX, OURARRAYSIZE*sizeof(cufftReal), cudaMemcpyDeviceToHost );
cudaMemcpy( forcesY, dev_forcesY, OURARRAYSIZE*sizeof(cufftReal), cudaMemcpyDeviceToHost );
diffuseVelocity(velocitiesX, velocitiesY, forcesX, forcesY);//MATH PART
cudaMemcpy( dev_forcesX, forcesX, OURARRAYSIZE*sizeof(cufftReal), cudaMemcpyHostToDevice );
cudaMemcpy( dev_forcesY, forcesY, OURARRAYSIZE*sizeof(cufftReal), cudaMemcpyHostToDevice );
cufftExecC2R(planC2R, (cufftComplex*) dev_forcesX, (cufftReal*)dev_forcesX);
cufftExecC2R(planC2R, (cufftComplex*) dev_forcesY, (cufftReal*)dev_forcesY);
cudaMemcpy( forcesX, dev_forcesX, OURARRAYSIZE*sizeof(cufftReal), cudaMemcpyDeviceToHost );
cudaMemcpy( forcesY, dev_forcesY, OURARRAYSIZE*sizeof(cufftReal), cudaMemcpyDeviceToHost );
After the math part both programs hold exactly the same data (matrix). Sadly after the reverse fourier transformation data in matrices differs. I noticed that corrupted data is that, which lies in bonus columns ( (DIM * 2*(DIM/2+1)) ) which are needed for in place transformation.
Does anybody have any idea, why? Is there something about CUFFT that i don't know?