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

応用数理工学特論
線形計算と
ハイパフォーマンスコンピューティング
2008年4月17日
計算理工学専攻張研究室
山本有作
計算科学・工学における大規模シミュレーション
流体計算
タンパク質の
シミュレーション
構造解析
電子状態計算
•数万～数億の自由度を持つ超多自由度の系
•解析のために多大な計算機パワーが必要
最新の計算機を効率的に使うには
• アーキテクチャの理解と最適化手法の習得
– 最新の計算機がどんな原理で高速化を実現しているかを理解
し，その高速性を引き出すための手法を学ぶことが必要
• 線形計算の高速化
– 科学技術計算では，計算の大部分が線形計算に帰着されること
が多い。
• 流体計算 → 連立一次方程式の解法，高速フーリエ変換
• 構造解析 → 連立一次方程式の解法，固有値計算
• 電子状態計算 → 固有値計算
– その高速化は重要な課題
授業の目標
• 最新の計算機で大規模な計算を高速に行うための技法
について学ぶ。
• 線形計算のための高速アルゴリズムについて学ぶ。
授業の構成
1.
2.
3.
4.
5.
6.
7.
8.
9.
高性能計算のための計算機アーキテクチャ
行列乗算の高性能化手法
連立一次方程式の高性能解法（密行列の場合）
連立一次方程式の高性能解法（疎行列の場合）
固有値計算の高性能化手法（非対称行列の場合）
固有値計算の高性能化手法（対称行列の場合）
高速フーリエ変換の高性能化手法
高性能計算に関する最新トピックス
グループ発表
テーマ間の関連
高性能計算のための計算機アーキテクチャ
行列乗算の高性能化手法
連立一次方程式
（密行列）
連立一次方程式
（疎行列）
固有値計算
グループ発表
高速フーリエ
変換
授業の進め方(1)
• レポート
– ２回出題予定
• グループ発表
– 最後の２回の授業を使って行う。
– ２～３人のグループで１つの発表を行う。
– 論文を読んでまとめるか，あるいは講義で勉強したことをもとに自分
たちで数値実験を行った結果を発表。
授業の進め方(2)
• 授業ノート
– 各回の授業内容を担当者（２名）が記録して提出
• TeX，テキスト形式，あるいは Word で
• 提出は授業後１週間以内に
– 山本のホームページに掲載予定
• http://www.na.cse.nagoya-u.ac.jp/~yamamoto
/lectures/hpc2008/hpc2008.html
– ノート担当はレポート提出１回分とみなす
授業の進め方(3)
• スーパーコンピューティング実習との連携
– この講義を受ける人で，OpenMP，MPIを使ったことがない人は，でき
るだけ石井克哉教授担当の講義「計算科学フロンティア特別講義：並
列計算特論」（5/7開講）も受講してください。
– これにより，OpenMP／MPIの使い方が学べるとともに，情報連携基盤
センターのスーパーコンピュータのアカウントが発行され，並列計算の
アルゴリズムなどを実地に試せるようになります。
授業の進め方(4)
• 連絡先
– ３号館南479号室山本まで
– あるいはメールで
[email protected] まで
参考書
• J. Demmel: Applied Numerical Linear Algebra, SIAM, 1997.
• G. H. Golub & C. F. van Loan: Matrix Computations, 3rd Edition,
Johns Hopkins Univ. Press, 1996.
• K. A. Gallivan et al.: Parallel Algorithms for Matrix Computations,
SIAM, 1994.
• K. R. Wadleigh and I. L. Crawford: Software Optimization for high
Performance Computing, Prentice-Hall, 2000.
• B. Wilkinson and M. Allen: Parallel Programming: Techniques and
Applications Using Network of Workstations and Parallel Computers,
Prentice-Hall, 1999.
ハイパフォーマンスコンピューティング（HPC）技術とは
• 大規模な計算を高速かつ高精度に行うための技術
• HPC技術の内容
–
–
–
–
高速・高精度な計算アルゴリズム
単体プロセッサ向けの性能最適化技術
並列化技術
ネットワーク・GRID技術
PC向けプロセッサの
クロック周波数の向上
• 周波数は10年間で約50倍も
向上
• たとえば AMD Opteronプロ
セッサ（1.6GHz）の場合，１
秒間に最高で32億回の浮動
小数点演算を実行可能
（3200MFLOPSのピーク性
能）
それでは，一般的な科学技術計算のプログラムで
は，ピーク性能の何％程度を発揮できているか？
• 例：連立一次方程式を解くためのガウスの消去法
do k=1, n
do i=k+1, n
a(i,k)=a(i,k)/a(k,k)
end do
do
j=k+1, n
do i=k+1, n
a(i,j)=a(i,j)–a(i,k)*a(k,j)
end do
end do
end do
• n=1000のとき，Opteronでの性能はピークの何％？
10%
30%
50%
80%
Opteronでの性能測定結果
• n=1000 での性能は225MFLOPS（ピークの7%）
• n が大きくなるにつれ，性能は低下
Performance (MFLOPS)
3500
3000
Gaussian elimination
peak performance
2500
2000
1500
1000
500
0
100
200
300
400
500
600
700
800
900
1000
n
ピークよりはるか下の性能しか得られない原因
は？
• 最大の要因は，データ転送速度のネック
– 計算を行うには，主メモリに格納されているデータをプロセッサ内
の演算器に転送する必要あり。
– 演算器は十分速いが，データを供給する速度が追い付かない。
• それでは，データ転送速度のネックを解消するにはどう
すればよいか？
まず，プロセッサのアーキテクチャを知る必要がある。
典型的なマイクロプロセッサのメモリ階層
データ転送速度
非常に大
レジスタ
演算器
8～128 ワード
データ転送速度大
キャッシュ
数100Kバイト～数Mバイト
データ転送速度小
ラインサイズ
主メモリ
数100Mバイト～数Gバイト
• 演算器に近い記憶装置ほど高速だが，容量は小さい。
• 演算器と主メモリの速度差は，年々大きくなっている。
性能最適化の原理
• データがレジスタ中にあれば，演算器は最高速度で計算が可能
いったんデータをレジスタに持ってきたら，そのデータに対して必要な
計算を集中して行うよう，計算の順序を変更する。
（データ参照の局所性を高める。）
• キャッシュと主メモリについても，同じ方針で最適化を行う。
性能最適化の具体例
• もっとも単純なアルゴリズムである行列乗算 C=AB に対して最適化
を行う。
最適化前のプログラム
do i=1, n
do j=1, n
sum=0.0d0
do k=1, n
sum=sum+a(i,k)*b(k,j)
end do
c(i,j)=sum
end do
end do
• 性能 = 77.7MFLOPS （Opteron 1.6GHz，n=500）
• ピーク性能の2.4%
レジスタの再利用性を高める最適化（レジスタブロッキング）
•
•
•
•
→
→
→
→
iのループを２倍展開
i, jの各ループを２倍展開
iを４倍，jを２倍展開
i, jの各ループを４倍展開
162.8MFLOPS
240.8MFLOPS
324.7MFLOPS
495.5MFLOPS
500
450
performance
400
350
300
250
200
150
100
50
0
1×1
2×1
2×2
4×2
4×4
キャッシュの再利用性を高める最適化（キャッシュブロッキング）
• 行列を部分行列（１個がL×L）に分割
• 部分行列単位で乗算を行う。
• Lは部分行列３個がキャッシュに格納できるサイズに取る。
do I=1, n/L
do J=1, n/L
CIJ=0
do K=1, n/L
CIJ=CIJ + AIKBKJ
end do
end do
end do
• 演算量は同じだが，主メモリへのアクセス回数が約1/Lに減少する。
キャッシュブロッキングの効果
800
performance
700
600
500
400
300
200
100
0
1
8
20
40
100
200
500
Block size
これまでに説明した最適化手法の限界
• レジスタブロッキングとキャッシュブロッキングの併用により，行列乗
算の性能は800MFLOSまで向上できた。
• しかし，ピーク性能に対しては，まだ25%に過ぎない。
• その理由
– 実際のプロセッサでは，キャッシュが２階層になっている。
– データ転送速度だけでなく，アクセス遅延時間も考慮した最適化が必
要，など。
• 性能を最大限に引き出すには，これらの点も考慮した最適化が必要
BLAS （Basic Linear Algebra Subprograms）
• 行列乗算をはじめとする基本的な行列演算については，BLASと呼
ばれる最適化済みのライブラリが各プロセッサに対して提供されて
いる。
• BLASの機能
– 行列乗算
– 行列とベクトルの積
– ベクトルの和，内積，など。
• ATLAS（Automatically Tuned Linear Algebra Subprograms）
– 自分で自分を最適化するBLAS
– インストール時に，ループ展開のサイズやキャッシュブロッキングのサイ
ズなどの最適値を自分で探し，設定。
ATLASの性能
• n=1000のとき，ピーク性能の95%以上を達成している。
3500
3000
2500
hand-coded
ATLAS
peak performance
2000
1500
1000
500
0
100
200
300
400
500
600
700
800
900
1000
その他の行列乗算の高速化手法
• Strassenのアルゴリズム
– A，B，Cをそれぞれ２×２に分割して乗算を行う。
– 乗算の回数を7/8に削減可能
– Strassenのアルゴリズムで現れる小行列の乗算に対し，再帰的に
Strassenのアルゴリズムを使うことにより，計算量を更に削減可能
– ただし，標準的なアルゴリズムに比べて加減算の回数が増える
ため，丸め誤差は増大
• 再帰的に使うと，さらに誤差が増大
– 誤差に敏感でない問題に限って使うべき。
Strassenのアルゴリズム
C11 C12   A11 A12   B11 B12 




C 21 C 22   A21 A22   B21 B22 
P1 = (A11+A22) (B11+B22)
P2 = (A21+A22) B11
P3 = A11 (B12 – B22)
P4 = A22 (B21 – B11)
P5 = (A11+A12) B22
P6 = (A21 – A11) (B11+B12)
P7 = (A12 – A22) (B21+B22)
C11 = P1 + P4 – P5 + P7
C12 = P3 + P5
C21 = P2 + P4
C22 = P1 + P3 – P2 + P6
まとめ
１. マイクロプロセッサのメモリ階層
– レジスタ，キャッシュ，主メモリ
２. 単体プロセッサ向けの高性能化手法
– データ参照の局所性の向上
• レジスタ有効利用のためのループ展開
• キャッシュ有効利用のためのブロック化
– 行列乗算の例
３. アルゴリズムの工夫による演算量の削減
– 行列乗算のための Strassen アルゴリズム

Download Report