分散メモリ型スパースソルバの開発と評価

固有値計算のための三重対角化処理の並列化と
キャッシュ向け最適化
2004年11月22日
名古屋大学計算理工学専攻
山本有作
目次
１. 固有値計算の概要
２. 並列化とキャッシュ向け最適化の基本的技法
３. 三重対角化処理の並列化
４. 三重対角化処理のキャッシュ向け最適化
５. 今後の展開
１. 固有値計算の概要
• 固有値計算とその応用
• 固有値・固有ベクトル計算の流れ
• 各部分の演算量
• 三重対角行列 T に対する固有値・固有ベクトル解法
固有値計算とその応用
• 対象とする問題
– 標準固有値問題 Ax = λx （A：実対称またはエルミート行列）
• 対象とする計算機
– 共有メモリ型並列計算機
– 分散メモリ型並列計算機
– SMPクラスタ
• 固有値計算の応用分野
– 分子計算，統計計算，構造解析など
固有値・固有ベクトル計算の流れ
実対称行列 A
計算内容
計算手法
三重対角化
Q*AQ = T
(Q: 直交行列)
ハウスホルダー法
三重対角行列 T
三重対角行列の
固有値・固有ﾍﾞｸﾄﾙ計算
Tyi =λi yi
Tの固有値 {λi }，
固有ベクトル {yi }
逆変換
Aの固有ベクトル {xi }
xi = Qyi
QR法
DC法
二分法・逆反復法
Dhillon のｱﾙｺﾞﾘｽﾞﾑ
逆変換
各部分の演算量
• N×N の実対称行列 A に対し，M 組の固有値・固有ベク
トルを求める場合
– ハウスホルダー法による三重対角化： (4/3)N3
– 三重対角行列の固有値・固有ベクトルの計算： O(NM) ～ O(N3)
（解法および固有値の分布により大きく異なる。）
– 逆変換： 2N2 M
M≪N の場合は，三重対角化の演算量が支配的
本報告ではこの部分に絞って並列化・最適化の手法を述べる。
三重対角行列 T に対する固有値・固有ベクトル解法
解法
概要
計算量
並列化一部の固有値
・固有ﾍﾞｸﾄﾙ
のみの計算
QR法
相似変換によりTを
対角行列に変換
O(N2)
+ O(N3)
容易
不可
DC法
（分割統治
法）
Tを再帰的に半分の
ｻｲｽﾞの行列に分割
O(N2)
～ O(N3)
容易
不可
二分法・
逆反復法
二分法で固有値，逆
反復法で固有ﾍﾞｸﾄﾙ
を計算
O(NM) 非自明
+ O(NM2)
Dhillon の
アルゴリズム
逆反復法の改良
（dqds，RRR，twisted
dqds = differential qd分解などの利用）
with shift
O(NM)
+ O(NM)
容易
RRR = Relatively Robust Representation; T をLLt の形で表現して計算を進める手法
可
可
２. 並列化とキャッシュ向け最適化の基本的技法
• 共有メモリ型並列計算機
• 分散メモリ型並列計算機
• キャッシュメモリの有効利用
• BLASの利用による高性能化
並列計算機のアーキテクチャ
共有メモリ型並列計算機（SMP）
SMPクラスタ
地球シミュレータ
Itaniumサーバ
Power Mac G5
分散メモリ型並列計算機
日立 SR11000
日立 SR2201
PCクラスタ
Power Mac G5 クラスタ
共有メモリ型並列計算機
• 構成
– 複数のプロセッサ（PU）がバスをキャッシュ PU0
通してメモリを共有
– PUはそれぞれキャッシュを持つ。
バス
PU1
PU2
PU3
メモリ
• 特徴
– メモリ空間が単一のためプログラミングが容易
– PUの数が多すぎる（＞20個）と，アクセス競合により性能が低下
• PU間同期
– 他PUの計算結果を使って計算を
行う場合は，同期が必要
PU0
PU1
PU2
PU3
時間
DO IP = 0, 3 内積計算の例
DO I = 1, 25
S(IP) = S(IP)
+ A(IP*25+I)*B(IP*25+I)
END DO
PU間同期
END DO
S = S(0)+S(1)+S(2)+S(3)
共有メモリ型並列計算機向けの高性能化手法
• PU間の負荷分散均等化
– 各PUの処理量が均等になるよう
処理を分割
– 密行列向けのアルゴリズムでは
比較的容易
キャッシュ PU0
PU1
PU2
PU3
バス
メモリ
• PU間同期の削減
– 同期には通常，数百サイクルの時間が必要
– 並列粒度の増大による同期回数の削減が性能向上の鍵
• キャッシュメモリの有効利用
– 演算順序の変更等により，キャッシュ中のデータの再利用性を向上
– PU間でのアクセス競合を防ぎ，アクセスを高速化
分散メモリ型並列計算機
• 構成
– 各々がメモリを持つ複数のPUをキャッシュ PU0
ネットワークで接続
メモリ
– PUはそれぞれキャッシュを持つ。
• 特徴
PU1
PU2
ネットワーク
– 数千～数万PU規模の並列が可能
– PU間へのデータ分散を意識したプログラミングが必要
• PU間でのデータ転送
– 他PUの計算結果，あるいは他PUの持つデータを使って計算を
行う場合は，PU間でのデータ転送が必要
PU3
分散メモリ型並列計算機向けの高性能化手法
• PU間の負荷分散均等化
– 共有メモリ型の場合と同様
• データ転送の削減
キャッシュ
PU0
PU1
PU2
PU3
メモリ
ネットワーク
– 転送には通常，数千サイクルのセットアップ時間が必要
– データ１個の転送には，演算１回の数十倍の時間が必要
– アルゴリズムとデータ分散方法の工夫により，データ転送量・転送
回数を削減することが性能向上の鍵
• キャッシュメモリの有効利用
– 演算順序の変更等により，キャッシュ中のデータの再利用性を向上
– 相対的に遅い主メモリへのアクセスを削減し，計算を高速化
SMPクラスタ
• 構成
– 複数の共有メモリ型並列計算機
（SMP）をネットワークで接続
• 特徴
PU0
PU1
PU2
メモリ
PU3
PU0
PU1
PU2
メモリ
PU3
PU0
PU1
PU2
メモリ
ネットワーク
– 各ノードの性能を高くできるため，比較的少ないノード数で高性能
を達成できる。
– プログラミングは，ノード内部の計算では共有メモリ型並列機とし
て，ノードをまたがる計算では分散メモリ型並列機として行う。
• 高性能化手法
– ノード内部の計算では共有メモリ型並列機向けの手法を適用
– ノードをまたがる計算では分散メモリ型並列機向けの手法を適用
PU3
キャッシュ向け最適化の基礎
データ転送速度
非常に大
レジスタ
演算器
8～128 ワード
データ転送速度大
キャッシュ
数100Kバイト～数Mバイト
データ転送速度小
主メモリ
数100Mバイト～数Gバイト
• 演算器と主メモリの速度差は，年々増大。
• キャッシュにデータがある間に演算をまとめて行う（データの再利用）
ことにより，主メモリのデータ転送速度の遅さをカバーできる。
BLASの利用による高性能化
• BLASとは
– Basic Linear Algebra Subprograms の略
– 個々のマシン向けにチューニングした基本行列演算のライブラリ
• BLASの種類
– BLAS1: ベクトルとベクトルの演算
• 内積 c := x t y
• AXPY演算 y: = ax + y など
– BLAS2: 行列とベクトルの演算
• 行列ベクトル積 y := Ax
• 行列のrank-1更新 A := A + xyt
– BLAS3: 行列と行列の演算
• 行列乗算 C := AB など
= A ×
A = A
＋
×
C = A × B
BLASにおけるデータ再利用性と並列粒度
• BLAS1
– 演算量： O(N)，データ量： O(N)
– データ再利用性：なし
– 並列粒度： O(N/p) （N：ベクトルの次元，p：プロセッサ台数）
• BLAS2
– 演算量： O(N2)，データ量： O(N2)
– ベクトルデータのみ再利用性あり
– 並列粒度： O(N2/p)
=
A ×
A = A
＋
×
• BLAS3
– 演算量： O(N3)，データ量： O(N2)
– O(N)回のデータ再利用が可能
– 並列粒度： O(N3/p)
C = A × B
演算をできる限り BLAS2，3で行うことが高性能化のポイント
現在利用可能な最適化BLAS
• プロセッサメーカーの提供するBLAS
– Intel MKL，AMD ACML，IBM ESSL など
– メーカーによっては共有メモリ向け並列化版もあり。
• ATLAS（Automatically Tuned Linear Algebra Subprograms）
– 対象プロセッサに合わせて自動的に性能を最適化するBLAS
– http://math-atlas.sourceforge.net/ より入手可能
– 共有メモリ向け並列化版もあり。
• GOTO BLAS
–
–
–
–
テキサス大学オースチン校の後藤和茂氏によるBLAS
http://www.cs.utexas.edu/users/flame/goto/ より入手可能
この３種のBLASの中では，多くの場合最高速
共有メモリ向け並列化版も（一部）あり。
３. 三重対角化処理の並列化
• ハウスホルダー法による三重対角化
• 共有メモリ向けの並列化
• サイクリック列分割による並列化
• Scattered Square 分割による並列化
• 性能評価
ハウスホルダー法による三重対角化
• 鏡像変換による列の消去
ベクトル
– 鏡像変換 H = I – a uut
• Hは対称な直交行列
• 与えられたベクトルの第１成分以外をゼロにする。
左からH
を乗算
• 第 k ステップでの処理
0
0
0
0
0
左からH
を乗算
0
右からH
を乗算
k
非ゼロ要素
ゼロにしたい部分
影響を受ける部分
各ステップでの処理の詳細
•鏡像変換ベクトルの作成
0
σ= sqrt(d td)
u = (d1 – sgn(d1)σ, d2, … , dN – k )
α= 2 /∥u∥2
0
C
d
•H による両側からの乗算
k
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 の対称性を利用し， 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
0
0
d (k )1
C (k )
d (k )
k
ピボット列
ピボット行
ハウスホルダー法の特徴
• 演算量
– 行列サイズが N のときの演算量：約 (4/3)N3
• 行列ベクトル積の演算量：
約 (2/3)N3
• rank-2更新の演算量：
約 (2/3)N3
• データの再利用性
– 演算はほとんどBLAS2のため，行列データの再利用性なし
– キャッシュマシンでは高い性能が期待できない。
ブロック化など，アルゴリズム上の工夫が必要
• 並列粒度
– BLAS2のため O(N2/p) であり，比較的高い。
共有メモリ向けの並列化
• 基本方針
– 行列ベクトル積および rank-2更新の部分について，BLAS2レベ
ルで並列化
– 並列BLASを使えば，並列化は極めて容易
• PU間でのアクセス競合の回避
– BLAS2はデータ再利用性がないため，
バスを通した共有メモリへのアクセスが
頻発
– アクセス競合による性能低下を回避
するには，後に述べるキャッシュ向け
最適化技法が有効
キャッシュ
PU0
PU1
PU2
バス
メモリ
PU3
分散メモリ向けの並列化 (1-1)
• ブロックサイクリック列分割による並列化
（Dongarra et al., 1992）
0
0
– 行列を幅 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)
(5) (3)と(4)の結果の集計，β(k) の計算，q(k) の計算
×
＝
部分和の集計
／ﾌﾞﾛｰﾄﾞｷｬｽﾄ
× ＝
全対全通信により
結果ﾍﾞｸﾄﾙを分配
(6) 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）
0
0
– 行列を b×b のブロックに分割し，ブロックに対して
大きさ √p× √ p のプロセッサテンプレートを周期的
に割り当てる。
鏡像変換
ベクトル
C (k )
k
• 各ステップの処理
(1) 鏡像変換ベクトル u(k) およびα(k) の作成
・第 k 列を担当するPU群が（通信しながら）実行
(2)
u(k)
およびα(k) のブロードキャスト
u (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) ）。
部分和を横方向に集計し，
結果を対角PUに集める。
(4) C (k) の上三角部分と u(k) との積
・各PUは自分の行列要素を使って部分和ベクトルを計算
・縦方向に集計し，対角PUに集める（O(N log p /√p) ）。
(5) (3)と(4)の結果の集計とブロードキャスト
×
＝
部分和を縦方向に集計し，
結果を対角PUに集める。
・対角PUにおいて，(3)と(4)の結果を集計
・結果の部分ベクトルを，横方向のPUにブロードキャスト（O(N log p /√p) ）
・結果の部分ベクトルを，縦方向のPUにブロードキャスト（O(N log p /√p) ）
分散メモリ向けの並列化 (2-3)
• 各ステップの処理（続き）
(6) β(k) の計算，q(k) の計算
(7) rank-2更新 C (k) = C (k) – u(k)q(k)t – q(k)u(k)t
・各PUで独立に計算
• 両並列化方式の比較
– 各PUの計算時間
• 各ステップの計算時間は両方式ともO(N2 / p) （p に反比例して減少）
– 通信時間
• ブロックサイクリック列分割は O(N log p) （pとともに増加）
• ブロックScattered Square 分割は O(N log p /√p) （pとともに減少）
PU台数が増えるほど，ブロックSS分割が有利
サイクリック列分割とSS分割の性能比較
10
5
性能予測モデルによる推定性能
SS分割
Mflops
10
4
ｻｲｸﾘｯｸ列分割
10
10
3
2
10
0
10
1
10
2
10
3
10
4
Num. of Proc.
超並列ではScattered Square分割が優位
SMPクラスタ上での性能評価
ブロックScattered Square 分割による並列化（MATRIX/MPP HZES1MDP使用）
SR8000/F1（SMPクラスタ）の１～16ノードで評価。キャッシュ向け最適化あり。
N=2000 ～ 32000 のエルミート行列での性能
100
90
Gflops
•
•
•
80
70
60
50
40
30
20
10
0
1000
1832[sec]
16-nodes
4-nodes
902[sec]
1-node
447[sec]
10000
Matrix size
100000
４. 三重対角化処理のキャッシュ向け最適化
• Dongarra のアルゴリズム
• Bischof のアルゴリズム
• Wu のアルゴリズム
• 性能・精度評価
Dongarraのアルゴリズム
• 原理（Dongarra et al., 1992）
– 複数（L’ 本）のピボット列・行を溜めておき，後でまとめて更新
– rank-2更新の部分を行列乗算（BLAS3）にでき，キャッシュマシン
で有効
– LAPACK，ScaLAPACK 等で採用され，広く使われる。
0
U(K*L’)
0
C(K*L’)
= C((K–
1)*L’)
d(k) d(k+1) d(k +2) d(k +3)
–
×
Q(K*L’)t
Q(K*L’)
–
U(K*L’)t
×
rank-2L’ 更新（BLAS3使用）
Dongarra のアルゴリズム
[Step 1] K = 1からN /L’ まで以下の[Step 2,3,12] を繰り返す。
[Step 2] U ((K–1)*L’) = φ，Q ((K–1)*L’) = φ（0行0列の行列）
[Step 3] k = (K–1)*L’+1からK*L’ まで以下の[Step 3-1～8]を繰り返す。
[Step 3-1] 部分ハウスホルダ変換：
d(k) := d(k) – U(k–1)(Q(k–1)t)k–(K–1)*L’ – Q(k–1)(U(k–1)t)k–(K–1)*L’
[Step 3-2] 内積 σ(k) = sqrt(d(k)t d(k))の計算
[Step 3-3] 鏡像変換ベクトル作成：
u(k) = (d(k)1 – sgn(d(k)1)σ(k), d(k)2, … , d(k)N– k )
[Step 3-4] 規格化定数の計算： α(k) = 2 /∥u(k)∥2
[Step 3-5] 行列ﾍﾞｸﾄﾙ積：p(k) =α(k) (C (k) – U(k–1)Q(k–1)t – Q(k–1)U(k–1)t) u(k)
[Step 3-6] ベクトルの内積： β(k) =α(k) u(k)tp(k) /2
[Step 3-7] ベクトルの加算： q(k) = p(k) –β(k)u(k)
[Step 3-8] U (k) = [U(k–1)| u(k) ]，Q (k) = [Q(k–1)| q(k) ]
[Step 12] 行列のrank-2L’ 更新：
C (K*L’) = C ((K–1)*L’) – U (K*L’)Q (K*L’)t – Q (K*L’)U (K*L’)t
BLAS-3使用
Dongarra のアルゴリズムの特徴と問題点
• データの再利用性
– 演算量の半分を占める rank-1 更新は BLAS3 （行列乗算）化
– しかし，残り半分を占める行列ベクトル積はそのまま
後者がネックとなり，キャッシュマシンではピークの10～25%程
度の性能しか達成できない場合が多い。
行列ベクトル積の部分も BLAS3に置き換えられるアルゴリズム
が必要
• 計算精度
– 広く利用され，計算精度については十分な確認が行われている。
Bischof のアルゴリズム
• 原理（Bishof et al., 1993, 1994）
– 密行列 A をまず半帯幅Lの帯行列 B に変換し，次にこの帯行列
を三重対角行列 T に変換する。
次数 N
半帯幅 L
帯行列化
約 (4/3)N3
A
0
村田法
O(N2L)
0
B
0
0
T
• 三重対角化を２段階で行うことの利点
– 帯行列への変換は，BLAS3のみを使って実行可能。
– 帯行列から三重対角行列への変換の演算量はO(N2L)であり，
前半部に比べてずっと小さい。
Bischof のアルゴリズム
ブロックベクトル
ブロック鏡像変換によるブロック列の消去
– ブロック鏡像変換 H = I – UαUt
• Hは直交行列
• 与えられたブロックベクトルの
第１ブロックを上三角行列にし，
それ以外をゼロにする。
左からH
を乗算
第kステップでの処理
0
0
左からH
を乗算
非ゼロ要素
0
0
0
ゼロにしたい部分
右からH
を乗算
0
影響を受ける部分
Bischof のアルゴリズム
[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 を
選べない可能性あり
Wu のアルゴリズム
• 原理（Wu et al., 1996）
– Bischofのアルゴリズムにおいて，複数のピボット列・行を溜めて
おき，後でまとめて更新
– L’ 本溜めることで，rank-2L 更新の部分を rank-2LL’ 更新にでき，
データ再利用性を更に向上可能
後半部（三重対角化）の演算量を増やさずに，前半部（帯行列
化）の性能を向上できる可能性あり
４. 性能・精度評価
• 性能評価
– テスト例題
• N = 480，960，1920，3840 のフランク行列の三重対角化
• Dongarra，Bischof，Wu の３つのアルゴリズムの性能を比較
– 評価マシン
• Opteron (1.6GHz)
• Alpha 21264A (750MHz)
・Xeon (2.0GHz)
・Ultra SPARC III (750MHz)
– 評価方法
• 演算量を(4/3)N3 とし，三重対角化に要した時間より各アルゴ
リズムのMFLOPS値を算出して比較
• 各アルゴリズムのブロックサイズ L，L’ は，性能が最大にな
るよう決定
• BLASルーチンとしては，全マシンで ATLAS 3.6.0 を使用
Opteron (1.6GHz) 上での性能
1800
Performance (GFLOPS)
1600
1400
L=24, L’=4
Dongarra
Bishof
Wu
L=48
1200
1000
800
L’=32
600
400
200
0
480
960
1920
3840
Matrix size
Wu の方法は Dongarra の方法に比べて約２倍の性能を達成
N = 3840 のとき，Wu の方法はピークの50%以上の性能を達成
Xeon (2.0GHz) 上での性能
1400
Performance (GFLOPS)
1200
L=24, L’=4
Dongarra
Bischof
Wu
L=48
1000
800
L’=32
600
400
200
0
480
960
1920
3840
Matrix size
Wu の方法は Dongarra の方法に比べて２倍近い性能を達成
Alpha 21264A (750MHz) 上での性能
900
Performance (GFLOPS)
800
700
L=24, L’=4
Dongarra
Bishof
Wu
L=48
600
500
400
L’=32
300
200
100
0
480
960
1920
3840
Matrix size
Wu の方法は Dongarra の方法に比べて約２倍の性能を達成
N = 3840 のとき，Wu の方法はピークの50%以上の性能を達成
Ultra SPARC III (750MHz) 上での性能
500
Performance (GFLOPS)
450
400
L=24, L’=2
Dongarra
Bishof
Wu
L=48
350
300
250
200
L’=32
150
100
50
0
480
960
1920
3840
Matrix size
Wu の方法は Dongarra の方法に比べて２倍以上の性能を達成
精度評価
• テスト例題
–
–
–
–
フランク行列 aij = min(i, j)
２次元ラプラス方程式の５点差分により得られる行列
Harwell-Boeing Sparse Matrix Collection の行列
乱数行列（要素が [0, 1] の一様乱数である対称行列）
• 評価マシン
– Opteron (1.6GHz)
• 評価方法
– 各解法に対し，計算した固有値λi’の相対誤差の最大値
max 1≦i≦N |λi’ – λi| / |λi| を評価
– 比較対象の固有値λiとしては，解析解あるいはオリジナルのハウ
スホルダー法により求めた固有値を使用
フランク行列に対する計算精度
Matrix size
1.00E-06
480
960
1920
3840
Maximum relative error
1.00E-07
1.00E-08
1.00E-09
Dongarra
Bischof
Wu
1.00E-10
1.00E-11
1.00E-12
Bischof 及び Wu の方法の誤差は，Dongarra の方法と同程度
２次元５点差分の行列に対する計算精度
Matrix size
1.00E-10
480
960
1920
3840
Maximum relative error
1.00E-11
1.00E-12
1.00E-13
Dongarra
Bischof
Wu
1.00E-14
1.00E-15
1.00E-16
Bischof / Wu の方法の誤差は Dongarra の方法に比べやや大きい。
しかし，増大の程度は１桁以内に収まっている。
Harwell-Boeing Collection の行列に対する計算精度
Matrix name
1.00E-03
BCSSTK8
BCSSTK12
BCSSTK13
BCSSTK23
BCSSTK24
Maximum relative error
1.00E-04
1.00E-05
1.00E-06
Dongarra
Bischof
Wu
1.00E-07
1.00E-08
1.00E-09
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’ の自動最適化

Download Report