多倍長精度演算の性能評価 日時 年 月 日 :3-: 場所工学院大学新宿校舎 8 階第 4 会議室 高エネルギー加速器研究機構 濱口信行 hgu@post.kek.jp // 第 回多倍長精度計算フォーラム
. はじめに 計算センター => ユーザプログラムの実行効率は何 % です よく出ています or 改善してください 実行性能 = 演算量 / 実行時間実行効率 = 実行性能 / 理論性能 ユーザ実行時間 = 演算量 / 実行性能に関心 例えば, フーリエ変換と FFT: フーリエ変換のほうが FFT より実行時間がかかるのに, 実行性能がよい. 両者の実感が一致する 演算量が確定, 理論性能の定義が明確 最近では実数の倍精度演算の場合のみ明確化されている. 今回の発表の動機 4 倍精度の性能がどうこうと言う話を聞くが そもそも 4 倍精度以上の多倍長精度演算の性能とは何? // 第 回多倍長精度計算フォーラム
. 多倍長精度演算 演算方式には, 浮動小数点演算方式と整数演算方式がある.. 浮動小数点演算方式 多倍長精度変数を複数個の倍精度変数の和で表す 表現できる数値範囲は,ieee754 の倍精度と同じ. : n 倍精度変数,,..., : 倍精度変数 有効ビット数 = 仮数部のビット数 n倍精度 5 n n n n 3 n 4... n 6ビット 59ビット ビット ビット 乗加算命令がある場合とない場合の処理に分かれる コンパイラの最適化方式に性能が大きく依存する // 第 回多倍長精度計算フォーラム 3
* *., : :,, ;,,,, 注意が必要とならないように 最適化オプションにより乗算乗加算命令ありならない様に注意が必要加算減算ともに最適化オプションで括弧をはずした演算順序に減算加算倍精度変数倍精度変数の加減算乗除算の結果を つの倍精度変数の和で表す 基本の演算 d c d c t t c d c t c t t c d c t c d c d c * * * * * * * 34779. 7 c d t t r t t t r t c d r 乗算乗加算命令なし // 4 第 回多倍長精度計算フォーラム
演算量 n 3 では 加減算 乗算 n に比例 3 n に比例 n 個の変数の絶対値が大きい順にならべるソート処理が必要.. 整数演算方式 多倍長精度変数は ieee754-8 の 4 倍精度の仮数部を変えたものであらわす 表現できる数値範囲は ieee754-8 の 4 倍精度と同じ. : p倍精度 p 3 変数 符号部 ビット 指数部 5ビット 仮数部 3 p 6ビット 有効ビット数 p 4 3ビット p 6 77ビット p 8 4ビット p 倍精度 演算量 加減算 乗算 pに比例 p に比例 符号部 指数部 仮数部 ビット 5ビット 3*p-6ビット p=4:,p=6: 76,p=8: 4 // 第 回多倍長精度計算フォーラム 5
これらから多倍長精度演算の性能に関しては, 以下の つの影響が大きい. 浮動小数点演算方式 dd 形式 か整数演算方式 ieee 形式 か. 乗加算命令 f を使用するかしないか. またこれらの性能には, 計算機アーキテクチャーとコンパイラの最適化手法が大きく左右する. これらを, 具体的な例題を用いて計算方法と精度の解説を加え, 精度上問題のないものに対して性能を評価した. // 第 回多倍長精度計算フォーラム 6
3.Rup s 例題 7767., 3396. f 333.75 5.5 6 より 6 4 5.5 8 f 54767.873659946... 669 倍精度演算 4倍精度演算 f f.85967743.7639453... 全く異なった結果となる.! // 第 回多倍長精度計算フォーラム 7
原因 A 333.75 B 5.5 8 6 6 797346689636347549485. 7973466896363475494848. 4 A+B の演算で ビットの桁落ちが発生する事 有効桁数 進 37 桁以上の多倍長精度演算が必要となる. 幾つかの計算機でこの計算を,, 回実行したときの実行時間 秒 コンパイラは Intel fortrn,ieee 形式は整数型 3 ビット演算. 計算機アーキテクチャーとコンパイラの最適化手法が大きく左右する事を示している. // 第 回多倍長精度計算フォーラム 8
xeon3.6ghz d75.ghz V8. V9. 精度 ieee dd 精度 ieee dd 5 倍精度 6 5 倍精度.53 6 倍精度 7.63 3.44 6 倍精度.79 3.4 7 倍精度 8.69 7 倍精度 3.8 8 倍精度.9 7.8 8 倍精度 3.97 7.7 d85.6ghz x5355.66ghz V9. V9. 精度 ieee dd 精度 ieee dd 5 倍精度.4 5 倍精度.75 6 倍精度.36.63 6 倍精度.99.95 7 倍精度.5 7 倍精度. 8 倍精度 3.34 6.8 8 倍精度.83 4.5 xeon と x5355 同じインテル系でもコンパイラ V8. と V9. で様子が異なる. 同じコンパイラ V9. でもインテル系とオプテロン系て様子が異なる. // 第 回多倍長精度計算フォーラム 9
C プログラム : 拡張倍精度変数を つつなげると, 有効ビット数 3, 有効桁数 進 39 桁. <= 正しく計算できる条件を満たす. gcc 標準でコンパイル.,, 回実行したときの実行時間 秒 は以下の表 dd 形式の C プログラム 精度 e543 x557 周波数.66GHz.93GHz long doule*.359.8446 doule*3.53868.47 doule*4.533.8675 // 第 回多倍長精度計算フォーラム
4. 数学関数 評価関数 : rctnx を,, 回実行したときの実行時間 秒 の比較 dd 形式 : 各精度に特化したものと一般化したもの Ieee 形式 : 整数型 3 ビット演算, 整数型 64 ビット演算 4. rctnx の計算概略 n n x rctn x n n x rctn x rctn x x rctn x sign これにより, rctn x x x[ rctn ] 4 x x の計算ができれば良い. j x w x に対して, w j,,...,5 として, v とすると, 6 6 xw rctn x rctn w rctn v から求める. v 6 x でのrctn x はテーラー展開で求める. 6 j 注 付近の値でのxに対し, vを求める際の桁落ち防止の処理が必要 6 // 第 回多倍長精度計算フォーラム
4. 性能評価に関して テーラー展開での項数は 4 n 有効ビット数 に対して, n となるnの最小値を取っている 8 倍精度までは 有効ビット数 は ieee 形式の値をとり, 倍精度以上では, dd 形式と ieee 形式で差が出るため, 性能測定では,dd 形式と ieee 形式の平均に設定した. 具体的項数 4 倍精度 5, 5 倍精度 9, 6 倍精度, 7 倍精度 6, 8 倍精度 9 倍精度 35, 倍精度 4, 4 倍精度 49, 6 倍精度 56 次ページ以降の結果では, 倍精度以上の演算精度が必要な場合は性能面より,ieee 形式を使用するのが良い事を示している. // 第 回多倍長精度計算フォーラム
4.3 性能測定結果 rctnx の,, 回実行時間 秒 d75.ghz+v9. 演算精度 dd 一般化 ieee3 dd 特化 ieee64 4 倍精度.94555.677 5 倍精度 6.648 5.7533 6 倍精度.94 7.78386 6.574 7.3688 7 倍精度.3637 9.45556 8 倍精度 9.4875 5.54964.9398 3.3498 倍精度 6.569 7.9676 倍精度 3.8389 48.9456 4 倍精度 96. 67.9367 6 倍精度 34.477 9.789 rctnx の,, 回実行時間 秒 e543.666ghz+v9. 演算精度 dd 一般化 ieee3 dd 特化 ieee64 4 倍精度.555764.483 5 倍精度 4.5973 3.77446 6 倍精度 6.846 5.864 3.559 4.798 7 倍精度 8.63688 6.483 8 倍精度 5.45765.644 9.94489 8.767 倍精度 34.8487 9.674 倍精度 63.7993 3.43 4 倍精度.9 45.469 6 倍精度 78.8458 65.5643 // 第 回多倍長精度計算フォーラム 3
rctnx の,, 回実行時間 秒 corei73.ghz+v. 演算精度 dd 一般化 ieee3 dd 特化 ieee64 4 倍精度.486.9846 5 倍精度 3.58456 3.9543 6 倍精度 4.83565 4.669.34648 3.83547 7 倍精度 6.7897 5.3888 8 倍精度.88 9.373575 6.94944 6.77977 倍精度 4.7584 5.96657 倍精度 48.73459 4.537 4 倍精度 86.98378 34.338 6 倍精度 35.8693 46.5797 rctnx の,, 回実行時間 秒 x557.93ghz+v. 演算精度 dd 一般化 ieee3 dd 特化 ieee64 4 倍精度.4848.963853 5 倍精度 3.39495.877563 6 倍精度 4.6897 4.9736.4663 3.597453 7 倍精度 6.585 5.798 8 倍精度.378 8.3573 6.65989 6.3534 倍精度 6.534 4.54879 倍精度 5.898.7654 4 倍精度 89.85934 3.746 6 倍精度 45.4339 4.3387 // 第 回多倍長精度計算フォーラム 4
5 行列積の計算 dd 形式での演算量 計算式 : 倍精度浮動小数点演算命令を単位 精度 加減算 乗算 fなし 乗算 fあり 4 倍精度 4 6 倍精度 7 78 364 8 倍精度 49 75 998 注 f 命令ありの場合,f 命令は flop となるので,flop を単位としたときの演算量は 内の値となる. SR6/M での性能モニターで測定した値 MFLOP n=4 精度 f なし 対倍精度 計算 f あり 対倍精度 計算 4 倍精度 39788 8.5 7.5 579.5 6 倍精度 7353 54.7 5.5 857 37.5 33.5 8 倍精度 353 9.6 6979 79 73.5 注 倍精度演算は 47 4 倍精度以上の多倍長精度演算の演算量はモニター値などを使用するのが良い. // 第 回多倍長精度計算フォーラム 5
// 第 回多倍長精度計算フォーラム 6 6.Two loop s C M se D u z y x u z y x M xy u z zu y x u y z x u z y x E u z y x u z y x u z y x C dudzdydx CD s S x y x z y x 5 4 3 5 4 3,,,, ; S の一般式今回の扱ったケースと解析近似解 i s y n y n s n n s S n n ] ln 3 [ 4 3 3!,,,, ;
精度検証結果 解析近似解 =3.844383e- 倍精度演算結果 =3.844383e- 4 倍精度演算結果 =3.844383e- s, 5, 3 4 s, 5, 3 4 解析近似解 =.667394555973 倍精度演算結果 =.667397783997 積分誤差.9D-8 積分法 : 二重指数関数型積分法加速法 : ε- 算法 s= の場合 5 回反復 // 第 回多倍長精度計算フォーラム 7
性能評価結果 s, 5, 3 4 二重指数関数型積分 n=576 GPGPU 理論性能 88GFLOPs 倍精度 3.83 秒 58GFLOPs ソースから倍精度演算が単位 4 倍精度 5.57 秒 35.6GFLOPs ソースから4 倍精度演算が単位 性能はプログラムでカウントしているので, 最適化による演算削減分は含まれていない. 倍精度 /4 倍精度 =6.3 倍 実行時間比 SR6/M 理論性能 98GFLOPs 倍精度.777 秒 性能モニター 346GFLOPs 倍精度演算が単位 4 倍精度 33.393 秒 性能モニター 7GFLOPs 倍精度演算が単位 倍精度 /4 倍精度 =8. 倍 実行時間比 倍精度演算数 =346*.777=475GFLOP 4 倍精度演算数 =7*33.393=44GFLOP 演算数比 4 倍精度演算数 / 倍精度演算数 =.4 最適化による演算削減より, 演算量はプログラムでカウントした場合より, 性能モニターでの値は少なくなっている. // 第 回多倍長精度計算フォーラム 8
4 倍精度以上の多倍長精度演算では, プログラムから演算量を算出して, 例えば 4 倍精度演算性能を QGFLOPs 等と表現する方法がある. 特に実行効率を重視するサイトで 4 倍精度演算 および多倍長精度演算を使用するユーザには. モニター等で倍精度演算量を採取し,4 倍精度演算および多倍長精度演算を同じベースで評価する方法が有効. // 第 回多倍長精度計算フォーラム 9
7. 特異点のない infr ox 積分式 I x x y D dzdydx D sxy tz x y z x y x y z x y e z x y f s 5, t 5, f 5, e.5, 3 解析近似解 I s t f 相対誤差 7 s ln ln 6 以下 t e f f // 第 回多倍長精度計算フォーラム
解析近似解の算出方法 // 第 回多倍長精度計算フォーラム
// 第 回多倍長精度計算フォーラム
性能測定結果 cpu :TK 筑波システム AMD qud-core Opteron 8 シリーズ Brcelon node: ピーク性能 47GFLOPs,6MPI/node 積分法 : 二重指数関数型積分 サイズ n=4 要求精度 : 解析近似解と 進 6 桁まで一致する事 実行時間一覧表 秒 精度 64MPI 8MPI 56MPI 形式 4 倍精度 8.57 4.39.9 dd 4 倍精度 7.569 3.784.89 ieee 6 倍精度.9.668 5.53 dd 8 倍精度 88.97 46.34 3.79 dd 8 倍精度 79.43 39.759.67 ieee 倍精度 37.786 9.357 59.69 ieee 6 倍精度 355.7 77.48 88.869 ieee 3 倍精度 37.44 65.35 36.99 ieee 注 この解析近似解と 進 6 桁まで一致する様にするため,dd 形式の4 倍精度のみ, 分点のテーブル作成部分を以下の様に変更している 解析近似解.3567368 D-6 修正前 x3i=expy/expy+exp-y, result=.35653377937485d-6 修正後 x3i=.q/.q+exp-.q*y,result=.356736644969533d-6 // 第 回多倍長精度計算フォーラム 3
8. 特異点のある infr ox I x x y dzdydx D D sxy tz x y z x y x y z x y z x y f e s 5, t 5, f 5, e.5, 3 特異点のある場合の計算では,s< で求めた式に s> の値を代入した値と積分結果を比較する. 3 よりこの解析近似解は十分な精度をもっているまた有理化は,. s s iで行い, を求める加速法は 算法を使用している. 理由 これはD D iで有理化すると虚数部の計算で異なる3つの値が得られ, プログラムは正しくても物理現象と合わない結果が出ることによる. // 第 回多倍長精度計算フォーラム 4
領域分割 実数部の値を求めるため, 積分領域を 4 つの領域に分割し, 各領域での二重指数関数型積分法で積分値を求めて最終結果を得ている. 分割の方法は積分領域を [,]*[,]*[,] にして z 軸を最外側にする. また.5 に対する対称性を利用して積分領域を [,] から [,.5] にする. 積分実行前に 4 e e s 4 s 4 4 Z の値がきまると Az の値が定まる. 4A z S 4A z s 4A z 4A z e e を求める を計算する. この値より, x, y [, ] [,] [,] [,.5] の領域を以下の4つに分割する. 領域 [, ] [,] 領域 [, ] [,] 領域 3[,] [, ] 領域 4[,] [,.5] 注 領域分割では, を使用している. 注 は領域の和を表している. // 第 回多倍長精度計算フォーラム 5
実数部の積分計算結果が要求される相対誤差.% 以下をみたすのに要する時間は 64 スレッド MPI で以下の様に変遷している. 実行時間一覧表 秒 CPU SR XM TK SR6 理論性能 537.6 844.8 588 98.49 GFLOPs 4 倍精度 68 589 44 4 6 倍精度 346 4 37 8 8 倍精度 38 458 78 34 infr ox 多倍長精度演算結果 実数部 解析近似解 -.3567368D-6 4 倍精度 -.3565943D-6 相対誤差.3% 6 倍精度 -.356683D-6 相対誤差.4% 8 倍精度 -.35644997D-6 相対誤差.8% 以上の結果をみるとこの問題には dd 形式 6 倍精度が最も適していると言える. // 第 回多倍長精度計算フォーラム 6
9. まとめ 4 倍精度, 多倍長精度演算ではモニター等により, 実行時の倍精度演算数を採取し それをもとに同じ基準で性能をみるのが良い dd 形式と ieee 形式では 4 倍精度から 8 倍精度まで性能面では特に優劣はない. 3 倍精度以上の演算精度が必要な場合は ieee 形式が良い. 44 倍精度演算で精度不足の場合,8 倍精度ではなく,dd 形式の 6 倍精度を使用するのが良い. // 第 回多倍長精度計算フォーラム 7