多体問題(2)

前回遅かったのは、グローバルな配列に何度もアクセスしていたためと思われる。共有メモリというものがあって、こちらは速いので、これを使うことを考える。
説明しやすくするために、4096個の質点に0〜4095のIDを振っておこう。共有メモリは、ブロックの中ですべてのスレッドで共有できるメモリである。ブロックサイズは256なので、一つのブロックでは256個の質点それぞれに作用する重力を計算する。そのため4096個の質点の座標が必要になる。
しかし、共有メモリは16kBという制限がある。floatなら4096個になる。座標が2つあるので、すべてをここに入れることはできない。2048個だとぎりぎり入りそうだが、多少は別に使っているらしく、無理らしい。というわけで、質点1024個ずつグローバルから共有メモリにコピーして、その質点との相互作用を計算する。
まず、0〜1023の質点の座標をコピーする。そのときに、256個のスレッドはダブらないように4つずつ分担する。すなわち、スレッド0は0〜3の質点、スレッド1は4〜7などと分担する。こうして共有メモリにコピーした座標を使って、0〜1023の質点から及ぼされる重力を計算する。次に、1024〜2047の質点と同様の計算をする。


このようにすると、CPUで235秒かかっていたものが、3.43秒となり、そんなに悪くない速度になった。



#include
#include
#include
#include

#define BLOCK_DIM 16
#define SHARED_DIM 1024

__device__ double diff_x(float x1, float x0) {
float dx = x1 - x0;
if(dx < -0.5)
dx += 1.0;
else if(0.5 < dx)
dx -= 1.0;
return dx;
}

__device__ void calcAcceleration(float *data,
int nParticles, int k, dim3 tid) {
float x = data[k*6];
float y = data[k*6+1];
float ax = 0.0;
float ay = 0.0;

__shared__ float sx[SHARED_DIM];
__shared__ float sy[SHARED_DIM];

for(int i0 = 0; i0 < nParticles; i0 += SHARED_DIM) {
// 共有メモリに自分の分担分のデータをコピーする
int iThread = tid.x + BLOCK_DIM * tid.y;
const int size = SHARED_DIM / BLOCK_DIM / BLOCK_DIM;
int j0 = iThread * size;
for(int j = j0; j < j0 + size; j++) {
sx[j] = data[(i0+j)*6];
sy[j] = data[(i0+j)*6+1];
}

__syncthreads(); // すべてコピーされるのを待つ

for(int i = 0; i < SHARED_DIM; i++) {
if(i0 + i == k)
continue;

float dx = diff_x(sx[i], x);
float dy = diff_x(sy[i], y);
float dist_square = dx * dx + dy * dy;
float d = dist_square * sqrt(dist_square);
ax += dx / d;
ay += dy / d;
}

__syncthreads(); // すべてのスレッドの計算が終わるのを待つ
}

data[k*6+4] = ax;
data[k*6+5] = ay;
}

// [ 0, 1 )に
__device__ void normalize(float& x) {
if(x < 0.0)
x += 1.0;
else if(x >= 1.0)
x -= 1.0;
}

__global__ void calc(float *data, int width, int height,
int nParticles, float dt) {
unsigned int xIndex = blockIdx.x * BLOCK_DIM + threadIdx.x;
unsigned int yIndex = blockIdx.y * BLOCK_DIM + threadIdx.y;
if( (xIndex < width) && (yIndex < height)) {
unsigned int index = yIndex * width + xIndex;
calcAcceleration(data, nParticles, index, threadIdx);
}
}

__global__ void update(float *data, int width, int height, float dt) {
unsigned int xIndex = blockIdx.x * BLOCK_DIM + threadIdx.x;
unsigned int yIndex = blockIdx.y * BLOCK_DIM + threadIdx.y;
if( (xIndex < width) && (yIndex < height)) {
unsigned int index = yIndex * width + xIndex;
int k = index * 6;
data[k+2] += data[k+4] * dt;
data[k+3] += data[k+5] * dt;
data[k ] += data[k+2] * dt;
data[k+1] += data[k+3] * dt;
normalize(data[k]);
normalize(data[k+1]);
}
}

void runTest( int argc, char** argv);

int main( int argc, char** argv) {
runTest( argc, argv);
cutilExit(argc, argv);
}

void runTest( int argc, char** argv) {
const unsigned int size_x = 64;
const unsigned int size_y = 64;
const unsigned int nParticles = size_x * size_y;
const unsigned int nData = nParticles * 6;
const unsigned int mem_size = sizeof(float) * nData;
const float dt = 0.0001;
const int nRep = 1000;

double x[nRep], y[nRep]; // ID0の粒子の位置

unsigned int timer;
cutCreateTimer(&timer);

if( cutCheckCmdLineFlag(argc, (const char**)argv, "device") )
cutilDeviceInit(argc, argv);
else
cudaSetDevice( cutGetMaxGflopsDeviceId() );

// 初期化
float *h_data = (float *)malloc(mem_size);
for(int i = 0; i < nParticles; i++) {
for(int k = 0; k < 2; k++) {
h_data[i*6+k] = (double)rand() / RAND_MAX; // 位置
h_data[i*6+k+2] = 0.0; // 速度
}
}

cutStartTimer(timer);

// データをGPUに転送
float *d_data;
cutilSafeCall(cudaMalloc( (void**)&d_data, mem_size));
cutilSafeCall(cudaMemcpy(d_data, h_data, mem_size, cudaMemcpyHostToDevice));

// setup execution parameters
dim3 grid(size_x / BLOCK_DIM, size_y / BLOCK_DIM, 1);
dim3 threads(BLOCK_DIM, BLOCK_DIM, 1);

for(int i = 0; i < nRep; i++) {
// GPUで計算実行
calc<<< grid, threads >>>(d_data, size_x, size_y, nParticles, dt);
cudaThreadSynchronize();

cutilCheckMsg("Kernel execution failed");

update<<< grid, threads >>>(d_data, size_x, size_y, dt);
cudaThreadSynchronize();

cutilCheckMsg("Kernel execution failed");

// データをGPUから転送
cutilSafeCall(cudaMemcpy(h_data, d_data, mem_size, cudaMemcpyDeviceToHost));
x[i] = h_data[0];
y[i] = h_data[1];
}
float optimizedTime = cutGetTimerValue(timer);
cutStopTimer(timer);

for(int i = 0; i < nRep; i++) {
printf("%5.2f %.9f %.9f\n", dt * i, x[i], y[i]);
}

printf("time: %0.3f ms\n\n", optimizedTime);

// cleanup memory
free(h_data);
cutilSafeCall(cudaFree(d_data));
cutilCheckError( cutDeleteTimer(timer));

cudaThreadExit();
}