今回の内容 CUDA 付属のライブラリ cublas 行列ベクトル積, 行列行列積 cusperse 行列格納形式 cufft 余弦波の FFT curand モンテカルロ法による円周率計算 Thrust 913

(1)

GPU最適化ライブラリ

(2)

今回の内容

2015/07/22 GPGPU実践プログラミング 913



_{CUDA付属のライブラリ}



_cuBLAS



行列‐ベクトル積，行列‐行列積



_cuSPERSE



行列格納形式



_cuFFT



余弦波のFFT



_cuRAND



モンテカルロ法による円周率計算



_Thrust

(3)

ライブラリ



特定の処理を行う複数のプログラムを再利用可能な形で

まとめた集合体



動画像処理やファイル圧縮，数値計算のライブラリが有名



数値計算ライブラリは自作のプログラムよりも性能が高い

ため，処理速度の向上に貢献

(4)

数値計算ライブラリ



_{FFT（Fast Fourier Transform）}



_FFTW



線形代数演算（ベクトル処理，行列処理）



_{BLAS（Basic Linear Algebra Subprogram）}



_{BLASを利用した線形代数演算ライブラリ}



_LAPACK



_LINPACK



_ScaLAPACK



_{BLASやLAPACKのメーカー別実装}



_{MKL Intel Math Kernel Library}



_{ACML AMD Core Math Library}



_{IMSL International Mathematics and Statistics Library}

(5)

CUDA付属のライブラリ



_{NVIDIA製GPUに最適化}



_cuBLAS



_{BLASのGPU向け実装（密行列）}



_cuSPARSE



_{BLASのGPU向け実装（疎行列）}



_cuFFT



高速フーリエ変換



_cuRAND



乱数生成

(6)

cuBLAS



_{BLASのGPU向け実装}



_Level 1



ベクトルに対する計算



_Level 2



行列－ベクトルの計算



_Level 3



行列－行列の計算



_{BLASはFORTRAN用ライブラリとして開発}



_{C言語から使用する場合には注意が必要}



配列の添字は1から



メモリの配置がC言語と異なる

(7)

cuBLAS



関数群はドキュメントを参照



_{http://docs.nvidia.com/cuda/cublas/}



グリッド，ブロック，スレッドの設定は自動



処理内容に対応した関数を呼び出す



_{cublas<>axpy()}



<>は用いる型に応じて選択



_{S 単精度実数}



_{D 倍精度実数}



_{C 単精度複素数型}



_{Z 倍精度複素数型}

(8)

cuBLAS



計算内容に対応した関数を呼ぶ事で計算を実行



_{cublas<>axpy()}

ベクトル和



_{cublas<>gemv()}

行列‐ベクトル積



_{cublas<>gemm()}

行列‐行列積



コンパイルの際はオプションとして

_–lcublas

を追加

y[] = x[] + y[]

y[] = OP(A[][])x[] + y[]

C[][] = OP(A[][])OP(B[][]) + C[][]

(9)

#include<stdio.h> #include<math.h> #include<cublas_v2.h> #define M 256 #define N 256 int main(){ int i,j; float A[M*N], B[N], C[M]; float *dev_A, *dev_B, *dev_C; float alpha, beta; for(j=0;j<N;j++){ for(i=0;i<N;i++){ A[i+M*j] = ((float) i)/256.0; } } for(i=0;i<M;i++){ C[i] = 0.0f; } for(i=0;i<N;i++){ B[i] = 1.0f; } cudaMalloc ((void**)&dev_A, N*N*sizeof(float)); cudaMalloc ((void**)&dev_B, N*sizeof(float)); cudaMalloc ((void**)&dev_C, M*sizeof(float)); alpha = 1.0f; beta = 0.0f;

cuBLASによる行列－ベクトル積

cublas_gemv.cu

(10)

cublasHandle_t handle; cublasCreate(&handle); cublasSetMatrix (M, N, sizeof(float), A, M, dev_A, M); cublasSetVector (N, sizeof(float), B, 1, dev_B, 1); cublasSetVector (M, sizeof(float), C, 1, dev_C, 1); cublasSgemv (handle, CUBLAS_OP_N, M, N, &alpha, dev_A, M, dev_B, 1, &beta, dev_A, 1); cublasGetVector (M, sizeof(float), dev_C, 1, C, 1); for(i=0;i<M;i++){ printf("C[%3d] = %f ¥n",j, C[j]); } cudaFree(dev_A); cudaFree(dev_B); cudaFree(dev_C); cublasDestroy(handle); return 0; }

cuBLASによる行列－ベクトル積

GPGPU実践プログラミング 921 2015/07/22

cublas_gemv.cu

(11)

cuBLASによる行列－ベクトル積



ヘッダファイルのインクルード



_{#include<cublas_v2.h>}



_{CUBLASのハンドルの作成}



_{cublasHandle_t handle;}



_{cublasCreate(&handle);}



ハンドル（行列やベクトルの情報を管理）の生成



返値はcublasStatus_t型



_{GPU上のメモリ確保}



_cudaMalloc

(12)

cuBLASによる行列－ベクトル積



データ転送



_{cublasSetMatrix(M,N,sizeof(float),A,M,dev_A,M);}



_{cublasSetVector(N,sizeof(float),B,1,dev_B,1);}



_{cublasSetVector(M,sizeof(float),C,1,dev_C,1);}



関数内部でcudaMemcpyが呼ばれる

アクセスの

ストライド

(13)

cuBLASによる行列－ベクトル積



演算を行う関数の呼び出し



_{cublasSgemv(handle,}

_{CUBLAS_OP_N}

_{, M, N, &alpha,}

dev_A, M, dev_B, 1, &beta, dev_C, 1);



_{CUBLAS_OP_N 行列Bに対する操作}



cublasOperation_t型



CUBLAS_OP_N OP(B)＝B

処理しない



_{CUBLAS_OP_T OP(B)＝B}

T

転置



_{CUBLAS_OP_C OP(B)＝B}

H

共役転置

C[M] = 

OP

(A[M][N])B[N] + C[M]

アクセスの

ストライド

(14)

cuBLASによる行列－ベクトル積



データの読み戻し



_{cublasGetVector(M,sizeof(float),dev_C,1,C,1);}



_{GPU上のメモリの解放}



_cudaFree



ハンドルの削除



_{cublasDestroy(handle);}

アクセスの

ストライド

(15)

cuBLASを使う際の注意点



行列を表現するときのメモリ配置



_A[i][j]



_{2次元配列でもメモリ上は1次元に配置}



_{i方向が先かj方向が先か}



_{C言語はj方向優先}



_{Fortranはi方向優先}

(16)

cuBLASを使う際の注意点



_C言語

におけるA[i][j]のメモリ上の配置



0/256,0/256,0/256･･･1/256,1/256,1/256･･･2/256,2/256,2/256 ･･･



_Fortran

におけるA(i,j)のメモリ上の配置 ←BLAS



0/256,1/256,2/256･･･0/256,1/256,2/256･･･0/256,1/256,2/256･･･



















256 /

2

256 /

2

256 /

2

256 /

1

256 /

1

256 /

1

256 /

0

256 /

0

256 /

0 j

i

(17)

cuBLASによる行列－行列積



_{cublasSgemm(cublasHandle_t handle,}

char

trans_a

, char

trans_b

,

int m, int n, int k,

float

alpha

, const float *A, int lda,

const float *B, int ldb,

float

beta

,

float *C, int ldc);



_{trans_a,trans_b}



CUBLAS_OP_N, CUBLAS_OP_T, CUBLAS_OP_C



_{lda, ldb, ldc 行数を指定（m, k, m）}



正方行列ならm,n,kのどれでも値が同じ

(18)

#include<stdio.h> #include<math.h> #include<cublas_v2.h> #define N 4096 int main(){ int i,j; float *A, *B, *C; float *dev_A, *dev_B, *dev_C; float alpha, beta; A=(float *)malloc(N*N*sizeof(float)); B=(float *)malloc(N*N*sizeof(float)); C=(float *)malloc(N*N*sizeof(float)); for(j=0;j<N;j++){ for(i=0;i<N;i++){ A[i+N*j] = (float)i; B[i+N*j] = (float)j; C[i+N*j] = 0.0f; } } cudaMalloc ((void**)&dev_A, N*N*sizeof(float)); cudaMalloc ((void**)&dev_B, N*N*sizeof(float)); cudaMalloc ((void**)&dev_C, N*N*sizeof(float)); alpha = 1.0f; beta = 0.0f;

cuBLASによる行列－行列積

cublas_gemm.cu

(19)

cublasHandle_t handle; cublasCreate(&handle); cublasSetMatrix(N, N, sizeof(float), A, N, dev_A, N); cublasSetMatrix(N, N, sizeof(float), B, N, dev_B, N); cublasSetMatrix(N, N, sizeof(float), C, N, dev_C, N); cublasSgemm(handle, CUBLAS_OP_N, CUBLAS_OP_N, N, N, N, &alpha, dev_A, N, dev_B, N, &beta, dev_C, N); cudaFree(dev_A); cudaFree(dev_B); cudaFree(dev_C); cublasDestroy(handle); return 0; }

cuBLASによる行列－行列積

cublas_gemm.cu

(20)

cuBLASによる行列－行列積の実効性能

2015/07/22 GPGPU実践プログラミング 931 2 2 2 2 2 2 2 2

行列サイズ

GFLOPS

Tesla M2050の理論演算性能

(21)

cuSPARSE



_{BLASのGPU向け実装（疎行列向け）}



_{Level 1, Level 2, Level 3が存在}



行列の格納形式の操作に関するHelper Function



疎行列



行列の要素の大半が0



非ゼロ要素がまばらに存在



非ゼロ要素のみを格納



そもそも疎行列を作るのが大変

疎行列

密行列

ゼロ要素

非ゼロ要素

(22)

疎行列格納形式



疎行列



行列の要素の大半が0



実際には8割程度の要素が0



全てを保持するのは効率が悪い



疎行列格納形式



非ゼロ要素だけを保持してメモリを有効に利用



_{COO形式(Coordinate Format)}



_{CSR形式(Compressed Sparse Row Format)}



_{CSC形式(Compressed Sparse Column Format)}



その他色々













9

0

8

0

7

6

0

5

4

0

3

0

2

0

1

(23)

疎行列格納形式（COO形式）



_{COO形式(Coordinate Format)}



非ゼロ要素の行番号，列番号，値を保持













9

0

8

0

7

6

0

5

4

0

3

0

2

0

1 Row = {0, 0, 1, 1, 2, 2, 2, 3, 3}

Column = {0, 2, 1, 3, 0, 2, 3, 1, 3}

Value = {1, 2, 3, 4, 5, 6, 7, 8, 9}

0

1

2

3 Column

0

1

2

3 Row

(24)

疎行列格納形式（CSR形式）



_{CSR形式(Compressed Sparse Row Format)}



_{COO形式のRowでは同じ数字が繰り返される}



_{Rowのデータも圧縮}













9

0

8

0

7

6

0

5

4

0

3

0

2

0

1 Value = {1, 2, 3, 4, 5, 6, 7, 8, 9}

Column = {0, 2, 1, 3, 0, 2, 3, 1, 3}

Row = {0, 2, 4, 7, 9}

0

1

2

3 Column

0

1

2

3 Row

(25)

疎行列格納形式（CSR形式）













9

0

8

0

7

6

0

5

4

0

3

0

2

0

1

0

1

2

3 Column

0

1

2

3 Row

Value = {1, 2, 3, 4, 5, 6, 7, 8, 9}

Column = {0, 2, 1, 3, 0, 2, 3, 1, 3}

Row = {0, 2, 4, 7, 9}

Column配列の0,1は0行目のデータ

i行目のデータの個数は Row[i+1]‐Row[i]

forループが簡単に書ける for(j=row[i]; j<row[i+1];j++)

(26)

疎行列格納形式（CSR形式）













9

0

8

0

7

6

0

5

4

0

3

0

2

0

1

0

1

2

3 Column

0

1

2

3 Row

Value = {1, 2, 3, 4, 5, 6, 7, 8, 9}

Column = {0, 2, 1, 3, 0, 2, 3, 1, 3}

Row = {0, 2, 4, 7, 9}

Column配列の2,3は1行目のデータ

i行目のデータの個数は Row[i+1]‐Row[i]

forループが簡単に書ける for(j=row[i]; j<row[i+1];j++)

(27)

疎行列格納形式（CSR形式）













9

0

8

0

7

6

0

5

4

0

3

0

2

0

1

0

1

2

3 Column

0

1

2

3 Row

Value = {1, 2, 3, 4, 5, 6, 7, 8, 9}

Column = {0, 2, 1, 3, 0, 2, 3, 1, 3}

Row = {0, 2, 4, 7, 9}

Column配列の4,5,6は2行目のデータ

i行目のデータの個数は Row[i+1]‐Row[i]

forループが簡単に書ける for(j=row[i]; j<row[i+1];j++)

(28)

疎行列格納形式（CSR形式）













9

0

8

0

7

6

0

5

4

0

3

0

2

0

1

0

1

2

3 Column

0

1

2

3 Row

Value = {1, 2, 3, 4, 5, 6, 7, 8, 9}

Column = {0, 2, 1, 3, 0, 2, 3, 1, 3}

Row = {0, 2, 4, 7, 9}

Column配列の7,8は3行目のデータ

i行目のデータの個数は Row[i+1]‐Row[i]

forループが簡単に書ける for(j=row[i]; j<row[i+1];j++)

(29)

cuSPARSEの行列－行列積



_{N×Nの正方行列同士のかけ算}



_{cuBLASの場合}



_[C] =



_OP

_([A])*

_OP

_([B]) +



_[C]



_{cublasSgemm(handle,}

_OPERATION

_,

_OPERATION

_,

N, N, N,



, [A], N, [B], N,



, [C], N);



_{cuSPERSEの場合}



_[C] =



_OP

_([A])*[B] +



_[C]



_{cusparseScsrmm(handle,}

_OPERATION

_{, N, N, N, NNZ,}



, descr[A], csrVal[A], csrRow[A], csrCol[A],

[B], N,



, [C], N);

行列サイズ

行数

行列転置の有無

非ゼロ

要素数

行列サイズ

(30)

cuFFT



_{FFTWと似たインターフェースを持つFFTライブラリ}



以下の各FFTをサポート



_1次元



_2次元



_3次元



複素数 → 複素数



実数

_{→ 複素数}



複素数 → 実数

(31)

cuFFTの利用手順

1. planの作成



_{FFTの次元と種類を指定}



_{1次元，2次元，3次元のいずれか}



複素数→複素数，実数→複素数，複素数→実数のいずれか

2. FFTの実行



_{planに沿ってFFTを実行}



順方向FFTか逆方向FFTかを実行時に指定

3. planの破棄

(32)

planの作成



_FFTの次元



_1次元

cufftPlan

1d

(cufftHandle *plan, int nx, cufftType type,int batch);



_2次元

cufftPlan

2d

(cufftHandle *plan, int nx,

int ny

, cufftType type);



_3次元

cufftPlan

3d

(cufftHandle *plan, int nx, int ny,

int nz

, cufftType

type);



_{FFTの種類（type）}



複素数 → 複素数

_{CUFFT_C2C}



実数

_{→ 複素数}

_{CUFFT_R2C}



複素数 → 実数

_{CUFFT_C2R}

cufftHandle型

ポインタ変数

FFTを行うデータの

x,y,z方向サイズ

1次元FFTを

何回行うか

(33)

FFTの実行



複素数 → 複素数

cufftExecC2C(cufftHandle plan, cufftComplex *idata,

cufftComplex *odata, int direction);



実数

_{→ 複素数}

cufftExecR2C(cufftHandle plan, cufftReal *idata,

cufftComplex *odata, int direction);



複素数 → 実数

cufftExecC2R(cufftHandle plan, cufftComplex *idata,

cufftReal *odata, int direction);



_FFTの方向



_{int direction}



_{CUFFT_FORWARD}

順FFT



_{CUFTT_INVERSE}

逆FFT

(34)

FFTの実行，planの破棄



倍精度利用する場合



_{cufftDoubleComplex, cufftDoubleRealを指定}

cufftExecZ2D(cufftHandle plan, cufftDoubleComplex *idata,

cufftDoubleReal *odata, int direction);



_{planの破棄}

(35)

#include<stdio.h> #include<math.h> #include<cufft.h> #define N 256 #define MAX_WAVE 5 int main(){ int i,wavenum; float real[N],imag[N]; cufftHandle plan; cufftComplex host[N]; cufftComplex *dev; cufftComplex re_host[N]; for(i=0;i<N;i++){ real[i] = 0.0; imag[i] = 0.0; } for(wavenum=1;wavenum<=MAX_WAVE; wavenum++){ for(i=0;i<N;i++){ real[i] += cos(2.0*M_PI*(float)i/N *(float)wavenum); } } for(i=0;i<N;i++){ printf("%f,%f¥n",(float)i,real[i]); } for(i=0;i<N;i++){ host[i] = make_cuComplex( real[i],imag[i]); }

cuFFTによる余弦波のFFT

cufft.cu

(36)

cudaMalloc((void**)&dev, N*sizeof(cufftComplex)); cudaMemcpy(dev, host, N*sizeof(cufftComplex), cudaMemcpyHostToDevice); cufftPlan1d(&plan, N, CUFFT_C2C,1); cufftExecC2C(plan, dev, dev, CUFFT_FORWARD); cudaMemcpy(re_host, dev, N*sizeof(cufftComplex), cudaMemcpyDeviceToHost); for(i=0;i<N/2;i++){ printf("dev[%3d] = %f %f ¥n", i,cuCrealf(re_host[i]), cuCimagf(re_host[i])); } cufftDestroy(plan); cudaFree(dev); return 0; }

cuFFTによる余弦波のFFT

cufft.cu

(37)

cuFFTによる余弦波のFFT



ヘッダファイルのインクルード



_{#include<cufft.h>}



二つの実数（実部と虚部）から複素数を作成



_{host[i] = make_cuComplex(real[i],imag[i]);}



_{GPU上のメモリ確保}



_cudaMalloc



_{1次元FFTのプランを作成}



_{cufftPlan1d(&plan, N, CUFFT_C2C, 1);}

N個のデータを使って

1次元FFTを1回行う

(38)

cuFFTによる余弦波のFFT



複素数から複素数へのFFTを実行



_{cufftExecC2C(plan, dev, dev,CUFFT_FORWARD);}



結果の表示（実部と虚部の取り出し）



_{cuCrealf(re_host[i])}

複素数の実部を取り出す



_{cuCimagf(re_host[i])}

複素数の虚部を取り出す



プランの破棄とGPU上のメモリの解放



_{cufftDestroy(plan);}



_cudaFree



コンパイルの際はオプションとして

_–lcufft

を追加

(39)

cuFFTによる余弦波のFFT



入力信号



波数1~5の余弦波の重ね合わせ

-2 -1 0 1 2 3 4 5 0 32 64 96 128 160 192 224 256

(40)

cuFFTによる余弦波のFFT



スペクトラムの実部

0 20 40 60 80 100 120 140 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16

入力した波の波数

波数

スペクトラム

(41)

cuRAND



乱数（列）を生成するライブラリ



既知の数列から次の数列の値が予測できない



一様乱数，正規乱数など



真の乱数列の生成は困難



乱数列と同じような性質を持った数列（疑似乱数）を生成



_{CPUから呼ばれ，GPUで乱数を生成}



_{GPUから呼ばれ，GPUで乱数を生成}

(42)

cuRAND



_{CPUから関数を呼ぶ場合}



ヘッダーファイルのインクルード



_{#include<curand.h>}



_{Generatorを宣言}



_{curandGenerator_t 変数名}

(43)

cuRAND



_{Generatorの生成}

curandStatus_t curandCreateGenerator

(curandGenerator_t* generator, curandRngType_t rng_type )



_{curandRngType_t（乱数の発生方法）の指定}



_{CURAND_RNG_PSEUDO_XORWOW}

_XORWOW



_{CURAND_RNG_PSEUDO_MRG32K3A}

_MRG32k3a



_{CURAND_RNG_PSEUDO_MTGP32}

_{メルセンヌ・ツイスタ}



_{CURAND_RNG_QUASI_SOBOL32}

_Sobol32



_{CURAND_RNG_QUASI_SCRAMBLED_SOBOL64}

_{Scrambled Sobol64}

(44)

cuRAND



乱数の種の設定

curandStatus_t curandSetPseudoRandomGeneratorSeed

(curandGenerator_t generator, unsigned long long seed )



乱数の生成



一様乱数を生成

curandStatus_t curandGenerate

Uniform

(curandGenerator_t generator, float* outputPtr,

size_t num )

(45)

#include <stdio.h> #include <stdlib.h> #include <math.h> #include <curand.h> #define N 1024 int main(){ int i; float *value, *value_d; curandGenerator_t gen; cudaMalloc((void**)&value_d, N*sizeof(float)); curandCreateGenerator (&gen, CURAND_RNG_PSEUDO_DEFAULT); curandGenerateUniform (gen, value_d, N); value = (float *)malloc(N*sizeof(float)); cudaMemcpy(value, value_d, N*sizeof(float), cudaMemcpyDeviceToHost); for(i=0;i<N;i++){ printf("%d¥n",value[i]); } curandDestoryGenerator(gen); cudaFree(value_d); free(value); return 0; }

cuRANDサンプル

curand.cu

(46)

cuRAND



_{GPUから乱数を生成するデバイス関数を呼ぶ場合}



ヘッダファイルのインクルード



_{#include<curand_kernel.h>}



乱数生成の初期化



_{device void curand_init(unsigned long long seed, unsigned}

long long subsequence, unsigned long long offset,

curandStateMRG32k3a_t* state )



乱数の生成



_{device float curand_uniform(curandStateMtgp32_t* state )}

GPGPU実践プログラミング

(47)

#include<stdio.h> #include<stdlib.h> #include<curand_kernel.h> #define N (1024) __global__ void random(float *value){ int i; unsigned int seed; curandState stat; seed = 1; curand_init(seed, 0, 0, &stat); for(i=0;i<N;i++){ value[i] = curand_uniform(&stat); } } int main(){ int i; float *value, *value_d; cudaMalloc((void**)&value_d, N*sizeof(float)); random<<<1,1>>>(value_d); value = (float *)malloc (N*sizeof(float)); cudaMemcpy(value, value_d, N*sizeof(float), cudaMemcpyDeviceToHost); for(i=0;i<N;i++){ printf("%f ¥n",value[i]); } cudaFree(value_d); free(value); return 0; }

cuRANDサンプル

curand_kernel.cu

(48)

Thrust



_{CUDA用の並列アルゴリズムライブラリ}



_{C++の標準テンプレートライブラリ(STL)とよく似た高水}

準なインタフェースを保有



_{CUDA4.0からCUDA本体に吸収}



_{https://developer.nvidia.com/thrust}

(49)

#include<iostream> #include<thrust/host_vector.h> #include<thrust/device_vector.h> #include<thrust/copy.h> #include<thrust/sort.h> int main(void) { thrust::host_vector < float > host_vec(3); thrust::device_vector < float > device_vec(3); host_vec[0] = 1.1; host_vec[1] = 3.3; host_vec[2] = 2.2; thrust::copy(host_vec.begin(), host_vec.end(), device_vec.begin()); thrust::sort(device_vec.begin(), device_vec.end()); thrust::copy(device_vec.begin(), device_vec.end(), host_vec.begin()); std::cout << host_vec[0] << std::endl; std::cout << host_vec[1] << std::endl; std::cout << host_vec[2] << std::endl; return 0; }

Thrustサンプル

thrust_sort.cu

(50)

Thrustサンプル



ヘッダのインクルード



_{#include<thrust/host_vector.h>}



_{#include<thrust/device_vector.h>}



_{#include<thrust/copy.h>}



メモリ確保



_{host_vector}

_{CPUにメモリを確保}



_{device_vector GPUにメモリを確保}



_{CPU側の配列の初期化とGPUへのコピー}



_{thrust::copy(host_vec.begin(), host_vec.end(),}

device_vec.begin());

(51)

Thrustサンプル



データのソート



_{thrust::sort(device_vec.begin(),}

device_vec.end());



_{GPUから読み戻し}



_{thrust::copy(device_vec.begin(),}

device_vec.end(), host_vec.begin());

(52)

モンテカルロ法



乱数を用いて数値シミュレーションや数値計算を行う手

法の総称



得られる解が必ずしも正しいとは保証されていない



円周率の計算



正方形の中に乱数を使って点を打ち，円の内側に入っている

点を数えることで円の面積を求める



点の数が無限に多く，点の座標が重複しないなら面積（や円

周率）が正しく求められる

(53)

モンテカルロ法による円周率計算



円の面積



r

2 

正方形の面積

4r

2 

面積比

a=



r

2 /4r

2 

円周率



=4a



円と正方形の面積比を点の数で近似

2r

r

(54)

#include<stdio.h> #include<stdlib.h> #include<math.h> #include<curand.h> #define N 1024 int main(){ int i,inside; float *x, *y, *x_d, *y_d; float pi; curandGenerator_t gen; cudaMalloc((void**)&x_d, N*sizeof(float)); cudaMalloc((void**)&y_d, N*sizeof(float)); curandCreateGenerator(&gen, CURAND_RNG_PSEUDO_DEFAULT); curandGenerateUniform(gen, x_d, N); //x座標 curandGenerateUniform(gen, y_d, N); //y座標 x = (float *)malloc(N*sizeof(float)); y = (float *)malloc(N*sizeof(float)); cudaMemcpy(x, x_d, N*sizeof(float), cudaMemcpyDeviceToHost); cudaMemcpy(y, y_d, N*sizeof(float), cudaMemcpyDeviceToHost); inside=0; for(i=0;i<N;i++){ if( (x[i]*x[i] + y[i]*y[i]) <=1.0f ) inside++; } pi = 4.0f*(float)inside/N; printf("%f¥n",pi); curandDestroyGenerator(gen); cudaFree(x_d); cudaFree(y_d); free(x); free(y); return 0; }

モンテカルロ法による円周率計算

montecarlo.cu

(55)

モンテカルロ法によって得られた円周率

点の数

円周率

相対誤差

2

1

_2.000000

_0.36338

2

_2.000000

_0.36338

2

3

_2.500000

_0.204225

2

4

_3.000000

_0.04507

2

5

_3.125000

_0.005282

2

6

_3.062500

_0.025176

2

7

_3.031250

_0.035123

2

8

_2.984375

_0.050044

2

9

_3.164062

_0.007152

2

10

_3.171875

_0.009639

点の数

円周率

相対誤差

2

11

_3.160156

_0.005909

2

12

_3.171875

_0.009639

2

13

_3.133301

_0.002639

2

14

_3.157227

_0.004976

2

15

_3.14917

_0.002412

2

16

_3.152649

_0.003519

2

17

_3.152374

_0.003432

2

18

_3.147675

_0.001936

2

19

_3.142609

_0.000323

2

20

_3.142635

_0.000332

(56)

モンテカルロ法によって得られた円周率

サンプル数

N

相対誤差

O(N

-0.5

₎

(57)

その他

GPU向けライブラリ



_MAGMA



_{http://icl.cs.utk.edu/magma/}



_{NVIDIA GPU向け（CUDA）の線形代数ライブラリ}



_{CPUも同時に使うハイブリッド型のライブラリであるため，}

今回の内容 CUDA 付属のライブラリ cublas 行列 ベクトル積, 行列 行列積 cusperse 行列格納形式 cufft 余弦波の FFT curand モンテカルロ法による円周率計算 Thrust 913

GPU最適化ライブラリ