ノート/CUDA/数値積分

サイトトップ

ノート
データ処理
├マイニング
├バスケット解析
├バスケット解析をＲで
├図書貸出をＲで
├モデルクラスタリング
└
テキストマイニング
├テキストマイニングTM
├TMとシソーラス
├PubMedをTM
├TMとMeSH
├TMとNLTK
├テキストマイニングとtagger
├医薬品添付文書ＤＢ
├論文の処理1
├剽窃１
├Bing検索
└
分子進化学
├分子進化学
└
Pythonと論文アクセス
├Pythonを使ってみる
├PythonでPubMed
├Pythonで...続き
└
Macメモ
├Script

&ref(): File not found: "favicon_story.txt" at page "ノート/ノート";

自信ないので
├Winのバッチ
│├スクリプト入門
│├バッチファイル
│├SETコマンド
│├自分のメモ
│└
└

メニュー編集

数値積分プログラムの並列実行その２（2010-11-04（木）) †

下のプログラムに(1)バグがあったのと、(2)CUDA2.0にして倍精度浮動小数doubleを使ってみたのと、(3)もう少し時間データを取ってみたかったのと、で再実験した。

プログラムは以下の通り

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

//#define FLOAT float
#define FLOAT double

#define BLOCK 64    /* Number of blocks */
#define THREAD 128  /* Number of threads per block */

#define STRIPS 65536*256  /* Number of total strips */
//#define STRIPS 65536  /* Number of total strips */

void Host(FLOAT *Result);
__global__ void Kernel(FLOAT *Result);

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

FLOAT h_result[THREAD*BLOCK];

int main(int argc, char** argv){
        unsigned int timer;
        FLOAT pi;

        /*CUT_DEVICE_INIT();*/
        CUT_DEVICE_INIT(argc, argv);

        FLOAT *d_result;
        cudaMalloc((void**) &d_result, sizeof(FLOAT)*THREAD*BLOCK);
        cudaMemset(d_result, 0, sizeof(FLOAT)*THREAD*BLOCK);

        SetTimer(&timer);
        cudaMemcpy(d_result, h_result, sizeof(FLOAT)*THREAD*BLOCK, cudaMemcpyHostToDevice);

        fprintf(stderr, "THREAD=%d, BLOCK=%d\n", THREAD, BLOCK);
        dim3 grid(BLOCK, 1, 1);
        dim3 threads(THREAD, 1, 1);
        Kernel<<< grid, threads >>>(d_result);

        cudaMemcpy(h_result, d_result, sizeof(FLOAT)*THREAD*BLOCK, cudaMemcpyDeviceToHost);

        printf("計算時間 =%f(ms)\n", EndTimer(&timer));
        pi=0.0;
        for (int k=0; k<THREAD*BLOCK; k++) {
                pi += h_result[k];
                //printf("Result[%d]=%15.12f\n", k, h_result[k]);
        }
        printf("GPU計算結果 =%15.12f\n",  (pi - (((FLOAT)0.5)/((FLOAT)STRIPS))) * (FLOAT)4.0);

        cudaFree(d_result);

        SetTimer(&timer);
        Host(h_result);
        printf("ホストの計算時間 =%f(ms)\n", EndTimer(&timer));

        pi=0.0;
        for (int k=0; k<THREAD*BLOCK; k++) {
                pi += h_result[k];
                //printf("Result[%d]=%15.12f\n", k, h_result[k]);
        }
        printf("host計算結果 =%15.12f\n",  (pi - (((FLOAT)0.5)/((FLOAT)STRIPS))) * (FLOAT)4.0);
        //getch();
}

__global__ void Kernel(FLOAT *Result)
{
        //GPUでの処理
        int ix=blockIdx.x*blockDim.x + threadIdx.x;
        FLOAT tmp=0.0;

        for(int k=0; k<(STRIPS/BLOCK/THREAD); k++){
                FLOAT x = ((FLOAT) ((ix*(STRIPS/BLOCK/THREAD))+k)) / ((FLOAT)STRIPS);
                tmp+=((FLOAT)sqrt(1.0-x*x)) / ((FLOAT)STRIPS);
        }

        Result[ix]=tmp;
}

void Host(FLOAT *Result){
        //ホスト側での計算
        int k, ix;
        FLOAT x, tmp;
        for(ix=0; ix<(THREAD*BLOCK); ix++){
                tmp = 0.0;
                for(k=0; k<(STRIPS/BLOCK/THREAD); k++){
                        x = ( (((FLOAT)ix)*((FLOAT)(STRIPS/BLOCK/THREAD))) + (FLOAT)k) / ((FLOAT)STRIPS);
                        tmp=tmp+(((FLOAT)sqrt(1-(x*x))) / ((FLOAT)STRIPS) );
                        //printf("xi=%15.12f\n", ((FLOAT)sqrt(1.0-x*x)) / ((FLOAT)STRIPS)) ;
                } 

               Result[ix]=tmp; 

        }
}

void SetTimer(unsigned int *t){
        cutCreateTimer(t);
        cutStartTimer(*t);
}

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

改変について若干の説明：
それぞれのTHREADに割当てた計算が、きちんと両端を含む（k=0からk<(STRIPS/BLOCK/THREAD) ）ようにすると、面積全体の両端をも含むが、台形公式からは、両端を2回含むことになってしまうので、両端（x=0とx=1）だけは1回分引かなければならない。領域分割をしない場合は、両端を含まない2回ずつの積上げ計算（2で割らない）をして最後に両端の半分を引いていたのだが、そのプログラムを援用してスレッドへ領域分割をしたとき、両端を含めるように直すのを忘れていた。

現在はこれらの点は正しく修正され、それらしい値が計算されている。たとえば

Using device 0: Tesla S2050
THREAD=8, BLOCK=8192
計算時間 =12.614000(ms)
GPU計算結果 = 3.141592653573
ホストの計算時間 =479.681000(ms)
host計算結果 = 3.141592653573

また、CUDA 1.3（倍精度が使えない）でなく、2.0を使うようにし（Makefile中で2.0のみに限定）、結果も上記のように、台形ストリップ数65536*256の状態で、10進で12桁の精度を得ている。

いちおう、正しいらしい計算ができるようになったところで、改めて所要時間を比較してみる。

BLOCK	THREAD	GPU	HOST
1	1	9591.1	479.1
1	8	1199	480.1
1	16	599.53	479.7
1	32	299.81	479.27
1	64	149.94	479.19
1	128	76.39	479.28
1	256	40.9	479.71
1	512	29.6	479.42
1	1024	?	479.6
2	8	599.52	479.72
2	16	299.8	479.18
2	32	149.94	479.17
2	64	75.01	479.15
2	128	38.23	479.16
2	256	20.42	479.3
2	512	14.85	480.08
4	8	288.8	479.18
4	16	149.94	479.27
4	32	75.01	480.52
4	64	37.55	479.23
4	128	19.16	501.82
4	256	10.28	479.78
4	512	7.47	501.65
4	1024	?
8	8	149.95	479.73
8	16	75.03	479.14
8	32	37.55	479.13
8	64	18.81	479.21
8	128	9.63	479.11
8	256	5.21	502.83
8	512	3.81	502.33
16	8	75.01	479.09
16	16	37.55	479.09
16	32	18.82	479.08
16	64	9.6	479.08
16	128	5.28	479.09
16	256	3.97	481.76
16	512	3.86	479.6
16	1024	?
32	8	37.54	501.85
32	16	18.82	479.18
32	32	9.45	479.56
32	64	5.04	479.17
32	128	3.19	479.08
32	256	3.16	479.15
32	512	3.04	479.26
64	8	19.39	501.55
64	16	9.75	502.02
64	32	4.93	479.21
64	64	2.95	479.15
64	128	2.84	480.16
64	256	2.64	479.28
64	512	2.77	487.31
128	8	19.8	479.13
128	16	9.86	479.34
128	32	4.99	479.12
128	64	3.04	479.17
128	128	2.69	502.22
128	256	2.75	479.52
128	512	3.16	502.64
256	8	15.04	501.98
256	16	7.58	479.14
256	32	3.89	501.77
256	64	2.64	479.65
256	128	2.67	479.49
256	256	2.99	502.15
256	512	3.87	480.73
256	1024	?
512	8	12.8	479.15
512	16	6.53	501.87
512	32	3.44	501.54
512	64	2.7	480
512	128	2.96	490.95
512	256	3.65	503
512	512	4.94	504.52
1024	8	12.74	479.2
1024	16	6.56	479.33
1024	32	3.61	502.26
1024	64	2.98	502.77
1024	128	3.68	485.68
1024	256	4.88	504.22
1024	512	7.26	484.14
2048	8	12.32	497.28
2048	16	6.75	485.9
2048	32	3.84	479.84
2048	64	3.69	480.34
2048	128	4.85	481.62
2048	256	7.16	484.59
2048	512	13.02	502.89
4096	8	12.3	501.9
4096	16	6.76	502.58
4096	32	4.55	480.43
4096	64	4.87	504.23
4096	128	7.13	506.76
8192	4	24.01	501.95
8192	8	12.61	479.68
8192	16	7.48	480.39
8192	32	5.75	504.34
8192	64	7.21	507.25
8192	128	11.52	479.18

GPU所要時間を表にすると

BLOCK\THREAD	1	8	16	32	64	128	256	512
1	9591.1	1199	599.53	299.81	149.94	76.39	40.9	29.6
2		599.52	299.8	149.94	75.01	38.23	20.42	14.85
4		288.8	149.94	75.01	37.55	19.16	10.28	7.47
8		149.95	75.03	37.55	18.81	9.63	5.21	3.81
16		75.01	37.55	18.82	9.6	5.28	3.97	3.86
32		37.54	18.82	9.45	5.04	3.19	3.16	3.04
64		19.39	9.75	4.93	2.95	2.84	2.64	2.77
128		19.8	9.86	4.99	3.04	2.69	2.75	3.16
256		15.04	7.58	3.89	2.64	2.67	2.99	3.87
512		12.8	6.53	3.44	2.7	2.96	3.65	4.94
1024		12.74	6.56	3.61	2.98	3.68	4.88	7.26
2048		12.32	6.75	3.84	3.69	4.85	7.16	13.02
4096		12.3	6.76	4.55	4.87	7.13
8192		12.61	7.48	5.75	7.21	11.52

上記の表で言えることをまとめると、

ブロック＝スレッド＝１の時に9591ミリ秒に対して、最速はブロック256、スレッド64、またはブロック64、スレッド256の2.64ミリ秒である。これは速度比で3633倍であり、プロセス数では256×64＝16,000である。
ストリーミングプロセッサ数は448と思われるので、上記の速度比3633倍は、およそ448×8＝3584に相当する。
並列化による加速は、1×1から32×32（＝1024）程度まではスレッド×ブロック数に比例している。32×64や64×32では加速率が比例に比べてやや鈍っている。
Xeonホスト(2.1GHz)の実行時間480ミリ秒に比べて、単一スレッド時は9591/480＝19.98と約20倍遅い。これをマルチスレッドで実行すると最速時で480/2.64＝182倍速くなるので、単純計算だと24時間かかる計算が8分で終わる。GPUが４台搭載されているので、それを全て使うと、単純計算では8/4＝2分になる。もし１ヶ月＝30日かかる計算が2分×30＝１時間で終わるのだとすれば、これは朗報だろう。

↑

数値積分プログラムの並列実行（2009-02-10 (火) 16:16:24） †

＞＞CUDAバージョンアップ（2.0から3.1）に伴い、一部ソースファイル修正（2010-08-23）

挑戦するのは、（たとえば）台形則で簡単な関数を数値積分する例
　　関数は　sqrt( 1 - x*x )、つまり原点を中心とする半径1の円（の上半分）
　　積分区間は、0から1まで。つまり1/4円になる。
　　だから、積分結果（＝面積）を4倍すると、半径1の円の面積＝πになるはず。
さし当って2つの問題に興味：
　1) 並列度
　　　元のプログラムの並列度はかなり高くできる（台形ごとの計算など）
　　　これをGPUのコア（CUDAだとスレッド・ブロック・グリッドの階層＋ワープ）にどう割り当てるか
　　　共有メモリをどう使うか、それによってどういう速度低下の影響があるか
　2) 単精度浮動小数の問題
　　　単純に確認を要するだけだが、気になる。

↑

単精度浮動小数問題 †

ステップ1）普通のCPU（Pentium）で単精度と倍精度を比較してみよう。
　Cygwinを利用。プログラムは下記（左端は行番号）。実は結構細かい工夫が要る。

001  #include <stdio.h>
002  #include <math.h>
003
004  //#define STRIPS 65536*32565
005  #define STRIPS 65536*1
006  //#define FLOAT double
007  #define FLOAT float
008
009  main() {
010    int i;
011    FLOAT xi, sum, pi;
012    for (i=1; i<STRIPS; i++) {
013      xi = ((FLOAT) i) / ((FLOAT) STRIPS);
014      sum = sum + ( ((FLOAT) sqrt(1-(xi*xi))) / ((FLOAT)STRIPS) );
015      if ((i % 4096)==0) printf("i=%d, xi=%f, sum=%f\n", i, xi,  sum);
016    }
017    pi = (sum + (((FLOAT)0.5) / ((FLOAT)STRIPS))) * (FLOAT)4.0;
018    printf("%15.12f\n", pi);
019  }

行006, 007は、文字列FLOATをfloatに置き換えるかdoubleに置き換えるかの切替えをしている。
行004, 005で定義しているSTRIPSは、台形の本数（区間0～1を何等分するか）の定数

積分は台形則で計算するので、両端（x=0とx=1）を除いて、外の点の関数値（yの値）を求めて、足し合わせる。
　　これを、行012～016のforループで行う。
　　xiは、i番目の分割点のx座標。
　　　　このとき、iは整数のままforループで増やし、それぞれのiに対してxiを計算。
　　　　xiの差分（台形の高さ）を浮動小数で求めておいてそれを次々に足すやり方をすると、xi差分の誤差が蓄積するのでよくないだろう。
　　行014のsumの累積は、関数値（高さy）に幅（台形の高さ）を掛けてから足し合わせる。幅を掛ける＝STRIPSで割る。
　　　　代案として、関数値のまま足し合わせて最後に１回だけ幅を掛ける（STRIPSで割る）のが考えられるが、
　　　　STRIPSをここで書いたぐらいの値にすると、累積値sumが大きくなりすぎてfloatに収まらなくなるので、具合が悪い。
　　　　doubleだと、sumを関数値のまま（STRIPSで割らずに）足し合わせても大丈夫であることは確かめた。
　　行015は適当な間隔でループ内の動作の様子を見るためのデバッグ用。
行017で、台形則の両端（x=0とx=1）を（1/2にして）加える。
　　但し、x=1は関数値が0なので不要、x=0は関数値が1なので（1/2にして0.5で、それに幅（台形の高さ1/STRIPS）を掛け）、0.5/STRIPS
　　最後に、1/4円を全円にする為に４倍
行018の出力印刷は、一応doubleに対応するために、%15.2f（15桁文のスペースで、小数点より上が１桁のみ）とした。

計算結果は次の通り。
floatの場合

STRIPS   65536*1         65536*2         65536*4         65536*8         65536*16
pi      3.141592264175  3.141610622406  3.141520738602  3.141903400421  3.140739917755

doubleの場合

STRIPS   65536*1         65536*2         65536*4         65536*8         65536*16
pi      3.141592583496  3.141592628808  3.141592644828  3.141592650492  3.141592652495

floatの場合は、分割数が65536の時がπに近く、分割数を増やすと却ってπの値からはずれてゆく。それに対して、doubleの場合はπの真値に徐々に近づいてゆく。
ちなみに、floatの場合は仮数部23ビット～およそ10進7桁相当、doubleの場合は仮数部52ビット～およそ10進15桁相当、のはずである。
また、分割数がおよその精度を決めるとすると、65536*1なら10進5～6桁、65536*16なら10進7桁のはずである。
上記のdoubleの場合は、内部の変数の保持精度は15桁あるので、台形則（分割）による近似精度が表面に現れている。
　　65536*1では8桁目があやしいが、*4で9桁目、*8で10桁目、*16で10桁目が真値に近づくぐらいの精度である。
　　この結果では、上記の「分割数がおよその精度を決める」は正しくなく、分割数よりはるかに多い桁数の結果が得られている。
他方、floatの結果は、分割数を増やすと、常識的にはある一定の桁数の精度でとどまると予想されるが、分割数を増やすと却って真値から外れている。
この原因を考えてみるに、１つ可能性があるのは、計算途中の精度低下が、STRIPSが増えるとより低下する仕組があることで、具体的には誤差の蓄積が考えられる。
　　行014でsumを足し合わせているが、ここに累積するであろう誤差は、足し合わせる台形の数が増えれば増加する。これが可能性として考えられる。
　　しかし、累積誤差の影響は、一般論としては累積誤差は累積個数＝台形の個数に比例するだろうが、*4、*8、*16では台形の個数が１桁増える程度なのに、精度は7桁から4桁に3桁分も劣化している。十分な説明がつかない。
　　この他に、たとえば引き算に伴う桁落ちなどがあり得るが、このプログラムでは考えにくい。sumの全ての項は正である。

↑

並列化による実行時間短縮 †

上記のプログラムを並列に実行することを考える。
台形を左から同数のかたまりに分けて、そのかたまり毎に「スレッド」に割り当てる。この辺りのことはマニュアルがよく理解出来ない部分があるので、正しく割り当てられているか自信がない。
動かすスレッド数（＝並列度）を変化させて、処理時間を測定する。
プログラムは下記の通り。但し、台形則の両端（x=0とx=1）は、本来半分しか入れないが、並列計算時にはフルに入れて、全て終わってから余分の0.5/THREADSを減算する。

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

#define BLOCK 4    /* Number of blocks */
#define THREAD 32  /* Number of threads per block */
#define STRIPS 65536*256  /* Number of total strips */

void Host(float *Result);
__global__ void Kernel(float *Result);

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

float h_result[THREAD*BLOCK];

int main(int argc, char** argv){    // modified from main() due to CUDA version up to 3.1 (2010-08-23)
	unsigned int timer;

	CUT_DEVICE_INIT(argc, argv);   // modified from CUT_DEVICE_INIT() due to CUDA version up to 3.1 (2010-08-23)

	float *d_result;
	cudaMalloc((void**) &d_result, sizeof(float)*THREAD*BLOCK);
	cudaMemset(d_result, 0, sizeof(float)*THREAD*BLOCK);

	SetTimer(&timer);
	cudaMemcpy(d_result, h_result, sizeof(float)*THREAD*BLOCK, cudaMemcpyHostToDevice);

	dim3 grid(BLOCK, 1, 1);
	dim3 threads(THREAD, 1, 1);
	Kernel<<< grid, threads >>>(d_result);

	cudaMemcpy(h_result, d_result, sizeof(float)*THREAD*BLOCK, cudaMemcpyDeviceToHost);

	printf("計算時間 =%f(ms)", EndTimer(&timer));
	float pi=0.0;
	for (int k=0; k<THREAD*BLOCK; k++) {
		pi += h_result[k];
		//printf("Result[%d]=%f\n", k, h_result[k]);
	}
	printf("　計算結果 =%f\n", (pi - ((float)0.5/(float)STRIPS)) * ((float)4.0));

	cudaFree(d_result);

	SetTimer(&timer);
	Host(h_result);
	printf("ホストの計算時間 =%f(ms)", EndTimer(&timer));

	pi=0.0;
	for (int k=0; k<THREAD*BLOCK; k++) {
		pi += h_result[k];
		//printf("Result[%d]=%f\n", k, h_result[k]);
	}
	printf("　計算結果 =%f\n", (pi - ((float)0.5/(float)STRIPS)) * ((float)4.0));
	//getch();
}

__global__ void Kernel(float *Result)
{
	//GPUでの処理
	int ix=blockIdx.x*blockDim.x + threadIdx.x;
	float tmp=0.0;

	for(int k=0; k<(STRIPS/BLOCK/THREAD); k++){
		float x = ((float) (ix*(STRIPS/BLOCK/THREAD)+k)) / ((float)STRIPS);
		tmp+=((float)sqrt(1.0-x*x)) / ((float)STRIPS);
	}

	Result[ix]=tmp;
}

void Host(float *Result){
	//ホスト側での計算
	int ix;
	float tmp;
	
	for(ix=0; ix<THREAD*BLOCK; ix++){
		tmp = 0.0;
 		for(int k=0; k<(STRIPS/BLOCK/THREAD); k++){
			float x = ((float) (ix*(STRIPS/BLOCK/THREAD)+k)) /  (float)STRIPS);
			tmp+=((float)sqrt(1.0-x*x)) / ((float)STRIPS);
		}

		Result[ix]=tmp;

	}
}

void SetTimer(unsigned int *t){
	cutCreateTimer(t);
    cutStartTimer(*t);    
}

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

これにおける実行結果は次の通りになった。いずれも５回測定し、各回でGPUでの実行時間とCPUでの実行時間を測定している。（CPUでの実行時間は殆ど変わらない）単位はmS。

65536*256
#Blocks	#Threads	GPU time1	CPU time1	GPU time2	CPU time2	GPU time3	CPU time3	GPU time4	CPU time4	GPU time5	GPU time5	pi value	1/time/#Blk/#Thr
1	1	2313.6	547.8	2299.5	548	2310.8	548.2	2301	547.6	2310.4	548.3	3.456289	0.43
1	4	575.8	547.9	575.7	548	575.7	548.5	575.7	548.1	575.6	548.5		0.43
1	8	288	548	288.5	548.3	288.6	547.6	287.9	548.3	288.5	547.8		0.43
1	32	80.5	547.9	80.4	548.1	80.6	548.3	80.4	548.5	80.6	548.4		0.39
1	64	41.4	547.9	41.2	547.8	41.2	548.4	41.2	548.7	41.3	548.3		0.38
1	128	22.1	547.7	23.2	548.4	21.7	548.4	22.1	547.8	22.7	547.9		0.35
1	256	15.2	547.8	14.7	547.8	15	548.3	14.6	548.2	14.6	548	3.133834	0.26
4	1	575.5	547.9	575.7	548	575.7	548.4	575.6	548	575.5	548.2	3.179302	0.43
4	32	20.5	548.4	20	548.1	20	547.9	20	548	20	548	3.134089	0.38
4	128	6.24	547.8	5.71	548	5.71	547.8	5.71	547.9	6.11	547.9	3.133776	0.31
4	256	3.71	547.9	3.71	548.2	3.71	547.8	3.71	547.8	3.71	547.8	3.133775	0.26
16	128	1.51	547.8	1.51	548	1.48	548.1	1.51	548.2	1.48	548	3.133774	0.32
16	256	0.99	548.3	0.99	547.8	1	548	0.99	548.1	0.99	548.2	3.133765	0.25

GPUの測定結果は、ブロック数指定#Blocksとブロック内スレッド数指定#Threadsにほぼ反比例して実行時間が減少している。その意味では並列化が行われているように見える。
しかし、物理仕様は、「マルチプロセッサ」（通常の意味のプロセッサ）当り８ストリーム（～スレッド）分の処理回路があり、かつ全体で１４個（GT8800の場合）のプロセッサがあるので、112並列が限界のはずである。
但し、各マルチプロセッサ内のストリーム（スレッド）実行はリソースの制限から４サイクル毎に実行するとちょうど全部利用でき、そのためには１プロセッサ当り８×４＝３２スレッドを割り当てるのがよい。（ワープと呼ぶ。プログラム上は４ブロックを作るのがこれに当るらしい）
従って、８ストリーム×４（ワープ）×１４個＝４４８以上の並列向上は見込めないはずである。
しかし、上記の測定範囲でも最も遅いBlk=1, Thr=1の2310mSから、最も速いBlk=16, Thr=256の0.99mSまで、2300倍を超える速度ゲインがあったことになる。この理由がまだ説明できない。

最新の5件

数値積分プログラムの並列実行その２（2010-11-04（木）) †

数値積分プログラムの並列実行（2009-02-10 (火) 16:16:24） †

単精度浮動小数問題 †

並列化による実行時間短縮 †