Hitachi SR8000/MPP フルノードによる並列有限要素法コード「GeoFEM」
における MPI/OpenMP ハイブリッド並列ソルバーの性能評価
中島 研吾（財団法人 高度情報科学技術研究機構，nakajima@tokyo.rist.or.jp）
奥田 洋司（東京大学大学院工学系研究科 システム量子工学専攻，okuda@q.t.u-tokyo.ac.jp）

1．はじめに
近年，超並列計算機の世界では Hitachi SR8000，「地球シミュレータ」 [1] などに代表されるような，複数のプロセッサが共 有メ モリ ユニ ット を構 成す る SMP（Symmetric Multiprocessor）クラスタ型のアーキテクチュアが一般的である。アメリカの ASCI（Accelerated Strategic Computing Initiative） 2] のハードウェア群もすべてこの SMP クラスタ型を採用している。このようなアーキテクチュアにおいては Loop Directive と Message Passing を組み合わせたいわゆる Hybrid 型の並列プログラミングが適していると考えられる。す なわち，SMP ノード（共有メモリユニット）間の通信には MPI [3] などのMessage Passing，SMP ノード内の並列化には OpenMP [4] などの Loop Directive を使用するというものである。このようなプログラミングモデルに関する研究 [5][6] はここ数年数多く実施されているが，多くは NAS Parallel Benchmarks（NPB） [7] に代表される構造格子を対象としたもので非構造格子を使用した例はほとんどない。
本稿では有限要素法，すなわち非構造格子における前処理付き並列反復法ソルバーに対して Hybrid型並列 プロ グ ラミ ン グモ デ ルを 適 用し ， 東京 大 学情報基盤センター Hitachi SR8000/MPP フルノード（128SMP ノード，1024 プロセッサ）を使用して最大 805,306,368 自由度までの計算を実施した例を示す。使用した有限要素法コードおよび線形ソルバーは財団法人高度情報科学技術研究機構によって開発されている「GeoFEM」 [8] の一部である。GeoFEM は文部科学省科学技術振興調整費「高精度の地球変動予測のための並列ソフトウエア開発に関する研究」の固体地球分野において，固体地球変動予測シミュレーションのための並列有限要素法プラットフォームとして開発されている。
以 下 ， （ 1） 反 復 法 を 使 用 し た 有 限 要 素 法 の 並 列 化 ， （ 2） 並 列 計 算 用 局 所 前 処 理 手 法 に つ い て ，「 GeoFEM」 に お け る 実 装 例 を 述 べ た 後 ， （3） SMPク ラ ス タ 型 並 列 計 算 機 向 け に 最 適 化 さ れ た ソ ル
バーの概要と計算結果について述べる。

2．反復法を使用した有限要素法の並列化 [8][9]並列計算で扱うデータのサイズ（メッシュ数）は非常に大きいため，全体領域を一括して取り扱うことは困難で，全体データを部分領域（局所データ）に分割する必要がある。有限要素法は差分法などと比較して並列化が困難であると考えられてきた。間接データ参照があるため，1 プロセッサ（Processing Element，PE） あ た り の計 算 効 率 は 差分法と比較して低いが，有限要素法の処理は基 本的に要素単位の局所的な処理であり，領域間の通信は図 1 に示すように線形ソルバーの部分だけで生じる。この特性を最大限利用し，適切なデータ構造を設定，並列計算に適した反復法を採用することによって後述するように 95%を越えるような並列化効率を達成することも可能である。
GeoFEM では領域間の通信の記述には MPI を使用している。差分法などに使用されている構造格子（Structured Grids）に関しては MPI 固有の領域間通信用のサブルーチンが準備されているが，有限要素法に代表される非構造格子（Unstructured Grids）では，プログラム開発者が独自にデータ構造と領域間通信を設計しなくてはならない。

PE#0

Local Data

FEM codeLinear Solvers
MPI

PE#1

Local Data

FEM codeLinear Solvers

PE#(M-2)

Local Data

FEM codeLinear Solvers
MPIPE#(M-1)

Local Data

FEM codeLinear Solvers

Total PE number is (M).
図１．並列有限要素法の処理
有限要素法の処理は要素単位の局所的な処理であり並列化が容易。通信は線形ソルバーでのみ発生。

GeoFEM では，領域間の負荷バランスを考慮して全体領域を「節点ベース（Node-Based）」に分割して い る 。 す な わ ち ，各 部 分領 域 ，各 プ ロセ ッ サ （ PE） で 扱 う 節 点 数 が 均 等 に な る よ う に 領 域 を 分割している。（図 2（a）参照）。
節点ベースの領域分割では各局所データは以下の情報を含んでいる：（1）

本来その領域に割り当てられた節点（2）（1）の節点を含む要素（3）（2）の要素に含まれる節点のうち（1）に含まれないもの（4）

領域間の通信テーブル（Communication Table）（5）

その他，節点／要素／面グループ等

上記のうち（2），（3）の情報は各領域において，要素単位のマトリクス生成を実施するために必要な情報である。そのため図 2（b）に示すように各領域境界では各領域によって共有されるオーバーラップ要素が生じる。
部分領域内の節点は領域間通信の見地から以下の 3 種類に分類される：・内点（Internal Nodes，上記の（1））

・

外点（External Nodes，上記の（3））

・境界点（Boundary Nodes，「内点」のうちで他領域の「外点」となっている点）

例えば図 2（b）の「PE#2」の領域に注目すると，内点：｛1，2，3，4，5，6｝，外点：｛7，8，9，10，11，12｝，境界点：｛1，4，5，6｝となる。
隣接する領域間の通信に関して記述しているデータが「通信テーブル（Communication Tables）」

である。「境界点」のデータが各隣接領域に「送信（send）」され，送信先の各領域では「外点」データとして「受信（receive）」される（図 3 参照）。

PE#121221716PE#023241813112514158710978911101252PE#311134510123891112109111245884714109666543151233162201265119

PE#0

PE#141

PE#2（a） 全 体 メ ッ シ ュ と 内 点

2PE#377123PE#2（b） 各 領 域 に 属 す る 内 点 ， 要 素 ， 外 点
領域間でオーバーラップする要素が生じる

図2．GeoFEMの局所データ：節点単位の領域分割（a）送信モード

156（b）受信モード

7PE#0PE#0111051296384PE#31214134510123891112109111287471113105963PE#21284PE#312767134512389111210911124510652141587126PE#23

図3．GeoFEMの局所データ：通信テーブル

3．並列計算用局所前処理手法 [8][9]
図 1 にも記述されているように，GeoFEM では領域間の通信は線形ソルバーの部分でのみ生じる。
前処理つき反復法における計算プロセスは以下の 4 種類に分類される：（1） 列ベクトル積（2） ベクトル～ベクトル内積（3） ベクトル（およびその実数倍）の加減（4） 前処理
このうち（3）を除く 各プ ロセスで は領 域間の通 信が発生する。（1）は計算前に 2．で述べた通信を[A]実施すれば局所的な処理が可能である。（2）{x}は MPI_ALLREDUCE などのサブルーチンを使用して容易に達成可能である。

1 A 2 3456PE#0（4）については前処理手法によって異なる。
例えば代表的な前処理手法である不完全 LU

PE#1分解（Incomplete LU Factorization, ILU），不完 全 コ レ ス キ ー 分 解 （ Incomplete Cholesky

PE#2Factorization, IC）などの手法は前進／後退代入により大域的な変数の依存性が生じるた

PE#3

め，並列化が困難である。単独プロセッサを使用した計算の場合，Fill-in のない ILU(0)法を前処理として使用すると，前処理計算部分が全体の 50%程度を占めるため[9]，前進／後退代入部分の並列化は計算効率の向上のために不可欠である。

図 4．局所 ILU（0）法における前進後退代入
た と え ば PE#0 の 要 素「A」に 注 目 す る と 6 個 の 非 対 角
成 分 が あ る が ，こ の う ち 4，5，6，は 他 領 域 に 属 す る「 外 点 」 で あ る た め ， 局 所 ILU（ 0） 法 に お け る 前 進
後退代入においては 0 とみなされる。

GeoFEM では局所前処理法（局所 ILU(0)z法，Localized ILU(0)） [4] を 使 用 し て い る 。

Uniform Distributed Force inz-direction @ z=Zmax

局所 ILU(0)法は一種の「擬似」ILU(0)法である。局所 ILU(0)法では前進／後退代入計

Uy=0 @ y=Ymin

算時に領域外からの影響（すなわち外点の影

Ux=0 @ x=Xmin響）を 0 とすることによって，前処理の局所化を行い，並列性の高いアルゴリズムを実現している（図 4 参照）。

Nz-1
Ny-1

図 5 に示す三次元固体力学の例題（弾性静解析，一様単純引張）を Hitachi SR2201（東京大学情報基盤センター，1024 PE，ピーク性能 300GFLOPS）を使用した計算結果を図6 に示す。各 PE における自由度数が一定とUz=0 @ z=ZminyNx-1x

図 5．三次元固体力学例題（弾性静解析）なるような問題設定となっている。この問題
では一節点あたり 3 成分の自由度（x,y,z 方向の変 位） が 存在し てい る ため， これ ら をブロ ック化して扱うブロック IC(0)前処理 [10] を使用している。
最大 1.97×10 8 自由度の問題を 1024 PE を使用して，68.7 GFLOPS の処理性能を達成している。これはピーク性能の約 22.4 %にあたる。図 6（b）は並列化効率（（通信およびそれに要する並替計算等を除いた演算時間／全演算時間））に関するデータである。PE 数が増加しても並列化効率 95%以上を維持しており，2．で述べたデータ構造と局所 IC(0)法によって高い並列化効率が得られている。

70.0

100.0

60.0

95.0
90.0

40.0%

GFLOPS

50.0

85.0

30.0
80.0

20.0

75.0

10.0
0.0

70.00128

25638451264076889610240128

256384PE #5126407688961024

PE #（a）PE数 とGFLOPS値 の 関 係：ほ ぼ 正 比 例（b）PE数 と 並 列 化 効 率 の 関 係：充 分 大 き い 問 題 規 模の 場 合 ， 並 列 化 効 率 は95%以 上

図 5．Hitachi SR2201 を使用した GeoFEM ソルバー計算例
三 次 元 弾 性 静 解 析 ： Block 局 所 IC(0)CG 法
最 大 1024 PE， 1.97×10 8 自 由 度 ， 68.7 GFLOPS， ピ ー ク 性 能 の 約 22.4%
各PEの 自 由 度 数 を 固 定 ● ： 1.23×10 4 ， ■ ：9.83×10 4 ， △ ：1.92×10 5

4．SMPクラスタ型アーキテクチュア向けソルバーの開発「GeoFEM」の SMP クラスタ用並列ソルバーは，もともとベクトル計算機である「地球シミュレータ」を主要なターゲットとして開発 さ れ て い る 。「 地 球 シ ミ ュ レ ータ」は，8 つのベクトルプロセッサから構成される SMP ノードが640 ノード，合計 5,120 プロセッ

Node-0Node-1

サであり，Hitachi SR8000 とは非常によく似た構成である。
「GeoFEM」の SMP クラスタ型

Node-2Node-3用並列ソルバーでは以下の 3 レベルの並列性が考慮されている：•SMP ノード間： MPI•SMP ノード内： OpenMP•

各プロセッサ： ベクトル化

Partitioning of analysis
domainNode-0Node-1

PE
PE
PE
PE
PE
PE
PE
PE

M
emoryPE
PE
PE
PE
PE
PE
PE
PE

M
emoryPE
PE
PE
PE
PE
PE
PE
PE

M
emoryPE
PE
PE
PE
PE
PE
PE
PE

M
emory

Node-2Node-3

図 7．SMP クラスタ型アーキテクチュアにおける領域分割
1 領 域 が 1SMP ノ ー ド に 対 応 し て い る

Hitachi SR8000 上で使用する場合は各プロセッサに関しては「疑似ベクトル」オプション [11] を適用している。並列計算にあたっては分割された各領域は各 SMP ノードに対応している（図 7 参照）。非構造格子を使用した計算において高いベクトル／並列性能を得るためには以下の 3 点が重要である：•

局所的な処理•

連続メモリアクセス•

充分なループ長

これらの要件を満たすように，鷲尾らによって考案されたベクトルプロセッサおよび SMP ノード向けのオーダリング手法 [12][13]（図 8，9 参照）と GeoFEM で使用されている局所前処理手法を統合した手法を適用し，特に ILU，IC 前処理における性能を向上させている。今回使用しているオーダリング（PDJDS/CM-RCM，Parallel Descending Order Jagged Diagonal Storage/Cyclic Multicolor-RCM）の詳
細については文献[14]を参照されたいが，以下の 4 段階から成っている：①

各領域に RCM（Reverse Cuthil-McKee）オーダリング [10] を適用し，依存性のない Hyperplane を生成する（図 8（a））。

②ループ長が均等となるように CM（Cyclic Multicolor）オーダリングを適用する（図 8（c））。

③DJDS（Descending Order Jagged Diagonal Storage）オーダリングによって各カラー内で非対角成分の多い順に並び替え，係数行列に関してメモリの連続した一次元配列を生成する（図 9）。

④SMP ノード内での負荷が均等となるようサイクリックなオーダリングを実施する。
8 9 10 11 12 13 14 15
7 8 9 10 11 12 13 14

3 4 3 4 3 4 3 41 2 1 2 1 2 1 2

4 1 2 3 4 1 2 3
3 4 1 2 3 4 1 2

6 7 8 9 10 11 12 13

3 4 3 4 3 4 3 42 3 4 1 2 3 4 15 6 7 8 9 10 11 12

1 2 1 2 1 2 1 2

1 2 3 4 1 2 3 44 5 6 7 8 9 10 11
3 4 5 6 7 8 9 10

3 4 3 4 3 4 3 41 2 1 2 1 2 1 24 1 2 3 4 1 2 3
3 4 1 2 3 4 1 2

2 3 4 5 6 7 8 9

3 4 3 4 3 4 3 4

2 3 4 1 2 3 4 1

1 2 3 4 5 6 7 8

1 2 1 2 1 2 1 2

1 2 3 4 1 2 3 4(a) Hyperplane/RCM(b) Multicoloring: 4 colors(c) CM-RCM: 4 colors

図 8．CM-RCM(Cyclic Multicolor+RCM)オーダリング(a) 非 対 角 成 分 の 多 い 順 番 に 並 び 替 え(b) 一 次 元 圧 縮 配 列

図 9． DJDS（Descending Order Jagged Diagonal Storage）オーダリング
上記①～④のうち，③はベクトル化向け，④は SMP ノード内並列化向けのオーダリングであり，①，②は両者に効果的なオーダリングである。
図 5 に示した三次元弾性静解析について，最大 8.05×10 9 自由度の問題を 128 ノード（1024 プロセッサ ） を使 用し て ，335 GFLOPS の 処 理性 能 を達成している（1-SMP ノードあたり 128 3 節点，6,291,456 自由度）。これはピーク性能の約 18.6 %にあたる。また 95%以上の高い並列性能を発揮している（図 10，11 参照）。
1000.0300
250

GFLOPSGFLOPS/Work Ratio(%)

350

200
150

10.0

95.9%

10081.2%50
1.E+06100.0

1.E+0797.5%

1.E+081.E+091.E+101.01DOF

101001000

SMP-Node #図 10. 三次元弾性静解析結果（SR8000/MPP128 ノード使用）

図 11. 三次元弾性静解析結果（SR8000/MPP）
SMP ノード数と GFLOPS 値の関係

● ：GFLOPS 値 ， ○ ： 並 列 性 能

1-SMP ノ ー ド あ た り 128 3 節 点 ，6,291,456 自 由 度

最 大 805,306,368 自 由 度 ， 335.2G FLOPS（ ピ ー ク 性 能 1.8 TFLOPS の 約 18.6%）

最 大 805,306,368 自 由 度 ， 335.2G FLOPS（ ピ ー ク 性 能 1.8 TFLOPS の 約 18.6%）

図 12 は 1-SMP ノードを使用して問題規模

100を変化 さ せ た 場 合 （ 最 小 ： 16 3 節 点 ，12,288の結果で Hitachi FORTRAN の XCLOCK 関数[11]を使用して，SMP ノード内での通信，同期オーバーヘッドを測定し，並列化効率を算出したものである。この図によると最小問題
では通信 ， 同 期 オー バ ー ヘ ッ ドは 20％ 以 上であるが，問題規模が大きくなるにつれて減少し，40 3 節点，192,000 自由度（1 プロセッ
サあたり 24,000 自由度）で 10%以下，64 3 節点 ， 786,432 自 由 度 （ 1 プ ロ セ ッ サ あ た り98,304 自由度）で 5%以下となっている。

95

Work Ratio (%)

自由度，最大：128 3 節点，6,291,456 自由度）

9085
80
75
1.E+04

1.E+05

1.E+06

1.E+07

DOF #図 12. SR8000/MPP 1-SMP ノード内
並列性能（XCLOCK サブルーチン使用）
40 3 節 点 ，192,000DOF（24,000 自 由 度/PE） で 90%以 上

続いて，オーダリングの影響を検討するため，図 13 に示すように，①PDJDS/CM-RCM（本稿で提案している手法），②PDCRS/CM-RCM（オーダリングは①と同様だが，最内ループ長が短い），③オーダリング無しの 3 つの場合について，1-SMPノードを使用し，問題規模を変更した計算を実施した。なお，この比較計算には東京大学情報基盤センターの Hitachi SR8000/128（ノードあたりピーク性能 8GFLOPS，SR8000/MPP は 14.4GFLOPS）を使用した。計算結果を図 14 に示す。

PDJDS/CM-RCM

2.00

PDCRS/CM-RCM CRS no re-ordering
short innermost loop

GFLOPS

1.50

1.00

0.50

0.00
1.E+04

1.E+05

1.E+06

1.E+07

DOF #図 13．各種のオーダリング，係数格納方法

図 14. SR8000/128 1-SMP ノード内並列性能（オーダリング，係数格納方法の効果）●：PDJDS/CM-RCM，○：PDCRS/CM-RCM，▲：No-Ordering

図 14 に示すように，問題規模が小さい場合は最内ループ長の短い PDCRS/CM-RCM が有利であるが問題規模が大きくなるにしたがって PDJDS/CM-RCM の方が性能が良くなり，もっとも規模が大きい ケ ー ス （ 128 3 節 点 ， 6,291,456 自 由 度 ） で は 1： 2 程 度 の 性 能 比 に な っ て い る 。 図 14 に お け る
PDJDS/CM-RCM の挙動は，ベクトルプロセッサの場合とよく似た傾向である。オーダリングをしない場合の性能は非常に低く，他のケースの 10 分の 1 程度である。SMP ノード内並列化，ベクトル化の効果が共に得られていないためであると考えられる。

5．結論および将来への展望
並列有限要素法プラットフォーム GeoFEM の前処理付き並列反復法ソルバーに対して Hybrid 型並列 プ ロ グ ラ ミ ン グ モ デ ル を 適 用 し ， 東 京 大 学 情 報 基 盤 セ ン タ ー Hitachi SR8000/MPP フ ル ノ ー ド（128SMP ノード，1024 プロセッサ）を使用して最大 805,306,368 自由度までの三次元弾性静解析を実施し，335GFLOPS（ピーク性能 1.8TFLOPS の 18.6%）の性能を得た。4 段階のオーダリングを適用する PDJDS/CM-RCM 法の効果は大きく，オーダリングを全く実施しない場合と比較して，10 倍以上の効 率 が 得 ら れ て い る 。 本 稿 の 手 法 は 本 来 ベ ク ト ル計 算機 向け に開 発さ れた 手法 であ るが ，SR8000の疑似ベクトルオプションを使用して，高いピーク性能比が得られた。
今回は単純な形状に関する線形問題の例を紹介した。「GeoFEM」の本来対象としている固体地球シミュレーションの代表的なものとして，地震発生サイクル予測のための断層接触シミュレーションがある。このようなシミュレーションは断層の複雑形状を扱う非線形問題であり，接触面の拘束条件に起因するペナルティ項の ために条件数 が増加する [15] 。このような問題に対しては，通常のILU(0)，IC(0)など の 前 処 理手 法 で は 収 束解 を 得 る こ とが できず，特殊な前処理手法が開発されている [15] 。今後はこれらの前処理手法に関してもSMPクラスタ型アーキテクチュア向けの改良を進め，より大規模な問題にチャレンジしていく予定である。

謝 辞
本研究は文部科学省科学技術振興調整費「高精度の地球変動予測のための並列ソフトウェア開発に関する研究」の一環として実施中の「GeoFEM」プロジェクトの成果の一部である。前処理手法に関して鷲尾 隆 氏（NEC） に 貴重 な助 言を 頂い た。 この場を借りて深甚なる謝意を表するものである。
今回の研究紹介の機会を与えてくださった金田康正教授（東京大学情報基盤センター）に篤く御礼申し上げます。

参考文献
[1]
[2]
[3]
[4]
[5]

http://www.es.jamstec.go.jp/
http://www.llnl.gov/asci/
http://www.mpi.org/
http//www.openmp.org/
Falgout, R. and Jones, J. : "Multigrid on Massively Parallel Architectures", Sixth European Multigrid
Conference, Ghent, Belgium, September 27-30, 1999.
[6] Cappelo, F. and Etiemble, D. :"MPI versus MPI+OpenMP on the IBM SP for the NAS Benchmarks",SC2000 Technical Paper, Dallas, Texas, 2000.
[7] http://www.nas.nasa.gov/Research/Software/swdescription.html#NPB/[8] http://geofem.tokyo.rist.or.jp/[9] K.Nakajima and H.Okuda, IJCFD, Vol.12, pp.315-322, 1999.
[10] Dongarra, J.J., Duff, I.S., Sorensen, D.C. and van der Vorst, H.A. : "Numerical Linear Algebra for
High-Performance Computers", SIAM, 1998.
[11] http://www.hitachi.co.jp/Prod/comp/hpc/foruser/sr8000/[12] Washio, T., Maruyama, K., Osoda, T., Shimizu, F. and Doi, S. : "Blocking and reordering to achieve
highly parallel robust ILU preconditioners", RIKEN Symposium on Linear Algebra and its Applications,
The Institute of Physical and Chemical Research, 1999, pp.42-49.
[13] Washio, T., Maruyama, K., Osoda, T., Shimizu, F. and Doi, S. : "Efficient implementations of blocksparse matrix operations on shared memory vector machines", SNA2000 : The Fourth International
Conference on Supercomputing in Nuclear Applications, 2000.
[14] Nakajima, K. and Okuda, H. : " Parallel Iterative Solvers for Unstructured Grids using an OpenMP/MPI
Hybrid Programming Model for the GeoFEM Platform on SMP Cluster Architectures", International
Workshop on OpenMP: Experiences and Implementations (WOMPEI 2002), Lecture Notes in Computer
Science 2327 (in press), May 2002.
[15] Nakajima, K. and Okuda, H, : "Parallel Iterative Solvers with the Selective Blocking Preconditioning
for Simulations of Fault-Zone Contact", GeoFEM 2001-010, RIST/Tokyo, 2001.