ペタフロップス環境における小規模行列用対称密行列固有値ソルバに向けて―逆変換の改良

全文

(1)情報処理学会論文誌. コンピューティングシステム. Vol. 3. No. 2. 1–8 (June 2010). する超並列計算機である．このことは，ペタフロップスの演算性能を達成する計算機環境（ペタフロップス計算機環境）では，10 万プロセッサもの超並列性を，1 つのアプリケーショ. ペタフロップス環境における小規模行列用対称密行列固有値ソルバに向けて——逆変換の改良片. 桐. 孝. 洋†1. ン内で達成しなくてはならないことを意味している．一方，従来開発されてきた並列計算機用の密行列数値計算ライブラリは，高効率並列実行のために超大規模行列での実行を前提とし，それに特化したアルゴリズムと実装方式が実現されている．ところが計算量が O(n3 ) で増加するため，大規模行列実行には並列化されているとはいえ時間的制約がある．加えて，スーパーコンピュータの運営上の制約で，1 ユーザが全系占有できる時間は数時間に限定される．. 超並列環境に向く固有値ソルバにおいて，従来法の問題は固有ベクトル逆変換部分の低い演算効率にあった．本論文では逆変換演算の特性を利用し，通信ブロック化，データアクセス局所化，および先行内積計算の手法を提案する．T2K オープンスパコンによる性能評価の結果，Weak Scaling で最大で 2.9 倍，Strong Scaling で最大で 2.3 倍の速度向上を達成した．. そこで本論文では，小規模行列に対してコア数を増加させても台数効果が見込める観点の超並列固有値ソルバの開発を目的にした．通信時間削減と単体性能向上の観点から従来方式の改良を行う．本論文の構成は以下のとおりである．2 章で小規模行列用の超並列アルゴリズムを用いた固有値ソルバの説明を行う．3 章は，逆変換方式の性能を改善する実装方式の提案を行う．. Towards Symmetric Dense Eigensolver for Small Sized Matrices on Petaflops Environment ——An Improvement of Inverse Transformation Takahiro. Katagiri†1. 4 章で，提案アルゴリズムを実際の並列計算機を用いて性能評価を行う．最後に，本論文で得られた知見を述べる．. 2. ABCLib DRSSED：小規模行列用超並列アルゴリズムを用いた固有値ソルバ 2.1 小規模行列用超並列ソルバに求められる特性. The conventional drawback of eigensolver on massively parallel processing was low efficiency of eigenvector computation for inverse transformation part. In this paper, communication blocking, data access localization, and previously computation for dot product methods are proposed with respect to characteristics of the inverse transformation computation. The result of the T2K open supercomputer indicated that we obtained the maximum speedup factors of 2.9× on weak scaling and 2.3× on strong scaling, respectively.. 2.1.1 小規模問題でも台数効果を維持行列の行/列サイズ np は，1 コアあたりの行列の行/列サイズを n1 で固定した場合の p コア使用時の行列の行/列サイズとする．行列は，行と列の 2 次元の方向を持つことから， √ p コア使用時に np = n1 × p である．密行列演算アルゴリズムは，行/列サイズ n に対し √ √ て O(n3 ) の計算量を必要するから，実行時間は n3p = (n1 × p)3 = n31 × p × p の関係から大まかに見積もることができる．. 1. はじめに. p コアを使うとして，並列化効率が 100%の場合でも np サイズの実行時には n3p /p = √ √ √ (n31 × p × p)/p = n31 × p より，1 コアあたりの実行時間が p 倍となる．一方，シミュ. 人類はペタフロップスコンピュータを持った．その並列数（コア数）は 10 万コアにも達. レーション全工程でソルバ呼び出しが 1 回の処理はまれである．通常は数十∼数百回，もしくは数千回呼び出される．この場合，運営上の都合から制約された実行時間内で，ソルバ 1. †1 東京大学情報基盤センター Information Technology Center, The University of Tokyo. 1. 回あたりの実行時間が規定される．たとえば，p を 1 万として，1 コアあたりの実行時間を. 4 時間に制限する．(n1 サイズの実行時間) × 100 = 4 時間であるから，並列化を考えない. c 2010 Information Processing Society of Japan .

(2) 2. ペタフロップス環境における小規模行列用対称密行列固有値ソルバに向けて——逆変換の改良. 場合の行列サイズを推定すると，(n1 サイズの実行時間) = 144 秒より，行列サイズとして. (Step1) 三重対角化：行列 A を相似変換で，三重対角行列 T に変換（QT A Q = T ）. は小規模になる．これが効率 100%で 1 万コアまで性能がスケールしないと実行時間内で処. (Step2) 三重対角行列 T の固有値行列 Λ を求解. 理が終わらない．つまり，小行列を超並列実行して効率が良いことが求められる．. (Step3) T の固有ベクトル行列 Y を求解. 2.1.2 データ領域が分散されていること. (Step4) 逆変換：Y を行列 A の固有ベクトル行列 X に変換（X = Q Y ）. 超並列実行において，必要データ領域が各コアに分散されていないと，全体問題サイズの大規模化にともないコアあたりのメモリサイズが増大する．小規模行列を対象としても，メモリ量もソルバ全体で O(n2 ) となることが知られている．なお，三重対角化に必要な行列. モリに関してスケーラブルでないと実行不可能となる．. 2.1.3 集団通信に関与するコア数の減少. Q は，各反復 k において，(n − k + 1) × (n − k + 1) の行列 Ak から計算される枢軸ベク. 超並列実行では通信時間の割合を減少させる必要がある．特に集団通信関数は並列実行に. トル uk から計算される Qk を用いて構成される．. 関わるコア数で実行時間が決まるが，全コアが通信に関連するようなアルゴリズムは好まし. Qk = (I − αk uk uTk ),. (3). ここで，αk は Ak から計算されるスカラ ∈ である．. くない．. 2.2 ABCLib DRSSED 概要. このとき. ABCLib DRSSED 2),4) は，実数対称密行列の固有値問題の任意の個数の固有値・固有ベクトルを計算できる機能を持つ並列数値計算ライブラリである．現在公開されている AB-. Q = Q1 Q2 · · · Qn−2 ,. (4). となる．. CLib DRSSED ver. 1.04 は，MPI-1 と Fortran90 で実装されている．数値アルゴリズム. 2.4 採用されている並列解法の特徴. として，対称密行列用 Householder 三重対角化，二分法，Gram-Schmidt 直交化付き逆反. 2.4.1 二次元サイクリック分散. 復法，Householder 逆変換の各ルーチンが利用できる．超並列実行に向くアルゴリズム1) が. (Step1)∼(Step4) を実現する解法として Householder 変換を用いた方法を用いる．この場合，負荷分散の観点から，行列 A を二次元サイクリック分散方式するのが有効である．. 採用されており，2.1 節の特性 2.1.1∼2.1.3 項を満たす．. 2.3 固有値解法の概略. A = (aij )（i, j = 1, 2, · · · , n）とするとき，. 対称密行列 A ∈ n×n ，実数 λ ∈ ，実数ベクトル x ∈ n とすると，以下の標準固有. aij → P (i % nprocx, j % nprocy). (5). への写像である．ここで，P (x, y) は，2 次元で表されたコア番号（x = 0, · · · , nprocx − 1，. 値問題. Ax = λx. (1). y = 0, · · · , nprocy − 1, nprocx ∗ nprocy = コア数）であり，%はモジュロ演算子である．. の解 λ を固有値，ベクトル x を固有ベクトルと呼ぶ．いま，式 (1) の固有値 n 個を対角要. (Step1) の三重対角化を上記の二次元サイクリック分散で行う場合，行列–ベクトル積と. 素に並べた行列 Λ を，Λ = diag(λ1 , λ2 , · · · , λn )，固有値 λi に対応する固有ベクトル xi を. 行列更新処理を並列化することで実現できる．ここで，対称行列を圧縮せずに行列要素すべ. 並べて構成された行列 X を，X = (x1 , x2 , · · · , xn ) とすると，式 (1) は. てを所有して二次元サイクリック分散を行うと，枢軸ベクトル計算のためのデータ転送量と √ 集団通信のための通信時間が 1/ p のオーダで削減できる超並列向きの三重対角化アルゴ. AX = XΛ,. (2). となる．いま，式 (2) の固有値行列 Λ，固有ベクトル行列 X を求めるが，その手順は以下である：. (Step1) および (Step4) は O(n3 ) である．(Step3) は解法と問題の性質に依存し，O(n)∼ 3. 3. リズムが構築できる1) ．一方，行列–ベクトル積と行列更新で必要なコアあたりの行列サイ √ ズが，二次元サイクリック分散の特性から O(n/ p) で減少していく．したがって，超並列実行でキャッシュミスヒットの削減が期待できる．. O(n ) である．したがって，ソルバ全体の演算量は O(n ) となる．メモリ量については，. 2.4.2 ブロック化なしアルゴリズムの採用. 入出力行列 A と X が密行列であるため O(n2 ) となるのは自明であるが，解法に必要なメ. 既存ライブラリの多くは，三重対角化 QT AQ の式レベルにおいて，ある幅 m の単位で. 情報処理学会論文誌. コンピューティングシステム. Vol. 3. No. 2. 1–8 (June 2010). c 2010 Information Processing Society of Japan .

(3) 3. ペタフロップス環境における小規模行列用対称密行列固有値ソルバに向けて——逆変換の改良. 1 度に変換するブロックアルゴリズムが採用されている．さらにデータ分散である二次元サ. (b) は，1 演算あたりのメモリバンド幅が低下している近年のマルチコア計算機で深刻な. イクリック分散においても，正方行列 m × m の単位で行う．この理由は，並列実装の容易さから m × m の単位で BLAS 演算を行う実装になっているからである．一方，コアあた. 問題となる．三重対角化では，式レベルのブロック化をしなくても二次元サイクリック分散 √ の特性から，行列サイズが O(n/ p) で減少していく．このため超並列実行時にキャッシュ. りの BLAS 演算の単体性能を高めるため，m は小さくできない．このことが，超並列実行. サイズまでアクセス領域を縮小できる．一方，逆変換では列ブロック分散のため，並列実行. で激しい負荷バランスの劣化を引き起こす．そこで文献 1) では，データ分散幅を 1 として. でもアクセス領域は減少せず O(n) となる．ゆえにキャッシュミスヒットが減少せず，単体. 超並列実行における負荷バランスの劣化を回避する．ただし，QT AQ の式レベルのブロッ. 実行性能がメモリバンド幅性能に律速されるという問題を生じる．. 3.2 通信ブロック化. ク化ができないため，コアあたりの単体性能を高めることができない．. 2.4.3 通信時間削減のための重複所有. 枢軸ベクトル uk を収集後，枢軸ベクトル 1 本を用いて行列を更新するシンプルな実装. 文献 1) の特徴として，三重対角化時に生成される枢軸ベクトル uk を各コアで重複し √ て所有する．ここで，枢軸ベクトル uk は p 個のコア間でデータ分散はサイクリック分 √ 散となる．枢軸ベクトル保存に必要なメモリサイズは O(n × n/ p) であり，完全分散の. は演算の最適化の余地がない．そこで，あるブロック幅 m 本の枢軸ベクトル uk を収集し，その後 m 本の枢軸ベクトルを 1 度に利用し逆変換の行列更新をする．演算時間を増加させ，通信時間の割合を減少できるため演算性能の向上が期待できる．この実装は，図 1 になる．. 場合の O(n × n/p) に対しメモリ量が増大するが，分散はされる．一方で (Step4) の逆変. ここで，各コアが担当する固有ベクトルの範囲は nstart∼nend で，ブロック幅 m の大域. 換時に，枢軸ベクトル収集時間が完全分散時では O(n × n log p) となるが，文献 1) では √ √ O((n × (n/ p) log p) に削減できる．. 的なインデックス範囲は jend∼j である．. 2.4.4 超並列向き三重対角行列固有値ソルバ. kk global = 0. (Step2)，(Step3) は，2 分法–逆反復法など古典的解法を採用する場合，並列実行時にボト. do jbk=j, jend, -1. ルネックとなる．本ソルバでは，超並列実行に向く先進解法 MRRR 法（Multiple Relatively. Robust Representations）3) を採用した．. local i = i - nstart + 1. 3. 逆変換の改良. norm = 0.0d0 do k=jbk, n. 3.1 単体性能の劣化問題. norm = norm + bufM(kk global+k-jbk+1) * X(k, local i). 逆変換は，式 (3) と (Step4) から以下となる．. enddo. X = (I − α1 u1 uT1 )(I − α2 u2 uT2 ) · · · (I − αn−2 un−2 uTn−2 )Y. (6). いま，X は固有ベクトル行列で，ユーザの利用の都合からベクトル要素を分散しない．列ブロック分散されて戻る．この要求を満たすため，三重対角行列 T の固有ベクトル行列 Y も列ブロック分散となる．先述のとおり文献 1) のアルゴリズムでは枢軸ベクトル uk は重複してサイクリック分散で保持されているが，uk を収集さえすれば自明な並列性が式 (6) に存在する．一方問題となるのは，(a) サイクリック分散されている uk の収集時間；(b) 単体演算性能の劣化である．. (a) は，uk を収集して，その後 (I − αk uk uTk )Yk の演算を交互に行うシンプルな実装では，通信時間が隠ぺいできない場合は演算効率劣化の要因となる．. 情報処理学会論文誌. do i=nstart, nend. コンピューティングシステム. Vol. 3. No. 2. 1–8 (June 2010). norm = ALP(jbk) * norm do k=jbk, n X(k, local i) = X(k, local i) - norm * bufM(kk global+k-jbk+1) enddo enddo kk global = kk global + (n-jbk+1) + 1 enddo 図 1 通信ブロック化法（CB）のカーネル Fig. 1 The kernel of communication blocking (CB) method.. c 2010 Information Processing Society of Japan .

(4) 4. ペタフロップス環境における小規模行列用対称密行列固有値ソルバに向けて——逆変換の改良. do i=nstart, nend. まる場合，キャッシュ上のデータを再利用し更新ができる．したがって，キャッシュに収ま. kk global = 0. る jbk∼n の範囲で劇的な速度向上が期待できる．一方，jbk∼n がキャッシュ範囲外である場合，最適化の効果は少ない．. local i = i - nstart + 1. 3.4 次ステップの内積演算を先行して行う方式. do jbk=j, jend, -1. 図 2 では，配列 X の連続アクセスについて次元数が大きい場合にキャッシュミスが避け. norm = 0.0d0. られないという問題がある．このキャッシュミスを避ける手法として，次のステップでの内. do k=jbk, n. 積計算必要な配列 X の要素が決まり次第，その要素を用いて先行して内積計算を行う方式. norm = norm + bufM(kk global+k-jbk+1)*X(k, local i). が考えられる．この値の決定範囲を，キャッシュサイズ内の長さ単位で行うことでデータの. enddo enddo. 局所化を行う．このことで，次元数 n に無関係に演算効率を高めることができる．この方. norm = ALP(jbk) * norm. 式を図 3 に示す．図 3 の先行内積計算法（PD）の kblksize は，先行して演算を行う長さの変数であり，. do k=jbk, n X(k, local i) = X(k, local i) - norm * bufM(kk global+k-jbk+1). キャッシュサイズに依存する性能パラメータである．. 3.5 LAPACK 方式との違い. enddo. LAPACK で採用されている逆変換方式は，前提として，ブロック化アルゴリズムを採用. kk global = kk global + (n-jbk+1) + 1. した三重対角化ルーチン中で蓄積された枢軸ベクトル行列を用いるものであり，本方式とは. enddo 図 2 データアクセス局所化法（AL）カーネル Fig. 2 The kernel of data access localization (AL) method.. 前提が異なる．アルゴリズムの違いにより，演算のデータ依存関係も異なることから直接比較ができない．本方式はすでに説明したように，負荷バランスの劣化からブロック幅が大きくとれないほど超並列実行する場合，および，小行列を処理する場合において，負荷バラン. ALP は，式 (3) のスカラ αk を保持している配列である．集積された枢軸ベクトル uk は配列 bufM に 1 次元配列として圧縮して収納されている．この圧縮方式とは，枢軸ベクトル収納用の配列 Q を n × n の 2 次元配列として持つとき，枢軸ベクトル uk の収納場所 Q(1,k). · · · Q(k-1,k) の値は 0 となるが，無駄な領域でかつ連続アクセス妨げとなり性能劣化要因となるので，uk−1 のデータの直後に uk のゼロ以外の値を詰め込むという方式である．図 1 の通信ブロック化法（CB）では，最内ループにおける配列 X，および枢軸ベクトルを集積した配列 bufM の連続アクセスが n に依存し，かつ，ノード内所有ベクトル数のループ i に依存するため，データアクセスの局所化が達成できない．. ス改善のメリットからブロック化を行わない三重対角化に基づく逆変換手法である．以上の前提を考慮したうえで，LAPACK の逆変換ルーチンである dormqr.f に対する外見上の違いについて表 1 にまとめる．. 4. 性能評価 4.1 計算機環境東京大学情報基盤センターに設置されている T2K オープンスパコン（東大版）を利用した．各ノードは，AMD Opteron 8356（2.3 GHz，4 コア）を 4 台（4 ソケット）搭載. 3.3 通信ブロック化にともなうデータアクセスの局所化. しており，ノードあたりのメモリ量は 32 GB である．理論最大演算性能は，ノードあたり. 通信ブロック化法（CB）により m 本の枢軸ベクトルを収集しているので，このデータ. 147.2 GFLOPS である．キャッシュサイズは，L1 命令キャッシュ，L1 データキャッシュと. を再利用し配列 X のアクセスを局所化することができる．これは，図 1 において jbk ルー. もに 64 Kbytes，L2 キャッシュは 512 Kbytes である．L1 キャッシュと L2 キャッシュはコ. プと i ループを交換することにほかならない．この方式を図 2 に示す．. アごとに独立している．L3 キャッシュ2048 Kbytes であり，ソケット内で共有されている．. 図 2 のデータアクセス局所化法（AL）では，配列 X のアクセスデータがキャッシュに収. 情報処理学会論文誌. コンピューティングシステム. Vol. 3. No. 2. 1–8 (June 2010). 各ソケットには，ローカルメモリ 8 GB（ソケットあたり 8 GB × 4 ソケット = 32 GB）が. c 2010 Information Processing Society of Japan .

(5) 5. ペタフロップス環境における小規模行列用対称密行列固有値ソルバに向けて——逆変換の改良表 1 LAPACK 方式と提案方式の違い Table 1 The differences between LAPACK method and the proposed method.. do i=nstart, nend kk global=0; local i=i-nstart+1; jbk=j;. 比較項目. NORMS1 = 0.0d0. 必要メモリ領域（枢軸ベクトル行列保存用以外）. do k=jbk, n. 固有ベクトル配列 X のコピー. dtemp=bufM(kk global+k-jbk+1); NORMS1=NORMS1+dtemp*X(k, local i);. 枢軸ベクトル行列 Q の圧縮収納. enddo. 先行内積計算主演算. do jbk=j, jend+1, -1. LAPACK 方式 n × m 領域 2 つ m × m 領域 1 つ n × m 領域を 1 回適用不可適用不可 BLAS2 および BLAS2-k. 提案方式なしなし適用適用 BLAS2 および BLAS1. dtemp2 = ALP(jbk) NORMS1=dtemp2*NORMS1; kk global2=kk global+(n-jbk+1)+1; k=jbk-1;. ある．通信性能は運用クラスタ群で異なる．ここでは Miri-10G が 4 本実装されており，最. dtemp=bufM(kk global2+k-jbk+2); NORMS2=dtemp*X(k, local i);. 大で 5 GB/sec の双方向性能を有する A 群を利用している．コンパイラは，日立最適化 Fortran90 V01-00-/B で，コンパイラオプションは -Oss. do kk=jbk, n, kblksize kkend = min(n, kk+kblksize-1). -noparallel である．ピュア MPI 実行を行っている．プロセスとプロセスが利用するメモ. do k=kk, kkend. リ領域についてハードウェア構成から高速となるように割り当てるため，numactl コマンドを MPI 起動時に実行している．. dtemp=bufM(kk global+k-jbk+1). 4.2 固有値ソルバ. X(k, local i)=X(k, local i)-NORMS(1)*dtemp. ABCLib DRSSED ver.1.04 4) に，LAPACK 3.1.1 における MRRR 法ルーチン（dstegr）. enddo. 相当のコードを実装し新規開発した．dstegr は分散並列版ではないため，各コアに均等な. do k=kk, kkend. 担当範囲を決め，その範囲について独立に固有ベクトル計算する実装となっている．ここ. dtemp=bufM(kk global2+k-jbk+2). で均等な担当範囲とは，コア数を p とするとき，n/p ずつの範囲で分担することであるが，. NORMS2=NORMS2+dtemp*X(k, local i). n/p が割り切れないとき，MPI ランクが小さいコアから順に n/p の切り捨て数に 1 を加算. enddo enddo. した数を担当する．自動チューニング機能を適用しておらず，デフォルトパラメータによる. NORMS1 = NORMS2. 実行である．アンローリングを適用することにより，さらなる速度向上が期待できる．. 4.3 テスト行列と演算精度. kk global = kk global + (n-jbk+1) + 1. MRRR 法を採用する場合，並列固有値ソルバのほとんどの時間は三重対角化と逆変換の. enddo jbk = jend; dtemp2 = ALP(jbk). 時間となる6) ．Frank 行列以外の行列では，三重対角行列の全固有値・全固有ベクトル計算. NORMS1 = dtemp2 * NORMS1. の実行時間が増加することがあるが，我々の経験では増加率はたかだか 10 倍程度である7) ．. do k=jbk, n. したがって，入力行列が変わってもほとんどの時間は三重対角化と逆変換の時間になるとい. dtemp=bufM(kk global+k-jbk+1); X(k, local i)=X(k, local i)-NORMS1*dtemp;. える．三重対角化と逆変換の時間は，入力行列の特性に依存しない．そこで，テスト行列には Frank 行列を利用した．Frank 行列の全固有値および全固有ベクトルの計算性能を評価. enddo. する．. enddo Fig. 3. 情報処理学会論文誌. 図 3 先行内積計算法（PD）のカーネル The kernel of precedence dot-products (PD) method.. コンピューティングシステム. Vol. 3. No. 2. 1–8 (June 2010). 式 (2) の標準固有値問題の演算精度について，各固有値における解析解との誤差最大値，. c 2010 Information Processing Society of Japan .

(6) 6. ペタフロップス環境における小規模行列用対称密行列固有値ソルバに向けて——逆変換の改良. 固有ベクトル直交性，元の行列 A に対する固有方程式の各残差ベクトル最大値で精度を確. 4.5 Strong Scaling 性能. 認した．. 図 5 に，小規模行列サイズとして N = 10,000 に固定し，1 ノード（16 コア）∼64 ノー. 4.4 Weak Scaling 性能. ド（1,024 コア）まで変化させる Strong Scaling の性能を載せる．. 図 4 に，行列 A について 1 ノードのメモリ量からほぼ限界となるノードあたり 9,600 ×. 図 5 から，シンプルな実装に対する通信ブロック化法，およびデータアクセス局所化法に. 9,600 に固定しコア数を増加させる Weak Scaling における逆変換の性能を示す．なお，逆. 対する先行内積計算法の効果がないことが分かる．データアクセス局所化法について，そ. 変換の FLOPS 値は (4/3)n3 + (5/2)n2 + (1/2)n とした．. の効果は 0∼2.3 倍の速度向上である．ノード数が増加すると BLK=16 に対する通信量が演. 図 4 では，シンプルな実装 ORG（BLK = 1），ここで BLK は通信ブロック化の個数，に対し，通信ブロック化法（BLK = 16）をすることで 1.18∼1.6 倍の速度向上を得る．. 算量に対し増大し性能が劣化する．しかし，1 コアに割り当てられる固有ベクトル数を勘案し，適切な通信ブロック幅 BLK を選ぶことで性能改善が期待できる．. さらにデータアクセス局所化法を導入すると，ORG（BLK = 1）に対し 1.64∼2.7 倍の. また，1,024 コア実行ですべての方式がほぼ同一の性能となっているが，これは並列化に. 速度向上を得る．先行内積計算法を導入することで，ORG（BLK = 1）に対し 1.76∼2.9. より演算時間が少なくなり，かつ通信時間がほとんどを占めるようになるため，通信ブロッ. 倍に達する．先行内積計算法の更新幅 kblksize は，1 ノード（16 コア），2000 次元の行. ク化をしても同期的な通信方式では通信時間の隠ぺいができなくなるからと考えられる．こ. 列で kblksize=64，128，256，512，および 1024 での実行結果から，kblksize=1024 で. の状況では，枢軸ベクトルの収集方式を非同期化するなどの通信時間削減方式が必要と考え. 固定した．. られる．. 図 4 逆変換の Weak Scaling 性能 Fig. 4 The weak scaling performance of the inverse transformation.. 情報処理学会論文誌. コンピューティングシステム. Vol. 3. No. 2. 1–8 (June 2010). 図 5 逆変換の Strong Scaling 性能 Fig. 5 The strong scaling performance of the inverse transformation.. c 2010 Information Processing Society of Japan .

(7) 7. ペタフロップス環境における小規模行列用対称密行列固有値ソルバに向けて——逆変換の改良表2. 小規模行列 N = 10,000 におけるソルバ内の各処理の実行時間の詳細．逆変換には，先行内積計算法（BLK = 8）を実装している Table 2 The break-down of execution time for each processes on the solver for small sized problem N = 10,000. PD (BLK = 8) is implemented on the inverse transformation. ノード数（コア数）. 三重対角化［秒］. 行列 T の再分散［秒］. 1 ノード（16 コア）. 372. 2 ノード（32 コア）. 逆変換［秒］. 全計算時間［秒］. 0.002. MRRR 法による全固有値・全固有ベクトル計算［秒］ 5.27. 131. 510. 190. 0.003. 2.86. 66.3. 260. 4 ノード（64 コア）. 97.0. 0.003. 1.76. 34.4. 133. 8 ノード（128 コア）. 53.0. 0.004. 0.838. 18.4. 72.7. 16 ノード（256 コア）. 27.1. 0.006. 0.473. 7.92. 35.8. 32 ノード（512 コア）. 20.1. 0.007. 0.286. 5.34. 26.0. 64 ノード（1,024 コア）. 10.1. 0.008. 0.199. 7.37. 18.0. 図 6 三重対角化の Strong Scaling 性能 Fig. 6 The strong scaling performance of the tri-diagonalization.. 化の時間に比べ 1/2.8 程度まで削減されており，もはやボトルネックではない．本改良を行う前，1 ノード（16 コア）でのシンプルな実装の逆変換の実行時間は，三重対角化の時間図 6 に，三重対角化における Strong Scaling の性能を載せる．. に対し 1/1.2 倍程度かかっていたことを考慮すると，劇的な性能改善を行うことができた．. 図 6 から，三重対角化においては 16 コア∼1,024 コアにおいてスケーラブルな性能が得. 一方，MRRR 法での三重対角行列の全固有値・全固有ベクトル計算部分は良好な台数効果. られていることが分かる．16 コアでの対ピーク実行性能は 5.4/147.2 ∗ 100 = 3.66%と低い. がある．さらに，全体の時間に占める割合は 1%程度に満たないので，無視できる時間とい. が，自動チューニング機能によるアンローリングの適用で，キャッシュ範囲外のデータサイ. える．. 5). ズで 5%程度，キャッシュ範囲内のデータサイズで 18%程度まで効率を改善できる．. 4.6 ソルバ全体処理の詳細. 演算精度に関し，理論固有値に対する最大誤差は 0.158E-008，フロベニウスノルムによる最悪の直交性は 0.110E-009，最大の残差ベクトルの 2 ノルムは 0.146E-006 であった．. 表 2 に，小規模行列 N = 10,000 におけるソルバ内の各処理の実行時間の詳細を載せる．ここで逆変換には，内積計算法（BLK = 8）を採用した．BLK = 8 にした理由は，64 ノード（1,024 コア）実行時には，1 コアあたりが保持する固有ベクトル数が最大で 10 個しか. 5. 関連研究式レベルのブロック化を行い演算効率を高めるが，ブロック対角行列の固有値・固有ベクトルを直接求めることで並列性を増加させる新しい固有値ソルバの研究が今村8) により開. ないため，BLK = 16 が実行できないという理由による．表 2 から，全固有値・全固有ベクトルを計算する場合において，最も時間を要する部分. 始されている．. は三重対角化で多くて約 70%を占めることが分かる．一方，逆変換の占める割合は多くて. また，本論文で提案する先行内積計算法は，三重対角化における「Wilkinson の技巧」と. 約 40%程度である．また，逆変換の時間は，1 ノード（16 コア）実行時において三重対角. 考え方が類似している．この技法は近年の CPU 上で，村上9) により有効性が検証されて. 情報処理学会論文誌. コンピューティングシステム. Vol. 3. No. 2. 1–8 (June 2010). c 2010 Information Processing Society of Japan .

(8) 8. ペタフロップス環境における小規模行列用対称密行列固有値ソルバに向けて——逆変換の改良. いる．. 6. おわりに本論文では，ブロック化しない超並列向き固有値ソルバでの逆変換処理において，従来法のシンプルな実装で問題となっていた性能劣化問題を改善するため，通信ブロック化法，データアクセス局所化法，先行内積計算法の提案を行った．性能評価の結果，Weak Scaling において最大で 2.9 倍，Strong Scaling において最大で 2.3 倍の速度向上を得た．今後の課題として，さらなる速度向上を得るため，通信隠ぺい手法の開発が必要である．また，適切な通信ブロック幅，キャッシュサイズに依存する更新サイズ幅，およびアンローリング段数など，単体性能を向上させるチューニング技法が必要である．提案手法に対する. 5) 片桐孝洋：T2K オープンスパコン（東大）における自動チューニング機能付き固有値ソルバの性能評価，東京大学情報基盤センタースーパーコンピューティングニュース， Vol.10, No.5, pp.52–65 (2008). 6) 片桐孝洋：ペタスケール環境に向けた固有値ソルバの開発，東京大学情報基盤センタースーパーコンピューティングニュース，Vol.11, No.1, pp.47–59 (2009). 7) Katagiri, T., Voemel, C. and Demmel, J.W.: Multi-section with Multiple Eigenvalues Method for Computing Eigenvalues in Symmetric Tridiagonal Eigensolvers, 情報処理学会研究報告，Vol.2006-HPC-105, pp.25–30 (2006). 8) 今村俊幸：10 万超コアを駆使する固有値ソルバについての検討，情報処理学会研究報告，Vol.2009-HPC-121 (2009). 9) 村上弘：両側ハウスホルダ変換に対する Wilkinson の著書 AEP 中の「技巧」について，情報処理学会研究報告，Vol.2008, No.125, pp.67–72 (2008). (平成 21 年 9 月 25 日受付). 自動チューニング方式の開発と性能評価が必要であろう．謝辞本研究は，科学技術研究費補助金，基盤研究（B），課題番号：21300007「メニー. (平成 22 年 3 月 22 日採録). コア・超並列時代に向けた自動チューニング記述言語の方式開発」，および，科学技術研究費補助金，特定領域研究（情報爆発），課題番号：21013014，「情報爆発時代のロバストな. 片桐孝洋（正会員）. 自動チューニングシステムに向けた数理的基盤技術の研究」の支援による．また，有益なコ. 東京大学情報基盤センター特任准教授．1994 年豊田工業高等専門学校情報工学科卒業．1996 年京都大学工学部情報工学科卒業．2001 年東京. メントをいただいた査読者の各位に感謝する．. 参. 考. 文. 大学大学院理学系研究科情報科学専攻博士課程修了．博士（理学）．2001. 献. 年 4 月日本学術振興会特別研究員 PD，12 月科学技術振興機構研究者，. 1) Katagiri, T. and Kanada, Y.: An Efficient Implementation of Parallel Eigenvalue Computation for Massively Parallel Processing，Parallel Computing, Vol.27, No.14, pp.1831–1845 (2001). 2) Katagiri, T., Kise, K., Honda, H. and Yuba, T.: ABCLib DRSSED: A Parallel Eigensolver with an Auto-tuning Facility，Parallel Computing, Vol.32, Issue 3, pp.231–250 (2006). 3) Parlett, B.N. and Dhillon, I.S.: Relatively Robust Representations of Symmetric Tridiagonals, Linear Algebra and Appl., Vol.309, No.1-3, pp.121–151 (2000). 4) ABClib DRSSED. http://www.abc-lib.org/. 情報処理学会論文誌. コンピューティングシステム. Vol. 3. No. 2. 1–8 (June 2010). 2002 年 6 月電気通信大学大学院情報システム学研究科助手，2005 年 3 月∼2006 年 1 月米国カリフォルニア大学バークレー校コンピュータサイエンス学科訪問学者を経て，2007 年 4 月より現職．超並列数値計算アルゴリズム，およびソフトウエア自動チューニングの研究に従事．2002 年情報処理学会山下記念研究賞受賞．2007 年 Microsoft. INNOVATION AWARD 2007 アカデミック部門最優秀賞受賞．日本ソフトウェア科学会，日本応用数理学会，ACM，IEEE-CS，SIAM 等各会員．. c 2010 Information Processing Society of Japan .

(9)