1
hp150246 「京」共用法第 12 条調査研究 Investigation Study Defined by the Public Utilization Act
高次精度離散化手法に関する調査研究
Performance Evaluation of High-Order Discretization Methods
富山 栄治Eiji TOMIYAMA
一般財団法人 高度情報科学技術研究機構
Research Organization for Information Science and Technology 要旨 本研究では，「京」を利用する研究課題の促進のため，コンピュータ・シミュレーションにて用 いられる有限要素法などにおいてホットスポットとなる疎行列ベクトル積に着目し，性能向上に 寄与すると期待される高次精度離散化手法の調査を実施した．高次精度化が行列を密にすること を踏まえ評価用プログラムを作成し，係数行列のデータサイズや連続性に生じる変化が実効性能 をどのように向上させるかを調査した． この結果，行列ベクトル積の内側ループ長の増大によって単体性能が最大で1.4 倍に向上する ことが明らかとなった．また，データの連続性増大によって1.5 倍に高速化することも示した． これらから，高次精度離散化手法は「京」コンピュータの特徴であるSIMD 性能とメモリバンド 幅を効率的に活用できる手法であることが分かった．一方で，高次精度離散化法を適用しても， 8 コア実行時にはメモリバンド幅が不足し性能が律速されることがわかった． キーワード：有限要素法，高次精度離散化，性能分析，演算性能，データの局所性 Abstract
In this project, we made a test code of the sparse matrix-vector product as a hot spot of the finite element method in order to investigate an impact of the density of sparse matrices in the high-order discretization method on the performance. Using this approach, we investigated the performance that is changed by the density of sparse matrices, the size of them, and their continuity/discreteness.
It was confirmed that the floating-point performance can be increased up to 1.4 times by expanding the inner most loop length. Also, the performance with higher continuity of the data became 1.5 times at the maximum. These results show that the higher-order discretization method can efficiently utilize the SIMD operation and memory bandwidth that are characteristics of the K computer. On the other hand, it was found that the performance is limited by the lack of memory bandwidth when the program was run with eight threads.
© 2021 Research Organization for Information Science and Technology All rights reserved. Received: 24 March 2020
Accepted: 2 April 2021 Available online: 28 April 2021
2
Keywords: Finite element method, high-order discretization, performance analysis, floating-point
performance, data locality
1. 研究の背景と目的 近年のプロセッサでは演算性能に対するメモリデータ転送性能が相対的に低下してきており， プログラムの実効性能がメモリ性能に律速される傾向が強まっている．計算機における高い実効 性能の実現にはプロセッサの理論ピーク演算性能（Flop/s）とメモリからのデータ転送性能 （Bytes/s）のバランスが重要であり，高い実効性能を示した事で知られる地球シミュレータでは メモリ転送性能とピーク演算性能の比Bytes/Flop（B/F 値）が 4.0 と非常に高かった[1,2]．しかし 「京」コンピュータにおけるB/F 値は 0.5 であるのに対し「富岳」では 0.30 ~ 0.33[3]，同様に富 士通 HPC システムでは FX10 は 0.36，FX100 は 0.48[4]，「富岳」と同じ CPU を搭載している FX1000 では 0.30[5]となっている．高次精度離散化法は，1 つのステンシル（要素，セル）内に多 数の自由度を持たせることで高精度近似を可能とし，かつデータと演算の局所性を高めることで プロセッサ演算性能のより一層の活用を可能とする事が期待されている[6]． そこで，本調査研究では，施設利用研究の促進のため，高次精度離散化法に関わる性能特性の 把握や，その手法を利用する技術を蓄積することを目的とする． 2. 計算モデル 2.1 高次精度離散化について コンピュータ・シミュレーションにおいては，物理現象をシミュレートするために様々な離 散化手法が用いられる．そして離散化の際，離散化の次数を調整することにより任意の計算精 度を得ることが出来る．しかしながら，従来は演算性能がボトルネックとなり高次精度離散化 法は利用されて来なかった．そして低次離散手法を用いてステンシル（要素，セル）の細分化 により精度を向上させ，かつ並列性を増加させる方向で対応していた．一方，近年のプロセッ サの傾向として演算性能に対するメモリデータ転送性能が相対的に低下しているためにプロ グラムの実行性能がメモリ性能に律速される傾向が強まっている状況から，従来のステンシル の細分化（空間分割法）による高精度化と性能向上がいきづまりつつある．そのため，高次精 度離散化法がデータと演算の局所性を高めプロセッサの演算性能のより一層の活用を可能と し，プログラムの実行性能を向上させる１つの有力な手法として期待されている． 高次精度離散化法がどのようにプログラムの実行性能を向上させるか，その機構の概要を以 下に示す． ある一つのステンシルを基準に高次精度離散化法と空間分割法で高精度化する場合を考え る．高次精度離散化法と空間分割法では，プログラムの実効性能を支配するデータと演算の局 所性に関して以下のような特徴がある．
3 高次精度離散化法では，1 つのステンシル内で節点を設け自由度を表し，それらを内挿点と した内挿関数を設定し自由度を増加させることで高精度近似を可能とする．例えば有限要素法 ２次元四辺形メッシュにおける次数の高次化は図1 のようになる． 図1 2 次元四辺形メッシュにおける 1 次要素（左），2 次要素（右） 1 次要素では，ある節点に接する要素に含まれる節点数は 9 であるのに対し，高次化された 2 次要素においては 21 となる．その結果，係数行列の計算においては高次関数の計算と積分点 数が増加する．そのときステンシル当たりの節点数（必要なデータ移動量に相当）に対する演 算量のオーダーの比は急激に増加し，演算が要求するB/F 値は減少することになる．また，ア プリケーションにおいてホットスポットとなることが多い疎行列ベクトル積においては，ステ ンンシル内の節点がステンシルの接続を通して関係する節点数は次数の増加とともに急増す る．その結果，各行の非ゼロ要素が増加するためプログラムの最内側ループ長が増加し，ベク トル長とデータの連続アクセス数が増加する．さらに，外側ループにおける変化としては，一 度ロードしたx ベクトルの再利用率が増加することによって B/F 値を低下させることが可能と なる．一方，空間分割法では要素のサイズを小さく変更した場合でも節点数と演算量の比には 変化が無い． 近年のシミュレーションにおいては1 次要素または 2 次要素を用いて離散化されることが多 いが，それ以上の高次精度離散化によってさらに節点数が増え（表1），演算負荷の割合が大き くなることで「京」を始めとする HPC に搭載されているプロセッサの性能をより効率的に活 かす事が期待出来る． 表1 高次化次数と演算に関わる節点数の対応(3 次元) 次数 1 2 3 4 5 6 節点数 27 81 135 189 243 297
4 2.2 行列ベクトル積計算について 有限要素法，有限差分法，有限体積法といった代表的な離散化手法においては，離散化の結 果いずれも連立一次方程式Ax = b を解く事に帰着する．ここで係数行列 A は疎行列となる． 計算の過程としては係数行列A の作成と行列ベクトル積の計算があるが，本研究では一般に計 算のホットスポットとなる行列ベクトル積に着目し高次精度離散化法の性能特性を調査した． そして，SIMD 性能とバンド幅性能に優れる「京」コンピュータの演算性能及び近年のプロセ ッサの演算性能を引き出すことの有効性を調査した． 高次精度化に伴う演算性能への影響を分析するためのプログラムとして，3 次元格子の場合 に対応する行列ベクトル積計算プログラムを作成した．プログラムでは，表1 に示した高次精 度化に相当する節点数の変化を再現出来るように内側非ゼロ要素ループのループ長を可変に するとともに，行列サイズやリストベクトルの連続性も制御可能にした．このプログラムは， 仮想的に高次精度離散化法のプログラムを模擬したものであるが，その本質を捉えたものとな っている． 行列ベクトル積計算部分のループ構造を図2 に示す． do j=1,n ！節点ループ b1 = 0.d0 b2 = 0.d0 b3 = 0.d0 do k=list(j-1)+1, list(j) ！非ゼロ要素ループ i=ia(k) x1 = x(3*i-2) x2 = x(3*i-1) x3 = x(3*i )
b1= b1 + a(1,1,k)*x1 + a(1,2,k)*x2 + a(1,3,k)*x3 b2= b2 + a(2,1,k)*x1 + a(2,2,k)*x2 + a(2,3,k)*x3 b3= b3 + a(3,1,k)*x1 + a(3,2,k)*x2 + a(3,3,k)*x3 end do b(3*j-2) = b1 b(3*j-1) = b2 b(3*j ) = b3 end do 図2 行列ベクトル積計算ループ 3. 並列計算の方法と効果（性能） 本研究においては自動並列によるスレッド並列化を適用した．コンパイル時のオプションは， ”-Kfast -Kparallel”を指定した．”-”-Kfast”は「京」における推奨最適化オプションであり，”-O3”をは じめとして複数の最適化オプションが誘導される．”-Kparallel”は自動並列によるスレッド並列化 を有効にするオプションである．外側の節点ループをブロック分割することで，「京」においてコ
5 ア数8 を用いたスレッド並列実行時においても充分なロードバランスとなった．8 スレッド実行 時の性能を計測したところ，演算の実効効率（ピーク性能比）は9.2%，演算で使用されたメモリ バンド幅はほぼ実効的な上限となる52 GB/s（理論性能の 80%）となった．またこの時のキャッ シュ効率はL1 キャッシュミス率 5.3 %，L2 ミス率 4.7 %であった． 次に，この演算性能の妥当性について以下に検討する．プログラムがメモリに対して要求する B/F 値と計算機の実効的な B/F 値を比較することでその計算機におけるおおよその性能上限を見 積もる事ができる[7]．本ループにおける要求 B/F 値を図２に基づき，ループがある程度回転しキ ャッシュの利用が十分になされている状態について算出する．最内のk ループ内にて，整数配列 ia および一次元配列 x はサイズが小さいためキャッシュに維持される．倍精度実数配列 a は再利 用性がないためメモリからのロードが9 要素発生する．演算子は 18 個あるため B/F = 9*8/18 = 4 である．一方，「京」の実効B/F 値は 0.36[7]であるので，0.36/4 = 0.09 となり性能上限は 9 %と見 積もる事ができる．従って，今回作成したテストプログラムは「京」の性能を充分引き出せてい ることがわかる． 次にキャッシュミス率について検討する．「京」における再利用性のない倍精度実数に連続アク セスする場合のL1 キャッシュミス率は，キャッシュのブロックロードが 128 Byte であることか ら128/8 = 16 となり必ず 16 要素に 1 要素のキャッシュミスが発生することとなり 1/16 = 0.0625 すなわち6.25 %である．行列ベクトル積においてはデータアクセスの多くは係数行列に対する再 利用性のない連続アクセスであることから，L1 キャッシュミス率 5.3 %は問題の無い値と判断し た．L2 キャッシュミス率 4.7 %についても，配列 a に対する L1 キャッシュミスは殆どが L2 キャ ッシュにおいてもミスとなると想定されるため L1 キャッシュミス率と同程度のミス率となるの が妥当と判断した． また，メモリに対する要求B/F 値が実効 B/F 値を大きく超えているために本計算はメモリデー タ転送性能によって律速されていると考えられ，このような計算においてはL1，L2 キャッシュ のミス率を低減させる最適化を行っても性能向上が難しい事を意味する． 4. 研究成果 「京」においてCPU 内 8 コア全てを用いて，L2 キャッシュサイズ 6MB を十分に超える問題 規模である節点数65,536 の条件で性能評価を実施した（図 3 左）．その結果、連続アクセスデー タ（図3 左の赤丸）およびランダムアクセスデータ（図 3 左の青三角）いずれにおいても横軸の
inner loop length（内側ループ長）の増大に伴う性能向上を確認した．内側ループ長の増大は高次 精度化に対応している．
6 図3 内側ループ長・アクセスパターンと実効効率の関係 8 コア実行（左），1 コア実行（右） 8 コア実行と同様に内側ループ長の増大によって性能が向上するとともに，8 コアの場合より も実効効率が向上していることを確認した．これは，使用コア数が減ったことで相対的にメモリ バンド幅が増大した効果によると考えられる．最も効果の高かった1 コア，連続アクセスの条件 において，内側ループ長256 での実効効率は 8 の場合と比較して 1.4 倍に高速化された（図 3 右 の赤矢印）． データアクセスの連続性による効果については，8 コア実行の場合に連続アクセスではランダ ムアクセスと比較して実効効率が1.5 倍程度になった（図 3 左，）．高次精度離散化手法において は，低次数でのメモリアクセスはランダム性が高く，高次数になると近接した節点が増える事で メモリアクセスの連続性が高い状態になると考えられるため，内側ループ長の増大とメモリアク セスの連続性向上の二つの効果により大きな性能向上が期待できる． 次に，問題規模による実効効率への影響を検証するため，8 コア実行・連続アクセスパターン において問題規模を変化させて調査した結果，内側ループ長が短い，即ち離散化の次数が小さい 時には問題規模が大きい程性能が高くなった．一方，内側ループ長が長い場合には問題規模の大 きさに依らず同じ性能を示した（図4）．これは，問題規模が小さい場合に次数の小さな離散化で はメモリバンド幅を充分に活用出来ていないのに対して，高次精度化することによってデータ転 送性能を向上させ演算性能を引き出せた事を示している．
7 図4 内側ループ長・問題規模と実効効率の関係 本調査により，高次精度離散化手法を用いる事によって「京」を含む計算機における実効効率 の向上が期待できることがわかった． 5. まとめと今後の課題 「京」における本調査研究の結果，疎行列ベクトル積の内側ループ長の増大によって，1 コア 実行・連続アクセスの場合に単体性能が約1.4 倍に向上することを確認した．また，利用データ の連続性が高いほど性能が向上することも確認した．これは，「京」の特徴である SIMD 演算機 能などを効率的に活用した演算が実現できるためであると考えられる．これらの効果は有限要素 法などにおいて高次精度離散化を適用することによって実現でき，現在最先端のHPC である「富 岳」やFX1000 のような「京」よりも小さな B/F 値を持つシステムを用いた場合にはより高い性 能向上が期待できる． 本研究は疎行列ベクトル積を対象としており，これは様々な数値解法に基づく実アプリケーシ ョンにおいても計算時間の大部分を占めるため，本研究で得られた知見は多くの実アプリケーシ ョンにおいても有効であると期待出来る．しかし，実アプリケーションにおいては係数行列作成 処理が必要であり，この処理は高次精度離散化によって負荷が高まると推測されるため，全体と してどの程度の効果となるかについては注意が必要である．また，既存のアプリケーションを高 次精度化しようとする場合には，数値安定性の問題などが生じないよう高次精度化する必要があ り実装は簡単ではないという問題もある．また，本研究では単体性能における評価であり，実ア プリケーションを並列実行する場合に経過時間にどの程度効果があるかについては明らかでは ない．今後，アルゴリズム全体，すなわち係数行列作成，疎行列の線形計算を通した全体の性能， また有限要素法，有限体積法等の具体的なアプリケーションに調査対象を広げて高次精度離散化 法の有効活用に向けた調査を進めたい．
8 参考文献
[1] L. Oliker, A. Canning, J. Carter, J. Shalf and S. Ethier, "Scientific Computations on Modern Parallel Vector Systems," SC '04: Proceedings of the 2004 ACM/IEEE Conference on Supercomputing, Pittsburgh, PA, USA, 2004, pp. 10-10.
[2] L. Oliker et al., "Leading Computational Methods on Scalar and Vector HEC Platforms," SC '05: Proceedings of the 2005 ACM/IEEE Conference on Supercomputing, Seattle, WA, USA, 2005, pp. 62-62. [3] https://www.r-ccs.riken.jp/fugaku/system/
[4] https://accc.riken.jp/wp-content/uploads/2015/06/chiba.pdf
[5] https://www.fujitsu.com/downloads/JP/jsuper/primehpc-fx1000-datasheet-ja.pdf
[6] V. Dobrev, T. Kolev, R. Rieben, B. Still, “Scalable multi-physics simulations will require new discretization and numerical methods research” DOE Workshop on Applied Mathematics Research for Exascale Computing, No.17 (2013).
[7] 南 一生, 井上 俊介, 堤 重信, 前田 拓人, 長谷川 幸弘, 黒田 明義, 寺井 優晃, 横川 三津
夫.：”「京」コンピュータにおける疎行列とベクトル積の性能チューニングと性能評価”ハイパ
