並列化

この節では得られた線形ソルバを並列化する方法について述べる．

6.4.1 領域分割

数値流体計算においては反復によって計算を解く場合，緩和法などで大規模問題の並列 処理するために計算領域を分割し，それぞれを独立に計算することで問題を処理する手法 が一般的である．このとき領域間データの依存性はある程度無視されるため，並列計算で 処理する場合と逐次計算で処理する場合の反復回数は，並列計算を行った方が多くなる．

分割された領域を適切に分割する．分割の方法は単純な短冊状の分割も可能であるが，

計算領域全体が複雑な形状をしている場合，単純に分割したときの分割領域間のサイズが まちまちになってしまい性能低下の原因となることがある．

今回の計算においてはメッシュの分割に自動分割エンジンであるミネソタ大学の開発し

たMETISを用いる[9]．分割領域間に，のりしろとなる領域としてゴーストノードを設定

し，計算ループを回す度に，ゴーストノードと領域境界点の値をコピーし，領域ごとに 新たな境界値を設定する．ゴーストノードの境界点との値のコピーはMPIの2点間通信 を用いて行われる．本研究のコードでは，MPI Isend, MPI Irecv と MPI Waitall の非ブ ロッキング通信を用い，通信時間とMPIメッセージ送出バッファを構築する計算時間の 重ね合わせが行われている．

緩和法のループ完了を判定する際の誤差の値は各領域間で最も大きい値を用いて判断 せねばならない．なぜならある領域での誤差がループを脱出する許容値より低かったとし ても，別の領域ではそうとは限らないためである．このような計算にはMPI の集合通信 である MPI Allreduceを用いて誤差の値をすべての領域上で同期させる．

MPI通信による並列化はコード上に概要で示すと次のようになる．

r₀ =b−Ax₀ p₀ =r₀

do k = 0, kmax q_k=Ap_k

!ゴーストノードにおける q_k の値の交換 α= (r_k,r_k) / (p_k, q_k)

!MPI領域間でαの値をAllreduceで足し合わせこれを新しいalphaの値とする x_k+1 =x_k+α_kp_k

r_k+1 =r_k−α_ks_k rtr= (r_k+1,r_k+1)

!MPI領域間で rtrの値をAllreduceで足し合わせこれを新しいrtrの値とする if (rtr < ε) exit

β_k = (r_k+1,r_k+1) / (r_k,r_k) p_k+1 =r_k+1+β_kp_k

end do

6.4.2 ループレベル分割

プログラム中のDOループをOpenMPを使って並列化する．

OpenMPでの並列化の手法はプログラム中に並列化ディレクティブを挿入することで

ある．CG法のアルゴリズムには初期化と内積計算の部分がそうであるが，第2章で述べ た2つのプログラミングモデルによって並列化部分が異なる．これらは以下のように表現 できる．

Hybrid マスターオンリーモデルCG法

...

do k = 0, kmax

!$OMP parallel do

!ここでスレッドを開く q_k=Ap_k

! ここでスレッドは閉じられる

!ゴーストノードにおける q_k の値の交換

!$OMP parallel do

! スレッドを開く α= (r_k,r_k) / (p_k, q_k)

! ここでスレッドは閉じられる . . .

end do

HybridMPI+OpenMP SPMDモデル

...

!$OMP parallel do

! 反復部分全体に渡ってスレッドを生成する ....

do k = 0, kmax

!$OMP do

!各々のスレッドが分割されたループを処理する q_k=Ap_k

!$OMP end do

!$OMP single

! MPI呼び出しの前に同期してから単一スレッドで実行

ゴーストノードにおける q_k の値の交換 (MPI通信)

!$OMP end single

!再び複数スレッドで実行 . . .

end do

!$OMP end parallel

!ここでスレッドを閉じる