2.6 流体解析プログラムのデータ&ループ構造に関して
宇宙航空研究開発機構 宇宙科学研究所
高木 亮治
2.6.1

はじめに一般的に流体解析プログラムは高いメモリバンド幅を要求するプログラムであり、高いメモリバンド幅を有するベクトル計算機との相性が良いと言われてきた。しかしながら、科学技術計算の分野においても専用計算機的なベクトル計算機から汎用計算機的なスカラー並列計算機へと計算機アーキテクチャの移行が行われ、しかも急激に増加する演算性能に比して、メモリ性能の伸びが追いつかず、B/F は確実に減少する傾向にある。例えば、現在検討が進められているエクサフロップス級計算機のうち、様々な計算に適用可能な汎用性を指向したものでは、B/F は 0.1 へと減少し次世代スーパーコンピュータ「京」の 1/5 と予想されている。そのため低い B/F においてもそれなりの実行性能を発揮する流体解析プログラムを開発する必要がある。
ここでは、低 B/F を前提とした圧縮性流体解析プログラムの実現を目指して、どの様なプログラム構造が適切かについての検討を試みた。まず、手始めに JAXA が開発している UPACS の主要部分を抜き出してカーネルプログラムを作成し、これを用いてループ構造の検討を行った。

2.6.2

ループ構造の検討

カーネルプログラムでは、一般曲線座標系で記述された支配方程式を対象として、右辺対流項の計算部分 (基本変数を用いた 2 次精度 MUSCL、van Albada のリミター、数値流速は SHUS) と左辺時間積分 (1 次精度 Euler 陽解法) を実装している。現時点では粘性項、陰解法による時間積分、境界条件は考慮していない。対象となるこれらの計算を行う際のループ構造として• ループ A：従来のループ• ループ B：局所性を意識した (空間スイープの 3 重ループの数をできるだけ減らした) ループを実装し、計算速度の違いを調べた。
ループ A は従来のプログラムに良く見られる構造で、圧縮性流体の離散方程式をプログラムとして実装する際に、対流項の計算に用いるセル面での物理量の外挿 (MUSCL+リミター)、セル面における対流項の数値流束の計算、(セル面における粘性流束の計算)、更新ベクトル ∆Q の計算、時間積分 (左辺の計算) のように、それぞれの計算を分割して、それぞれに対して空間の多重ループで計算を実行する。そのために 1 ステップの計算を実行するのに、何度も空間スイープを行うことになる。ループ A を模式的に書くと以下のようになる。ここで、dir=1,3のループは 3 次元のインデックス方向 (i、j 、k 方向) のループである。

do dir=1,3
do k=1,kmax
do j=1,jmax
do i=1,imax
MUSCL による外挿
enddo
enddo
enddo

1

do k=1,kmax
do j=1,jmax
do i=1,imax
数値流束の計算
enddo
enddo
enddo
do k=1,kmax
do j=1,jmax
do i=1,imax
更新ベクトルの計算
enddo
enddo
enddo
enddo
do k=1,kmax
do j=1,jmax
do i=1,imax
時間積分enddo
enddo
enddo
この様にループ A では 3(MUSCL、数値流束、∆Q) × 3 方向 (i,j,k 方向)+1(時間積分) の計 10 回の空間スイープを実行することになる。空間を何度もスイープすることはそれだけメモリアクセスが増加し、キャッシュを有効に活用することができなくなる。そのため、アルゴリズムの観点からメモリアクセスを減らすことを意図して、空間のスイープを極力減らし、データの再利用性を心掛けるループ B を考える。ループ B では空間のスイープは必要最低限な 2 回 (右辺の計算と左辺の計算) とした。つまり、あるセル (i,j,k) に着目し、そのセルで必要な右辺の計算を全て行い、それが終わると次のセルに移動する。この様にして 1 つの空間スイープで右辺を全て計算する。
次に左辺に関して、2 つ目の空間スイープで計算を済ませることとする。ここで注意すべき点として、一般には数値流束の計算はセル面のループで回すが、これをセルのループで回すと何も考慮しないと 1 つの面での数値流束を 2 回計算することになる。これを避けるためにはフラグを設定して計算したかどうかを判別する必要があるが、ここでは構造格子である利点を活かして、個々のセルでは各方向でインデックスが増える方向のセル面での数値流束の計算を実施することとする。そのためインデックスの始点側境界 (i = 1 or j = 1 or k = 1) での処理が必要となる。

do k=1,kmax
do j=1,jmax
do i=1,imax
do dir=1,3
MUSCL による外挿、数値流束の計算、更新ベクトルの計算
境界での処理 (MUSCL、数値流束、更新ベクトルの計算)
enddo

enddo
enddo
2

enddodo k=1,kmax
do j=1,jmax
do i=1,imax
時間積分
enddo

enddo
enddo
ループ A、B ともに最外ループである k のループを対象に OpenMP でスレッド並列化を行った。

数値実験

2.6.3

ループ A および B を実装したカーネルプログラムを幾つかの計算機上で実行し、それぞれの性能測定を実施した。計算対象は単純な立方体格子であり、格子ブロックサイズ (=ループ長) を変えて測定を行った。性能測定に用いた計算機環境を表 1 に示す。コンパイル時の最適化オプションは富士通コンパイラー (JSS)、インテルコンパイラー (PC-S)、GNU コンパイラー (PC-S, PC-N) でそれぞれ、-O5、-fast、-O3 を用いた。
表 1: 計測環境
CPU Clock
GFLOPS
[GHz]

Name

CPU

# of CoresJSS

SPARC64 VII

42.5

PC-S

Core i7-3960X
Sandy Bridge

63.3

Core i7-965
Nehalem

4

PC-N

2.6.3.1

3.2

Memory bandwidth[GB/s]

Compiler

4040

Fujitsu

158.4

51.2 (DDR3-1600)
orIntelor42.7 (DDR3-1333)

GNU

31.2 (DDR3-1333)

GNU

51.2

JSS での測定

JSS(富士通 FX1、CPU は SPARC64VII) 上でループ A とループ B の比較を行った。計算は 1 プロセス、4 スレッドである。図 1(a) にループ A、B それぞれの計算時間と L2 キャッシュのミス率の傾向を示す。図より、ループ B は狙い通り L2 キャッシュのミス率が半減していることがわかる。しかしながら計算速度としてはループ A の方がループ B よりも速い結果となった。この原因であるが、ループ B はキャッシュのミス率は改善されたが、ループの中身が大きくなった分、レジスター溢れや、パイプライン処理の最適化など、他の性能要因の影響によって性能が悪化した可能性がある。コンパイラーの最適化能力と関連するので、引き続き詳細な検討が必要である。
両方に共通する傾向として、経過時間に細かな振動が見られるが、これはスレッド数 4 の周期となっており、格子ブロックサイズがスレッド数で割りきれる場合が局所的に経過時間が短く、余りが 3 の場合に局所的に経過時間が長くなるためである。また、格子ブロックサイズが小さい場合はデータがキャッシュに収まるため、演算ネックとなり格子ブロックサイズの増加とともに経過時間が増加している。一方、格子ブロックサイズが 60 を越える辺りからデータがキャッシュから溢れるため、メモリバンド幅ネックとなり、格子ブロックサイズが増加しても経過時間は殆ど変化しなくなると考えられる。
図 1(b) に仮想的にメモりバンド幅を変化させた時の影響を示す。ここでは 1 つの CPU 内に 1 プロセス× 1 ス
レッド (1P/CPU) で計算を行った場合と、4 プロセス× 1 スレッド (4P/CPU) で計算を行った場合を比較することで仮想的にメモリバンド幅が変化した場合の計算性能の変化を調べた。詳細に関しては 2.6.A を参照のこと。
1P/CPU は CPU のメモリバンド幅をほぼ占有できるが、4P/CPU は 4 プロセスでメモリバンド幅を共有するた3

め、1P/CPU のケースに比べて 1/4 のメモリバンド幅とみなせる。図より、ループ A はメモリバンド幅が減少すると 10%程度性能が下がるが、ループ B は 2%程度しか下がらず、この範囲ではメモリバンド幅にあまり影響を受けないことがわかる。
4 Threads

4 Threads

0.40.4
B (Elaps)
B (L2 MISS)

A (1P/CPU)A (4P/CPU)

0.3

0.3

0.2

0.25

0.1

0.2
020406080100
Block size

120140B (1P/CPU)
B (4P/CPU)

1.1
Elaps time/grid point [µ sec.]

0.35

L2 MISS [%]Elaps time/grid point [µ sec.]

A (Elaps)
A (L2 MISS)

1.2

10.9

0160

0.8
020(a) 経過時間と L2 ミス率

406080100
Block size

120140160(b) メモリバンド幅の影響

図 1: JSS におけるループ性能

2.6.3.2

インテル系 CPU での測定インテル系 CPU でも PC-S(Core i7-3960X) を中心にループ A および B の比較を行った。図 2 に 1 スレッドと6 スレッドの場合 (どちらも 1 プロセス実行) の測定結果を示す。横軸は格子ブロックサイズ (=ループ長)、縦軸は1 格子点あたりの計算時間 (経過時間) である。ここで、図中の 1333 および 1600 はメモリクロックの値を示し、大きい方 (1600) がデータ転送能力が高いため、ループ A、B ともメモリクロックが高い方が性能が高いことがわかる。
ループ A は格子ブロックサイズが 23 および 78 の前後で不連続な特性、また 1 スレッドのケースで格子ブロックサイズが 63 の時に局所的な性能悪化が見られる。キャッシュやメモリバンク競合などの原因が考えられるが詳細は不明である。一方ループ B は比較的素直な特性を示しており、一般的な傾向 (格子ブロックサイズが小さい範囲では演算器ネックのため演算量の増加に伴って計算時間が増加し、格子ブロックサイズが大きくなりキャッシュが溢れる様になるとメモリネックになり計算時間がほぼ一定となる) が見られる。スレッド並列の場合は格子ブロッ
クサイズが小さい範囲ではスレッド並列のオーバーヘッドが顕著になり、結果的に計算時間が増加していると考えられる。スレッド並列の場合、プログラム中の各ループはスレッド数で分割されるため、ループ長がスレッド数で割った余りに応じて計算時間が変動している様子が観察できる。
ループ A とループ B との比較では、スレッド数の増加および格子ブロックサイズの増加など、メモリアクセスの負荷が大きくなると、ループ A とループ B で計算性能の逆転現象が見られ、ループ B がループ A に比べて 18%程度良い性能を示している。この事は、それぞれのループでメモリ周波数の違いによる計算時間の差を見ても同様のことが言える。つまり、それぞれのループで 1333 と 1600 の違いを見ると、メモリ性能が低くなった場合に計算速度がどの程度悪化するかがわかる。ループ B はメモリ性能が悪化してもほとんど計算時間が悪化していないが、ループ A はメモリ性能の悪化に対して計算時間が 7%程度 (6 スレッドの場合) 増加している。この結果からもループ B はループ A に比べてメモリ性能にあまり依存しないと考えられる。

41 Thread

6 Threads

0.40.15
A (1333)
B (1333)
A (1600)
B (1600)

0.14

0.38

0.13
Elaps time/grid point [µ sec.]

Elaps time/grid point [µ sec.]

0.36

0.34

0.32

0.12

0.11

0.1

0.3
0.09
A (1333)
B (1333)
A (1600)
B (1600)

0.28

0.08

0.26

0.07
020406080
100
Block size

1201401600(a) 1 スレッド

20406080
100
Block size

120140160(b) 6 スレッド

図 2: PC-S(Core i7-3960X) におけるループ性能

2.6.3.3

B/F による整理

理論メモリ性能 (Byte/s) と理論演算性能 (FLOPS) の比を B/F と呼ぶが、アプリケーションの特性を議論する際に重要な指標となる。一般に圧縮性流体解析プログラムは高い B/F が必要と言われている。1CPU で実行するプロセス数を変化させることで、1 プロセス当たりのメモリバンド幅を仮想的に変化させることを考えた。この手法を用いてそれぞれのループのメモリバンド幅が減少した際の性能特性を調査した。メモリバンド幅を仮想的に変化させる手法の詳細については 2.6.A を参照の事。なお、2.6.A の結果より、JSS および Core i7-965 は比較的この手法の精度が良いと考えられる。
測定結果を図 3 に示す。この図で「3960X」、「965」は表 1 の PC-S (Core i7-3960X)、PC-N (Core i7-965) を示す。また、「Intel」、「GNU」はそれぞれインテルコンパイラー、GNU コンパイラーを示す。
格子ブロックサイズは 80 である。横軸は B/F であるが、ここでは理論性能の値を用いた。理論性能に対して実際に出る性能 (実行性能) は計算機システムによって異なるため、異る CPU の結果を比較する際は注意が必要である。
図 3(a) の縦軸は「性能」を示す。ここで「性能」は単位理論性能あたりの計算速度 (経過時間の逆数) とした。

JSS を始めとして全ての CPU でループ A の方が良い性能を示している。しかしながら、B/F が減少するとループ A は急速に性能が悪化している。一方、ループ B は B/F の減少にともなう性能悪化はそれほど酷くないことがわかる。1 ケースだけではあるが、B/F が最低の時にループ B の方が速い結果が得られている。更に B/F が減少した場合に 2 つのループの特性がどうなるかは現時点では不明ではあるが、低 B/F の領域ではループ B の方が良い性能を示す可能性がある。
図 3(b) の縦軸は「計算効率」を示す。ここで、「計算効率」は前述の「性能」を最大 B/F の時の「性能」で正規化した値である。つまりある CPU に対して B/F が最大の場合 (1CPU に対して 1 プロセスを実行し、メモリバンド幅を占有した場合) の「性能」を 100 とした時の「性能」比である。この図からも、B/F が減少した時にループ A はループ B に比べて急激に性能が悪化していることがわかる。

51 Thread

1Thread, Block size=80

0.2

100

Performance (1./Elaps time/GFLOPS)

0.18

0.16

Performance / Performance@(Band_Width)max [%]A (3960X Intel)
B (3960X Intel)
A (3960X GNU)B (3960X GNU)
A (JSS)
B (JSS)A (965 GNU)B (965 GNU)

0.14

0.12

0.1

0.08

0.06

959085A (3960X Intel)
B (3960X Intel)
A (JSS)
B (JSS)A (965 GNU)B (965 GNU)

8075
00.5

11.522.53
Byte/Flops

3.544.5

500.51(a) 計算性能

1.522.5
Byte/Flops

33.54(b) 計算効率

図 3: B/F の変化による計算性能の変化

2.6.4

まとめ

UPACS の主要部を切り出したカーネルプログラムを作成し、2 種類のループ構造に対して計算性能の比較を行った。従来のループ構造よりも局所性を意識し空間スイープを極力減らしたループは、狙い通りキャッシュミス率を従来のループ構造に比べてほぼ半減させることができた。計算性能を支配する他の最適化要因が複雑に関係するため、新しいループ構造が常に良いという結果は得られなかったが、低メモリバンド幅のシステムでは有利になる可能性が示された。今回の測定では低 B/F 領域を十分に設定することができなかったため、更なる低 B/F領域での測定を実施する予定である。更には、粘性項、陰解法を含めた評価や、ループ構造だけではなく、キャッシュの有効利用など低 B/F を前提としたプログラミングモデルの検討を今後進めていく予定である。

2.6.A

STREAM と B/F の制御メモり性能を測定するベンチマークテストである STREAM を用いて、今回性能測定を行った CPU のメモリ性能を測定した。STREAM では COPY、SCALE、ADD、TRIAD の性能を測定できるが、ここでは COPY とTRIAD の性能をブロックサイズ N を変化させて測定した。ここで COPY は
do i=1,Nc(i) = a(i)enddoとなる。図 4 に測定結果を示す。
ブロックサイズが小さい場合はキャッシュを有効に活用することができ、一般的にメモリ性能が良いことがわかる。最大のブロックサイズ（4 × 107 ）での COPY 性能は Core i7-3960X が 25.3GB/s、Core i7-965 が 14.1GB/s、JSS が 9.35GB/s となった。ちなみに、この性能は OpenMP を用いて CPU 内の全コアを使ったスレッド並列での性能である。それぞれの理論メモリバンド幅は 51.2GB/s、31.2GB/s、40GB/s であるため、実行効率は 49.4%、45.2%、23.4%となる。
それぞれの CPU でコア数分プロセス× 1 スレッド、1 プロセス× 1 スレッド、1 プロセス×コア数分スレッドの比較を行った。ちなみに 1 プロセス×コア数分スレッドはマルチコア CPU を利用する際の標準的な手法であ6

STREAM COPY

Performance [GB/s]

1000

100

25.3
14.1
9.35

10

110000100000

1e+06

1e+07

1e+08

Block size
965 GNU (4P1T)
965 GNU (1P1T)
965 GNU (1P4T)
3960X GNU (6P1T)
3960X GNU (1P1T)
3960X GNU (1P6T)

3960X Intel (6P1T)
3960X Intel (1P1T)
3960X Intel (1P6T)
JSS (4P1T)
JSS (1P1T)
JSS (1P4T)

図 4: STREAM を用いたメモリ性能の測定結果り、ハイブリッド並列の基本となる。ここで注目したいのは 1 プロセス× 1 スレッドと 1 プロセス×コア数分スレッドとの比較である。特に JSS においては、1 プロセス× 1 スレッド (8.22GB/s) と 1 プロセス× 4 スレッド(9.35GB/s) の COPY 性能の差が小さい。これは 1 プロセス× 1 スレッドの場合は 1 スレッドが CPU のメモリバンド幅をほぼ使い切ることができることを意味する。一方、4 プロセス× 1 スレッドはメモリバンド幅の観点から
は 1 プロセス× 1 スレッドの場合と比べて 4 倍のメモリアクセスが発生するため、相対的に 4 分の 1 のメモリバンド幅と考えることができる。実際、4 プロセス× 1 スレッドでは 2.39GB/s となり、1 プロセス× 1 スレッドに比べて 1/3.44 倍である。この結果を踏まえて以下の様に考える。JSS の場合は、CPU あたり、4 コアを有し、1コアの演算性能は 10[GFLOPS]、CPU 全体でのメモリバンド幅は 40[GB/s] である。ここで 1CPU に 1 プロセス× 1 スレッドを実行すると、1 コアを使った演算であるため演算性能は 10[GFLOPS]、また CPU のメモリバンド幅をほぼ占有できると考えると、メモリバンド幅は 40[GB/s] となるため、B/F は 40/10 = 4 と考えられる。次に1CPU に対して複数のプロセス (各プロセスは 1 スレッド) を実行する。各プロセスは 1 コアで実行されるので演算性能はプロセス数がコア数を越えない範囲では常に一定で、プロセス当たりは 10[GFLOPS] となる。一方でプロセス間でメモリバンド幅を共有することになるため、プロセス当たりのメモリバンド幅は大体プロセス数分の1 と考えられる。つまり、2 プロセスであれば 1/2、4 プロセスであれば 1/4 と考えられるので、B/F はそれぞれ40/2/10 = 2、40/4/10 = 1 と考えられる。複数プロセスの実行によるオーバーヘッドの影響など、厳密にはこの通りにはならないが、この方法で B/F が変化した時の性能特性の傾向を見ることは可能と考える。
インテル系の CPU では、1 プロセス× 1 スレッドと 1 プロセス×コア数分スレッドとの性能差は、JSS に比べると大きく Core i7-3960X は 14.0GB/s 対 25.3GB/s、Core i7-965 は 11.1GB/s 対 14.1GB/s となり、特に Corei7-3960X は 1CPU に 1 プロセス× 1 スレッドの場合に CPU のメモリバンド幅を占有できるとは言えない。また、
コア数分プロセス× 1 スレッドと、1 プロセス× 1 スレッドの比は Core i7-3960X が 14.0GB/s 対 4.23GB/s で1/3.31(6 コア)、Core i7-965 が 11.1GB/s 対 3.58GB/s で 1/3.10(4 コア) となり、Core i7-3960X は想定からの乖離が大きい。
以上の結果より、CPU 内のプロセス数を変化させて仮想的に B/F を変化させるという方法は、JSS および Corei7-965 に対しては定量的にもほぼ適用できると思われる。しかしながら、Core i7-3960X ではこの考え方は定量的には問題があるが、定性的な議論には利用できると思われる。

7