4 倍精度基本線形代数ルーチン群 QPBLAS の紹介 [index] 1. Introduction 2. Double-double algorithm 3. QPBLAS 4. QPBLAS-GPU 5. Summary 佐々成正 1, 山田進 1, 町田昌彦 1, 今村俊幸 2, 奥田洋司

Size: px

Start display at page:

Download "4 倍精度基本線形代数ルーチン群 QPBLAS の紹介 [index] 1. Introduction 2. Double-double algorithm 3. QPBLAS 4. QPBLAS-GPU 5. Summary 佐々成正 1, 山田進 1, 町田昌彦 1, 今村俊幸 2, 奥田洋司"

しなつてらわ
6 years ago
Views:

1 4 倍精度基本線形代数ルーチン群 QPBLAS の紹介 [index] 1. Introduction 2. Double-double algorithm 3. QPBLAS 4. QPBLAS-GPU 5. Summary 佐々成正 1, 山田進 1, 町田昌彦 1, 今村俊幸 2, 奥田洋司日本原子力研究開発機構システム計算科学センター 2 理科学研究所計算科学研究機構 3 東京大学新領域創成科学研究科

2 1. Introduction はじめに京コンピュータのような超大規模な並列計算機が一般ユーザでも利用できる環境大規模並列計算機の性能をフルに利用するためシミュレーション手法の大規模化演算回数が増大するため本来無限桁である実数を有限桁で打ち切って計算することによる誤差の累積の影響が大きくなる倍精度演算では有効な精度の結果が得られない可能性あり地球シミュレータで 375,000 次元の行列の全固有値固有ベクトルを直接法で計算精度は数ケタ (SC06 Yamada et. al.)

3 4 倍精度化の方法 [ 対応策多倍長 (4 倍精度 ) 計算 ] real*16 変数を利用計算時間がかかるため実用的でない多倍長計算用のライブラリ Bailey の double-double アルゴリズム 2 つの倍精度実数 (real*8) を組み合わせることで 4 倍精度演算を実現倍精度演算の組み合わせなので高速に計算可能 a=a.hi+a.lo a.hi : 上位データ a.lo : 下位データ倍精度の仮数部は 52bit で表現しているためこの方法だと 104bit で表現 (real*16 の場合仮数部は 112bit)

4 2. double-double algorithm 2 つの倍精度実数 (real*8) を組み合わせて 4 倍精度演算を実現 double-double 型実数 a (ah,al) ah: 上位データ (real*8), al: 下位データ (real*8) 52bit( 仮数部 ) 52bit( 仮数部 ) 104bit( 仮数部 ) [10 進約 32 桁 ] REAL*16 型実数 112bit( 仮数部 ) [10 進約 34 桁 ]

5 Dailey の 4 倍精度演算アルゴリズム [Bailey の 4 倍精度演算アルゴリズム ] Bailey の double-double アルゴリズムによる加算 C=A+B Bailey の double-double アルゴリズムによる乗算 C=A*B 11 回の演算で実現可能コンパイラの最適化で計算順序変更不可 24 回の演算で実現可能 (= ) 上位 26bit, 下位 26bit のデータに分割する際に利用

6 3. 4 倍精度化 BLAS BLAS (Basic Linear Algebra Subprograms) 線形基本演算のルーチン群 (40 個 ) double-double アルゴリズムで 4 倍精度化 QPBLAS (Quadrature Precision Basic Linear Algebra Subprograms) 倍精度版 [ ルーチン名の接頭辞にdをつけて倍精度を表現 ] call dgemm( TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LCD) 4 倍精度版 [ ルーチン名の接頭辞にddをつけて4 倍精度を表現 ] call ddgemm( TRANSA, TRANSB, M, N, K, ALPHAH, ALPHAL, AH, AL, LDA, BH, BL, LDB, BETAH, BETAL, CH, CL, LCD) (*H : 4 倍精度の上位データ *L : 4 倍精度の下位データ )

7 3.1 QPBLAS QPBLAS の性能評価の性能評価 DDDOT ( 内積 ) DDGEMV( 行列ベクトル積 ) DDGEMM( 行列行列積 ) の性能評価使用する計算機 Intel/ Windows Processor : Intel Core 2 Duo E8400 (3.0GHz) OS : Windows XP Professional Compiler : Intel Fortran 10.0 IA32 AMD/ Linux Processor : Dual Core AMD Opteron Processor 2800 (2.4GHz) OS : Cent OS 4.4 Compiler : gfortran 4.1.0

8 3.1 DDDOT および DDGEMV の性能評価 DDDOT( 内積計算 ) [α x T y] 計算時間は次元サイズに比例 AMD の方が若干早い DDGEMV( 行列ベクトル積 ) 計算時間 ( 秒 ) [y αax + βy] TRANS= N 計算時間は次元サイズの 2 乗に比例最適化をした Intel が早い計算時間 ( 秒 ) Intel O0 Intel O2 Intel O AMD O0 AMD O2 AMD O Intel O0 Intel O2 Intel O3 配列次元 AMD O0 AMD O2 AMD O 配列次元

9 3.1 DDGEMM の性能評価計算時間 ( 秒 ) Intel O0 Intel O2 Intel O3 500 AMD O0 AMD O2 AMD O 配列次元 DDGEMM( 行列行列積 ) [C αab + βc] TRANSA= N, TRANSB= N 計算時間は次元サイズの 3 乗に比例 Intel の方が早い DGEMM と DDGEMM の性能比較 (DDGEMM の計算時間 /DGEM の計算時間 ) 最適化 O2,O3(Intel) DGEMM と DDGEMM の計算時間比は最大 10 倍程度計算時間の増加率配列次元

3.2 QPBLAS の公開 HP QPBLAS の公開 HP http://ccse.jaea.go.jp/ja/download/qpblas.

10 3.2 QPBLAS の公開 HP QPBLAS の公開 HP ライセンスオープンソース (2 条項 BSDライセンス ) [subrou^nes] Level 1 Level 2 Level 3 ddswap ddscal ddcopy ddaxpy dddot ddnorm2 ddsum ddidmax ddrot ddrotg ddrotm ddrotmg ddzdotc ddzdotu ddgemv ddsymv ddtrmv ddtrsv ddsyr ddsyr2 ddgbmv ddger ddsmbv ddtbmv ddtbsv ddzgerc ddzgeru ddzhbmv ddzhemv ddzher ddzher2 ddgemm ddsymm ddsyr2k ddsyrk ddtrmm ddtrsm ddzhemm ddzher2k ddzherk

11 4. QPBLAS-GPU BLAS (Basic Linear Algebra Subprograms) 線形基本演算のルーチン群 (40 個 ) double-double アルゴリズムで 4 倍精度化 CUDA4.0 で GPGPU 用に実装 QPBLAS-GPU (Quadrature Precision Basic Linear Algebra Subprograms on GPUs) 倍精度版 [ ルーチン名の接頭辞にdをつけて倍精度を表現 ] call dgemm( TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LCD) GPU4 倍精度版 [ ルーチン名の接頭辞にgddをつけて4 倍精度を表現 ] call gddgemm( TRANSA, TRANSB, M, N, K, ALPHAH, ALPHAL, AH, AL, LDA, BH, BL, LDB, BETAH, BETAL, CH, CL, LCD) (*H : 4 倍精度の上位データ *L : 4 倍精度の下位データ )

12 4.1 QPBLAS-GPU の性能評価 DDDOT ( 内積 ) DDGEMV( 行列ベクトル積 ) DDGEMM( 行列行列積 ) の性能評価計算機環境 Intel/ Linux Processor : Intel Xeon W3565 (3.2GHz) GPU : Tesla C2075(448core, 1.15GHz) OS : RedHat Fortran Compiler : gfortran C Compiler : gcc GPGPU 開発環境 : CUDA4.0 計算時間測定は純粋な GPU の計算時間ではなくホストメモリーデバイスメモリ間の転送時間等を含んだもの

13 4.1 GDDDOT の性能評価 [α x T y] 1.0E E-01 CPU GPU 計算時間 (sec) 1.0E E E E E E E E E E+07 データサイズ 3 GPU/CPU 速度比最大 5 倍程度

14 4.1 GDDGEMV の性能評価 [y αax + βy] 1.0E E+00 CPU GPU 1.0E-01 計算時間 (sec) 1.0E E E E E E E E E+04 データサイズ 3 GPU/CPU 速度比最大 10 倍程度最大演算速度 0.16GFlops( メモリ間の転送時間を含む )

15 4.1 GDDGEMM の性能評価 [C αab + βc] 1.0E E+02 CPU GPU 1.0E+01 計算時間 (sec) 1.0E E E E E E E E E E E E+04 データサイズ 3 GPU/CPU 速度比最大 500 倍程度最大演算速度 85.9GFlops( メモリ間の転送時間を含む )

4.2 QPBLAS-GPU の公開 HP QPBLAS- GPU の公開 HP http://ccse.jaea.go.jp/ja/download/ qpblas_gpu.

16 4.2 QPBLAS-GPU の公開 HP QPBLAS- GPU の公開 HP qpblas_gpu.html ライセンスオープンソース (2 条項 BSD ライセンス ) [subrou^nes] Level 1 Level 2 Level 3 gddswap gddscal gddcopy gddaxpy gdddot gddnorm2 gddsum gddidmax gddrot gddrotg gddrotm gddrotmg gddzdotc gddzdotu gddgemv gddsymv gddtrmv gddtrsv gddsyr gddsyr2 gddgbmv gddger gddsmbv gddtbmv gddtbsv gddzgerc gddzgeru gddzhbmv gddzhemv gddzher gddzher2 gddgemm gddsymm gddsyr2k gddsyrk gddtrmm gddtrsm gddzhemm gddzher2k gddzherk

5 公開ソフトウェア公開用ソフトウェア一覧 HP http://ccse.jaea.go.jp/ja/download/ software.html 1. 京コンピュータ用固有値計算ライブラリ : EigenK 2.

17 5 公開ソフトウェア公開用ソフトウェア一覧 HP software.html 1. 京コンピュータ用固有値計算ライブラリ : EigenK 2. 4 倍精度 Basic Linear Algebra Subprograms: QPBLAS 3. 4 倍精度 Basic Linear Algebra Subprograms on GPU: QPBLAS- GPU 今後 4 倍精度ソフトウェアを中心に拡充を予定

18 6. Summary double-double アルゴリズムを用いて 4 倍精度化 BLAS (QPBLAS) の作成と公開 1 主要 40 routine 2 倍精度版 gemm との速度比は約 10 倍 GPU 用 4 倍精度化 BLAS (QPBLAS-GPU) の作成と公開 1 主要 40 routine 2 CUDA4.0 3 CPU 版 gemm との速度比は最大 500 倍高速化 8

untitled

untitled A = QΛQ T A n n Λ Q A = XΛX 1 A n n Λ X GPGPU A 3 T Q T AQ = T (Q: ) T u i = λ i u i T {λ i } {u i } QR MR 3 v i = Q u i A {v i } A n = 9000 Quad Core Xeon 2 LAPACK (4/3) n 3 O(n 2 ) O(n 3 ) A {v i }