発表の流れ 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 = 34779. * 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
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.6 8.6 8.7.
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
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.34.4.646.49.3.79.786.6356.488.33.577.79.66783.585.3364.8734.985 6.8664 5.987 3.466 8.44 99.388 68.9387 53.4777 34.9778 算術平均.86557 44.3665 5.9864.85944 7.74985 pteron 次数 倍精度 FRTRAN _SSE _SSE.37.86.669.457.36.357.837.6687.456.37.44.7888.68794.4766.337.67599 7.95484 7.3384 4.8489 3.4769 7.955 79.85377 7.953 48.4959 34.934 算術平均.58384 7.737 5.59738.767 7.749 スピードテストのまとめ _SSEの性能 比較対象 pteon.9.46 _SSEの性能 比較対象 pteon..7 FRTRAN 6.4.4 倍精度.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. () スピードテストのまとめ ( 続 ) 4. 3.5 3..5..5..5...4 倍精度倍精度 ( 比例 ) _SSE _SSE( 比例 ).3.93 4 6 8 次数 倍精度演算は, から徐々に大きく増加, で LD 8KB, で L 5KB から外れる 4 倍精度演算は次数にほぼ比例して増加 4 倍精度演算はデータのロードとストアよりも演算の割合が大きい 収束テスト 上でテスト行列 A( 次数,) の結果 (, 反復回数, 残差ノルム ) 倍精度 FRTRAN _SSE γ 実行時間反復回数 残差 実行時間 反復回数 残差 実行時間反復回数残差.7866 58.84E- 8.69 58.84E- 3.78737 58.84E-..9436 7.3E-.7676 7.3E- 4.583 7.3E-..744 86 3.3E- 7.9793 86 3.3E- 5.6647 86 3.3E-.3 36.534 3.47E- 7.445 3.47E-.4 7.887 55.85E-.3 55.85E- γ =. までは両方ともに収束している γ =.3 以上では倍精度は停滞.4 倍精度は収束 完全な 4 倍精度と _SSE は同等の収束性 4
新演算精度を変えたリスタート 倍精度で解いた解 ( あるいは途中の値 ) を初期値として 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 倍精度 リスタート 7.4 4.43 ( 倍 :.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