PowerPoint プレゼンテーション

Transcription

1 スーパーコンピュータ有効活用の手引き ( 平成 年度 - 平成 7 年度 平成 年度から高性能計算の扉に記載しました各種計算機基本性能調査 各種計算機アプリケーション性能比較 多倍長計算手法のまとめとして性能関連事項を各種計算機性能調査まとめと 精度関連事項を多倍長計算手法まとめの 冊にしました 性能関連事項と精度関連事項は異なる視点からの検討が必要な部分が多くあるため この 冊を列挙する体裁を取っています 論文リストはスライド 8 以降にまとめて掲載しています スーパーコンピュータ使用するケースは主に以下の つのケースに分類できます ( 大学等での計算機教育や学問としての計算機論の展開 ( 研究所等での論文作成の手段として ( 業務を行う公的機関や民間会社での手段として 一般的には ( と ( のサイトでは目的と手段が明確になっていますが ( のサイトでは目的と手段が不明確で逆転している事が多々見られます 一例として過去ある大学で逆行列を求めるサブルーチン群を学内で公開した所, 参照回数の 5% 以上がクラメルの公式で求めるサブルーチンだった事があります ( や ( の分野の人から見れば, 時間もかかり, 計算精度からみて問題のあるものを使用するのはおかしく見えますが 教育という面からみると, 連続変数を扱う数学と離散変数を扱う計算機の差や, 計算機の演算量や精度を知る上では目的にかなっていると言えます

2 今回は業務を行う場合から見る事にし シュミレーションにおいて許容される精度で正しい結果を速く実行するためにスーパーコンピュータを使用するという観点からまとめています 例えばスーパーコンコンピュータの性能と言う場合性能 (FLOPs= 演算量 (FLOP/ 実行時間 (sec を重視しますが ここでは実行時間 (sec で中心に見ます また実ジョブでは使用されない様な条件での性能は検討から外しています すなわち中心におくのは 計算機学の研究ではなく, シュミレーションを実施する側にすると言う事です 理由としては幾つか例をあげますが ( フーリエ変換を使用するアプリケーションでは FFT を使用すれば, 実行時間は短縮されますが 性能の定義では FFT の FLOPs 値は低下します ( 最適化オプションを上げれば 一般的には実行時間は短縮されますが FLOPs 値は低下します ( スーパーコンピュータを使用するアプリケーションでは必ずしも浮動小数点演算だけとはかぎりません などです また実行時間に関しては PC やサーバーから移植する場合には互換性などの問題も含め経過時間が問題となりますので測定結果はすべて経過時間としています たとえ CPU での測定でも CPU 時間ではなく経過時間としています

3 各種計算機性能調査まとめ 平成 年度 - 平成 7 年度 目次. はじめに. 使用計算機. スーパーコンピュータ. アクセラレータ. サーバー. 基本演算. メモリバンド幅測定. rup s 例題.. ieee75-8 系.. power 系.. 超多倍長演算. 行列積, ストラッセン行列積.. 行列積計算... SR6/M システム... システム... BG/Q システム... 各種演算精度の行列積計算.. ストラッセンの行列積. 実アプリケーション. 重力多体問題 (N 体問題.. SR6/M システム.. BG/Q システム.. 各種計算機性能比較... 演算量可変ケース... 演算量一定ケース

4 各種計算機性能調査まとめ 平成 年度 - 平成 7 年度. 分子動力学計算.. SR6/M システム詳細.. 各種計算機性能比較. 大規模疎行列の反復計算.. 対称問題.. 非対称問題. QCD 計算.5 量子モンテカルロ法による物性スペクトル計算.5. SR6/M システム.5. ieee75-8 形式.6 ファインマンループ積分計算.6. infra vt 計算.6. infra bo 計算.6. assless 計算 次元積分計算.6.. s 計算 次元積分計算 次元積分計算.6.5 loop 計算.6.6 loop 計算

5 . はじめに 性能を考える場合, 実行時間をもとにしても, 一般的に以下の事が言えます ( 扱う問題がある程度特定される場合とそうでない場合の性能は大きく異なる場合がある ( カタログ性能と実行時間は必ずしも 比例しない ( 最適化オプションの効果には注意する必要がある 具体例を挙げると以下の様なものがあります ( 複素変数の絶対値計算 c a ib c a b a, bが絶対値がある範囲である事がわかっている場合は, 実数演算でつの乗算,つの加算,つの平方根計算 で済む.c 6 i 6 の場合,cは表現可能な数値 で, 結果も表現可能な値となるが 同じ計算をすると演算途中でオーバーフローが発生します これを防ぐため, 関数ルーチンでは大小判定や除算等が使用されるのが一般的で数倍の性能差が出ます

6 ( カタログ性能は計算機に備わる演算器が全て同時に動作する事を前提としているが, 実際のアプリケーションでは必ずしもそうはならない 倍精度演算と 倍精度演算の性能比は SR では約 倍と言われてきたが SR6 では 倍精度演算には SIMD 命令が適用されないためその差は大きくなり 5- 倍となっている ( 最適化オプションでは結果が変わらない範囲での最適化と 結果が変わる可能性がある最適化のレベルがある 特に, 並列実行の際には, 精度に敏感なプログラムで総和, 内積演算を含む場合には性能以前に結果の検証が必要になります

7 . 使用計算機. スーパーコンピュータ ( SR/K( ノード プロセッサ :power5 周波数 :.GHz CPU コア数 6 論理コア数 6 理論最大性能. GFLOPs メモリ容量 GB メモリアーキテクチャー Flat Meory Interleave L キャッシュ Off-Chip 6MB/ コア演算器 / 物理コア乗加算器 つ ( SR6/M プロセッサ :power7 周波数 :.8GHz CPU コア数 ( 物理的,6( 論理的 理論最大性能 98.8 GFLOPs メモリ容量 56GB メモリアーキテクチャー NUMA,(6 論理コア単位で flat SIMD(Single Instruction Multiple Data をサポートする VSX 機構付き L キャッシュ On-Chip MB/8 コア演算器 / 物理コア乗加算器 つ

8 ( SR6/ SR6/ は周波数が.GHz で他は SR6/M と同じです 演算性能だけみれば,SR6/M ノードは SR6/ の 6% 性能向上版ともいえます またメモリ性能は SR6/M の半分です この 機種では,( の SR6/M がメインで ( の SR/K は互換性, 実行性能比を, SR6/X はメモリ性能の影響と, st の効果を見ています (BG/L 周波数.7GHz ノード 論理性能 57. GFLOPs L キャシュ KB( データ, 命令,L キャシュ KB,L キャシュ MB (5 BG/Q 周波数.6GHz ノード 6core 論理性能.8GFLOPs L キャシュ 6/6KB (Core L MB (node Main storage 6GB (Core St=,, (,(5 に関してはメインは (5 で ( は互換性を見ています

9 (6TK 筑波システム AMD quad-core Opteron 8 シリーズ (Barcelona node: ピーク性能 7GFLOPs,6MPI/node

10 . アクセラレータ (HD587 GPU カード型番 :ATI RadeonHD587 メモリ : GDDR5, GB, 5.6 GB/s ホストインタフェース : PCI Epress. 6strea processing unit: 個 ( 演算プロセッサ 動作周波数 : 85 MHz ピーク性能 ( 倍精度 : 88 Gflops (HD697 GPU カード型番 :ATI RadeonHD697 メモリ : GDDR5, GB, 76 GB/s processing unit: 6 個 ( 演算プロセッサ 動作周波数 : 88 MHz ピーク性能 ( 倍精度 :7 Gflops

11 ( グラフィックボード HOST E5 67.6GHz cpu 8core キャシュ MB.6GHz GFLOPs cpu cpu.8gflops Xeon Phi5P.5GHz 6コア,スレッド/core.5GHz6.88GFLOPs. サーバー (557 8 コア.9GHz キャシュ 8MB/ コア (e5 6 コア.66GHz キャシュ MB/ コア ( 上記グラフィックボードの HOST E5-67 もここに含めました

12 . 基本演算. メモリバンド幅測定 メモリバンドテスト copy do j, N do i, N a (i, j b(i, j end do end do trans do do a (i, end end j j, N i, N do do b( j, i 性能測定結果 (MB/sec N 演算 SR6 E5-67 Phi5P core core 6sp sp 89 copy trans copy trans

13 詳細 trans は参照非連続 N=89 非連続 sr6 sr6 E5-67 Phi5P 要素 core 6sp 6sp sp 連続 なし 非連続 参照 非連続 格納 非連続 両方 N=89 非連続 sr6 sr6 E5-67 Phi5P 要素 core 6sp 6sp sp 連続 なし 非連続 参照 非連続 格納 非連続 両方 N=89 非連続 st st st st st 要素 sp 6sp sp 6sp 8sp 連続 なし 非連続 参照 N=89 非連続 st st st st st 要素 sp 6sp sp 6sp 8sp 連続 なし 非連続 参照

14 演算を含めた場合の結果 (MMAX=89 D Copy D Scale D Add do j =,MMAX do j =,MMAX do j =,MMAX do i =,MMAX do i =,MMAX do i =,MMAX r(i,j=p(i,j p(i,j=scalar*q(i,j q(i,j=p(i,j+r(i,j enddo enddo enddo enddo enddo enddo D Triad do j=,mmax do i=,mmax r(i,j = p(i,j + scalar*q(i,j enddo enddo Transpose do j=,mmax do i=,mmax r(i,j=p(j,i enddo enddo SR6/XM と SR6/M の 機種で測定していますが, 約 倍の差があり, 演算性能の差 6% とは大きな差があります 特に,Transpose 6 スレッドでの差が顕著です

15 メモリバンド幅性能測定結果 (MB/sec SR6/XM MB/sec function スレッド スレッド スレッド 8スレッド 6スレッド スレッド 6スレッド D Copy D Scale D ADD D Triad Transpose SR6/M MB/sec function スレッド スレッド スレッド 8スレッド 6スレッド スレッド 6スレッド D Copy D Scale D ADD D Triad Transpose

16 . rup s 例題 rup's例題とは, 四則演算からなる簡単な 計算で, 有効ビット数がビット以上を要求 する計算のため, 精度 性能の検証に よく使用されるものです a f a a 7767., b a b 5.5b (a より,f.75b b 6 b a 6 b a b 96. (a b b b b 5.5b.75b

17 .. ieee75-8 系 拡張 倍精度 : 拡張倍精度変数を つつなげたもの 6 倍精度 : 倍精度変数を つつなげたもの 8 倍精度 : 倍精度変数を つつなげたもの rup 例題テスト結果 e5 実行結果一覧表 実行回数,, 演算量.9GFLOP ( 各精度 精度 言語 実行時間 sp 実行時間 台数効果 ( 秒 8sp( 秒 拡張 倍精度 icc(strict 拡張 倍精度 icc(etended 倍精度 icc 倍精度 ifort 倍精度 icc 倍精度 ifort tk 実行回数 6,, 演算量.GFLOP ( 各精度 精度 言語 実行時間 sp 実行時間 台数効果 ( 秒 6sp( 秒 拡張 倍精度 icc(strict 拡張 倍精度 icc(etended 倍精度 icc 倍精度 icc

18 拡張 倍精度は etended オプション効果あり ( 両 CPU で Fortran と c++ での有意な性能差はない 拡張 倍精度と 6 倍精度の性能は, e5 は拡張 倍精度 Tk では 6 倍精度が良いという逆の傾向がでています

19 rup 's 例題実行結果一覧,, 回実行 実行時間 ( 秒 C FORTRAN 精度 E5-67 Phi5P E5-67 Phi5P 6sp core 6sp core 拡張倍精度 * 倍精度 倍精度 倍精度 rup 例題実行時間 ( 秒 一覧表 整数演算方式 精度 E5-67 Phi5P 5 倍精度 倍精度 倍精度 倍精度 倍精度 倍精度 倍精度 倍精度 Phi5P の threads/core 実行時間 ( 秒 改良版 threads 倍精度 倍精度 DD 形式 ( 一般サブルーチン使用 精度 E5-67 Phi5P 6sp core 倍精度 倍精度 倍精度 倍精度

20 E5 67とPhi5Pの性能に関しては, DD形式では,8倍精度演算で性能が逆転する またDD形式と整数演算方式では E5 67において倍精度演算で性能が逆転しています Phi5Pでは, 最適化オプションで optとoptsで使用メモリ量が大きく変わるので, 基本的にはoptで実行するのが良い またレジスタの問題で整数演算方式では演算量が多くなると, 最適化が妨げられる事があります

21 E5-66 rup's 例題実行時間一覧表,, 回実行の実行時間 ( 実行時間 : 秒 sp 数 拡張倍精度 * 6 倍精度 6 倍精度 8 倍精度 8 倍精度 c c F c F 拡張倍精度 * 及び 6 倍精度 FORTRAN で sp の場合の性能が良い

22 .. POWER 系 スレッド実行 fa,nofa SR6/M,, 回実行 コンパイラ f9 f9 lf lf 精度 fa nofa fa nofa ( 秒 ( 秒 ( 秒 ( 秒 倍精度 倍精度 倍精度

23 8 倍精度演算方式 QD 方式 : 倍精度変数を つつなげた方式 DQ 方式 : 倍精度変数を つつなげた方式 ( 注 (DQ 方式は QD 方式より精度が悪くなる 最悪の場合 TD( 倍精度変数を つつなげたもの と同じ精度になる (ieee75-8 での DQ 方式の有効ビット数は 6 ビットで,DD 方式の QD 方式より精度が良い Rup's 例題 8 倍精度演算,, 回実行の時間 ( 秒 SR6/M ノード実行 sp 数 QD 方式 DQ 方式

24 .. 超多倍長演算 多倍長 rup's例題一覧表,, 回実行 E5-67 実行時間一覧表 ( 秒 演算精度 6sp 8sp sp sp sp 68 倍精度 倍精度 倍精度 倍精度 倍精度 倍精度 倍精度 倍精度 倍精度 8.87 Phi5P sp/core 演算精度 sp 8sp sp 6sp 68 倍精度 倍精度 倍精度 倍精度 倍精度 倍精度 倍精度 8 倍精度 sp sp 88 倍精度 倍精度.59 ともに並列化効果は出ていますが,E5-67 の効果が大きくなっています

25 . 行列積, ストラッセン行列積 計算機の基本性能を見る演算としては行列積計算が適しています ここでは倍精度演算を SR6/M SR6/XM,BG/Q でその特徴を検証し倍精度以外の行列積計算をサーバー系を含めてその性能を見ています ストラッセン行列積は行列積演算の演算量削減方法と知られていて, 連立一次方程式の求解や後に示すルジャンドル陪関数を用いたいた計算にも適用が試みられています このためその演算量削減効果と実行時間短縮の効果を SR6/M を使用して検証しています

26 .. 行列積計算... SR6/M システム 行列積計算を行う場合, 行列のサイズ N が小さい場合, つの行列すべて合わせてもキャシュに収まる範囲にあります T 行列積計算でC A Bを実行する場合の総メモリ所要量 N 実数型 複素数型 KB 8KB 6 96KB 9KB 96 6KB KB 8 8KB 768KB 56 56KB MB 5 6MB MB 768.5MB 7MB 赤色 LD 内に収まる 緑色 L 内に収まる 紫色 L 内に収まる N> ( 桁 でキャシュチューニングが重要となる

27 キャシュチューニングの際, 内積型計算を使用する場合主要演算時の所要メモリ量 N(N+*8or6 が MB 以下に抑える事が重要になります 主要演算時の所要メモリ量 実数型 複素数型 N ( N 8バイト N ( N 6バイト N 実数型 複素数型 8KB 6KB 6 KB 6KB 6 KB 66KB 8KB 6KB 9 9KB 58KB 56 56KB KB KB KB 8 58KB 6KB

28 SR6 / M ノードで内積型行列積計算で性能が出る条件 ( sid命令が適用される ( 最内側ループ( k 内の演算での乗加算命令が適用される ( 最外側ループ( j のループ長が 8の倍数 ( 中間のループ( i のループ長が6 の倍数 (5 主要演算時の所要メモリ量がMB以下 n 実数型では N=8 が上記すべての条件を満たす 複素数型ではすべてをみたすサイズはないので N に応じて最適なものを探します これらの事から N が 桁でも実数型では N=5, 複素数型では N=6 を境にチューニング手法を変える必要が出てきます

29 実数型行列積計算 実数型行列積計算性能一覧表 内積タイプとは, 以下の重 DO DO j, n, lj DO i, n, li DO k, n ループで lj * liで示しています. ノード st=off,sp= 実行結果一覧 N 内積タイプ GFLOPs * * * * * * * * 定義式どうりのソースでの結果 (lj=,li= N sp 6sp

30 性能向上要件 (5 のテスト N ( N 8バイト MB での実数型行列積計算の性能 ノード st=off,sp 内積タイプ N=8 N=5 * * * *

31 複素数変数行列積計算 複素数型行列積計算性能一覧表 内積タイプとは, 以下の重 DO DO j, n, lj DO i, n, li DO k, n ループで lj * liで示しています. ノード st=off,sp= 実行結果一覧 N 内積形式 GFLOPs * * * * * * * * * * * 定義式どうりのソースでの結果 (lj=,li= N sp 6sp

32 ... システム st= の効果検証 MULTIN N N B A C nosid sid RMULT RMULT CMULT A B A C nosid sid T T T プログラムは中間アンローリング数は最外側回回実行回数はを計算する 実行列積で計算は実数変数でて定義どうり計算するプログラム定義式どうりに計算する複素数変数プログラム回計算する を複素行列積で行いました 検証は複素数変数行列積と実数変数行列積計算で 6,, 8,5, 88,8 (, 8 ( (,

33 実数変数行列積計算 C A T B ultin *6 タイプ sid N=8, 演算量 =5GFLOP 実行時間 ( 秒 一覧表 sp 数 実行時間 GFLOPs N=88 演算量 =78GFLOP 実行時間 ( 秒 一覧表 sp 数 実行時間 GFLOPs N=8 演算量 =GFLOP N=5 演算量 =56GFLOP 実行時間 ( 秒 一覧表 sp n=8 n=8 n=5 n=5 実行時間 GFLOPs 実行時間 GFLOPs 最高性能一覧表 N GFLOPs 実行効率 (% 最高性能は SR6/M と同じ N=8 で達成

34 sp での st=,st= の性能比較 st とst でのsp性能比較 内積タイプ *6 N st= st= 演算量 実行時間 性能 演算量 実行時間 性能 (GFLOP ( 秒 (GFLOPs (GFLOP ( 秒 (GFLOPs 内積タイプ * N st= st= 演算量 実行時間 性能 演算量 実行時間 性能 (GFLOP ( 秒 (GFLOPs (GFLOP ( 秒 (GFLOPs

35 複素変数行列積計算 C ( A T B T A N=89, 演算量 =8GFLOP 実行時間 ( 秒 一覧表 sp 数 rult rult cult cult sid on off on off N=88, 演算量 =7GFLOP 実行時間 ( 秒 一覧表 sp 数 rult rult cult cult sid on off on off N=, 演算量 =.668GFLOP 実行時間 ( 秒 一覧表 sp 数 rult rult cult cult on off on off Sp 数,96 が良い値となっています

36 ... BG/Q システム 複素数変数行列積計算 C T A T B の計算を複素数変数のままと 実数変数に直した場合を実行 SIMD命令を適用させるため, qstrict オプションは外しています 測定条件は ノード openpによる並列化

37 bg 複素変数行列積 C T A T B - O5- qarch qp - qtune qp - qcache auto - qhot level - qipa level - qsid auto - qsp op 複素数型 実行時間一覧表 ( 秒 n 演算量 タイプ 6sp sp 8sp 6sp (GFLOP 7 * * * * * * * ( 注 8sp では,96,9 はタイプ * にしている 最高性能は 8GFLOPs

38 bg 複素変数行列積 C T A T B - 複素数型 O- qarch qp - qtune qp - qcache - qipa level - qsid auto ー O - qsp op auto - qhot level 複素数型 実行時間一覧表 ( 秒 n 演算量 タイプ 6sp sp 8sp 6sp (GFLOP * * * * * * ( 注 8sp では,96,9 はタイプ * にしている 性能は 8GFLOPs,8GFLOPs と -O5 の場合より少し良くなっています ( 誤差の範囲

39 bg 複素変数行列積 C T A T B - O5- qarch - qhot level - qipa level - qsid auto - qsp op qp - qtune qp - qcache auto 実数型 実行時間一覧表 ( 秒 n 演算量 タイプ 6sp sp 8sp 6sp (GFLOP 7 * * * * * * * ( 注 8sp では,96,9 はタイプ * にしている ソースは SR6/M,XM と同じものを使用 性能は 6GFLOPs,6GFLOPs,68GFLOPs と複素変数の場合に比べて有意に劣っています

40 ... 各種演算精度の行列積計算 比較は SR6 ノード (core,6sp,es-67 6sp, Phi5P core で実施しています 最も多くのケースを実施したのは行列サイズ N=,N=8 の 回実行です C は,ES-67,Phi5P で倍精度, 拡張倍精度, 拡張倍精度 + 拡張倍精度のケースを行っています 名称はそれぞれ Double,edouble,ddedouble としています FORTRAN は SR6,ES-67,Phi5P で倍精度 ultd, 倍精度 (ieee75-8 データ形式 ultq,dd 形式 倍精度 ddultd,dd 形式の 6 倍精度,8 倍精度, 倍精度, それぞれ ult6,ult8,ult です.SR6 では特に ultd,ultq は実施していません また比較という事で, テストプログラムは定義式どうりのコーデイングでアンローリング, キャシュチューニング等は実施していません 演算量は以下の様になっています 演算量 (GFLOP 一覧 プログラム N= N=8 ultd 8 ddultd 7 76 ult ult ult

41 行列サイズ N= の結果は以下の様になっています 行列積 N=, 回実行の性能比較 実行時間 ( 秒 プログラム SR6 phi5p E5-67 core 6sp sp 6sp ultd.8. ddultd ultq double.9. edouble ddedouble ult ult ult (dd 形式 倍精度は phi5p は高速 (dd 形式 6 倍精度,8 倍精度, 倍精度は SR6 core と phi5p の性能が同等 (C とソフトウエアサポートルーチンは E5-67 が高速

42 行列サイズ N=8 の結果は以下の様になっています 行列積 N=8, 回実行の性能比較 実行時間 ( 秒 プログラム SR6 phi5p E5-67 core 6sp sp 6sp ultd.5. ddultd ultq double edouble.6.76 ddedouble ult ult ult 性能比較の傾向は N= の場合と全く同じ

43 各種行列積計算の性能は以下の様になっています 各種行列積計算性能比較表 実行時間一覧表 ( 秒 N CPU double edouble ddedouble ultd ddultd E Phi5P E Phi5P E Phi5P E Phi5P E Phi5P ( 倍精度では c,fortran とも N=7 で Phi5 P の性能が E5-67 を上回ります ( 拡張倍精度はすべて E5-67 の性能が上回っています (SR6 の 倍精度では,Phi5P が非常に効力を発揮する事がわかります

44 これに関連した計算は以下の様になっています 倍精度行列積計算 N=56 回の実行時間 SR6 single ジョブ 9.68 sec ieeequad sec 557 double+double sec etend double +etend double 7.6 sec e5 gcc double+double 65.9 sec etend double +etend double sec e5 icc double+double 9.57 sec etend double +etend double 6.66 sec Phi5P sp が SR6 ノードより高速という結果と合わせ Phi5P が有効な事がわかります

45 6 倍精度,8 倍精度, 倍精度行列積計算 SR6/M ノード実行コンパイルオプションは -W, opt(o(s,disbracket( opt(o(ss にしたり disbracket( を外すと結果不正となります 行列積テスト結果一覧表 N=, 回実行時間 ( 秒 精度 single core 6sp 6 倍精度 倍精度 倍精度 ( 並列化効果はよくでている (8 倍精度は SIMD がよく効いているので他の つの精度に比べ演算数比率より性能が良い (8 倍精度は つの倍精度変数を扱うため

46 に比例すると言えます の行列積の演算量は数学的にはとなりの行列積の演算量は通常の行列積の 演算量比演算量通常の行列積段目はもとの行列のサイズ行列和の回数段目の行列積の回数 N,strassen, 8 7 ( strassen 7 8 ( (N 7 N 8 N (7 N 8 7 ( N N (7 N 8 7 ( strassen N N N, (7 A, 7 M A 8M A, 7M M A, M A, M P P P P C P P C P P C P P P P C B (B A (A P B (B A (A P B A (A P B (B A P B (B A P B A (A P B (B A (A P C C C C,C B B B B B, A A A A A AB C.. ストラッセンの行列積代表的な例でその計算方法と演算量を示しました

47 演算量に関して strassenの行列積性能一覧表 SR6/M ー Oss N=6556 st on 6sp/core プログラムの性能 (GFLOPs 演算量 N (FLOP / 実行時間 ( 秒 段数 プログラムモニター 実行時間 (GFLOPs (GFLOPs ( 秒

48 見かけ上の性能値と論理最大性能 ストラッセン行列積性能一覧表 SR6/M st=off 測定値 (GFLOPs 演算量 N (FLOP / 実行時間 ( 秒 N 段数 sp 測定値 論理最大性能 (GFLOPs (GFLOPs

49 . 実アプリケーション. 重力多体問題 (N 体問題 重力多体問題は F i G r j ij r G : 万有引力定数, i : 粒子 iの質量 r : 粒子 iの位置, F : 粒子 iにかかる力 i i ji で計算します ij i r ij r i r j 演算精度は倍精度でSIMD命令が適用されます 演算量は粒子数 Nの二乗に比例します

50 .. SR6/M システム ノードでフラット MPI で SMT ON SMT OFF の場合の結果です 演算量が N=**n の二乗に比例しますので 実行時間は log スケールで表しています タイムステップ数は です Log( 実行時間 ( 秒 st on st off 粒子数 N=**n n SMT ON,OFF での実行時間に有意差は見られません

51 また フラット MPI とハイブリッドの差を N=**8,** で調べた結果は以下の様になりました サイズ**n sp 数 task 数 ノード数 実行時間 ( 秒 サイズ**n sp 数 task 数 ノード数 実行時間 ( 秒 フラット MPI とハイブリッドでも実行時間に有意差は見られません

52 .. BG/Q システム 領域分割をせずに並列実行する場合 粒子数 N=**5=768 から N=**=8576 とすると スレッド当たりの所要メモリ量が.8MB から MB ため SMP と MPI のハイブリッド方式がフラット MPI より効果がでます 演算量が N の二乗に比例するため 実行時間の軸は log を取っています 実行結果は以下の様になりました タイムステップ数は です Log( 実行時間 ( 秒 フラット MPI ハイブリッド n 粒子数 N=**n

53 .. 各種計算機性能比較 使用した計算機は以下のものです 計算機 A BG/Q 8node sp 計算機 B BG/Q 8node 89MPI A C BG/Q node 8node 8MPI sp D SR6/M node 8MPI B E SR6/M BG/Q 8node node 89MPI 56MPI C F HD587 BG/Q node core 8MPI G HD697 6core D SR6/M node 8MPI E HD697 SR6/M node 56MPI 88MHz 56core 68GFLOPs 76GB/sec F HD587 85MHz HD587 6core core 5GFLOPs 5.6GB/sec G HD697 6core HD697 88MHz 56core 68GFLOPs 76GB/sec HD587 85MHz 6core 5GFLOPs 5.6GB/sec 扱ったケースは, 演算量は可変で粒子数を変化させた場合と演算量は一定で粒子数を変化させた場合を扱っています

54 ... 演算量可変ケース 重力多体問題実行時間 log( 実行時間 ( 秒 6 実行性能 BG/Q ノード SR6/M HD697 ノード 8 6 A B C D E F G n 粒子数 N=**n

55 5... 演算量一定ケース 粒子数 N,,で タイムステップ数はそれぞれ,, 65,として演算量を同じにして ( 平方根計算あり7GFLOP, なし58GFLOP 性能を比較. 演算は倍精度演算 N の場合, 並列化オーバーヘッドにより SR6の性能が勝っているが, N,ではGPUの性能が勝っています また ieee75-8 系の特徴でもある拡張倍精度演算に関する結果も扱っています

56 N 体問題実行結果一覧表 ( 秒 粒子数 =, タイムステップ = 演算量 ( ソース =5GFLOP SR6 ノード core 6sp 最適化 OPT SOPT E5-67 6sp 6.9 各種加速装置実行時間一覧表 GPU ボード ボード ボード ボード (sp/core (sp/core (sp/core (sp/core HD HD HD W8.6 HD Phi5P

57 N 体問題実行結果一覧表 ( 秒 粒子数 =, タイムステップ =65 演算量 ( ソース =5GFLOP SR6 ノード core 6sp 最適化 OPT SOPT E5-67 6sp.69 各種加速装置実行時間一覧表 GPU ボード ボード ボード ボード (sp/core (sp/core (sp/core (sp/core HD HD HD W8 8.8 HD Phi5P

58 N 体問題実行結果一覧表 ( 秒 粒子数 N=, タイムステップ = 演算量 ( ソース =5GFLOP SR6 ノード core 6sp 最適化 OPT SOPT E5-67 6sp.55 各種加速装置実行時間一覧表 GPU ボード ボード ボード ボード (sp/core (sp/core (sp/core (sp/core HD HD HD W8.796 HD Phi5P

59 拡張倍精度演算 N 体問題倍拡張倍精度実行時間一覧表 ( 秒 N 反復回数 E5-67 Phi5P E5-66 E5-66 6sp sp 6sp sp 倍精度の場合とは逆に E5-67 の性能が E5-66 の性能を上回っています 拡張倍精度に対する処理で E5-67 のコンパイラが e5 のコンパイラより非常に適している事によります

60 . 分子動力学計算 分子動力学計算は以下のつの計算からなり クーロン力計算は重力多体問題と同様の計算となりますが 異なる計算として ワンデルワールス力計算があります 逆二乗則と逆六乗則の差により 演算する範囲が異なり 演算量も単純に粒子数の何乗に比例するとは言えません 演算精度は倍精度で SIMD命令は一部適用されます クーロン力 q i F : 電荷量, i q i ji q r ij r : 真空の誘電率 j ij r ij r i r j ワンデルワールス力分子間に働く分散力で 等方向性で原子間距離の6 乗に反比例する力 a F k 6 r b : 分極率

61 .. SR6/M システム詳細 SR6/M 複数ノードでの実行時間は以下の様になっています 実行はフラット MPI を使用しています 所要メモリはそれぞれ.8MB,8.MB,5.7MB となっています 実行時間 ( 単位秒 N=8 MPI 数 ノード数 VdW Coulob 全体 N=6 MPI 数 ノード数 VdW Coulob 全体 N=8 MPI 数 ノード数 VdW Coulob 全体 粒子の数は n=n** で演算量は Coulb( クーロン力 は n(n-/ に比例,VdW( ファンデルワールス力 は n(n-// 近傍にある粒子の数に比例します VdW ではすべてのケースで SMT OFF の効果があり クーロン力は N=8,node では SMT ON がそれ以外は SMT OFF が効果があります これは転送量が MPI 数に比例するため 演算量との比で SMT ON,OFF の効果が異なる様になります

62 .. 各種計算機性能比較 N 所要メモリ (MB 実行時間 ( 秒 一覧表 N=8 SR6/M MPI 数 ノード数 VdW Coulob 全体 BG/Q MPI 数 ノード数 VdW Coulob 全体 N=6 SR6/M MPI 数 ノード数 VdW Coulob 全体 BG/Q MPI 数 ノード数 VdW Coulob 全体 N=8 SR6/M MPI 数 ノード数 VdW Coulob 全体 BG/Q MPI 数 ノード数 VdW Coulob 全体

63 メモリの制限により,N=,N=8 では BG/Q 8MPI は実行不可 N=96 SR6/M MPI 数 ノード数 VdW Coulob 全体 BG/Q MPI 数 ノード数 VdW Coulob 全体 N= SR6/M MPI 数 ノード数 VdW Coulob 全体 BG/Q MPI 数 ノード数 VdW Coulob 全体 N=8 SR6/M MPI 数 ノード数 VdW Coulob 全体 BG/Q MPI 数 ノード数 VdW Coulob 全体

64 BG/Q dcore sp 数 ノード数 VdW Coulob 全体 対 st= dcore sp 数 ノード数 VdW Coulob 全体 対 st= dcore sp 数 ノード数 VdW Coulob 全体 対 st= dcore sp 数 ノード数 VdW Coulob 全体 対 st= MPI オーバーヘッド削減の効果は BG/Q ででています X では ノード並列のためほとんど変化はありません

65 st= に対する効果一覧表 プログラム cpu st= st= st= dcore BG/Q dcore BG/Q dcore ocore 論理コア 8(new,6(old の比較 dcore sp 数 ノード数 VdW Coulob 全体 new old new old new dcore sp 数 ノード数 VdW Coulob 全体 new old new old new BG/Q MPIオーバーヘッド削減の効果が大きい (st= X 論理コア8(st= の拡大の効果が大きい

66 . 大規模疎行列の反復計算.. 対称問題倍精度演算 ポアソン方程式 u f (,u ( [,] [,] [,] N 収束判定値 : 共役残差 u f sin( sin( y sin( z sin( sin( y sin( z ( 注 初期値は(i cgs, cgsilu, cgsilは収束しない ので(i. : (, の一様乱数値. cgsはcgsiluを,cgsはcgsilを並列化して初期値 (i.としたもの

67 対称行列反復法実行時間一覧表 ( 秒 プログラム 557 E5-67 Phi5P E5-66 E5-66 SR6 SR6 sp 6sp sp 6sp sp sp 6sp dsor dsor dadi dbcg dcg dcgs dscg dbicgs dcgsilu dcgsil dgpbicg dcgs dcgs E5-67,E5-66,Phi5P は SR6 と比較して実行効率ではほぼ同等の性能がでています

68 .. 非対称問題 y,z u(, 9 nz ny, n R square gradient cgs(conjugate 9 nz ny, n R gradient e bcg(biconjugat z sin( y sin( sin( e y,z u(, [,] [,] [,] f u R u yz 初期値収束判定値は共役残差法法解法非対称問題 解析解領域問題は次元ポアソン方程式演算は 倍精度演算

69 非対称問題実行時間一覧表 ( 秒 プログラム E5-67 Phi5P E5-66 E5-66 SR6 SR6 6sp sp 6sp sp sp 6sp bcg cgs 倍精度演算では SR6 の性能がカタログ性能比を考慮しても非常に良い性能を示しています 今回の cgs 法は前処理付きだが以前の方式ではサイズが ** で 6 倍精度 SR 6 6sp 9.5 秒 E5-67 6sp 7.6 秒 8 倍精度 SR 6 6sp.55 秒 E5-67 6sp 7.8 秒 と倍精度変数をつなげた形式では E5-67 の実行効率が勝っていました

70 . QCD 計算 テストケース case data 6666 case data case data 測定条件 BG / Q, SR6 / M, すべてノード st数 BG / Q, st,, SR6 / M st, スレッド並列のソース BG/Q ppgenで作成したソース -qsp=auto でコンパイル. SR6/M, ともに自動並列でコンパイル

71 QCD 性能測定結果一覧表 ( 対 スレッドに対する台数効果 data BG/Q SR6/M st= st= st= st= st= 最高性能 sp 数 case case case 最高性能 BG/Q st= =>sp SR6/M st= =>sp =>sp 扱うデータがおおきくなると SR6/M は大きな性能低下が見られます BG/Q はデータサイズによる性能低下は見られません BG/Q st= (6sp で台数効果 と高い

72 性能比較 QCD ***6 並列実行性能測定結果 実行時間一覧表 ( 秒 ノード sp BG/Q ソースは Parallel Progra Generetor で作成 SR6 9.5 BG/Q.9 SR6/BG 実行性能比 6.55 SR6/BG カタログ性能比.79 pi,sp 並列 SR6 BG/Q node total node total SR6 ( 対 BG/Q node 実行性能比 カタログ性能比 ソースチューニングがない場合,SR6/M が BG/Q に対してカタログ性能比以上に実行性能比がでています

73 .5 量子モンテカルロ法による物性スペクトル計算 プログラムは高エネルギー加速器研究機構, 物構研岩野氏より提供いただいたのをベースにしています 数値計算からみた場合の計算内容 G,A A A B C G i パラメータ : β L,dt L C の絶対値の最大値, 最小値 e i i Tr(G C,B,C I ( 単位行列 L C C : n n 行列...C...C A (I B L i n,l, β,u dtu を計算し, 結果のチェックは Tr(G,e L dtu で行います 最初の実行結果 ( 倍精度 n=,β=,l=,u=5 n=,β=,l=8,u=5 結果は正確結果は不正確 目標 n=,l=8 を固定 β=, u=5,6,7,8,9, での結果を 進 桁一致させる

74 使用されている計算手法 行列の QDR 分解 行列 V を V=QR (Q: ユニタリ行列,R: 上三角行列 に分解後,R の対角要素を D に入れ,R の対角要素を にする 扱ったケース : 設定したパラメータより結果の精度が見積もりやすい P e dtul,e( Tr(G プログラム実行中に現れる数値の最大値 P 絶対値最小値 /P 行列やベクトル要素に現れる数値の最大値 P, 以外の 絶対値最小値 /P E(,E(L, L L でE( 演算に必要な最小ビット数 log P, 以外の で最小になるので, P logp 精度検証結果 ( 進 桁一致したもの 6 倍精度 β=,u=5,6,7,8,9 8 倍精度 β=,u=5,6,7,8,9, β=,u=5,6,7,8 倍精度 β=,u=9,u=

75 .5. SR6/M システム 性能測定条件 倍精度 β=,u=5 結果は 進 桁までしか一致しないが他との比較のため. 6 倍精度 β=,u=6 8 倍精度 β=,u=8 倍精度 β=,u= チューニング手順 (6 倍精度,8 倍精度での加減算, 乗算, 除算, 平方根計算をサブルーチン化して, uinline オプションを適用 ( 演算量の多い行列積計算,QDR 積計算部分をサブルーチン化し 並列実行する ( 加減算, 乗算において,if 文削除のアルゴリズムに変更 ((~( の方式を取り入れて 倍精度演算ルーチンを作成

76 性能測定結果 量子モンテカルロ実行時間 ( 秒 SR6/M ノード 倍精度 倍精度 st on,off 実行時間 QDR 積計算 6sp core core 9.7 並列化不可 sp. 分割並列化 6 56 数値制限解除 5 5 uinlineオプション 精度 サブルーチン化 6sp core 6 倍精度 倍精度 行列積,QDR 積計算 IF 文削除 精度 6sp core 6sp core 6 倍精度 倍精度 倍精度演算の問題 加減算, 乗算での演算数の増加 ( 最適化,uinline の効果が減少 ( 行列積計算での並列化効果の減少 (QDR 積計算が並列化不可となり, つのサブルーチンに分割した事による並列化効果の減少表現できる数値範囲の制限大きな数の平方根計算では アルゴリズムの変更が必要

77 .5. ieee75-8 形式 E5-67 6sp Phi5P core モンテカルロシュミレーション E5-67 Phi5P 精度 β U L 実行時間 実行時間 ( 秒 ( 秒 6 倍精度 倍精度 倍精度 ieee8 倍精度 DQ(8 倍精度 (β=,u= を計算する場合 dd 形式では 倍精度演算が必要で ieee 形式では 8 倍精度で事足ります (SR6 で scope オプションで実行すると 倍精度は 69 秒. このため, 計算機種の差ではなく, コンパイラの最適化能力の差が効いていると言えます

78 DD 形式 8 倍精度と ieee75-8 形式の 8 倍精度演算では, 正しく計算出来るパラメータ領域に差がある 一概にどちらが良いとは言えない β=,u=8 では DD 形式 8 倍精度が良い. β=,u= では DD 形式では 倍精度を使用する必要があり,ieee75-8 形式の 8 倍精度が良い 倍精度演算では ieee75-8 形式の方が性能が良いと考えられ,GPU 等のアクセラレータで 8 倍精度演算までの DD 形式を適用している理由とも合う結果となっている

79 .6 ファインマンループ積分の多項式実定数ループの多重度ループ内の素粒子の数時空次元 N nl N n ( L N N N i i nl N,...,,,C,D : :,L : N :, (n n : C i (D C... ( d nl (N ( ( I( ファインマンループ積分の一般式は以下の様になる SMP Phi5P 6SMP 67 E5 SMP HD587.5, 5 s y ( y ( ys D dyd D I vt ra inf.6. 5 実行条件で実施しています 倍精度演算で計算

80 GFLOPs. 次元積分 (VTX 性能測定結果 ( 演算量はソース上から算出 HD587 E5-67 Phi5P. 5.. N サイズ

81 .6. Infra bo 計算 ( ( ( ( ( f e y y z y z y y z y tz sy D dzdyd D I 計算式 (, (.5 5,, 5, 5 倍精度倍精度 e f t s テストデータ

82 倍精度 N=89, 演算量 7GFLOP SR6/M ノード 自動並列 BG/Q ノード ノード間 MPI ノード内 スレッド並列 計算機 sp 数 st 数 性能 (GFLOPs 実行効率 (% SR6/M 68 7 SR6/M BG/Q 6 86 BG/Q 9 BG/Q HD587 6 HD 実行効率は 5%-%

83 GFLOPs 6. 次元積分 (BOX 性能測定結果 ( 演算量はソース上から算出 5... HD587 E5-67 Phi5P サイズ N

84 倍精度 サイズ N=8 演算量 68GFLOP テーブルサイズ KB SR6/M,8ノード フラットMPI BG/Q ノード フラットMPI 倍精度性能測定結果 (GFLOPs 計算機 ノード数 st 数 性能 実行効率 (% SR6/M 8 6. SR6/M 8.7 SR6/M SR6/M BG/Q 5.9 BG/Q BG/Q 75.8 HD587 ボード 8 6 HD697 ボード SR6/M,BG/Q では SIMD が適用されないので,GPGPU 系の計算機の実行効率が良い

85 ノード性能比較 実行時間 ( 秒 サイズN= BG/Q( 換算 6sp SR6/M 6sp BG/Q BG/Q( 換算 SR6/M 5 SMT 数 BG/Q SR6/M SMT= 6 sp sp = sp 6 sp = 8 sp 96 sp = 6 sp 8 sp BG/Q( 換算.8GFLOPs => 98.8GFLOPs にしたときの値

86 実行時間 ( 秒 サイズ N=8 5 5 BG/Q( 換算 6sp SR6/M 6sp BG/Q BG/Q( 換算 SR6/M 5 SMT 数 BG/Q SR6/M SMT= 6 sp sp = sp 6 sp = 8 sp 96 sp = 6 sp 8 sp BG/Q( 換算.8GFLOPs => 98.8GFLOPs にしたときの値

87 サイズ N=96 実行時間 ( 秒 5 5 BG/Q( 換算 6sp SR6/M 6sp BG/Q BG/Q( 換算 SR6/M 5 SMT 数 BG/Q SR6/M SMT= 6 sp sp = sp 6 sp = 8 sp 96 sp = 6 sp 8 sp BG/Q( 換算.8GFLOPs => 98.8GFLOPs にしたときの値

88 BG/Q node MPI SR 6PE GPGPU XM 6SMP XM MPI*6SMP TK 6MPI TK 8MPI TK 56MPI BG/L node MPI 倍精度 INFRA BOX N= 各種計算機実行時間 ( 実行時間 ( 秒 論理最大性能 BG/Q GFLOPs SR 57.6GFLOPs GPGPU 88GFLOPs XM 8.8GFLOPs TK 6MPI 588GFLOPs TK 8MPI 76GFLOPs TK 56MPI 5GFLOPs BG/L 57.GFLOPs 実行時間 ( 秒

89 倍精度 INFRA BOX N= 各種計算機実行時間 (..8.6 実行時間 ( 秒 論理最大性能 BG/Q GFLOPs SR6/M node 9.9GFLOPs 8node 78.8GFLOPs.. 実行時間 ( 秒

90 BG/Q node 5MPI BG/Q node MPI BG/Q node 8MPI BG/L node 8MPI GPGPU SR6/M node 8MPI SR6/M node 56MPI SR6/M 8node 56MPI SR6/M 8node 5MPI 倍精度 INFRA BOX N=8 各種計算機実行時間 実行時間 ( 秒 論理最大性能 BG/Q GFLOPs SR6/M node 9.9GFLOPs 8node 78.8GFLOPs GPGPU 88GFLOPs BG/L 57. GFLOPs 実行時間 ( 秒

91 倍精度 INFRA BOX N=8 と N=96 各種計算機実行時間比率 演算量比率 (N=96/N=8=8 8 6 実行時間 (N=96/N=8 比率 倍率

92 .6. assless 計算を計算する必要があります 倍精度演算で精度良く計算するにはとなり n n... n n n... n n n... n n...d d d... (... n...(n n ( n (,,...d d d... (... n...(n n ( n (...d d d... (... n n n

93 倍精度演算による実行結果 条件一覧表 n h N.5^6= ^6= ^= ^= 実行時間一覧表 ( 秒 st=on,6sp #NAME? n SR n Phi5P sp E5-67 6sp

94 Phi5Pでのチューニング例 n... n ( ( ( n ( (... ( n( n...(n n n n n において., n 7,n 8の場合を計算した 倍精度で計算するため, 右辺の式を使用 n 7,h.6, N n 8,h., N で実行した結果は下記のとうり n n d n d d n n d......d n ( n...d n d n d n...d n=7,n=8 assless 計算実行時間 ( 秒 一覧表 SR6 6sp,-Os 自動並列 XM 6sp,-Os 自動並列 E5-67 6sp -Os -openp Phi5P sp -Os -openp n SR6 XM E5-67 Phi5P n=7,n=8 の場合の Phi5P の値が悪すぎる

95 n 7 実行時間 78秒 68秒 n 8 実行時間 79秒 6秒 n 7だとSR6とほぼ同等,n SR6の.倍の性能 8だと 結果 n 7 n 8 解析解 実測結果 解析解 実測結果 N の場合の 実測結果 結果は倍精度演算でも十分な精度となっている n=7 と n=8 の差は h=.6(n=,h=.(n= の違いによる 一般に DE では倍精度演算では N=-5 で十分な精度が得られる

96 このため,y を以下の様に一次元化 do i=,n =(i by=.d- Do i=,n yy=(i*by bz=.d--yy do i=,n*n =list(i yy=ylist(i bz=.d--yy gw(i*by*gw(i 前処理 l= do ip=,n =(ip b=.d- do iq=,n =(iq*b l=l+ list(l= ylist(i= wlist(l=gw(ip*b*gw(iq end do end do wlist(i

97 次元積分計算.6.. S 計算 dudzdyd DC s S y z y 5,,,, ; ( C M se D u z y u z y M y u z zu y u y z u z y E u z y u z y u z y C 5 ( ( ( ( ( ( ( ( ( ( :.887 :,,,, ( 解析近似解 s サイズ N=8, 演算量 =86GFLOP テーブルサイズ KB SR6/M ノード BG/Q ノード計算機 SR6/M SR6/M BG/Q BG/Q 並列化 st= st= flat sp 実行時間 ( 秒 性能 (GFLOPs

98 S (5;,,,,. 5 N GPU HD sec HD sec HD587 HD697 HOST sec sec

99 S ( ;,,,, の性能測定結果 サイズ N=576, 演算量 76GFLOP ( 性能モニター X557,SR 自動並列,E5 MPI スレッド数 X557 E5 SR 実行時間 性能 実行時間 性能 実行時間 性能 秒 GFLOPs 秒 GFLOPs 秒 GFLOPs MPI MPI MPI BG/L MPI スレッド数 N 演算量 実行時間 性能 GFLOP 秒 GFLOPs S (;,,, の性能測定結果 算法反復回数 5回 CPU N 演算量 スレッド数実行時間 性能 GFLOP 秒 GFLOPs BG/L X E

100 s (5;,,,,. 5 N 演算量 76GFLOP 実行時間 ( 秒 一覧表 ボード数 HD587 HD697 gpu gpu HD587 GFLOPs,HD697 6GFLOPs SR6/M システムに匹敵する性能がでています HD587 HD697 HD587 HD697 GPU.55 sec sec HOST sec sec GPU が HOST に比べて非常に良い性能がでています

101 HOST コンピュータでの性能が低い事の原因追究 (GPU 用に 重 DO ループを 重 DO ループにしたのが原因?( 次元積分, 次元積分 ( 使用コンパイラの問題では? 実行時間一覧表 ( 秒 cpu コンパイラ 次元積分 次元積分 557 gcc ifort cpu コンパイラ 次元積分 次元積分 e5 icc e5 icc parallel sp ソースで二次元化したためキャシュミスが発生 コンパイラの最適化処理能力 GPU でのチューニングは 容量の大きいベクトルレジスタを持ったベクトル機でのチューニングと同じ傾向を示します

102 スーパーコンピュータでの 次元積分と 次元積分の比較 実行時間一覧表 ( 秒 6MPI CPU 次元積分 次元積分 SR6/M TK 次元積分の方が性能が良いという傾向はサーバー系 557,e5 と同じ

103 f9 と lf_r コンパイラの比較 S ( ;,,,, 演算量 755 N 5 GFLOP SR6/M, ノード cpu sp 数 f9 lf_r SR6/M SR6/M 7.87 f9 実行効率 7.%, SR6/M.9% と非常に良い値がでています lf_r は f9 以上に実行効率.% とよく, キャシュの有効利用と SIMD が非常によく効いています

104 SR6,E5-67,Phi5P の比較の結果は以下の様になっています 次元積分実行時間一覧表 ( 秒 case 精度 SR6 E5-67 Phi5P 演算量 core 6sp ((GFLOP case 倍精度 case 拡張倍精度 case 倍精度 case 拡張倍精度 case 倍精度 ( 多重 DO ループの一重化は E5-67, Phi5P は若干速くなりますが SR6 は ~ 倍程度遅くなっています (ieee75-8 形式の 倍精度は非常に遅くなっています

105 サイズを大きくした場合と HD587 との比較は以下の様になっています 次元積分計算詳細テスト結果 (s case N= 実行時間 ( 秒 精度 E5-67 Phi5P 6sp core 倍精度 拡張倍精度 case 実行時間 ( 秒 精度 HD587 Phi5P 面 面 9core core 倍精度 拡張倍精度 倍精度 Phi5P は倍精度は高速で E5-67 は拡張倍精度が高速 HD587 の 倍精度は dd 形式のため高速です

106 S テスト S (, N 576,,,, 演算量 GFLOP 倍精度 568 倍精度 87 SR6 実行時間 ( 秒 精度 core 6sp 倍精度 倍精度 E5-67 6SMP 精度 実行時間 ( 秒 倍精度 拡張倍精度 倍精度 8.76 HD587 実行時間 ( 秒 精度 面 面 倍精度 倍精度 Phi5P 実行時間 ( 秒 精度 sp /core sp/core sp/core sp/core 倍精度 拡張倍精度 倍精度 HD697 実行時間 ( 秒 精度 一面 二面 三面 四面 倍精度 倍精度 グラフィックボード性能比較表 精度 HD797 W8 HD798 倍精度 倍精度

107 次元積分計算 laportad I 6 5 D d 5 d d d d D = &-***-***-***-***6-***-** &-.d*** &-**5-**6-***-.d***-**5- **6 &-*** &-**5-.d***6-*5*6-*6**- ***- ***5 &-** &-**5-***-.d***5-**6-*5**- *5*6 &-** * &-***5-***-.d***5-**6-*5**- *5*6 &-***5 &-***6-*5**-.d**5*6-*6**-5***6-5*6** 解析近似解 =

108 laportaf I 6 5 D d 5 d d d d D=.-**-**-**5-**6-**-**.-**.-.d***5-.d***6-**-**5-**6.-**-.d***5.-.d***6-*5*5-*5*6- *6*6-**.-**5-**6.-**-**-**5-.d***6-**5.-**6-*5*5-*5*6-*6*6-**-**5.-**6-**-.d***5-.d***6- *5*5-*5*6.-*6*6-**5-**6-*5*5-*5*6-*6*6 解析解は算出できなったため, つの異なる積分方法で結果を確認

109 Laporta (d SR6/M, ノード サイズ N=8 演算量 668GFLOP 実行時間 ( 秒 一覧表 sp 数 sr6/m SR6/M f9 lf_r f9 lf_r f9,sr6/m の実行効率 6.% と突出して良い値となっています

110 laporta (f サイズ N=8 演算量 6GFLOP ソースは二次元積分 スーパーコンピュータと GPU 実行時間 ( 秒 CPU ノード数 sp 数 pi 数 実行時間 性能 実行効率 ( 秒 (GFLOPs (% BG/Q SR HD HD どの機種も高い実行効率をしめしています GPU と HOST 実行時間 ( 秒 cpu gpu host host/gpu HD HD

111 各種計算機詳細 GPU 実行時間 ( 秒 一覧表 ボード数 HD587 HD697 gpu gpu BG/Q case -O5 -qhot=level= case -O -qstrict -qhot=level= ノード実行時間 ( 秒 case flat hybrid GPU: ボード数の効果がよくでています BG/Q: コンパイルオプションの効果が大きい. オプションにより flat MPI とハイブリッド MPI の効果が変わっています

112 SR6/M ノード実行時間 ( 秒 st sp 実行時間 性能 実行効率 ( 秒 (GFLOPs (% on off.67.8 TK 実行時間 ( 秒 一覧表 MPI 数 実行時間 性能 実行効率 ( 秒 (GFLOPs (% XM 実行時間 ( 秒 一覧表論理コア数 6 f9 st 6 on off lf_r eecution 実行時間 ( 秒 一覧表論理コア数 8 sp 6 8 f lf_r SR6/M,TK ともに実行効率は高くなっています X: 論理コア数増加の効果は f9 に大きくでています.

113 演算量 (N= case GFLOP laporta d 5 laporta f 68 実行時間一覧表 ( 秒 CPU SR6 E5-67 Phi5P HD587 プログラム core 6sp c fortran c fortran 面 laporta d laporta f 演算量が約 5GFLOP 程度だと SR6 と Phi5P,HD587 の性能はほぼ等しい. これは 次元積分でも同様でした

114 次元積分計算 laportag I C D d 7 d 6 d 5 d d d C=*+*5+*6+*+*5+*6+*+*5+*6+ *5 & + *6+*7+5*7+6*7 D = &- (**+**+**+7**+*+*+*7+*+*7+ *7 &*(+5+6 &-***(+++5+X6+7 &-(5**+6**+5*6*(++++7 &-.**(*5+6*7 &-.*((++**6+(++7**5 解析近似解 =

115 laportah C=(+++*( * cc=*+*+*+*+5*5+6*6 +7*7 D=-c*cc.+s*(**( *6*(++++* *6+**5.+t***7.+p*(**( **5.+p*(**( **6.+p*(5*7*(++++**7.+p*(6*7*(++++** d d d d d d D C I

116 テストデータ =.d =.d =.d =.d 5=.d 6=.d 7=.d p=.d p=.d p=.d p=.d s =.d t =.d 解析近似解 =.67989

117 aportai C=(++++5*( (++** cc=*+*+*+*+5*5+6*6 +7*7 D=-c*cc.+s*(**( *5*6+**7- **6.+t**(-*6+5*7.+p*(**( **(6+7.+p*(**( *6*(+5.+p*(*5*( *6*(++*5 *7.+p*(6*7*(++++5+*6*(++*5 * d d d d d d D C I

118 テストデータ =.d =.d =.d =.d 5=.d 6=.d 7=.d p=.d p=.d p=.d p=.d s =.d t =.d 解析近似解 =

119 6 次元積分実行結果一覧表 演算量 (N= case GFLOP laporta g 6 laporat h 8 laporta i 896 実行時間一覧表 ( 秒 cpu SR6 E5-67 Phi5P HD587 プログラム core 6sp c fortran c fortran 二面 laporta g laporat h laporta i 演算量がこれだけ大きくなると,Phi5P, HD587 の性能が SR6 を大きく上回ります また,Phi5P,HD587 のカタログ性能はほぼ同じですが HD587 の性能が Phi5P の倍近く良くなっています これは Phi5P は 6CPU, core という構成によるものと考えられます

120 laporta 計算結果一覧表 演算量 (N= プログラム GFLOP LAPORTA D 5 LAPORTA F 68 LAPORTA G 6 LAPORTA H 8 LAPORTA I 896 実行時間一覧表 ( 秒 SR6 プログラム core 6sp LAPORTA D LAPORTA F LAPORTA G LAPORTA H LAPORTA I E5-67 HOST 6sp プログラム c fortran LAPORTA D LAPORTA F LAPORTA G LAPORTA H LAPORTA I Phi sp(6core プログラム c fortran LAPORTA D LAPORTA F LAPORTA G LAPORTA H LAPORTA I

121 laporta Phi5P 実行時間 ( 秒 一覧表 プログラム sp /core sp /core sp /core sp /core laporta d laporta f laporta g laporta h laporta i グラフィックボード性能比較表 実行時間 ( 秒 プログラム HD797 W8 HD798 laporta d laporta f laporta g laporta h laporta i

122 SR6/M 複数ノード BG/Q LAPORTA (D LAPORTA(D MPI 数 ノード数 実行時間 MPI 数 ノード数 実行時間 ( 秒 ( 秒 LAPORTA(G LAPORTA(G MPI 数 ノード数 実行時間 MPI 数 ノード数 実行時間 ( 秒 ( 秒 LAPORTA ((H LAPORTA(H MPI 数 ノード数 実行時間 MPI 数 ノード数 実行時間 ( 秒 ( 秒 LAPORTA (I LAPORTA(I MPI 数 ノード数 実行時間 MPI 数 ノード数 実行時間 ( 秒 ( 秒 S N=8 S N=8 MPI 数 ノード数 実行時間 MPI 数 ノード数 実行時間 ( 秒 ( 秒 GPGPU LAPORTA プログラム実行時間 ( 秒 D 5.85 G H I 6.88 S 理論最大性能 (GFLOPs BG/Q SR6 node 9.9 8node 78.8 GPGPU 88

123 .6.5 loop 積分 今回扱ったloop積分は d d 7 d 6 でDが次のケース(N, L, L の場合です 解析解は d 5 5 d d d d D d loop N 5 刻み幅 N L L 実行結果一覧表 h.7897 / h / h / SR6 6sp 倍精度 プログラム実行時間 ( 秒 N 9.69 L 8.7 L 67.6

124 N 8= =+ =+ 7=++7 5=+++5 6= = D = ***5**-6*7*5**-***6*5678 +****5678-5****5678 -***7*5678+5*6*7*5678 +***5*6-**7*5*6+***6** -5*7*6**+**6*5*7 -***5*7+***6*7 -*5**6*7+***7**-5*6*7**

125 L 8= =++ 67=+6+7 8= = D = -67*578***-67****8-***578*8 +*67*578*8-***8*7+***7** -*8*7**-*67***8-***7*8 +***8**-*67*8** L 8= =++ 56= = =++7 D = ***6**-*7*6**-56****678 -***7*678+*56*7*678-****678 -****678-***6*7-***6*7 +***7**-*56*7**

126 .6.6 loop 積分今回扱ったloop積分は d d d d d d CD でC,Dが次のケース(M 解析解は M M 5 は 8 は6 7 d,m d d の場合です 演算精度に関しては, 解析解と進 桁以上合わせるには M M 5 では倍精度演算で事足ります では, 拡張倍精度,倍精度演算が必要

127 M 実行結果一覧表 条件 : 分点数 5 刻み幅 h.8 / t a.666, E5-67 6sp ー O -OPENMP 言語 精度 実行時間 ( 秒 備考 FORTRAN 倍精度 6.96 C 倍精度 7.7 C 拡張倍精度 9. FORTRAN 倍精度 FORTRAN 倍精度 ー O Hpi5P SMP ー O -OPENMP 言語 精度 実行時間 ( 秒 FORTRAN 倍精度 C 倍精度 C 拡張倍精度.97 e5 8sp ー O -OPENMP 言語 精度 実行時間 ( 秒 C 倍精度

128 M 5 実行結果一覧表 条件 : 分点数 5 刻み幅 h / t a.97977, E5-67 6sp ー O -OPENMP 言語 精度 実行時間 ( 秒 備考 C 拡張倍精度 FORTRAN 倍精度 FORTRAN 倍精度 ー O FORTRAN 倍精度 75.6 ー O SR6 6sp ー O -OMP 言語 精度 実行時間 ( 秒 FORTRAN 倍精度 Phi5P sp ー O -OPENMP 言語 精度 実行時間 ( 秒 C 拡張倍精度 性能的には拡張倍精度演算を使用するのが良い事を示しています

129 M 9 = = = = = = C = 67* 78* 9* 589-9* 589* 7** - 67* 9* 8** - 67* 78* 9** + 7*** 9** D = - (67* 78* *** * *** 9* *** 78* 9* * 67* 78* 9* * * * 9* 589* 7 + *** 589* 7** - * 9* 589* 7** + 67* *** 8** + *** 9* 8** - * 67* 9* 8** - * 67* * * 8* 9 - * * * 7* 8* * *** 9** + *** 78* 9** - * 67* 78* 9** + * * * 7* 9** + * 7*** 9**

130 M 5 9 = = = = = = c = -79* 9 * 5 ** * 79 * 9 * * 589 * 7 ** - * 9 * 5 * 7 * * 79 * 9 ** * 9 * 9 ** + 5 ** * 9 ** * 589 * 9 ** + * 5 * 7 * 9 ** + 7 ** * 9 ** + * 567 * 9 ** d = 79 * ** * 5 ** + ** * 9 * 5 ** - * 79 * 9 * 5 ** * 79 * ** * * ** * 9 * ** * 79 * 9 * * 567 * 79 * 9 * * * * 9 * 589 * 7 + ** * 589 * 7 ** - * 9 * 589 * 7 ** + * * 79 * * 5 * 9 - * * * 9 * 5 * 9 + * * * 5 ** * 9 - * 567 * * * 589 * 9 + * * * 5 * 7 * 9 + * ** * 5 * 7 * 9 - * * 9 * 5* 7 * 9 - * * * 589 * 7 * * ** * 9 ** + ** * 79 * 9 ** - * 567 * 79 * 9 ** + * 567 * * * 9 ** * ** * 9 ** + ** * 9 * 9 ** - * 567 * 9 * 9 ** + * * * 5 * 9 ** - * * * 5 * 9 ** + * 5 ** * 9 ** + ** * 589 * 9 ** - * 567 * 589 * 9 ** + * * * 7 * 9 ** + * * * 7 * 9 ** + * * 5 * 7 * 9 ** + * 7 ** * 9 ** - * ** * 9 ** + * * 567 * 9 **

131 多倍長計算手法まとめ 平成 年度 - 平成 7 年度. はじめに. 多倍長精度演算. 浮動小数点演算方式. 整数演算方式. ルーチン作成前に決定すべき事項. ルーチン作成の要点.5 注意点詳細.6 整数演算方式の並列化に関して. 実アプリケーション. ヒルベルト行列 H による連立一次方程式 H=b の求解.. ヒルベルト行列の条件数.. ヒルベルト行列実測結果... 6 倍精度演算までの結果 倍精度演算での結果 倍精度演算での結果 倍精度演算での結果. 対称疎行列 A による連立一次方程式 A=b の求解. 非対称疎行列 A による連立一次方程式 A=b の求解.. n=ny=nz=65 の場合... bcg 法の収束状況... cgs 法の収束状況.. サイズを拡大した場合.. 精度改善例

132 多倍長計算手法まとめ. ファインマンループ積分.. 数値積分法.. assless 計算.. bsgaa.. 次元積分..5 two loop verte..6 loop 積分..7 loop 積分..8 loop,loop 積分作業手順まとめ.5 量子モンテカルロ法による物性スペクトル計算.5. 指数部 ビットの限界について.5. 必要演算精度見積り

133 . はじめに 倍精度演算で正しい結果が得られない場合の原因としては以下の つの場合が考えられます ( 近接する 変数の減算による桁落ち ( 演算中に表現できる数値範囲を超える場合 対策としては ( 倍精度浮動小数点数を複数個つなげた演算を行う (ieee75-8 形式の 倍精度変数を拡張する があり,( を DD 形式 ( 浮動型,( を ieee 形式 ( 整数型 と称します 場合によっては,( の形式を複数個つなげる方式 (DQ 形式 を使用します どちらの方式も一長一短がありすべてのケースを網羅する方式は現在はないと言えます また並列化により発生する問題として, 指数部の差が大きい場合の加算での演算順序や部分和の採り方により逐次実行と結果が大きくと異なる場合があります それもかならずしも並列化時の精度が悪いわけでもないと言う厄介な問題もあります 簡単な問題でその例を示します

134 c iplicit real*(a-h,o-z s=. a=. do i=, s=s+a end do write(6, s forat(h,,f.5 stop end 正しい結果は ですが次ページ以降に示すように最適化オプションや並列化で様子が異なります

135 SR6/ f9 並列化をしない場合 % f9 a.f -O -o a f9: copile start : a.f *OFORT9 V--/B entered. *progra nae = MAIN *end of copilation : MAIN %./a & % % f9 a.f -Os -noloopepand - noparallel -o a f9: copile start : a.f *OFORT9 V--/B entered. *progra nae = MAIN *end of copilation : MAIN %./a

136 % f9 a.f -Os -noparallel -o a f9: copile start : a.f *OFORT9 V--/B entered. *progra nae = MAIN *end of copilation : MAIN %./a. SR6/ f9 並列化をした場合 % f9 a.f -Os -o a f9: copile start : a.f *OFORT9 V--/B entered. *progra nae = MAIN *end of copilation : MAIN *end of copilation : _parallel_func MAIN %./a. %

137 次にスレッド数指定で異なる例を以下の例で示します 先のプログラムの を に変更しただけのものです c iplicit real*(a-h,o-z s=.d a=.d C$OMP PARALLEL DO C$OMP& reduction(+:s do i=,loop s=s+a end do C$OMP END PARALLEL DO write(6, s forat(h,,f.5 stop end

138 E5 67 6sp より結果は sp より正しい結果となる となる Phi5P sp sp.. ここで は単精度加算演算で正確に計算できる 個数となっています

139 . 多倍長精度演算 演算方式には, 浮動小数点演算方式と整数演算方式がある.. 浮動小数点演算方式 多倍長精度変数を複数個の倍精度変数の和で表す 表現できる数値範囲は,ieee75 の倍精度と同じ. a : n 倍精度変数 a, a a a,..., a a n : 倍精度変数... a 有効ビット数 = 仮数部のビット数 n倍精度 (5 n ビット n 6ビット n 59ビット n ビット n 乗加算命令がある場合とない場合の処理に分かれる コンパイラの最適化方式に性能が大きく依存する 次に示す基本演算例よりわかります 9

140 基本の演算倍精度変数の加減算, 乗除算の結果をつの倍精度変数の和で表す a, b,c,d倍精度変数 ; a b,a b,a b c d 加算 減算 : c t c a d (a (c t (b t : c t c a a b a b d (a (c t (b t 加算, 減算ともに最適化オプションで括弧をはずした演算順序に ならない様に注意が必要 乗算 ( 乗加算命令あり. c a * b d a * b c ( 最適化オプションによりd とならないように注意が必要 乗算 ( 乗加算命令なし r 779.d ( c a * b t a * r a t (t a a a a t b b * r t (t b b b b d ((a* b c a* b a* b a* b 7

141 演算量 加減算 乗算 n に比例 n に比例 n では n 個の変数の絶対値が大きい順にならべるソート処理が必要.. 整数演算方式 多倍長精度変数は ieee75-8 の 倍精度の仮数部を変えたものであらわす 表現できる数値範囲は ieee75-8 の 倍精度と同じ. a : p倍精度 ( p 変数 符号部 ビット 指数部 5ビット 仮数部 p 6ビット p 倍精度 演算量 加減算 乗算 pに比例 p に比例 有効ビット数 p p 6 p 8 ビット 77ビット ビット 符号部 指数部 仮数部 ビット 5ビット *p-6ビット (p=:,p=6: 76,p=8:

142 . ルーチン作成前に決定すべき事項 ルーチン作成にあたって以下の項目を決めて置く必要があります ( 最終ビット処理現在まで一般的に使用されていた丸め処理を適用 ( アンダーフロー, オーバーフロー処理浮動小数点演算方式では既存の倍精度演算と同じとし, 整数演算方式ではともに結果を とし, オーバーフローは使用者責任とした ( ゼロ割り処理浮動小数点演算方式では既存の倍精度演算と同じとし, 整数演算方式ではすべて結果を としこれも使用者責任とした /=,a/=.

143 . ルーチン作成の要点 最初に加減算, 乗算ルーチンを作成し平方根は既存の倍精度演算を使用して多倍長演算の初期値として反復法により求める 除算に関しては, 浮動小数点演算方式では既存の倍精度除算を複数回使用する 6 倍精度以上 ( 使用する倍精度変数が 以上の場合, 絶対値が大きい順になる様にする注意が必要 整数演算方式では逆数の初期値を既存の 倍精度除算を使用して反復法により逆数を求め, 被除数との積をとり求める 数学関数等の作成は上記の 5 演算の組み合わせで行う 具体的に立方根, 対数関数, 指数関数, 正弦関数計算でその注意事項を記述する

144 引数の取り得る範囲 ( ここでは表現可能な数値の絶対値 最大値をMとする は立方根, 指数関数, 正弦関数は M. 対数関数は 立方根, 対数関数は y M, 指数関数は 正弦関数は y ここで, 指数関数は より log(m,e 数値,zをある標準区間 とすると, 分解によって誤差は発生しにくい 立方根に関してはこの関数 が準備されていないシステムもあり, される事があるが 判れば簡単な事ではあるが, 多くの人が同じ間違いを繰り返した 問題であるので注意が必要である 指数関数及び正弦関数は テーラー展開式があるが, 引数範囲すべてで計算すると 性能, 精度の問題が発生する また対数関数では引数範囲の テーラー展開はなく で計算してを求めるが,e M ここで関数値 yのとる範囲は, 立方根に関しては より log(m となる. (yの立方根は容易に計算出来る ではエラーとなるので注意が必要である, e y z e e e 発生しにくい 対数関数と正弦関数を比較すると, 結果の範囲 が y M, と y から, 正弦関数の計算の方が誤差が出易い 事がわかる 正弦関数の計算方法の一例をあげる I [ / ], t I とし t, t より sin( sin(t, sin( t で計算する ここで問題なのは, が大きい場合, とtの持つ精度が異なる ことである ( tを求める際の桁落ちのため. また, が小さくても がnに近い場合, tは桁落ちが発生するがこの場合は sin( のため, 相対誤算は大きいが絶対誤差は小さくて済む 精度上の 問題が発生した場合には, 作成者により精度に関する思想 扱い が異なるという意識を持つことが必要となる y z y M でy e y zより, log( ( べき乗 で計算 y zの様に分解してある近似区間でz なので指数関数は誤差が y zの

145 .5 注意点詳細 ( 浮動小数点演算方式 ( ア 小さな数は精度が保持できない場合がある / 5 などでは6倍精度,8倍精度演算では 精度分の結果を得る事が出来ない ( イ 性能やソース作成の手間を考慮すると, 倍精度演算までにする必要がある ( ウ 数値表現の制限にかかる場合は, 拡張倍精度, ieee75 8の倍精度の変数をつなげた演算を行う必要がありこの場合はつなげる個数は つまで ( 整数演算方式 ( ア プログラム作成でもっとも間違えやすいのは, 近接する つの数の減算での桁落ちの数を求める部分 ( イ倍精度演算で一付近の数 に対し の桁落ちが発生する場合には,変数の8倍精度演算より倍精度変数を つつなげた8倍精度演算を使用した方が良い 拡張倍精度変数をつつなげた演算で行えれば 性能的にはより有効な方法となる ( ウ( イ と同じ様なケースが 倍超精度演算で発生する場合 変数のさらに高精度な演算ではなく,倍超精度変数を つつなげた高精度な演算を使用した方が良い

146 .6 整数演算方式の並列化に関して 整数演算方式では, 浮動小数点数 a, bを正整数 ia,ibに変換して,iaとibの加減算, 乗算, 除算を行なう. 加減算ではia,ibに 変換する(ia 一組の数に対する操作では問題ないが sidの様に core内で複数組の数を並列に操作する場合には,core内で各組毎に異なるシフト処理が出来る様にする必要がある 複数 coreで複数組の数を並列に操作する場合は各 core毎に異なるシフト処理ができるので問題はない これはsidと言ってもSR6, BG / QとHD587の様なGPUでは, 適用範囲が異なる事を意味する 乗算, 除算に関しては複数組の数で異なる処理はない 演算後の処理では 加算, 乗算, 除算は桁上げ ( 加算, 乗算, 桁下げ ( 除算 があるかないかのケースしかないので,core内でのsidが適用できる様にするには, マスク処理が出来る様にする必要がある 減算の場合, 桁下げの数を求めるのに, 左から検索して最初にが現れるビット位置をもとめる処理が必要になる core内でのsidが適用できる様にするには, この処理が出来る様にする必要がある 多倍長整数演算の場合, nビット(n 演算するが, 加減算, 乗算は市販本などでよく知られている方法で容易に演算出来るが, 除算では特別なアルゴリズムが必要となる また除算では逆数 ( a の初期近似を既存の浮動 小数点除算で求め, ib 際,a, bの値によりibのシフト数が定まる n で除算を行う方法がある a n ( a n を複数ビット毎に分けて で/ aを求めて, b / a b / a なので a の桁下がりは あるかないかのケースとなり,core内でのsid処理には, 他の演算と同じ様にマスク処理があれば良い これまでの実行では除算に関しては,倍精度から8倍精度演算ではビット数が固定なら特別なアルゴリズムの使用した除算, 可変なら任意の整数除算を使用した除算,8倍精度超精度演算では 既存の除算命令と反復法で逆数を求めて除算を行うのが良いという結果となっている n n

147 これを以下の Rup s 例題を使用して検証しました Rup's例題 a b 96. f.75b 6 a (a b 8 a 5.5b b 5767 理論解 [.75b a (a b b b ] 5.5b の演算でビットの桁落ちが発生 b b 6 8 桁落ちが影響を受けると, 受けないと となる 倍精度演算では有効ビット数がのため桁落ちの影響を受けるため.7695となり, 有効ビット数がより大きい演算 (5倍超精度演算 では桁落ちの影響を 受けないので となる

148 Core 内 sid は適用できないので core 間 sp で実行している 多倍長 rup's例題一覧表,, 回実行 実行時間一覧表 ( 秒 E5-67 Phi5P sp/core 演算精度 6sp sp 68 倍精度 倍精度 倍精度 倍精度 倍精度 倍精度 倍精度 倍精度 倍精度

149 Core 間のスレッド並列の効果が出ている 多倍長 rup's例題一覧表,, 回実行 E5-67 実行時間一覧表 ( 秒 演算精度 6sp 8sp sp sp sp 68 倍精度 倍精度 倍精度 倍精度 Phi5P sp/core 演算精度 sp 8sp sp 6sp 68 倍精度 倍精度 倍精度 倍精度 演算精度 sp sp 68 倍精度

150 a a ln( f '(a ] a a ( a a ln[( f (a a a ( a a ( (an a ( a P b,a b a P b b ( a a ( b b ( a a ( bn b((a a( ( b b a a ab bn! an!(( (( b[(an!(bn!] a( ( bn ((a bn! an!(( (( n b a P, bn j an, i j! i!(n (n [i!j!] j j(i i(n (n j j!i i!(n (n P, P n e n!~ n H H ( ln(n, n n... H j! i!(n (n!]!(j [(i j (i! j!(n i (n ( H n j i, ( j i H.. a a ( an an ( bn an bn ( an ( n n ji ji ij j i ji ij とするととすると に対して対称よりはととするとは変形して逆行列の絶対値を使用する は大として スターリングの公式またの最大値とする のノルムを逆行列オイラー数が大とするとのノルムは逆行列ヒルベルト行列ヒルベルト行列の条件数. 実アプリケーション. ヒルベルト行列 H による連立一次方程式 H=b の求解

151 n n ( /( ( ( ( /( ( ( ( n ( P.5 a ( 7 7 ( 7 (.n P.7 a.n. 7 ( 7 7 (.n. P..7,a bn j an, i a f (a.75 a.5.5n n ( log (.n log n n ( log n n n ( (.n log n.8n.6n.5979n.6n.8n a となる近似的には条件数条件数比率のより詳細にすると以下の様になる するととするが整数となる事からで最大となるはでは

152 条件数はA bを(a A( (b b で計算する場合のの大きさに関する目安となるものである ヒルベルト行列 H( ここではAとする の係数 を8倍精度演算で求めるとします. ここで解 (,,..., となる様にbを8倍精度演算で計算します. すなわち A, b の場合,8倍精度演算で計算すると, DIMENSION SIZE int 8bai seido gosa E - 6 DIMENSION SIZE 6 int 8bai seido gosa となる. ここでAの係数を倍精度演算でもとめると, A, b となり DIMENSION SIZE int 8bai seido gosa e - 6 Aの係数を倍精度演算で求め, 他の演算を 倍精度で行った場合でも DIMENSION SIZE GOSAMAX E - 7 となる T この事は, 入力データは低精度 ( もともとデータで進 桁の精度がある事は稀である. で演算のみを高精度で行い, 計算結果を低精度でも出力するという場合に適用できる 多倍長精度演算のテストでは, A, b の場合で実施

153 .. ヒルベルト行列実測結果... 6 倍精度演算までの結果 ヒルベルト行列最大誤差一覧表 次元数 倍精度 倍精度 8 倍精度 5.66.E-7.E E 次元数 6 倍精度 倍精度.7E-99.E-76.9E-86.E- 6.66E-6.6E- 8.65E-.8E-78.58E+.E-5 ヒルベルト行列の精度 N 演算精度 最大誤差 条件数 ( ビット数表示 7 倍精度.E 倍精度.E- 9 9 倍精度 7.E 倍精度.7E- 6 倍精度.7E 倍精度 9.6E 倍精度 7.98E 倍精度 7.6E 倍精度 7.E+

154 倍精度演算結果 ヒルベルト行列の精度 N 最大誤差精度.E-6 8 倍精度 5.68E-67 8 倍精度 6.E- 8 倍精度 7.9E-6 8 倍精度 8 7.6E-9 8 倍精度.E-8 88 倍精度.E-7 88 倍精度.E-5 8 倍精度 最大誤差のビット数 有効ビット数 条件数のビット数となっています 8倍精度演算ではN が使用できるメモリ容量の上限となっています

155 倍精度演算での結果 N= でのヒルベルト行列計算 条件数.7E+5 ビット数 5 精度 有効ビット数 誤差 実測 理論値 E-8.8E E-9.5E E-97.66E E-98.8E-56 ( 注 E-57.8E-56 ( 注 逆数計算で 倍精度の初期値から 7 回反復した場合の値 8,68,8 倍精度は 7 回反復,88 倍精度は 8 回反復が必要

156 倍精度演算での結果 誤差は解 (,,..., T と比較していますので 小さなNに対しては差がと正しく計算しているかの検証とは言えません そこでNを変えながら正しい事を検証しました 557 最大誤差 実行時間 条件数 有効ビット- 条件数 ( 秒 ( ビット数.D D-.D-.57.D-.67D D D D-666.5D D-5 誤差が出るサイズ 条件数 8 ビットとなる N N 条件数 ( ビット数

157 . 対称疎行列 A による連立一次方程式 A=b の求解 扱った問題 ポアソン方程式 u f (,u ( [,] [,] [,] N 収束判定値 u f : 共役残差 sin( sin( y sin( z sin( sin( y sin( z 初期値は(i. cgs, cgsilu, cgsilの非対称疎行列用反復法では収束しない この様な場合, 初期値は (i : (, の一様乱数値としている

158 倍精度演算実行での反復回数 事になる 反復回数 解法 反復回数 備考 sor 9 陰解法 sor 799 陽解法 sor 75 陰解法 +odd-even 法 adi 95 陽解法 cg 69 bcg 7 bicgs 95 cgs 6 初期値 (i=.では収束しない cgs 88 cgsiluの並列版 前処理付き scg 69 前処理付きの最も簡単なケース cgsil 初期値 (i=.では収束しない cgsilu 7 初期値 (i=.では収束しない cgs 6 cgsilの並列版 gpbicg 並列実行にすると, 初期値 (i.でも 収束するように改善される場合がある 反復回数の削減には, 初期値の変更, 並列化などがある

159 . 非対称疎行列 A による連立一次方程式 A=b の求解サイズは精度改善により拡大しています 初期値収束判定値は共役残差法法解法非対称問題 解析解領域次元ポアソン方程式問題は z y, u(, 9 nz ny n, R nz ny n, R 65 nz ny n, R square gradient cgs(conjugate 9 nz ny n, R 65 nz ny n, R gradient bcg(biconjugate z sin( y sin( sin( e z y, u(, [,] [,] [,] f u R u yz

160 .. n=ny=nz=65 の場合... BCG 法の収束状況 反復回数 5 bcg 法 倍精度演算では反復回数 でも収束せず 反復回数 ヒ ット数

161 -log( 残差 倍精度 反復回数

162 BCG 法の収束状況 n=ny=nz=65 収束判定値 =.q- -log( 残差 ヒ ット 8 ヒ ット 反復回数

163 -log( 残差 倍精度 5 倍精度 6 倍精度 7 倍精度 8 倍精度 - 反復回数

164 CGS 法の収束状況 n=ny=nz=65 収束判定値 =.q- 倍精度演算では収束しない 反復回数 5 cgs 法 反復回数 ヒ ット数

165 ... CGS 法の収束状況 -log( 残差 5 倍精度 5 倍精度 (n= 反復回数

166 -log( 残差 反復回数 6 倍精度,8 倍精度 6 倍精度 (n= 8 倍精度 (n=

167 cgs 法反復回数 倍精度は 回で収束せず 反復回数 精度 ieee ieee sp dd 倍精度 倍精度 6 55 none 6 倍精度 7 倍精度 8 8 none 8 倍精度 qsu,q5su 使用 精度 反復回数 倍精度 78 5 倍精度 55 (ieee 形式,dd 形式ともに 6 倍精度演算が最も効率が良い事を示しています (qsu,q5su の様に総和演算を無限精度演算を使用すると,5 倍精度の収束が良くなり,ieee での sp 実行と同じ反復回数となっています

168 詳細な収束状況 dd 形式 cgs 反復回数 ビット数 反復回数

169 .. サイズを拡大した場合 通常サイズを大きくすると収束に要する反復回数は増えますが その状況を以下の条件で検証しました 解法 bcg(biconjugate gradient 法 R, n cgs(conjugate R, n 収束判定値は共役残差 初期値 u(, y,z の場合 ny nz 9 gradient square 法 ny nz 反復回数 まででは bcg 法は倍精度演算, cgs 法は 倍精度演算まで収束しません

170 bcg 法反復回数一覧 演算精度 SR6 557 E5-67 Phi5P 倍精度 8 倍精度 倍精度 7 6 倍精度 7 7 倍精度 7 8 倍精度 7 cgs 法反復回数一覧表 演算精度 SR6 557 E5-67 Phi5P 5 倍精度 倍精度 倍精度 倍精度 ( 注. 赤字はその計算機で最も性能がよかったもの.6倍精度の有効ビット数はSR6は 59ビット, 他は77ビット.E5 67, Phi5Pでは内積部分に無限精度演算を適用しています

171 bcg 法反復回数 精度 反復回数 正規化 実測 / 正規化 56ビット 5 none none 6ビット 5 node none 6ビット ビット ビット ビット 倍精度 ビット ビット ビット ビット ビット ビット ビット.96 倍精度 8 8 ( 注 5 回は 5 回反復で収束しなかった事を示しています 正規化 P 有効ビット数 = 精度のビット数 + 正規化反復回数 8 ( P 性能的には 倍精度が最も効率が良いと言えます

172 .. 精度改善例 以下の様に精度改善を行う事により, これまでのサイズおよびサイズをさらに拡大しても 倍精度演算で収束するようになりました 前処理付きcgs法により, 並列処理で精度改善 cgs(conjugate gradient square 法 R, n ny nz 9 収束判定値は共役残差 初期値 u(, y, z 倍精度演算 CGS (PARALLEL SIZE NX, NY, NZ = 9 EPS =.E - ITERATION = 9 9

173 . ファインマンループ積分 ファインマンループ積分の一般式は以下の様になる I( ( N ( nl n : 時空次元 (n (N : 実定数,C,D :, N :, nl ループ内の素粒子の数,L : ループの多重度,..., N N d ( i i の多項式... N N n ( L C (D ic NnL 演算精度が問題となるのは, 以下の例の様に仮想光子 λ が含まれる場合である ( Infra vt I(s dyd D (D ys ( y e ( y ( Infra bo y I(s dzdyd D (D ys tz( y z ( ( y z( y e y z( y f

174 .. 数値積分法 数値積分法としては ( 誤差を評価しながら領域を二分割して行く Gauss-Kronrod 積分法 ( 端点特異点を持つ被積分関数の積分に強く, 基本的には領域は分割せず, 標本点をふやして積分計算結果の精度をあげていく二重指数関数型積分法がある これ以外にも一重指数関数型積分法もあるが, 二重指数関数型積分法に比べ標本点数を大きくとる必要がある 端点と領域内部で同時に特異点をもつ 場合の計算を, sin c求積法によるhadaard 有限部分計算で行う方法もあるがファインマン ループ積分への適用出来る例は少ない 数値積分で積分領域内に特異点がある場合 微小量 を使用してD D iと有理化して で得た値を求める積分値とする場合に 加速法を使用するが主にRicardson加速法, エイトケン加速法 Wynnの 算法等を用いる

175 Gauss Kronrod積分法と二重指数関数型積分法 の選択に関しては以下の差があります 分点と重み係数のテーブル作成方法では Gauss Kronrod法では, まずGauss 積分法で n の分点と重み係数を求めます これにn を求めます これは反復法により求めるため, 積分計算の 前に分点と重み係数の計算が必要となります これに対し, 二重指数関数型積分は分点と重み係数は 関数計算 (sinh, cosh, ep で求まるため, 積分計算のプログラム に入れる事が可能です 二重指数関数型積分の誤差評価 I 誤差 は, N h の分点を追加して計 n n で求めた値 w n i w ep( c... i とI n h n i n, w w,..., w N で評価ができます log(n 個の分点と重み係数 を比較して求められますし, 内部特異点がある場合は領域分割によって計算が可能です i n n n で... n 幾つかの例を次ページに記しました

176 例題 倍精度演算, 分点数 96, 変数変換区間 [,] で計算. 分母 dでd それぞれの領域でd d i( 算法で計算 となる点で領域を分割し, d ( e 計算結果 ,...,. (sin cos d 計算結果 sin.cos. 計算結果 d d log log( 計算結果 ( 注 この積分値は約 年間 と考えられていました

177 ( ( ( ] 8 ( [ I ( ( 8 ( dy y y ( y dy y ( y dy y ( y dy y y ( y J ( ( dy y y ( y dy y ( y dy y ( y dy y y ( y f (d I y ( y, dy d ( y ( ( ( ( 6 ( ( ( ( ( ( ( f ( ln(n N C ep(.76 C N., h ( ( ( f ( f (d I 倍精度演算実測値の相対誤差

178 被積分関数が内部特異点を持つ場合, D->D-iε で有理化して, ε-> として積分値を求めるのが一般的であるが収束が遅くなったり, 異なる積分値を得る場合がある ( inf では s の場合の解析近似解は ( [ln( s でs s ra dyd D verte e (D sy ( y e ln( s ln e ( s ( y ] 6 のときは,s s iで解析接続して求める 5 倍精度演算で計算した時の実行時間 ( 秒 は D D i s s i e e y,.5でgauss Kronrod法 以下の様になっています