Shared Memoryは各ブロックあたり16kBしかないので、floatの場合で、4096個しか使用できない。
一方、Constant Memoryは64kBあるので、その4倍の16384個（floatの場合）を使うことができる。
Constant Memoryはカーネル内からは書き換えられないが、CPU側からはいつでも書き換えることができる。
また、アクセス速度もGlobal Memoryより速い。
今回は、以前の2次元拡散方程式を、Constant Memoryを使って解いてみる。
まず、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_idata2;
CUDA_SAFE_CALL( cudaMalloc( (void**) &d_idata2, sizeof( float) * X*Y ));
//ブロック数を増やして並列度を上げる
dim3 grid( 16, 1, 1);
dim3 threads(256, 1, 1);
//メインメモリからデバイスのConstantメモリにデータを転送する
CUDA_SAFE_CALL( cudaMemcpyToSymbol( g_idata, h_idata, sizeof( float) * X*Y) );
for (int t=0;t<100;t++)
{
for (int n=0;n<10;n++)
{
//ここでGPUを使った計算が行われる
diff2dKernel<<< grid, threads>>>( d_idata2,X,Y);
//作業用領域から書き戻す
CUDA_SAFE_CALL( cudaMemcpyToSymbol( g_idata, d_idata2, sizeof( float) * X*Y,0, cudaMemcpyDeviceToDevice) );
}
//デバイスからメインメモリ上に実行結果をコピー
CUDA_SAFE_CALL( cudaMemcpy( h_idata, d_idata2, 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);
CUDA_SAFE_CALL(cudaFree(d_idata2));
fclose(fp);
//終了処理
CUT_EXIT(argc, argv);
return 0;
}
diff2d_kernel.cuは、
__constant__ float g_idata[10000];
__global__ void diff2dKernel(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: 561.552124 (ms)
Press ENTER to exit...
Global Memoryを使ったときの結果、762.077759 msよりも高速化されていることがわかる。
Shared Memoryよりも大きいとはいえ、たかだか16384個しかないので、いつも使える方法というわけではないが、Global Memoryを使った場合と比べて、有為な差が出るので必要に応じて使える方法である。