CUDA技術を利用したGPUコンピューティングの実際（後編） ―― FFTを利用した光波の伝播（フレネル回折）をGPUで高速計算

～　途中省略　～

#include <cutil.h>　//CUTIL用のヘッダ・ファイル
#include <cufft.h>　//CUFFT用のヘッダ・ファイル

～　途中省略　～

  int N=512;
  int N2=N*2;

  double p=10.0e-6; //サンプリング間隔
  double lambda=633.0e-9; //赤色レーザ波長
  double d=0.5; //開口面と観察面の間の距離
  double PI=3.1415926535897932384626433832795; //円周率

  //メモリ確保
  int mem_size = sizeof(cufftComplex) * N2 * N2;
  cufftComplex* h_buf = (cufftComplex*) malloc( mem_size);

  cufftComplex* d_buf1, *d_buf2, *d_buf3;
  cudaMalloc( (void**) &d_buf1, mem_size);
  cudaMalloc( (void**) &d_buf2, mem_size);
  cudaMalloc( (void**) &d_buf3, mem_size);

  //開口データa(x,y)の作成と拡張　←(1)
  for(int i=0;i<N2;i++){
    for(int j=0;j<N2;j++){
      int adr=j+i*N2;
      if(j>256-20 && j<256+20 && i>256-10 && i<256+10){
        h_buf[adr]=make_cuComplex(1.0, 0.0);
      }
      else{
        h_buf[adr]=make_cuComplex(0.0, 0.0);
    }
  }
  }

  //開口面a(x,y)のフーリエ変換　←(2)
  cufftHandle fftplan;
  cufftPlan2d(&fftplan, N2, N2, CUFFT_C2C);
  cudaMemcpy( d_buf1, h_buf, mem_size,cudaMemcpyHostToDevice);
  cufftExecC2C(fftplan, d_buf1, d_buf2, CUFFT_FORWARD);

  //p(x,y)を算出　←(3)
  for(int i=0;i<N2;i++){
    for(int j=0;j<N2;j++){
      int adr=j+i*N2;
      double dx=((double)( (j-N2/2) )) * p;
      double dy=((double)( (i-N2/2) )) * p;
      double tmp= 1.0 * PI* ( dx*dx + dy*dy ) / (lambda * d);
      h_buf[adr]=make_cuComplex(cos(tmp),sin(tmp));
    }
  }

  FFTShift(h_buf, N2, N2); //象限の交換

  //p(x,y)のFFT　←(4)
  cudaMemcpy( d_buf3, h_buf, mem_size,cudaMemcpyHostToDevice);
  cufftExecC2C(fftplan, d_buf3, d_buf1, CUFFT_FORWARD);

  //複素乗算　←(5)
  dim3 block(16, 16, 1);
  dim3 grid(N2 / block.x, N2 / block.y, 1);
  KernelMult<<< grid, block>>>(d_buf2, d_buf1, d_buf3, N2, N2);

  //逆FFT　←(6)
  cufftExecC2C(fftplan, d_buf3, d_buf1, CUFFT_INVERSE);

  //計算結果を転送　←(7)
  cudaMemcpy( h_buf, d_buf1, mem_size,cudaMemcpyDeviceToHost);

  //結果をビットマップ・ファイルとして保存　←(8)
  Save("diffraction.bmp", N2, N2, h_buf);

  //開放　←(9)
  cufftDestroy(fftplan);
  cudaFree(d_buf1);
  cudaFree(d_buf2);
  cudaFree(d_buf3);
  free(h_buf);