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

Size: px

Start display at page:

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

しょうこふしはら
5 years ago
Views:

1 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 deﬁnite 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 eﬃciently on modern multicore rocessors. However, the conventional methods show on such environment since they contain many cache ineﬃcient 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

2 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 に等しくなるように ﬁll- 本稿ではこれらについては取り扱わない in を消去するような直交行列であるまた同時に本稿の構成は以下のとおりである第節では標準固 Z S P を計算する DSBGST ルーチン有値問題を経由する従来法について述べその演算量演 (iii) 帯行列 C を直交行列 Q によって T Q CQ と三算の種類について分析する第節では一般化固有値問重対角行列に変換し同時に X ZQ を計算する 5 Information Processing Society of Jaan

3 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 が存在するた

4 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= とランク行列がつ少ない式に変換できる同様の変換

5 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 をとお

6 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 行目の行列積として実行されることになる

7 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

8 HPCS5 5/5/9 5年ハイパフォーマンスコンピューティングと計算科学シンポジウム High Performance Comuting Symosium 5 Relative residual B-orthogonality Maximum error.e-9.e-.e-.e-.e-.e- TOTAL 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

9 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, 日本応用数理学会年度年会予稿集 (). 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 eﬃcient algorithm for the ran-one modiﬁcation of the symmetric eigenroblem, SIAM Journal on Matrix Analysis and Alications, Vol. 5, No., (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

09.pptx

09.pptx 講義内容数値解析第 9 回 5 年 6 月 7 日水理学部物理学科情報理学コース. 非線形方程式の数値解法. はじめに. 分法. 補間法.4 ニュートン法.4. 多変数問題への応用.4. ニュートン法の収束性. 連立次方程式の解法. 序論と行列計算の基礎. ガウスの消去法. 重対角行列の場合の解法項目を変更しました.4 LU 分解法.5 特異値分解法.6 共役勾配法.7 反復法.7. ヤコビ法.7.