ここでは、CUDAを使って拡散方程式を解いてみる。
解き方はもっとも単純な差分法を使っている。
まずは単純に、diff1d.cuというテキストファイルに、
#include <stdio.h> #include <cutil.h> #include "diff1d_kernel.cu" int main( int argc, char** argv) { //デバイスの初期化 CUT_DEVICE_INIT(argc, argv); //結果書き込み用ファイルのオープン FILE *fp=fopen("result.txt","w"); //タイマーを作成して計測開始 unsigned int timer = 0; CUT_SAFE_CALL( cutCreateTimer( &timer)); CUT_SAFE_CALL( cutStartTimer( timer)); //メインメモリ上にfloat型のデータを100個生成する float* h_idata = (float*) malloc(sizeof( float) * 100); //初期条件をセット for( int i = 0; i < 50; i++) h_idata[i] = 0; for( int i = 50; i < 100; i++) h_idata[i] = 1; //デバイス上(ビデオカードのこと)にも同じくfloat型100個分のメモリを確保する float* d_idata; CUDA_SAFE_CALL( cudaMalloc( (void**) &d_idata, sizeof( float) * 100 )); dim3 grid( 1, 1, 1); //100は100個並列であるため dim3 threads(100, 1, 1); //メインメモリからデバイスのメモリにデータを転送する CUDA_SAFE_CALL( cudaMemcpy( d_idata, h_idata, sizeof( float) * 100 , cudaMemcpyHostToDevice) ); for (int t=0;t<100;t++) { //ここでGPUを使った計算が行われる diff1dKernel<<< grid, threads, sizeof( float) * 200 >>>( d_idata); //デバイスからメインメモリ上に実行結果をコピー CUDA_SAFE_CALL( cudaMemcpy( h_idata, d_idata, sizeof( float) * 100, cudaMemcpyDeviceToHost) ); //実行結果を表示 for (int i=0;i<100;i++) { fprintf(fp,"%f\t",h_idata[i]); } fprintf(fp,"\n"); } //タイマーを停止しかかった時間を表示 CUT_SAFE_CALL( cutStopTimer( timer)); printf("Processing time: %f (ms)\n", cutGetTimerValue( timer)); CUT_SAFE_CALL( cutDeleteTimer( timer)); //各種メモリを解放 free( h_idata); CUDA_SAFE_CALL(cudaFree(d_idata)); fclose(fp); //終了処理 CUT_EXIT(argc, argv); return 0; }
diff1d_kernel.cuの方はというと、
__global__ void diff1dKernel( float* g_idata) { extern __shared__ char sharedmem[]; float * sdata = (float *) sharedmem; float * sdata2 = (float *) (sharedmem+sizeof(float)*100); // スレッドIDを取得 const unsigned int tid = threadIdx.x; //グローバルメモリから入力データの読み込み sdata[tid] = g_idata[tid]; __syncthreads(); //ここで計算を行う const float Dfu=1; const float dt=0.2; const float dx=1; float Dfudtdx2=Dfu*dt/(dx*dx); for (int n=0;n<100;n++) { if (tid==0) sdata2[tid]=sdata[tid]+(2*sdata[tid+1]-2*sdata[tid])*Dfudtdx2; else if (tid==99) sdata2[tid]=sdata[tid]+(2*sdata[tid-1]-2*sdata[tid])*Dfudtdx2; else sdata2[tid]=sdata[tid]+(sdata[tid-1]+sdata[tid+1]-2*sdata[tid])*Dfudtdx2; __syncthreads(); sdata[tid]=sdata2[tid]; } //グローバルメモリに結果を書き込む g_idata[tid] = sdata[tid]; }
これを、
nvcc diff1d.cu -lcutil32 -lkernel32 -o diff1d.exe
とコンパイルし、できたdiff1d.exeを実行する。
>diff1d.exe
Processing time: 45.219727 (ms)
Press ENTER to exit...
Mac OS Xでコンパイルする場合は、
$ nvcc -L/Developer/CUDA/lib diff1d.cu -lcutil -o diff1d
のようにコンパイルする。
拡散方程式を2次元で解く場合、1次元の時とは違い系が大きいので1次元の時ような形でのShared Memoryは使えない(Shared Memoryは各ブロックあたり16kBしかないので、floatの場合で、4096個まで使用できる)。
また、前回は、ブロック数は1であったので、GPUのプロセッサをすべて使い切っていない。ここでは、複数のブロックを使ってより並列度を上げてみる。
ブロックを複数使った場合の問題は、カーネル内でブロック間の同期をとる方法が存在しない点である。そのため、下記のプログラムでは、1回計算するたびに、カーネルを終了し同期をとっている。
まず、diff2d.cuであるが、
#include <stdio.h> #include <cutil.h> #include "diff2d_kernel.cu" //系のサイズ 100 X 100のグリッド上で拡散方程式を解く const int X=100; const int Y=100; int main( int argc, char** argv) { //デバイスの初期化 CUT_DEVICE_INIT(argc, argv); //結果書き込み用ファイルのオープン FILE *fp=fopen("result.txt","w"); //タイマーを作成して計測開始 unsigned int timer = 0; CUT_SAFE_CALL( cutCreateTimer( &timer)); CUT_SAFE_CALL( cutStartTimer( timer)); //メインメモリ上にfloat型のデータをX*Y個生成する float* h_idata = (float*) malloc(sizeof( float) * X*Y); //初期条件をセット for( int i = 0; i < X; i++) for( int j = 0; j < X; j++) if((i-X/2)*(i-X/2)+(j-Y/2)*(j-Y/2)<10*10) h_idata[i*Y+j] = 1; else h_idata[i*Y+j] = 0; //デバイス上(ビデオカードのこと)にも同じくfloat型X*Y個分のメモリを確保する float* d_idata; CUDA_SAFE_CALL( cudaMalloc( (void**) &d_idata, sizeof( float) * X*Y )); //デバイス上(ビデオカードのこと)にfloat型X*Y個分の作業用メモリを確保する float* d_idata2; CUDA_SAFE_CALL( cudaMalloc( (void**) &d_idata2, sizeof( float) * X*Y )); //ブロック数を増やして並列度を上げる dim3 grid( 16, 1, 1); dim3 threads(256, 1, 1); //メインメモリからデバイスのメモリにデータを転送する CUDA_SAFE_CALL( cudaMemcpy( d_idata, h_idata, sizeof( float) * X*Y , cudaMemcpyHostToDevice) ); for (int t=0;t<100;t++) { for (int n=0;n<10;n++) { //ここでGPUを使った計算が行われる diff2dKernel<<< grid, threads>>>( d_idata, d_idata2,X,Y); //作業用領域から書き戻す CUDA_SAFE_CALL( cudaMemcpy( d_idata, d_idata2, sizeof( float) * X*Y, cudaMemcpyDeviceToDevice) ); } //デバイスからメインメモリ上に実行結果をコピー CUDA_SAFE_CALL( cudaMemcpy( h_idata, d_idata, sizeof( float) * X*Y, cudaMemcpyDeviceToHost) ); //実行結果を表示 for (int i=0;i<X;i+=2) { for (int j=0;j<Y;j+=2) { fprintf(fp,"%f\t",h_idata[i*Y+j]); } } fprintf(fp,"\n"); } //タイマーを停止しかかった時間を表示 CUT_SAFE_CALL( cutStopTimer( timer)); printf("Processing time: %f (ms)\n", cutGetTimerValue( timer)); CUT_SAFE_CALL( cutDeleteTimer( timer)); //各種メモリを解放 free( h_idata); free( h_idata2); CUDA_SAFE_CALL(cudaFree(d_idata)); fclose(fp); //終了処理 CUT_EXIT(argc, argv); return 0; }
diff2d_kernel.cuは、
__global__ void diff2dKernel(float* g_idata, float* g_idata2,int X, int Y) { const unsigned int tid = threadIdx.x; const unsigned int bid = blockIdx.x; const unsigned int bdim = blockDim.x; const unsigned int gdim = gridDim.x; int step=bdim*gdim; int num=X*Y; //ここで計算を行う const float Dfu=1; const float dt=0.2; const float dx=1; float Dfudtdx2=Dfu*dt/(dx*dx); for (int id=bid * bdim + tid;id<num;id+=step) { //境界は、nonflux境界条件 if ((id)==0)//四隅 g_idata2[id]=g_idata[id]+(g_idata[id+1]+g_idata[id+1]+g_idata[id+Y]+g_idata[id+Y]-4*g_idata[id])*Dfudtdx2; else if ((id)==Y-1)//四隅 g_idata2[id]=g_idata[id]+(g_idata[id-1]+g_idata[id-1]+g_idata[id+Y]+g_idata[id+Y]-4*g_idata[id])*Dfudtdx2; else if ((id)==X*Y-Y)//四隅 g_idata2[id]=g_idata[id]+(g_idata[id+1]+g_idata[id+1]+g_idata[id-Y]+g_idata[id-Y]-4*g_idata[id])*Dfudtdx2; else if ((id)==X*Y-1)//四隅 g_idata2[id]=g_idata[id]+(g_idata[id-1]+g_idata[id-1]+g_idata[id-Y]+g_idata[id-Y]-4*g_idata[id])*Dfudtdx2; else if ((id)<Y)//四辺 g_idata2[id]=g_idata[id]+(g_idata[id-1]+g_idata[id+1]+g_idata[id+Y]+g_idata[id+Y]-4*g_idata[id])*Dfudtdx2; else if((id)>X*Y-Y)//四辺 g_idata2[id]=g_idata[id]+(g_idata[id-1]+g_idata[id+1]+g_idata[id-Y]+g_idata[id-Y]-4*g_idata[id])*Dfudtdx2; else if ((id)%Y==0)//四辺 g_idata2[id]=g_idata[id]+(g_idata[id+1]+g_idata[id+1]+g_idata[id-Y]+g_idata[id+Y]-4*g_idata[id])*Dfudtdx2; else if ((id)%Y==Y-1)//四辺 g_idata2[id]=g_idata[id]+(g_idata[id-1]+g_idata[id-1]+g_idata[id-Y]+g_idata[id+Y]-4*g_idata[id])*Dfudtdx2; else g_idata2[id]=g_idata[id]+(g_idata[id-1]+g_idata[id+1]+g_idata[id-Y]+g_idata[id+Y]-4*g_idata[id])*Dfudtdx2; } }
上記のプログラムを下記のようにしてコンパイルする。
>nvcc diff2d.cu -lcutil32 -lcuda -o diff2d.exe
そして実行する。当方の環境では、下記のような実行時間となった。
>diff2d.exe
Processing time: 762.077759 (ms)
Press ENTER to exit...
結果は下記の通り。
この結果の表示にはMatlabを使った。一応そのコード。
data=load('result.txt'); for lp=1:9 subplot(3,3,lp) z=data(lp*5,:); z=reshape(z,50,50); surf(z) shading interp daspect([50 50 1]) axis([1,50,1,50,0,1]) caxis([0 1]) end
最終更新日