Microsoft PowerPoint - 先端GPGPUシミュレーション工学特論(web).pptx

Transcription

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

2 今日の内容 シミュレーションの歴史と進歩 差分法 1 階微分 階微分に対する差分法 1 次関数の差分 次元拡散方程式 付録 共有メモリの典型的な使い方 49 先端 GPGPUシミュレーション工学特論

3 数値計算 計算機を利用して数学 物理学的問題の解を計算 微積分を計算機で扱える形に変換 処理自体はあまり複雑ではない 精度を上げるために計算量が増加 493 先端 GPGPUシミュレーション工学特論

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

5 シミュレーションの歴史と進歩 1985 年 次元でのシミュレーション 実験結果と 傾向は 一致 495

6 シミュレーションの歴史と進歩 1987 年 3 次元化し 実車に近い形状で計算 496

7 シミュレーションの歴史と進歩 1988 年 ~1990 年 車輪や床下も含めたモデル化 497

8 シミュレーションの歴史と進歩 1993 年 空気抵抗を誤差 1% で予測可能 空力解析以外にも利用 498

9 シミュレーションの歴史と進歩 199 年 ~1993 年 車体から発生する騒音のシミュレーション ドアミラーやピラーの形状の改良 エンジンルームの冷却 10mm 以上の部品は全て含めて解析 499

10 シミュレーションの歴史と進歩 衝突解析 500

11 現象を支配する方程式 支配方程式 ) 現象は微分方程式によって記述 微分と積分が計算できれば現象を明らかにできる 支配方程式 場所や時間によって方程式は変わらない 流体の支配方程式 0 t p t q p E t E 質量保存式 ) エネルギ保存式 ) 運動量保存式 ) 501

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

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

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

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

16 差分法 理論的なお話 ) 誤差の主要項 空間打ち切り誤差の中で最も大きい誤差は第 1 項 は小さい の高次項はさらに小さく 無視できる Δ! d d Δ 3! d 3 d 3 O Δ) 直感的に導いた微分の数値計算法の誤差は O) を 1/10 にすれば誤差も 1/10 になる 506

17 差分法 理論的なお話 ) 507 差分の取り方 関数値の選び方にいくつか選択肢がある テイラー展開で整理 Δ Δ Δ Δ d d Δ ) ) ) ) lm !! ) ) d d Δ d d Δ d d Δ Δ !! 1 ) ) d d Δ d d Δ Δ Δ Δ d d

18 差分法 理論的なお話 ) 508 差分の取り方 テイラー展開で整理 Δ Δ Δ Δ Δ Δ d d Δ ) ) ) ) lm 0 Δ Δ Δ Δ ) ) ) ) ! 1 ) ) d d Δ Δ Δ Δ Δ d d

19 差分法の概念図 ) 中心差分 Δ) Δ Δ) 後退差分 ) Δ) Δ 前進差分 Δ) ) Δ =

20 差分法の概念図 ) サンプリングされた関数値を配列 [] で保持 = =0 1) +1) 510

21 差分法の概念図 [] 中心差分 [+1] [ 1]) /*d) サンプリングされた関数値を配列 [] で保持 [+1] [] [ 1] d =

22 階微分の離散化 テイラー展開を応用 方向に 離れた 点で展開 3 3 Δ Δ Δ) ) Δ 3! 3! 3 3 Δ Δ Δ) ) Δ 3! 3! 式を足すと 1 階微分が消滅 Δ Δ) Δ) )! 51 先端 GPGPUシミュレーション工学特論

23 階微分の離散化 階微分の式に整理 Δ) ) Δ Δ) [] 階の中心差分 [+1] [+1] *[]+[ 1]) /d*d) [ 1] [] = 先端 GPGPUシミュレーション工学特論 d サンプリングされた関数値を配列 [] で保持

24 差分法の実装 階微分の中心差分近似 d d Δ) ) Δ Δ) 1 Δ 1 [] Δ dd[] 514

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

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

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

28 CPU プログラム #nclde<stdlb.h> #nclde<math.h>/* lmオプションが必要*/ #defne L.0*M_PI) #defne N 56) #defne d L/N 1)) #defne Nbytes N*szeofdoble)) #defne dd d*d) vod ntdoble *){ nt ; for=0; <N; ++){ [] = sn*d); vod dfferentatedoble * doble *dd){ nt ; dd[0] =.0*[0] 5.0*[1] +4.0*[] [3])/dd; for=1; <N 1; ++) dd[] = [+1].0*[ ] + [ 1])/dd; dd[n 1]= [N 4] +4.0*[N 3] 5.0*[N ] +.0*[N 1])/dd; nt manvod){ doble **dd; = doble *)mallocnbytes); dd = doble *)mallocnbytes); nt); dfferentatedd); retrn 0; dfferentate.c 518

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

30 実行結果 d/d 50 先端 GPGPUシミュレーション工学特論

31 GPU への移植 計算領域内部を計算するスレッド dd[]=[ 1] *[]+[+1])/dd; 境界を計算するスレッド dd[0 ]=*[0 ] 5*[1 ]+4*[ ] [3 ])/dd; dd[n 1]=*[N 1] 5*[N ]+4*[N 3] [N 4])/dd; [] Δ dd[] 51

32 GPU プログラム #nclde<stdo.h> #nclde<stdlb.h> #nclde<math.h>/* lmオプションが必要*/ #defne L.0*M_PI) #defne N 56) #defne d L/N 1)) #defne Nbytes N*szeofdoble)) #defne dd d*d) #defne NT 18) #defne NB N/NT) vod ntdoble *){ nt ; for=0; <N; ++){ [] = sn*d); global vod dfferentate doble * doble *dd){ nt = blockid.*blockdm. + threadid.; f==0) dd[] =.0*[ ] 5.0*[+1] +4.0*[+] [+3])/dd; f0< && <N 1) dd[] = [+1].0*[ ] + [ 1])/dd; f==n 1) dd[]= [ 3] +4.0*[ ] 5.0*[ 1] +.0*[ ])/dd; dfferentate.c 5

33 GPU プログラム nt manvod){ doble *host_*host_dd; doble **dd; host_ =doble *)mallocnbytes); cdamallocvod **)&Nbytes); cdamallocvod **)&ddnbytes); nthost_); cdamemcpy host_ Nbytes cdamemcpyhosttodevce); dfferentate<<<nb NT>>> dd); //host_dd=doble *)mallocnbytes); //cdamemcpyhost_dd dd Nbytes cdamemcpydevcetohost); //fornt =0; <N; ++)prntf"%f%f%f n"*dhost_[]host_dd[]); freehost_); freehost_dd); cdafree); cdafreedd); retrn 0; dfferentate.c 53

34 次元への拡張 54 拡散方程式 1 次元 次元 ) ) t t t ) ) ) y t y t y t t y dfferentate.c で計算どのように 次元に拡張するかどのように離散化するか

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

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

37 離散化された方程式の記述 簡略化した表現 配列との対応をとるため下付き添字 を利用 y) y) 1 y y) 1 時間は 上付き添字 n を利用 n y t) y t t) n1 57 先端 GPGPUシミュレーション工学特論

38 離散化された拡散方程式 連続系 離散系 t 秒後の値 ) ) ) y t y t y t t y y t n n n n n n n n y t n n n n n n n n 58

39 拡散方程式 熱伝導方程式 ) y y y y t t+t y t n n n n n n n n

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

41 計算手順 1. 計算条件の決定 計算領域の大きさ L L y 計算領域の分割数 離散点の個数 )N N y 離散点同士の間隔 格子間隔 ) y 計算時間間隔 t. 初期値の決定 の初期分布の決定 3. 差分値の計算 の分布から y 方向の 階微分値を計算 境界条件に基づいて境界の値を決定 t 秒後のを計算 531 先端 GPGPUシミュレーション工学特論

42 CPU プログラム 計算条件の決定 計算領域の大きさ L L y 計算領域の分割数 離散点の個数 )N N y 離散点同士の間隔 格子間隔 ) y #nclde<stdo.h> #nclde<stdlb.h> #nclde<math.h> #defne L.0*M_PI) #defne Ly.0*M_PI) #defne N 51 #defne Ny 51 #defne d L/N 1)) #defne dy Ly/Ny 1)) #defne dt 計算時間間隔 t 53 先端 GPGPUシミュレーション工学特論

43 CPU プログラム 初期条件 sn sn y 境界条件 sn sn y 0 for=0;<ny;++){ for=0;<n;++){ = doble)*d; y = doble)*dy; [*Ny+]=sn)*sny); 533 先端 GPGPUシミュレーション工学特論

44 CPU プログラム #nclde<stdo.h> #nclde<stdlb.h> #nclde<math.h> #defne L.0*M_PI) #defne Ly.0*M_PI) #defne N 18 #defne Ny 18 #defne d L/N 1)) #defne dy Ly/Ny 1)) #defne dt #defne endt 1.0) #defne Nt nt)endt/dt) #defne DIFF 1.0) #defne dd d*d) #defne dydy dy*dy) #defne Nbytes N*Ny*szeofdoble)) vod laplacandoble * doble *); vod ntegratedoble * doble * doble *); vod pdatedoble * doble *); nt manvod){ doble **new*lapy; nt p1m1p1m1n; = doble *)mallocn*ny*szeofdoble)); new = doble *)mallocn*ny*szeofdoble)); lap = doble *)mallocn*ny*szeofdoble)); for=0;<ny;++){ for=0;<n;++){ = doble)*d; y = doble)*dy; [*Ny+]=sn)*sny); new[*ny+]=0.0f; lap[*ny+]=0.0f; forn=0;n<nt;n++){ laplacanlap); ntegratelapnew); pdatenew); retrn 0; dffson.c 534

45 CPU プログラム vod laplacandoble * doble *lap){ nt p1m1p1m1; for=1;<ny 1;++){ for=1;<n 1;++){ = *Ny+; p1 = +1)*Ny+; m1 = 1)*Ny+; p1 = *Ny++1; m1 = *Ny+ 1; lap[] = [p1].0*[]+[m1])/dd +[p1].0*[]+[m1])/dydy; vod ntegrate doble * doble *lap doble *new){ nt ; for=0;<ny;++){ for=0;<n;++){ = *Ny+; new[] = [] + dt*diff*lap[]; vod pdatedoble * doble *new){ nt ; for=0;<ny;++){ for=0;<n;++){ = *Ny+; [] = new[]; dffson.c 535

46 CPU プログラム 差分計算 方向 y 方向偏微分を個別に計算して加算 vod laplacandoble * doble *lap){ nt p1m1p1m1; for=1;<ny 1;++){ for=1;<n 1;++){ = *Ny+; p1 = +1)*Ny+; m1 = 1)*Ny+; p1 = *Ny++1; m1 = *Ny+ 1; lap[] = [p1].0*[]+[m1])/dd +[p1].0*[]+[m1])/dydy; 536 先端 GPGPUシミュレーション工学特論

47 ラプラシアン計算のメモリ参照 ある 1 点 のラプラシアンを計算するために 周囲 5 点の を参照 [] lap[] 先端 GPGPUシミュレーション工学特論

48 ラプラシアン計算のメモリ参照 ある 1 点 のラプラシアンを計算するために 周囲 5 点の を参照 [] lap[] 先端 GPGPUシミュレーション工学特論

49 ラプラシアン計算のメモリ参照 ある 1 点 のラプラシアンを計算するために 周囲 5 点の を参照 [] lap[] 先端 GPGPUシミュレーション工学特論

50 ラプラシアン計算のメモリ参照 ある 1 点 のラプラシアンを計算するために 周囲 5 点の を参照 [] lap[] 先端 GPGPUシミュレーション工学特論

51 ラプラシアン計算のメモリ参照 ある 1 点 のラプラシアンを計算するために 周囲 5 点の を参照 全ての を参照し 領域内部の lap を計算 [] lap[] 541 先端 GPGPUシミュレーション工学特論

52 CPU プログラム 境界条件 のラプラシアン sn sn y lap[] 境界ではどの時刻においても 0 =0 0 y = 0 y 0 y=0 0 y= 変数 lap を 0 で初期化すれば 計算しなくてもよい 54 先端 GPGPUシミュレーション工学特論

53 CPU プログラム の積分 vod ntegrate doble * doble *lap doble *new){ nt ; for=0;<ny;++){ for=0;<n;++){ = *Ny+; new[] = [] + dt*diff*lap[]; y t n n n n n n n n

54 CPU プログラム の更新 n から n+1 を計算 n+1 から n+ を計算 同じアルゴリズム 今の時刻から次の時刻を求める 求められた次の時刻を今の時刻と見なし 次の時刻を求める vod pdatedoble * doble *new){ nt ; for=0;<ny;++){ for=0;<n;++){ = *Ny+; [] = new[]; 544 先端 GPGPUシミュレーション工学特論

55 GPU プログラム CPU 処理用共通部分 ) #nclde<stdo.h> #nclde<stdlb.h> #nclde<math.h> #defne L.0*M_PI) #defne Ly.0*M_PI) #defne N 18 #defne Ny N #defne d L/N 1)) #defne dy Ly/Ny 1)) #defne dt #defne endt 1.0) #defne Nt nt)endt/dt) #defne DIFF 1.0) #defne dd d*d) #defne dydy dy*dy) #defne Nbytes N*Ny*szeofdoble)) #nclde "df1.c" //#nclde "df.c" //#nclde "df3.c" nt manvod){ nt n; doble *dev_*dev_new*dev_lap; dm3 Thread Block; fdiff*dt/dd > 0.5){ prntf"confgraton error n"); et1); cdamalloc vod**)&dev_ Nbytes ); cdamalloc vod**)&dev_new Nbytes ); cdamalloc vod**)&dev_lap Nbytes ); Thread = dm3threadxthready1); Block = dm3blockx BLOCKY 1); nt<<<block Thread>>> dev_ dev_lap dev_new); forn=0;n<nt;n++){ laplacan<<<block Thread>>> dev_dev_lap); ntegrate<<<block Thread>>> dev_dev_lapdev_new); pdate<<<block Thread>>>dev_dev_new); retrn 0; dff.c 545

56 GPU プログラム 1 スレッド版 ) #defne THREADX 1 #defne THREADY 1 #defne BLOCKX 1 #defne BLOCKY 1 global vod nt doble * doble *lap doble *new){ nt ; doble y; for=0;<ny;++){ for=0;<n;++){ = *Ny+; = doble)*d; y = doble)*dy; []=sn)*sny); new[]=0.0; lap[]=0.0; global vod laplacandoble * doble *lap){ nt p1m1p1m1; for=1;<ny 1;++){ for=1;<n 1;++){ = *Ny+; p1 = +1)*Ny+; m1 = 1)*Ny+; p1 = *Ny++1; m1 = *Ny+ 1; lap[] = [p1].0f*[]+[m1])/dd +[p1].0f*[]+[m1])/dydy; global vod ntegrate doble * doble *lap doble *new){ nt ; for=0;<ny;++){ for=0;<n;++){ = *Ny+; new[] = [] + dt*diff*lap[]; global vod pdatedoble * doble *new){ nt ; for=0;<ny;++){ for=0;<n;++){ = *Ny+; [] = new[]; df1.c 546

57 次元ブロック分割 1 スレッドが 1 点のラプラシアンを計算 grddm.= blockid.=0 blockid.=1 grddm.y= blockid.y=0 blockid.y=1 blockdm.=4 blockdm.y=4 threadid.= threadid.y= 547

58 次元ブロック分割 N=8 Ny=8 y 方向スレッド数 4 ブロック数 = blockid.*blockdm. + threadid. = blockid.y*blockdm.y + threadid.y block00) block10) = )10)0)30)00) 01)11)1)31) 0)1))3) 03)13)3)33) 33) 00) 00) 33) 33) = threadid. threadid.y 548 block01) block11)

59 次元的な配列アクセスの優先方向 CPU) 次元配列の 1 次元配列的表現 for=0;<n;++) for=0;<ny;++) ot[][]=n[][]; n[]ot[] for=0;<n;++) for=0;<ny;++) ot[*ny+]= n[*ny+]; Ny N 549

60 次元的な配列アクセスの優先方向 GPU) CUDA で 次元的に並列化してアクセスする場合 for=0;<n;++) for=0;<ny;++) ot[*ny+]= n[*ny+]; = blockid.*blockdm. + threadid.; = blockid.y*blockdm.y + threadid.y; ot[*n+]=n[*n+]; n[]ot[] threadid. threadid.y Ny N 550

61 GPU プログラム 1 スレッドが 1 点を計算 ) #defne THREADX 16 #defne THREADY 16 #defne BLOCKX N/THREADX) #defne BLOCKY Ny/THREADY) global vod laplacandoble * doble *lap){ nt p1m1p1m1; = blockid.*blockdm. + threadid.; = blockid.y*blockdm.y + threadid.y; f 0< && <N 1) && 0< && <Ny 1) ){ = +N*; p1 = +1+N*; m1 = 1+N*; p1 = +N*+1); m1 = +N* 1); lap[] = [p1].0*[]+[m1])/dd +[p1].0*[]+[m1])/dydy; global vod ntegrate doble * doble *lap doble *new){ nt ; = blockid.*blockdm. + threadid.; = blockid.y*blockdm.y + threadid.y; = +N*; new[] = [] + dt*lap[]; global vod pdatedoble * doble *new){ nt ; = blockid.*blockdm. + threadid.; = blockid.y*blockdm.y + threadid.y; = +N*; [] = new[]; global vod nt doble * doble *lap doble *new){ nt ; doble y; = blockid.*blockdm. + threadid.; = blockid.y*blockdm.y + threadid.y; = +N*; = doble)*d; y = doble)*dy; []=sn)*sny); new[]=0.0; lap[]=0.0; df.c 551

62 境界でラプラシアンが 0 にならない場合 共有メモリを利用してキャッシュを模擬 共有メモリに付加的な領域を追加 [] lap[] 55 先端 GPGPUシミュレーション工学特論

63 境界でラプラシアンが 0 にならない場合 共有メモリを利用してキャッシュを模擬 共有メモリに付加的な領域を追加 [] lap[] 553 先端 GPGPUシミュレーション工学特論

64 境界でラプラシアンが 0 にならない場合 共有メモリを利用してキャッシュを模擬 共有メモリに付加的な領域を追加 s[][] [] lap[] 554 先端 GPGPUシミュレーション工学特論

65 付加的な領域 袖領域 ) の取り扱い データがグローバルメモリに存在する場合は グローバルメモリから読み込み s[][] [] lap[] 555 先端 GPGPUシミュレーション工学特論

66 付加的な領域 袖領域 ) の取り扱い グローバルメモリに無い場合は境界条件から決定 s[][] [] lap[] 556 先端 GPGPUシミュレーション工学特論

67 境界条件 557 階微分が Δ Δy y

68 境界条件 558 階微分を片側差分で計算 Δy y Δ

69 GPU プログラム ラプラシアン kernel) #defne THREADX 16 #defne THREADY 16 #defne BLOCKX N/THREADX) #defne BLOCKY Ny/THREADY) global vod laplacandoble * doble *lap){ nt ; nt tty; shared float s[1+threadx+1][1+thready+1]; = blockid.*blockdm. + threadid.; = blockid.y*blockdm.y + threadid.y; t = threadid. + 1; ty = threadid.y + 1; = +N*; s[t][ty] = []; syncthreads); f blockid. == 0 && threadid. == 0 ) s[t 1][ty] =.0*s[t][ty] s[t+1][ty]; f blockid.!= 0 && threadid. == 0 ) s[t 1][ty] = [*N+ 1]; f blockid. == grddm. 1 && threadid. == blockdm. 1) s[t+1][ty] =.0*s[t][ty] s[t 1][ty]; f blockid.!= grddm. 1 && threadid. == blockdm. 1) s[t+1][ty] = [*N++1]; f blockid.y == 0 && threadid.y == 0) s[t][ty 1] =.0*s[t][ty] s[t][ty+1]; f blockid.y!= 0 && threadid.y == 0) s[t][ty 1] = [ 1)*N+]; f blockid.y == grddm.y 1 && threadid.y == blockdm.y 1) s[t][ty+1] =.0*s[t][ty] s[t][ty 1]; f blockid.y!= grddm.y 1 && threadid.y == blockdm.y 1) s[t][ty+1] = [+1)*N+]; syncthreads); lap[] = s[t 1][ty ].0*s[t][ty]+s[t+1][ty ])/dd +s[t ][ty 1].0*s[t][ty]+s[t ][ty+1])/dydy; df3.c 559

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

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

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

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

74 GPU プログラム ラプラシアン kernel) lap[] = s[t 1][ty ].0*s[t][ty]+s[t+1][ty ])/dd +s[t ][ty 1].0*s[t][ty]+s[t ][ty+1])/dydy; 共有メモリのデータを利用してラプラシアンを計算 f 分岐を排除 lap[] s[][] 564

75 GPU プログラムの評価 共有メモリを使用した df3.c df.c 共有メモリ不使用 ) より速いこともあれば遅いこともある Ferm 世代以降の GPU はキャッシュを搭載 共有メモリを使っても速くならない 共有メモリへ明示的にデータを移動 余分な処理により負荷が増加 565

76 その他の処理の高速化 値の更新 newのデータをにコピーしているだけ cdamemcpyで代用可能 global vod pdatedoble * doble *new){ nt ; = blockid.*blockdm. + threadid.; = blockid.y*blockdm.y + threadid.y; = +N*; [] = new[]; 566

77 その他の処理の高速化 値の更新 newのデータをにコピーしているだけ cdamemcpyで代用可能 forn=0;n<nt;n++){ laplacan<<<block Thread>>>dev_dev_lap); ntegrate<<<block Thread>>>dev_dev_lapdev_new); //pdate<<<block Thread>>>dev_dev_new); cdamemcpydev_ dev_new Nbytes cdamemcpydevcetodevce); retrn 0; 567

78 その他の処理の高速化 初期値の計算 三角関数がそれらしい値で求められればいい場合 有効数字 4 桁程度 global vod ntdoble * doble *lap doble *new){ nt ; doble y; = blockid.*blockdm. + threadid.; = blockid.y*blockdm.y + threadid.y; = +N*; = doble)*d;y = doble)*dy; []=sn)*sny); new[]=0.0; lap[]=0.0; 568

79 その他の処理の高速化 初期値の計算 三角関数がそれらしい値で求められればいい場合 有効数字 4 桁程度 global vod ntdoble * doble *lap doble *new){ nt ; doble y; = blockid.*blockdm. + threadid.; = blockid.y*blockdm.y + threadid.y; = +N*; = doble)*d;y = doble)*dy; []= snf)* snfy); new[]=0.0; lap[]=0.0; 569

80 FastMath 関数 三角関数など数学関数を高速に計算 精度は落ちる 手動で書き換え snf cosf powf など se_fast_math オプション sncos を利用していても se fast math オプションを付けることで snf cosf に置き換えられる 570

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

82 共有メモリの典型的な使い方 差分法 あるスレッドが中心差分を計算するために配列の要素 1 +1 を参照 配列要素は複数回参照される Ferm 世代以降はキャッシュが利用可能 [] 参照される回数 1 Δ dd[] 57

83 共有メモリの典型的な使い方 境界で異なる処理を行うために f 分岐が必要 キャッシュを搭載していてもグローバルメモリへのアクセスを伴う f 分岐は高負荷 データの再利用と f 文の排除に共有メモリを利用 [] 3 3 参照される回数 Δ dd[] 573

84 共有メモリによる明示的なキャッシュ グローバルメモリから共有メモリにデータをキャッシュ 共有メモリ上で境界条件を処理 中心差分の計算から f を排除 blockid.=0 blockid.=1 [] 共有メモリ 計算のために必要になる余分な領域 袖領域 ) dd[] 574

85 共有メモリによる明示的なキャッシュ グローバルメモリから共有メモリにデータをキャッシュ 共有メモリ上で境界条件を処理 中心差分の計算から f を排除 blockid.=0 blockid.=1 [] 共有メモリ 何らかの方法で境界条件を反映 何らかの方法で境界条件を反映 dd[] 575

86 共有メモリによる明示的なキャッシュ グローバルメモリから共有メモリにデータをキャッシュ 共有メモリ上で境界条件を処理 中心差分の計算から f を排除 blockid.=0 blockid.=1 [] 共有メモリ 全スレッドが同じ式で中心差分を計算 dd[] 576

87 GPU プログラム 単純な実装 ) #nclde<stdlb.h> #nclde<math.h>/* lm オプションが必要 */ #defne L.0*M_PI) #defne N 104*104) #defne d L/N 1)) #defne Nbytes N*szeofdoble)) #defne NT 56) #defne NB N/NT) vod ntdoble *){ nt ; for=0; <N; ++){ [] = sn*d); global vod dfferentate doble * doble *dd){ nt = blockid.*blockdm. + threadid.; f==0) dd[]= 3.0*[ ] +4.0*[+1] [+])/.0*d); f0< && <N 1) dd[] = [+1] [ 1])/.0*d); f==n 1) dd[]= [ ] 4.0*[ 1] +3.0*[ ])/.0*d); dfferentate1d.c 577

88 GPU プログラム 単純な実装 ) nt manvod){ doble *host_*host_dd; doble **dd; host_ =doble *)mallocnbytes); host_dd=doble *)mallocnbytes); cdamallocvod **)&Nbytes); cdamallocvod **)&ddnbytes); nthost_); cdamemcpyhost_nbytes cdamemcpyhosttodevce); dfferentate<<<nb NT>>>dd); cdamemcpyhost_dd dd Nbytes cdamemcpydevcetohost); freehost_); freehost_dd); cdafree); cdafreedd); retrn 0; dfferentate1d.c 578

89 共有メモリを用いた書き換え global vod dfferentatedoble * doble *dd){ nt = blockid.*blockdm. + threadid.; shared doble s[1+nt+1]; nt t = threadid.+1; s[t] = []; syncthreads); fblockid.> 0 && threadid.==0 ) s[t 1] = [ 1]; fblockid.< grddm. 1 && threadid.==blockdm. 1) s[t+1] = [+1]; fblockid.==0 && threadid.==0 ) s[t 1] = 3.0*s[t] 3.0*s[t+1]+s[t+]; fblockid.==grddm. 1 && threadid.==blockdm. 1) s[t+1] = 3.0*s[t] 3.0*s[t 1]+s[t ]; syncthreads); dd[] = s[t+1] s[t 1])/.0*d); dfferentate1d_shared.c 579

90 共有メモリの宣言と代入 nt = blockid.*blockdm. + threadid.; shared doble s[1+nt+1]; // 右と左の袖領域を追加して宣言 nt t = threadid.+1; s[t] = []; syncthreads); [] = NT=3 NT=3 threadid.= s[nt+] [0] [1] [] [3] [4] [5] t=

91 袖領域の処理 fblockid.> 0 && threadid.==0 ) s[t 1] = [ 1]; fblockid.< grddm. 1 && threadid.==blockdm. 1) s[t+1] = [+1]; [] blockid.=0 blockid.=1 = =threadid. s[nt+] [0] [1] [] [3] [4] [5] t=

92 袖領域の処理 fblockid.> 0 && threadid.==0 ) s[t 1] = [ 1]; fblockid.< grddm. 1 && threadid.==blockdm. 1) s[t+1] = [+1]; [] blockid.=0 blockid.=1 = =threadId. s[nt+] [0] [1] [] [] [3] [4] [5] t=

93 袖領域の処理 fblockid.> 0 && threadid.==0 ) s[t 1] = [ 1]; fblockid.< grddm. 1 && threadid.==blockdm. 1) s[t+1] = [+1]; [] blockid.=0 blockid.=1 = threadid.= 0 1 s[nt+] [0] [1] [] [] [3] [4] [5] t=

94 袖領域の処理 fblockid.> 0 && threadid.==0 ) s[t 1] = [ 1]; fblockid.< grddm. 1 && threadid.==blockdm. 1) s[t+1] = [+1]; [] blockid.=0 blockid.=1 = threadid.= s[nt+] [0] [1] [] [3] [] [3] [4] [5] t=

95 境界条件の処理 fblockid.==0 && threadid.==0 ) s[t 1] = 3.0*s[t] 3.0*s[t+1]+s[t+]; fblockid.==grddm. 1 && threadid.==blockdm. 1) s[t+1] = 3.0*s[t] 3.0*s[t 1]+s[t ]; blockid.=0 blockid.=1 s[nt+] [ 1] [0] [1] [] [3] [] [3] [4] [5] t= d d 3 1 境界での差分式と中心差分式が一致するように 1 を決定 Δ Δ 0 585

96 境界条件の処理 fblockid.==0 && threadid.==0 ) s[t 1] = 3.0*s[t] 3.0*s[t+1]+s[t+]; fblockid.==grddm. 1 && threadid.==blockdm. 1) s[t+1] = 3.0*s[t] 3.0*s[t 1]+s[t ]; blockid.=0 blockid.=1 s[nt+] [ 1] [0] [1] [] [3] [] [3] [4] [5] [6] t= N 3 N 1 3 N N 3 境界での差分式と中心差分式が一致するように N を決定 d d 3 N 1 4N N 3 N N Δ Δ N 1 586

97 中心差分の計算 syncthreads); dd[] = s[t+1] s[t 1])/.0*d); blockid.=0 blockid.=1 t= s[nt+] [ 1] [0] [1] [] [3] [] [3] [4] [5] [6] 1 Δ 全スレッドが同じ式 で中心差分を計算 dd[] 587