Project EulerをGPUで（Problem 1）

Project Eulerの問題をGPU（CUDA）で解いてみるシリーズ。卒論等が本格的に始まる前に、CUDAの演習がてら始めてみる。
(とはいえ、アルゴリズム的にGPUで解ける問題は多くないので悩みどころ。)

とりあえずProblem 1から。

初期の問題は、特に計算スピードが要求されることもない（CPUでナイーブなアルゴリズムを使っても、一瞬で終わる）ので、特に速度の比較等はしていない。体感的にはGPUの方が明らかに遅い。というか、この領域だとGPUカーネル呼び出しのオーバーヘッドが大きくて、GPUでやるとかえって遅くなる。

というわけでここら辺は完全に手習いの領域。

アルゴリズムの説明

CPUで計算する場合と違い、GPUはSIMDなので、特殊なアルゴリズムを使う必要がある。具体的には、多数のスレッド（数百〜数千）で、配列の上をデータ並列に計算する。

今回は、長さ1000の配列を用意し、

1. 1000個のスレッドで、第i要素を 0 または i に初期化(3または5で割れるかどうかのチェック)
2. 配列の要素を、2個ごとに足し合わせてリダクションしていく
3. 第0要素に答えが入っている！

というアルゴリズムにしてみた。

アルゴリズムやGPUの原理の詳細は、青木・額田本を参考にしてほしい。

はじめてのCUDAプログラミング―驚異の開発環境[GPU+CUDA]を使いこなす! (I・O BOOKS)
	青木 尊之 額田 彰 第二I O編集部 工学社 2009-11 売り上げランキング : 27832 おすすめ平均 CUDAプログラミングの概要を知るには良い本だと思います。 Amazonで詳しく見る by G-Tools

ソースコードと結果

// Project Euler Problem 001 with cuda
// Written by Keisuke Fukuda, All rights reserved.
// This program is licensed under the The BSD License.

#include <cstdio>
#include <cstdlib>
#include <cmath>
#include <cassert>

#define CUDA_CALL(stmt)  assert( (stmt) == cudaSuccess )

#ifndef NUM
#  define NUM 1000
#endif

__global__ void kernel_init_array(int* array) {
    int id = blockIdx.x * 100 + threadIdx.x;

    // 配列の id 番目要素を、3または5の倍数なら id 、それ以外なら０にする
    if (id % 3 == 0 || id % 5 == 0) {
        array[id] = id;
    }
    else {
        array[id] = 0;
    }
}

// NUM以下の3または5の倍数の和を求める1ループを実行する kernel
__global__ void kernel_prob_001(int *array, int base) {

    int id = blockIdx.x * 100 + threadIdx.x;

    int j = id * base;

    if( j + base/2 < NUM ) {
        array[j] = array[j] + array[j + base/2];
    }
}

int main(int argc, char* argv[]) {
    int answer = 0, *gpu_mem = NULL;

    CUDA_CALL( cudaMalloc((void**)&gpu_mem, sizeof(int) * NUM) );
    assert(gpu_mem);

    // kernelの呼び出し
    dim3 gridDim(NUM > 100 ? NUM/100 : 1);
    dim3 blockDim(NUM > 10 ? 100 : NUM);

    kernel_init_array<<<gridDim, blockDim>>>(gpu_mem);
    cudaThreadSynchronize();

    int base = 1;
    do {
        base *= 2;
        kernel_prob_001<<<gridDim, blockDim>>>(gpu_mem, base);
        cudaThreadSynchronize();
    } while(base < NUM);


    CUDA_CALL( cudaMemcpy(&answer, gpu_mem, sizeof(int),       cudaMemcpyDeviceToHost) );

    CUDA_CALL( cudaFree(gpu_mem) );

    printf("Answer : %d\n", answer);
}

実行結果：

$ nvcc euler001.cu && ./a.out
Answer : 233168