私の論文では、特別な MPI-Navier Stokes-Solver プログラムを CUDA で最適化する必要があります。元のプログラムは、いくつかの PDE を解くために FFTW を使用します。詳細には、いくつかの上三角行列を 2 次元でフーリエ変換しますが、1 次元配列として扱います。今のところ、元のコードの一部に苦労しています: (N は常に 64 に設定されています)
オリジナル:
//Does the complex to real in place fft an normalizes
void fftC2R(double complex *arr) {
fftw_execute_dft_c2r(plan_c2r, (fftw_complex*)arr, (double*)arr);
//Currently ignored: Normalization
/* for(int i=0; i<N*(N/2+1); i++)
arr[i] /= (double complex)sqrt((double complex)(N*N));*/
}
void doTimeStepETDRK2_nonlin_original() {
//calc velocity
ux[0] = 0;
uy[0] = 0;
for(int i=1; i<N*(N/2+1); i++) {
ux[i] = I*kvec[1][i]*qvec[i] / kvec[2][i];
uy[i] = -I*kvec[0][i]*qvec[i] / kvec[2][i];
}
fftC2R(ux);
fftC2R(uy);
//do some stuff here...
//...
return;
}
ux と uy は次のように割り当てられます (double complex arrays):
ux = (double complex*)fftw_malloc(N*(N/2+1) * sizeof(double complex));
uy = (double complex*)fftw_malloc(N*(N/2+1) * sizeof(double complex));
fft プランは次のように作成されます。
plan_c2r = fftw_plan_dft_c2r_2d(N, N,(fftw_complex*) qvec, (double*)qvec, FFTW_ESTIMATE);
qvec は ux や uy と同じ方法で割り当てられ、データ型は double complex です。
CUDA コードの関連部分は次のとおりです。
NN2_VecSetZero_and_init<<<block_size,grid_size>>>();
cudaSafeCall(cudaDeviceSynchronize());
cudaSafeCall(cudaGetLastError());
int err = (int)cufftExecZ2D(cu_plan_c2r,(cufftDoubleComplex*)sym_ux,(cufftDoubleReal*)sym_ux);
if (err != CUFFT_SUCCESS ) {
exit(EXIT_FAILURE);
return;
}
err = (int)cufftExecZ2D(cu_plan_c2r,(cufftDoubleComplex*)sym_uy,(cufftDoubleReal*)sym_uy);
if (err != CUFFT_SUCCESS ) {
exit(EXIT_FAILURE);
return;
}
//do some stuff here...
//...
return;
sim_ux と sim_uy は次のように割り当てられます。
cudaMalloc((void**)&sym_ux, N*(N/2+1)*sizeof(cufftDoubleComplex));
cudaMalloc((void**)&sym_uy, N*(N/2+1)*sizeof(cufftDoubleComplex));
関連するカフパーツの初期化は次のようになります
if (cufftPlan2d(&cu_plan_c2r,N,N, CUFFT_Z2D) != CUFFT_SUCCESS){
exit(EXIT_FAILURE);
return -1;
}
if (cufftPlan2d(&cu_plan_r2c,N,N, CUFFT_D2Z) != CUFFT_SUCCESS){
exit(EXIT_FAILURE);
return -1;
}
if ( cufftSetCompatibilityMode ( cu_plan_c2r , CUFFT_COMPATIBILITY_FFTW_ALL) != CUFFT_SUCCESS ) {
exit(EXIT_FAILURE);
return -1;
}
if ( cufftSetCompatibilityMode ( cu_plan_r2c , CUFFT_COMPATIBILITY_FFTW_ALL) != CUFFT_SUCCESS ) {
exit(EXIT_FAILURE);
return -1;
}
そのため、FFTW の完全な互換性を使用し、すべての関数を FFTW 呼び出しパターンで呼び出します。両方のバージョンを実行すると、ux と uy (sim_ux と sim_uy) でほぼ同じ結果が得られます。しかし、配列の定期的な位置では、Cufft はこれらの要素を無視しているように見えます。ここで、FFTW はこれらの要素の実部をゼロに設定し、複素数部分を計算します (配列が大きすぎてここに表示できません)。これが発生する歩数は N/2+1 です。ですから、Cufft と FFTW の間の fft パディング理論を完全には理解していなかったと思います。Cufft-executions が呼び出されるまで、これらの配列間の以前の不一致を除外できます。したがって、上記のコードの他の配列はここでは関係ありません。
私の質問は、FFTW 呼び出しスタイルをほぼ 100% 使用することに楽観的すぎるのでしょうか? fft の前に配列を処理する必要がありますか? Cufft のドキュメントによると、データの入力配列と出力配列のサイズを変更する必要があります。しかし、インプレース変換を実行しているときに、どうすればこれを行うことができますか? 元のコードから離れすぎたくないし、各 fft 呼び出しにこれ以上コピー命令を使用したくありません。これは、メモリが限られているため、配列は gpu でできるだけ長く保持して処理する必要があるためです。
すべてのヒント、批判的な声明、またはアイデアに感謝します!
私の構成:
- コンパイラ: gcc 4.6 (C99 標準)
- MPI-Package: mvapich2-1.5.1p1 (単一処理デバッグ モードが縮小されているため、役割を果たさないはずです)
- CUDA バージョン: 4.2
- GPU: CUDA-arch-compute_20 (NVIDIA GeForce GTX 570)
- FFTW3.3