偏微分方程式の差分計算 長岡技術科学大学電気電子情報工学専攻出川智啓

偏微分方程式の差分計算 長岡技術科学大学電気電子情報工学専攻出川智啓

今回の内容 差分法 階微分 階微分に対する差分法 次元拡散方程式 guplot による結果の表示 分岐の書き方による実行時間の変化 高速化に利用できるいくつかのテクニック 7

前回授業 ビットマップを使った画像処理 配列の 要素が物理的な配置に対応 配列の 要素に物理的なデータが定義 B G R 7

数値計算 ( 差分法 ) 計算機を利用して数学 物理学的問題の解を計算 微積分を計算機で扱える形に変換 差分法では微分を差分に置き換え 処理自体は単純 精度を上げるために計算量が増加 73

シミュレーションの歴史と進歩 このスライドは諸事情により空白です 74

数値計算 ( 差分法 ) 物理空間に観測点を設置 観測点での物理量の変化を計算 観測点群を配列として確保し その点で物理量を定義 75

偏微分方程式 ( 拡散方程式 ) 物質の拡散を表す方程式 水の中に落ちたインクの拡散 金属中の熱伝導等 時刻 t=0 における の分布 ( 初期値 ) が既知 時間進行に伴い がどのように変化するかを計算 時間積分しながら の分布を求める ) ( ) ( x t x t t x 76 ) ( ) ( ) ( y t y x x t y x t t y x

差分法 計算機で微分を計算する方法の一つ 微分の定義 x の関数 について x だけ離れた 点間の傾きを計算し 点の間隔を無限小に近づけたときの極限 差分近似 d dx lm 0 Δx 関数をある間隔でサンプリング ( x Δx) Δx ( x) その間隔 x が の変化に対して十分小さいと仮定 d ( x Δx) ( x) dx Δx 77

差分法 ( 理論的なお話 ) 差分の誤差 lm を排除したことでどの程度の誤差が入るのか d dx lm Δx0 ( x Δx) Δx ( x) ( x Δx) Δx ( x) 関数 のテイラー展開を利用 ( x 3 3 d Δx d Δx d Δx) ( x) Δx 3 dx! dx 3! dx d dx ( x Δx) Δx ( x) Δx Δx! d dx Δx 3! 3 d 3 dx 3 78

差分法 ( 理論的なお話 ) 空間打ち切り誤差 定義とテイラー展開の比較 d dx ( x Δx) ( x) ( x Δx) ( x) lm Δx 0 Δx Δx 誤差 d dx ( x Δx) Δx ( x) Δx Δx! d dx Δx 3! 3 d 3 dx 3 テイラー展開を有限項で打ち切ったことによる誤差 3 Δx d Δx d 3! dx 3! dx 79

階微分の離散化 テイラー展開を応用 x 方向に x 離れた 点で展開 ( x 3 3 Δx Δx Δx) ( x) Δx 3 x! x 3! x ( x 3 3 Δx Δx Δx) ( x) Δx 3 x! x 3! x 式を足すと 階微分が消滅 ( x Δx Δx) ( x Δx) ( x)! x 730

階微分の離散化 階微分の式に整理 x ( x Δx) ( x) Δx ( x Δx) (x) (x+x) (x) (x x) x x=0 x x x x+x x 73

階微分の離散化 階微分の式に整理 x ( x Δx) ( x) Δx ( x Δx) [] 階の中心差分 [+] ([+] *[]+[ ]) /(dx*dx) [ ] [] 73 dx =0 + サンプリングされた関数値を配列 [] で保持

階微分の離散化 サンプリングされた関数値 の保持 x 方向長さ ( 計算領域の長さ ) を等間隔 xに分割 間隔 xの両端に離散点 ( 観測点 ) を置き そこでの値を定義 原点から順に離散点に番号を付与し で表現 離散点 におけるの値を配列 [] に保存 [] 配列 ( メモリ ) 配列 [] 離散点 =0 + 733

差分法の実装 計算領域内部 ddx[]=([ ] *[]+[+])/dxdx; 境界条件 ( 関数値が無いため処理を変更 ) ddx[0 ]=(*[0 ] 5*[ ]+4*[ ] [3 ])/dxdx; ddx[n ]=(*[N ] 5*[N ]+4*[N 3] [N 4])/dxdx; [] + + + + + + Δx ddx[] 734

差分法の実装 計算領域内部 ddx[]=([ ] *[]+[+])/dxdx; 境界条件 ( 関数値が無いため処理を変更 ) ddx[0 ]=(*[0 ] 5*[ ]+4*[ ] [3 ])/dxdx; ddx[n ]=(*[N ] 5*[N ]+4*[N 3] [N 4])/dxdx; [] + Δx ddx[] 735

差分法の実装 計算領域内部 ddx[]=([ ] *[]+[+])/dxdx; 境界条件 ( 関数値が無いため処理を変更 ) ddx[0 ]=(*[0 ] 5*[ ]+4*[ ] [3 ])/dxdx; ddx[n ]=(*[N ] 5*[N ]+4*[N 3] [N 4])/dxdx; [] + Δx ddx[] 736

CPU プログラム #clude<stdlb.h> #clude<math.h>/* lmオプションが必要*/ #dee Lx (.0*M_PI) #dee Nx (56) #dee dx (Lx/(Nx )) #dee Nbytes (Nx*szeo(double)) #dee dxdx (dx*dx) vod t(double *){ t ; or(=0; <Nx; ++){ [] = s(*dx); vod deretate(double * double *ddx){ t ; ddx[0] = (.0*[0] 5.0*[] +4.0*[] [3])/dxdx; or(=; <Nx ; ++) ddx[] = ( [+].0*[ ] + [ ])/dxdx; ddx[nx ]=( [Nx 4] +4.0*[Nx 3] 5.0*[Nx ] +.0*[Nx ])/dxdx; t ma(vod){ double **ddx; = (double *)malloc(nbytes); ddx = (double *)malloc(nbytes); t(); deretate(ddx); retur 0; deretate.c 737

関数の離散化 計算領域の長さと離散点の数 離散点の間隔の関係 (x) L x x 0 x 0 から L x の間に設けられた点の数 N x =L x /x + 738

実行結果 d/dx 739 x

GPU への移植 計算領域内部を計算するスレッド ddx[]=([ ] *[]+[+])/dxdx; 境界を計算するスレッド ddx[0 ]=(*[0 ] 5*[ ]+4*[ ] [3 ])/dxdx; ddx[n ]=(*[N ] 5*[N ]+4*[N 3] [N 4])/dxdx; [] + + + + + + + + Δx ddx[] 740

GPU プログラム #clude<stdo.h> #clude<stdlb.h> #clude<math.h>/* lmオプションが必要*/ #dee Lx (.0*M_PI) #dee Nx (56) #dee dx (Lx/(Nx )) #dee Nbytes (Nx*szeo(double)) #dee dxdx (dx*dx) #dee NT (8) #dee NB (Nx/NT) vod t(double *){ t ; or(=0; <Nx; ++){ [] = s(*dx); global vod deretate (double * double *ddx){ t = blockidx.x*blockdm.x + threadidx.x; (==0) ddx[] = (.0*[ ] 5.0*[+] +4.0*[+] [+3])/dxdx; (0< && <Nx ) ddx[] = ( [+].0*[ ] + [ ])/dxdx; (==Nx ) ddx[]=( [ 3] +4.0*[ ] 5.0*[ ] +.0*[ ])/dxdx; deretate.cu 74

GPU プログラム t ma(vod){ double *host_*host_ddx; double **ddx; host_ =(double *)malloc(nbytes); cudamalloc((vod **)&Nbytes); cudamalloc((vod **)&ddxnbytes); t(host_); cudamemcpy( host_ Nbytes cudamemcpyhosttodevce); deretate<<<nb NT>>>( ddx); //host_ddx=(double *)malloc(nbytes); //cudamemcpy(host_ddx ddx Nbytes cudamemcpydevcetohost); //or(t =0; <Nx; ++)prt("%%% "*dxhost_[]host_ddx[]); ree(host_); ree(host_ddx); cudafree(); cudafree(ddx); retur 0; deretate.cu 74

次元への拡張 743 拡散方程式 次元 次元 ) ( ) ( x t x t t x ) ( ) ( ) ( y t y x x t y x t t y x deretate.cu で計算どのように 次元に拡張するかどのように離散化するか

次元への拡張 x 方向 階偏微分 y 方向を固定して x 方向に偏微分 y 方向 階偏微分 x 方向を固定して y 方向に偏微分 ) ( ) ( ) ( Δx y Δx x y x y Δx x x ) ( ) ( ) ( Δy Δy y x y x Δy y x y 744

時間積分 時間微分項の離散化 時間微分項を前進差分で離散化 ( x y t Δt) ( x y t) t Δt 右辺のt+tの項を移行 ( x y t Δt) ( x y t) Δt t t 拡散方程式を代入 ( x y t Δt) ( x y t) Δt ( x y t) (x y t) の 階微分を計算できれば (x y t+t) が求められる 745

離散化された方程式の記述 簡略化した表現 配列との対応をとるため下付き添字 を利用 ( x y) ( x Δx y) ( x y Δy) 時間は 上付き添字 を利用 ( x y t) ( x y t Δt) 746

離散化された拡散方程式 連続系 離散系 t 秒後の値 ) ( ) ( ) ( y t y x x t y x t t y x Δy Δx Δt Δy Δx Δt 747

離散化された関数値の保持 各方向を等間隔 x yに分割し 離散点を配置 離散点で関数値を定義し 配列に保存 配列 [][] =0 + =0 + 離散点 配列 ( メモリ ) 748

拡散方程式 ( 熱伝導方程式 ) x y + + x y x y x y t t+t 749 Δy Δx Δt + +

拡散方程式の安定性 プログラムを正しく作成しても正常な計算結果が得られない場合がある 安定条件 拡散の強さを表す係数 ( 拡散係数 ) を使った形 Δt 二つの条件を満たすことが必要条件 ( 結果が正しいかは別 ) Δt 0.5 Δt v v 5 Δx Δy 0. Δx Δy 750

計算手順. 計算条件の決定 計算領域の大きさ L x L y 計算領域の分割数 ( 離散点の個数 )N x N y 離散点同士の間隔 ( 格子間隔 )x y 計算時間間隔 t. 初期値の決定 の初期分布の決定 3. 差分値の計算 の分布からx y 方向の 階微分値を計算 境界条件に基づいて境界の値を決定 t 秒後のを計算 75

CPU プログラム 計算条件の決定 計算領域の大きさ L x L y 計算領域の分割数 ( 離散点の個数 )N x N y 離散点同士の間隔 ( 格子間隔 )x y #clude<stdo.h> #clude<stdlb.h> #clude<math.h> #dee Lx (.0) #dee Ly (.0) #dee Nx 8 #dee Ny 8 #dee dx (Lx/(Nx )) #dee dy (Ly/(Ny )) #dee dt 0.000 計算時間間隔 t 75

CPU プログラム 初期条件 境界条件 753 0.) ( 0.) ( 0 r r x y r 0 : 外向き法線ベクトル Δx 0 0 Δy y x r

CPU プログラム #clude<stdlb.h> #clude<math.h> #dee Lx (.0) #dee Ly (.0) #dee Nx 8 #dee Ny 8 #dee dx (Lx/(Nx )) #dee dy (Ly/(Ny )) #dee dt 0.000 #dee edt (.0) #dee Nt (t)(edt/dt) #dee DIFF (0.0) #dee dxdx (dx*dx) #dee dydy (dy*dy) #dee Nbytes (Nx*Ny*szeo(double)) vod t(double * double * double *); vod laplaca(double * double *); vod tegrate(double * double * double *); vod update(double * double *); t ma(vod){ double **_ew*_lapxy; t ; = (double *)malloc(nx*ny*szeo(double)); _ew = (double *)malloc(nx*ny*szeo(double)); _lap = (double *)malloc(nx*ny*szeo(double)); t( _lap _ew); or(=0;<nt;++){ laplaca(_lap); tegrate(_lap_ew); update(_ew); retur 0; duso.c 754

CPU プログラム vod t(double * double *_lap double *_ew){ t ; double xy; or(=0;<ny;++){ or(=0;<nx;++){ [+Nx*] = 0.0; _lap[+nx*] = 0.0; _ew[+nx*] = 0.0; x=*dx Lx/.0; y=*dy Ly/.0; (sqrt(x*x+y*y)<=0.) [+Nx*] =.0; vod laplaca(double * double *_lap){ t pmpm; or(=0;<ny;++){ or(=0;<nx;++){ = +Nx*; m= +Nx* ; ( < 0 )m=++nx*; p=++nx* ; (+>=Nx)p= +Nx*; m= +Nx*( );( < 0 )m=+nx*(+); p= +Nx*(+);(+>=Ny)p=+Nx*( ); _lap[] = ([p].0*[]+[m])/dxdx +([p].0*[]+[m])/dydy; vod tegrate (double * double *_lap double *_ew){ t ; or(=0;<ny;++){ or(=0;<nx;++){ = +Nx*; _ew[] = [] + dt*diff*_lap[]; vod update(double * double *_ew){ t ; or(=0;<ny;++){ or(=0;<nx;++){ = +Nx*; [] = _ew[]; duso.c 755

CPU プログラム 差分計算 x 方向 y 方向偏微分を個別に計算して加算 vod laplaca(double * double *_lap){ t pmpm; or(=0;<ny;++){ or(=0;<nx;++){ = +Nx* ; m= +Nx* ; ( < 0 )m=++nx*; p=++nx* ; (+>=Nx)p= +Nx*; m= +Nx*( );( < 0 )m=+nx*(+); p= +Nx*(+);(+>=Ny)p=+Nx*( ); _lap[] = ([p].0*[]+[m])/dxdx +([p].0*[]+[m])/dydy; 756

ラプラシアン計算のメモリ参照 ある 点 のラプラシアンを計算するために 周囲 5 点の を参照 [] _lap[] + + 757

ラプラシアン計算のメモリ参照 ある 点 のラプラシアンを計算するために 周囲 5 点の を参照 [] _lap[] + + 758

ラプラシアン計算のメモリ参照 ある 点 のラプラシアンを計算するために 周囲 5 点の を参照 [] _lap[] + + 759

ラプラシアン計算のメモリ参照 ある 点 のラプラシアンを計算するために 周囲 5 点の を参照 [] _lap[] + + 760

ラプラシアン計算のメモリ参照 ある 点 のラプラシアンを計算するために 周囲 5 点の を参照 全ての を参照し 領域内部の _lap を計算 [] _lap[] 76

CPU プログラム 境界条件 法線方向の 階微分が 0 境界での 階の差分式 76 _lap[] 0 : 外向き法線ベクトル Δx 0 0 Δy or Δx Δx Δx or Δy Δy Δy

CPU プログラム 境界条件 参照する点が境界からはみ出している場合は添字を修正 階微分を計算する式は変更なし _lap[] ( < 0 )m=++nx*; (+>=Nx)p= +Nx*; ( < 0 )m= +Nx*(+); (+>=Ny)p= +Nx*( ); 右端では+を に修正 左端では を+に修正 上端では+を に修正 下端では を+に修正 763

CPU プログラム の積分 vod tegrate(double * double *_lap double *_ew){ t ; or(=0;<ny;++){ or(=0;<nx;++){ = +Nx*; _ew[] = [] + dt*diff*_lap[]; 764 Δy Δx Δt

CPU プログラム の更新 から + を計算 + から + を計算 同じアルゴリズム 今の時刻から次の時刻を求める 求められた次の時刻を今の時刻と見なし 次の時刻を求める vod update(double * double *_ew){ t ; or(=0;<ny;++){ or(=0;<nx;++){ = +Nx*; [] = _ew[]; 765

GPU プログラム (CPU 処理用共通部分 ) #clude<stdo.h> #clude<stdlb.h> #clude<math.h> #dee Lx (.0) #dee Ly (.0) #dee Nx 048 #dee Ny Nx #dee dx (Lx/(Nx )) #dee dy (Ly/(Ny )) #dee dt 0.00000 #dee edt (.0) #dee Nt (t)(edt/dt) #dee DIFF (0.0) #dee dxdx (dx*dx) #dee dydy (dy*dy) #dee Nbytes (Nx*Ny*szeo(double)) #clude "d.cu" //#clude "d.cu" t ma(vod){ t ; double *dev_*dev ew*dev lap; dm3 Thread Block; (DIFF*dt/dxdx > 0.5){ prt("cogurato error "); ext(); cudamalloc( (vod**)&dev_ Nbytes ); cudamalloc( (vod**)&dev ew Nbytes ); cudamalloc( (vod**)&dev lap Nbytes ); Thread = dm3(threadxthready); Block = dm3(blockx BLOCKY ); t<<<block Thread>>> (dev_ dev lap dev ew); or(=;<=nt;++){ laplaca<<<block Thread>>>(dev_dev lap); tegrate<<<block Thread>>> (dev_dev lapdev ew); update<<<block Thread>>>(dev_dev ew); retur 0; duso.cu 766

GPU プログラム ( スレッドが 点を計算 ) #dee THREADX 6 #dee THREADY 6 #dee BLOCKX (Nx/THREADX) #dee BLOCKY (Ny/THREADY) global vod laplaca(double * double *_lap){ t mpmp; = blockidx.x*blockdm.x + threadidx.x; = blockidx.y*blockdm.y + threadidx.y; = +Nx* ; m= +Nx* ; ( < 0 )m=++nx*; p=++nx* ; (+>=Nx)p= +Nx*; m= +Nx*( );( < 0 )m=+nx*(+); p= +Nx*(+);(+>=Ny)p=+Nx*( ); _lap[] = ([p].0*[]+[m])/dxdx +([p].0*[]+[m])/dydy; global vod tegrate (double * double *_lap double *_ew){ t ; = blockidx.x*blockdm.x + threadidx.x; = blockidx.y*blockdm.y + threadidx.y; = +Nx*; _ew[] = [] + dt*diff*_lap[]; global vod update(double * double *_ew){ t ; = blockidx.x*blockdm.x + threadidx.x; = blockidx.y*blockdm.y + threadidx.y; = +Nx*; [] = _ew[]; global vod t (double * double *_lap double *_ew){ t ; double xy; = blockidx.x*blockdm.x + threadidx.x; = blockidx.y*blockdm.y + threadidx.y; = +Nx*; [] = 0.0; _lap[] = 0.0; _ew[] = 0.0; x=*dx Lx/.0; y=*dy Ly/.0; (sqrt(x*x+y*y)<=0.) [] =.0; d.cu 767

guplot による結果の表示 次元 3 次元データをプロットするアプリケーション コマンドラインで命令を実行してグラフを描画 関数の描画 ファイルから読み込んだデータの表示が可能 tesla?? では正しく動作しないため grouse で実行すること 768

guplot による結果の表示 表示可能なファイルフォーマット データはスペース区切り 3 次元表示では列 ( または行 ) の区切りとして空白行が必要 0.500000 0.500000 0.000000e+00 0.4990 0.500000 0.000000e+00... 0.4990 0.500000 0.000000e+00 0.500000 0.500000 0.000000e+00 0.500000 0.4990 0.000000e+00 0.4990 0.4990 0.000000e+00... 0.4990 0.4990 0.000000e+00 0.500000 0.4990 0.000000e+00 0.500000 0.498045 0.000000e+00 0.4990 0.498045 0.000000e+00 < 改行 < 改行 769

guplot による結果の表示 スクリプトの実行 guplotのコマンドを記述したファイルを実行 load ' スクリプトファイル名 ' am.gpl 結果を単色の等値面で 3 次元表示 am_color.gpl 結果をカラーの等値面で 3 次元表示 am_d.gpl 結果をカラーの等高線で 次元表示 770

guplot による結果の表示 スクリプト am.gpl の内容 set xrage [ 0.5:0.5] x 軸の表示範囲を 0.5~0.5に固定 set yrage [ 0.5:0.5] y 軸の表示範囲を 0.5~0.5に固定 set zrage [0:] z 軸の表示範囲を0~に固定 set tcslevel 0 xy 平面とz 軸の最小値表示位置の差 set hdde3d 隠線消去をする set vew 4530 視点をx 方向 45 y 方向 30 に設定 set sze 0.75 グラフの縦横比を:0.75 uset cotour 次元等高線は表示しない set surace 3 次元で等値面を表示 uset pm3d カラー表示はしない 以下でファイルを読み込み 3 次元表示 splot '0.00.txt' usg ::3 wth le ttle "t=0.00s"... splot '.00.txt' usg ::3 wth le ttle "t=.00s" 77

guplot による結果の表示 スクリプト am.gpl による表示結果 77

guplot による結果の表示 スクリプト am_color.gpl の内容 set xrage [ 0.5:0.5] x 軸の表示範囲を 0.5~0.5に固定 set yrage [ 0.5:0.5] y 軸の表示範囲を 0.5~0.5に固定 set zrage [0:] z 軸の表示範囲を0~に固定 set tcslevel 0 xy 平面とz 軸の最小値表示位置の差 set hdde3d 隠線消去をする set vew 4530 視点をx 方向 45 y 方向 30 に設定 set sze 0.75 グラフの縦横比を:0.75 uset cotour 次元等高線は表示しない set surace 3 次元で等値面を表示 set pm3d カラー表示する set cbrage[0:] カラーバーの表示範囲を0~に固定 以下でファイルを読み込み 3 次元表示 splot '0.00.txt' usg ::3 wth le ttle "t=0.00s"... splot '.00.txt' usg ::3 wth le ttle "t=.00s" 773

guplot による結果の表示 スクリプト am_color.gpl による表示結果 774

guplot による結果の表示 スクリプト am_d.gpl の内容 set xrage [ 0.5:0.5] x 軸の表示範囲を 0.5~0.5に固定 set yrage [ 0.5:0.5] y 軸の表示範囲を 0.5~0.5に固定 set zrage [0:] z 軸の表示範囲を0~に固定 set vew 00 真上から見下ろす set sze グラフの縦横比を: set cotour 次元等高線を表示 uset surace 3 次元で等値面を表示しない uset pm3d カラー表示しない set ctrparam levels cremetal 00. 等高線を0からまで0. 刻みに設定以下でファイルを読み込み 3 次元表示 splot '0.00.txt' usg ::3 wth le ttle "t=0.00s"... splot '.00.txt' usg ::3 wth le ttle "t=.00s" 775

guplot による結果の表示 スクリプト am_d.gpl による表示結果 776

様々な境界条件の計算 共有メモリを利用してキャッシュを模擬 共有メモリに付加的な領域を追加 袖領域 [] _lap[] 777

様々な境界条件の計算 共有メモリを利用してキャッシュを模擬 共有メモリに付加的な領域を追加 [] _lap[] 778

様々な境界条件の計算 共有メモリを利用してキャッシュを模擬 共有メモリに付加的な領域を追加 s[][] [] _lap[] 779

付加的な領域 ( 袖領域 ) の取り扱い データがグローバルメモリに存在する場合は グローバルメモリから読み込み s[][] [] _lap[] 780

付加的な領域 ( 袖領域 ) の取り扱い グローバルメモリに無い場合は境界条件から決定 s[][] [] _lap[] 78

境界条件 78 階微分が 0 0 Δy y 0 Δx x u

境界条件 783 階微分が 0 0 Δx x 0 Δy y

境界条件 784 階微分を片側差分で計算 3 4 6 4 3 4 5 Δy y 3 4 5 Δx x 3 4 6 4

GPU プログラム ( ラプラシアン kerel) #dee THREADX 6 #dee THREADY 6 #dee BLOCKX (Nx/THREADX) #dee BLOCKY (Ny/THREADY) global vod laplaca(double * double *_lap){ t ; t txty; shared loat s[+threadx+][+thready+]; = blockidx.x*blockdm.x + threadidx.x; = blockidx.y*blockdm.y + threadidx.y; tx = threadidx.x + ; ty = threadidx.y + ; = +Nx*; s[tx][ty] = []; sycthreads(); (blockidx.x == 0 && threadidx.x == 0) s[tx ][ty] = s[tx+][ty]; (blockidx.x!= 0 && threadidx.x == 0) s[tx ][ty] = [ +Nx*]; (blockidx.x == grddm.x && threadidx.x == blockdm.x ) s[tx+][ty] = s[tx ][ty]; (blockidx.x!= grddm.x && threadidx.x == blockdm.x ) s[tx+][ty] = [++Nx*]; (blockidx.y == 0 && threadidx.y == 0) s[tx][ty ] = s[tx][ty+]; (blockidx.y!= 0 && threadidx.y == 0) s[tx][ty ] = [+Nx*( )]; (blockidx.y == grddm.y && threadidx.y == blockdm.y ) s[tx][ty+] = s[tx][ty ]; (blockidx.y!= grddm.y && threadidx.y == blockdm.y ) s[tx][ty+] = [+Nx*(+)]; sycthreads(); _lap[] = (s[tx ][ty ].0*s[tx][ty]+s[tx+][ty ])/dxdx +(s[tx ][ty ].0*s[tx][ty]+s[tx ][ty+])/dydy; d.cu 785

GPU プログラム ( ラプラシアン kerel) shared loat s[+threadx+][+thready+]; = blockidx.x*blockdm.x + threadidx.x; = blockidx.y*blockdm.y + threadidx.y; tx = threadidx.x; ty = threadidx.y; [] = +Ny*; s[tx][ty] = []; sycthreads(); ブロック内のスレッド数 + 袖領域分の共有メモリを確保袖領域があるために添字の対応が変化 添字の対応を考えないと 必要なデータを袖領域に置いてしまう sycthreads() を呼んでスレッドを同期 共有メモリにデータが正しく書き込まれた事を保証 s[][] 786

GPU プログラム ( ラプラシアン kerel) shared loat s[+threadx+][+thready+]; = blockidx.x*blockdm.x + threadidx.x; = blockidx.y*blockdm.y + threadidx.y; tx = threadidx.x + ; ty = threadidx.y; [] = +Ny*; s[tx][ty] = []; sycthreads(); ブロック内のスレッド数 + 袖領域分の共有メモリを確保袖領域があるために添字の対応が変化 添字の対応を考えないと 必要なデータを袖領域に置いてしまう sycthreads() を呼んでスレッドを同期 787 共有メモリにデータが正しく書き込まれた事を保証 s[][] +

GPU プログラム ( ラプラシアン kerel) shared loat s[+threadx+][+thready+]; = blockidx.x*blockdm.x + threadidx.x; = blockidx.y*blockdm.y + threadidx.y; tx = threadidx.x + ; ty = threadidx.y + ; [] = +Ny*; s[tx][ty] = []; sycthreads(); ブロック内のスレッド数 + 袖領域分の共有メモリを確保袖領域があるために添字の対応が変化 添字の対応を考えないと 必要なデータを袖領域に置いてしまう sycthreads() を呼んでスレッドを同期 788 共有メモリにデータが正しく書き込まれた事を保証 s[][] + +

GPU プログラム ( ラプラシアン kerel) 袖領域の設定 (blockidx.x == 0 && threadidx.x == 0 ) s[tx ][ty] = s[tx+][ty]; (blockidx.x!= 0 && threadidx.x == 0 ) s[tx ][ty] = [ +Nx*]; (blockidx.x == grddm.x && threadidx.x == blockdm.x ) s[tx+][ty] = s[tx ][ty]; (blockidx.x!= grddm.x && threadidx.x == blockdm.x ) s[tx+][ty] = [++Nx*]; (blockidx.y == 0 && threadidx.y == 0 ) s[tx][ty ] = s[tx][ty+]; (blockidx.y!= 0 && threadidx.y == 0 ) s[tx][ty ] = [+Nx*( )]; (blockidx.y == grddm.y && threadidx.y == blockdm.y ) s[tx][ty+] = s[tx][ty ]; (blockidx.y!= grddm.y && threadidx.y == blockdm.y ) s[tx][ty+] = [+Nx*(+)]; sycthreads(); グローバルメモリにデータがある箇所はグローバルメモリから読み込みグローバルメモリにデータがない箇所は 階微分が0になるように袖領域の値を決定 ブロックが境界に接しているか否かで処理を切替 789

GPU プログラム ( ラプラシアン kerel) _lap[] = (s[tx ][ty ].0*s[tx][ty]+s[tx+][ty ])/dxdx +(s[tx ][ty ].0*s[tx][ty]+s[tx ][ty+])/dydy; 共有メモリのデータを利用してラプラシアンを計算 分岐を排除 _lap[] s[][] 790

GPU プログラムの評価 共有メモリを使用した d.cu d.cu( 共有メモリ不使用 ) と比較して遅い Ferm 世代以降の GPU はキャッシュを搭載 共有メモリを使っても速くならない 共有メモリへ明示的にデータを移動 余分な処理により負荷が増加 せっかくなので 文の書き方について評価 79

GPU プログラムの評価 分岐の書き方による実行速度の変化 ブロックとスレッドの条件を同時に記述 (blockidx.x == 0 && threadidx.x == 0 ) s[tx ][ty] = s[tx+][ty]; (blockidx.x!= 0 && threadidx.x == 0 ) s[tx ][ty] = [ +Nx*]; (blockidx.x == grddm.x && threadidx.x == blockdm.x ) s[tx+][ty] = s[tx ][ty]; (blockidx.x!= grddm.x && threadidx.x == blockdm.x ) s[tx+][ty] = [++Nx*]; (blockidx.y == 0 && threadidx.y == 0 ) s[tx][ty ] = s[tx][ty+]; (blockidx.y!= 0 && threadidx.y == 0 ) s[tx][ty ] = [+Nx*( )]; (blockidx.y == grddm.y && threadidx.y == blockdm.y ) s[tx][ty+] = s[tx][ty ]; (blockidx.y!= grddm.y && threadidx.y == blockdm.y ) s[tx][ty+] = [+Nx*(+)]; Warp 内のスレッドが分岐すると実行速度が著しく低下 GPU はブロック単位で分岐 Warp 単位で分岐しても実行速度の低下は抑えられる 上の書き方は最適ではない? ブロック単位の分岐とスレッド単位の分岐を同時に記述 79

GPU プログラムの評価 分岐の書き方による実行速度の変化 ブロックとスレッドの条件を分け を入れ子に (blockidx.x == 0 )(threadidx.x == 0 ) s[tx ][ty] = s[tx+][ty]; (blockidx.x!= 0 )(threadidx.x == 0 ) s[tx ][ty] = [ +Nx*]; (blockidx.x == grddm.x )(threadidx.x == blockdm.x ) s[tx+][ty] = s[tx ][ty]; (blockidx.x!= grddm.x )(threadidx.x == blockdm.x ) s[tx+][ty] = [++Nx*]; (blockidx.y == 0 )(threadidx.y == 0 ) s[tx][ty ] = s[tx][ty+]; (blockidx.y!= 0 )(threadidx.y == 0 ) s[tx][ty ] = [+Nx*( )]; (blockidx.y == grddm.y )(threadidx.y == blockdm.y ) s[tx][ty+] = s[tx][ty ]; (blockidx.y!= grddm.y )(threadidx.y == blockdm.y ) s[tx][ty+] = [+Nx*(+)]; ブロック単位で分岐した後 スレッド単位で分岐 ブロック単位での分岐による速度低下が抑えられそうな気がする ブロック単位の分岐を記述する 8 個の は 個一組 ( 同じ色の行 ) で排他的に分岐 blockidx.x==0 が成立すると blockidx!=0 は絶対に成立しない 793

GPU プログラムの評価 分岐の書き方による実行速度の変化 個一組の の片方を else に書き換え (blockidx.x == 0 ){(threadidx.x == 0 ) s[tx ][ty] = s[tx+][ty]; else {(threadidx.x == 0 ) s[tx ][ty] = [ +Nx*]; (blockidx.x == grddm.x ){(threadidx.x == blockdm.x ) s[tx+][ty] = s[tx ][ty]; else {(threadidx.x == blockdm.x ) s[tx+][ty] = [++Nx*]; (blockidx.y == 0 ){(threadidx.y == 0 ) s[tx][ty ] = s[tx][ty+]; else {(threadidx.y == 0 ) s[tx][ty ] = [+Nx*( )]; (blockidx.y == grddm.y ){(threadidx.y == blockdm.y ) s[tx][ty+] = s[tx][ty ]; else {(threadidx.y == blockdm.y ) s[tx][ty+] = [+Nx*(+)]; 個一組の における条件判断は 回限り 回 の条件判断を行うよりは速いような気がする 794

分岐の書き方の違いによる実行速度の変化 格子分割数 N x N y =048 048 ブロックあたりのスレッド数 6 6 文の書き方 ブロックとスレッドの条件を同時に記述ブロックとスレッドの条件を分け を入れ子に 個一組のの片方を elseに書き換えカーネル ( 共有メモリ不使用 ) 実行時間 [ms].9.49.33.00 795

その他の処理の高速化 laplaca カーネルと tegrate カーネルの融合 カーネルフュージョン (kerel uso) 異なる処理を行うカーネルをまとめることでデータを再利用 グローバルメモリへの書き込み 読み出しを回避 global vod tegrate(double * double *_ew){ t ; double _lap; // ある 点のラプラシアンを計算してレジスタに保存 = blockidx.x*blockdm.x + threadidx.x; = blockidx.y*blockdm.y + threadidx.y; = +Nx*; _lap =... _ew[] = [] + dt*diff*_lap; //[]=[]+... とすると 全ての点の _lap を計算する前に を更新する可能性がある 796

その他の処理の高速化 値の更新 _ewのデータをにコピーしているだけ cudamemcpyで代用可能 global vod update(double * double *_ew){ t ; = blockidx.x*blockdm.x + threadidx.x; = blockidx.y*blockdm.y + threadidx.y; = +Nx*; [] = _ew[]; 797

その他の処理の高速化 値の更新 _ewのデータをにコピーしているだけ cudamemcpyで代用可能 or(=;<=nt;++){ tegrate<<<block Thread>>>(dev_dev ew); //update<<<block Thread>>>(dev_dev ew); cudamemcpy(dev_ dev ew Nbytes cudamemcpydevcetodevce); retur 0; 798

その他の処理の高速化 値の更新 ダブルバッファリング ( 第 7 回授業参照 ) でコピーを回避することも可能 double *dev_[]; t =0;out= ; cudamalloc((vod **)&dev_[ ]Nbytes); cudamalloc((vod **)&dev_[out]nbytes); t<<<block Thread>>>(dev_[(=0)] dev_[out(=)]); or(=;<=nt;++){//= tegrate<<<block Thread>>>(dev_[(=0)]dev_[out(=)]); = out; out = ; //= out=0 799

その他の処理の高速化 値の更新 ダブルバッファリング ( 第 7 回授業参照 ) でコピーを回避することも可能 double *dev_[]; t =0;out= ; cudamalloc((vod **)&dev_[ ]Nbytes); cudamalloc((vod **)&dev_[out]nbytes); t<<<block Thread>>>(dev_[(=0)] dev_[out(=)]); or(=;<=nt;++){//= tegrate<<<block Thread>>>(dev_[(=)]dev_[out(=0)]); = out; out = ; //=0 out= 800

その他の処理の高速化 値の更新 ダブルバッファリング ( 第 7 回授業参照 ) でコピーを回避することも可能 tegrate<<<block Thread>>>(dev_[(=0)]dev_[out(=)]); dev_[0] 時刻 dev_[] 時刻 + 80

その他の処理の高速化 値の更新 ダブルバッファリング ( 第 7 回授業参照 ) でコピーを回避することも可能 tegrate<<<block Thread>>>(dev_[(=)]dev_[out(=0)]); dev_[0] 時刻 + dev_[] 時刻 80