応用数理工学特論 線形計算と ハイパフォーマンスコンピューティング 応用数理工学特論 線形計算と ハイパフォーマンスコンピューティング 第7回 名古屋大学　計算理工学専攻 山本有作 日立製作所の山本有作です。 「～」について発表いたします。

目次 １. 固有値計算の概要 ２. 三重対角化処理の並列化 ３. 三重対角化処理のキャッシュ向け最適化 １.　固有値計算の概要 ２.　三重対角化処理の並列化 ３.　三重対角化処理のキャッシュ向け最適化 本発表では，はじめに，研究の背景を述べてから，スパースソルバの概要，並列化手法，そして本研究で工夫した点の一つであるRISCプロセッサ向けの最適化についてご説明します。最後に，並列計算機SR2201上での性能評価とまとめを述べます。

１. 固有値計算の概要 固有値計算とその応用 固有値・固有ベクトル計算の流れ 各部分の演算量 １.　固有値計算の概要 固有値計算とその応用 固有値・固有ベクトル計算の流れ 各部分の演算量 三重対角行列 T に対する固有値・固有ベクトル解法 はじめに，研究の背景として，有限要素法とスパースソルバの必要性についてご説明します。

固有値計算とその応用 対象とする問題 対象とする計算機 固有値計算の応用分野 標準固有値問題　Ax = λx　（A： 実対称またはエルミート行列） 対象とする計算機 共有メモリ型並列計算機 分散メモリ型並列計算機 SMPクラスタ 固有値計算の応用分野 分子計算，統計計算，構造解析など ここでは特に，Aが実対称またはエルミートの密行列の場合を考える。

固有値・固有ベクトル計算の流れ Q*AQ = T (Q: 直交行列) Tyi =λi yi xi = Qyi 実対称行列 A 計算内容 計算手法 三重対角化 Q*AQ = T (Q: 直交行列) ハウスホルダー法 三重対角行列 T QR法 DC法 二分法・逆反復法 Dhillon のｱﾙｺﾞﾘｽﾞﾑ 三重対角行列の 固有値・固有ﾍﾞｸﾄﾙ計算 Tyi =λi yi Tの固有値 {λi }， 　 固有ベクトル {yi } 密行列Aをまず三重対角行列Tに相似変換してからTの固有値・固有ベクトルを求めるのが最も一般的な計算法。 三重対角化には，後に述べるハウスホルダー法を使う場合がほとんど。 三重対角行列の固有値・固有ベクトルの計算には，色々なアルゴリズムがある。 xi = Qyi 逆変換 逆変換 Aの固有ベクトル {xi }

各部分の演算量 N×N の実対称行列 A に対し，M 組の固有値・固有ベクトルを求める場合 三重対角行列の固有値・固有ベクトルの計算：　O(NM) ～ O(N3) 　　（解法および固有値の分布により大きく異なる。） 逆変換：　2N2 M 　M≪N の場合は，三重対角化の演算量が支配的 　本報告ではこの部分に絞って並列化・最適化の手法を述べる。

三重対角行列 T に対する固有値・固有ベクトル解法 概要 計算量 並列化 一部の固有値 ・固有ﾍﾞｸﾄﾙ のみの計算 QR法 相似変換によりTを対角行列に変換 O(N2) + O(N3) 容易 不可 DC法 （分割統治法） Tを再帰的に半分のｻｲｽﾞの行列に分割 ～ O(N3) 二分法・ 逆反復法 二分法で固有値，逆反復法で固有ﾍﾞｸﾄﾙを計算 O(NM) + O(NM2) 非自明 可 Dhillon の アルゴリズム 逆反復法の改良 （dqds，RRR，twisted分解などの利用） + O(NM) ・三重対角行列の固有値・固有ベクトル計算の部分をもう少し詳しく見ると，次のようなアルゴリズムがある。本報告ではこの部分は主題ではないので，簡単にまとめてある。やや誤解を招く部分もある。 ・QR法が並列化容易というのは，O(N^3)の固有ベクトル計算部分のこと。また，一部の固有ベクトルのみを計算する手法もないことはない。 ・Dhillonのアルゴリズムは７年前のDhillonの博士論文で出てきたものであり，ここ数年大変注目され，論文もたくさん出てきている。 これについては，11/28の応用数理学会研究会で簡単なレビューを行う予定。 dqds = differential qd with shift RRR = Relatively Robust Representation; T をLLt の形で表現して計算を進める手法

２. 三重対角化処理の並列化 ハウスホルダー法による三重対角化 共有メモリ向けの並列化 サイクリック列分割による並列化 ２.　三重対角化処理の並列化 ハウスホルダー法による三重対角化 共有メモリ向けの並列化 サイクリック列分割による並列化 Scattered Square 分割による並列化 性能評価 はじめに，研究の背景として，有限要素法とスパースソルバの必要性についてご説明します。

ハウスホルダー法による三重対角化 鏡像変換による列の消去 第 k ステップでの処理 鏡像変換 H = I – a uut Hは対称な直交行列 与えられたベクトルの第１成分以外をゼロにする。 第 k ステップでの処理 ベクトル 左からH を乗算 左からH を乗算 右からH を乗算 k 非ゼロ要素 ゼロにしたい部分 影響を受ける部分

各ステップでの処理の詳細 鏡像変換ベクトルの作成 σ= sqrt(d td) σ= sqrt(d td) u = (d1 – sgn(d1)σ, d2, … , dN – k ) α= 2 /∥u∥2 HCH = (I – a uu t)C (I – a uu t) 　 　 = C – a uu tC – aC uu t + a2 uu tC uu t 　 　 = C – up t – pu t + (1/2) a uu tpu t + (1/2) a up tuu t　 　　( p ≡ aC u ) 　　 = C – uq t – qu t 　 　 　( q ≡ p – (1/2) a (p tu) u ) C d H による両側からの乗算 k 実際には，C の対称性を利用し， HCH は下三角部分のみ計算する。 p ≡ aC u の計算でも，C の上三角部分は下三角部分で代用する。

ハウスホルダー法のアルゴリズム [Step 1] k = 1からN –2まで以下の[Step 2] ～ [Step 8]を繰り返す。 　[Step 2]　内積 σ(k) = sqrt(d (k)t d (k)) の計算 　[Step 3]　鏡像変換ベクトル作成： 　　u(k) = (d (k)1 – sgn(d (k)1)σ(k), d (k)2, … , d (k)N – k ) 　[Step 4]　規格化定数の計算：　α(k) = 2 /∥u(k)∥2 　[Step 5]　行列ベクトル積：　p(k) =α(k)C (k)u(k) 　[Step 6]　ベクトルの内積：　β(k) =α(k) u(k)tp(k) /2 　[Step 7]　ベクトルの加算：　q(k) = p(k) –β(k)u(k) 　[Step 8]　行列のrank-2更新：　C (k) = C (k) – u(k)q(k)t – q(k)u(k)t ピボット列 ピボット行 d (k)1 C (k) d (k) k

ハウスホルダー法の特徴 演算量 データの再利用性 並列粒度 行列サイズが N のときの演算量： 約 (4/3)N3 rank-2更新の演算量： 約 (2/3)N3 データの再利用性 演算はほとんどBLAS2のため，行列データの再利用性なし キャッシュマシンでは高い性能が期待できない。 　　ブロック化など，アルゴリズム上の工夫が必要 並列粒度 BLAS2のため O(N2/p) であり，比較的高い。

共有メモリ向けの並列化 基本方針 PU間でのアクセス競合の回避 行列ベクトル積 および rank-2更新 の部分について，BLAS2レベルで並列化 並列BLASを使えば，並列化は極めて容易 PU間でのアクセス競合の回避 BLAS2はデータ再利用性がないため， 　 バスを通した共有メモリへのアクセスが 　 頻発 アクセス競合による性能低下を回避 　 するには，後に述べるキャッシュ向け 　 最適化技法が有効 キャッシュ PU0 PU1 PU2 PU3 バス メモリ

分散メモリ向けの並列化 (1-1) ブロックサイクリック列分割による並列化 （Dongarra et al., 1992） 各ステップの処理 行列を幅 b の列ブロックに分割し，ブロックを周期的にプロセッサに割り当てる。 各ステップの処理 C (k) 鏡像変換 ベクトル k (1) 鏡像変換ベクトル u(k) およびα(k) の作成 ・第 k 列を担当するPUが実行 ﾌﾞﾛｰﾄﾞｷｬｽﾄ (2) u(k) およびα(k) のブロードキャスト ・第 k 列を担当するPUが他の全PUに送信 ・二進木を使った場合，通信にかかる時間は O(N log p) 　（ベクトル u(k) の長さは平均 N/2） u(k)

分散メモリ向けの並列化 (1-2) 各ステップの処理（続き） (3) C (k) の下三角部分と u(k) との積 ・各PUは自分の行列要素を使って部分和ベクトルを計算 ・結果の集計／ブロードキャストの通信時間はO(N log p) (4) C (k) の上三角部分と u(k) との積 ・各PUは自分の行列要素を使って結果ﾍﾞｸﾄﾙの一部を計算 ・全対全通信により全PUに結果ﾍﾞｸﾄﾙの全要素を分配する 　通信時間はO(N log p) (3)と(4)の結果の集計，β(k) の計算，q(k) の計算 rank-2更新 C (k) = C (k) – u(k)q(k)t – q(k)u(k)t ・各PUで独立に計算 × ＝ 部分和の集計 ／ﾌﾞﾛｰﾄﾞｷｬｽﾄ × ＝ 全対全通信により 結果ﾍﾞｸﾄﾙを分配 以上より，ブロックサイクリック分割の通信時間は O(N log p)

分散メモリ向けの並列化 (2-1) Scattered Square 分割による並列化 （Choi et al., 1995） 行列を b×b のブロックに分割し，ブロックに対して大きさ √p× √ p のプロセッサテンプレートを周期的に割り当てる。 各ステップの処理 C (k) 鏡像変換 ベクトル k (1) 鏡像変換ベクトル u(k) およびα(k) の作成 ・第 k 列を担当するPU群が（通信しながら）実行 ﾌﾞﾛｰﾄﾞｷｬｽﾄ u(k) (2) u(k) およびα(k) のブロードキャスト ・第 k 列担当のPUが，まず横方向のPUにﾌﾞﾛｰﾄﾞｷｬｽﾄ ・次に，対角ブロックを担当するPUが縦方向のPUに 　ﾌﾞﾛｰﾄﾞｷｬｽﾄ ・通信にかかる時間は O(N log p /√p) ﾌﾞﾛｰﾄﾞｷｬｽﾄ

分散メモリ向けの並列化 (2-2) 各ステップの処理（続き） (3) C (k) の下三角部分と u(k) との積 ・各PUは自分の行列要素を使って部分和ベクトルを計算 ・横方向に集計し，対角PUに集める（O(N log p /√p) ）。 (4) C (k) の上三角部分と u(k) との積 ・縦方向に集計し，対角PUに集める（O(N log p /√p) ）。 (3)と(4)の結果の集計とブロードキャスト ・対角PUにおいて，(3)と(4)の結果を集計 ・結果の部分ベクトルを，横方向のPUにブロードキャスト（O(N log p /√p) ） ・結果の部分ベクトルを，縦方向のPUにブロードキャスト（O(N log p /√p) ） × ＝ 部分和を横方向に集計し， 結果を対角PUに集める。 × ＝ 部分和を縦方向に集計し， 結果を対角PUに集める。

分散メモリ向けの並列化 (2-3) 各ステップの処理（続き） 両並列化方式の比較 (6) β(k) の計算，q(k) の計算 各PUの計算時間 各ステップの計算時間は両方式ともO(N2 / p) 　（p に反比例して減少） 通信時間 ブロックサイクリック列分割は O(N log p) 　（pとともに増加） ブロックScattered Square 分割は O(N log p /√p) 　（pとともに減少） (6) β(k) の計算，q(k) の計算 (7) rank-2更新 C (k) = C (k) – u(k)q(k)t – q(k)u(k)t ・各PUで独立に計算 PU台数が増えるほど， ブロックSS分割が有利

超並列ではScattered Square分割が優位 サイクリック列分割とSS分割の性能比較 性能予測モデルによる推定性能 10 5 SS分割 10 4 ｻｲｸﾘｯｸ列分割 Mflops 10 3 10 2 10 10 1 10 2 10 3 10 4 Num. of Proc. 超並列ではScattered Square分割が優位

SMPクラスタ上での性能評価 16-nodes 4-nodes 1-node 1832[sec] 902[sec] 447[sec] ブロックScattered Square 分割による並列化（MATRIX/MPP　HZES1MDP使用） SR8000/F1（SMPクラスタ）の１～16ノードで評価。キャッシュ向け最適化あり。 N=2000 ～ 32000 のエルミート行列での性能 16-nodes 4-nodes 1-node 1832[sec] 902[sec] 447[sec]

４. 三重対角化処理のキャッシュ向け最適化 従来のハウスホルダー法の問題点 Bischof のアルゴリズム 性能・精度評価 ４.　三重対角化処理のキャッシュ向け最適化 従来のハウスホルダー法の問題点 Bischof のアルゴリズム 性能・精度評価 はじめに，研究の背景として，有限要素法とスパースソルバの必要性についてご説明します。

従来のハウスホルダー法の問題点 演算量 データの再利用性 行列サイズが N のときの演算量： 約 (4/3)N3 rank-2更新の演算量： 約 (2/3)N3 データの再利用性 演算はほとんどBLAS2のため，行列データの再利用性なし キャッシュマシンでは高い性能が期待できない。 　　ブロック化など，アルゴリズム上の工夫が必要

Bischof のアルゴリズム 原理（Bishof et al., 1993, 1994） 三重対角化を２段階で行うことの利点 密行列 A をまず半帯幅Lの帯行列 B に変換し，次にこの帯行列を三重対角行列 T に変換する。 三重対角化を２段階で行うことの利点 帯行列への変換は，BLAS3のみを使って実行可能。 帯行列から三重対角行列への変換の演算量はO(N2L)であり， 　　前半部に比べてずっと小さい。 次数 N 半帯幅 L 帯行列化 村田法 約 (4/3)N3 O(N2L) A B T

Bischof のアルゴリズム ブロック鏡像変換によるブロック列の消去 ブロック鏡像変換 H = I – UαUt Hは直交行列 ブロックベクトル ブロック鏡像変換によるブロック列の消去 ブロック鏡像変換 H = I – UαUt Hは直交行列 与えられたブロックベクトルの 第１ブロックを上三角行列にし， それ以外をゼロにする。 第kステップでの処理 左からH を乗算 右からH を乗算 左からH を乗算 非ゼロ要素 ゼロにしたい部分 影響を受ける部分

Bischof のアルゴリズム すべてBLAS3（行列乗算） [Step 1]　K = 1からN /L–1まで以下の[Step 2] ～ [Step 6]を繰り返す。 [Step 2]　 D(K) を第１ブロックが上三角行列で第２ブロック以下が ゼロ行列であるブロックベクトルに変換するブロック鏡像変換 I – U(K)α(K)U(K)tを求める。 [Step 3]　行列・ブロックベクトル積：　P(K) = C (K)U(K)α(K) [Step 4]　ブロックベクトルの内積：　β(K) = α(K)tU(K)tP(K) / 2 [Step 5]　ブロックベクトルの加算：　Q(K) = P(K) –U(K)β(K) [Step 6]　行列のrank-2L更新：　C (K) = C (K) – U(K)Q(K)t – Q(K)U(K)t すべてBLAS3（行列乗算）

Bischof のアルゴリズムの特徴と問題点 データの再利用性 演算のうち，O(N3) の部分はすべて BLAS3 化 キャッシュサイズに応じて L を選ぶことで，データ再利用性を最大化することが可能 計算精度 多様な例題を用いて体系的に評価した結果は，まだ報告されていない。 演算量 後半部（村田法）の演算量は L に比例して増加 全体としての実行時間が最小になるよう，L を適切に選ぶ必要がある。

４. 性能・精度評価 性能評価 テスト例題 N = 480，960，1920，3840 のフランク行列の三重対角化 ４.　性能・精度評価 性能評価 テスト例題 N = 480，960，1920，3840 のフランク行列の三重対角化 Bischof のアルゴリズムの性能を LAPACK と比較 評価マシン Opteron (1.6GHz) Alpha 21264A (750MHz) Power 5 (1.9GHz) 評価方法 演算量を(4/3)N3 とし，三重対角化に要した時間より各アルゴリズムのMFLOPS値を算出して比較 ブロックサイズ L は，全実行時間が最小になるよう決定 BLASルーチンとしては，全マシンで ATLAS 3.6.0 を使用

Opteron (1.6GHz) 上での性能 Wu の方法は LAPACK に比べて約２倍の性能を達成 Performance (GFLOPS) L’=32 Matrix size Wu の方法は LAPACK に比べて約２倍の性能を達成 N = 3840 のとき，Bischof の方法はピークの50%以上の性能を達成

Alpha 21264A (750MHz) 上での性能 Wu の方法は LAPACK に比べて約２倍の性能を達成 L=24, L’=4 L=48 Performance (GFLOPS) L’=32 Matrix size Wu の方法は LAPACK に比べて約２倍の性能を達成 N = 3840 のとき，Bischof の方法はピークの50%以上の性能を達成

Power5 (1.9GHz) 上での性能 Wu の方法は LAPACK に比べて２倍以上の性能を達成 ピークの 60% L=48, L’=2 L=96 Performance (GFLOPS) L’=64 Matrix size Wu の方法は LAPACK に比べて２倍以上の性能を達成 N = 11520 のとき，Bischof の方法はピークの60%の性能を達成

Power5 (1.9GHz) 上でのSMP性能 帯行列化部分のみの性能（村田法は含まず） L=48, L’=2 L=24, L’=4 Performance (GFLOPS) L=24, L’=2 L=24, L’=2 L=12, L’=2 Number of processors 帯行列化部分のみの性能（村田法は含まず） N = 7680 のとき， Bischof の方法は16PEで10倍の加速を達成

精度評価 テスト例題 評価マシン 評価方法 フランク行列 aij = min(i, j) ２次元ラプラス方程式の５点差分により得られる行列 Harwell-Boeing Sparse Matrix Collection の行列 乱数行列 （要素が [0, 1] の一様乱数である対称行列） 評価マシン Opteron (1.6GHz) 評価方法 各解法に対し，計算した固有値λi’の相対誤差の最大値 　　max 1≦i≦N |λi’ – λi| / |λi| を評価 比較対象の固有値λiとしては，解析解あるいはオリジナルのハウスホルダー法により求めた固有値を使用

Maximum relative error フランク行列に対する計算精度 Matrix size Maximum relative error Bischof 及び Wu の方法の誤差は，Dongarra の方法と同程度

Maximum relative error ２次元５点差分の行列に対する計算精度 Matrix size Maximum relative error Bischof / Wu の方法の誤差は Dongarra の方法に比べやや大きい。 しかし，増大の程度は１桁以内に収まっている。

Harwell-Boeing Collection の行列に対する計算精度 Matrix name Maximum relative error Bischof / Wu の方法の誤差は Dongarra の方法に比べやや大きい。 しかし，増大の程度は２～３倍以内に収まっている。

乱数行列に対する計算精度 Bischof / Wu の方法の誤差は Dongarra の方法に比べやや大きい。 しかし，増大の程度は１桁以内に収まっている。

キャッシュ向け最適化に関するまとめ 性能・精度評価の結論 キャッシュマシン向け三重対角化アルゴリズムとしては Wu のアルゴリズムが最高速であり，大規模問題では Dongarraのアルゴリズムの約２倍の性能が達成できる。 特に N = 3840の場合，Wu のアルゴリズムでは Opteron，Alpha 21264Aでピークの50%以上の性能が達成できる。 Bischof / Wuのアルゴリズムによる固有値の相対誤差は，Dongarraのアルゴリズムに比べてやや大きい。しかし，増大の程度は多くの場合２～３倍以内，最大でも１桁程度である。 はじめに，研究の背景として，有限要素法とスパースソルバの必要性についてご説明します。 より詳細なデータについては， 　　http://www.na.cse.nagoya-u.ac.jp/~yamamoto/work.html　を参照

５. 今後の展開 より多様な例題での精度評価 固有ベクトル計算部分の実装と性能・精度評価 共有／分散メモリ型計算機上での並列化と性能評価 ５.　今後の展開 より多様な例題での精度評価 固有ベクトル計算部分の実装と性能・精度評価 共有／分散メモリ型計算機上での並列化と性能評価 性能パラメータ L，L’ の自動最適化