Microsoft PowerPoint - omp-f-01.ppt [互換モード]

Transcription

1 科学技術計算のための マルチコアプログラミング入門 Fortr 編第 Ⅰ 部 : 概要 対象アプリケーション OpeMP 中島研吾 東京大学情報基盤センター

2 OMP- 本編の背景 マイクロプロセッサのマルチコア化 メニーコア化 低消費電力 様々なプログラミングモデル OpeMP 指示行 ディレクティヴ を挿入するだけで手軽に 並列化 ができるため 広く使用されている 様々な解説書 データ依存性 t epeecy メモリへの書き込みと参照が同時に発生 並列化を実施するには 適切なデータの並べ替えを施す必要がある このような対策はOpeMP 向けの解説書でも詳しく取り上げられることは余りない : とても面倒くさい Hybr 並列プログラミングモデル

3 OMP- 3 Ft MPI vs. Hybr Ft-MPI:Ech PE -> Iepeet core core core memory core core memory core core memory core core core core core Hybr:Herrch Structure memory core core core core memory core core core core memory core core core core

4 4 We re ow Post-Pet-Sce Er GFLOPS PF PF TF PFLOPS: Pet = 5 Fotg OPertos per Sec. E-FLOPS = 8 w be tte -

5 5 Key-Issues towrs App./Agorthms o E-Sce Systems Jc Dogrr ORNL/U. Teessee t SIAM/PP Hybr/Heterogeeous Archtecture Mutcore + GPU Mutcore + Mycore more teget Me Precso Computto Auto-Tug/Sef-Aptg Fut Toert Commucto Reucg Agorthms

6 Heterogeeous Archtecture by 6 CPU+GPU or CPU+Mycore w be geer ess th 5 yers NVIDIA Ferm Ite MIC

7 7 CPU+Acceertor/Co-Processor GPU Mycore 高いメモリーバンド幅 GPU プログラミング環境 :CUDAOpeCL 一部のアプリケーションでは高効率 : 陽的 FDMBEM メニーコア Mycores Ite My Itegrte Core Archtecture MIC GPU より賢い : 軽い OS コンパイラが使える Ite eo Ph

8 8 Hybr 並列プログラミングモデルは必須 Messge Pssg MPI Mut Threg OpeMP CUDA OpeCL OpeACC 京 でも Hybr 並列プログラミングモデルが推奨されている 但し MPI+ 自動並列化 ノード内

9 OMP- 9 本編の目的 有限体積法から導かれる疎行列を対象としたICCG 法 を題材とした データ配置 reorergなど 科学技術計算のためのマルチコアプログラミングにおいて重要なアルゴリズムについての講習 更に理解を深めるための Fを利用した実習 単一のアプリケーションに特化した内容であるが 基本的な考え方は様々な分野に適用可能である 実はこの方法は意外に効果的である

10 OMP- 本編の目的 続き 単一のアプリケーションに特化した内容であるが 基本的な考え方は様々な分野に適用可能である 実はこの方法は意外に効果的である いわゆる 並列化講習会 とはだいぶ趣が異なる SMASH: Scece 無き科学技術計算はあり得ない! Scece Moeg Agorthm Softwre Hrwre

11 OMP- ファイルの用意 o your PC コピー 展開 >$ c <$CUR> >$ tr vf mutcore-c.tr >$ tr vf mutcore-f.tr >$ c mutcore 以下のディレクトリが出来ていることを確認 L L これらを以降 <$P-L><$P-L> Your PC F

12 OMP- 可視化には PrVew を使用 フリーソフトウェア Wows 版 Mc 版がある UNI 版もあり

13 OMP- 3 資料は Web 上にもあります

14 OMP- 4 背景 有限体積法 前処理付反復法 ICCG 法によるポアソン方程式法ソルバーについて 実行方法 データ構造 プログラムの説明 OpeMP 初期化 係数マトリクス生成 ICCG 法

15 OMP- 5 本編の目的より 有限体積法から導かれる疎行列を対象とした ICCG 法 を題材とした データ配置 reorerg など 科学技術計算のためのマルチコアプログラミングにおいて重要なアルゴリズムについての講習 有限体積法 疎行列 ICCG 法

16 OMP- 6 対象とするアプリケーションの概要 支配方程式 : 三次元ポアソン方程式 f y z 有限体積法 Fte Voume MethoFVM による空間離散化 任意形状の要素 要素中心で変数を定義 直接差分法 Drect Fte Dfferece Metho とも呼ばれる 境界条件 ディリクレ 体積フラックス 反復法による連立一次方程式解法 共役勾配法 CG+ 前処理

17 OMP- 7 解いている問題 : 三次元ポアソン方程式 変数 : 要素中心で定義 ポアソン方程式 y z f = t Z=Z m 境界条件 各要素で体積フラックス Z=Z m 面で= z y 各要素の体積フラックス = fot+j+ 要素体積 j=yzce3

18 OMP- 8 対象 : 規則正しい三次元差分格子半非構造的に扱う NZ Z z y N NY Y

19 OMP- 9 体積フラックス f の内容 f fot j y f z j YZ ce YZ ce YZ ce 3 YZce =3 は YZ 方向の差分格子のインデックス各メッシュが YZ 方向の何番目にあるかを示している NZ Z z y N NY Y

20 OMP- 有限体積法 Fte Voume Metho FVM 面を通過するフラックスの保存に着目 隣接要素との拡散 b S V Q 体積フラックス S e S b b c b V : 要素体積 S : 表面面積 j : 要素中心から表面までの距離 Q : 体積フラックス S c c c

21 OMP- 一次元ポアソン方程式 : 中央差分 / / Q - +

22 OMP- 有限体積法 Fte Voume Metho FVM 面を通過するフラックスの保存に着目 隣接要素との拡散 b S V Q 体積フラックス S e S b b c b V : 要素体積 S : 表面面積 j : 要素中心から表面までの距離 Q : 体積フラックス S c c c

23 OMP- 3 一次元差分法との比較 /3 h S b b S b b 一辺の長さhの正方形メッシュ接触面積 : S = h 要素体積 : V = h 接触面までの距離 : j =h/ 厚さ : この面を通過するフラックス :Qs b Qs b b b b S b フーリエの法則面を通過するフラックス = ー ポテンシャル勾配

24 OMP- 4 一次元差分法との比較 /3 h S b b S b b 一辺の長さhの正方形メッシュ接触面積 : S = h 要素体積 : V = h 接触面までの距離 : j =h/ 厚さ : S V Q 両辺を V で割る : この部分に注目すると V S Q

25 OMP- 5 一次元差分法との比較 3/3 h h h b b S S b b b b b h h h h S V h b b b h h h h h h h h 一辺の長さ h の正方形メッシュ接触面積 : S = h 要素体積 : V = h 接触面までの距離 : j =h/ 厚さ : 要素 について成立連立一次方程式

26 OMP- 6 三次元では NEIBcece NEIBcece NEIBcece3 NEIBcece5 NEIBcece4 NEIBcece6 NEIBcece NEIBcece NEIBcece3 NEIBcece5 NEIBcece4 NEIBcece6 y z y z z y f y z y z z y z y z y z y ce ce ce eb ce ce eb ce ce eb ce ce eb ce ce eb ce ce eb

27 OMP- 7 整理すると : 連立一次方程式 ce ce ce ce V f S z y f y z y z z y z y z y z y ce ce ce eb ce ce eb ce ce eb ce ce eb ce ce eb ce ce eb f A 対角項非対角項 N ce V f S S ce ce ce ce ce ce

28 8 有限要素法の係数マトリクスゼロが多い : 疎行列 } { } ]{ [ F F F F F F F F F F F F F F F F D D D D D D D D D D D D D D D D F K

29 OMP- 9 FVM の係数行列も疎行列面を通過するフラックスの保存に着目周囲の要素とのみ関係がある 隣接要素との拡散 b S V Q 体積フラックス S e S b b c b V : 要素体積 S : 表面面積 j : 要素中心から表面までの距離 Q : 体積フラックス S c c c

30 D-Prt 3 疎行列 : が多い Aj のように正方行列の全成分を記憶することは疎行列では非効率的 密 行列向け F F F F F F F F F F F F F F F F D D D D D D D D D D D D D D D D 有限体積法 : 非零非対角成分の数は高々 数百 規模 例えば未知数が 8 個あるとすると記憶容量 ワード数 は 正方行列 :O 6 非零非対角成分数 :O 非零成分のみ記憶するのが効率的

31 D-Prt 3 行列ベクトル積への適用 非零 非対角成分のみを格納 疎行列向け方法 Compresse Row Storge CRS Dg 対角成分 実数 =N Ie 非対角成分数に関する一次元配列 通し番号 整数 =N Item 非対角成分の要素 列 番号 整数 = en AMt 非対角成分 実数 = en {Y}= [A]{} o = N Y= Dg* o = Ie-+ Ie Y= Y + Amt*Item eo eo F F F F F F F F F F F F F F F F D D D D D D D D D D D D D D D D

32 行列ベクトル積への適用 非零 非対角成分のみを格納 疎行列向け方法 Compresse Row Storge CRS 3 {Q}=[A]{P} for=;<n;++{ W[Q][] = Dg[] * W[P][]; for=ie[];<ie[+];++{ W[Q][] += AMt[]*W[P][Item[]]; } } D-Prt

33 D-Prt 33 行列ベクトル積 : 密行列 とても簡単 N N N N N N N N N N N N N N N N N N N N y y y y {Y}= [A]{} o j= N Yj=. o = N Yj= Yj + Aj* eo eo

34 TK-FVM- 34 Compresse Row Storge CRS

35 TK-FVM- 35 Compresse Row Storge CRS NODE_tot= 8 対角成分 D=. D= 3.6 D3= 5.7 D4= 9.8 D5=.5 D6=.4 D7= 3. D8= 5.3

36 TK-FVM- 36 Compresse Row Storge CRS

37 TK-FVM- 37 Compresse Row Storge CRS 非対角 e= 成分数 e= 4 e= 6 e3= 8 3 e4= 4 e5= 5 e6= 7 4 e7= 4 e8= 5 e-+~e 番目が 行目の非対角成分 NPLU= 5 =en

38 TK-FVM- 38 Compresse Row Storge CRS 非対角 e= 成分数.4 3. e= 5 4 e= 6 e3= e4= e5= e6= e7= e8= e-+~e 番目が 行目の非対角成分 NPLU= 5 =en

39 TK-FVM- 39 Compresse Row Storge CRS 例 : tem 7= 5 AMAT 7=.5 tem9= 3 AMAT9=.5

40 TK-FVM- 4 Compresse Row Storge CRS D 対角成分 実数 =N e 非対角成分数に関する一次元配列 通し番号 整数 =N tem 非対角成分の要素 列 番号 整数 = en AMAT 非対角成分 実数 = en {Y}= [A]{} o = N Y= D* o = e-+ e Y= Y + AMAT*tem eo eo

41 4 疎行列 : 非零成分のみ記憶 メモリへの負担大 memory-bou: 間接参照 差分 FEMFVM {Y}= [A]{} o = N Y= D* o = e-+ e = tem Y= Y + AMAT* eo eo

42 4 行列ベクトル積 : 密行列 とても簡単メモリへの負担も小さい N N N N N N N N N N N N N N N N N N N N y y y y {Y}= [A]{} o j= N Yj=. o = N Yj= Yj + Aj* eo eo

43 OMP- 43 背景 有限体積法 前処理付反復法 ICCG 法によるポアソン方程式法ソルバーについて 実行方法 データ構造 プログラムの説明 OpeMP 初期化 係数マトリクス生成 ICCG 法

44 44 科学技術計算における 大規模線形方程式の解法 多くの科学技術計算は 最終的に大規模線形方程式 A=b を解くことに帰着される mportt epesve アプリケーションに応じて様々な手法が提案されている 疎行列 sprse 密行列 ese 直接法 rect 反復法 tertve 密行列 ese グローバルな相互作用 :BEM スペクトル法 MOMD 気液 疎行列 sprse ローカルな相互作用 :FEMFDMMD 固 高速多重極展開付 BEM

45 45 直接法 Drect Metho Guss の消去法 完全 LU 分解 逆行列 A - を直接求める 利点 安定 幅広いアプリケーションに適用可能 Prt Pvotg 疎行列 密行列いずれにも適用可能 欠点 反復法よりもメモリ 計算時間を必要とする 密行列の場合 ON 3 の計算量 大規模な計算向けではない ON の記憶容量 ON 3 の計算量

46 46 反復法とは 適当な初期解 から始めて 繰り返し計算によって真の解に収束 coverge させていく A b It Souto 初期解 Ler Equtos 連立一次方程式 b b b

47 47 反復法 Itertve Metho 定常 sttory 法 反復計算中 解ベクトル以外の変数は変化せず SORGuss-SeeJcobなど A b 概して遅い M 非定常 osttory 法 拘束 最適化条件が加わる Kryov 部分空間 subspce への写像を基底として使用するため Kryov 部分空間法とも呼ばれる CGCojugte Gret: 共役勾配法 BCGSTABB-Cojugte Gret Stbze GMRESGeerze Mm Resu Nb

48 48 反復法 Itertve Metho 続き 利点 直接法と比較して メモリ使用量 計算量が少ない 並列計算には適している 欠点 収束性が アプリケーション 境界条件の影響を受けやすい 前処理 precotog が重要

49 Sover-Itertve 49 非定常反復法 : クリロフ部分空間法 / Kryov Subspce Metho A I b b A 以下の反復式を導入し... を求める : r A b A I b where A b r : 残差ベクトル resu r A r I Ar r Ar A b A r b A b r

50 Sover-Itertve 5 非定常反復法 : クリロフ部分空間法 / Kryov Subspce Metho z は 次のクリロフ部分空間 Kryov Subspce に属するベクトル 問題はクリロフ部分空間からどのようにして解の近似ベクトル を求めるかにある : r A I I r A I r z r A I r A r I r r r A r A Ar r

51 Sover-Itertve 5 代表的な非定常反復法 : 共役勾配法 Cojugte Gret 法 略して CG 法 最も代表的な 非定常 反復法 対称正定値行列 Symmetrc Postve Defte:SPD 任意のベクトル {} に対して {} T [A]{}> 全対角成分 > 全固有値 > 全部分行列式 主小行列式 首座行列式 > と同値 アルゴリズム 最急降下法 Steepest Descet Metho の変種 = - + p : 反復解 p : 探索方向 : 定数 et 厳密解を y とするとき {-y} T [A]{-y} を最小とするような {} を求める 詳細は参考文献参照 例えば : 森正武 数値解析 第 版 共立出版

52 Sover-Itertve 5 共役勾配法 CG 法 のアルゴリズム Compute r = b-[a] for = - = r - r - f = p = r ese - = - / - p = r ef q = [A]p = - /p q = - + p r = r - - q chec covergece r e p - 行列ベクトル積 ベクトル内積 ベクトル定数倍の加減 DAPY : Vector : Scr

53 Sover-Itertve 53 共役勾配法 CG 法 のアルゴリズム Compute r = b-[a] for = - = r - r - f = p = r ese - = - / - p = r ef q = [A]p = - /p q = - + p r = r - - q chec covergece r e p - 行列ベクトル積 ベクトル内積 ベクトル定数倍の加減 DAPY : Vector : Scr

54 Sover-Itertve 54 共役勾配法 CG 法 のアルゴリズム Compute r = b-[a] for = - = r - r - f = p = r ese - = - / - p = r ef q = [A]p = - /p q = - + p r = r - - q chec covergece r e p - 行列ベクトル積 ベクトル内積 ベクトル定数倍の加減 DAPY : Vector : Scr

55 Sover-Itertve 55 共役勾配法 CG 法 のアルゴリズム Compute r = b-[a] for = - = r - r - f = p = r ese - = - / - p = r ef q = [A]p = - /p q = - + p r = r - - q chec covergece r e p - 行列ベクトル積 ベクトル内積 ベクトル定数倍の加減 DAPY Doube {y}= {} + {y} : Vector : Scr

56 Sover-Itertve 56 共役勾配法 CG 法 のアルゴリズム Compute r = b-[a] for = - = r - r - f = p = r ese - = - / - p = r ef q = [A]p = - /p q = - + p r = r - - q chec covergece r e p - : Vector : Scr

57 Sover-Itertve 57 CG 法アルゴリズムの導出 /5 y を厳密解 Ay=b とするとき 下式を最小にする を求める : T y A y T y A y A y A Ay y Ay A Ay y Ay A b y b 定数 従って 下記 f を最小にする を求めればよい : f A b f h f h A b h Ah 任意のベクトル h

58 58 b A f Ah h b A h f h f 任意のベクトル h Ah h b A h f Ah h b h A h b A b h b Ah h Ah A h A b h b Ah h A h b h h A h h f Sover-Itertve

59 59 CG 法アルゴリズムの導出 /5 p Sover-Itertve CG 法は任意の から始めて f の最小値を逐次探索する 今 番目の近似値 と探索方向 p が決まったとすると : f A b p Ap p p f f + を最小にするためには : Ap p r p Ap p A b p p f A b r は第 近似に対する残差

60 Sover-Itertve 6 CG 法アルゴリズムの導出 3/5 残差 r も以下の式によって計算できる : r r Ap r r b r A A r b A A Ap 探索方向を以下の漸化式によって求める : p r p r p 3 本当のところは下記のように + 回目に厳密解 y が求まれば良いのであるが 解がわかっていない場合は困難 y p

61 6 CG 法アルゴリズムの導出 4/5 ところで 下式のような都合の良い直交関係がある : 従って以下が成立する : Ap p p Ap y Ap Ap p r p Ap r p Ap A b p p A b p A b p A Ay p y Ap Ap p r p y Ap Sover-Itertve

62 Sover-Itertve 6 CG 法アルゴリズムの導出 5/5 p Ap r p Ap r Ap p Ap r Ap p Ap 4 p Ap p と p + が行列 A に関して共役 cojugte Compute p =r = b-[a] for = cc. - = p - r = r - - [A]p - chec covergece r f ot coverge p p r p r Ap Ap Ap e cc. - p = r + - p -

63 Sover-Itertve 63 CG 法アルゴリズム 任意の j に対して以下の共役関係が得られる : j p Ap j 探索方向 p 残差ベクトル r についても以下の関係が成立する : j r r j p r r r N 次元空間で互いに直交で一次独立な残差ベクトル r は N 個しか存在しない 従って共役勾配法は未知数が N 個のときに N 回以内に収束する 実際は丸め誤差の影響がある 条件数が大きい場合 Top Agorthms the th Cetury SIAM モンテカルロ法 シンプレックス法 クリロフ部分空間法 行列分解法 最適化 Fortr コンパイラ QR 法 クイックソート FFT 整数関係アルゴリズム FMM 高速多重極法

64 64 Proof /3 Mthemtc Iucto 数学的帰納法 Sover-Itertve j Ap p j r r j j Ap p r p 3 4 Ap r r p r p r p Ap p Ap r 直交性共役性

65 65 Proof /3 Mthemtc Iucto 数学的帰納法 Sover-Itertve * s stsfe for j j where Ap p Ap p Ap p p Ap r Ap r r r r r r j Ap p j r r j j f < f = Ap p r p Ap r p r p r r p r r r p r r Ap p r r Ap p p r r Ap r r r r r * * * * Ap p r p 3 4 Ap r r p r p Ap p Ap r j または

66 Sover-Itertve Proof 3/3 Mthemtc Iucto 数学的帰納法 j r r j j p Ap j * 66 * s stsfe for f < f = j where j または j 3 p Ap r p Ap r Ap * p r p Ap r r r r r Ap 3 p r p r Ap 4 p Ap 3 p Ap r Ap p Ap 4

67 67 Sover-Itertve r p r r Ap p r r Ap p p r r Ap r r r r r r r 3 * Ap p r p 3 4 Ap r r p r p Ap p Ap r r p r r

68 68 r r r r r Ap r r r r r Ap p Ap r 実際は はもうちょっと簡単な形に変形できる : r r r p Ap p r r Ap p r p Ap p A b p Sover-Itertve

69 Sover-Itertve 69 共役勾配法 CG 法 のアルゴリズム Compute r = b-[a] for = - = r - r - f = p = r ese - = - / - p = r ef q = [A]p = - /p q = - + p r = r - - q chec covergece r e p - : Vector : Scr r r r p r r r Ap

70 OMP- 7 前処理 precotog とは? 反復法の収束は係数行列の固有値分布に依存 固有値分布が少なく かつに近いほど収束が早い 単位行列 条件数 coto umber 対称正定 = 最大最小固有値比 条件数がに近いほど収束しやすい もとの係数行列 [A] に良く似た前処理行列 [M] を適用することによって固有値分布を改善する 前処理行列 [M] によって元の方程式 [A]{}={b} を [A ]{}={b } へと変換する ここで [A ]=[M] - [A] {b }=[M] - {b} である [A ]=[M] - [A] が単位行列に近ければ良いということになる [A ]=[A][M] - のように右からかけることもある 前処理 は密行列 疎行列ともに使用するが 普通は疎行列を対象にすることが多い

71 OMP- 7 前処理付共役勾配法 Precotoe Cojugte Gret Metho PCG Compute r = b-[a] for = sove [M]z - = r - - = r - z - f = p = z ese - = - / - p = z ef q = [A]p = - /p q = - + p r = r - - q chec covergece r e p - [M]= [M ][M ] [A ] =b [A ]=[M ] - [A][M ] - =[M ] b =[M ] - b p =>[M ]p r =>[M ] - r p = r p - [M ]p = [M ] - r [M ]p - p = [M ] - [M ] - r p - p = [M] - r p - - = [M] - r - r - / [M] - r - r - - = [M] - r - r - / p - [A]p -

72 OMP- 7 前処理付共役勾配法 Precotoe Cojugte Gret Metho PCG Compute r = b-[a] for = sove [M]z - = r - - = r - z - f = p = z ese - = - / - p = z ef q = [A]p = - /p q = - + p r = r - - q chec covergece r e p - 実際にやるべき計算は : z M r 近似逆行列 の計算が必要 : M A M A 究極の前処理 : 本当の逆行列 M A M A 対角スケーリング : 簡単 = 弱い M D M D

73 OMP- 73 対角スケーリング 点ヤコビ前処理 前処理行列として もとの行列の対角成分のみを取り出した行列を前処理行列 [M] とする 対角スケーリング 点ヤコビ pot-jcob 前処理 M D... D sove [M]z - = r - という場合に逆行列を簡単に求めることができる 簡単な問題では収束する D N... D N

74 OMP- 74 ILU IC 最もよく使用されている前処理 疎行列用 不完全 LU 分解 Icompete LU Fctorzto 不完全コレスキー分解 Icompete Choesy Fctorzto 対称行列 不完全な直接法 もとの行列が疎でも 逆行列は疎とは限らない f- もとの行列と同じ非ゼロパターン f- 無し を持っているのが ILUIC

75 LU 分解法 : 完全 LU 分解法 直接法の一種 逆行列を直接求める手法 逆行列 に相当するものを保存しておけるので 右辺が変わったときに計算時間を節約できる 逆行列を求める際にF- もとの行列ではであったところに値が入る が生じる LU fctorzto OMP- 75

76 不 完全 LU 分解法 ILU fctorzto Icompete LU fctorzto F- の発生を制限して 前処理に使う手法 不完全な逆行列 少し弱い直接法 F- を許さないとき :ILU OMP- 76

77 OMP- 77 LU 分解による連立一次方程式の解法 A が 行列のとき A を次式のように表すことを あるいは そのような L と U そのものを A の LU 分解という. u u u u u u u u u u LU A L:Lower trgur prt of mtr A U:Upper trgur prt of mtr A

78 OMP- 78 連立一次方程式の行列表現 元の連立一次方程式の一般形 b b b 行列表現 b b b A b b A

79 LU 分解を用いた A=b の解法 A LU となるAのLU 分解 LとUを求める. Ly b の解 yを求める. 簡単! 3 U y の解 を求める. 簡単! この が A b の解となる A LU Ly b OMP- 79

80 OMP- 8 Ly=b の解法 : 前進代入 b Ly b b b y y y b y y y b y y b y y b y y b y y b y b y 芋づる式に oe fter other 解が求まる.

81 OMP- 8 U=y の解法 : 後退代入 y U y y y u u u u u u y u u u y u u y u / / / u u y u u y u y j j 芋づる式に oe fter other 解が求まる.

82 OMP- 8 LU 分解の求め方 u u u u u u u u u u u u u u u u u u u u u u u u u u u u

83 OMP- 83 数値例 u u u u u u u u u u A 第 行 u u u u 第 列 / / / u u u u u u 第 行 u u u u u u u u u 第 列 u u u u

84 数値例 続き A u u u u u u u u u u 第 3 行 u u u u 3 4 u u u u 第 3 列 u u u u 第 4 行 第 4 列 u u u u u 行 列 行 列 の順に求める式を作っていく. OMP- 84

85 OMP- 85 数値例 続き 結局 A L U

86 OMP- 86 実例 :5 点差分

87 OMP- 87 実例 :5 点差分

88 OMP- 88 実例 : 係数マトリクス =

89 OMP- 89 実例 : 解 =

90 OMP- 9 完全 LU 分解したマトリクス もとのマトリクス LU 分解したマトリクス [L][U] 同時に表示 [L] 対角成分 = 省略 f- が生じている もともと だった成分が非ゼロになっている

91 OMP- 9 不完全 LU 分解したマトリクス f- 無し 不完全 LU 分解したマトリクス f- 無し [L][U] 同時に表示 [L] 対角成分 = 省略 完全 LU 分解したマトリクス [L][U] 同時に表示 [L] 対角成分 = 省略 f- が生じている もともと だった成分が非ゼロになっている

92 OMP- 9 解の比較 : ちょっと違う 不完全 LU 分解 完全 LU 分解

93 OMP- 93 ILU IC 前処理 F- を全く考慮しない 不完全な 分解 記憶容量 計算量削減 これを解くと 不完全な 解が得られるが 本来の解とそれほどずれているわけではない 問題に依存する

94 OMP- 94 前処理の分類 :Tre-off We Strog Pot Jcob Dgo Bocg ILU ILU ILU Guss Emto Smpe Esy to be Preze Chep Compcte Gob Depeecy Epesve

95 OMP- 95 背景 有限体積法 前処理付反復法 ICCG 法によるポアソン方程式法ソルバーについて 実行方法 データ構造 プログラムの説明 初期化 係数マトリクス生成 ICCG 法 OpeMP 超 入門

96 OMP- 96 対象とするアプリケーションの概要 支配方程式 : 三次元ポアソン方程式 有限体積法 Fte Voume MethoFVM による空間離散化 任意形状の要素 要素中心で変数を定義 直接差分法 Drect Fte Dfferece Metho とも呼ばれる 境界条件 y ディリクレ 体積フラックス 反復法による連立一次方程式解法 共役勾配法 CG+ 前処理 z f

97 OMP- 97 対象 : 規則正しい三次元差分格子半非構造的に扱う NZ Z z y N NY Y

98 OMP- 98 解いている問題 : 三次元ポアソン方程式 変数 : 要素中心で定義 ポアソン方程式 y z f = t Z=Z m 境界条件 各要素で体積フラックス Z=Z m 面で= z y 各要素の体積フラックス = fot+j+ 要素体積 j=yzce3

99 OMP- 99 有限体積法 Fte Voume Metho FVM 面を通過するフラックスの保存に着目 隣接要素との拡散 b S V Q 体積フラックス S e S b b c b V : 要素体積 S : 表面面積 j : 要素中心から表面までの距離 Q : 体積フラックス S c c c

100 OMP- プログラムの実行プログラム 必要ファイル 実行ディレクトリ :<$P-L>/ru mg メッシュジェネレータ mesh.t メッシュファイル L-so ポアソン方程式ソルバー test.p 結果ファイル PrVew 実際は mesh.t は使わない INPUT.DAT 制御ファイル

101 OMP- $> c <$P-L>/ru プログラムの実行コンパイル $> gfortr -O mg.f -o mg or cc O mg.c o mg $> s mg メッシュジェネレータ : mg mg $> c../src $> me $> s../ru/l-so L-so ポアソン方程式ソルバー : L-so gfortr が無い場合は g95 などをご利用ください Mefe も書き換える

102 OMP- プログラムの実行メッシュ生成 $> c../ru $>./mg 4 3 $> s mesh.t mesh.t 下図の NNYNZ を入力すると mesh.t が生成される NZ N NY

103 OMP mesh.t/5 re '' N NY NZ re '' ICELTOT o = ICELTOT re NEIBce = 6 YZj j= 3 eo

104 OMP mesh.t/5 NZ NY N YZ 方向の要素数 re '' N NY NZ re '' ICELTOT o = ICELTOT re NEIBce = 6 YZj j= 3 eo

105 OMP mesh.t3/5 要素数 :N NY NZ re '' N NY NZ re '' ICELTOT o = ICELTOT re NEIBce = 6 YZj j= 3 eo

106 OMP re '' N NY NZ re '' ICELTOT mesh.t4/5 隣接要素 :NEIBce o = ICELTOT re NEIBce = 6 YZj j= 3 eo z y :4 3: 項目目は通し番号です 読み飛ばし 6: 5 5:3 4:9 :6

107 OMP- 7 NEIBce: 隣接している要素番号 境界面の場合は NEIBcece6 NEIBcece4 NEIBcece NEIBcece NEIBcece3 NEIBcece5 NEIBcece= ce NEIBcece= ce + NEIBcece3= ce N NEIBcece4= ce + N NEIBcece5= ce N*NY NEIBcece6= ce + N*NY

108 OMP mesh.t5/5 YZ 方向の位置 :YZj z y re '' N NY NZ re '' ICELTOT o = ICELTOT re NEIBce = 6 YZj j= 3 eo

109 OMP- 9 NEIBce: 隣接している要素番号 境界面の場合は NEIBcece6 NEIBcece4 NEIBcece NEIBcece NEIBcece3 = YZce j= YZce = YZce3 ce= -*N*NY + j-*n + z y NEIBcece5 NEIBcece= ce NEIBcece= ce + NEIBcece3= ce N NEIBcece4= ce + N NEIBcece5= ce N*NY NEIBcece6= ce + N*NY

110 OMP- プログラムの実行制御データ <$P-L>/ru/INPUT.DAT の作成 N/NY/NZ MEHOD :.e-.e-.e- D/DY/DZ.e-8 EPSICCG N NY NZ 各方向のメッシュ数 METHOD 前処理行列の作成方法 : 次ページ D DY DZ 各要素の YZ 方向辺長さ EPSICCG ICCG 法の収束判定値 NZ N NY

111 OMP- 前処理手法の選択 N/NY/NZ MEHOD :.e-.e-.e- D/DY/DZ.e-8 EPSICCG METHOD= 不完全修正コレスキー分解 非対角項保存 METHOD= 不完全修正コレスキー分解 METHOD=3 対角スケーリング 点ヤコビ METHOD=3 について計算してみよ!

112 OMP- $> c <$P-L>/ru $>./L-so $> s test.p test.p プログラムの実行計算実行 ポスト処理

113 FEM3D-Prt 3 PrVew ファイルを開く 図の表示 イメージファイルの保存

114 4 UCD Formt /3 Ustructure Ce Dt 要素の種類点線三角形四角形四面体角錐三角柱六面体二次要素線 三角形 四角形 四面体 角錐 三角柱 六面体 キーワード pt e tr qu tet pyr prsm he e tr qu tet pyr prsm he 4

115 5 UCD Formt /3 Orgy for AVS mcroavs Eteso of the UCD fe s p There re two types of formts. Oy o type c be re by PrVew. 5

116 6 UCD Formt 3/3: O Formt 全節点数 全要素数 各節点のデータ数 各要素のデータ数 モデルのデータ数 節点番号 座標 Y 座標 Z 座標 節点番号 座標 Y 座標 Z 座標 要素番号 材料番号 要素の種類 要素を構成する節点のつながり 要素番号 材料番号 要素の種類 要素を構成する節点のつながり 要素のデータ成分数 成分 の構成数 成分 の構成数 各成分の構成数 要素データ成分 のラベル 単位 要素データ成分 のラベル 単位 各要素データ成分のラベル 単位 要素番号 要素データ 要素データ 要素番号 要素データ 要素データ 節点のデータ成分数 成分 の構成数 成分 の構成数 各成分の構成数 節点データ成分 のラベル 単位 節点データ成分 のラベル 単位 各節点データ成分のラベル 単位 節点番号 節点データ 節点データ 節点番号 節点データ 節点データ 6

117 OMP- 7 背景 有限体積法 前処理付反復法 ICCG 法によるポアソン方程式法ソルバーについて 実行方法 データ構造 プログラムの説明 OpeMP 初期化 係数マトリクス生成 ICCG 法

118 OMP- 8 プログラムの構成 progrm MAIN use STRUCT use PCG use sover_iccg use sover_iccg use sover_pcg mpct REAL*8 A-HO-Z c INPUT c POINTER_INIT c BOUNDARY_CELL c CELL_METRICS c POI_GEN PHI=. f METHOD.eq. c sove_iccg f METHOD.eq. c sove_iccg f METHOD.eq.3 c sove_pcg c OUTUCD stop e MAIN メインルーチン INPUT 制御ファイル読込 INPUT.DAT POINTER_INIT メッシュファイル読込 mesh.t BOUNDARY_CELL = を設定する要素の探索 CELL_METRICS 表面積 体積等の計算 POI_GEN 行列コネクティビティ生成 各成分の計算 境界条件 SOLVER_ICCG ICCG 法ソルバー METHOD= SOLVER_ICCG ICCG 法ソルバー METHOD= j FORM_ILU ~ j j ~ ~ ~ j SOLVER_PCG ICCG 法ソルバー METHOD=3

119 OMP- 9 moue STRUCT moue STRUCT cue 'precso.c'!c!c-- METRICs & FLU teger =t :: ICELTOT ICELTOTp N teger =t :: N NY NZ NP NYP NZP IBNODTOT teger =t :: Nc NYc NZc re =re :: & & D DY DZ AREA YAREA ZAREA RD RDY RDZ & & RD RDY RDZ RD RDY RDZ re =re meso: octbe :: & & VOLCEL VOLNOD RVC RVN teger =t meso:: octbe :: & & YZ NEIBce ICELTOT: 要素数 N NY NZ N: 節点数 NNYNZ: yz 方向要素数 NPNYPNZP: yz 方向節点数 IBNODTOT: NP NYP YZICELTOT3: 要素座標 NEIBceICELTOT6: 隣接要素!C!C-- BOUNDARYs teger =t :: ZmCELtot teger =t meso: octbe :: BC_INDE BC_NOD teger =t meso: octbe :: ZmCEL!C!C-- WORK teger =t meso:: octbe :: IWK re=re meso:: octbe :: FCV e moue STRUCT

120 OMP- moue PCG/5 moue PCG teger prmeter :: N= 56 teger :: NUm NLm NCOLORtot NCOLOR NU NL METHOD re=8 :: EPSICCG re=8 meso: octbe :: D PHI BFORCE re=8 meso: octbe :: AL AU teger meso: octbe :: INL INU COLORe teger meso: octbe :: OLDtoNEW NEWtoOLD teger meso:: octbe :: IAL IAU teger meso: octbe :: el teml teger meso: octbe :: eu temu e moue PCG 扱う行列 : 疎行列 自分の周辺のみと接続 非ゼロ成分のみを記憶する 上下三角成分を別々に記憶 L U

121 OMP- moue PCG/5 moue PCG teger prmeter :: N= 56 teger :: NUm NLm NCOLORtot NCOLOR NU NL METHOD teger :: NPL NPU re=8 :: EPSICCG re=8 meso: octbe :: D PHI BFORCE re=8 meso: octbe :: AL AU teger meso: octbe :: INL INU COLORe teger meso: octbe :: OLDtoNEW NEWtoOLD teger meso:: octbe :: IAL IAU teger meso: octbe :: el teml teger meso: octbe :: eu temu e moue PCG 補助配列 下三角成分 列番号 : 非対角成分で自分より要素番号が小さい IALcou < 上三角成分 列番号 : 非対角成分で自分より要素番号が大きい IAUcou > INLICELTOT 非零下三角成分の数 IALNLICELTOT 非零下三角成分 列番号 INUICELTOT 非零上三角成分の数 IAUNUICELTOT 非零上三角成分 列番号 NUNL 非零上下三角成分の最大数 ここでは6 U el:iceltot eu:iceltot NPLNPU temlnpltemunpu 各行の非零下三角成分数 CRS 各行の非零上三角成分数 CRS 非零上下三角成分数の合計 CRS 比零上下三角成分 列番号 CRS L

122 OMP- moue PCG3/5 moue PCG teger prmeter :: N= 56 teger :: NUm NLm NCOLORtot NCOLOR NU NL METHOD teger :: NPL NPU re=8 :: EPSICCG re=8 meso: octbe :: D PHI BFORCE re=8 meso: octbe :: AL AU teger meso: octbe :: INL INU COLORe teger meso: octbe :: OLDtoNEW NEWtoOLD teger meso:: octbe :: IAL IAU teger meso: octbe :: el teml teger meso: octbe :: eu temu 補助配列 下三角成分 列番号 : IALcou < その個数が INL 上三角成分 列番号 : IAUcou > その個数が INU e moue PCG INLICELTOT 非零下三角成分の数 IALNLICELTOT 非零下三角成分 列番号 INUICELTOT 非零上三角成分の数 IAUNUICELTOT 非零上三角成分 列番号 NUNL 非零上下三角成分の最大数 ここでは6 U el:iceltot eu:iceltot NPLNPU temlnpltemunpu 各行の非零下三角成分数 CRS 各行の非零上三角成分数 CRS 非零上下三角成分数の合計 CRS 比零上下三角成分 列番号 CRS L

123 OMP- 3 moue PCG4/5 moue PCG teger prmeter :: N= 56 teger :: NUm NLm NCOLORtot NCOLOR NU NL METHOD teger :: NPL NPU re=8 :: EPSICCG re=8 meso: octbe :: D PHI BFORCE re=8 meso: octbe :: AL AU teger meso: octbe :: INL INU COLORe teger meso: octbe :: OLDtoNEW NEWtoOLD teger meso:: octbe :: IAL IAU teger meso: octbe :: el teml teger meso: octbe :: eu temu 補助配列 下三角成分 列番号 : IALcou < その個数が INL 上三角成分 列番号 : IAUcou > その個数が INU e moue PCG INLICELTOT 非零下三角成分の数 IALNLICELTOT 非零下三角成分 列番号 INUICELTOT 非零上三角成分の数 IAUNUICELTOT 非零上三角成分 列番号 NUNL 非零上下三角成分の最大数 ここでは6 U el:iceltot eu:iceltot NPLNPU temlnpltemunpu 各行の非零下三角成分数 CRS 各行の非零上三角成分数 CRS 非零上下三角成分数の合計 CRS 比零上下三角成分 列番号 CRS L

124 OMP- 4 moue PCG 5/5 moue PCG teger prmeter :: N= 56 teger :: NUm NLm NCOLORtot NCOLOR NU NL METHOD teger :: NPL NPU re=8 :: EPSICCG re=8 meso: octbe :: D PHI BFORCE re=8 meso: octbe :: AL AU teger meso: octbe :: INL INU COLORe teger meso: octbe :: OLDtoNEW NEWtoOLD teger meso:: octbe :: IAL IAU teger meso: octbe :: el teml teger meso: octbe :: eu temu e moue PCG METHOD 前処理手法 ===3 EPSICCG D ICELTOT PHI ICLETOT BFORCEICELTOT ALNPLAUNPU 収束打切残差係数行列の対角成分従属変数連立一次方程式の右辺ベクトル係数行列の比零上下三角成分 CRS

125 OMP- 5 行列関係変数 : まとめ 配列 変数名 型 内容 DN R 対角成分 N: 全メッシュ数 =ICELTOT BFORCEN R 右辺ベクトル PHIN R 未知数ベクトル el:n I 各行の非零下三角成分数 CRS eu:n I 各行の非零上三角成分数 CRS NPL I 非零下三角成分総数 CRS NPU I 非零上三角成分総数 CRS temlnpl I 非零下三角成分 列番号 CRS temunpu I 非零下三角成分 列番号 CRS ALNPL R 非零下三角成分 係数 CRS AUNPL R 非零上三角成分 係数 CRS

126 OMP- 6 行列関係変数 : まとめ 補助配列 配列 変数名 型 内容 NLNU I 各行の非零上下三角成分の最大数 ここでは6 INLN I 各行の非零下三角成分数 INUN I 各行の非零上三角成分数 IALNLN I 各行の非零下三角成分に対応する列番号 IAUNUN I 各行の非零上三角成分に対応する列番号 補助配列を使う理由 NPLNPUの値が前以てわからない後掲の並び替え orergreorerg のとき CRS 形式ではやりにくい

127 OMP- 7 行列ベクトル積 :{q}=[a]{p} o = N VAL= D*p o = el-+ el VAL= VAL + AL*ptemL eo o = eu-+ eu VAL= VAL + AU*ptemU eo q= VAL eo

128 OMP- 8 プログラムの構成 progrm MAIN use STRUCT use PCG use sover_iccg use sover_iccg use sover_pcg mpct REAL*8 A-HO-Z c INPUT c POINTER_INIT c BOUNDARY_CELL c CELL_METRICS c POI_GEN MAIN メインルーチン INPUT 制御ファイル読込 INPUT.DAT POINTER_INIT メッシュファイル読込 mesh.t BOUNDARY_CELL = を設定する要素の探索 CELL_METRICS 表面積 体積等の計算 PHI=. f METHOD.eq. c sove_iccg f METHOD.eq. c sove_iccg f METHOD.eq.3 c sove_pcg c OUTUCD stop e POI_GEN 行列コネクティビティ生成 各成分の計算 境界条件 SOLVER_ICCG ICCG 法ソルバー METHOD= SOLVER_ICCG ICCG 法ソルバー METHOD= SOLVER_PCG ICCG 法ソルバー METHOD=3 j FORM_ILU ~ j j ~ ~ ~ j

129 OMP- 9 put: INPUT.DAT の読み込み!C!C***!C*** INPUT!C***!C!C INPUT CONTROL DATA!C subroute INPUT use STRUCT use PCG mpct REAL*8 A-HO-Z chrcter*8 CNTFIL!C!C-- CNTL. fe ope fe='input.dat' sttus='uow' re * N NY NZ re * METHOD re * D DY DZ re * EPSICCG cose!c=== retur e N/NY/NZ MEHOD -3.e-.e-.e- D/DY/DZ.e-8 EPSICCG

130 OMP- 3 poter_t/3: mesh.t の作成!C!C***!C*** POINTER_INIT!C***!C subroute POINTER_INIT use STRUCT use PCG mpct REAL*8 A-HO-Z!C!C !C Geertg MESH fo.!c !C=== ICELTOT= N * NY * NZ NP= N + NYP= NY + NZP= NZ + octe NEIBceICELTOT6 YZICELTOT3 NEIBce= NNYNZ: yz 方向要素数 NPNYPNZP: yz 節点数 可視化用 ICELTOT: 要素数 N NY NZ NEIBceICELTOT6: 隣接要素 octe/eocte for C YZICELTOT3: 要素座標

131 OMP- 3 poter_t/3: mesh.t の作成 o = NZ o j= NY o = N ce= -*N*NY + j-*n + NEIBcece= ce - NEIBcece= ce + NEIBcece3= ce - N NEIBcece4= ce + N NEIBcece5= ce - N*NY NEIBcece6= ce + N*NY f.eq. NEIBcece= f.eq.n NEIBcece= f j.eq. NEIBcece3= f j.eq.ny NEIBcece4= f.eq. NEIBcece5= f.eq.nz NEIBcece6= YZce= YZce= j YZce3= NEIBcece6 NEIBcece NEIBcece3 NEIBcece5 NEIBcece4 NEIBcece!C=== eo eo eo = YZce j= YZce = YZce3 ce= -*N*NY + j-*n + NEIBcece= ce NEIBcece= ce + NEIBcece3= ce N NEIBcece4= ce + N NEIBcece5= ce N*NY NEIBcece6= ce + N*NY

132 OMP- 3 poter_t3/3!c!c !C Prmeters!C !C=== f D.e..e the D=. / fotn DY=. / fotny DZ=. / fotnz ef D=. となっていた場合のみ DDYDZ をこのように指定 NP= N + NYP= NY + NZP= NZ +!C=== IBNODTOT= NP * NYP N = NP * NYP * NZP retur e

133 OMP- 33 poter_t3/3: 節点 可視化用!C!C !C Prmeters!C !C=== f D.e..e the D=. / fotn DY=. / fotny DZ=. / fotnz ef NP= N + NYP= NY + NZP= NZ !C=== IBNODTOT= NP * NYP N = NP * NYP * NZP retur e NPNYPNZP: yz 方向節点数 IBNODTOT: NP NYP N: 節点数 可視化に使用

134 OMP- 34 boury_ce!c!c***!c*** BOUNDARY_CELL!C***!C subroute BOUNDARY_CELL use STRUCT mpct REAL*8 A-HO-Z Z=Z m の要素の定義 総数 : ZmCELtot 要素番号 : ZmCEL:!C!C !C Zm!C !C=== IFACTOT= N * NY ZmCELtot= IFACTOT octe ZmCELZmCELtot!C=== cou= = NZ o j= NY o = N ce= -*IFACTOT + j-*n + cou= cou + ZmCELcou= ce eo eo retur e Z

135 OMP- 35 ce_metrcs!c!c***!c*** CELL_METRICS!C***!C subroute CELL_METRICS use STRUCT use PCG mpct REAL*8 A-HO-Z!C!C-- ALLOCATE octe VOLCELICELTOT octe RVCICELTOT!C!C-- VOLUME AREA PROJECTION etc. AREA= DY * DZ YAREA= D * DZ ZAREA= D * DY RD=. / D RDY=. / DY RDZ=. / DZ RD=. / D** RDY=. / DY** RDZ=. / DZ** RD=. /.5*D RDY=. /.5*DY RDZ=. /.5*DZ 計算に必要な諸パラメータ DZ z y AREA DY D AREA Y Z YAREA Z ZAREA Y RD RDY RDZ Y Z V= D * DY * DZ RV=./V VOLCEL= V RVC = RV retur e

136 OMP- 36 ce_metrcs!c!c***!c*** CELL_METRICS!C***!C subroute CELL_METRICS use STRUCT use PCG mpct REAL*8 A-HO-Z!C!C-- ALLOCATE octe VOLCELICELTOT octe RVCICELTOT!C!C-- VOLUME AREA PROJECTION etc. AREA= DY * DZ YAREA= D * DZ ZAREA= D * DY RD=. / D RDY=. / DY RDZ=. / DZ RD=. / D** RDY=. / DY** RDZ=. / DZ** RD=. /.5*D RDY=. /.5*DY RDZ=. /.5*DZ V= D * DY * DZ RV=./V VOLCEL= V RVC = RV retur e 計算に必要な諸パラメータ DZ z y RD AREA DY D RDY Y RDZ RD RDY.5.5 Y RDZ.5 Z Z

137 OMP- 37 ce_metrcs!c!c***!c*** CELL_METRICS!C***!C subroute CELL_METRICS use STRUCT use PCG mpct REAL*8 A-HO-Z!C!C-- ALLOCATE octe VOLCELICELTOT octe RVCICELTOT!C!C-- VOLUME AREA PROJECTION etc. AREA= DY * DZ YAREA= D * DZ ZAREA= D * DY RD=. / D RDY=. / DY RDZ=. / DZ 計算に必要な諸パラメータ DZ z y AREA DY D RD=. / D** RDY=. / DY** RDZ=. / DZ** RD=. /.5*D RDY=. /.5*DY RDZ=. /.5*DZ V= D * DY * DZ RV=./V VOLCEL= V RVC = RV retur e VOLCEL V Y RV RVC VOLCEL Z

138 OMP- 38 プログラムの構成 progrm MAIN use STRUCT use PCG use sover_iccg use sover_iccg use sover_pcg mpct REAL*8 A-HO-Z c INPUT c POINTER_INIT c BOUNDARY_CELL c CELL_METRICS c POI_GEN MAIN メインルーチン INPUT 制御ファイル読込 INPUT.DAT POINTER_INIT メッシュファイル読込 mesh.t BOUNDARY_CELL = を設定する要素の探索 CELL_METRICS 表面積 体積等の計算 PHI=. f METHOD.eq. c sove_iccg f METHOD.eq. c sove_iccg f METHOD.eq.3 c sove_pcg c OUTUCD stop e POI_GEN 行列コネクティビティ生成 各成分の計算 境界条件 SOLVER_ICCG ICCG 法ソルバー METHOD= SOLVER_ICCG ICCG 法ソルバー METHOD= SOLVER_PCG ICCG 法ソルバー METHOD=3 j FORM_ILU ~ j j ~ ~ ~ j

139 OMP- 39 po_ge/7 subroute POI_GEN use STRUCT use PCG mpct REAL*8 A-HO-Z!C!C-- INIT. = ICELTOT NU= 6 NL= 6 octe BFORCE D PHI octe INL INU IALNL IAUNU PHI=. D=. INL= INU= IAL= IAU=

140 OMP- 4 配列の宣言 配列 変数名 型 内容 DN R 対角成分 N: 全メッシュ数 =ICELTOT BFORCEN R 右辺ベクトル PHIN R 未知数ベクトル el:n I 各行の非零下三角成分数 CRS eu:n I 各行の非零上三角成分数 CRS NPL I 非零下三角成分総数 CRS NPU I 非零上三角成分総数 CRS temlnpl I 非零下三角成分 列番号 CRS temunpu I 非零下三角成分 列番号 CRS ALNPL R 非零下三角成分 係数 CRS AUNPL R 非零上三角成分 係数 CRS

141 OMP- 4!C!C !C CONNECTIVITY!C !C=== o ce= ICELTOT cn= NEIBcece cn= NEIBcece cn3= NEIBcece3 cn4= NEIBcece4 cn5= NEIBcece5 cn6= NEIBcece6 coug= f cn5.e. the cou= INLce + IALcouce= cn5 INL ce= cou ef f cn3.e. the cou= INLce + IALcouce= cn3 INL ce= cou ef f cn.e. the cou= INLce + IALcouce= cn INL ce= cou ef NEIBcece po_ge/7 NEIBcece6 NEIBcece3 NEIBcece5 NEIBcece4 下三角成分 NEIBcece5= ce N*NY NEIBcece3= ce N NEIBcece= ce NEIBcece

142 OMP- 4 f cn.e. the cou= INUce + IAUcouce= cn INU ce= cou ef po_ge3/7 f cn4.e. the cou= INUce + IAUcouce= cn4 INU ce= cou ef NEIBcece6 NEIBcece4!C=== f cn6.e. the cou= INUce + IAUcouce= cn6 INU ce= cou ef eo NEIBcece NEIBcece3 NEIBcece NEIBcece5 上三角成分 NEIBcece= ce + NEIBcece4= ce + N NEIBcece6= ce + N*NY

143 OMP- 43!C!C-- D rry octe el: eu: el= eu= o ce= ICELTOT elce= INLce euce= INUce eo o ce= ICELTOT elce= elce + elce- euce= euce + euce- eo NPL= eliceltot NPU= euiceltot octe temlnpl ALNPL octe temunpu AUNPU teml= temu= po_ge4/7 配列の宣言 配列 変数名 型 内容 DN R 対角成分 N: 全メッシュ数 =ICELTOT BFORCEN R 右辺ベクトル PHIN R 未知数ベクトル el:n I 各行の非零下三角成分数 CRS eu:n I 各行の非零上三角成分数 CRS NPL I 非零下三角成分総数 CRS NPU I 非零上三角成分総数 CRS temlnpl I 非零下三角成分 列番号 CRS temunpu I 非零下三角成分 列番号 CRS ALNPL R 非零下三角成分 係数 CRS AUNPL R 非零上三角成分 係数 CRS!C=== AL=. AU=.

144 OMP- 44 有限体積法 Fte Voume Metho FVM 面を通過するフラックスの保存に着目 隣接要素との拡散 b S V Q 体積フラックス S e S b b c b V : 要素体積 S : 表面面積 j : 要素中心から表面までの距離 Q : 体積フラックス S c c c

145 OMP- 45 全体マトリクスの生成要素 に関する釣り合い S V Q S S V Q S S V Q D 対角成分 AL AU 非対角成分 BFORCE 右辺

146 OMP- 46!C!C !C INTERIOR & NEUMANN BOUNDARY CELLs!C !C=== coug= o ce= ICELTOT cn= NEIBcece cn= NEIBcece cn3= NEIBcece3 cn4= NEIBcece4 cn5= NEIBcece5 cn6= NEIBcece6 VOL= VOLCELce cou= f cn5.e. the coef =RDZ * ZAREA Dce= Dce - coef cou= cou + = cou + elce- teml= cn5 AL= coef ef f cn3.e. the coef =RDY * YAREA Dce= Dce - coef cou= cou + = cou + elce- teml= cn3 AL= coef ef f cn.e. the coef =RD * AREA Dce= Dce - coef cou= cou + = cou + elce- teml= cn AL= coef ef po_ge5/7 係数の計算 境界面以外 W D N DY DY D E W y y N S y y S f c E y

147 OMP- 47 三次元では NEIBcece6 NEIBcece NEIBcece4 NEIBcece f cn5.e. the coef = RDZ * ZAREA Dce= Dce coef z y NEIBcece3 cou= cou + = cou + elce- NEIBcece5 eb ce eb ce3 y eb ce5 z ce ce ce yz z y eb ce eb ce4 y eb ce6 z ce ce ce yz z y teml= cn5 AL= coef ef f ce yz

148 OMP- 48 po_ge6/7 cou= f cn.e. the coef = RD * AREA Dce= Dce - coef cou= cou + = cou + euce- temu= cn AU= coef ef 係数の計算 境界面以外 N f cn4.e. the coef = RDY * YAREA Dce= Dce - coef cou= cou + = cou + euce- temu= cn4 AU= coef ef W D DY D E!C=== f cn6.e. the coef = RDZ * ZAREA Dce= Dce - coef cou= cou + = cou + euce- temu= cn6 AU= coef ef = YZce jj= YZce = YZce3 BFORCEce= -fot+jj+ * VOL eo DY S E W y y N S y y f c y

149 OMP- 49 po_ge6/7 cou= f cn.e. the coef = RD * AREA Dce= Dce - coef cou= cou + = cou + euce- temu= cn AU= coef ef f cn4.e. the coef = RDY * YAREA Dce= Dce - coef cou= cou + = cou + euce- temu= cn4 AU= coef ef f cn6.e. the coef = RDZ * ZAREA Dce= Dce - coef cou= cou + = cou + euce- temu= cn6 AU= coef ef = YZce jj= YZce = YZce3 体積発熱 f fot j j YZ ce YZ ce YZ ce 3 YZce =3 は YZ 方向の差分格子のインデックス 各メッシュが YZ 方向の何番目にあるかを示している!C=== BFORCEce= -fot+jj+ * VOL eo

150 OMP- 5!C!C !C DIRICHLET BOUNDARY CELLs!C !C TOP SURFACE!C=== o b= ZmCELtot ce= ZmCELb coef=. * RDZ * ZAREA Dce= Dce - coef eo!c=== retur e po_ge7/7 係数の計算 境界面 =- DZ Z=Zm = 境界面の外側に 大きさが同じで 値が =- となるような要素があると仮定 境界面で丁度 = となる : 一次近似

151 OMP- 5 ディリクレ条件 S V Q W D DZ S D E S D 対角成分 AL AU BFORCE 非対角成分 右辺 S V Q

152 OMP- 5 ディリクレ条件 D D DZ DZ W E S N Q V S D 対角成分 BFORCE 右辺 AL AU 非対角成分 Q V S S N N Q V y z S S

153 OMP- 53 ディリクレ条件 D D DZ DZ W E S N Q V S D 対角成分 BFORCE 右辺 AL AU 非対角成分 Q V S S N N Q V y z S S Q V y z S S

154 OMP- 54 ディリクレ条件 D D DZ DZ W E S N Q V S D 対角成分 BFORCE 右辺 AL AU 非対角成分 Q V S S Q V y z S S

155 OMP- 55 ディリクレ条件 D D DZ DZ W E S N Q V S D 対角成分 BFORCE 右辺 AL AU 非対角成分 Q V S S Q V y z S S Q V S y z S

156 OMP- 56 ディリクレ条件 W D N DZ DZ S D E S S S y z S V Q D 対角成分 AL AU BFORCE 非対角成分 右辺 o b= ZmCELtot ce= ZmCELb coef=. * RDZ * y V Q ZAREA Dce= Dce - coef eo z S S VQ S V Q

157 OMP- 57 プログラムの構成 progrm MAIN use STRUCT use PCG use sover_iccg use sover_iccg use sover_pcg mpct REAL*8 A-HO-Z c INPUT c POINTER_INIT c BOUNDARY_CELL c CELL_METRICS c POI_GEN MAIN メインルーチン INPUT 制御ファイル読込 INPUT.DAT POINTER_INIT メッシュファイル読込 mesh.t BOUNDARY_CELL = を設定する要素の探索 CELL_METRICS 表面積 体積等の計算 PHI=. f METHOD.eq. c sove_iccg f METHOD.eq. c sove_iccg f METHOD.eq.3 c sove_pcg c OUTUCD stop e POI_GEN 行列コネクティビティ生成 各成分の計算 境界条件 SOLVER_ICCG ICCG 法ソルバー METHOD= SOLVER_ICCG ICCG 法ソルバー METHOD= SOLVER_PCG ICCG 法ソルバー METHOD=3 j FORM_ILU ~ j j ~ ~ ~ j

158 OMP- 58 背景 有限体積法 前処理付反復法 ICCG 法によるポアソン方程式法ソルバーについて 実行方法 データ構造 プログラムの説明 OpeMP 初期化 係数マトリクス生成 ICCG 法

159 OMP- 59 あとは線形方程式を解けば良い 共役勾配法 Cojugte GretCG 前処理 不完全コレスキー分解 Icompete Choesy FctorztoIC 実は不完全 修正 コレスキー分解 ICCG

160 OMP- 6 修正コレスキー分解 対称行列 A の LU 分解 対称行列 A は 対角行列 D を利用して [A]= [L][D][L] T のような形に分解することができる この分解を LDL T 分解または修正コレスキー分解 mofe Choesy ecomposto と呼ぶ [A]= [L][L] T とするような分解法もある コレスキー分解 N=5 の場合の例

161 OMP- 6 修正コレスキー分解 j ここで j j j とすると以下が導かれる j j j j j

162 j j j j

163 OMP- 63 不完全修正コレスキー分解 j ここで j j j とすると以下が導かれる j j j j j 実際には 不完全 な分解を実施し このような形を用いることが多い j j

164 OMP- 64 プログラムの実行制御データ <$P-L>/ru/INPUT.DAT の作成 N/NY/NZ MEHOD :.e-.e-.e- D/DY/DZ..e-8 OMEGA EPSICCG METHOD 前処理行列の作成方法 : 不完全修正コレスキー分解 METHOD= 対角成分のみ変更 METHOD= 非対角成分変更 F-は無し: j の場合のみ METHOD=3 対角スケーリング 点ヤコビ j j j j j

165 OMP- 65 不完全修正コレスキー分解 j ここで j j j とすると以下が導かれる j j j j j j j 対角成分のみがもとの行列と変わる

166 OMP- 66 不完全修正コレスキー分解を使用した前進後退代入 r LDL z T r z LDL z M T r y L y z DL T / / / / / / / / / /

167 OMP- 67 不完全修正コレスキー分解を使用した 前進後退代入!C!C !C {z}= [Mv]{r}!C !C=== o = N WY= WR eo W DD / Ly r o = N WVAL= WY o = el-+ el WVAL= WVAL - AL * WtemLY eo WY= WVAL * WDD eo DL T z y!c=== o = N - SW =. o = eu-+ eu SW= SW + AU * WtemUZ eo WZ= WY - WDD * SW eo / 3 3 / / / / / / / / / 33 44

168 OMP- 68 不完全修正コレスキー分解を使用した 前進後退代入 : 計算順序考慮!C!C !C {z}= [Mv]{r}!C !C=== o = N WZ= WR eo W DD / Lz z o = N WVAL= WZ o = el-+ el WVAL= WVAL - AL * WtemLZ eo WZ= WVAL * WDD eo DL T z z!c=== o = N - SW =. o = eu-+ eu SW= SW + AU * WtemUZ eo WZ= WZ - WDD * SW eo / 3 3 / / / / / / / / / 33 44

169 OMP- 69 sove_iccg/7:method=!c***!c*** moue sover_iccg!c***! moue sover_iccg cots!c!c*** sove_iccg!c subroute sove_iccg & & N NPL NPU el teml eu temu D B & & AL AU EPS ITR IER mpct REAL*8 A-HO-Z re=8 meson :: D re=8 meson :: B re=8 meson :: re=8 mesonpl :: AL re=8 mesonpu :: AU teger meso:n:: el eu teger mesonpl:: teml teger mesonpu:: temu re=8 meso:: octbe :: W teger prmeter :: R= teger prmeter :: Z= teger prmeter :: Q= teger prmeter :: P= 3 teger prmeter :: DD= 4 変数名 : ICELTOT N BFORCE B PHI EPSICCG EPS W= WR {r} W= WZ {z} W= WQ {q} W3= WP {p} W4= WDD {}

170 OMP- 7 sove_iccg/7:method=!c!c !C INIT!C !C=== octe WN4 o = N =. W=.D W=.D W3=.D W4=.D eo!c=== o = N VAL= D o = el-+ el VAL= VAL - AL** * WtemLDD eo WDD=./VAL eo 不完全修正コレスキー分解 WDD= WDD: D: temlj: ALj: j j j j j j j 対角成分のみがもとの行列と変わる

171 OMP- 7 sove_iccg3/7:method=!c!c !C {r}= {b} - [A]{}!C !C=== o = N VAL= D* o = el-+ el VAL= VAL + AL*temL eo o = eu-+ eu VAL= VAL + AU*temU eo WR= B - VAL eo!c=== BNRM=.D o = N BNRM = BNRM + B ** eo BNRM= b あとで収束判定に使用 Compute r = b-[a] for = sove [M]z - = r - - = r - z - f = p = z ese - = - / - p = z ef q = [A]p = - /p q = - + p r = r - - q chec covergece r e p -

172 OMP- 7 sove_iccg4/7:method=!c!c************************************************* ITERATION ITR= N o L= ITR!C!C !C {z}= [Mv]{r}!C !C=== o = N WZ= WR eo!c=== o = N WVAL= WZ o = el-+ el WVAL= WVAL - AL * WtemLZ eo WZ= WVAL * WDD eo o = N - SW =. o = eu-+ eu SW= SW + AU * WtemUZ eo WZ= WZ - WDD*SW eo Compute r = b-[a] for = sove [M]z - = r - - = r - z - f = p = z ese - = - / - p = z p ef q = [A]p = - /p q = - + p r = r - - q chec covergece r e

173 OMP- 73 sove_iccg4/7:method=!c!c************************************************* ITERATION ITR= N o L= ITR!C!C !C {z}= [Mv]{r}!C !C=== o = N WZ= WR eo!c=== T M z LDL z r Lz o = N WVAL= WZ o = el-+ el WVAL= WVAL - AL * WtemLZ eo WZ= WVAL * WDD eo o = N - SW =. DL T o = eu-+ eu SW= SW + AU * WtemUZ eo WZ= WZ - WDD*SW eo r z z 前進代入 Forwr Substtuto 後退代入 Bcwr Substtuto

174 OMP- 74 sove_iccg5/7:method=!c!c !C RHO= {r}{z}!c !C=== RHO=. o = N RHO= RHO + WR*WZ eo!c=== Compute r = b-[a] for = sove [M]z - = r - - = r - z - f = p = z ese e - = - / - p = z p ef q = [A]p = - /p q = - + p r = r - - q chec covergece r

175 OMP- 75 sove_iccg6/7:method=!c!c !C {p} = {z} f ITER=!C BETA= RHO / RHO otherwse!c !C=== f L.eq. the!c=== o = N WP= WZ eo ese BETA= RHO / RHO o = N WP= WZ + BETA*WP eo ef!c!c !C {q}= [A]{p}!C !C=== o = N VAL= D*WP o = el-+ el VAL= VAL + AL*WtemLP eo o = eu-+ eu VAL= VAL + AU*WtemUP eo WQ= VAL eo!c=== Compute r = b-[a] for = sove [M]z - = r - - = r - z - f = p = z ese - = - / - p = z p ef q = [A]p = - /p q = - + p r = r - - q chec covergece r e

176 OMP- 76 sove_iccg7/7: METHOD=!C!C !C ALPHA= RHO / {p}{q}!c !C=== C=. o = N C= C + WP*WQ eo ALPHA= RHO / C!C===!C!C !C {}= {} + ALPHA*{p}!C {r}= {r} - ALPHA*{q}!C !C=== o = N = + ALPHA * WP WR= WR - ALPHA * WQ eo DNRM=. o = N DNRM= DNRM + WR** eo!c=== ERR = sqrtdnrm/bnrm f ERR.t. EPS the IER = goto 9 ese RHO = RHO ef eo IER = Compute r = b-[a] for = sove [M]z - = r - - = r - z - f = p = z ese - = - / - p = z p ef q = [A]p = - /p q = - + p r = r - - q chec covergece r e

177 OMP- 77 sove_iccg7/7: METHOD=!C!C !C ALPHA= RHO / {p}{q}!c !C=== C=. o = N C= C + WP*WQ eo ALPHA= RHO / C!C===!C!C !C {}= {} + ALPHA*{p}!C {r}= {r} - ALPHA*{q}!C !C=== o = N = + ALPHA * WP WR= WR - ALPHA * WQ eo DNRM=. o = N DNRM= DNRM + WR** eo!c=== ERR = sqrtdnrm/bnrm f ERR.t. EPS the IER = goto 9 ese RHO = RHO ef eo IER = r= b-[a] DNRM= r BNRM= b ERR= r / b Compute r = b-[a] for = sove [M]z - = r - - = r - z - f = p = z ese - = - / - p = z p ef q = [A]p = - /p q = - + p r = r - - q chec covergece r e

178 OMP- 78 sove_iccg/3: METHOD=!C!C***!C*** moue sover_iccg!c***! moue sover_iccg cots!c!c*** sove_iccg!c subroute sove_iccg & & N NPL NPU el teml eu temu D B & & AL AU EPS ITR IER mpct REAL*8 A-HO-Z re=8 meson :: D re=8 meson :: B re=8 meson :: re=8 mesonpl :: AL re=8 mesonpu :: AU teger meso:n :: el eu teger mesonpl:: teml teger mesonpu:: temu re=8 meso:: octbe :: W teger prmeter :: R= teger prmeter :: Z= teger prmeter :: Q= teger prmeter :: P= 3 teger prmeter :: DD= 4 Compute r = b-[a] for = sove [M]z - = r - - = r - z - f = p = z ese - = - / - p = z ef q = [A]p = - /p q = - + p r = r - - q chec covergece r e p - re=8 meso: octbe :: ALu AUu re=8 meso: octbe :: Du

179 OMP- 79 sove_iccg/3: METHOD=!C!C !C INIT!C !C=== octe WN4 o = N =. W=.D W=.D W3=.D W4=.D eo!c=== c FORM_ILU Du ALuAUu には ILU 分解された対角 下三角 上三角成分が入る 行列 [M]

180 OMP- 8 FORM_ILU/ 不完全修正コレスキー分解 : 正確には不完全修正 LU 分解 sover_iccg.f に附属 cots!c!c***!c*** FORM_ILU!C***!C!C form ILU mtr!c subroute FORM_ILU mpct REAL*8 A-HO-Z teger meso: octbe :: IW IW teger meso: octbe :: IWsL IWsU re =8:: RHS_Aj DINV A Aj!C!C !C INIT.!C !C=== octe ALuNPL AUuNPU octe DuN!C=== o = N Du= D o = INL ALu= AL eo o = INU AUu= AU eo eo j j j DuALuAUu 初期値として DALAU の値を代入する j j Du ALuAUu には ILU 分解された対角 下三角 上三角成分が入る 行列 [M]

181 OMP- 8 FORM_ILU /!C!C !C ILU fctorzto!c !C=== octe IWN IWN IW= IW= o = N o = el-+ el IWtemL= eo o = eu-+ eu IWtemU= eo o co= el-+ el = temlco D= Du DINV=./D A= ALuco o co= eu-+ eu j= temuco f j.eq. the Aj= AUuco RHS_Aj= A * DINV * Aj Du= Du - RHS_Aj ef f j.t... IWj.e. the Aj= AUuco RHS_Aj= A * DINV * Aj j = IWj ALuj= ALuj - RHS_Aj ef f j.gt... IWj.e. the Aj= AUuco RHS_Aj= A * DINV * Aj eo eo j = IWj AUuj= AUuj - RHS_Aj ef o = el-+ el IWtemL= eo o = eu-+ eu IWtemU= eo eo o = N Du=. / Du eo eocte IW IW!C=== e subroute FORM_ILU o = N o = - f A s o-zero the o j= + N f Aj s o-zero the Aj= Aj & -A*A - *Aj ef eo ef eo eo

182 OMP- 8 FORM_ILU /!C!C !C ILU fctorzto!c !C=== octe IWN IWN IW= IW= o = N o = el-+ el IWtemL= eo o = eu-+ eu IWtemU= eo o co= el-+ el = temlco D= Du DINV=./D A= ALuco o co= eu-+ eu j= temuco f j.eq. the Aj= AUuco RHS_Aj= A * DINV * Aj Du= Du - RHS_Aj ef f j.t... IWj.e. the Aj= AUuco RHS_Aj= A * DINV * Aj j = IWj ALuj= ALuj - RHS_Aj ef f j.gt... IWj.e. the Aj= AUuco RHS_Aj= A * DINV * Aj eo eo j = IWj AUuj= AUuj - RHS_Aj ef o = el-+ el IWtemL= eo o = eu-+ eu IWtemU= eo eo o = N Du=. / Du eo eocte IW IW!C=== e subroute FORM_ILU o = N o = - f A s o-zero the o j= + N f Aj s o-zero the Aj= Aj & -A*A - *Aj ef eo ef eo eo

183 OMP- 83 FORM_ILU /!C!C !C ILU fctorzto!c !C=== octe IWN IWN IW= IW= o = N o = el-+ el IWtemL= eo o = eu-+ eu IWtemU= eo o co= el-+ el = temlco D= Du DINV=./D A= ALuco o co= eu-+ eu j= temuco f j.eq. the Aj= AUuco RHS_Aj= A * DINV * Aj Du= Du - RHS_Aj ef f j.t... IWj.e. the Aj= AUuco RHS_Aj= A * DINV * Aj j = IWj ALuj= ALuj - RHS_Aj ef f j.gt... IWj.e. the Aj= AUuco RHS_Aj= A * DINV * Aj eo eo j = IWj AUuj= AUuj - RHS_Aj ef o = el-+ el IWtemL= eo o = eu-+ eu IWtemU= eo eo o = N Du=. / Du eo eocte IW IW!C=== e subroute FORM_ILU o = N o = - f A s o-zero the o j= + N f Aj s o-zero the Aj= Aj & -A*A - *Aj ef eo ef eo eo

184 OMP- 84 FORM_ILU /!C!C !C ILU fctorzto!c !C=== octe IWN IWN IW= IW= o = j N o = el-+ j el IWtemL= j j j eo o = eu-+ eu IWtemU= eo o co= el-+ el = temlco D= Du DINV=./D A= ALuco o co= eu-+ eu j= temuco f j.eq. the Aj= AUuco RHS_Aj= A * DINV * Aj Du= Du - RHS_Aj ef f j.t... IWj.e. the Aj= AUuco RHS_Aj= A * DINV * Aj j = IWj ALuj= ALuj - RHS_Aj ef j= f j.gt... IWj.e. the Aj= AUuco RHS_Aj= A * DINV * Aj eo eo j = IWj AUuj= AUuj - RHS_Aj ef o = el-+ el IWtemL= eo o = eu-+ eu IWtemU= eo eo o = N Du=. / Du eo eocte IW IW!C=== e subroute FORM_ILU o = N o = - f A s o-zero the o j= + N f Aj s o-zero the Aj= Aj & -A*A - *Aj ef eo ef eo eo

185 OMP- 85 FORM_ILU /!C!C !C ILU fctorzto!c !C=== j octe IWN j IWN IW= j j j IW= o = N o = el-+ el IWtemL= eo o = eu-+ eu IWtemU= eo o co= el-+ el = temlco D= Du DINV=./D A= ALuco o co= eu-+ eu j= temuco f j.eq. the Aj= AUuco RHS_Aj= A * DINV * Aj Du= Du - RHS_Aj ef f j.t... IWj.e. the Aj= AUuco RHS_Aj= A * DINV * Aj j = IWj ALuj= ALuj - RHS_Aj ef j< f j.gt... IWj.e. the Aj= AUuco RHS_Aj= A * DINV * Aj eo eo j = IWj AUuj= AUuj - RHS_Aj ef o = el-+ el IWtemL= eo o = eu-+ eu IWtemU= eo eo o = N Du=. / Du eo eocte IW IW!C=== e subroute FORM_ILU o = N o = - f A s o-zero the o j= + N f Aj s o-zero the Aj= Aj & -A*A - *Aj ef eo ef eo eo

186 OMP- 86 FORM_ILU /!C!C !C ILU fctorzto!c !C=== j octe IWN j IWN IW= j j j IW= o = N o = el-+ el IWtemL= eo o = eu-+ eu IWtemU= eo o co= el-+ el = temlco D= Du DINV=./D A= ALuco o co= eu-+ eu j= temuco f j.eq. the Aj= AUuco RHS_Aj= A * DINV * Aj Du= Du - RHS_Aj ef f j.t... IWj.e. the Aj= AUuco RHS_Aj= A * DINV * Aj j = IWj ALuj= ALuj - RHS_Aj ef f j.gt... IWj.e. the Aj= AUuco RHS_Aj= A * DINV * Aj eo eo j = IWj AUuj= AUuj - RHS_Aj ef o = el-+ el IWtemL= eo o = eu-+ eu IWtemU= eo eo o = N Du=. / Du eo eocte IW IW!C=== e subroute FORM_ILU j> o = N o = - f A s o-zero the o j= + N f Aj s o-zero the Aj= Aj & -A*A - *Aj ef eo ef eo eo

187 OMP- 87 FORM_ILU /!C!C !C ILU fctorzto!c !C=== j octe IWN j IWN IW= j j j IW= o = N o = el-+ el IWtemL= eo o = eu-+ eu IWtemU= eo o co= el-+ el = temlco D= Du DINV=./D A= ALuco o co= eu-+ eu j= temuco f j.eq. the Aj= AUuco RHS_Aj= A * DINV * Aj Du= Du - RHS_Aj ef f j.t... IWj.e. the Aj= AUuco RHS_Aj= A * DINV * Aj j = IWj ALuj= ALuj - RHS_Aj ef j< f j.gt... IWj.e. the Aj= AUuco RHS_Aj= A * DINV * Aj eo eo j = IWj AUuj= AUuj - RHS_Aj ef o = el-+ el IWtemL= eo o = eu-+ eu IWtemU= eo eo o = N Du=. / Du eo eocte IW IW!C=== e subroute FORM_ILU 実はこの IF 文の中は通らない 理由は後述 したがって ALu= AL AUu= AU j>

188 OMP- 88 j= j< j> /3 j ある要素 に接続する下三角成分 の上三角成分 j が : j= 自身である場合Du 更新 j j j j j

189 OMP- 89 j= j< j> /3 j ある要素 に接続する下三角成分 の上三角成分 j が : j< の下三角成分である場合 ALu-j 更新 j j j j j

190 OMP- 9 j= j< j> 3/3 j ある要素 に接続する下三角成分 の上三角成分 j が : 3 j> の上三角成分である場合 AUu-j 更新 実際は 3 に該当する場合は無し j j j j j

191 OMP- 9 sove_iccg3/3: METHOD=!C!C***************************************************** ITERATION ITR= N o L= ITR!C!C !C {z}= [Mv]{r}!C !C=== o = N WZ= WR eo!c=== o = N WVAL= WZ o = el-+ el WVAL= WVAL - ALu * WtemLZ eo WZ= WVAL * Du eo o = N - SW =. o = eu-+ eu SW= SW + AUu * WtemUZ eo WZ= WZ - Du*SW eo これ以下の処理は sove_iccg と全く同じ Compute r = b-[a] for = sove [M]z - = r - - = r - z - f = p = z ese - = - / - e p = z ef q = [A]p = - /p q = - + p r = r - - q chec covergece r p -