1 / 19

応用数理工学特論 線形計算と ハイパフォーマンスコンピューティング

応用数理工学特論 線形計算と ハイパフォーマンスコンピューティング. 第 8 回 名古屋大学　計算理工学専攻 山本有作. １ . 　はじめに. 本研究で対象とする問題 標準固有値問題　 A x = l x A ： n × n 非対称密行列 応用分野 MHD 化学工学 量子化学 流体力学 Cf. Bai and Demmel: A test matrix collection for non-Hermitian eigenvalue problems. 高い信頼性. 非対称行列の固有値計算の流れ. 計算内容. 計算手法. 密行列 A.

evangelia-chloe
Download Presentation

応用数理工学特論 線形計算と ハイパフォーマンスコンピューティング

An Image/Link below is provided (as is) to download presentation Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author. Content is provided to you AS IS for your information and personal use only. Download presentation by click this link. While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server. During download, if you can't get a presentation, the file might be deleted by the publisher.

E N D

Presentation Transcript

  1. 応用数理工学特論線形計算とハイパフォーマンスコンピューティング応用数理工学特論線形計算とハイパフォーマンスコンピューティング 第8回 名古屋大学　計算理工学専攻 山本有作

  2. １.　はじめに • 本研究で対象とする問題 • 標準固有値問題　Ax = lx • A： n×n 非対称密行列 • 応用分野 • MHD • 化学工学 • 量子化学 • 流体力学 • Cf. Bai and Demmel: A test matrix collection for non-Hermitian eigenvalue problems.

  3. 高い信頼性 非対称行列の固有値計算の流れ 計算内容 計算手法 密行列 A QTAQ = H (Q: 直交行列) ハウスホルダー法 ヘッセンベルグ化 ヘッセンベルグ行列H 固有値の計算 QR法 分割統治法 A の固有値 {li} 固有ベクトルの計算 Hui =li ui Hの固有ﾍﾞｸﾄﾙ {ui } 逆変換 逆変換 vi = Qui Aの固有ベクトル {vi }

  4. 各部分の実行時間 Execution time (min) • 全固有値を求める場合の演算量 • ヘッセンベルグ 化：　(10/3) n3 • QR法： 　10n3（経験値） • 演算量（時間）の大部分をQR法が占める。 • QR法の高速化が必要 Origin2000 1PU （R12000，400MHz）上での実行時間 K. Braman et al: “The Multishift QR Algorithm I” より抜粋 Hessenberg 化の時間は推定値

  5. QR法の特徴と高性能計算 高性能計算の条件 QR法の特徴 • 計算の逐次性 　　（bulge-chasing 型演算） 計算の並列性 低いデータ再利用性 • 高いデータ再利用性 　　（キャッシュの有効利用） オリジナルのQR法のままでは高性能化が困難 高性能計算の条件を満たす新しいアルゴリズムが必要 マルチシフトQR法

  6. QR法の基礎 • アルゴリズム（Francis, 1961 & Kublanovskaya, 1961） • 行列 A0から始めて次のように QR 分解と相似変換を繰り返す。 • A0 = Q1R1 • A1 = R1Q1 (= Q1–1A0 Q1 ) • A1 = Q2R2 • A2 = R2Q2 (= Q2–1A1 Q2 = Q2–1Q1–1A0 Q1Q2 ) • 収束定理 • 適当な条件の下で，Ak は（ブロック）上三角行列に収束 • A0 の固有値を絶対値の大きい順に l1, l2, … , lnとすると，対角要素 aijは li に１次収束 • 非対角要素 aij （j < i）は収束率 rij≡ |li| / |lj| で 0 に１次収束

  7. ダブルシフトQR法 • シフトの導入 • Akから固有値 li の近似値 sを引いた行列に対して QR 法を適用 • Ak– s I = Qk Rk • Ak+1 = Rk Qk+ sI (= Qk–1AkQk) rij≡ |li – s| / |lj – s| ≒ 0 より収束が加速 • ダブルシフトQR法 • 共役複素数の固有値ペアに対し，シフト s, sによる2反復をまとめて行う • (Ak– s I)(Ak– s I) = Qk Rk • Ak+2 = Qk–1AkQk • シフトは右下隅の 2×2 行列の固有値に取る 　　　　→　局所的に2次収束 • 複素数の固有値を持つ場合でも，実数演算だけで QR 法を実行可能

  8. 0 Hessenberg 形と Implicit Q定理 • Hessenberg 形の利用 • Akが Hessenberg 形のとき， • QR法の1ステップは O(n2) で実行可能 • Ak+1 も Hessenberg 形 • Implicit Q定理の利用により，QR 分解を陽に行う 　 ことなく1ステップの計算を実行可能 A0を Hessenberg 形に相似変換してからQR法を適用 • Implicit Q定理 • U，Vが直交行列で，G = UtAU， H = VtAV が共に既約な Hessenberg 行列であるとする。 このとき，もし Uと Vの第1列が等しいならば，±1 を要素に持つ対角行列 Dが存在して，V = UD，H = DGD が成り立つ • すなわち，行列 Aを既約な Hessenberg 行列に相似変換する直交行列は，第1列目だけが与えられれば（実質的に）一意に定まる Hessenberg 形 aij = 0 if i > j+1

  9. Implicit shift QR法 Ak 0 • 1ステップの計算 • (Ak– s2I)(Ak– s1I)の第1列を e1の定数倍にするハウスホル 　　ダー変換H0を求める • Ak’ = H0tAkH0 • 直交行列による相似変換を繰り返すことにより， Ak’ を再び Hessenberg 行列に変形する（bulge-chasing） • この変換の性質 • ある直交行列 H による相似変換 • H の第1列は，Qk の第1列に等しい • 上記第3のステップは第2行・第2列以降のみに影響 • H は Akを Hessenberg 形に変換 Implicit Q 定理より，Ak+1 = H tAk H はQR法の1ステップと等価 H0 AkH0 0 0 0 Ak+1 0

  10. 陰的ダブルシフトQR法の演算パターン • バルジ追跡における演算 • 3×3のハウスホルダー変換を左右からかけることにより，bulge を1つ右下に動かす • 演算の特徴 • 並列粒度は O(n) • 並列アルゴリズム：　R. Suda et al.（1999），G. Henry et al. (2002) • QR法の1ステップで，各行列要素は3回のみ更新 データ再利用性が低く，キャッシュの有効利用が困難 左から Hlを乗算 右から Hlを乗算 0 0 0 更新される要素 0にしたい要素

  11. ２.　マルチシフトQR法 • 原理 （Bai & Demmel, 1989） • Ak の右下の m×m行列の固有値 s1，s2 , … , smをシフトとして用い，QR法の m ステップを一度に実行 • (Ak– smI) ・・・ (Ak– s2I)(Ak– s1I) = Qk Rk • Ak+m = Qk–1AkQk • シフト数増加の効果 • 並列粒度は O(m) 倍に増加 • 行列の各要素に対する更新回数も O(m) 倍に増加 • データ再利用性の向上により，キャッシュの有効利用が可能 　　（WY representation を用いた BLAS3化）

  12. マルチシフトQR法におけるバルジ追跡 m = 6 の場合 • 方式１：　大きなバルジ1個を追跡 （Bai & Demmel, 1989） • シフト s1, s2, …, smを含む (m+1)×(m+1) のバルジを導入 • (m+1)×(m+1) のハウスホルダー変換を用いてこれを追跡 • 方式２： 小さなバルジ複数個を追跡 （Braman et al., 2003） • シフトを2つずつ組にし， 3×3 の小さなバルジ m/2 個を導入 • ダブルシフトQR法と同様にして，これらを追跡 • バルジ同士は3行離れていれば，干渉なしに計算が可能 • 密に詰めることで，データ参照の局所性を向上 • 得られる行列は，無限精度演算では方式１と同じ 0 0

  13. 方式１と方式２の収束性比較 • シフトの数と総演算量との関係 • DHSEQR： 方式１（LAPACK） • TTQR： 方式２ K. Braman et al.: “The Multishift QR Algorithm I” (2002) より引用 有限精度演算では，方式２が有利 （Cf. Watkins による理論的解析） 以下では方式２について考える

  14. レベル3 BLAS の利用 • 更新処理の分割 • m/2個のバルジをそれぞれ r行追跡する際， 　 まず対角ブロックのみを更新 • 非対角ブロックは，更新に使ったハウスホルダー変換を1個の行列に蓄積し，後でまとめて更新 　　非対角ブロックの更新で レベル3 BLASが使用可能 • rの決め方 • 演算量の点から，r～3mとするのが最適 • このとき，演算量は方式１の約2倍 　　（非ゼロ構造を利用した場合） 0 バルジ（3×3） 最初に更新 まとめてBLAS3で更新 蓄積されたハウスホルダー変換の非ゼロ構造

  15. マルチシフトQR法の性能 • 流体力学の問題に対する計算時間の例 • PowerPC G5 2.0GHz （1CPU） • 行列は金田研究室 水野氏提供 N=3645 N=11232

  16. ３.　シフト数の最適化 • シフト数mに関するトレードオフ • mを大きくすると，BLAS3 部分の性能は向上 • しかし，シフトの計算量は増加（O(m3)） • 最適な mの値は計算機と問題サイズに依存 • r（1回のバルジ追跡行数）の最適化 • 演算量最小化のためには r～3mが最適 • BLAS3 の実行時間は必ずしも演算量に比例しないため，実行時間の最適値はこれとは異なる m，rに対する自動チューニングの必要性 実験的に求めた mの最適値 （Origin2000上，Braman et al. より引用） そこで，性能予測モデルを用いてシフト数 mの自動最適化を行うことを考える

  17. PowerPC G5上での予測結果（1CPU） • m を変えたときの相対実行時間（最小値を1に規格化） • モデルは mを変えたときの相対実行時間の変化を定性的に再現 　　→　シフト数の自動最適化に使える可能性あり • ただし，絶対時間では20～40%の誤差 相対実行時間 相対実行時間 予測結果 実測結果

  18. PowerPC G5上での予測結果（2CPU） • m を変えたときの相対実行時間（最小値を1に規格化） • モデルは mを変えたときの相対実行時間の変化を定性的に再現 • 絶対時間の誤差は1CPUの場合と同程度 相対実行時間 相対実行時間 予測結果 実測結果

  19. Opteron 上での予測結果（4CPU） • m を変えたときの相対実行時間（最小値を1に規格化） • モデルは mを変えたときの相対実行時間の変化を定性的に再現 • 絶対時間の誤差は PowerPC G5 の場合と同程度 相対実行時間 相対実行時間 予測結果 実測結果

More Related