LAPACK/BLAS入門

Transcription

1 LAPACK/BLAS 入門 静岡理工科大学情報学部コンピュータシステム学科幸谷智紀

2 本日のメニュー. LAPACK/BLAS 入門 について 2. LAPACK/BLAS とは? 2. LAPACK/BLAS の概略 2.2 LAPACKE/CBLAS 2.3 ベクトル, 行列のデータ型 2.4 LAPACK/BLAS を使うメリット デメリット 3. BLAS の機能と例題 3. BLAS Level, 2, GEMV ベンチマーク 3.3 GEMM ベンチマーク 3.4 [ 例題 ] べき乗法 4. LAPACK の機能と例題 4. 連立一次方程式の直接解法 4.2 行列の固有値問題 5. [ 応用 ] 積分方程式 2

3 . LAPACK/BLAS 入門 について 幸谷著 LAPACK/BLAS 入門 森北出版, 206 年 2 月刊行 BLAS, 連立一次方程式, 標準固有値問題の解説に限定 LAPACK/BLAS 関連の知識を広く浅く実例 ( プログラム ) と実行結果で示す サンプル C プログラムは別配布 [ サポートページ ] 反省点 解説が浅すぎる A5 版 35 ページ 目次第 章 LAPACK/BLAS って何? 第 2 章 LAPACK/BLAS, 最初の一歩第 3 章 BLAS を極める第 4 章 LAPACK ドライバルーチンひとめぐり第 5 章疎行列用の線型計算ライブラリ第 6 章並列化の方法 (Pthread, OpenMP のみ ) 第 7 章 GPU 上の LAPACK/BLAS--- cublas と MAGMA, cusparse 第 8 章非線型問題にもチャレンジ! 3

4 2. LAPACK/BLAS とは? 2. LAPACK/BLAS の概略 2.2 ベクトル, 行列のデータ型 2.3 LAPACK/BLAS では出来ないこと 4

5 2. LAPACK/BLAS の概略 (/4) Fortran で記述された, 線型計算ライブラリ (C ライブラリも同梱 ) IEEE754 単精度, 倍精度の実数, 複素数を要素とするベクトル, 行列をサポート Version 3.8.0(207 年 月リリース ) が最新版 LAPACK ドライバルーチン計算ルーチン BLAS(Level /2/3) aux ハードウェア (CPU, GPU) BLAS(Basic Linear Algebra Subprograms) ベクトル, 行列の基本演算を担当 LAPACK(Linear Algebra PACKage) BLASを基盤とし, より複雑な線型計算を担当 テネシー大学, カリフォルニア大学バークレイ校, コロラド大学デンバー校が 著作権保持 Julie&Julian Langou, Dimmel, Dongarraらによるサポート メンテナンス 5

6 2. LAPACK/BLAS の概略 (2/4) 西暦 カテゴリ 統合型数学ソフトウェア Matlab R Scilab LAPACK/BLAS EISPACK LINPACK ScaLAPACK LAPACK ATLAS BLAS MAGMA cublas Intel MKL プログラミング言語 C++ C FORTRAN66 FORTRAN77 Fortran 並列 分散処理用技術標準 MPI OpenMP PThreads CUDA ハードウェアのトレンド CPU 周波数の増加による高速化 マルチコア メニーコアによる高速化 LINPACK( 連立一次方程式 )+EISPACK( 固有値問題 ) LAPACK/BLAS 様々な派生形ライブラリも登場 APIとして現役 SLATE, MAGMA 等が後継っぽい が,LAPACK/BLASのバージョンアップも継続中

7 2. LAPACK/BLAS の概略 (3/4) BLAS Basic Linear Algebra Subprograms CBLAS BLAS の C バージョン Level ベクトル同士の演算 ( スカラー倍, 内積など ) Level 2 行列 & ベクトル演算 ( 行列 ベクトル積など ) Level 3 行列演算 7

8 2. LAPACK/BLAS の概略 (4/4) LAPACK Linear Algebra PACKage LAPACKE LAPACK の C バージョン ( CLAPACK) ドライバルーチン 問題のタイプ毎に存在 連立一次方程式 線型最小二乗 (LLS, Linear Least Squares) 問題 一般化線型最小二乗問題 標準固有値問題 対称行列の固有値 固有ベクトル計算 非対称行列の固有値 固有ベクトル計算 特異値分解 一般化固有値問題および特異値問題 一般化対称固有値問題 一般化非対称固有値問題 一般化特異値分解 計算ルーチン ドライバルーチンの下支え アクセサリ (aux) 8

9 2.2 ベクトル, 行列のデータ型 (/3) 行目 2 行目 3 行目 *mat_a 行優先 (Row-major) 方式 行目 2 行目 3 行目 列目 2 列目 3 列目 *mat_a 列目 2 列目 3 列目 列優先 (Column-major) 方式 LAPACK/BLAS 全て Fortran90 で記述 列優先 (Column-major) 方式による密行列格納 サブルーチン内の一時変数は全て引数に指定 LAPACKE/CBLAS C/C++ 用 LAPACK/BLAS IF + α 行優先 (Row-major) か列優先を選択利用可 関数内の一時変数を引数に指定しない方式も採用 9

10 API( 関数名 ) の基本命名法 昔のFORTRANの6 文字制限 xyyzzz x 計算精度と実数, 虚数の指定 S... 単精度実数 D... 倍精度実数 C... 単精度複素数 Z... 倍精度複素数 yy 使用する行列のタイプの指定 一般の密行列は xgezzz( 一般行列 ) を使用 zzz 実行される計算の内容を示す文字列 ドライバルーチン xyysv( 連立一次方程式 ) xyye, xyyev( 固有値 固有ベクトル ) 例 ) 倍精度行列ベクトル積 DGEMV (BLAS) cblas_dgemv 倍精度一般行列の連立一次方程式 DGESV (LAPACK) LAPACKE_dgesv 計算ルーチン xyytrf(lu 分解 ), xyytrs( 前進, 後退代入 ) xyyqrf(qr 分解 ) 0

11 LAPACK/BLAS で扱える行列タイプ 一般の非対称 密行列は GE (GEneral) で指定 行列のタイプに合わせて処理時間, 使用メモリの削減ができるものを選択 例 ) 実対称密行列 :SY メモリ量半減 固有値 固有ベクトルは実数

12 派生ライブラリ :Intel Math Kernel CPU コア (Core) レジスタ (Register) RAM CPU Core 0 Core Core 2 Core 3 L (Level) キャッシュ (Cache) L2 キャッシュ L3 キャッシュ 外部記憶 SSD, HDD Thread 0 Thread Thread 2 プログラム (Process) Thread 3 Intel による商用高性能計算ライブラリ LAPACK/BLAS と LAPACKE/CBLAS IF 利用可 疎行列計算機能もサポート キャッシュメモリチューニング,SIMD 命令利用による高速化 マルチコア環境を利用した並列計算機能 // LAPACK(Row-Major or Column-Major) info = LAPACKE_dgeev(LAPACK_COL_MAJOR, 'N', 'V', dim, ma, dim, re_eig, im_eig, NULL, dim, revec, dim); RAM A 2

13 派生ライブラリ :cublasとmagma Shared memory Thread 2 Thread Thread 0 Local memory Block 0 L L L Grid 0 Block S L L L S Block 0 L L L Grid Block S L L L S Block 2 L L L Block 0 L L L Block 2 L L L Block 3 S L L L Grid 2 Block S L L L Block 3 S L L L S S S Block 2 Block 3 S L L L L L L Grid 3 Block 0 Block S L L L L L L Block 2 Block 3 S L L L L L L Device(GPU) S S S Global memory Texture memory Constant memory LAPACK BLAS Level, 2 and 3 MAGMA MAGMA BLAS CUDA cublas PCIe bus RAM Host(CPU) NVIDIA 社提供 GPGPU 開発環境 :CUDA cublas BLAS 互換ライブラリ cusparse 疎行列計算ライブラリ MAGMA CPU & GPU 用高性能 LAPACK/BLAS MAGMA BLAS MAGMA // MAGMA(Column-Major) magma_dgeev('n', 'V', dim, magma_ma, dim, re_eig, im_eig, NULL, dim, magma_revec, dim, h_work, (magma_int_t)lwork_num, &info); 3

14 2.3 LAPACK/BLAS では出来ないこと メモリに入り切れないサイズの行列 ベクトルを扱う問題 オンメモリが基本 単精度計算で 0 進約 7 桁, 倍精度計算で 0 進約 6 桁以上の精度を求める問題 連立一次方程式 条件数の大きな問題 固有値問題 条件数の大きな問題 2 次以上の Jordan ブロックを持つ対角化不可能な行列 4

15 2.4 LAPACK/BLAS を使うメリット デメリット 線型計算をしたいだけなら統合型数値計算環境や動的言語を使うのが better 試行錯誤が楽 昔ほど処理が遅くならない グラフィックスや他のパッケージと連携しやすい Fortran, C/C++ から直接 LAPACK/BLAS を使わざるを得ないケースは有効 C/C++ が好き ( 年寄り?), ブラックボックス嫌い できうる限り高速化したい (MPI, GPU の利用 ) LAPACK/BLAS API を使う他のライブラリを使いたい 線型計算を利用する組み込み機器を開発したい 5

16 3. BLAS の機能と例題 3. BLAS Level, 2, GEMVベンチマーク 3.3 GEMMベンチマーク 3.4 [ 例題 ] べき乗法 6

17 BLAS Level ベクトル ( 行列 ) のコピー 内積 スカラー 各種ノルムの計算 7

18 BLAS Level 2, Level 3 Level 2 行列とベクトルの演算 演算結果が行列になる場合もこちら Level 3 行列積 和 差 スカラー倍は Level で十分 8

19 BLAS Level2(xGEMV) ベンチマーク 9

20 BLAS Level 3(xGEMM) ベンチマーク 20

21 xgemm の比較 (cublas, MAGMA BLAS, IMKL) 2

22 3.4 [ 例題 ] べき乗法 : power_eig.c 22

23 [ 余談 ] Top500.org と LAPACK/BLAS 207 年 月の Top500 4, 9, 0 位に日本スパコン 23

24 反復法で日本は世界一に!=HPCG 疎行列の例 理研の京が No DRIVCAV/cavity04: 37 x 37, 非零成分数 7327 ( 全要素数の約 7.3%) CG 法のアルゴリズム CG 法 BLAS Level, 2 のみで構成 疎行列 ベクトルの演算 (SpMV) の強化が必要 IMKL, cusparse の疎行列演算機能を使おう! ベンチマークは大変 24

25 4. LAPACK の機能と例題 連立一次方程式の直接解法 行列の固有値問題 25

26 連立一次方程式のベンチマーク 連立一次方程式を直接法で解く 自作プログラム LAPACK(DGESV 関数 ) IMKL(DGESV 関数 ) 単 - 倍混合精度反復改良法 (DSGESV) IMKL MAGMA (on CUDA) 26

27 linear_eq.c : DGESV ルーチン使用 // vec_b :=.0 * mat_a * vec_x * vec_b alpha =.0; beta = 0.0; cblas_dgemv(cblasrowmajor, CblasNoTrans, dim, dim, alpha, mat_a, dim, vec_x, inc_vec_x, beta, vec_b, inc_vec_b); // ピボット初期化 pivot = (lapack_int *)calloc(dim, sizeof(lapack_int)); // solve A * X = C -> C := X info = LAPACKE_dgesv(LAPACK_ROW_MAJOR, dim,, mat_a, dim, pivot, vec_b, ); // print printf("calculated x = n"); for(i = 0; i < dim; i++){ printf("%3d -> %3d: ", i, pivot[i]); printf("%25.7e ", vec_b[i]); printf(" n"); } 27

28 直接法とは?=LU 分解と前進 後退代入 xx xx 2 = 6xx + 3xx 2 = 2 2 6xx 2xx 2 = 2 +) 6xx + 3xx 2 = 2 xx 2 = 0 xx 2 = 0 を に代入 3xx = xx = 3 ( 答 ) xx = 3 xx 2 = 0 左辺の未知数が 次式で記述できる複数式の方程式を連立 一次 方程式と呼ぶ 実用上重要な方程式 微分方程式に基づくシミュレーション 建築物の構造解析等々 未知数は大量にあるケースが多い 短い計算時間で求めるにはどうすればよいか? 28

29 LU 分解と前進 後退代入 (2 2 の場合 ) xx xx = 係数行列 A 2 2 解 x. LU 分解 3. 後退代入 3 6/3 3 ( 6/3) ( ) 3 2 LL = 0 3, UU = 2 0 AA = LLLL LL(UUx) = b 2. 前進代入 yy LLy = b 0 2 yy = 2 2 yy = 2yy + yy 2 = 2 yy = yy 2 = 2 + 2yy = 0 UUx = y 3 xx 0 xx = 2 0 3xx xx 2 = xx 2 = 0 xx = + xx 2 = 3 3 xx 2 = 0 ( 答 ) 定数ベクトル b xx = 3 xx 2 = 0 29

30 規模の大きな LU 分解と前進 後退代入 AAxx = bb aa aa 2 aa nn aa 2 aa 22 aa 2nn aa nnn aa nnn aa nnnn 係数行列 A xx xx 2 xx nn = 解 x bb bb 2 bb nn 定数ベクトルb A b U L b L U x LU 分解 nn 3 に比例した計算量 前進 後退代入 nn 2 に比例した計算量 30

31 LU 分解の並列化 Thread 0 Thread Thread 2 Thread 0 Thread Thread 行 ( 横方向 ) 単位でスレッドに計算を割り当てる 計算が進むにつれて, 並列化できる余地が減る L U 3 6 b

32 前進 後退代入の並列化 前進代入 L U b L Solve Ly = b b 後退代入 U y Solve Ux = y スレッド割り当ては市松模様に LU 分解より計算量が少ないので並列化の効果も少ない L U 32 x

33 DGESV 関数ベンチマーク IMKL,2LAPACKE( ソースをそのままコンパイルして利用 ),3 自作ルーチンの順に高速 ブロック化アルゴリズム SIMD 命令利用 大規模問題になるとさらにその差が増す LAPACK/BLAS 相当の線型計算を自作するな! (by 福井大 細田 ) 33

34 混合精度反復改良法 (DSGESV) LAPACK Version 3.2 (2008 年 ) から追加 単精度と倍精度の混合精度バージョン 34

35 混合精度反復改良法 (DSGESV) 計算時間 L 桁計算直接法 LU 分解に要する時間 前進 後退代入に要する時間 S 桁計算直接法 LU 分解 前進後退代入 S-L 桁計算混合精度反復改良法 LU 分解 行列 ベクトル積 前進後退代入 行列 ベクトル積 前進後退代入 反復処理 ( この場合は 2 回反復 ) 単精度計算 (S 桁 ) が高速な GPU では効果的 倍精度計算 (L 桁 ) が遅い環境 & 良条件問題に対しては実用性が高い 35

36 CPU と GPU の処理速度の比較 CPU GPU GPU GPU は CPU より 4 倍近い高速化を達成 DSGESV を使用することで更に倍以上の高速化を達成 36

37 ( 標準 ) 固有値問題 右固有ベクトル x : 左固有ベクトル y : ここで yy HH = yy TT [ 質問 ] 固有値問題と連立一次方程式の解を求める問題との本質的な違いは? [ 質問 2] 対角化不可能な正方行列と対角化可能な正方行列の違いは? 37

38 DSYEV( 実対称 ), DGEEV( 一般 ) 複素数が必要となる計算 実数で閉じている計算 38

39 5. [ 応用 ] 積分方程式 39

40 積分方程式を LAPACK で解く!(/5) 40

41 積分方程式を LAPACK で解く!(2/5) 4

42 積分方程式を LAPACK で解く!(3/5). Gauss-Legendre 積分公式を DSTEQR で求める 2. 離散化した非線型方程式を Derivative-free 解法で求める 42

43 積分方程式を LAPACK で解く!(4/5) 43

44 積分方程式を LAPACK で解く!(5/5) 反復回数 回分, Derivative-free 解法が高速 (n>28) GPU で高速になるかどうかは研究課題 44

45 まとめ 線型計算のグローバルスタンダードである LAPACK/BLAS の概要と機能 ( のごく一部 ) 派生型の高性能ライブラリの性能評価 非線型問題への応用事例 [ サポートページ ] GitHub にてサンプルプログラムは公開中 45