4 倍精度基本線形代数ルーチン群 QPBLAS の紹介 [index] 1. Introduction 2. Double-double algorithm 3. QPBLAS 4. QPBLAS-GPU 5. Summary 佐々成正 1, 山田進 1, 町田昌彦 1, 今村俊幸 2, 奥田洋司 3 1 1 日本原子力研究開発機構システム計算科学センター 2 理科学研究所計算科学研究機構 3 東京大学新領域創成科学研究科
1. Introduction はじめに 京 コンピュータのような超大規模な並列計算機が一般ユーザでも利用できる環境 大規模並列計算機の性能をフルに利用するため シミュレーション手法の大規模化 演算回数が増大するため 本来無限桁である実数を有限桁で打ち切って計算することによる誤差の累積の影響が大きくなる 倍精度演算では有効な精度の結果が得られない可能性あり 地球シミュレータで 375,000 次元の行列の全固有値 固有ベクトルを直接法で計算 精度は数ケタ (SC06 Yamada et. al.)
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)
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 桁 ]
Dailey の 4 倍精度演算アルゴリズム [Bailey の 4 倍精度演算アルゴリズム ] Bailey の double-double アルゴリズムによる加算 C=A+B Bailey の double-double アルゴリズムによる乗算 C=A*B 11 回の演算で実現可能 コンパイラの最適化で 計算順序変更不可 24 回の演算で実現可能 134217729(= 2 27 +1) 上位 26bit, 下位 26bit のデータに分割する際に利用
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 倍精度の下位データ )
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
3.1 DDDOT および DDGEMV の性能評価 DDDOT( 内積計算 ) [α x T y] 計算時間は次元サイズに比例 AMD の方が若干早い DDGEMV( 行列ベクトル積 ) 計算時間 ( 秒 ) [y αax + βy] TRANS= N 計算時間は次元サイズの 2 乗に比例 最適化をした Intel が早い 計算時間 ( 秒 ) 0.14 0.12 0.10 0.08 0.06 0.04 0.02 0.00 0 0.6 0.5 0.4 0.3 0.2 0.1 0 0 Intel O0 Intel O2 Intel O3 200000 500 AMD O0 AMD O2 AMD O3 400000 600000 Intel O0 Intel O2 Intel O3 配列次元 AMD O0 AMD O2 AMD O3 1000 配列次元 800000 1500 1000000 2000
3.1 DDGEMM の性能評価 140 120 100 80 60 40 20 0 0 計算時間 ( 秒 ) Intel O0 Intel O2 Intel O3 500 AMD O0 AMD O2 AMD O3 1000 配列次元 1500 2000 DDGEMM( 行列行列積 ) [C αab + βc] TRANSA= N, TRANSB= N 計算時間は次元サイズの 3 乗に比例 Intel の方が早い DGEMM と DDGEMM の性能比較 (DDGEMM の計算時間 /DGEM の計算時間 ) 最適化 O2,O3(Intel) DGEMM と DDGEMM の計算時間比は最大 10 倍程度 計算時間の増加率 30 25 20 15 10 5 0 10 250 750 1250 1750 配列次元
3.2 QPBLAS の公開 HP QPBLAS の公開 HP http://ccse.jaea.go.jp/ja/download/qpblas.html ライセンスオープンソース (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
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 倍精度の下位データ )
4.1 QPBLAS-GPU の性能評価 DDDOT ( 内積 ) DDGEMV( 行列ベクトル積 ) DDGEMM( 行列行列積 ) の性能評価 計算機環境 Intel/ Linux Processor : Intel Xeon W3565 (3.2GHz) GPU : Tesla C2075(448core, 1.15GHz) OS : RedHat 4.1.2 Fortran Compiler : gfortran 4.1.2 C Compiler : gcc 4.1.2 GPGPU 開発環境 : CUDA4.0 計算時間測定は純粋な GPU の計算時間ではなく ホストメモリーデバイスメモリ間の転送時間等を含んだもの
4.1 GDDDOT の性能評価 [α x T y] 1.0E+00 1.0E-01 CPU GPU 計算時間 (sec) 1.0E-02 1.0E-03 1.0E-04 1.0E-05 1.0E-06 1.0E+03 1.0E+04 1.0E+05 1.0E+06 1.0E+07 データサイズ 3 GPU/CPU 速度比最大 5 倍程度
4.1 GDDGEMV の性能評価 [y αax + βy] 1.0E+01 1.0E+00 CPU GPU 1.0E-01 計算時間 (sec) 1.0E-02 1.0E-03 1.0E-04 1.0E-05 1.0E-06 1.0E+01 1.0E+02 1.0E+03 1.0E+04 データサイズ 3 GPU/CPU 速度比最大 10 倍程度 最大演算速度 0.16GFlops( メモリ間の転送時間を含む )
4.1 GDDGEMM の性能評価 [C αab + βc] 1.0E+03 1.0E+02 CPU GPU 1.0E+01 計算時間 (sec) 1.0E+00 1.0E-01 1.0E-02 1.0E-03 1.0E-04 1.0E-05 1.0E-06 1.0E+00 1.0E+01 1.0E+02 1.0E+03 1.0E+04 データサイズ 3 GPU/CPU 速度比最大 500 倍程度最大演算速度 85.9GFlops( メモリ間の転送時間を含む )
4.2 QPBLAS-GPU の公開 HP QPBLAS- GPU の公開 HP http://ccse.jaea.go.jp/ja/download/ 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. 4 倍精度 Basic Linear Algebra Subprograms: QPBLAS 3. 4 倍精度 Basic Linear Algebra Subprograms on GPU: QPBLAS- GPU 今後 4 倍精度ソフトウェアを中心に拡充を予定
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