階微分が 0 になるように袖領域の値を決定

GPU プログラム（ラプラシアン kernel ）

 袖領域の設定



グローバルメモリにデータがある箇所はグローバルメモリから読み込み



グローバルメモリにデータがない箇所は 2 階微

GPU プログラム（ラプラシアン kernel ）



共有メモリのデータを利用してラプラシアンを計算



if 分岐を排除

ulap[]

su[][]

ulap[ij] = (su[tx‐1][ty ]‐2.0*su[tx][ty]+su[tx+1][ty ])/dxdx +(su[tx ][ty‐1]‐2.0*su[tx][ty]+su[tx ][ty+1])/dydy;

GPU プログラムの評価

 共有メモリを使用した dif3.cu

 dif2.cu （共有メモリ不使用）より速いこともあれば遅いこともある

 Fermi 世代以降の GPU はキャッシュを搭載

 共有メモリを使っても速くならない

 共有メモリへ明示的にデータを移動

 余分な処理により負荷が増加

global void update(double u, double unew){

int i,j,ij;

i = blockIdx.x*blockDim.x + threadIdx.x;

j = blockIdx.y*blockDim.y + threadIdx.y;

ij = i+Nx*j;

u[ij] = unew[ij];

}

その他の処理の高速化

 値の更新

 unew のデータを u にコピーしているだけ

 cudaMemcpy で代用可能

その他の処理の高速化

 値の更新

 unew のデータを u にコピーしているだけ

 cudaMemcpy で代用可能

for(n=0;n<Nt;n++){

laplacian<<<Block, Thread>>>(dev_u,dev_ulap);

integrate<<<Block, Thread>>>(dev_u,dev_ulap,dev_unew);

//update<<<Block, Thread>>>(dev_u,dev_unew);

cudaMemcpy(dev_u, dev_new, Nbytes, cudaMemcpyDeviceToDevice);

}

return 0;

その他の処理の高速化

 初期値の計算

 三角関数がそれらしい値で求められればいい場合

 有効数字 4 桁程度

global void init(double u, double ulap, double *unew){

int i,j,ij;

double x,y;

i = blockIdx.x*blockDim.x + threadIdx.x;

j = blockIdx.y*blockDim.y + threadIdx.y;

ij = i+Nx*j;

x = (double)idx;y = (double)jdy;

u[ij]=sin(x)*sin(y);

unew[ij]=0.0;

ulap[ij]=0.0;

}

その他の処理の高速化

 初期値の計算

 三角関数がそれらしい値で求められればいい場合

 有効数字 4 桁程度

global void init(double u, double ulap, double *unew){

int i,j,ij;

double x,y;

i = blockIdx.x*blockDim.x + threadIdx.x;

j = blockIdx.y*blockDim.y + threadIdx.y;

ij = i+Nx*j;

x = (double)idx;y = (double)jdy;

u[ij]=sinf(x)*sinf(y);

unew[ij]=0.0;

ulap[ij]=0.0;

}

FastMath 関数

 三角関数など数学関数を高速に計算

 精度は落ちる

 手動で書き換え

 __sinf

 __cosf

 __powf など

 ‐use_fast_math オプション

 sin,cos を利用していても， ‐use‐fast‐math オプションを

付けることで sinf, cosf に置き換えられる

付録

共有メモリの典型的な使い方

 差分法

 あるスレッドが中心差分を計算するために配列の要素 i‐1, i, i+1 を参照

 配列要素は複数回参照される

 Fermi 世代以降はキャッシュが利用可能

1 1 2 2 1 1

dudx[i]

u[i]

＋＋＋＋

Δx 2

参照される回数

共有メモリの典型的な使い方

 境界で異なる処理を行うために if 分岐が必要

 キャッシュを搭載していてもグローバルメモリへのアクセスを伴う if 分岐は高負荷

 データの再利用と if 文の排除に共有メモリを利用

2 2 3 3 2 2

dudx[i]

u[i]

＋＋＋＋＋＋

Δx 2

参照される回数

共有メモリによる明示的なキャッシュ

 グローバルメモリから共有メモリにデータをキャッシュ

 共有メモリ上で境界条件を処理

 中心差分の計算から if を排除

dudx[i]

u[i]

共有メモリ

blockIdx.x=0 blockIdx.x=1

計算のために必要になる余分な領域（袖領域）

共有メモリによる明示的なキャッシュ

 グローバルメモリから共有メモリにデータをキャッシュ

 共有メモリ上で境界条件を処理

 中心差分の計算から if を排除

dudx[i]

u[i]

共有メモリ

blockIdx.x=0 blockIdx.x=1

何らかの方法で境界条件を反映

何らかの方法で

境界条件を反映

共有メモリによる明示的なキャッシュ

 グローバルメモリから共有メモリにデータをキャッシュ

 共有メモリ上で境界条件を処理

 中心差分の計算から if を排除

dudx[i]

u[i]

共有メモリ

＋

＋＋＋＋＋

blockIdx.x=0 blockIdx.x=1

全スレッドが同じ式

で中心差分を計算

#include<stdlib.h>

#include<math.h>/*‐lm

オプションが必要

#define Lx (2.0*M_PI)

#define Nx (1024*1024)

#define dx (Lx/(Nx‐1))

#define Nbytes (Nx*sizeof(double))

#define NT (256)

#define NB (Nx/NT) void init(double *u){

int i;

for(i=0; i<Nx; i++){

u[i] = sin(i*dx);

} }

__global__ void differentiate

(double *u, double *dudx){

int i = blockIdx.x*blockDim.x

+ threadIdx.x;

if(i==0)

dudx[i]=(‐3.0*u[i ] +4.0*u[i+1]

‐ u[i+2])/(2.0*dx);

if(0<i && i<Nx‐1)

dudx[i] = ( u[i+1]

‐u[i‐1])/(2.0*dx);

if(i==Nx‐1)

dudx[i]=( u[i‐2]

‐4.0*u[i‐1]

+3.0*u[i ])/(2.0*dx);

}

GPU プログラム（単純な実装）

differentiate1d.cu

int main(void){

double *host_u,*host_dudx;

double *u,*dudx;

host_u =(double *)malloc(Nbytes);

host_dudx=(double *)malloc(Nbytes);

cudaMalloc((void **)&u,Nbytes);

cudaMalloc((void **)&dudx,Nbytes);

init(host_u);

cudaMemcpy(u,host_u,Nbytes,

cudaMemcpyHostToDevice);

differentiate<<<NB, NT>>>(u,dudx);

cudaMemcpy(host_dudx, dudx, Nbytes, cudaMemcpyDeviceToHost);

free(host_u);

free(host_dudx);

cudaFree(u);

cudaFree(dudx);

return 0;

}

GPU プログラム（単純な実装）

differentiate1d.cu

__global__ void differentiate(double *u, double *dudx){

int i = blockIdx.x*blockDim.x + threadIdx.x;

__shared__ double su[1+NT+1];

int tx = threadIdx.x+1;

su[tx] = u[i];

__syncthreads();

if(blockIdx.x> 0 && threadIdx.x==0 ) su[tx‐1] = u[i‐1];

if(blockIdx.x< gridDim.x‐1 && threadIdx.x==blockDim.x‐1) su[tx+1] = u[i+1];

if(blockIdx.x==0 && threadIdx.x==0 )

su[tx‐1] = 3.0*su[tx]‐3.0*su[tx+1]+su[tx+2];

if(blockIdx.x==gridDim.x‐1 && threadIdx.x==blockDim.x‐1)

su[tx+1] = 3.0*su[tx]‐3.0*su[tx‐1]+su[tx‐2];

__syncthreads();

dudx[i] = ( su[tx+1]‐su[tx‐1])/(2.0*dx);

}

共有メモリを用いた書き換え

differentiate1d_shared.cu

共有メモリの宣言と代入

int i = blockIdx.x*blockDim.x + threadIdx.x;

__shared__ double su[1+NT+1]; //

右と左の袖領域を追加して宣言

int tx = threadIdx.x+1;

su[tx] = u[i];

__syncthreads();

u[i]

su[NT+2]

i= 0 1 2 3 4 5

threadIdx.x= 0 1 2 0 1 2

tx= 0 1 2 3 4 0 1 2 3 4 NT=3 NT=3

u[0] u[1] u[2] u[3] u[4] u[5]

袖領域の処理

if(blockIdx.x> 0 && threadIdx.x==0 ) su[tx‐1] = u[i‐1];

if(blockIdx.x< gridDim.x‐1 && threadIdx.x==blockDim.x‐1) su[tx+1] = u[i+1];

u[i]

su[NT+2]

i= 0 1 2 3 4 5

0 1 2=threadIdx.x

tx= 0 1 2 3 4 0 1 2 3 4 blockIdx.x=0 blockIdx.x=1

u[0] u[1] u[2] u[3] u[4] u[5]

袖領域の処理

if(blockIdx.x> 0 && threadIdx.x==0 ) su[tx‐1] = u[i‐1];

if(blockIdx.x< gridDim.x‐1 && threadIdx.x==blockDim.x‐1) su[tx+1] = u[i+1];

u[i]

su[NT+2]

i= 0 1 2 3 4 5

0=threadIdx.x

tx= 0 1 2 3 4 0 1 2 3 4 blockIdx.x=0 blockIdx.x=1

u[0] u[1] u[2] u[2] u[3] u[4] u[5]

袖領域の処理

if(blockIdx.x> 0 && threadIdx.x==0 ) su[tx‐1] = u[i‐1];

if(blockIdx.x< gridDim.x‐1 && threadIdx.x==blockDim.x‐1) su[tx+1] = u[i+1];

u[i]

su[NT+2]

i= 0 1 2 3 4 5

tx= 0 1 2 3 4 0 1 2 3 4

blockIdx.x=0 blockIdx.x=1

u[0] u[1] u[2] u[2] u[3] u[4] u[5]

threadIdx.x= 0 1 2

袖領域の処理

if(blockIdx.x> 0 && threadIdx.x==0 ) su[tx‐1] = u[i‐1];

if(blockIdx.x< gridDim.x‐1 && threadIdx.x==blockDim.x‐1) su[tx+1] = u[i+1];

u[i]

su[NT+2]

i= 0 1 2 3 4 5

tx= 0 1 2 3 4 0 1 2 3 4

blockIdx.x=0 blockIdx.x=1

u[0] u[1] u[2] u[3] u[2] u[3] u[4] u[5]

threadIdx.x=2

境界条件の処理

if(blockIdx.x==0 && threadIdx.x==0 )

su[tx‐1] = 3.0*su[tx]‐3.0*su[tx+1]+su[tx+2];

if(blockIdx.x==gridDim.x‐1 && threadIdx.x==blockDim.x‐1)

su[tx+1] = 3.0*su[tx]‐3.0*su[tx‐1]+su[tx‐2];

su[NT+2]

tx= 0 1 2 3 4 0 1 2 3 4

blockIdx.x=0 blockIdx.x=1

u[0] u[1] u[2] u[3] u[2] u[3] u[4] u[5]

2 1

3 u 3 u u

u

_

  

Δx u u

Δx

u u

u dx

du

2 2

4

3

₀ ₁ ₂ ₁ ₁





 





 

境界での差分式と中心差分式が一致するように u

₋₁

を決定

u[‐1]

境界条件の処理

if(blockIdx.x==0 && threadIdx.x==0 )

su[tx‐1] = 3.0*su[tx]‐3.0*su[tx+1]+su[tx+2];

if(blockIdx.x==gridDim.x‐1 && threadIdx.x==blockDim.x‐1)

su[tx+1] = 3.0*su[tx]‐3.0*su[tx‐1]+su[tx‐2];

su[NT+2]

tx= 0 1 2 3 4 0 1 2 3 4 blockIdx.x=0 blockIdx.x=1

u[0] u[1] u[2] u[3] u[2] u[3] u[4] u[5]

3 2

3

_



_



_



_N _N _N

u u u

u

Δx u u

Δx

u u

u dx

du

_N _N _N _N _N

2 2

4

3

₁ ₂ ₃ ₂





 



 

境界での差分式と中心差分式が一致するように u

を決定

u[‐1] u[6]

ドキュメント内 Microsoft PowerPoint - 先端GPGPUシミュレーション工学特論(web).pptx (ページ 73-97)

階微 分が 0 になるように袖領域の値を決定

GPU プログラム（ラプラシアン kernel ）

 袖領域の設定

グローバルメモリにデータがある箇所はグロー バルメモリから読み込み

グローバルメモリにデータがない箇所は 2 階微

GPU プログラム（ラプラシアン kernel ）

共有メモリのデータを利用してラプラシアンを計算

if 分岐を排除

ulap[]

su[][]

GPU プログラムの評価

 共有メモリを使用した dif3.cu

 dif2.cu （共有メモリ不使用）より速いこともあれば遅いこ ともある

 Fermi 世代以降の GPU はキャッシュを搭載

 共有メモリを使っても速くならない

 共有メモリへ明示的にデータを移動

 余分な処理により負荷が増加

__global__ void update(double *u, double *unew){

int i,j,ij;

i = blockIdx.x*blockDim.x + threadIdx.x;

j = blockIdx.y*blockDim.y + threadIdx.y;

ij = i+Nx*j;

u[ij] = unew[ij];

}

その他の処理の高速化

 値の更新

 unew のデータを u にコピーしているだけ

 cudaMemcpy で代用可能

その他の処理の高速化

 値の更新

 unew のデータを u にコピーしているだけ

 cudaMemcpy で代用可能

for(n=0;n<Nt;n++){

laplacian<<<Block, Thread>>>(dev_u,dev_ulap);

integrate<<<Block, Thread>>>(dev_u,dev_ulap,dev_unew);

//update<<<Block, Thread>>>(dev_u,dev_unew);

cudaMemcpy(dev_u, dev_new, Nbytes, cudaMemcpyDeviceToDevice);

}

return 0;

その他の処理の高速化

 初期値の計算

 三角関数がそれらしい値で求められればいい場合

 有効数字 4 桁程度

__global__ void init(double *u, double *ulap, double *unew){

int i,j,ij;

double x,y;

i = blockIdx.x*blockDim.x + threadIdx.x;

j = blockIdx.y*blockDim.y + threadIdx.y;

ij = i+Nx*j;

x = (double)i*dx;y = (double)j*dy;

u[ij]=sin(x)*sin(y);

unew[ij]=0.0;

ulap[ij]=0.0;

}

その他の処理の高速化

 初期値の計算

 三角関数がそれらしい値で求められればいい場合

 有効数字 4 桁程度

__global__ void init(double *u, double *ulap, double *unew){

int i,j,ij;

double x,y;

i = blockIdx.x*blockDim.x + threadIdx.x;

j = blockIdx.y*blockDim.y + threadIdx.y;

ij = i+Nx*j;

x = (double)i*dx;y = (double)j*dy;

u[ij]=__sinf(x)*__sinf(y);

unew[ij]=0.0;

ulap[ij]=0.0;

}

FastMath 関数

 三角関数など数学関数を高速に計算

 精度は落ちる

 手動で書き換え

 __sinf

 __cosf

 __powf など

 ‐use_fast_math オプション

 sin,cos を利用していても， ‐use‐fast‐math オプションを

付けることで __sinf, __cosf に置き換えられる

付録

グローバルメモリにデータがある箇所はグローバルメモリから読み込み

 dif2.cu （共有メモリ不使用）より速いこともあれば遅いこともある

global void update(double u, double unew){

global void init(double u, double ulap, double *unew){

x = (double)idx;y = (double)jdy;

global void init(double u, double ulap, double *unew){

x = (double)idx;y = (double)jdy;

u[ij]=sinf(x)*sinf(y);

付けることで sinf, cosf に置き換えられる

＋＋＋＋

 キャッシュを搭載していてもグローバルメモリへのアクセスを伴う if 分岐は高負荷

＋＋＋＋＋＋

何らかの方法で境界条件を反映