ガウスの故事にならって、1000個ずつ区切って素数の個数を数える。1001〜2000に135個、2001〜3000に127個ある。nを1500、2500、…として、横軸を1/log n、縦軸を素数の個数としてプロットした。
素数の密度は、ほぼ1/log nになっていることがわかる。これを素数定理という（厳密には違うが）。
素数の個数を数えるプログラムをそのままGPUに使うとあまり速度が出ない。
__global__ void count_prime(int *data, int width, int height, int nRep) {
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 x = data[index];
int result = 0;
if(x > 0) {
for(int i = x; i <= x + nRep; i += 2) {
for(int p = 3; p * p <= i; p += 2) {
if(i % p == 0) {
result--;
break;
}
}
result++;
}
}
data[index] = result;
}
}
一般にGPUは分岐に弱い。これは、GPUによって違うが、ある程度のまとまったスレッドは同じ命令を実行するという制約があるためである。多少の分岐はなんとかしてくれるが、特に上のようにループが入れ子になると遅くなるようである。これを多少工夫して、一重のループにすると速くなる。
if(x > 0) {
int i = x;
int p = 3;
while(true) {
if(i >= x + nRep)
break;
if(p * p > i) {
result++;
i += 2;
p = 3;
}
else {
if(i % p == 0) {
i += 2;
p = 3;
}
else {
p += 2;
}
}
}
}
CPU GPU(上) GPU(下)
〜1e7 2856ms 1063ms 285ms
〜2e7 7626ms 2829ms 636ms
〜4e7 20654ms 6877ms 1565ms
CPUに対して、単純にGPUで計算させると3倍程度にしか速くならないが、下のようにすると14倍程度になった。
// prime.cpp
#include
#includeint main() {
const unsigned int size_x = 256;
const unsigned int size_y = 256;
const unsigned int nData = size_x * size_y;
const unsigned int nIntervals = 10000;
const unsigned int mem_size = sizeof(int) * nData;
const int nRep = 1000;
int *h_data = (int *)malloc(nData * sizeof(int));
for(int i = 1; i < nData; i++) {
h_data[i] = i < nIntervals ? i * 1000 + 1 : 0;
}
LARGE_INTEGER nFreq, nStart, nEnd;
QueryPerformanceFrequency(&nFreq);
QueryPerformanceCounter(&nStart);
for(int i = 1; i < nData; i++) {
int x = h_data[i];
int result = 0;
if(x > 0) {
for(int i = x; i < x + nRep; i += 2) {
for(int p = 3; p * p <= i; p += 2) {
if(i % p == 0) {
result--;
break;
}
}
result++;
}
}
h_data[i] = result;
}
QueryPerformanceCounter(&nEnd);
for(int i = 1; i <= nIntervals; i++) {
printf("%7d 〜 %8d : %3d\n", i * 1000 + 1, (i + 1) * 1000, h_data[i-1]);
}
printf("time: %0.3f ms\n\n",
(nEnd.QuadPart - nStart.QuadPart) * 1e3 / nFreq.QuadPart);
free(h_data);
}
// prime2.cu
#include
#include
#include
#include#define BLOCK_DIM 16
__global__ void count_prime(int *data, int width, int height, int nRep) {
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 x = data[index];
int result = 0;
if(x > 0) {
int i = x;
int p = 3;
while(true) {
if(i >= x + nRep)
break;
if(p * p > i) {
result++;
i += 2;
p = 3;
}
else {
if(i % p == 0) {
i += 2;
p = 3;
}
else {
p += 2;
}
}
}
}
data[index] = result;
}
}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 = 256;
const unsigned int size_y = 256;
const unsigned int nData = size_x * size_y;
const unsigned int nIntervals = 10000;
const unsigned int mem_size = sizeof(int) * nData;
const int nRep = 1000;
unsigned int timer;
cutCreateTimer(&timer);
if( cutCheckCmdLineFlag(argc, (const char**)argv, "device") )
cutilDeviceInit(argc, argv);
else
cudaSetDevice( cutGetMaxGflopsDeviceId() );
// 初期化
int *h_data = (int *)malloc(nData * sizeof(int));
h_data[0] = 0;
for(int i = 1; i < nData; i++) {
h_data[i] = i < nIntervals ? i * 1000 + 1 : 0;
}
cutStartTimer(timer);
// データをGPUに転送
int *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);
// GPUで計算実行
count_prime<<< grid, threads >>>(d_data, size_x, size_y, nRep);
cudaThreadSynchronize();
cutilCheckMsg("Kernel execution failed");
// データをGPUから転送
cutilSafeCall(cudaMemcpy(h_data, d_data, mem_size, cudaMemcpyDeviceToHost));
float optimizedTime = cutGetTimerValue(timer);
cutStopTimer(timer);
for(int i = 1; i < nIntervals; i++) {
printf("%7d 〜 %8d : %3d\n", i * 1000 + 1, (i + 1) * 1000, h_data[i]);
}
printf("time: %0.3f ms\n\n", optimizedTime);
// cleanup memory
free(h_data);
cutilSafeCall(cudaFree(d_data));
cutilCheckError( cutDeleteTimer(timer));
cudaThreadExit();
}
