hirayama

Size: px

Start display at page:

Download "hirayama"

ともあきみょうだに
5 years ago
Views:

1 128 ビット 4 倍精度と 160 ビット拡張 4 倍精度演算プログラムの作成平山弘神奈川工科大学自動車システム開発工学科 hirayama@sd.kanagawa-it.ac.jp 工学院大学新宿校舎 28 階第 4 会議室 2013 年 3 月 8 日 ( 金 )

2 なぜ 4 倍精度か 4 倍精度程度の精度では多倍長計算はあまり速くない精度の小さい計算の方が計算精度が大きい計算より需要は多い少ない精度で済む計算を大きな精度で計算するのはムダが大きい

3 4 倍精度数

4 4 倍精度の加算 4 倍精度数を構造体で表す Fortran 90 C++ 言語 Type quad class quad { Real (8) :: high double high ; Real (8) :: low double low ; End Type quad } 負の数は high と low の符号を変える

5 4 倍精度数の加算 Bailey の方法によって非常に簡単に書くことができる Function qq_add(a, b)! Fortran 90! Type (quad), Intent (In) :: a, b Type (quad) :: qq_add Real (8) :: sh, eh, v, ss sh = a%high + b%high v = sh - a%high eh = (a%high-(sh-v)) + (b%high-v) eh = eh + a%low + b%low ss = sh + eh v = ss - sh qq_add%low = (sh-(ss-v)) + (eh-v) qq_add%high = ss End Function qq_add 条件判断の文が含まないので高速に計算できる可能性がある

6 4 倍精度の乗算乗算も短いプログラムで記述できる } quad operator*( const quad &b ) const // C++ 言語 { quad c ; double u1, p1, p2, r1, r2, s1 ; c.high = high * b.high ; u1 = * high ; // = pow(2,27)+1 p1 = u1 - (u1 - high) ; p2 = high - p1 ; u1 = * b.high ; r1 = u1 - (u1 - b.high) ; r2 = b.high - r1; c.low = ((p1*r1-c.high) + p1*r2 + p2*r1) + p2*r2 ; c.low = c.low + (high * b.low +low * b.high) ; s1 = c.high +c.low ; c.low =c.low - (s1 - c.high) ; c.high =s1 ; return c ;

7 4 倍精度数の 2 乗関数乗算のプログラムからやや高速な 2 乗関数が作成できる Function q_square(a) Implicit None Type (quad), Intent (In) :: a Type (quad) :: q_square Real (8) :: u1, p1, p2, s1 q_square%high = a%high*a%high u1 = _dp*a%high p1 = u1 - (u1-a%high) p2 = a%high - p1 s1 = p1*p2 q_square%low = ((p1*p1-q_square%high)+s1+s1) + p2*p2 s1 = a%high*a%low q_square%low = q_square%low + (s1+s1) s1 = q_square%high + q_square%low q_square%low = q_square%low - (s1-q_square%high) q_square%high = s1 End Function q_square

8 4 倍精度数の逆数関数 ( 除算 )

9 4 倍精度数の逆数関数プログラム非常に簡単に書ける Function q_invers(a) Implicit None Type (quad), Intent (In) :: a Type (quad) :: q_invers, x x = to_quad_d(1.0_dp/a%high)! 倍精度で逆数の近似値を計算 q_invers = x + x - a*q_square(x)!newton 法 End Function q_invers これを使えば除算ができる

10 4 倍精度数の平方根逆数関数平方根の逆数も逆数と同様に簡単に求められるこの場合解くべき方程式は 1 a 2 x = 0 近似値は倍精度で計算しそれから Newton 法を使う 3 x = n 1 x + n 0.5( x + n axn) 計算は加減乗算でできる

11 4 倍精度数の関数絶対値等 abs(x) aint(x) 指数対数関数 exp(x) log(x) log10(x) 三角関数 sin(x) cos(x) tan(x) 逆三角関数 asin(x) acos(x) atan(x) 双曲線関数 sinh(x) cosh(x) tanh(x) 逆双曲線関数 asinh(x) acosh(x) atanh(x)

12 4 倍精度数プログラムへ変換 Fortran 90 subroutine add(a,b,c) c =a+b end 4 倍精度への変換倍精度 subroutine add(a,b,c) subroutine add(a,b,c) use def_quad type(quad) a, b,c c =a+b end real(8) a,b,c c=a+b end

13 FORTRAN77 から Fortran90 へ変換コメントの C や * を! に変更継続行を & を使う継続行に変更を a=p+q+r+s * +t+u a=p+q+r+s & +t+u

14 4 倍精度数プログラムへ変換 C++ 言語 void add(double a, double b,double &c){ } c =a+b ; 4 倍精度への変換倍精度 #include def_quad.h void add(quad a,quad b, quad &c){ } c =a+b ;

15 4 倍精度数独自の関数 4 倍精度数への変換 a=to_quad( ) a=to_quad(123) 4 倍精度数から文字列へ変換 write(6,*) to_string(a, e15.7 ) cout << to_string(a, %15.7e ) プログラムの大きさ Fortran 3300 行 C 行

16 4 倍精度数独自の関数 4 倍精度数への変換 a=to_quad( ) a=to_quad(123) 4 倍精度数から文字列へ変換 write(6,*) to_string(a, e15.7 ) cout << to_string(a, %15.7e ) プログラムの大きさ Fortran 3300 行 C 行 ( 使用した多倍長プログラムの違い )

17 4 倍精度プログラム作成に使ったコンパイラ NAG Fortran 32 ビット 64 ビットの Fortran コンパイラ 4 倍精度実数 (double-double) が組み込まれている ( 講演では IEEE 型と述べたが間違い ) G95 32 ビットの無料コンパイラ 80 ビット実数と 4 倍精度 ( 四則と平方根 ) 実数が利用できる Visual C++ 64 ビット C++ コンパイラ Borland C++ 32 ビット C++ コンパイラ 80 ビット実数が使える

18 簡単な拡張 4 倍精度プログラム例 2 次方程式の計算 ( 拡張 4 倍精度 g95) 約 38 桁の精度で計算されている

19 簡単な 4 倍精度プログラム例二重指数型数値積分の計算 (4 倍精度 g95) 大浦氏によるプログラムを変更 I_1=int_0^1 1/sqrt(x) dx I_1= E+0000 err= E-0030 n= 155 I_2=int_0^2 sqrt(4-x*x) dx I_2= E+0000 err= E-0030 n= 131 平方根双曲線関数指数関数を含むプログラムも容易に 4 倍精度に変換できる

20 Fortran 用 Blas の 4 倍精度プログラム Blas 関数の absm 関数の変更例 Fortran 90 FORTRAN 77 function wabsm(n,dx,incx) type(quad), dimension(:) :: dx DOUBLE PRECISION FUNCTION WABSM(N,DX,INCX) DOUBLE PRECISION DX(1),DTEMP type(quad) :: dtemp, wabsm INTEGER I,INCX,M,MP1,N,NINCX integer i,incx,m,mp1,n,nincx WABSM = 0.0D0 wabsm = 0.0d0 DTEMP = 0.0D0 dtemp = 0.0d0 IF(N.LE.0)RETURN if(n.le.0)return IF(INCX.EQ.1)GO TO 20 if(incx.ne.1) then NINCX = N*INCX nincx = n*incx do i = 1,nincx,incx dtemp = dtemp + abs(dx(i)) end do wabsm = dtemp return DO 10 I = 1,NINCX,INCX DTEMP = DTEMP + DABS(DX(I)) 10 CONTINUE WABSM = DTEMP RETURN 20 M = MOD(N,6) end if IF( M.EQ. 0 ) GO TO 40 m = mod(n,6) DO 30 I = 1,M if( m.eq. 0 ) go to 40 do i = 1,m DTEMP = DTEMP + DABS(DX(I)) 30 CONTINUE dtemp = dtemp + abs(dx(i)) IF( N.LT. 6 ) GO TO 60 end do 40 MP1 = M + 1 if( n.lt. 6 ) go to 60 DO 50 I = MP1,N,6 40 mp1 = m + 1 DTEMP = DTEMP + DABS(DX(I)) + DABS(DX(I + 1)) + DABS(DX(I + 2)) do i = mp1,n,6 * + DABS(DX(I + 3)) + DABS(DX(I + 4)) + DABS(DX(I + 5)) dtemp = dtemp + abs(dx(i)) + abs(dx(i + 1)) + abs(dx(i + 2)) & 50 CONTINUE + abs(dx(i + 3)) + abs(dx(i + 4)) + abs(dx(i + 5)) 60 WABSM = DTEMP end do RETURN 60 wabsm = dtemp END end function wabsm

21 Blas の 4 倍精度数プログラム C++ 言語 template 化した Blas を使用 complex の template を使用すれば 4 倍精度の複素数計算が可能

22 4 倍精度数プログラム性能次のような性質の良い行列の逆行列を計算し元の行列と掛けて単位行列を引く各行列要素の絶対値最大の数値を求める行列は次のように与える a i, j 10 = 1 + i i = j i j

23 4 倍精度数プログラム性能 (2) いろいろな大きさの行列について倍精度で計算したその実行時間を示す Intel i7-2700k(3.5ghz) 倍精度計算 (R8) 単位 : 秒 Nag Borland C+ g95 Visual C++ Fortran + CPU 32bit 64bit 32bit 32bit 64bit

24 4 倍精度数プログラム性能 (3) いろいろな大きさの行列について計算したその実行時間を示す Intel i7-2700k(3.5ghz) Fortran 4 倍精度計算 (R16) 単位 : 秒行列のサイ IEEE4 倍精 Double-Double 型ズ度 160bit 型コンパイラ Nag Fortran g95 g95 g95 CPU 32bit 64bit 32bit 32bit 64bit

25 4 倍精度数プログラム性能 (4) いろいろな大きさの行列について計算したその実行時間を示す Intel i7-2700k(3.5ghz) C++ 言語単位 : 秒 4 倍精度計算 (Double-Double) 行列のサイズ Double-Double 型 160bit 型コンパイラ Visual C++ Borland 6.0 Borland 6.0 CPU 64bit 32bit 64bit

26 結論 4 倍精度を持たない 64 ビットの C++ コンパイラを使用する場合 double-double 型の 4 倍精度は非常に高性能で有用である C++ 言語ではインライン展開されるためか Fortran より高性能であった 64 ビット C++ で double-double 型の 4 倍精度数は double 型より 11~14 倍 32 ビットなら 15~18 倍時間がかかる 64 ビット Fortran で 12~30 倍 32 ビットで 35~38 倍であった拡張倍精度を利用した 160 ビット 4 倍精度は double-double の 4 倍精度の 1.5 倍程度の計算時間を要する Fortran 90 や C++ 言語を使うとその変更は単精度を倍精度に変換する程度の作業でできる

27 主な参考文献山田, 佐々, 今村, 町田, 4 倍精度基本線形代数ルーチン群 QPBLAS の紹介とアプリケーションへの応用, 情報処理学会研究報告, Vol.2012-HPC-137,No.23,1-6(2012) 小武守恒, 藤井昭宏, 長谷川秀彦, 西田晃 : 反復法ライブラリ向け 4 倍精度演算の実装と SSE2 を用いた高速化, 情報処理学会論文誌コンピューティングシステム,Vol.1, No. 1, pp (2008)

<4D F736F F F696E74202D F A282BD94BD959C89F A4C E682528D652E707074>

<4D F736F F F696E74202D F A282BD94BD959C89F A4C E682528D652E707074> 発表の流れ SSE を用いた反復解法ライブラリ Lis 4 倍精度版の高速化小武守恒 (JST 東京大学 ) 藤井昭宏 ( 工学院大学 ) 長谷川秀彦 ( 筑波大学 ) 西田晃 ( 中央大学 JST) はじめに 4 倍精度演算について Lisへの実装 SSEによる高速化性能評価スピード収束まとめはじめにクリロフ部分空間法たとえば CG 法は, 理論的には高々 n 回 (n は係数行列の次元数