CUDAで遊んでみた

nVIDIAのGPUで並列計算を行うプログラムをCで書けるっていうCUDAで遊んでみました。
ということで、みしょさんの11以下の素因数を持つ携帯電話番号には価値がない。でのお題をやってみました。
いくつかあるのですが、このお題。

僕の携帯電話番号は113未満の素因数を持ちません。これはどの程度すばらしい電話番号なのでしょう？つまり，80から始まる10桁の自然数について，113未満の素因数を持たないものはどれくらいあるでしょうか？

と思ったら、なんだかCUDAでは大きい数字がうまく扱えない様子。
なので、10から始まる9桁の自然数について、113未満の素因数を持たないものを数えてみることにしました。
実行結果は次のとおり
ホストが、CPUでやった結果です。
Kernel1が共有メモリを使わずにGPUでやった場合。Kernel2が共有メモリ使った場合。

100000000から109999999までで113より小さい素因数を持たないものの数
ホスト 計算時間=1326.090088(ms) 結果 1157751
Kernel1 計算時間=6846.341309(ms) 計算結果=1157751
Kernel2 計算時間=5039.050293(ms) 計算結果=1157751

ぐぬぬ！GPU使ったほうが遅くなってます。共有メモリ使うとかなりスピードが改善してることもわかりますね。
ということで、整数演算だと、CPUよりも遅いかもしれない。
使ったのがGeForce8400GSという、最廉価のバージョンなので、8800GTとかを使うとCPUよりも速くなるのかもしれません。
ちなみに、Javaだと1556msになってます。Cよりちょっと遅いですが、数回繰り返すようにして最適化がかかればCと同じくらいの速さになりそうですね。

一応、80で始まる10桁を求めると、11579542になりました。Javaだと31秒、Cだと15秒程度でした。64ビット整数の演算は、Javaは苦手なんでしょか？

ということでソース。
Java版

import java.util.*;

public class CountPrime {
    public static void main(String[] args){
        //素数を求める
        boolean[] primeFlag = new boolean[113];
        for(int i = 0; i < primeFlag.length; ++i) primeFlag[i] = true;
        primeFlag[1] = false;
        for(int i = 1; i < primeFlag.length; ++i){
            if(!primeFlag[i]) continue;
            for(int j = i * 2; j < primeFlag.length; j += i){
                primeFlag[j] = false;
            }
        }
        //素数の配列を作成
        List<Integer> primsList = new ArrayList<Integer>();
        for(int i = 1; i < primeFlag.length; ++i){
            if(primeFlag[i]) primsList.add(i);
        }
        int[] prims = new int[primsList.size()];//Listのままだと65秒 配列だと31秒
        for(int i = 0; i < primsList.size(); ++i){
            prims[i] = primsList.get(i);
        }
        
        //
        long start = System.nanoTime();
        int count = 0;
        LOOP:
        //for(long lg = 8000000000L; lg < 8100000000L; ++lg){
        for(int lg = 100000000; lg < 110000000; ++lg){
            for(int i : prims){
                if(lg % i == 0) continue LOOP;
            }
            count++;
        }
        System.out.println((System.nanoTime() - start) / 1000 / 1000);
        System.out.println(count);
        // 090の場合 11579520
        // 080の場合 11579542
        // 070の場合 11579562
    }
}

で、CUDA版

#include <stdio.h>
#include <conio.h>
#include <cutil.h>


#define BLOCK 32
#define START 100000000L
#define RANGE  10000000L
#define MOD(X,Y) (X)%(Y)
#define PRIM 113
#define LONG unsigned __int32

//プロトタイプ
void Host(long *a, int count);
__global__ void Kernel1(LONG *a, LONG *b, int count);
__global__ void Kernel2(LONG *a, LONG *b, int count);

void SetTimer(unsigned int *t);
float EndTimer(unsigned int *t);


int main()
{
	printf("%ldから%ldまでで%dより小さい素因数を持たないものの数\n", START, START + RANGE - 1, PRIM);
	//素数を求める
	bool primflag[PRIM];
	for(int i = 1; i < PRIM; ++i){
		primflag[i] = true;
	}
	primflag[1] = false;
	for(int i = 1; i < PRIM; ++i){
		if(!primflag[i]) continue;
		for(int j = i * 2; j < PRIM; j += i){
			primflag[j] = false;
		}
	}
	//配列に格納
	int count = 0;
	LONG prims[PRIM];
	for(int i = 1; i < PRIM; ++i){
		if(!primflag[i]) continue;
		prims[count++] = i;
	}
	
	unsigned int timer;
	CUT_DEVICE_INIT();
	LONG pcount = 0;
	//ホスト側で計算
	
	SetTimer(&timer);

	pcount = 0;
	for(LONG lg = START; lg < START + RANGE; ++lg){
		bool flag = true;
		for(int i = 0; i < count; ++i){
			if(MOD(lg, prims[i]) == 0){
				flag = false;
				break;
			}
		}
		if(flag) pcount++;
	}


	printf("ホスト 計算時間=%f(ms)", EndTimer(&timer));
	printf(" 結果 %d\n", pcount);
	

	//初期化
	LONG *d_a;
	LONG *d_b;
	LONG ret[PRIM * BLOCK];
	dim3 grid(BLOCK, 1, 1);
	dim3 threads(count, 1, 1);

	//=== GPU ===
	for(int k = 0; k < 2; ++k){
	SetTimer(&timer);

	//CUDAグローバルメモリの確保
	cudaMalloc((void**) &d_a, sizeof(LONG) * count);
	cudaMalloc((void**) &d_b, sizeof(LONG) * count * BLOCK);
	cudaMemset(d_b, 0, sizeof(LONG) * count * BLOCK);

	//CUDAへコピー
	cudaMemcpy(d_a, prims, sizeof(LONG) * count,
		cudaMemcpyHostToDevice);

	//実行
	switch(k){
		case 0:
			Kernel1 <<< grid, threads >>>(d_a, d_b, count);
			break;
		case 1:
			Kernel2 <<< grid, threads >>>(d_a, d_b, count);
			break;
	}

	//CUDAからコピー
	cudaMemcpy(ret, d_b,
		sizeof(LONG) * count * BLOCK, cudaMemcpyDeviceToHost);
	//集計
	pcount = 0;
	for(int i = 0; i < count * BLOCK; ++i){
		pcount += ret[i];
	}

	printf("Kernel%d 計算時間=%f(ms)", k + 1, EndTimer(&timer));
	printf(" 計算結果=%d\n", pcount);

	//メモリ解放
	cudaFree(d_a);
	cudaFree(d_b);
	}

	_getch();
}

__global__ void Kernel1(LONG *a, LONG *b, int count)
{
	int idx = blockIdx.x * count + threadIdx.x;
	LONG split = RANGE / count / BLOCK;
	LONG s = START + idx * split;
	//LONG e = START + (idx + 1) * split;
	LONG e = (idx + 1 == BLOCK * count) ? START + RANGE : START + split * (idx + 1);
	//printf("%d:%d %lld - %lld\n", blockIdx.x, threadIdx.x, s, e);
	LONG pcount = 0;
	for(LONG lg = s; lg < e; ++lg){
		bool flag = true;
		for(int i = 0; i < count; ++i){
			if(MOD(lg, a[i]) == 0){
				flag = false;
				break;
			}
		}
		if(flag) pcount++;
	}
	b[idx] = pcount;

}

__global__ void Kernel2(LONG *a, LONG *b, int count)
{
	int idx = blockIdx.x * count + threadIdx.x;
	LONG split = RANGE / count / BLOCK;
	LONG s = START + idx * split;
	LONG e = (idx + 1 == BLOCK * count) ? START + RANGE : START + split * (idx + 1);
	//printf("%d:%d %lld - %lld\n", blockIdx.x, threadIdx.x, s, e);

	__shared__ LONG prim[PRIM];
	prim[threadIdx.x] = a[threadIdx.x];
	__syncthreads();

	LONG pcount = 0;
	for(LONG lg = s; lg < e; ++lg){
		bool flag = true;
		for(int i = 0; i < count; ++i){
			if(MOD(lg, prim[i]) == 0){
				flag = false;
				break;
			}
		}
		if(flag) pcount++;
	}
	b[idx] = pcount;

}


void SetTimer(unsigned int *t){
	cutCreateTimer(t);
	cutStartTimer(*t);
}
float EndTimer(unsigned int *t){
	cutStopTimer(*t);
	float tmp = cutGetTimerValue(*t);
	cutDeleteTimer(*t);
	return tmp;
}