<4D F736F F F696E74202D F A282BD94BD959C89F A4C E682528D652E707074>

Size: px

Start display at page:

Download "<4D F736F F F696E74202D F A282BD94BD959C89F A4C E682528D652E707074>"

ふさこいくのや
5 years ago
Views:

1 発表の流れ SSE を用いた反復解法ライブラリ Lis 4 倍精度版の高速化小武守恒 (JST 東京大学 ) 藤井昭宏 ( 工学院大学 ) 長谷川秀彦 ( 筑波大学 ) 西田晃 ( 中央大学 JST) はじめに 4 倍精度演算について Lisへの実装 SSEによる高速化性能評価スピード収束まとめはじめにクリロフ部分空間法たとえば CG 法は, 理論的には高々 n 回 (n は係数行列の次元数 ) の反復で収束丸め誤差の影響収束までに多くの反復が必要停滞収束の改善には高精度演算, 例えば 4 倍精度演算が有効であるが計算コスト大反復解法ライブラリ Lis へ 4 倍精度演算を実装 4 倍精度演算に SSE を活用して高速化 double-double 精度演算 FRTRAN は 4 倍精度浮動小数をサポート 5 ビット Bailey 倍精度浮動小数を個用いて 4 倍精度を実現 double-double 精度 a = a.hi + a.lo, a.hi > a.lo ビット IEEE754 準拠の 4 倍精度の表現形式 double-double 精度の表現形式仮数部 5 ビット仮数部ビットビット仮数部 5 ビット基本演算 ( 表記 ) fl(a+ b) :a + bの倍精度浮動小数加算 err(a + b): a + b = fl(a + b) + err(a + b) 乗算 a bについても同様 4 倍精度四則演算のための倍精度基本演算 s=fl(a+b), e=err(a+b) を計算. a >= b FAST_TW_SUM(a,b,s,e) { s = a + b e = b - (s - a) s=fl(a+b), e=err(a+b) を計算. TW_SUM(a,b,s,e) { s = a + b v = s - a e = (a - (s - v)) + (b - v) a を a=h+l に分割 ( 仮数部を 6bit ずつに ) SPLIT(a,h,l) { t = * a h = t - (t - a) l = a - h p=fl(a b),e=err(a b) を計算 TW_PRD(a,b) { p = a * b SPLIT(a,ah,al) SPLIT(b,bh,bl) e = ((ah*bh-p)+ah*bl+al*bh)+al*bl

2 4 倍精度加算乗算 Lis への実装加算 a = b + c. a=(a.hi,a.lo), b=(b.hi,b.lo), c=(c.hi,c.lo) ADD(a,b,c) { TW_SUM(b.hi,c.hi,sh,eh) TW_SUM(b.lo,c.lo,sl,el) FAST_TW_SUM(sh,eh,sh,eh) FAST_TW_SUM(sh,eh,a.hi,a.lo) 乗算 a = b c MUL(a,b,c) { TW_PRD(b.hi,c.hi,p,p) p = p + (b.hi * c.lo) p = p + (b.lo * c.hi) FAST_TW_SUM(p,p,a.hi,a.lo) Lis 反復解法ライブラリさまざまな反復法前処理が組合わせて使える反復解法を 4 倍精度演算に置き換え行列ベクトル積 (matvec) ベクトルの内積 (dot) ベクトルおよびその実数倍の加減 (axpy) Lis への実装の条件ユーザーインタフェイスに影響を与えることなく 4 倍精度演算が利用できるよう以下のようにする係数行列 A, 右辺ベクトル b は倍精度解ベクトル x の入出力は倍精度, 反復解法中では 4 倍精度反復解法中のベクトル, スカラーは 4 倍精度混合精度演算関数が必要 matvec, dot, axpy 主な処理は積和演算. (Floating Multiply-Add) 4 倍精度の積和演算関数 dot と axpy で利用. D 混合精度 (4 倍精度と倍精度 ) の積和演算関数 matvec で利用 D の実装 4 倍精度 vs 倍精度 a = a + b c (a,b,c) { TW_PRD(b.hi,c.hi,p,p) p = p + (b.hi * c.lo) p = p + (b.lo * c.hi) FAST_TW_SUM(p,p,p,p) TW_SUM(a.hi,p,sh,eh) TW_SUM(a.lo,p,sl,el) FAST_TW_SUM(sh,eh,sh,eh) FAST_TW_SUM(sh,eh,a.hi,a.lo) D a = a + b c D(a,b,c) { TW_PRD(b,c.hi,p,p) p = p + (b * c.lo) FAST_TW_SUM(p,p,p,p) TW_SUM(a.hi,p,sh,eh) TW_SUM(a.lo,p,sl,el) FAST_TW_SUM(sh,eh,sh,eh) FAST_TW_SUM(sh,eh,a.hi,a.lo) ポアソン方程式を 5 点中心差分で離散化した行列 ( 次数,,) 行列格納形式は CRS 前処理なし BiCG 法を 5 回反復 FRTRAN4 倍精度 matvecでを利用 matvec で D を利用倍精度実行時間 ( 比 )

3 SSE(Streaming SIMD Extension) SSE は Intel Pentium4 に搭載された x87 命令に代わる高速化命令 8bit のデータに対して SIMD 処理 64bit 倍精度浮動小数なら同時につの演算 SSE の組込み関数を利用 SSE 化のつのアプローチ A: 四則演算関数の SSE 化演算の依存関係が存在関数 _SSE, D_SSE を用意 B: dot, axpy, matvec を含めた SSE 化段のループアンローリングを行うとが回個同時に積和演算をする関数 _SSE, D_SSE を用意 A: 四則演算関数の SSE 化依存関係のため約 5% 程度しか SSE の pd 命令で処理できない _SSEのADD:pdで処理できる部分 xb = _mm_loadu_pd((double*)b); xc = _mm_loadu_pd((double*)c); xs = _mm_add_pd(xb,xc); xe = _mm_sub_pd(xs,xb); xt = _mm_sub_pd(xs,xe); xc = _mm_sub_pd(xc,xe); xb = _mm_sub_pd(xb,xt); xb = _mm_add_pd(xb,xc); xe = _mm_unpackhi_pd(xb,xb); xc = _mm_unpackhi_pd(xs,xs); _SSEのADD :pdで処理できない部分 xt = xs; xb = _mm_add_sd(xb,xc); xt = _mm_add_sd(xt,xb); xs = _mm_sub_sd(xt,xs); xb = _mm_sub_sd(xb,xs); xb = _mm_add_sd(xb,xe); xc = _mm_add_sd(xt,xb); _mm_store_sd(&a->hi,xc); xc = _mm_sub_sd(xc,xt); xb = _mm_sub_sd(xb,xc); _mm_store_sd(&a->lo,xb); B:dot,axpy,matvec を含めた SSE 化 %SSE の pd 命令で処理できる bh = _mm_loadu_pd(&y[i]); bl = _mm_loadu_pd(&yl[i]); sh = _mm_add_pd(bh,ch); th = _mm_sub_pd(sh,bh); t = _mm_sub_pd(sh,th); ch = _mm_sub_pd(ch,th); bh = _mm_sub_pd(bh,t); bh = _mm_add_pd(bh,ch); sl = _mm_add_pd(bl,p); th = _mm_sub_pd(sl,bl); t = _mm_sub_pd(sl,th); p = _mm_sub_pd(p,th); _SSE の ADD の部分 bl = _mm_sub_pd(bl,t); bl = _mm_add_pd(bl,p); bh = _mm_add_pd(bh,sl); th = sh; th = _mm_add_pd(th,bh); sh = _mm_sub_pd(th,sh); bh = _mm_sub_pd(bh,sh); bh = _mm_add_pd(bh,bl); sh = _mm_add_pd(th,bh); _mm_storeu_pd(&y[i],sh); sh = _mm_sub_pd(sh,th); bh = _mm_sub_pd(bh,sh); _mm_storeu_pd(&yl[i],bh); 性能評価テスト行列疎行列反復解法において前処理なしの BiCG 法を用いて比較 Lis の倍精度 FRTRAN で作成した 4 倍精度 Lis の 4 倍精度 ( を利用 ) Lis の 4 倍精度 (_SSE を利用 A) Lis の 4 倍精度 (_SSE を利用 B) 行列 A ポアソン方程式を5 点中心差分で離散化行列 A Toeplitz 行列 γ 右辺ベクトルbはすべてを初期ベクトルx はすべてを収束判定基準は相対残差ノルム - γ γ 3

4 CPU Clock LD Cache L Cache Memory S Compiler 計算環境 Linux.4.smp.8GHz 8KB 5KB GB Intel C++ 9. Intel Fortran 9. pteron Linux.6.4smp.GHz 64KB MB GB 最適化オプションは -3 が記述されている C ファイルは浮動小数の最適化を抑制 (-mp) スピードテストテスト行列 A で 5 回反復した結果次数倍精度 FRTRAN _SSE _SSE 算術平均 pteron 次数倍精度 FRTRAN _SSE _SSE 算術平均スピードテストのまとめ _SSEの性能比較対象 pteon.9.46 _SSEの性能比較対象 pteon..7 FRTRAN 倍精度.7.5 スピードテスト (N= 6 ) pteron FRTRAN 99. (4) 79.8 () 68.9 (8.) 7. (9.9) _SSE 53.4 (6.4) 48.4 (6.8) _SSE 34.9 (4.) 34.9 (4.9) 倍精度 8.4 () 7. () スピードテストのまとめ ( 続 ) 倍精度倍精度 ( 比例 ) _SSE _SSE( 比例 ) 次数倍精度演算は, から徐々に大きく増加, で LD 8KB, で L 5KB から外れる 4 倍精度演算は次数にほぼ比例して増加 4 倍精度演算はデータのロードとストアよりも演算の割合が大きい収束テスト上でテスト行列 A( 次数,) の結果 (, 反復回数, 残差ノルム ) 倍精度 FRTRAN _SSE γ 実行時間反復回数残差実行時間反復回数残差実行時間反復回数残差 E E E E E E E E E E E E E- γ =. までは両方ともに収束している γ =.3 以上では倍精度は停滞.4 倍精度は収束完全な 4 倍精度と _SSE は同等の収束性 4

5 新演算精度を変えたリスタート倍精度で解いた解 ( あるいは途中の値 ) を初期値として 4 倍精度で解く for(k=;k< 最大反復回数 ;k++) { 倍精度演算の反復解法 if( nrm<restart_tol ) break; for(k=k+;k< 最大反復回数 ;k++) { 4 倍精度演算の反復解法 nrmは相対残差ノルム restart_tolはリスタートするための判定基準精度を変えたリスタートの結果上でテスト行列 A( 次数,) γ=.3, restart_tol を -6 とした場合 4 倍精度リスタート ( 倍 :.45,4 倍 :3.98) 反復回数 3 4 ( 倍 :35,4 倍 :69) 4 倍精度で解いた場合と比べて.67 倍速い倍精度で解いた場合と比べて倍まとめ SSE を用いた 4 倍精度演算を実装 Lis へ double-double 精度演算を実装係数行列 A, 右辺ベクトル b は倍精度解ベクトル x の入出力は倍精度, 反復解法中では 4 倍精度反復解法中のベクトル, スカラーは 4 倍精度四則演算関数の SSE 化 dot, axpy, matvec を含めた SSE 化まとめ ( 続き ) スピードテストより完全な 4 倍精度と比較して平均 4.33 倍段のアンローリングにより約倍の高速化倍精度の約 4.5 倍程度収束テストより完全な 4 倍精度とほぼ同等の収束特性収束の改善を図る新たな手法が可能となった倍精度と 4 倍精度の混合精度の反復解法まとめ ( 続き ) 高精度演算の展望 ILU 前処理などの逐次的な前処理を利用しなくとも少ない反復回数で収束する可能性データアクセスよりも演算の割合が大きく並列化には適している並列化によって倍精度との性能差は縮まりより現実的になる今後の課題 4 倍精度をうまく活用して収束性を向上, 安定して解けるロバストなライブラリの開発 5

tabaicho3mukunoki.pptx

tabaicho3mukunoki.pptx 1 2 はじめに n 目的 4倍精度演算より高速な3倍精度演算を実現する l 倍精度では足りないが4倍精度は必要ないケースに欲しい l 4倍精度に比べてデータサイズが小さい Ø 少なくともメモリ律速な計算では4倍精度よりデータ転送時間を減らすことが可能 Ø PCIeやノード間通信がボトルネックとなりやすい GPUクラスタ環境に有効か n 研究概要 l DD型4倍精度演算 DD演算に基づく3倍精度演算