ノート/CUDA/数値積分

サイトトップ

ノート
マイニング
├マイニング
├バスケット解析
├バスケット解析をＲで
├図書貸出をＲで
├モデルクラスタリング
└
テキストマイニング
├テキストマイニングTM
├TMとシソーラス
├PubMedをTM
├TMとMeSH
├TMとNLTK
├テキストマイニングとtagger
├医薬品添付文書ＤＢ
├論文の処理1
├剽窃１
├Bing検索
└
分子進化学
├分子進化学
└
Pythonと論文アクセス
├Pythonを使ってみる
├PythonでPubMed
├Pythonで...続き
└
CUDA
├GPU
├ProgrammingGuide
├Selandメモ
├数値積分
└
ACS、PHPからLDAP
├ACS導入
├新規ホスト移行
├再度やり直し
├OpenLDAP導入
└

Linuxでビデオ
├Linuxでビデオデータを作る
├AVIフォーマット
├x264ソースを見る
├Macでフレーム抜き
└

MPIと並列化
├並列処理の実用化
├Fedora10でMPICH2
├EPMRを見る
├EPMRを見る２
├EPMRを見る〜Ｃ＋＋２
├EPMRを見る〜Ｃ＋＋２
├EPMRを見る〜Ｃ＋＋３
├EPMRを見る〜Ｃ＋＋４
├CASP7と日本人
└
Macメモ
├Script

favicon物語

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

メニュー編集

最新の15件

CUDA
訪問者数　4524　　　　　　最終更新　2010-11-05 (金) 12:41:54

数値積分プログラムの並列実行その２（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桁の精度を得ている。

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

GPU所要時間を表にすると

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

1. ブロック＝スレッド＝１の時に9591ミリ秒に対して、最速はブロック256、スレッド64、またはブロック64、スレッド256の2.64ミリ秒である。これは速度比で3633倍であり、プロセス数では256×64＝16,000である。
2. ストリーミングプロセッサ数は448と思われるので、上記の速度比3633倍は、およそ448×8＝3584に相当する。
3. 並列化による加速は、1×1から32×32（＝1024）程度まではスレッド×ブロック数に比例している。32×64や64×32では加速率が比例に比べてやや鈍っている。
4. 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。

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