HPCS5 5/5/9 5年ハイパフォーマンスコンピューティングと計算科学シンポジウム High Performance Comuting Symosium 5 などの行列とベクトルの演算 Level- 演算 は 演算回数 に対して必要となるデータ量が多く マルチコア計算機に おいて高い実行性能を実

HPCS5 5/5/9 5年ハイパフォーマンスコンピューティングと計算科学シンポジウム High Performance Comuting Symosium 5 帯行列の一般化固有値問題向け分割統治法 廣田 悠輔,,a) 今村 俊幸, 概要 本稿では 実対称正定値帯行列向けの一般化固有値解法を提案する 提案法は Elsner らによって 提案された三重対角行列の一般化固有値問題の分割統治法の拡張であり 三重対角行列向け解法の統治 フェーズを繰り返し適用することで一般の帯幅の帯行列の固有値問題を解く 近年のマルチコア CPU の 普及と性能向上により マルチコア計算機に適した数値解法の重要性はますます高くなっているが 問題 を標準固有値問題に変換して解く従来法はデータ再利用性の低い演算を多く含むため マルチコア計算機 上で高い性能を実現することが難しい 一方 提案法では演算の殆どが行列積として実行され 従来法に 比べて高い実行性能が実現できる Intel Xeon E5-66 ソケットを備えるマルチコア計算機における性 能評価では 次数 の五重対角行列の一般化固有値問題を解くとき 提案法は 従来法の. 倍高速 であり 9 GFLOPS ピーク性能比 77.6% の高い性能を示すことが確認された キーワード 一般化固有値問題 三重対角行列 帯行列 分割統治法 マルチコア Divide-and-Conquer Method for Banded Generalized Eigenvalue Problems Yusue Hirota,,a) Toshiyui Imamura, Abstract: In this aer, we resent a new solution method for symmetric-ositive definite generalized eigenvalue roblems of banded matrices. The roosed method is an extension of the divide and conquer method roosed by Elsner et al., which reeats the conquer hase of the divide and conquer method for a roblem of tridiagonal matrices. Recently, numerical solution methods are required to wor efficiently on modern multicore rocessors. However, the conventional methods show on such environment since they contain many cache inefficient oerations. On the hand, the roosed method is dominated by matrix roducts thus it shows higher erformance than the conventional methods. The roosed method is. times faster than the conventional method, achieving 9 GFLOPS (77.6% of the ea erformance) on a multicore environment (two octa-core Intel Xeon E5-66 CPUs). Keywords: generalized eigenvalue roblem, tridiagonal matrix, banded matrix, divide and conquer method, multicore. はじめに ピーク性能で動作する CPU に対してデータを供給し続け るだけのメモリ帯域をもたないことが一般的である した 近年 多くの計算機においてマルチコア CPU が利用さ がって マルチコア計算機において高い性能で演算を実行 れている マルチコア CPU を備える計算機 マルチコア するためには 度メモリから読み込んだデータを CPU の 計算機 では CPU は高いピーク演算性能をもつ一方 キャッシュメモリに蓄えて再利用し メモリからのデータ のロード回数をできるだけ削減する必要がある a) 理化学研究所 計算科学研究機構 RIKEN Advanced Institute for Comutational Science, Kobe, Jaan 科学技術推進機構 戦略的創造研究推進事業 Jaan Science and Technology Agency CREST yusue.hirota@rien.j 5 Information Processing Society of Jaan 行列やベクトルの計算のデータ再利用性について考え ると ベクトルの内積や加算などのベクトル同士の演算 Level- 演算 や 行列ベクトル積 行列のランク 更新 9

HPCS5 5/5/9 5年ハイパフォーマンスコンピューティングと計算科学シンポジウム High Performance Comuting Symosium 5 などの行列とベクトルの演算 Level- 演算 は 演算回数 に対して必要となるデータ量が多く マルチコア計算機に おいて高い実行性能を実現することが難しい 一方 行列 積などの行列同士の演算 演算 は 演算回数に 対して必要となるデータ量が少なく 適切にキャッシュメ モリを利用すれば 高い実行性能が実現できる したがっ て 数値計算アルゴリズムを基本行列演算の組み合わせと して構築する場合 できるだけ 演算が中心的とな るようにアルゴリズムを構築することが マルチコア計算 図 帯行列一般化固有値問題に対する解法アプローチ Fig. Solution aroaches for banded generalized eigenvalue roblems. 機で高い性能を実現するために必須となる 本稿では 半帯幅 が小さな値の実対称帯行列 A Rn n および 同じ半帯幅の実対称正定値帯行列 B R n n の一 題向けの分割統治法について述べる その中で 帯行列向 け分割統治法を提案する また 演算量 演算の種類につ いて分析し 従来法との比較を行う 第 節では 従来法 般化固有値問題 および提案法の精度および性能についてマルチコア計算機 Ax = λbx の固有値 λ 固有ベクトル x をすべて求める数値解法につ いて考える は n 組の固有値 固有ベクトル 固有 対 をもつ したがって の全固有対を求めることは を満たす対角行列 Λ R. 標準固有値問題を経由する解法 なく 両辺に左から S B-直交行列 X R n n を求 めることに等しく Λ の対角項 X の各列ベクトルがそれ ぞれ についてまとめる 一般化固有値問題 X (A λb)x = Λ λi n n 上で評価し 評価結果をもとに議論を行う 第 5 節で本稿 は A, B が帯行列か否かに関係 をかけることで (S AS )y = λy, y = S x の固有値 固有ベクトルとなる このような問 という標準固有値問題に変換することができる ただし 題に対する解法は 帯化前処理と組み合わせた密行列向け S は B = SS を満たす任意の行列である したがって 解法の部品 [] として応用可能であるほか 電子状態計算 一般化固有値問題 は コレスキー分解などにより に利用できる B SS と分解し C S AS を構成し C の標 問題 に対する解法は 図 に示されるように 様々 準固有値 固有ベクトルを求め 固有ベクトルを逆変換 なアプローチが考えられる 従来法では 赤や緑の線で示 することで解くことができる この原理に基づく解法は されるように 与えられた一般化固有値問題を標準固有値 数値計算ライブラリ LAPACK[6] に採用され DSBGV 問題に変換し 標準固有値問題を解いた結果を 一般化固 DSBGVD DSYGVD ルーチンとして実装されている 有値問題の固有ベクトルに逆変換するという手順が取ら 本節では を標準固有値問題を経由して解く つの れる しかしながら これらの解法は Level- Level- 演 解法について述べ その演算量および演算の種類について 算を多く含み マルチコア計算機で十分な性能を引き出す 分析する ことが難しい そこで 青の線で示される 中間形を経ず に直接 の一般化固有値 固有ベクトルを求めること. 行列の帯構造を利用する解法 を考える このような方法は = すなわち三重対角 本副節では 上述の原理に基づく解法のうち 行列 A, B 行列 の場合には Elsner らによって解法が提案されてお の帯構造を利用する解法について述べる このような帯構 り [] そのアルゴリズムは 演算が支配的となる 造を利用する解法は LAPACK では DSBGV DSBGVD 本研究では Elsner らの解法を の場合に拡張するこ などとして実装されている DSBGVD の解法は以下のよ とで 演算が支配的となる数値解法を新たに提案 うになる する なお Elsner らの方法とは別に [], [] で = の (i) 帯行列 B を B = SS と slit Cholesy 分解する 場合の解法が提案されているが 固有ベクトルの高精度化 手段 [5] の適用手段が確立されておらず また その解法 DPBSTF ルーチン (ii) A を C Z AZ と合同変換する ただし Z = S P の性質上 への拡張が困難であるなどの理由により であり P は C の帯幅が A に等しくなるように fill- 本稿ではこれらについては取り扱わない in を消去するような直交行列である また 同時に 本稿の構成は以下のとおりである 第 節では 標準固 Z S P を計算する DSBGST ルーチン 有値問題を経由する従来法について述べ その演算量 演 (iii) 帯行列 C を直交行列 Q によって T Q CQ と三 算の種類について分析する 第 節では 一般化固有値問 重対角行列に変換し 同時に X ZQ を計算する 5 Information Processing Society of Jaan

HPCS5 5/5/9 5年ハイパフォーマンスコンピューティングと計算科学シンポジウム High Performance Comuting Symosium 5 表 行列の帯構造を利用する解法の演算量内訳および演算の種類 表 密行列として扱う解法の演算量内訳および演算の種類 Table The number of FLOPs and the comutational attern Table The number of FLOPs and the comutational attern in the conventional method which exloits the band in the conventional method which does not exloit the structure of the matrix. band structure of the matrix. 演算量 種類 演算量 種類 ( + ) n Level- (i) /n A の変換 6n Level- および Level- (ii) n Z の計算 Level- および Level- (iii) 三重対角化 /n Level- が半分ずつ 6n Level- 分割統治法 /n Level- 逆変換 n (/)n n (i) (ii) (iii) 三重対角化 X の計算 (iv) 分割統治法 行列積 (/)n n n (iv) 変換する DTRSM ルーチン (i) から (iv) の演算量および支配的となる演算の種類を表 DSBTRD ルーチン (iv) T Q DQ と分割統治法により標準固有値問題を 解き DSTEDC ルーチン その後 行列積ルーチン DGEMM ルーチン により X X Q を計算する ことで の固有ベクトルに変換する ただし = すなわち A, B が三重対角行列 である場 合にはステップ (iii) はスキップされる (i) から (iv) の演 算量および支配的となる演算の種類を表 に示す 総演算 量は = の場合には (9/6)n + O(n ) の場合 には (/6)n + O(n ) となる そのうち 演算 は (/)n であり いずれの半帯幅 でも多くの Leve- Level- 演算が含まれる このため マルチコア計算機上 において高い性能 FLOPS 値 が得られない可能性があ る なお DSBGV は DSBGVD と同様に帯構造を利用 するが DSTEDC および DGEMM ルーチンの代わりに QR 法ルーチン DSTEQR を用いて (iv) のステップを実 行している したがって DSBGVD と同様に (ii) (iii) に おいて Level- Level- 演算が必要となり これらの部分 は同様の性能を示すと考えられる に示す 行列の帯構造を利用する解法と比べると演算量が 増大するが 標準固有値解法の三重対角化以外が 演算となるためマルチコア計算機などでより高い性能が得 られると予想され 結果的に帯構造を利用する解法より高 速に問題を解ける可能性がある. 一般化固有値問題向け分割統治法 本節では まず Elsner らによって提案された = の 帯行列 三重対角行列 向けの分割統治法アルゴリズムに ついて説明する 続いて その拡張である半帯幅が の帯行列向けの分割統治法アルゴリズムを提案し 提案法 および従来法のアルゴリズムの性質について高性能計算の 観点から議論する. 三重対角行列向け分割統治法 本副節では Elsner らの分割統治法について述べる.. 原理 三 重 対 角 行 列 A, B は 任 意 の 分 割 点 m を 定 め て bm,m+ = であれば A λb. 行列を密行列として扱う解法 問題 の A, B が帯行列であっても その疎構造を無 視して密行列として問題を解くことも可能である そのよ うな方法は LAPACK では DSYGVD などとして実装さ れており 以下の手順で の固有対を求める (i) 実対称正定値行列 B を密行列として B LL コレ スキー分解する DPOTRF ルーチン (ii) B の分解結果を A に作用させて C L AL を計 算する DSYGST ルーチン (iii) C QDQ と標準固有値問題を解く DSYGVD で 用いられる解法 DSYEVD ルーチン は C をブロッ ク化されたハウスホルダーによって三重対角化し 三 重対角行列を標準固有値問題向け分割統治法によって 固有値分解し 逆変換によって C の固有値分解に戻す という手順が取られる (iv) X = L Q を計算することで の固有ベクトルに 5 Information Processing Society of Jaan = (A A ρvv ) λ(b B vv ) () と ブ ロ ッ ク 対 角 行 列 と 共 通 の ベ ク ト ル v に よ っ て 表 現 さ れ る ラ ン ク 行 列 に 分 解 で き る た だ し A, B Rm m, A, B R(n m) (n m) であり ρ = am+,m /bm+,m v = bm+,m em sign(bm+,m ) bm+,m em+ である このような分解によって得られる B, B は正定 値行列であり A, A, B, B はそれぞれ対称三重対角行 列となる A, A, B, B が実対称かつ B, B は正定値行列である ので Yi (Ai λbi )Yi = Di λi (i =, ) () を満たす B -直交行列 Y B -直交行列 Y が存在する た

HPCS5 5/5/9 5年ハイパフォーマンスコンピューティングと計算科学シンポジウム High Performance Comuting Symosium 5 A λb = だし D, D は対角行列である したがって () の右 辺は Y = Y Y の合同変換によって (A A V EV ) λ(b B V V ) Y [(A A ρvv ) λ(b B vv )]Y = (D ρww ) λ(i ww ) を満たす 正整数 半帯幅 の実対称帯行列 A, B () Rm m A, B R(n m) (n m) 行列 V Rn 対角行 列 E R が存在する 具体的な, A, A, B, B, V, E と対角行列とランク 摂動の和に変換できる ただし D = D D w = Y v である の構成法については後述する 分解 (6) () の右辺を を満たす B, B はいずれも正定値行列である ことを示す ブロック対角行列 B B は B B = B + W [(D ρww ) λ(i ww )]W = D λi と対角化すれば () () (5) (5) より 正定値行列である B, B のいずれかが正定値行列でない と仮定する 行列 B が正定値でない場合 z B z を満 と定義すれば z (B B )z = z B z となり の解が X = Y W, Λ = D として表さ れる なお 一般化固有値問題 (5) V V T と正定値行列 B と半正定値行列 V V の和であるので たす z Rm が存在する このとき z := [z, ] Rn (Y W ) (A λb)(y W ) = D λi が成り立ち (6) は 必要に応じて減 次 デフレーション を行った後 一変数非線型方程式 secular 方程式 を解いて固有値を求め その後に対応す る固有ベクトルを計算することで解くことができる 具体 的な方法については [] を参照されたい.. アルゴリズム B B が正定値であることに矛盾する 行列 B が正定 値でない場合も同様である したがって B, B はいずれ も正定値行列である A, A, B, B が実対称かつ B, B は正定値行列である ので (Xi ) (Ai λbi )Xi = Di λi (i =, ) 以上の原理に基づく Elsner らの分割統治法は以下の手 順にまとめられる (i) 行列 A, B を () の形に分解する (ii) もとの行列よりも次数の小さな固有値問題 () を Elsner らの解法によって再帰的に解く (iii) 小さな固有値問題を解いた結果をもちいて w = Y v を計算する を満たす B -直交行列 X B -直交行列 X る ただし D, D (6) の右辺は X は対角行列である したがって := X X による合同変換で λ(b B V V )]X ei,i ui (ui ) ] = [D λ[i を求める ui (ui ) ] () (v) の固有ベクトル X = Y W を計算する と対角行列と 個のランク 行列の和に変換できる た の演算量は Y のブロック対角性を考慮する場合には だし U = (X ) V ul nm + n(n m) である したがって 常に m n/ D = D D である ここで として行列を中心付近で分割して再帰的に問題を解く場 合 総演算量は (/)n + O(n ) となる なお (iv) でデ フレーションが行われる場合 W を陽に計算せずに低次 が存在す (X ) [(A A V EV ) (iv) secular 方程式を解くことにより (5) を満たす W, D (i), (iii), (iv) の演算量はいずれも O(n ) であり (v) (7) は U の第 l 列ベクトル (W ) {[D e, u (u ) ] λ[i u (u ) ]}W = D λi (9) 元の密行列と特殊な構造をもつ疎行列の積として陰的に求 と [D e, u (u ) ] λ[i u (u ) ] を対角化 めることで (v) の行列積の演算量の削減が可能であるが する W により () の右辺を合同変換すると 本稿では W を陽に計算する場合について考えている (W ) {[D. 帯行列向け分割統治法 本副節では 一般の半帯幅 の帯行列に適用可能な分割 統治法を提案する 最初の副々節で提案法の原理について 述べ 次の副々節でアルゴリズムについて述べる = [D.. 原理 帯行列 A, B に対して m n を満たす分割点 m を定めたとき 5 Information Processing Society of Jaan ei,i ui (ui ) ] λ[i ui (ui ) ]}W ei,i ui (ui ) ] λ[i ui (ui ) ] i= i= と ランク 行列が つ少ない式に変換できる 同様の 変換

HPCS5 5/5/9 5年ハイパフォーマンスコンピューティングと計算科学シンポジウム High Performance Comuting Symosium 5 (W ) {[D ej,j uj λ[i uj ) ] ) ]}W = D λi, (uj (W ) {[D ei,i uj = X D X, D, X R (uj () が存在する ただし M (i : j, : l) は行列 M の第 i 行から (ui ) ] 第 j 行 第 列から第 l 列を取り出した (j i+) (l +) i=j λ[i ui 部分行列を意味する また D は対角行列 X は正則行列 ) ]}W (ui である このとき i=j = [D ei,i ui (ui ) ] i=j+ λ[i ui (ui ) ] V := i=j+ を j =,,..., と繰り返すことで 最終的に (X W W () W () O(m ) V S V S E := D,, V = X, V = X O(n m ) ) [(A A V EV ) λ(b B V V )](X W W () W () ) とおき A, A をそれぞれ A + V EV の B, B をそれ λi ぞれ B + V V の対角ブロックとおけば = を満たす =D という関係が得られ () (6) が成り立つことがわかる ただし S は任意の正則な の固有値 固有ベクトルは 実対角行列である Λ = D(), X = X W W () W () () 行列 S の決定方法の一つとして A, B のブロック対角 要素の修正量 と表されることがわかる.. つの帯行列の同時分割法 (6) を満たす, A, A, B, B V E の構成法につい f (S) := A( : m, : m) A F て述べる + A(m + : n, m + : n) A F 帯行列 A の分解 + B( : m, : m) B F A = A A VA VA, VA Rn () + B(m + : n, m + : n) B F = V S EV F + V () S EV () F は常に存在し 帯行列の標準固有値問題の分割統治法で用 いられており [7], [], [9] などで言及されている + V S V F + V () S V () F () B + VA VA = B B VB VB, VB Rn の自然な拡張として得られる こ をできるだけ小さくすることを考える 帯行列の標準固有 のとき A, A, B, B V = [VA, VB ] E = I O 値問題の分割統治法では 帯行列のブロック対角行列と摂 を満たす ただし I, O はそれぞれ 動行列への分解において ブロック対角要素への修正量が の単位行列 ゼロ行列を意味する したがって 上記の分 大きくなる場合に解の精度が悪化することが経験的に知ら 解を行うことで = である分解を構成できる れている 我々は予備実験で 帯行列の一般化固有値問題 を満たす分解は () は (6) 分解 () は一意ではなく また (6) 一意ではない また () を満たす分解も に示されるとおり固有ベクトル の分割統治法においても () が増大するほど解の残差ノ ルム X は + 個の行列の積として表現されるため 実際に固有 ベクトルを計算する際の演算量を減らすことを考えると ができるだけ小さな値となる分解が望まれる そこで 実 AX BXΛ F 上三角行列 := B(m+ : m+, m+ : m) の対角要 素がすべて非ゼロ すなわち が正則 であり の逆 が増大する傾向を確認しており 一般化固有値問題にお 行列と実上三角行列 := A(m + : m +, m + : m) いても修正量 () の積 の標準固有値が重複しない すなわち を防ぐために有効であると考えている そこで () が異なる対角要素をもつ という仮定の元では ヒューリスティック最小化する S の決定方法を考える は を小さく抑えることが解の精度悪化 () 非対称行列だが 固有値の重複しない実三角行列であるの V, V () の i 番目の列ベクトルをそれぞれ vi, vi で 分解 くと f (S) について 5 Information Processing Society of Jaan を とお

HPCS5 5/5/9 5年ハイパフォーマンスコンピューティングと計算科学シンポジウム High Performance Comuting Symosium 5 f (S) () () ( ei,i si,i vi (vi ) F + ei,i s i,i vi (vi ) F ) +( si,i vi (vi ) F + = s i,i vi (vi ) F ) () () ( ) () () ( ei,i + ) si,i vi (vi ) F +s i,i vi (vi ) F Algorithm Divide-and-conquer algorithm for banded generalized eigenvalue roblems : : : : 5: 6: = ( ei,i + ) ( () si,i vi F +s i,i vi F ) が成り立つ ここで 右辺が厳密に最小化されるように S uj (uj ) ]}W = D λi X X W ui (W ) ui (i = j +, j +,..., ) end for X := X (), Λ := D() A A A V EV, B B B V V Solve (Xi ) (Ai λbi )Xi = Di λi (i =, ) X := X X, D := D D, E := E ui (X ) vi (i =,,..., ) for j =,,..., do Solve (W ) {[D ej,j uj (uj ) ] λ[i 7: : 9: : を選ぶことで f (S) をヒューリスティックに最小化する 右辺の最小化は 最右辺の総和の各項を個別に最小化する (i =,..., ) を計算する ように si,i を決定すれば達成できる したがって () si,i = vi / vi (i =,,..., ) と S の対角要素を選べば 修正量 () (iv) 以下の手順を j =,,..., について繰り返すことに より 順番に W,..., W () を求め () はヒューリスティッ 繰り返し行い () の変換を に示される + 個の行列の積 の計算を進める ク最小化される が対角にゼロ要素を持つ場合や に重複固有値 ( a ) 一般化固有値問題 (9) () を [] に示され た反復法により解き D, W を求める が存在する場合でも ゼロ要素の個数や固有値の重複度に応 じたランクの分解が得られる 今 の対角要素のうち 第 ( b ) 行列積 X X W を計算する i 対角要素のみがゼロであるとする このとき B と同じ帯 ( c ) ui 構造をもつ B := B +α(em +i +em+i )(em +i +em+i ) を考えると := B (m + : m, m + + αei e i となり の第 i 対角要素に α を加えたも のになっている したがって α が非ゼロかつ が を満たす行列 V R (i = j +, j +,..., ) を 上記の手順により 最終的に Λ=D () X = X () の固有値 固有ベクトル が計算できる 以上をまとめたもの を alg. に示す なお = の場合 三重対角行列の場 味で提案法は Elsner らの分割統治法の拡張となっている A = A A V E (V ), B = B B V (V ) 次に アルゴリズムの演算量 演算の種類について考え, 対角行列 E R が存在す る したがって V = [V, α (em +i + sign(α)em+i )] n 合 に 提案法は Elsner らの分割統治法と一致し その意 重複固有値をもたない値にとられていれば (W ) ui 計算する : m + ) = る 行目は の計算 の固有値問題の求 解 A, A, B, B の計算から構成され いずれも E = E とおけば = + を満たす (6) が成り 型の演算によって求められ 演算量は O( ) となる 立つことがわかる が対角に 数のゼロ要素を持つ 行目は行列積として実行され 演算量はそれぞれ O(n ) 場合には 個だけ同様のランク 行列を加えることで となる 6 行目では 反復法による secular 方程式の求解が を満たす分解が可能である また の 支配的な演算となり [] で述べられる反復法が少ない反復 第 i 第 i 対角要素のみが重複する 重複固有値が存在す 回数で収束すれば 演算量は O(n ) となる 7 行目の行 る 場合 の第 i 対角要素がゼロである場合と同じ方 列積は j = のとき X のブロック対角性を考慮すれば 法で = + の分解を得ることができる 複数の重複が 演算量は nm + n(n m) で計算でき j のときは 存在する場合にはその重複度に応じた個数のランク 行列 密行列同士の積となり演算量は n となる = + 常に m n/ として行列を中心付近で分割して を加えることで 同様の分解が実現できる 行目の次数の低い固有値問題を提案法によって再帰的.. アルゴリズム.. で述べた原理に基づく 提案法の手順を以下に示す (i).. で述べた原理に基づき の計算を計算し の固有値問題を解いて E, V を求め (6) の A, A, B, B を計算する X, D を計算する (iii) 得られた X により () 右辺の 5 Information Processing Society of Jaan なる演算の種類をまとめたものを表 に示す このとき の総演算量は (/)( )n + O( n ) となる = となるように 行目の分解を行う場合には 総演算量は (ii) も と の 行 列 よ り 小 さ な 固 有 値 問 題 (7) を 解 き に問題を解く場合の各ステップの演算量および支配的と (/)( )n + O( n ) となる n である場合 総 演算量の第 項は無視できるため 演算の殆どすべてが 7 ui = (X ) vj 行目の行列積として実行されることになる

HPCS5 5/5/9 5年ハイパフォーマンスコンピューティングと計算科学シンポジウム High Performance Comuting Symosium 5 表 提案法の演算量内訳および演算の種類 表 Table The number of FLOPs and the comutational attern in the roosed method. Intel Xeon E5-66 コア. GFLOPS CPU 演算量 計算機環境 Table Comutational environment. Hyer-Threading 無効 ソケット 種類 の計算 O( ) メモリ の固有値問題 O( ) コンパイラ O( ) LAPACK Intel Math Kernel Library. 行目 O( ) 6 行目 O( ) 反復法 BLAS Intel Math Kernel Library. 7 行目 (/)( )n 行目 O( ) 行目 A, A, B, B の計算 6 GB Intel Fortran Comiler.. および LAPACK.5. 装で内部的に呼び出される LAPACK ルーチン DPBSTF や DSBGST など は Intel による実装 Intel を使 用している また 精度評価を行う際には 問題を帯行列. アルゴリズムの比較 表 にまとめられた つのアルゴリズムの演算量 の標準固有値問題に変換し標準固有値問題を QR 法によっ と演算の種類について比較する まず 提案法で = を て解く解法の Intel 実装 DSBGV も使用した 提 満たす分解が可能であると仮定すると 帯構造を用いる従来 案法は 行列積などの基本行列計算については BLAS ライ 法 密行列として扱う従来法 提案法の演算量はそれぞれ ブラリを用いて実装し alg. は Fortran および OenMP (/6)n +O(n ) = の場合のみ (9/6)n +O(n ) によって独自に実装した ただし alg. の 6 行目は W (/)n (/)( )n + O( n ) となる 帯構造を用 が陽に計算されるように実装した また 分割点は m と問 いる従来法の演算量は 最高次の項は 三重対角化ステッ 題を二等分するように選び alg. の 行目の小問題につ プが不要な = の場合を例外として に対して一定で いては 提案法を再帰的に適用して解いた ただし 行列 あり 低次の項のみが に伴って増大する また 密行列 の次数が 未満になった問題については LAPACK の として扱う従来法の演算量は の影響を受けない これに DSBGVD を適用して解いた 対して 提案法の演算量は最高次の項が に対して線形に テスト行列は A が半帯幅 の実対称行列 B が半帯幅 増大しており 従来法と比べて半帯幅 に対して強い依存 の実対称正定値行列を満たすようにするため 以下のよ 性がある n を仮定して演算量の最高次の項のみを比 うに乱数を用いて生成を行った { [, ) 乱数 ( i j ) ai,j =, (wise) (i = j) 較すると では提案法の演算量が最小になり では最大となる また それぞれ 演算として実行 される演算量は (/)n (/)n (/)( )n であ り 提案法のみ演算量の殆どすべてが 型演算とし て実行されることがわかる 半帯幅 では 提案法がもっとも演算量が少なく bi,j = [, ) 乱数 ( i j ). (wise) 性能面でも有利であるため 提案法の実行時間は つの従 このように生成される問題は 固有値分布がクラスタを持 来法よりも短くなると考えられる 一方 では 提 ちにくく 絶対値の極めて小さな固有値をもつ確率が低い 案法は 問題を密行列として問題を扱う従来法と比較し ため 分割統治法を適用する際に精度面で有利に働く可能 て Level- 演算の演算量が (/)n 少なく 演算 性がある しかしながら 任意の固有値分布をもつ帯構造 が (/)n n だけ多い したがって ある値以上の半 正定値性を備えたテスト行列 A, B を生成する方法が確立 帯幅の行列に対しては 従来法の実行時間がより短くなる されていないため 本実験では上記の方法で生成されたテ と予想される スト行列を使用する. 数値実験 マルチコア CPU 上で従来法および提案法の精度および 実験で使用する計算機環境は表 のとおりである. 精度評価 性能を評価する 標準固有値問題を経由する従来法の実装 本副節では Intel のみを使用した従来法の実装 には それぞれ Intel による LAPACK 実装 Intel [] DSBGV DSBGVD DSYGVD および提案法の各実装 に含まれる DSBGVD DSYGVD ルーチンを使用した について精度評価を行う n = = および = また 実行時間の内訳を評価するため 上記の実装とは としてテスト行列を作成し 求めた近似固有対の精度を 別に による LAPACK 実装の DSBGVD および 以下に定義される相対残差 近似固有ベクトルの B-直交 DSYGVD 各ルーチンのソースコードを修正して時間計測 性 QR 法を利用する従来法 DSBGV を基準にした解法 機能を付加した実装を作成した ただし 時間計測付き実 間の近似固有値の最大相対誤差 5 Information Processing Society of Jaan 5

HPCS5 5/5/9 5年ハイパフォーマンスコンピューティングと計算科学シンポジウム High Performance Comuting Symosium 5 Relative residual B-orthogonality Maximum error.e-9.e-.e-.e-.e-.e- TOTAL 5 5 5 5 5.E-5.E-6 DSBGV 図 DSBGVD DSYGVD Proosed method 固有値 固有ベクトルの精度 n =, = DPBSTF DSBGST DSTEDC DGEMM 6 図 DSBGVD の実行時間 n =, = Fig. The execution time of DSBGVD (n =, = ). Fig. The accuracy of the comuted eigenvalues and eigenvectors (n =, = ). TOTAL 6 DPOTRF DSYGST DSYEVD DTRSM Relative residual B-orthogonality Maximum error 5.E-.E-9.E-.E-.E-.E-.E-5 図 5.E-6 DSBGV DSYGVD Proosed method 固有値 固有ベクトルの精度 n =, = Fig. The accuracy of the comuted eigenvalues and eigenvectors (n =, = ). (DSBGV) λj λj max j (DSBGV) λj (DSBGV) によって評価する ここで λj (j =,,..., n) は DSBGV によって求められた近似固有値である 6 DSYGVD の実行時間 n =, = DGEMM solution of seqular equations AX BXΛ F X BX I F,, A F n Fig. 5 The execution time of DSYGVD (n =, = ). 図 DSBGVD 6 図 6 6 提案法の実行時間 n =, = Fig. 6 The execution time of the roosed method (n =, = ). く増大しており 実行時間の強い 依存性が確認できる 次に スレッド数の増大による加速についてみてみる 評価結果を図 に示す いずれの解法も同程度の相 と 帯構造を用いる従来法の実装は殆ど加速されないこと 対残差 B-直交性 固有値の最大相対誤差を示しており が確認できる また DSYGVD では スレッド数の増加 提案法は 三重対角行列 五重対角行列のいずれに対して によって性能は向上しているものの Level- 演算を含む も実用的な精度の解が得られていることが確認できる DSYEVD が性能上のボトルネックになり 6 スレッド実 行時の逐次実行時に対する加速率は = の場合に 9. 倍. 性能評価 となっている 一方 提案法は順調にスケールし 6 ス 半帯幅を =, 行列の次数 n = として 各実 レッド実行時の加速率は = の場合に. 倍となって 装を 6 スレッドで実行し 実行時間およびそ いる 従来法に対する加速率は 並列実行時に逐次実行時 の内訳を調べる = の場合の各実装の評価結果を図, よりも高く 提案法のマルチコア計算機における優位性を 5, 6 に = の場合の結果を図 7,, 9 に示す 確認できる 半帯幅の異なる つのテスト問題の実行結果を比較する また いずれのスレッド数で比較した場合でも 提案法 と DSBGVD 図 7 では = の場合に比べて = は =, の両問題で実行時間が従来法より短く での実行時間が増大しているが これは = の場合のみ では従来法より高速であるという副節. の予想に合致す スキップされる三重対角化ステップ DSBTRD ルーチン る結果となった の実行時間が = では加わったことが大きく影響してい る また DSYGVD 図 5 の実行時間は殆ど半帯幅 5. おわりに の影響を受けないことが確認できる 提案法 図 6 9 三重対角行列の一般化固有値問題向け分割統治法をもと は = のときの実行時間が = の場合に比べて大き に 帯行列の一般化固有値問題向け分割統治法を提案した 5 Information Processing Society of Jaan 6

HPCS5 5/5/9 5年ハイパフォーマンスコンピューティングと計算科学シンポジウム High Performance Comuting Symosium 5 TOTAL DPBSTF DSBGST DSBTRD DSTEDC DGEMM ンジンの開発 の援助を受けている スケールに対応した階層モデルによる超並列固有値解析エ 6 参考文献 [] 6 図 7 DSBGVD の実行時間 n =, = [] Fig. 7 The execution time of DSBGVD (n =, = ). TOTAL 6 DPOTRF DSYGST DSYEVD DTRSM [] 5 [] 図 6 DSYGVD の実行時間 n =, = Fig. The execution time of DSYGVD (n =, = ). DGEMM solution of seqular equations [6] 5 [5] 5 5 図 9 6 [7] 提案法の実行時間 n =, = Fig. 9 The execution time of the roosed method (n = [], = ). [9] 提案法の演算量は問題の半帯幅 に比例して増大するが では 帯行列の標準固有値問題を経由して解く従来 法と比べても演算量が少ない また 演算の殆どが行列積 として実行される 次数 の三重対角行列 五重対角 行列の一般化固有値問題をマルチコア計算機上で解いて性 能を評価し 提案法は従来法に対して三重対角行列では 6.6 [] Du, L. and Imaura, A.: Reducing Two Symmetric Matrices to Band Form by Congruence Transformations, 日 本応用数理学会 年度年会 予稿集. 66 67 (). Elsner, L., Fasse, A. and Langmann, E.: A divide-andconquer method for the tridiagonal generalized eigenvalue roblem, Journal of comutational and alied mathematics, Vol. 6, No.,. (997). Beattie, C., Ribbens, C. J., Dongarra, J., Kennedy, K., Mesina, P., Sorensen, D. and Voight, R.: Parallel solution of a generalized symmetric matrix eigenvalue roblem, Proceedings of the Fifth SIAM Conference on Parallel Processing for Scientific Comuting, Society for Industrial and Alied Mathematics,. 6 (99). Borges, C. F. and Gragg, W. B.: A arallel divide and conquer algorithm for the generalized real symmetric definite tridiagonal eigenroblem, Technical reort, DTIC Document (99). Gu, M. and Eisenstat, S. C.: A stable and efficient algorithm for the ran-one modification of the symmetric eigenroblem, SIAM Journal on Matrix Analysis and Alications, Vol. 5, No.,. 66 76 (99). Anderson, E., Bai, Z., Bischof, C., Blacford, S., Demmel, J., Dongarra, J., Du Croz, J., Greenbaum, A., Hammarling, S., McKenney, A. and Sorensen, D.: LAPACK Users Guide, Society for Industrial and Alied Mathematics, Philadelhia, PA, third edition (999). Arbenz, P.: Divide and conquer algorithms for the bandsymmetric eigenvalue roblem, Parallel comuting, Vol., No.,. 5 (99). Gansterer, W. N., Schneid, J. and Ueberhuber, C. W.: A Divide-and-Conquer Method for Symmetric Banded Eigenroblems-Part I: Theoretical Results (999). Pham, H. P., Imamura, T., Yamada, S. and Machida, M.: Novel aroach in a divide and conquer algorithm for eigenvalue roblems of real symmetric band matrices, Proc. Joint Int. Conf. Suecomuting in Nuclear Alications + Monte Carlo (SNA+MC),. 7 (). Intel Math Kernel Library (online): htts://software. intel.com/en-us/intel-ml (5..5). 倍 五重対角行列では約. 倍高速であり その優位性が 確認できた 今後の課題としては 任意の固有値分布のテスト行列生 成手法の確立後 固有値が密集する問題や 絶対値の小さ な固有値が存在する問題に対する提案法の精度評価があげ られる 謝辞 本論文に対して数々の貴重なコメントを頂いた匿 名の査読者に深い感謝の意を表す 本稿執筆にあたり多く の有用な意見を頂いた深谷猛氏 北海道大学 に深く感謝 する 本稿の図の作成の一部を支援して頂いた椋木大地氏 理化学研究所計算科学研究機構 にお礼申し上げる 本研究は 科学技術振興機構戦略的創造研究推進事業研 究領域 ポストペタスケール高性能計算に資するシステム ソフトウェア技術の創出 における研究課題 ポストペタ 5 Information Processing Society of Jaan 7