第9回 配信講義　計算科学技術特論A （2023）

行列計算における高速アルゴリズム ― エクサフロップス時代における 線形計算アルゴリズムの課題と研究動向 ― 2023年6月15日 計算科学技術特論A 電気通信大学 情報理工学研究科 情報・ネットワーク工学専攻 山本 有作

本講義の構成 • 第1回（6月8日） 「大規模連立1次方程式の反復解法」 – 代表的な行列計算である連立1次方程式の解法について，現在最も 広く利用されているクリロフ部分空間法の基礎を紹介する • 第2回（6月15日） 「エクサフロップス時代における線形計算アルゴリズムの課 題と研究動向」 – 様々な行列計算について，計算機上での実装における課題と，その 解決のための取り組みを紹介する – また，量子アニーリングマシンの行列計算への適用可能性について も簡単に触れる 2

はじめに： スーパーコンピュータの性能動向 https://samhatfield.co.uk/2019/04/10/supercomputer-trend-data/ 3

本講義の目的 • 本講義の目的 – エクサフロップス時代における数値計算アルゴリズムの課題を， 線形計算に焦点を当てて考える – 最近の研究動向について，簡単なサーベイを行う • 目次 – エクサフロップスマシンのハードウェア特性 – 線形計算アルゴリズムの技術課題 – 最近の研究動向 4

ポストペタ時代のスーパーコンピュータ： 4つのタイプ • 汎用型 「京」の路線 – メモリ容量・帯域・演算性能をバランス良く向上 – 「京」のように汎用的に様々な問題に適用可能 • 容量・帯域重視型 ベクトル計算機タイプ – 汎用型から演算性能を落として，メモリ性能により多くの資源を割く • 演算重視型 メニーコアタイプ – メモリ性能を落とし，演算性能により多くの資源を割く • メモリ容量削減型 – メモリ容量を極限まで削減し，オンチップメモリですべての計算を完結 出典： HPCI技術ロードマップ白書 5

各タイプの性能諸元（予測値） • 条件 – 「京」と同程度の消費電力（20MW） – 「京」と同程度の設置面積（2000～3000m2） 最も総演算性能が高いのは演算重視型 「富岳」は汎用型 出典： HPCI技術ロードマップ白書 6

エクサフロップスマシンの特徴（1） • 多階層の並列性 – 107 ～109 個程度の演算コア（動作周波数はGHzオーダー） – コアレベル，チップレベル，ノードレベル，システムレベルの並列性 • 複雑なメモリ階層 – – – – – コア毎のレジスタ コア毎のキャッシュメモリ オンチップメモリ ノード内共有メモリ 他ノードのメモリ 1～10チップ / ノード 共有メモリ 102～103コア / チップ メモリ メモリ 105～106ノード / システム 7

エクサフロップスマシンの特徴（2） • データ移動コストの増大： スループット – 主メモリのByte/Flop（データ転送性能と演算性能の比）による比較 総演算性能 総メモリ帯域 比 演算重視型 1000～2000PFLOPS 5～10PB/s 0.005Byte/Flop 富岳*1) 400PFLOPS 160PB/s 0.37Byte/Flop 京 10PFLOPS 5PB/s 0.5Byte/Flop 主メモリアクセスの相対的コストは徐々に増加 レジスタ 演算コア データ転送速度大 キャッシュ データ転送速度小 主メモリ *1) https://www.v-t.co.jp/product/primehpc-a64fx/ 8

エクサフロップスマシンの特徴（3） • データ移動コストの増大： レイテンシ（遅延時間） – コア間同期・通信レイテンシ： 100サイクル – ノード間通信レイテンシ： 80～200サイクル AllReduce のような全ノード同期型通信で特に影響大 80～200サイクル 9

エクサフロップスマシンの特徴（4） • 電力消費の問題 – 発熱抑制と節電の両面から，電力消費の抑制が重要 – オフチップのデータアクセスが大きな電力を消費 • 倍精度演算： 1.1pJ/FLOP • オンチップデータアクセス： 2.1pJ/Word • オフチップデータアクセス： 205pJ/Word 100倍 • 部品数増加に伴う故障確率上昇 – ハードエラー（熱などによる部品の故障） – ソフトエラー（α線などによるビット反転） ハード／システムソフトレベルでのエラー耐性強化は可能だが， 部品数・速度・消費電力の面で大きなコストが必要 出典： HPCI技術ロードマップ白書 10

アプリケーション並列化における要求の変化 • これまでの並列化研究： 弱スケーリングに重点 – ノード当たりの問題サイズが一定という条件下で並列化 – ノード数を増やすとともに，問題サイズもどんどん大きくする – 比較的，並列性能を出しやすい • エクサフロップス時代における問題点 – すべてのアプリで問題サイズ拡大が求められているわけではない • 問題サイズは固定し，モデルを精密化したい場合 • 問題サイズは固定し，時間ステップを増やしたい場合（分子動力学など） – 弱スケーリングでは，ノード数増加につれて計算時間も増加 • 多くの場合，演算量は問題サイズに対して線形より速く増加 • エクサフロップスマシンの性能を引き出せる問題サイズでは，計算時間 が長くなり過ぎて実用的でなくなる場合も 11

線形計算アルゴリズムの課題(1) • 107～109個のコアを活用できる並列性 • データ移動の削減（I）： データ移動量の削減 – データ移動量（階層間／ノード間）に比例してかかるコストの影響を削減 データ再利用性の向上 データが上位のメモリにある間に できるだけ集中して演算を行う 消費電力節減の面からも有効 演算コア キャッシュ 主メモリ • データ移動の削減（II）： データ移動回数の削減 – データ移動／同期1回ごとにかかるコストの影響を削減 計算粒度 計算粒度の大きいアルゴリズム 同期やデータ移動をできるだけまとめて行う 通信と通信の間で，できるだけ多くの演算を行う 12

線形計算アルゴリズムの課題(2) • アルゴリズムレベルでの耐故障性 – 結果不正のノード，結果を返さないノードがあっても，計算が破綻せ ずに進行 • 計算量のオーダーの低減 – ある程度の（確率的）誤差を許容することでオーダーを低減 テンソルの近似に基づくアルゴリズム 確率的アルゴリズム • 強スケーリングの意味で効率的なアルゴリズム – 問題サイズを固定したとき，ノード数を増やすほど実行時間が短縮 強スケーリングの条件下では，通信・同期時間が支配的 大粒度のアルゴリズムはこの場合にも有効 13

109個のコアを活用できる並列性： 密行列の場合 • 正方行列に対するアルゴリズム – LU分解，QR分解，3重対角化，Hessenberg化 並列性は十分 – 全体の計算量： O(n3)，並列性： O(n2) – ネックとなりうる部分： ピボット生成，通信 これらを隠蔽するためのスケジューリングが重要 • DAGスケジューリング • 帯行列・縦長行列に対するアルゴリズム b b – 帯行列のLU分解，帯行列の3重対角化， n – 縦長行列のQR分解，ベクトル逐次添加型の直交化 – 全体の計算量： O(n2b) or O(nb2) ，並列性： O(b2), O(nb) or O(n) 新たな並列性を導入できるアルゴリズムが必要 • • 帯行列に対する分割統治型のLU分解 Compact WY 表現に基づくベクトル逐次添加型直交化法 14

109個のコアを活用できる並列性： 疎行列の場合 • 部分空間への射影に基づく反復解法 – 連立1次方程式の解法： Krylov部分空間法全般 – 固有値問題の解法： Lanczos法，Arnoldi法，Jacobi-Davidson法 部分空間の拡大操作は本質的に逐次的 並列性： 行列ベクトル積の並列度 O(nz) • 並列性を向上させる手法 – ブロックKrylov部分空間法 • 複数の初期ベクトルから出発し，複数の部分空間を同時に生成 • 部分空間の数だけ並列性が向上 – 数値積分を用いたフィルタで部分空間を生成する手法 • 櫻井・杉浦法 • 積分点ごとの並列性を新たに利用可能 15

データ移動の削減： 密行列の場合 • 想定する状況と目標 コア – コア：1個，キャッシュ： M ワード，行列サイズ： n – キャッシュと主メモリとの間のデータ移動を最小化 • ノード間通信最小化でも同様の手法を適用可能 • ブロック化（タイル）アルゴリズム レジスタ キャッシュ データ移動 最小化 主メモリ (M/3)1/2 – 行列を，サイズ (M/3)1/2×(M/3)1/2 のブロックに分割 • ブロック3個がキャッシュに入る – 各ブロックを要素と見て行列計算を行う • 行列乗算 • LU分解，コレスキー分解，QR分解 • ブロック3重対角化，ブロック Hessenberg 化 n – 演算の主要部分はブロックどうしの乗算であり，これはキャッシュ上 で実行可能 16

ブロック化アルゴリズムの例 • ブロック化コレスキー分解 – ブロックサイズを L= (M/3)1/2，第 (I, J) ブロックを AIJ とする do K=1, n/L AKK := Cholesky(AKK) do I=K+1, n/L AIK := AIKAKK-T end do do J=K+1, n/L do I=J, n/L AIJ := AIK AJKT end do end do end do 対角ブロックのコレスキー分解 ブロックピボット列の作成 ブロックどうしの行列乗算 （演算の主要部） – LAPACKで採用されているアルゴリズム 17

ブロック化アルゴリズムの最適性 • 定理（Ballard, Demmel, Holtz and Schwartz, 2009） – ブロックサイズを (M/3)1/2 としたブロック化コレスキー分解は，前記の 仮定の下で，キャッシュと主メモリの間のデータ移動量をオーダーの 意味で最小化する • 証明 – 行列乗算についてはデータ移動量の下界がわかっていることに着目 – コレスキー分解を用いて行列乗算を計算するアルゴリズムを構築 – これより，定数倍の差を除いて， コレスキー分解における データ移動量の下界 ≧ 行列乗算における データ移動量の下界 が言える – ブロック化コレスキー分解が右辺の下界を達成することを示す 18

Communication optimal なアルゴリズム • データ移動回数の最小化 – 1回のデータ移動では，主メモリの連続した 領域しか持ってこられないと仮定 連続でないため，1回で 持ってこられない – このとき，通常の行列格納形式では，ブロッ ク化コレスキーは移動回数最小にならない – ブロック単位の格納順を使うことで，データ 移動回数も最小化 最適でない 最適 Communication optimal なアルゴリズム • データ移動量・回数の下界が知られているアルゴリズムの例 – – – – 行列乗算（O(n3)／Strassen） LU分解，コレスキー分解 QR分解，最小2乗法 固有値分解，特異値分解 両方について下界を達成する （Communication optimal な） アルゴリズムの開発が，現在 でも活発な研究テーマ 19

データ移動の削減： 疎行列の場合 • Krylov部分空間法 – Kk(A; b) = span{b, Ab, A2b, …, Ak–1b} の中で近似解を求めていく解法 – 演算の主要部分は疎行列ベクトル積 y = Ax – 1ステップ中に複数回の内積（またはノルム計算）が存在 • データ移動の観点からの問題点 – 行列ベクトル積は行列データの再利用性が低い • y = Ax において，A の各要素は1回しか計算に利用しない – 内積は全ノードでの AllReduce を必要とし，レイテンシの影響大 • 課題 – 行列ベクトル積におけるデータ再利用性の向上 – 複数の内積をまとめて同期回数を削減 以下ではGMRES法を例として手法を紹介 20

GMRES法 • 原理 – A をかける操作と直交化を繰り返し，部分空間を広げながら正規直 交基底 {q1, q2, q3, …}を生成 – 各ステップにおいて，||r(n)||2 = ||Ax(n) – b||2 が最小になるよう解を更新 • アルゴリズム（正規直交基底の生成部分） 新たなベクトルの生成（空間の拡張） 直交化 正規化 21

手法（I）： ブロック GMRES 法 • アイディア – 複数（b本）の初期ベクトル R(0) = [r1(0), …, rb(0)] から出発し，ブロック Krylov 部分空間 Kb(m)(A; R(0)) 内で解を求める • 普通の GMRES 法との比較（1ステップあたり） GMRES ブロックGMRES 演算量（行列ベクトル積） 1（y=Ax） b （Y=AX） 行列データのアクセス回数 1（y=Ax） 1 （Y=AX） 1 1 同期回数* • 効果 * 直交化に compact WY 表現またはCGS法を用いた場合 再利用性 b 倍向上 同期回数 実質 1 / b – 行列アクセス回数・同期回数は同じため，実行時間増加は b 倍以下 – ブロック Krylov 部分空間を使うことによる収束性向上の効果も存在 右辺が複数本の場合は有利．1本でも高速化できる場合もある 22

手法（II）： k-step GMRES 法 • アイディア – 行列ベクトル積 Ar(m), A2r(m) …, Akr(m) を一度に行って Krylov 部分空間 を一度に k 次元拡大し，その後に正規直交基底を生成する • 普通の GMRES 法との比較（1ステップあたり） GMRES k-step GMRES 演算量（行列ベクトル積） 1 1+a 行列データのアクセス回数 1 1/k 同期回数* 1 1/k • 効果 * 直交化に compact WY 表現またはCGS法を用いた場合 再利用性 k 倍向上 同期回数 1/k – 再利用性向上と同期の削減により，大きな性能向上が期待できる – ただし，（直交化前の）基底が線形従属に近い場合は，収束性が悪化 線形独立性を高めるため，A の単項式の代わりに直交多項式を 使い，3項間漸化式で計算するなどの工夫が検討されている 23

アルゴリズムレベルでの耐故障性 • 想定する状況と目標 – 複数のノードが通信を行いつつ協調して計算 – そのうち1個のが不正な結果を返すか，あるいは 結果を返さなくても，計算は破綻せずに進行 – その場合，精度劣化，収束性劣化は許容 – 計算の一部を，高信頼モード／高信頼ハード ウェア（ただし計算コスト大）で行ってもよい × • 考察 – 複数のノードの結果を集めて部分空間を改良するタイプのアルゴリ ズムは，耐故障性と親和性が高い – 大粒度並列性は，耐故障性にとっても有利 • 高信頼モードの使用頻度を削減できる 24

MERAM（Multiple Explicitly Restarted Arnoldi Method） • アルゴリズム（Emad，Petiton & Edjlali, 2005） (1) P 個のノードで，異なる初期ベクトルを用いて独立にArnoldi法を実行 (2) 各ノードで作ったKrylov部分空間を合わせ，大きな部分空間を生成 (3) その中で P 本の初期ベクトルを新たに生成し，Arnoldi法をリスタート • 耐故障性 – 上記(1)を低信頼モード，(2)，(3)を高信頼モードで実行 • 計算量の大部分は(1)のステップ – (1) において，結果不正のノード／結果を返さないノードがあっても， 収束性が落ちるだけで，計算は破綻しない ノード0 Arnoldi ノード1 Arnoldi ノード2 Arnoldi 低信頼モード 部分空間合成 P本のリスタート ベクトルの生成 Arnoldi Arnoldi Arnoldi 部分空間合成 P本のリスタート ベクトルの生成 Arnoldi Arnoldi Arnoldi 高信頼モード 25

耐故障性を持つその他のアルゴリズム • Fault-Tolerant GMRES法（Hoemmen & Heroux, 2011） – Flexible GMRES法の枠組みを利用 – 結果不正を，不適切な前処理と解釈して計算を続行 • そのステップでは無駄に部分空間が大きくなるが，精度には悪影響なし Im z • 櫻井・杉浦法 – ある積分点を担当するノードが結果を返さない場合， その点を使わない積分公式を新たに構築して使用 O × Re z • 密行列向けの解法 – 直接解法では，1箇所の計算間違いで結果に壊滅的な影響 – アルゴリズムレベルでの耐故障性の実現は，原理的に難しそう 26

確率的アルゴリズムによる計算量のオーダー低減 • 特異値分解： 行列の最良の低ランク近似を与える分解 – 画像処理，信号処理，情報検索など広い応用 – m×n 行列（m≧n）の場合の演算量は O(mn2) と大きい • CX分解： 特異値分解の代替手法 – C： A の列ベクトルを k 本（1≦ k ≦ n）選んでできる m×k 行列 – X： 適当な k×n 行列 – このとき，∥A – CX∥F をできるだけ小さくする C と X を求める A の特徴を最もよく表す s 本の列ベクトルを選ぶことに相当 Fノルムに対しては，X = C+A と選ぶのが最適 したがって，問題は C を選ぶことに帰着 27

Leverage score に基づく確率的アルゴリズム • Leverage score（Drineas et al., 2008） – A の打ち切り型特異値分解を Ak = Uk Sk VkT とする – このとき，次の量を，VkT の第 j 列の Leverage score と呼ぶ • サンプリングと誤差評価 – 確率 pj に従い，A の列をサンプリングする – このとき，確率 0.9 以上で次の誤差評価が成り立つ • 問題点 サンプリング 相対評価 回数に依存 – Leverage score の計算には，VkT が必要（このままでは非実用的） 28

Leverage score の近似 • 近似アルゴリズム（Clarkson et al., 2012） – O(mn log m) の計算量で，Leverage score を近似的に計算するアルゴ リズムが存在する • アイディア – A に左から適当な直交行列をかけることで，一様な Leverage score を 持つ行列 A’ に変換 • Johnson & Lindenstrauss の補題と高速 Walsh-Hadamard 変換を使う – A’ に対して一様なサンプリングを行い，CX 分解を求める – A’ のCX 分解から，A のCX 分解を求める 相対誤差の意味での（確率的）誤差評価を持つCX分解の計算が， O(mn log m) で可能に 29

強スケーリングの意味で効率的なアルゴリズム • 実対称行列の全固有値・固有ベクトル計算 – 分子軌道法をはじめ，様々な問題で重要 – n=10,000 程度の中規模問題に対する需要も多い • 想定する状況と目標 – 行列サイズは n=10,000 に固定 – ノードは何個使ってもよいから，できるだけ高速に解きたい • 分子軌道法では，多電子積分の計算量が O(n4) 以上で並列性も高い • この部分を並列化するため，1万ノード以上を確保している例もある • 標準的な行列計算ライブラリ（ScaLAPACK）での結果 – 「京」上では，ノード数を400以上にしても，計算が加速しない • 実行時間の70%以上を占める通信オーバーヘッドのため – 確保した1万ノードの大部分は，固有値計算部分では遊んでしまう 30

強スケーリングの意味で効率的なアルゴリズム • ブロックヤコビ法による全固有値・固有ベクトル計算 – 行列をブロックに分割し，直交変換により複数の非対角ブロックを並 列に消去 – これを繰り返すことで，行列を対角行列に近づける = • アルゴリズムの特性 – 演算量は3重対角化に基づく解法の10倍以上 – Pノードで2次元分割を行った場合，通信回数は O(P1/2 log P * Iter) n=P=10,000 ならば3重対角化の通信回数 O(N log P) に比べて小さい 通信OHの大きい状況では，3重対角化法より高速となる可能性 31

強スケーリングの意味で効率的なアルゴリズム • ScaLAPACK とブロックヤコビ法の性能比較 – 「京」上で1万ノードまでを使用して実行時間を測定 – 行列サイズは n=10,000 に固定 35 強スケーリングの条件下では ブロックヤコビ法が優位 消去順序の最適化，前処理に より，更に高速化の可能性 実 行 時 間 （ 秒 ） ScaLAPACK Block Jacobi 30 25 20 15 10 8.3s 5 5.5s 0 100 400 625 1600 2500 10000 ノード数 • 強スケーリングの条件下でのアルゴリズム設計 – 通信・同期オーバーヘッドの削減が重要 – そのためならば，演算量をある程度増やしてもよい • 中務氏の新しい固有値・特異値計算アルゴリズム 32

大粒度並列性を持つQR分解手法 • A ∈ Rm×n （m≧n）のQR分解 – full QR分解 A = QR A = Q A = Q R ~~ – thin QR分解 A = QR ~ ~ R • QR分解に対する2つの見方 ~ – Q に着目 → A の列ベクトルの正規直交化 – R に着目 → 直交行列による A の上三角化 様々な応用 • QR分解のアルゴリズム – ハウスホルダー法 – 古典／修正グラム・シュミット法 高性能計算の観点からは課題が多い 33

CholeskyQR2法による大粒度並列QR分解 • Cholesky QR 法による行列 X∈Rm×n （m≧n）のQR分解 – アルゴリズム （コレスキー分解： A = RTR） – 長所： 大粒度並列，すべての計算がブロック化可能 – 短所： X の条件数が大きいとき，直交行列 Y の直交性が悪化 • CholeskyQR2 法*1 – Cholesky QR 法で得られた Y を再直交化 • Y を X と見て，上記のアルゴリズムをもう一度繰り返す – X の条件数が 108 以下ならば，極めて精度の良い結果が得られる • 直交性 ||QTQ – I|| も残差 ||X – YR|| も丸め誤差レベル *1 Y. Yamamoto et al.: “Roundoff Error Analysis of the CholeskyQR2 Algorithm”, Electronic Transactions on Numerical Analysis, Vol. 44, pp. 306-326 (2015). 34

「京」上でのCholeskyQR2法の性能 ScaLAPACK (Householder QR) (sec.) m=4,194,304 n=16 x2.0 speedup over TSQR P: #processes (= #nodes) (sec.) TSQR m=4,194,304 n=64 x2.9 speedup over TSQR P: #processes (= #nodes) CholeskyQR2 (sec.) m=4,194,304 n=256 x7.6 speedup over TSQR P: #processes (= #nodes) • TSQRとCholeskyQR2はScaLAPACK（ハウスホルダーQR）より高速 （n=256でPが大きい場合を除く） • CholeskyQR2はTSQRよりも高速 （Pが大きい場合だけでなく，Pが小さい場合でも） コレスキーQR分解を2回繰り返してもまだTSQRより高速！ 35

CholeskyQR2法の応用 • 第一原理分子動力学法における波動関数の時間発展 • 時間離散化を行った場合 – Xt = [ψ 1(r, t), …, ψ N(r, t)] とおくと， ~ Xt+Dt = Xt +CDt ~ Xt+Dt = Xt+Dt R 時間発展 直交化（QR分解） ~ – Xt は直交化されているから，Dt が小さければ， Xt+Dt も直交に近い – すなわち，条件数は小さい（列が直交する場合，条件数は1） CholeskyQR2法が適用可能のはず 36

悪条件の行列のための shifted CholeskyQR3法 • Shifted CholeskyQR3 法*1 – 行列 X が悪条件の場合，Cholesky QR2 法において A = X TX が特異 に近くなり，コレスキー分解が破綻することがある – A に微小なシフト sI を加えることで，条件数 1016 程度まで対応可能 – 大粒度並列性など，HPC面から見た利点は Cholesky QR2 法と同じ *1 T. Fukaya et al.: “Shifted CholeskyQR for computing the QR factorization of ill-conditioned matrices”, SIAM Journal on Scientific Computing, Vol. 42, , No. 1, pp. A477-A503 (2020). 37

不完全LU分解型前処理の超並列化 • クリロフ部分空間法と前処理（前回の復習） – 連立1次方程式 K–1Ax = K–1b を考える – いま，行列 K が A の近似行列であるとする – このとき，方程式 K–1Ax = K–1b を考えると，その係数行列 K–1A は A と比べて単位行列に近くなっていると考えられる．これにより，クリロ フ部分空間法の収束はより速くなると期待できる – このような方程式の変形を前処理（左前処理）と呼ぶ • 不完全LU分解型前処理（ILU(0)前処理） – 行列 A のLU分解の際に，フィルインを無視して近似的に A≒LU と 分解し，K = LU とおく．これを ILU(0) 前処理と呼ぶ – ILU(0) 前処理は，幅広いクラスの行列に対して効果があるが，一般 に並列化は困難 • K–1 をかける計算で，前進消去と後退代入が必要となるため 38

ブロック赤黒順序付けによる並列化（1） • 概要 – 偏微分方程式を格子点上で離散化して得られる行列を考える – この場合，格子点の番号の付け替え（リオーダリング）により，ILU(0) 分解の並列性を抽出することが可能 – 代表的なリオーダリング手法として，赤黒ブロック順序付けがある 4×4 格子の ブロック分割数 x 方向: 2， y 方向: 2 の順序付けの例 ・ 青字： 自然な格子点の番号 ・ 赤黒字： リオーダリング後の格子点 の番号 岩下武史, 島崎眞昭: 同期点の少ない並列化 ICCG 法のためのブロック化赤-黒 順序付け 情報処理学会論文誌, Vol.43, No.4，pp.893-904，2002． 39

ブロック赤黒順序付けによる並列化（2） ~ • リオーダリングは係数行列 A の行と列の置換 A = PAPT に相当する ~ • ブロック赤黒順序付けによって得られる行列 A は下図のようになる • 同色のブロック間に依存関係がないため，各プロセスに1ブロックを担当 させることで，ILU(0) 分解と前進消去・後退代入が並列化できる 40

ブロック赤黒順序付けによる並列化（3） • ブロック数を増やした場合の収束性 – 120×120×60 格子でのポアソン方程式の求解 – 並列数 = 各色のブロック数 = 全ブロック数 / 2 全 ブロック数が10000程度までは，あまり反復回数が増えない 超並列向き 41

低精度計算の利用による高性能化（1） • 近年のCPU/GPUの低精度演算性能 富士通 A64FX（富岳） NVIDIA Quadro GV100 倍精度 3.072TFLOPS 7.40TFLOPS 単精度 6.144TFLOPS 14.80TFLOPS 半精度 12.288TFLOPS 29.60TFLOPS 低精度演算の活用により高速化を達成できる可能性あり • 前処理付きクリロフ部分空間法と低精度演算 – 前処理付きクリロフ部分空間法では，残差の正確性（r(n) = Ax(n) – b が どの程度正確に成り立つか）が最も重要 – この関係さえ保持しておけば，一部で低精度演算を使ったとしても， 解の精度は保たれる • ただし，収束性は低下する可能性あり 42

低精度計算の利用による高性能化（2） • 右前処理に低精度計算を使った BiCGSTAB 法1) 右前処理： AM y = b, x = M y 低精度計算による uk, sk の誤差 低精度演算 低精度演算 低精度計算による xk+1 の誤差 低精度計算による rk+1 の誤差 xk+1, rk+1 はそれぞれ誤差を含むが，関係式 rk+1 = Axk+1– b は保たれている 1) H. Tadano et al., “On Single Precision Preconditioners for Krylov Subspace Iterative Methods”, Proceedings of LSSC 2007, LNCS 4818, pp. 721-728, 2008. 43

低精度計算の利用による高性能化（3） • 右前処理に低精度計算を使った BiCGSTAB 法の数値結果 – 前処理としてブロック赤黒順序付けのMILU(0) 分解を使用 – ベクトルに前処理行列を適用する部分を単精度で計算 – 倍精度とほぼ同じ結果が得られている A. Shioya and Y. Yamamoto: “Block red–black MILU(0) preconditioner with relaxation on GPU”, Parallel Computing, Vol. 103, No. 1, pp. 102760-102772 (2021). 44

ブロック赤黒順序付けに基づくソルバのGPU性能 • ソルバの概要 – ブロック赤黒順序付けによる並列 MILU(0) 分解 + BiCGSTAB法 – コアレスアクセスのためにデータ格納形式を最適化 – 前処理部分で単精度計算を利用 • NVIDIA Tesla K40t 上での性能 – 例題： 3次元ポアソン方程式 59×59×29 119×119×59 239×239×119 45

おわりに • 本発表のまとめ – エクサフロップスマシンでは，多階層の超並列性，データ移動コスト の増大，故障確率の上昇などが課題となる．また，強スケーリングで の並列化効率がより重視される． – そのため，線形計算アルゴリズムの研究では次の点が重要になる． • 109レベルの並列性 • データ移動量・移動回数の削減 • アルゴリズムレベルでの耐故障性 • （確率的）近似による計算量のオーダーの削減 • 強スケーリング性 – これらの方向に沿った最近の研究をいくつか紹介した． 46