スーパーコンピュータ有効活用の手引き(平成23年度-平成27年度)
平成23年度から高性能計算の扉に記載しました各種計算機基本性能調査、各種計算機アプリケーション性能比較、多倍長計
算手法のまとめとして性能関連事項を各種計算機性能調査まと
めと、精度関連事項を多倍長計算手法まとめの2冊にしました。
性能関連事項と精度関連事項は異なる視点からの検討が必要な部分が多くあるため、この2冊を列挙する体裁を取っています。
論文リストはスライド238以降にまとめて掲載しています。

スーパーコンピュータ使用するケースは主に以下の3つのケースに分類できます。
(1)大学等での計算機教育や学問としての計算機論の展開(2)研究所等での論文作成の手段として(3)業務を行う公的機関や民間会社での手段として一般的には(1)と(3)のサイトでは目的と手段が明確になっていますが(2)のサイトでは目的と手段が不明確で逆転している事が多々見られます
一例として過去ある大学で逆行列を求めるサブルーチン群を学内で公開した所,参照回数の50%以上がクラメルの公式で求めるサブルーチンだった事があります。
(2)や(3)の分野の人から見れば,時間もかかり,計算精度からみて問題のあるものを使用するのはおかしく見えますが、教育という面からみると,連続変数を扱う数学と離散変数を扱う計算機の差や,計算機の演算量や精度を知る上では目的にかなっていると言えます。

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

また実行時間に関しては、PCやサーバーから移植する場合には互換性などの問題も含め経過時間が問題となりますので測定結果はすべて経過時間としています。たとえ1CPUでの測定でもCPU時間ではなく経過時間としています。

各種計算機性能調査まとめ
平成２3年度-平成27年度
目次
1.はじめに
2.使用計算機
2.1 スーパーコンピュータ
2.2 アクセラレータ
2.3 サーバー
3.基本演算3.1 メモリバンド幅測定
3.2 rump’s例題
3.2.1 ieee754-2008系
3.2.2 power系
3.2.3 超多倍長演算
3.3 行列積,ストラッセン行列積
3.3.1 行列積計算
3.3.1.1 SR16000/M1 システム
3.3.1.2 xm1 システム
3.3.1.3 BG/Qシステム
3.3.1.4 各種演算精度の行列積計算
3.3.2 ストラッセンの行列積
4.実アプリケーション
4.1 重力多体問題(N体問題)
4.1.1 SR16000/M1 システム
4.1.2 BG/Qシステム
4.1.3 各種計算機性能比較
4.1.3.1 演算量可変ケース
4.1.3.2 演算量一定ケース

各種計算機性能調査まとめ
平成２3年度-平成27年度
4.2 分子動力学計算
4.2.1 SR16000/M1 システム詳細
4.2.2 各種計算機性能比較
4.3 大規模疎行列の反復計算
4.3.1 対称問題
4.3.2 非対称問題
4.4 QCD計算
4.5 量子モンテカルロ法による物性スペクトル計算
4.5.1 SR16000/M1 システム
4.5.2 ieee754-2008形式
4.6 ファインマンループ積分計算
4.6.1 infra vtx計算
4.6.2 infra box計算
4.6.3 massless計算
4.6.4 4-6次元積分計算
4.6.4.1 s221計算
4.6.4.2 5次元積分計算
4.6.4.3 6次元積分計算
4.6.5 3loop計算
4.6.6 4loop計算

1. はじめに性能を考える場合,実行時間をもとにしても,一般的に以下の事が言えます。
(1)扱う問題がある程度特定される場合とそうでない場合の性能は大きく異なる場合がある。
(2)カタログ性能と実行時間は必ずしも比例しない。
(3)最適化オプションの効果には注意する必要がある。
具体例を挙げると以下の様なものがあります。
（１）複素変数の絶対値計算c  a  ib
| c |

a 2  b2

a , bが絶対値がある範囲で ある事がわかっている 場合
は, 実数演算で 2つの乗算 ,1つの加算 ,1つの平方根計算で済む.c  10160  i10160 の場合 , cは表現可能な数値
で, 結果も表現可能な値と なるが、同じ計算をす ると演算途中でオーバーフ ローが発生します。
これを防ぐため , 関数ルーチンでは大小 判定や除算等が使用されるのが一般 的で数倍の性能差が出 ます。

(2)カタログ性能は計算機に備わる演算器が全て同時に動作する事を前提としているが,実際のアプリケーションでは必ずしもそうはならない。倍精度演算と4倍精度演算の性能
比はSR11000では約10倍と言われてきたがSR16000では4倍精度演算にはSIMD命令が適用されないためその差は大きくなり、
15-20倍となっている。

(3)最適化オプションでは結果が変わらない範囲
での最適化と、結果が変わる可能性がある
最適化のレベルがある。
特に,並列実行の際には,精度に敏感なプログラムで総和,内積演算を含む場合には性能以前に結果の検証が必要になります。

2.使用計算機
2.1 スーパーコンピュータ(1) SR11000/K1(1ノード)
プロセッサ:power5
周波数:2.1GHz
CPUコア数 16
論理コア数 16
理論最大性能 134.4 GFLOPs
メモリ容量 32GB
メモリアーキテクチャー Flat Memory Interleave
L3キャッシュ Off-Chip 36MB/2コア
演算器/物理コア 乗加算器2つ(2) SR16000/M1
プロセッサ:power7
周波数:3.83GHz
CPUコア数 32(物理的),64(論理的)理論最大性能 980.48 GFLOPs
メモリ容量 256GB
メモリアーキテクチャー NUMA,(16論理コア単位でflat)
SIMD(Single Instruction Multiple Data)をサポートするVSX機構付き
L3キャッシュ On-Chip 32MB/8コア
演算器/物理コア 乗加算器４つ

(3) SR16000/xm1
SR16000/xm1は周波数が3.3GHzで他はSR16000/M1と同じです。 演算性能だけみれば,SR16000/M1 1ノードはSR16000/xm1の１６％性能向上版ともいえます。
またメモリ性能はSR16000/M1 の半分です。

この3機種では,(2)のSR16000/M1 がメインで(1)のSR11000/K1 は互換性,実行性能比を,
SR16000/Xm1 はメモリ性能の影響と,smtの効果を見ています。
(4)BG/L
周波数 0.7GHz
1 ノード 2 論理性能 5734.4 GFLOPs
L1キャシュ 32KB(データ,命令),L2キャシュ 2KB,L3キャシュ4MB(5) BG/Q
周波数 1.6GHz
1 ノード 16core 論理性能 204.8GFLOPs
L1 キャシュ
16/16KB (Core)
L2 32MB (node)
Main storage 16GB (Core)
Smt=1,2,4(4),(5)に関してはメインは(5)で(4)は互換性を見ています。

(6)T2K 筑波システム
AMD quad-core Opteron 8000 シリーズ
(Barcelona)
1 node:ピーク性能
147GFLOPs,16MPI/node2.2 アクセラレータ(1)HD5870
GPU カード型番:ATI RadeonHD5870
メモリ: GDDR5, 1 GB, 153.6 GB/s
ホストインタフェース: PCI Express 2.1
x16stream
processing unit: 3200個(演算プロセッサ)
動作周波数: 850 MHz ピーク性能(倍精度): 1088 Gflops(2)HD6970
GPU カード型番:ATI RadeonHD6970
メモリ: GDDR5, 2 GB, 176 GB/s
processing unit: 6144個(演算プロセッサ)
動作周波数: 880 MHzピーク性能(倍精度):2703 Gflops

(3)グラフィックボード
HOST
E5  2670
キャシュ

2.60GHz

1cpu  8core

20MB

2.6GHz  8  8  166.4GFLOPs

2cpu

2cpu  332.8GFLOPs
Xeon

Ph i5110P

1.053GHz

60コア,4 スレッド/1core
1.053GHz  60  4  4  1010.88GFLOPs

2.3 サーバー(1)x5570 8コア 2.93GHz
キャシュ 8MB/コア(2)e5430 16コア 2.66GHz
キャシュ 12MB/コア(3)上記グラフィックボードのＨＯＳＴ
E5-2670もここに含めました。

3.基本演算
3.1 メモリバンド幅測定

メモリバンドテスト
copy
do j  1, N

trans
do j  1, N

do i  1, Na ( i , j)  b ( i , j)

do

end

do

end

do

end

do

end

do

i  1, Na ( i , j)  b ( j, i )性能測定結果(MB/sec)
N演算8192 copytrans8193 copytransSR16000 xm1
E5-2670 Phi5110P
32core32core16smp
240smp216248
84717
25687
80130
3223
958
4692
8508
221165
78572
20112
78970
31387
7216
17959
10018

詳細

transは参照非連続

N=8192

連続
非連続
非連続
非連続

非連続
要素
なし
参照
格納
両方

sr16000 sr16000 E5-2670 Phi5110P
32core64smp
16smp
240smp216248189155
25687
80130
3223
3082
4692
8508
2282
3461
5103
6927
871
1849
24912
10354

非連続
要素
なし
参照
格納
両方

sr16000 sr16000 E5-2670 Phi5110P
32core64smp
16smp
240smp221165182590
20112
78970
31387
30856
17959
10018
56743
52338
20496
16833
54351
43623
18743
12710

非連続
要素
なし
参照

smt2
smt2
smt4
smt4
smt4
32smp
64smp
32smp
64smp
128smp
84717
86216
84754
84639
71859
958
961
960
948
937

非連続
要素
なし
参照

smt2
smt2
smt4
smt4
smt4
32smp
64smp
32smp
64smp
128smp
78572
68711
78484
70198
79770
7216
6328
6297
6078
6173

N=8193

連続
非連続
非連続
非連続
xm1
N=8192連続
非連続
N=8193連続
非連続

演算を含めた場合の結果(MMAX=8193)

2D Copy
do j = 1,MMAX
do i = 1,MMAXr(i,j)=p(i,j)
enddo
enddo

2D Scale

2D Adddo j = 1,MMAX
do i = 1,MMAXp(i,j)=scalar*q(i,j)
enddo
enddo

do ｊ = 1,MMAX
do i = 1,MMAXq(i,j)=p(i,j)+r(i,j)
enddo
enddo

2D Triad

Transpose

do j=1,MMAX
do i=1,MMAXr(i,j) = p(i,j) + scalar*q(i,j)
enddo
enddo

do j=1,MMAX
do i=1,MMAXr(i,j)=p(j,i)
enddo
enddo

SR16000/XM1とSR16000/M1の2機種で測定していますが,約３倍の差があり,演算性能の差16%とは大きな差があります。
特に,Transpose 64スレッドでの差が顕著です。

メモリバンド幅性能測定結果(MB/sec)
SR16000/XM1 MB/sec
function 1スレッド ２スレッド ４スレッド ８スレッド １６スレッド ３２スレッド ６４スレッド
2D Copy
5070 8004 13786 27134 43057 68565 62532
2D Scale4682 7982 13707 27011 41294 61917 58637
2D ADD6563 11054 16604 32834 48544 78820 70376
2D Triad
7307 12097 16954 33525 49752 77751 70680
Transpose
426
725 1470 2649 4883 3801 2866
SR16000/M1 MB/sec
function 1スレッド ２スレッド ４スレッド ８スレッド １６スレッド ３２スレッド ６４スレッド
2D Copy
11661 23356 46292 92051 165536 138902 148588
2D Scale 12528 25012 49218 97630 170661 140657 163474
2D ADD
12918 25797 49839 97796 176648 151006 174514
2D Triad
12680 25216 48697 96034 174808 149709 171588
Transpose 1052 2055 4072 7971 15921 16231 30540

3.2 rump’s

例題

rump' s例題とは , 四則演算からなる簡単 な
計算で , 有効ビット数が121ビット以上を要求する計算のため , 精度、性能の検証によく使用されるもので す。
a  77617.0, b  33096.0
f  333.75b 6  a 2 (11a 2 b 2  b 6  121b 4  2)  5.5b 8a2ba 2  5.5b 2  1  6024398689


a 2 (11a 2 b 2  b 6  121b 4  2)  5.5b 2  333.75b 6  2
より, f  2 

a
54767

 0.827396059946
2b
66192

3.2.1 ieee754-2008系
拡張4倍精度：拡張倍精度変数を２つつなげたもの6倍精度：倍精度変数を３つつなげたもの
8倍精度：倍精度変数を4つつなげたもの
rump例題テスト結果
e5430

実行結果一覧表
実行回数
演算量

100,000,000
1.9GFLOP
(各精度)

精度

言語

拡張４倍精度拡張４倍精度
6倍精度
6倍精度
8倍精度
8倍精度

実行時間
smp実行時間 台数効果(秒)
8smp(秒)icc(strict)
45.5899
5.817464
7.84
icc(extended)
37.5245
4.79337.83
icc
54.6475
6.94647.87
ifort
55.40236.98077.94
icc
125.8814
15.8532
7.94
ifort
126.4694
15.7271
8.04

t2k
実行回数60,000,000
演算量
1.14GFLOP
精度
拡張４倍精度拡張４倍精度
6倍精度
8倍精度

言語

(各精度)

実行時間
smp実行時間 台数効果(秒)
16smp(秒)
icc(strict)
58.88582
3.766515.63
icc(extended) 50.962057
3.250428
15.68
icc
48.549858
3.132938
15.5
icc
113.161652
7.336832
15.42拡張4倍精度はextended オプション効果あり。
(両CPUで)Fortranとc++ での有意な性能差はない。
拡張４倍精度と6倍精度の性能は,
e5430 は拡張４倍精度、Ｔ２ｋでは6倍精度が良いという逆の傾向がでています。

rump ' s

例題実行結果一覧
100,000,000回実行

実行時間(秒)CFORTRAN精度
E5-2670 Phi5110P E5-2670 Phi5110P
16smp
240core 16smp
240core拡張倍精度*2 1.949568 2.993835
6倍精度
2.391048 3.373063
2.35453.33678倍精度
5.514834 5.949745
5.39945.864610倍精度
17.7158
15.0406

rump例題実行時間(秒)一覧表
整数演算方式
精度
5倍精度
6倍精度
7倍精度
8倍精度
10倍精度
12倍精度
１４倍精度
16倍精度

E5-2670 Phi5110P
6.765920.3040
8.270624.5392
10.0500
28.9972
12.9452
35.601717.0668
52.6718
25.220578.2484
32.5928 102.1592
39.9565 127.1153

Phi5110Pのthreads/core実行時間(秒)改良版
threads
12
3
4

10倍精度 12倍精度
81.3013 124.3947
60.672392.756257.749486.7546
52.6718
78.2484

DD形式(一般サブルーチン使用)

精度
10倍精度
12倍精度
14倍精度
16倍精度

E5-2670 Phi5110P
16smp
240core15.8159
27.930227.741946.732046.821170.783172.6414 107.0264

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

E5-2660 rump's 例題実行時間一覧表
100,000,000回実行の実行時間(実行時間:秒)
smp数

拡張倍精度*2 6倍精度
6倍精度
8倍精度
8倍精度
c
c
F
c
F
1
27.61538.20533.16296.05778.7232
13.85219.391
16.62648.81539.5684
7.248
10.081
8.607
25.221
20.414
8
4.315
5.367
4.643
13.54610.986
16
3.038
3.729
3.517
7.131
8.288
32
2.257
3.562
2.755
7.803
6.412

拡張倍精度*2 及び6倍精度FORTRANで32smpの場合の性能が良い。

3.2.2 ＰＯＷＥＲ系

1スレッド実行

fma,nofma

SR16000/M1
1,000,000 回実行

コンパイラ f90
f90
xlf
xlf
精度
fma
nofma
fma
nofma
（秒）
（秒）
（秒）
（秒）
4倍精度
0.241
0.249
0.15
0.21
6倍精度
0.514
0.453
0.31
0.39
8倍精度
0.937
0.849
0.56
0.73

8倍精度演算方式
QD方式:倍精度変数を4つつなげた方式
DQ方式:4倍精度変数を2つつなげた方式(注)(1)DQ方式はQD方式より精度が悪くなる。
最悪の場合ＴＤ(倍精度変数を3つつなげたもの)と同じ精度になる。
(2)ieee754-2008でのDQ方式の有効ビット数は226ビットで,DD方式の
ＱＤ方式より精度が良い。
Rump's例題８倍精度演算
100,000,000回実行の時間(秒)SR16000/M1
smp数

１ノード実行

QD方式
ＤＱ方式
1
99.22
334.012
52.95
168.78
4
31.49
86.83
8
36.55
47.55
16
45.68
42.06
32
54.52
47.67
64
49.09
43.3

4.2.3 超多倍長演算

多倍長rump' s例題一覧表
1,000,000回実行
実行時間一覧表(秒)E5-2670
演算精度
68倍精度
128倍精度188倍精度248倍精度308倍精度368倍精度428倍精度488倍精度548倍精度

16smp
8smp
4smp
2smp
1smp
11.9698
23.8582
47.858595.5189191.5384
42.3808
84.7253
169.0772
338.4619
677.3167
89.0890
177.8675
355.6048
711.6836 1425.8922
152.0543
303.7099
607.2627 1214.3922 2438.3971
248.2340
355.4690
476.9385
617.0996
833.8722

Phi5110P

4smp/core

演算精度
68倍精度
128倍精度188倍精度248倍精度308倍精度368倍精度428倍精度488倍精度548倍精度

240smp180smp120smp60smp
46.3840
48.7955
50.8138
69.1542156.0239
178.3118
181.1152244.4755
331.2860
378.7431
384.5338
525.7136
585.7436
669.6312694.6440
938.0996
966.5418
1378.1645
68倍精度
1859.8102
4smp
1smp
2410.4360
868.7135 3463.0626
3244.1159

ともに並列化効果は出ていますが,E5-2670の効果が大きくなっています。

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

3.3.1 行列積計算
3.3.1.1 SR16000/M1 システム

行列積計算を行う場合,行列のサイズNが小さい場合,
3つの行列すべて合わせてもキャシュに収まる範囲にあります。
行列積計算でC  AT Bを実行する場合の総メモリ所要量N

実数型
32 24KB
64 96KB
96 216KB
128 384KB
256 1536KB
512 6MB768 13.5MB

複素数型
48KB
192KB
432KB
768KB
3MB
12MB
27MB

赤色 L1D内に収まる。
緑色 L2内に収まる。
紫色 L3内に収まる。

N>100 (3桁)でキャシュチューニングが重要となる。

キャシュチューニングの際,内積型計算を使用する場合
主要演算時の所要メモリ量N(N+2)*8or16が2MB以下に抑える事が重要になります。

主要演算時の所要メモリ量

実数型
複素数型
N

N  ( N  2)  8バイトN  ( N  2)  16バイト

実数型
31 8KB63 32KB
64 33KB
100 80KB
192 291KB
256 516KB
289 657KB
384 1158KB

複素数型
16KB
64KB
66KB
160KB
582KB
1032KB
1314KB
2316KB

SR16000 / M 1 1ノードで内積型行列積計算で性能が出る条件(1) simd命令が適用される。
(2)最内側ループ (k )内の演算で 24の乗加算命令が適用される(3)最外側ループ ( j )のループ長が 4  32  128の倍数(4)中間のループ (i )のループ長が 6  2 n の倍数(5)主要演算時の所要メモリ量が 2 MB以下

実数型ではN=384が上記すべての条件を満たす。
複素数型ではすべてをみたすサイズはないのでＮに応じて最適なものを探します。
これらの事からＮが３桁でも実数型ではＮ＝５０４，
複素数型ではＮ＝３６０を境にチューニング手法を変える必要が出てきます。

実数型行列積計算
実数型行列積計算性能一覧表

内積タイプとは, 以下の3重DO ループでDO
j  1, n, lj
DOi  1, n, li
DOk  1, nlj * li で示しています.
1ノード smt=off,smp=32 実行結果一覧
N

内積タイプ
32 1*1664 2*896 3*8128 4*4
128 4*8
192 6*6
288 3*6
384 4*6

GFLOPs425.109645.747725.319807.875
807.924827.583
800.614860.561

定義式どうりのソースでの結果(lj=1,li=1)
N

32smp
64smp
31
105.7891.08061145.624 221.050
67
210.445 211.297
97
253.611 320.202
127
333.003 424.066

性能向上要件(5)のテスト

N  ( N  2 )  8 バイト  2 MB
での実数型行列積計算の性能

1ノードsmt=off,32smp
内積タイプ N=480 N=504
4*6
819 833.002
6*6
739 780.324
4*8
813 849.898
4*4
827 831.878

複素数変数行列積計算
複素数型行列積計算性能一覧表

内積タイプとは, 以下の3重DO ループで
DO j  1, n, lj
DO i  1, n, li
DO k  1, n
lj * liで示しています.
1ノード smt=off,smp=32 実行結果一覧
N

内積形式 GFLOPs
32 1*16227.79664 2*4454.764
64 2*8451.572
96 3*4602.258
128 4*4
470.995
192 3*6
677.963
192 3*4
719.704
192 6*3
817.149288 9*1
774.796288 3*6
739.758
288 3*4
809.752

定義式どうりのソースでの結果(lj=1,li=1)
N

32smp
64smp
31 86.505
76.37761 254.475 305.272
67 191.442 182.165
97 206.989 262.580
127 273.559 348.667

3.3.1.2 xm１システム

smt=４の効果検証
検証は複素数変数行列積と実数変数行列積計算でsimd , nosimdで行いました。
(1)複素行列積C  ( A B ) A を28000回計算する。
T

T

定義式どうりに計算する複素数変数プログラムCMULT
実数変数でて定義どうり計算するプログラムRMULT
RMULTはsimd , nosimdで計算(2)実行列積C  AT Bを計算する。
実行回数は N  288,384

100,000回

N  480,504

10,000回

アンローリング数 は最外側4, 中間6プログラムは MULTIN

実数変数行列積計算C  AT B
multin4*6タイプsimd

N=384,演算量=11325GFLOP
実行時間（秒)一覧表
smp数

実行時間 GFLOPs
32
16.436689
6434.195331
9624.207468
128
25.639442

N=288

演算量＝4478ＧＦＬＯＰ
実行時間（秒)一覧表

smp数
32
64
72
96
128

N=480
N=504

smp

実行時間 GFLOPs
10.116443
18.239246
14.379311
13.047343
17.722253

演算量＝2212ＧＦＬＯＰ
演算量＝2560ＧＦＬＯＰ
実行時間（秒)一覧表
n=480
n=480
n=504
n=504
実行時間 GFLOPs
実行時間 GFLOPs
32
3.477
636
3.802
673
64
5.626
393
5.816
440
96
6.898
321
7.566
338
128
6.055
365
6.650
385

最高性能一覧表N

GFLOPs
実行効率(%)
288
443
52
384
689
82
480
636
75
504
673
80

最高性能はSR16000/M1と同じN=384で達成

32smpでのsmt=2,smt=4の性能比較

smt  2とsmt  4での32smp性能比較
内積タイプ 4*6
N

smt=2
smt=4
演算量
実行時間 性能
演算量
実行時間 性能(GFLOP) (秒)(GFLOPs) (GFLOP) (秒)(GFLOPs)
288
4478
11.301396
4478
10.116
443
384
11325
17.645643
11325
16.436
689
408
1358
2.651
512
1358
2.771
490
432
1612
2.838
568
1612
2.927
488
456
1896
3.626
523
1896
3.394
559
480
2212
3.638
608
2212
3.599
615
504
2560
4.02
637
2560
3.802
673
528
2944
5.39
546
2944
5.351
560
768
9060
66.18
137
9060
58.292155
内積タイプ 3*4

N

smt=2
smt=4
演算量
実行時間 性能
演算量
実行時間 性能(GFLOP) (秒)(GFLOPs) (GFLOP) (秒)(GFLOPs)
192
1416
3.269
433
1416
3.008
471
288
4478
12.048
372
4478
8.869
505



複素変数行列積計算C  ( A B) A
T

T

N=289,演算量=10814GFLOP
実行時間（秒)一覧表
smp数
simd

rmult
rmult
cmult
cmult
on
off
on
off
32
44.887193.428
71.38150.7866439.372206.027
58.25460.9959746.45584.65262.35
51.837128
64.08684.47672.36
52.887

N=288,演算量=10702GFLOP
実行時間（秒)一覧表
smp数
simd

rmult
rmult
cmult
cmult
on
off
on
off
32
31.29987.59456.23587.0816436.029107.32458.96467.3629630.183
159.467
46.90248.036
12834.372110.54455.66
86.615

N=31,演算量=13.346368GFLOP
実行時間（秒)一覧表
smp数

rmult
rmult
on
off
32
0.124
64
0.155
96
0.195
128
0.321

cmult
cmult
on
off
0.176
0.189
0.193
0.205
0.212
0.212
0.236
0.395
0.219
0.219
0.434
0.261

Smp数32,96が良い値となっています。

3.3.1.3 BG/Qシステム
複素数変数行列積計算

CT  A T Bの計算を複素数変数のままと、実数変数に直した場合を実行。
SIMD命令を適用させるため,
 qstrict

オプションは外し ています。

測定条件は
1ノード

openmpによる並列化

bg

複素変数行列積

C T  AT B
- O5 - qarch  qp - qtune  qp - qcache  auto
- qhot  level  2 - qipa  level  2 - qsimd  auto
- qsmp  omp
複素数型
n

実行時間一覧表（秒)
演算量 タイプ
16smp 32smp 48smp 64smp(GFLOP)32
734 1*16
78.359 97.998
64
587 2*825.179 22.174
96
198 3*4
4.094
3.930
3.831
128
470 2*4
10.283 11.250
8.200
192
1585 3*4
32.378 25.078 22.773
22.917288
5351 3*4
102.125 64.950 70.568
384
12684 3*4
231.934 154.786
156.162(注)48smpでは,96,192はタイプ2*4にしている。

最高性能は82GFLOPsbg

複素変数行列積

C T  AT B
-複素数型
O3 - qarch  qp - qtune  qp - qcache  auto - qhot  level  1- qipa  level  1

ーＯ３- qsimd  auto
- qsmp  omp
複素数型
n

実行時間一覧表（秒)
演算量 タイプ 16smp 32smp 48smp 64smp(GFLOP)96
198 3*44.075 3.448
3.673
128
470 2*4
10.155 7.928
8.123
192
1585 3*4
32.373 20.034 22.314 20.523
2563758 2*4
81.983 49.812
57.069288
5351 3*4
100.029 63.432 70.078
384 12684 3*4
238.699 153.374
155.539(注)48smpでは,96,192はタイプ2*4にしている。

性能は84GFLOPs,83GFLOPsと-O5の場合より少し良くなっています。（誤差の範囲)

bg

複素変数行列積

C T  AT B
- O5 - qarch  qp - qtune  qp - qcache  auto
- qhot  level  2 - qipa  level  2 - qsimd  auto
- qsmp  omp
実数型
n

実行時間一覧表（秒)
演算量 タイプ
16smp 32smp 48smp 64smp(GFLOP)32
734 1*16
125.247 124.539
64
587 2*8
19.187 20.635
96
198 3*44.894 4.909
5.200
128
470 2*4
10.513 10.522
10.205192
1585 3*4
43.199 33.340
24.860288
5351 3*4
123.406 106.704 87.407
384 12684 3*4
263.06 256.902
187.866(注)48smpでは,96,192はタイプ2*4にしている。

ソースはSR16000/M1,XM1と同じものを使用。
性能は64GFLOPs,61GFLOPs,68GFLOPsと複素変数の場合に比べて有意に劣っています。

3.3.1.4 各種演算精度の 行列積計算
比較はSR16000 1ノード(32core,64smp),ES-2670 16smp,
Phi5110P 240core で実施しています。最も多くのケースを実施したのは行列サイズN=240,N=480の1000回実行です。
Cは,ES-2670,Phi5110Ｐで倍精度,拡張倍精度,拡張倍精度＋拡張倍精度のケースを行っています。名称はそれぞれ
Double,exdouble,ddexdoubleとしています。FORTRANは
SR16000,ES-2670,Phi5110Pで倍精度multd,4倍精度(ieee754-2008データ形式)multq,dd形式4倍精度ddmultd,dd形式の6倍精度,8倍精度,10倍精度,それぞれmult6,mult8,mult10
です.SR16000では特にmultd,multqは実施していません。
また比較という事で,テストプログラムは定義式どうりのコーデイングでアンローリング,キャシュチューニング等は実施していません。演算量は以下の様になっています。

演算量(GFLOP)一覧
プログラム N=240
N=480
multd
28
221
ddmultd
470
3760
mult6
1849
14789
mult8
3866
30928
mult108657
69265

行列サイズN=240の結果は以下の様になっています。
行列積N=240,1000回実行の性能比較
実行時間(秒)プログラム
multd
ddmultd
multq
doubleexdouble
ddexdouble
mult6
mult8
mult10

SR16000
32core64smp
7.360

21.27437.59979.776

phi5110P E5-2670
240smp
16smp
1.824
0.403
3.753
2.952
5.618
54.77126.297
1.941
0.421
5.735
2.663
71.82055.96614.78920.97226.14424.57533.49656.35153.65886.748 210.930(1)ｄd形式4倍精度はphi5110Pは高速(2)dｄ形式6倍精度,8倍精度,１０倍精度はSR16000 32coreとphi5110Pの性能が同等。
(3)Cとソフトウエアサポートルーチンは
E5-2670が高速。

行列サイズN=480の結果は以下の様になっています。
行列積N=480,1000回実行の性能比較
実行時間(秒)プログラム
multd
ddmultd
multq
doubleexdouble
ddexdouble
mult6
mult8
mult10

SR16000
32core 64smp
55.214164.514289.634605.880

phi5110P E5-2670
240smp
16smp
3.445
3.204
29.658
11.21244.669469.772 209.251
3.728
3.196
31.24621.076256.872 191.585
119.144
159.443 205.080
196.611285.541 460.643
436.391688.837 1686.164

性能比較の傾向はN=240の場合と全く同じ。

各種行列積計算の性能は以下の様になっています。
各種行列積計算性能比較表
実行時間一覧表(秒)NCPU
240 E5-2670
Phi5110P
480 E5-2670
Phi5110P
720 E5-2670
Phi5110P
960 E5-2670
Phi5110P
1200 E5-2670
Phi5110P

double exdouble ddexdouble multd
ddmultd
0.421
2.663
55.966
0.403
5.618
1.941
5.735
71.8201.824
2.952
3.196 21.076 191.585
3.204 44.669
3.728 31.246 256.872
3.445 11.212
17.364 70.782 485.278 19.230 150.859
11.141 103.516 948.543 10.315 35.615
40.647 204.819 1035.046 45.322 357.364
25.932 253.483 2284.650 25.001 116.462
79.576 369.371 1890.814 88.222 720.564
56.584 564.900 4642.316 57.053 241.207(1)倍精度ではc,FORTRANともN=720でPhi511
0P の性能がE5-2670を上回ります。
(2)拡張倍精度はすべてE5-2670の性能が上回っています。
(3)SR16000の4倍精度では,Phi5110Pが非常に効力を発揮する事がわかります。

これに関連した計算は以下の様になっています。
4倍精度行列積計算

N=256

1回の実行時間

SR16000 singleジョブ
92.638 msec
ieeequad 8704.767 msec
x5570
double+double
457.259 msec
extend double +extend double 734.146 msec
e5430 gcc
double+double
652.392 msec
extend double +extend double 759.277 msec

e5430 icc
double+double
309.527 msec
extend double +extend double 621.636 msec

Phi5110P 240smpがSR16000 １ノードより高速という結果と合わせPhi5110Pが有効な事がわかります。

6倍精度,8倍精度,10倍精度行列積計算
SR16000/M1
1ノード実行
コンパイルオプションは-W0,’opt(o(s),disbracket(0))‘
opt(o(ss))にしたり、disbracket(0)を外すと結果不正となります。
行列積テスト結果一覧表
N=100,1000回実行時間(秒)精度
single32core64smp
6倍精度
46.315
1.877
1.304
8倍精度
80.729
3.285
2.14
10倍精度
188.7617.58
5.275(1)並列化効果はよくでている。
(2)8倍精度はSIMDがよく効いているので他の２つの精度に比べ演算数比率より性能が良い。
（8倍精度は４つの倍精度変数を扱うため。）

4.3.2 ストラッセンの行列積
代表的な例でその計算方法と演算量を示しました。
C  AB A 11 A 12 
 B11 B12  C11 C12 




A
|
,
B

|
,
C

A

B

C | C 

A
B
22 
22 
22 
 21
 21
 21
P1  ( A 11  A 22 )(B11  B 22 )P2  ( A 21  A 22 ) B11
P3  A 11 ( B12  B 22 )P4  A 22 ( B 21  B11 )P5  ( A 11  A 12 ) B 22
P6  ( A 21  A 11 )(B11  B12 )P7  ( A 12  A 22 )(B 21  B 22 )C11  P1  P4  P5  P7
C12  P3  P5
C 21  P2  P4
C 22  P1  P3  P2  P6
m段目の行列積の回数 M m , 行列和の回数 A m
M 0  1, A 0  0
M m  7 M m 1 , A m  18M m 1  A m 1
M m  7 m , A m  3(7 m  1)
もとの行列のサイズ N, m段目は
演算量

N
2m

通常の行列 積 2 N 3

73(7 m  1) 2
2( ) m N 3 
N
8
4m
N
8m N
8
演算量比 
≒≒ ( )mm
m
m 1
7
7
3 (7  1) 2
7 (N  3  2)( )m N3 
N
m
8
2
4
7
strassenの行列積の演算量は通 常の行列積の ( ) m
8となり, 数学的には , strassenの行列積の演算量はstrassen

N 2.807 に比例すると言えます 。

演算量に関してstrassenの行列積 性能一覧表
SR16000/M1
ーＯｓｓ Ｎ＝６５５３６ smt on 64smp/32core

プログラムの性能(GFLOPs)
演算量2  N 3 (FLOP) / 実行時間(秒)
段数

プログラム モニター
実行時間
(GFLOPs) (GFLOPs) (秒）0304.764305.839 1847.162
1
343.957
302.132 1636.684
2
377.379
290.097 1491.735
3
407.541274.253 1381.331
4
454.6
267.3811238.345
556.783
287.836 1011.074
6461.391
210.794 1220.112
7300.128121.916 1875.696

見かけ上の性能値と論理最大性能
ストラッセン行列積性能一覧表
ＳＲ１６０００／Ｍ１

smt=off

測定値(GFLOPs)
演算量2  N 3 (FLOP) / 実行時間(秒)
Ｎ

段数
38400
38400
46080
46080
61440
61440
61440
61440
61440
61440
61440
61440
61440

smp
4
4
3
3
4
4
4
4
4
4
4
5
5

測定値
論理最大性能(GFLOPs) (GFLOPs)
8 254.12
245.1232 536.70
980.488 255.76
245.1232 726.88
980.488 307.47
245.12
12 408.76
367.6816 589.95
490.2420 586.77
612.8024 600.06
735.3628 578.24
857.9232 1027.31
980.488 292.72
245.1216 510.21
490.24

4.実アプリケーション
4.1 重力多体問題(N体問題)
重力多体問題は
mj 
Fi  Gmi  3 rij
j  i rij 
rij | ri  rj |G : 万有引力定数, mi : 粒子iの質量ri : 粒子iの位置, Fi : 粒子iにかかる力で計算します。
演算精度は倍精度でSIMD命令が適用されます。
演算量は粒子数Nの二乗に比例します。

4.1.1 ＳＲ１６０００/M1 システム
４ノードでフラットMPIでSMT ON、SMT OFFの場合の結果です。
演算量がN=2**nの二乗に比例しますので、実行時間はlog2スケールで表しています。タイムステップ数は100です。
Log2(実行時間（秒））
14

7103

12
1778

10
446

8

113

smt on6smt off

28
4

7
2
0
15

16

17

18

19

20

n

粒子数N=2**n

SMT ON,OFF での実行時間に有意差は見られません。

また、フラットMPIとハイブリッドの差をN=2**18,2**20で調べた結果は以下の様になりました。

サイズ2**n smp数
task数
ノード数 実行時間（秒）
20
1
256
8
3599
20
1
512
8
3612
20
32
8
8
3603
20
64
8
8
3604
20
8
32
8
3593
20
8
64
8
3607
サイズ2**n smp数
task数
ノード数 実行時間（秒）
18
1
256
8
224
18
1
512
8
231
18
32
8
8
222
18
64
8
8
228
18
8
32
8
222
18
8
64
8
223

フラットMPIとハイブリッドでも実行時間に有意差は見られません。

4.1.2 ＢＧ／Ｑ システム
領域分割をせずに並列実行する場合、粒子数N=2**15=32768
からN=2**20=1048576 とすると1スレッド当たりの所要メモリ量
が3.328MBから104MBため、SMPとMPIのハイブリッド方式がフラットMPIより効果がでます。演算量がNの二乗に比例する
ため、実行時間の軸はlog2を取っています。実行結果は以下の様になりました。タイムステップ数は１００です。

Log2(実行時間(秒))
12
2179

10

1922

595
482

184

8

121

58

6
21.1

フラットMPI
ハイブリッド

31

4
7.9

6.1

2

2.1

0
15

16

17

18

19

粒子数N=2**n20n4.1.3 各種計算機性能比較
使用した計算機は以下のものです。
計算機
A
計算機
B
AC
D
BE
CF
G
D

BG/Q 128node smp
BG/Q 128node 8192MPI
BG/Q
BG/Q 32node
128node2048MPI
smp
SR16000/M1 4node 128MPI
BG/Q 128node
8192MPI
SR16000/M1
4node
256MPI
HD58703200core
BG/Q 32node
2048MPI
HD69706144core
SR16000/M1 4node 128MPI

E
HD6970
FHD5870

SR16000/M1 4node 256MPI
880MHz1536core 683GFLOPs 176GB/sec
HD58703200core 544GFLOPs 153.6GB/sec
850MHz1600coreG

HD69706144core

HD6970880MHz1536core

683GFLOPs 176GB/sec

HD5870850MHz1600core

544GFLOPs 153.6GB/sec

扱ったケースは,演算量は可変で粒子数を変化させた場合と演算量は一定で粒子数を変化させた場合を扱っています。

4.1.3.1 演算量可変ケース
重力多体問題実行時間
log2(実行時間(秒))

16
14

実行性能
BG/Q 32ノード≒SR16000/M1 4ノード
≒HD6970

12
A

10

B
C

8

D

E

6

F
G

4
2

n

0
15

16

17
粒子数

18

19

Ｎ＝2**n20

5.1.3.2 演算量一定ケース

粒子数N  1000,4000,10000でタイムステップ数はそれぞれ,10000,
625,100として演算量を同じにして(平方根計算あり270GFLOP, なし580GFLOP)
性能を比較.演算は倍精度演算。
N  1000の場合, 並列化オーバーヘッド により
SR16000の性能が勝っているが,N  4000,10000ではGPUの性能が勝っています。
またieee754-2008系の特徴でもある
拡張倍精度演算に関する結果も扱っています。

N体問題実行結果一覧表(秒)粒子数=1000,タイムステップ=10000
演算量(ソース)=500GFLOP
SR16000

1ノード

32core64smp
4.73322.6379

E5-2670

最適化
4.7668 OPT3
2.6471 SOPT

16smp

16.0409

各種加速装置実行時間一覧表
GPU
HD5870HD6970HD7970
W8000
HD7980Phi5110P

1ボード
２ボード
３ボード
4ボード
(1smp/core) (2smp/core) (3smp/core) (4smp/core)
24.0556
32.9047
30.434242.083952.272964.291141.0905
41.061447.1802
4.21373.18912.82922.6806

N体問題実行結果一覧表(秒)粒子数=4000,タイムステップ=625
演算量(ソース)=500GFLOP
SR16000

1ノード

32core64smp
4.25972.6304E5-2670

最適化
4.3365 OPT3
2.7826 SOPT

16smp

4.6922

各種加速装置実行時間一覧表
GPU
HD5870HD6970HD7970
W8000
HD7980Phi5110P

1ボード
２ボード
３ボード
4ボード
(1smp/core) (2smp/core) (3smp/core) (4smp/core)
4.61265.24075.25416.0501
6.86957.72157.89258.04388.78003.1736
1.92311.6821
1.3974

N体問題実行結果一覧表(秒)粒子数N=10000,タイムステップ=100
演算量(ソース)=500GFLOP
SR16000

1ノード

32core64smp
4.16372.5880E5-2670

最適化
4.5797 OPT3
2.6867 SOPT

16smp

4.5523

各種加速装置実行時間一覧表
GPU
HD5870HD6970HD7970
W8000
HD7980Phi5110P

1ボード
２ボード
３ボード
4ボード
(1smp/core) (2smp/core) (3smp/core) (4smp/core)
2.74801.99173.0483
2.18852.37252.59313.00593.0796
3.24653.17361.86531.56181.3542拡張倍精度演算

拡張倍精度実行時間一覧表(秒)Ｎ体問題 倍精度実行時間一覧表(秒)N

反復回数 E5-2670 Phi5110P E5-2660 E5-2660
16smp 240smp 16smp 32smp
1000 10000 15.8442 32.0745 20.6893 22.7324
4000 625 10.1019 28.8471 18.4697 14.5187
10000 100 10.0968 28.5530 11.1613 15.1299

倍精度の場合とは逆にE5-2670の性能がE5-2660の性能を上回っています。

拡張倍精度に対する処理でE5-2670のコンパイラがe5430のコンパイラより非常に適している事によります。

4.2 分子動力学計算
分子動力学計算は以下 の２つの計算からなり 、
クーロン力計算は重力 多体問題と同様の計算 となりますが、異なる計 算として、
ワンデルワールス力計 算があります。逆二乗 則と逆六乗則の差により、 演算する範囲が異なり 、演算量も単純に粒子数 の何乗に比例するとは言えません。
演算精度は倍精度で SIMD命令は一部適用されま す。

qj 
qi
 
クーロン力 Fi rr
|ri  rj |ij
3 ij
4 0 ji rij
q i : 電荷量 ,  0 : 真空の誘電率
ワンデルワールス 力
分子間に働 く分散力で、等方向性 で原子間距離の６乗に 反比例する力
ab
F  k 6r

 : 分極率

4.2.1 SR16000/M1 システム詳細
SR16000/M1 複数ノードでの実行時間は以下の様になっています。
実行はフラットMPIを使用しています。所要メモリはそれぞれ、
11.8MB,28.0MB,54.7MBとなっています。
実行時間（単位秒）
N=48
MPI数
ノード数
VdW
Coulomb
全体
128
4
1.81984.44435
6.26732256
42.52998
8.61865
11.15244256
81.718892.66658
4.38859512
83.35877
3.98576
7.34833
N=64
MPI数
ノード数
VdW
Coulomb
全体
128
4 12.54277
63.21624
75.76465
256
4 16.55753
97.01006 113.57417
256
88.7749335.4702244.25077512
8 20.90441
63.69971
84.61077
N=80
MPI数
ノード数
VdW
Coulomb
全体
128
4 44.93113 412.63567 457.57689
256
4 63.25216 362.52039 425.78437
256
8 33.06468 259.85763 292.93251
512
8 91.30541 281.84471 373.16153

粒子の数はn=N**3 で演算量はCoulmb(クーロン力)はn(n-1)/2に比例,VdW(ファンデルワールス力)はn(n-1)/2/
近傍にある粒子の数に比例します。
VdWではすべてのケースでSMT OFF の効果があり、クーロン力
はN=80,4node ではSMT ONがそれ以外はSMT OFFが効果があります。これは転送量がMPI数に比例するため、演算量との比でSMT ON,OFFの効果が異なる様になります。

4.2.2 各種計算機性能比較
N
48
64
80
96
112
128

所要メモリ
(MB)
11.812528
54.6875
94.5
150.625224
実行時間（秒)一覧表
N=48

SR16000/M1
MPI数
ノード数
256
256
512
BG/Q
MPI数

VdW

ノード数

512
1024
2048

2.52998
1.71889
3.35877

Coulomb
全体
8.61865
11.15244
2.66658
4.38859
3.98576
7.34833

8.16086
11.82922
24.60092

Coulomb
全体
8.8191616.98838
9.9951121.83329
17.70815
42.32042

4
8
8

VdW
32
32
32

N=64
SR16000/M1
MPI数
ノード数
128
256
256
512
BG/Q
MPI数

4
4
8
8

VdW
Coulomb
全体
12.54277
63.21624
75.76465
16.55753
97.01006
113.57417
8.7749335.4702244.25077
20.90441
63.69971
84.61077

32
32
32

VdW
Coulomb
全体
29.56434
37.76663
67.34834
32.48871
32.148564.65611
57.66445.01047
102.69909

ノード数

512
1024
2048

N=80
SR16000/M1
MPI数
ノード数
128
256
256
512
BG/Q
MPI数
512
1024
2048

4
4
8
8

VdW
Coulomb
全体
44.93113
412.63567
457.57689
63.25216
362.52039
425.78437
33.06468
259.85763
292.93251
91.30541
281.84471
373.16153

32
32
32

VdW
Coulomb
全体
92.13976
128.87558
221.04942
80.97082
91.95302
172.96071
121.47084
106.98893
228.50708

ノード数

メモリの制限により,N=112,N=128ではBG/Q 2048MPIは実行不可
N=96
SR16000/M1
全体
Coulomb
VdW
ノード数
MPI数
4 135.30229 1018.25815 1153.57737
128
4 170.37763 1035.44061 1205.83632
256
96.82403 592.82008 689.65997
8
256
906.2826
8 217.66733 688.59726
512
BG/Q
MPI数

ノード数

512
1024
2048

全体
Coulomb
VdW
32 247.50348 392.02168 639.58082
32 185.35274 274.96719 460.38124
32 231.82941 340.60404 572.51359
N=112

SR16000/M1
全体
Coulomb
VdW
ノード数
MPI数
4 333.35061 2736.97492 3070.34942
128
4 411.22115 2556.15124 2967.40039
256
8 244.40453 1447.09706 1691.52595
256
8 369.37199 1905.33538 2274.73516
512
BG/Q
MPI数

ノード数
32
32

512
1024

全体
Coulomb
VdW
591.7365 1262.05753 1853.88598
394.3729 1181.64034 1576.11295
N=128

SR16000/M1
全体
Coulomb
VdW
ノード数
MPI数
755.4252 5581.71674 6337.17718
4128
4 924.26017 6382.94117 7307.24141
256
539.1269 3196.45186 3735.61507
8256
8 750.14716 4215.49161 4965.67888
512
BG/Q
MPI数
512
1024

ノード数

全体
Coulomb
VdW
32 1317.87477 3370.73599 4688.74093
32 805.80939 3364.92311 4170.8758

BG/Q
mdcoresmp数

mdcore2
smp数

ノード数
VdW
Coulomb
全体
16
32
219.958
326.51232
32
133.103
198.09948
32
113.692
166.06764
32
110.894158.311対smt=1546.518331.25279.81
269.255

1
1.65
1.95
2.03

ノード数
VdW
Coulomb全体
対smt=116
32 232.7985
361.71276
594.61322
1
32
32 119.5483
186.49663
306.14688
1.94
48
32
80.8
140.882221.784
2.68
64
32 61.61527
109.81348
171.53069
3.47

xm1
mdcore
smp数

ノード数
32
64
96
128

mdcore2
smp数

1
1
1
1

VdW
Coulomb全体
対smt=1
370.355
1031.531
1402.0251
246.078
611.013857.267
1.64
223.676446.46670.3282.09
173.189
340.836514.222
2.73

1
1
1
1

VdW
Coulomb全体
対smt=1353.256930.879
1254.2761
306.571605.227911.972
1.38
229.255
442.445671.891.87
177.715
339.45517.357
2.42

ノード数
32
64
96
128

ＭＰＩオーバーヘッド削減の効果はBG/Qででています。
Xm1では１ノード並列のためほとんど変化はありません。

smt=1に対する効果一覧表
プログラム
mdcoremdcore2
mdcoremocore2

ｃpu
BG/Q
BG/Q
xm1
xm1

smt=2

smt=3
1.65
1.94
1.64
1.38

smt=4
1.95
2.68
2.09
1.87

2.03
3.47
2.73
2.42

xm1
論理コア128(new),64(old)の比較
mdcoresmp数
32 new
32 old64 new
64 old128 new
mdcore2smp数
32 new
32 old64 new
64 old128 new

ノード数

VdW
Coulomb
全体
1 370.355
1031.531
1402.025
1 370.239
1051.356
1421.733
1 246.078
611.013857.267
1 250.312
554.729805.976
1 173.189
340.836514.222

ノード数

VdW
Coulomb
全体
1 353.256
930.879
1254.276
1 352.106
931.449
1283.693
1 306.571
605.227911.972
1 297.089
556.781854.376
1 177.715
339.45517.357

BG/Q MPIオーバーヘッド削減の効果が大きい(smt=4)
Xm1 論理コア128(smt=4)の拡大の効果が大きい。

4.3 大規模疎行列の反復計算
4.3.1 対称問題
倍精度演算
ポ アソン方 程 式
u  f (), u  0 ()  [0,1]  [0,1]  [0,1]N  200  200  200
収 束 判 定 値: 共 役 残 差10 1 2
u  sin( x ) sin( y) sin( z )f  3   2  sin( x ) sin( y) sin( z )（注）
初期値は x (i)  1.0
cgs, cgsilu, cgsmilは収束しない
ので x (i) : (0,1)の一様乱数値 .
cgsmは cgsiluを, cgsmmは cgsmilを並列化して初期値
x (i)  1.0としたもの。

対称行列反復法
実行時間一覧表(秒)プログラム x5570
E5-2670 Phi5110P E5-2660 E5-2660 SR16000 SR16000
1smp
16smp
240smp
16smp
32smp 32smp 64smp
d3sor3305.8545 107.2121
45.4492
81.5021 85.1895 18.4587 15.9792
d3sor25157.8583 2106.8935 845.0379 1523.1281 1658.299 285.6998 300.2200
d3adi6484.9054
976.999 732.4666 774.8730 1005.332 291.1860 281.3099
d3bcg
137.1272
54.8386
19.8924
41.5801 38.4800 12.9068 13.1204
d3cg
90.8600
31.1266
11.2524
20.9102 19.2947 7.8468 7.7593
d3cgs
127.0035
45.2665
18.884737.1193 34.9020 14.0279 15.7210
d3scg
89.9558
36.7480
13.218124.6104 23.9617 9.1910 11.4737
d3bicgs
112.0541
31.3391
14.602025.8810 24.5952 10.2689 11.7728
d3cgsilu
148.3735
92.7729 572.2768
92.6816 93.9379 60.3351 61.4112
d3cgsmil 118.5756
82.1866 511.6339
81.9387 82.9746 53.7503 53.3771
d3gpbicg 122.8012
46.647118.274530.8516 33.7076 13.3188 12.8016
d3cgsm456.177748.207843.7151
45.3634 42.8455 27.0258 27.2273
d3cgsmm 487.1845
58.6426
77.231573.7438 62.3460 34.3474 42.1951

E5-2670,E5-2660,Phi5110PはSR16000と比較して実行効率ではほぼ同等の性能がでています。

4.3.2 非対称問題

演算は4倍精度演算
問題は3次元ポアソン方程式
uu  R fx
領域[0,1]  [0,1]  [0,1]解析解u ( x , y, z)  e xyz  sin( x )  sin( y)  sin( z)
非対称問題。
解法
bcg ( biconjugat e

gradient ) 法

R  100, nx  ny  nz  129
cgs (conjugate

gradient

square)法R  100, nx  ny  nz  129
収束判定値は共役残差10 12
初期値u ( x , y, z)  0

非対称問題実行時間一覧表(秒)プログラム E5-2670 Phi5110P E5-2660 E5-2660 SR16000 SR16000
16smp 240smp 16smp 32smp 32smp 64smp
bcg
47.8155 92.7816 75.5405 54.6095 5.6664 4.1678
cgs
50.6370 99.7317 77.4960 59.7968 9.3170 7.4054
4倍精度演算ではSR16000の性能がカタログ
性能比を考慮しても非常に良い性能を示しています。
今回のｃｇｓ法は前処理付きだが以前の方式
ではサイズが113*113*113で6倍精度
SR1 6000 64smp 29.1450秒
E5-2670 16smp 47.0610秒
8倍精度
SR1 6000 64smp 33.5245秒
E5-2670 16smp 70.0802秒と倍精度変数をつなげた形式ではE5-2670の実行効率が勝っていました。

4.4 QCD計算
テストケース
case1

data 16  16  16  16

case2

data

24  24  24  24

case3

data

32  32  32  32

測定条件
BG / Q, SR16000 / M 1, xm1
smt数

すべて1ノード

BG / Q, xm1

smt  1,2,4

SR16000 / M 1

smt  1,2

スレッド並列のソース
BG/Q

ppgenで作成したソース。-qsmp=auto
でコンパイル.
SR16000/M1,xm1 ともに自動並列でコンパイル。

QCD 性能測定結果一覧表（対1スレッドに対する台数効果)
data

BG/Q
SR16000/M1 xm1
smt=1 smt=2 smt=4 smt=1 smt=2 最高性能smp数
case1
13.63 15.84 11.22 18.61 14.85
32case2
12.85 13.28 12.75 9.81 11.18
32
case3
12.33 13.03 12.88 11.09 9.48
32
最高性能
BG/Q
smt=2 =>32smp
SR16000/M1 smt=1 =>32smp
xm1=>32smp扱うデータがおおきくなるとSR16000/M1
は大きな性能低下が見られます。
BG/Qはデータサイズによる性能低下は見られません。
BG/Q smt=1 (16smp)で台数効果13と高い。

性能比較
QCD 32*32*32*64 並列実行性能測定結果
実行時間一覧表(秒)1 ノードsmp BG/Q ソースはParallel Program Generetor で作成
SR16000
339.15 BG/Q
SR16000/BG実行性能比 6.55
SR16000/BG カタログ性能比 4.79

2221.9mpi,ｓｍｐ並列
SR16000
node
total
1
4
8

BG/Q
node
217.4
46.37
18.8

total
32

56.34

SR16000 (対BG/Q)
node
実行性能比
カタログ性能比
1
0.259
0.15
4
1.215
0.598
8
3.029
1.197

ソースチューニングがない場合,SR16000/M1
がBG/Qに対してカタログ性能比以上に実行性能比がでています。

4.5 量子モンテカルロ法による物性スペクトル計算
プログラムは高エネルギー加速器研究機構,
物構研岩野氏より提供いただいたのをベースにしています。
数値計算からみた場合の計算内容
G i , A i , B, C i : n  n行列
A m  C m C m 1 .......C1
A0  I(単位行列)B  C L C L 1 .........C 1
G m  A m(I  B)1
パラメータ : n, L,β, u
β0  m  L, dt L
C i の絶対値の最大値, 最小値

e dtu , e  dtu

Tr(G m )を計算し, 結果のチェックはTr(G L ) 1で行います。

最初の実行結果(4倍精度)
n=40,β=10,L=112,u=5 結果は正確
n=100,β=10,L=448,u=5 結果は不正確
目標 n=100,L=448を固定
β=10,20 u=5,6,7,8,9,10での結果を10進10桁一致させる。

使用されている計算手法
行列のＱＤＲ分解。
行列VをV=QR (Q:ユニタリ行列,R:上三角行列)に分解後,Rの対角要素をDに入れ,Rの対角要素を1にする。
扱ったケース：設定したパラメータより結果の精度が見積もりやすい。
P  e dtu L , E(m)| Tr(G m )|プログラム実行中に現れる数値の最大値P 2 ,0以外の絶対値最小値1/P 2
行列やベクトル要素に現れる数値の最大値P,0以外の絶対値最小値1/P
LL
1でE( )
で最小になるので,
22P1演算に必要な最小ビット数log 2 P  log 2 P。
2E(0) 1, E(L) 1, m 

精度検証結果(10進10桁一致したもの)
6倍精度 β=10,u=5,6,7,8,9
8倍精度 β=10,u=5,6,7,8,9,10
β=20,u=5,6,7,8
10倍精度 β=20,u=9,u=10

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

性能測定結果
量子モンテカルロ実行時間(秒)SR16000/M1 1ノード
4倍精度
smt on,off 実行時間
32core
9.327
64smp
11.04

１０倍精度
QDR積計算
64smp
32core
並列化不可
733
783
分割並列化
426
456
数値制限解除
451
435

精度
6倍精度
8倍精度

uinlineオプション
サブルーチン化 64smp
32core
612
263
1190
463

精度
6倍精度
8倍精度

行列積,QDR積計算
64smp
32core
112
192

243
386

ＩＦ文削除
64smp
32core10523.74328.837174
42.07451.73910倍精度演算の問題

加減算,乗算での演算数の増加(1)最適化,uinlineの効果が減少(2)行列積計算での並列化効果の減少(3)QDR積計算が並列化不可となり,2つのサブルーチンに分割した事による並列化効果の減少表現できる数値範囲の制限
大きな数の平方根計算では、アルゴリズムの変更が必要。

4.5.2 ieee754-2008形式

E5-2670

16smp

Phi5110P

240coreモンテカルロシュミレーション
精度
6倍精度
8倍精度
10倍精度
ieee8倍精度
DQ(8倍精度)

U

β
10
20
20
20
20

L
6
8
10
10
9

E5-2670 Phi5110P
実行時間 実行時間(秒)(秒)
448 127.5865 1088.6992
448 261.7628 1762.0854
448 619.8694 3727.0650
448 505.0661 4487.2604
448 1961.2144 16287.9011(1)β=20,U=10を計算する場合dd形式では10倍精度演算が必要でieee形式では8倍精度で事足ります。
(2)SR16000でscopeオプションで実行すると10倍精度は649秒.このため,計算機種の差ではなく,コンパイラの最適化能力の差が効いていると言えます。

DD形式8倍精度とieee754-2008形式の8倍精度演算では,正しく計算出来るパラメータ領域に差がある。
一概にどちらが良いとは言えない

β=20,U=8ではDD形式8倍精度が良い.
β=20,U=10ではDD形式では10倍精度を使用する必要があり,ieee754-2008形式の8倍精度が良い。

10倍精度演算ではieee754-2008形式の方が性能が良いと考えられ,GPU等のアクセラレータで8倍精度演算までのDD形式を適用している理由とも合う結果となっている。

4.6 ファインマンループ積分
ファインマンループ積分の一般式は以下の様になる。
N  n ( L 1)

2
1
N
1 nL
C
I()  (1) N ( ) 2 ( N  nL )   dx i (1  x 1  x 2  .....  x N )N  nL2 0 i1
4
(D  iC) 2

n : 時空次元 (n  4), N : ループ内の素粒子の数 , L : ループの多重度
 : 実定数 , C, D : x 1 , x 2 ,..., x N の多項式

4.6.1 inf ra

vtx

計算

1 1 x

1
dydx
D
0 0

I 

D   xys  ( x  y) 2 m 2  (1  x  y)2s  500 2 , m  0.0005
倍精度演算で   10 5 で実施しています。
実行条件
HD5870 3200SMP
E5  2670 16SMP
Phi5110P 240SMPGFLOPs

2次元積分(VTX)性能測定結果（演算量はソース上から算出)

400.000

350.000
300.000
250.000
200.000

HD5870E5-2670

150.000

Phi5110P

100.000
50.000
0.000
N

サイズ

4.6.2 Infra box 計算
計算式
1 1 x 1 x  y

1
I     2 dzdydx
0 0 0 DD  sxy  tz (1  x  y  z )  ( x  y)2  (1  x  y  z )(1  x  y)me z (1  x  y)m f

2

テストデータs  500 , t  150 , m  150,
22fm  0.0005e  10 (倍精度),   10 (4倍精度)
12

30

2

倍精度

N=8192,演算量10447GFLOP

SR16000/M1 １ノード 自動並列
BG/Q
３２ノード ノード間
ノード内

MPI
スレッド並列計算機
smp数 smt数
性能(GFLOPs) 実行効率(%)
SR16000/M1
32
1
268
27
SR16000/M1
64
2
263
27
BG/Q
16
1
846
13
BG/Q
32
2
1493
23
BG/Q
64
4
1766
27
HD5870
1
3200
326
30
HD6970
1
6144
588
22
実行効率は２５％－３０％

GFLOPs

3次元積分(BOX)性能測定結果（演算量はソース上から算出)

600.0000

500.0000

400.0000
HD5870300.0000

E5-2670
Phi5110P

200.0000

100.0000

0.0000400 120020002800360044005200
サイズ

N

４倍精度
サイズ
N=2048
テーブルサイズ 32KB

演算量

618GFLOP テーブルサイズ 32KB

SR16000/M1 4,８ノード フラットMPI
BG/Q
32ノード フラットMPI
４倍精度性能測定結果(GFLOPs)
計算機
ノード数 smt数
性能
SR16000/M1
4
1
SR16000/M1
4
2
SR16000/M1
8
1
SR16000/M1
8
2
BG/Q
32
1
BG/Q
32
2
BG/Q
32
4
HD58702ボード
3200
HD6970
4ボード
6144

実行効率(%)
248
6.3
418
10.7
495
6.3
836
10.7
253
3.9
475
7.2
705
10.8
283
26
560
20.7

SR16000/M1,BG/QではＳＩＭＤが適用されないので,GPGPU系の計算機の実行効率が良い

１ ノード性能比較

サイズN=１０２４

実行時間（秒）
3530

BG/Q(換算) 64smp≒SR16000/M1 64smp

25

20

BG/Q
BG/Q(換算）

15

SR16000/M1

10
5

0
1

2

3

4

SMT数

BG/Q
SR16000/M1
SMT=1 16 smp
32 smp=2 32 smp
64 smp=3 48 smp
96 smp=4 64 smp
128 smpBG/Q(換算) 204.8GFLOPs => 980.48GFLOPs にしたときの値

サイズN=2048
実行時間（秒）300
250

BG/Q(換算) 64smp≒SR16000/M1 64smp
200
BG/Q

150

BG/Q(換算）
SR16000/M1

100
50
0
1

2

3

4

SMT数

BG/Q
SR16000/M1
SMT=1 16 smp
32 smp=2 32 smp
64 smp=3 48 smp
96 smp=4 64 smp
128 smpBG/Q(換算) 204.8GFLOPs => 980.48GFLOPs にしたときの値

サイズN=4096
実行時間（秒）
2500
2000
BG/Q(換算) 64smp≒SR16000/M1 64smp
1500

BG/Q
BG/Q(換算）

1000

SR16000/M1

500
0
1

2

3

4

SMT数

BG/Q
SR16000/M1
SMT=1 16 smp
32 smp=2 32 smp
64 smp=3 48 smp
96 smp=4 64 smp
128 smpBG/Q(換算) 204.8GFLOPs => 980.48GFLOPs にしたときの値

４倍精度INFRA BOX N=1024 各種計算機実行時間（１）実行時間（秒）8

論理最大性能BG/Q 6553.6 GFLOPs
SR11000 537.6GFLOPs
GPGPU 1088GFLOPs
XM1 844.8GFLOPs
T2K 64MPI 588GFLOPs
T2K 128MPI 1176GFLOPs
T2K 256MPI 2352GFLOPs
BG/L 5734.4GFLOPs

7
6
5

4
3
2

実行時間（秒）1
BG/L 1024node 1024MPIT2K 256MPIT2K 128MPIT2K 64MPI

XM1 4MPI*16SMP

XM1 64SMP

GPGPU

SR11000 64PE

BG/Q 32node 1024MPI

0４倍精度INFRA BOX N=1024 各種計算機実行時間（２）実行時間（秒）1.2
10.8

論理最大性能BG/Q 6553.6 GFLOPs
SR16000/M1
4node 3921.92GFLOPs
8node 7843.84GFLOPs

0.6
0.4

0.2
0

実行時間（秒）

9
8
7
6
5
4
3
2
1
0

SR16000/M1 8node 512MPI

SR16000/M1 8node 256MPI

SR16000/M1 4node 256MPI

SR16000/M1 4node 128MPI

GPGPU

BG/L 1024node 2048MPI

BG/Q 32node 2048MPI

BG/Q 32node 1024MPI

BG/Q 32node 512MPI

４倍精度INFRA BOX N=2048 各種計算機実行時間

実行時間（秒）

論理最大性能BG/Q 6553.6 GFLOPs
SR16000/M1
4node 3921.92GFLOPs
8node 7843.84GFLOPs
GPGPU 1088GFLOPs
BG/L 5734.4 GFLOPs

実行時間（秒）４倍精度INFRA BOX N=2048とN=4096
各種計算機実行時間比率

演算量比率(N=4096/N=2048)=8

実行時間（N=4096/N=2048)比率
12
10
8
6
4
2
0

倍率

4.6.3 massless計算

1 1 x1 x1 ........ x n 1

    ........
n

0

0


0

1
dx n dx n 1 .....dx 1
1 ( x 1 x 2 .....x n ) (1  n)(2  n)....(n  n) 
1 1 x1 x1 ........ x n 1

  ........
0

0


(
x
x
.....
x
)
dx n dx n 1 .....dx 1
 1 2
n
0

となり, 倍精度演算で精度良く 計算するには ,(1  n)(2  n)....(n  n) 
1 1 x1 x1 ........ x n 1

  ........
0

0


(
x
x
.....
x
)
dx n dx n 1 .....dx 1
 1 2
n
0

を計算する必要があり ます。

倍精度演算による実行結果
条件一覧表
n
h
3 0.5^6=0.015625
4 0.5^6=0.015625
5 0.5^4=0.0625
6 0.5^3=0.125
7
8

N
503
503
103
51
41
33

0.16
0.20

実行時間一覧表(秒)smt=on,64smpnSR16000
3
4
5
6
7
8

n

#NAME?
xm1

0.148960
74.977050
16.623410
25.476964
362.011287
3199.456509

0.278902
112.659721
23.287427
38.277622
417.984915
3635.130558

Phi5110P 240smp
E5-2670 16smp
3
0.60690.17774
29.1766
81.6845547.393115.64066
20.8052
27.881471478.4449
265.0752
813738.5545
2255.2062

Phi5110Pでのチューニング例
1 1 x1 x1 ........ x n 1

  ........
0



0

0

1
dx n dx n 1 .....dx 1( x 1 x 2 .....x n )1

1 ((1  )) n n (1  n)
1 1 x1 x1 ........ x n 1

   ........
n

0

0


0

1
dx n dx n 1 .....dx 1( x 1 x 2 .....x n )1
1 1 x1 x1 ........ x n 1 (1  n)(2  n)....(n  n)   ........
00 ( x x .....x ) dx dx
1



2

n

n

n 1

.....dx 10において   0.01, n  7, n  8の場合を計算した。
倍精度で計算するため , 右辺の式を使用。
n  7, h  0.16, N  41
n  8, h  0.2, N  33
で実行した結果は下記 のとうり。

n=7,n=8 massless計算実行時間(秒)一覧表
SR16000 64smp,-Os 自動並列
XM1 64smp,-Os 自動並列
E5-2670 16smp -Os -openmp
Phi5110P 240smp -Os -openmpnSR16000 XM1
E5-2670 Phi5110P
6
25
38
21
27
7
362
418
265
1478
8
3199
3635
2255
13739

n=7,n=8の場合のPhi5110Pの値が悪すぎる。

n7

実行時間 1478秒  368秒

n 8

実行時間 13739秒  2623秒

n  7だと SR16000とほぼ同等 , n  8だと
SR16000の1.2倍の性能。
結果
n7
n 8

解析解

0 .996679645317184

実測結果

0.996679645317182

解析解

0 .9955954352742048

実測結果

0.995595435126120

N  41の場合の実測結果

0.995595435274190

結果は倍精度演算でも十分な精度となっている。
n=7とn=8の差はh=0.16(N=41),h=0.2(N=33)の違いによる。
一般にＤＥでは倍精度演算ではN=40-50で十分な精度が得られる。

このためx,yを以下の様に一次元化

do i1=1,nxx=x30(i1)by=1.0d0-xx
Do i2=1,nyy=x30(i2)*by
bz=1.0d0-xx-yy

前処理
l=0
do ip=1,nx1=x30(ip)b2=1.0d0-x1
do iq=1,nx2=x30(iq)*b2
l=l+1
xlist(l)=x1
ylist(i)=x2
wlist(l)=gw30(ip)*b2*gw30(iq)end doend dodo i1=1,n*nxx=xlist(i1)yy=ylist(i1)bz=1.0d-xx-yy
gw30(i1)*by*gw30(i2)wlist(i1)4.6.4 4-6次元積分計算
4.6.4.1 S221計算
1 1 x 1 x  y 1 x  y  z

S 221( s; m1 , m2 , m3 , m4 , m5 )   
2

2

2

2

2

0 0

 
0

0

1
dudzdydxDCC  ( x  y  z  u )(1  x  y  z  u )  ( x  y )( z  u )E  (1  x  y  z  u )( x  z )( y  u )  ( x  y ) zu  ( z  u ) xy
2

2

2

2

M 2  xm1  ym2  zm3  um4  (1  x  y  z  u )m52

D   sE  M 2C

s 221(1 : 100,100,0,0,100)
解析近似解 : 0.0380004438127

サイズ N=2048,演算量=281622GFLOP
テーブルサイズ 32KB
SR16000/M1 １ノード

BG/Q 32ノード

計算機
SR16000/M1 SR16000/M1 BG/Q BG/Q
並列化
smt=1
smt=2
flat
smp
実行時間(秒)
1872
1187
244
207
性能(GFLOPs)
150
237
1054
1360

S (5;100,100,0,0,100)
 1.2221

15

N 1024
HD5870HD6970GPU
122.551 sec
55.667 sec

HD5870HD6970HOST79903.86 sec
96984.31 sec

S 221 (1;100,100,0,0,100)の性能測定結果
サイズ N=576,演算量 1762GFLOP (性能モニター)
X5570,SR11000 自動並列,E5430 MPI
スレッド数 X5570
E5430
SR11000
実行時間 性能
実行時間 性能
実行時間 性能秒
GFLOPs 秒
GFLOPs 秒
GFLOPs
1 367.441
4.795 460.579
3.826 1477.16
1.193
2 204.562
8.614
232.52
7.578 738.884
2.385
498.03617.973 116.234
15.159 369.904
4.763
848.40236.40390.20619.533 185.167
9.516
1657.32630.73658.12830.31292.87318.97232
28.79
61.202
MPI16
175.644
10.032
MPI32
87.31520.18
MPI64
43.78240.245
BG/L MPI
スレッド数 N
256
512
1024
2048

演算量
実行時間 性能
GFLOP
秒
GFLOPs1024
17601 273.879
64.2661024
17601 136.939 128.532
1024
17601
68.469 257.065
2048
281622 557.638 505.027

S 221 (1;100,100,0,100)の性能測定結果

  算法 反復回数 15回
CPU
BG/L
X5570
E5430

N

演算量
スレッド数 実行時間 性能
GFLOP
秒
GFLOPs1024
264015
1024
1061
249
548
21660
8
609
36
576
26430
32
516
51

s 221(5;100,100,0,0,100)演算量  1.215 N  1024

17601GFLOP

実行時間（秒)一覧表
ボード数 HD5870
HD6970
gpu0
gpu1
1 242.313 243.377 221.059
2 122.455
110.715374.5194
55.638

HD5870 143GFLOPs,HD6970 316GFLOPs
SR16000/M1 システムに匹敵する性能がでています。
HD5870HD6970GPU
122.551 sec
55.667 sec

HD5870HD6970HOST79903.86 sec
96984.31 sec

GPU がHOSTに比べて非常に良い性能がでていますHOSTコンピュータでの性能が低い事の原因追究(1)GPU用に４重DOループを２重DOループにしたのが原因 ?(４次元積分,2次元積分)(2)使用コンパイラの問題では？
実行時間一覧表（秒）
cpu
x5570
x5570

コンパイラ
gcc
ifort

４次元積分 ２次元積分
60808
57308
4227
6255

cpu
e5430
e5430

コンパイラ
icc
icc parallel

４次元積分 ２次元積分
5918
12616
928
2785 8smp

ソースで二次元化したためキャシュミスが発生。
コンパイラの最適化処理能力。

ＧＰＵでのチューニングは、容量の大きいベクトル
レジスタを持ったベクトル機でのチューニングと同じ傾向を示します。

スーパーコンピュータでの４次元積分と２次元積分の比較
実行時間一覧表（秒）

CPU

64MPI

４次元積分

２次元積分

xm1

15.00185

25.03721

SR16000/M1

8.8044217.75754

T2K

11.88899

31.80749

４次元積分の方が性能が良いという傾向はサーバー系x5570,e5430と同じ。

ｆ90とxlf_r コンパイラの比較
S 221 (1;100,100,0,0,100)
演算量  47505

N  1152
GFLOP

SR16000/M1,xm1 １ノード

cpu
smp数
f90
xlf_r
xm1
32 247.354
133.49xm1
64 279.959
131.89xm1
96 274.067
133.65xm1
128 205.521
130.01
SR16000/M1
64 178.875
SR16000/M1
32 147.387
f90 xm1 実行効率27.4%,
SR16000/M1 32.9%と非常に良い値がでています。
xlf_r はf90以上に実行効率４３．３％とよく,キャシュの有効利用とＳＩＭＤが非常によく効いています。

SR16000,E5-2670,Phi5110Pの比較の結果は以下の様になっています。
4次元積分実行時間一覧表(秒)
case
case1
case1
case2
case2
case2

精度

SR16000
E5-2670 Phi5110P 演算量
32core64smp((GFLOP)
倍精度
11.698903 11.645035 27.618532 17.836520
2976
拡張倍精度
108.657508 257.571165
倍精度
24.126208 30.395169 26.958704 17.815714
4568
拡張倍精度
108.635537 256.691088
4倍精度
4538.965000 9762.636500(1)多重DOループの一重化はE5-2670,
Phi5110Pは若干速くなりますがSR16000は2~3倍程度遅くなっています。
(2)ieee754-2008形式の4倍精度は非常に遅くなっています。

サイズを大きくした場合とHD5870との比較は以下の様になっています。

４次元積分計算詳細テスト結果 (s221)
case2

N=1200精度

E5-2670
Phi5110P
16smp
240core506.127663 287.737872
2047.347488 4682.559193

倍精度
拡張倍精度
case2
精度
倍精度
拡張倍精度
4倍精度

実行時間(秒)

実行時間(秒)
HD5870Phi5110P
1面
2面
192core240core27.190329
13.814706
19.255160
17.815714
302.454400 256.691088
453.623999 232.079049

Phi5110Pは倍精度は高速でE5-2670は
拡張倍精度が高速。HD5870の4倍精度はdd形式のため高速です。

S 221 テスト
S 221 (  1,100 ,100 , 0 , 0 ,100 )N  576
演算量
倍精度
4倍精度

GFLOP

SR16000

実行時間(秒)

精度
倍精度
4倍精度

32core64smp
30.395169
24.126208
550.758631
330.429147

E5-2670

16SMP

4568
42187

精度
実行時間(秒)
倍精度
26.9587
拡張倍精度
108.6355
4倍精度
4008.4746
HD5870

実行時間(秒)

精度
倍精度
4倍精度

1面

2面
27.1903453.6240

13.8147232.0790

Phi5110P 実行時間(秒)精度
1smp /core
2 smp/core
3smp/core4smp/core
倍精度
42.0126
27.611623.3237
20.6132
拡張倍精度
335.2602
258.8425
255.4137
249.6453
4倍精度
3364.0728
1822.1579
1683.9084
1488.7508
HD6970 実行時間(秒)
精度
倍精度
4倍精度

一面
二面
三面
四面
11.443008
5.923182
4.036383
3.196938
322.827294
162.131639
107.290817
82.046807

グラフィックボード性能比較表
精度
倍精度
4倍精度

HD7970
W8000
HD79809.164966
9.455362
9.963973
267.937949
305.351604
304.181516

4.6.4.2 ５次元積分計算
laportad1 1 x1 1 x1  x2 1 x1  x2  x3 1 x1  x2  x3  x4

I  
0 0


0


0


0

1
dx5 dx4 dx3dx2 dx1
2Dx6  1  x1  x2  x3  x4  x5

D=&-x1**2*x2-x1**2*x3-x1**2*x4-x1**2*x6-x1*x2**2-x1*x2*x3&-2.d0*x1*x2*x4&-x1*x2*x5-x1*x2*x6-x1*x3**2-2.d0*x1*x3*x4-x1*x3*x5x1*x3*x6&-x1*x4**2&-x1*x4*x5-2.d0*x1*x4*x6-x1*x5*x6-x1*x6**2- x2**2*x4x2**2*x5&-x2*x3*x4&-x2*x3*x5-x2*x4**2-2.d0*x2*x4*x5-x2*x4*x6-x2*x5**2x2*x5*x6&-x3**2 *x4&-x3**2*x5-x3*x4**2-2.d0*x3*x4*x5-x3*x4*x6-x3*x5**2x3*x5*x6&-x4**2*x5&-x4**2*x6-x4*x5**2-3.d0*x4*x5*x6-x4*x6**2-x5**2*x6x5*x6**2

解析近似解=0.2762092253588

laportaf
1 1 x1 1 x1  x2 1 x1  x2  x3 1 x1  x2  x3  x4

I  
0 0


0


0


0

1
dx5 dx4 dx3dx2 dx1
2Dx6  1  x1  x2  x3  x4  x5

D=.-x1*x1*x2-x1*x1*x4-x1*x1*x5-x1*x1*x6-x1*x2*x2-x1*x2*x3.-x1*x2*x4.-2.0d0*x1*x2*x5-2.0d0*x1*x2*x6-x1*x3*x4-x1*x3*x5-x1*x3*x6.-x1*x4*x4-3.0d0*x1*x4*x5.-2.0d0*x1*x4*x6-x1*x5*x5-x1*x5*x6- x1*x6*x6-x2*x2*x3.-x2*x2*x5-x2*x2*x6.-x2*x3*x3-x2*x3*x4-x2*x3*x5-2.0d0*x2*x3*x6-x2*x4*x5.-x2*x4*x6-x2*x5*x5-x2*x5*x6-x2*x6*x6-x3*x3*x4-x3*x3*x5.-x3*x3*x6-x3*x4*x4-2.0d0*x3*x4*x5-2.0d0*x3*x4*x6x3*x5*x5-x3*x5*x6.-x3*x6*x6-x4*x4*x5-x4*x4*x6-x4*x5*x5-x4*x5*x6-x4*x6*x6

解析解は算出できなったため,2つの異なる積分方法で結果を確認。

Laporta (d)SR16000/M1,xm1 １ノード。
サイズＮ＝１２８、演算量 2668GFLOP

実行時間（秒）一覧表

smp数

xm1

xm1

sr16000/M1

SR16000/M1

f90

xlf_r

f90

xlf_r

32

30.00552.86

11.61538.46

64

31.151

54.52

7.498

38.51

128

19.36954.08

f90,SR16000/M1 の実行効率36.2%と突出して良い値となっています。

laporta (f)
サイズＮ＝１２８、演算量 3642GFLOP
ソースは二次元積分

スーパーコンピュータとＧＰＵ
実行時間(秒)
CPU

ノード数 smp数

BG/Q
SR16000
HD5870HD6970

32
1
2
4

mpi数

1
64
3200
6144

実行時間 性能
実行効率(秒)(GFLOPs) (%)
2048 1.289
2825
43.1
1 8.808413
42.2
1 8.094450
41.4
1 4.252857
31.7

どの機種も高い実行効率をしめしています。
ＧＰＵとＨＯＳＴ
実行時間(秒)cpu
HD5870HD6970

gpu

host
8.094
4.252

host/gpu
2459
304
3009
708

各種計算機詳細GPU
実行時間（秒)一覧表
ボード数

HD5870HD6970
gpu0
gpu1
1
15.67815.70215.428
2
8.08
7.894
3
5.448
4
4.251

BG/Q
case1 -O5 -qhot=level=2
case2 -O3 -qstrict -qhot=level=1
32ノード実行時間(秒)
case

flat
1
2

hybrid
1.289
1.324
4.855
4.106

GPU:ボード数の効果がよくでています。
BG/Q: コンパイルオプションの効果が大きい.
オプションによりflat ＭＰＩとハイブリッドMPIの効果が変わっています。

SR16000/M1
1ノード実行時間(秒)
smt
on
off

smp

実行時間 性能
実行効率(秒)(GFLOPs) (%)
64
7.909
460
47
32
11.674
312
31.8

T2K
実行時間（秒)一覧表
MPI数

実行時間 性能
実行効率（秒)(GFLOPs) (%)
64
17.275
211
35.6
128
8.646
421
35.8
256
4.33
841
35.8
512
2.199
1656
35.2

XM1
実行時間（秒)一覧表 論理コア数64
f90

xlf_r

smt
on
off
execution

32
88.45857.75629.01

6455.06318.95

実行時間（秒)一覧表 論理コア数 128
smp
f90
xlf_r

32
40.911
26.62

64
38.31
22.22

128
24.67720.82

SR16000/M1,T2K ともに実行効率は高くなっています。
Xm1: 論理コア数増加の効果はf90に大きくでています.

演算量 (N=120)
case GFLOP
laporta d2514
laporta f
2638
実行時間一覧表(秒)CPU SR16000
E5-2670
Phi5110P
HD5870プログラム 32core 64smp c
fortran cfortran 2面
laporta d 8.404 5.426 27.756 27.621 7.474 7.665 5.868
laporta f
9.841 6.373 28.275 27.921 7.324 7.462 5.928

演算量が約2500GFLOP程度だとSR16000とPhi5110P,HD5８70の性能はほぼ等しい.
これは4次元積分でも同様でした。

4.6.4.3 6次元積分計算
laportag
1 1 x1 1 x1  x2 1 x1  x2  x3 1 x1  x2  x3  x5 1 x1  x2  x3  x5  x6

I  
0 0


0


0


0


0

C
dx dx dx dx dx dx
3 7 6 5 3 2 1Dx4  1  x1  x2  x3  x5  x6  x7

C=x1*x4+x1*x5+x1*x6+x2*x4+x2*x5+x2*x6+x3*x4+x3*x5+x3*x6+x
4*x5& + x4*x6+x4*x7+x5*x7+x6*x7
D=&(x1**2+x2**2+x3**2+x7**2+x1*x2+x1*x3+x1*x7+x2*x3+x2*x7+x3*x7)&*(x4+x5+x6)&-x4**2*(x1+x2+x3+x5+X6+x7)&-(x5**2+x6**2+x5*x6)*(x1+x2+x3+x4+x7)&-3.0*x4*(x1*x5+x6*x7)&-2.0*((x1+x2+x3)*x4*x6+(x2+x3+x7)*x4*x5)

解析近似解=0.1723367907503

laportah
1 1 x1 1 x1  x3 1 x1  x3  x2 1 x1  x3  x2  x7 1 x1  x3  x2  x7  x6

I  
0 0

 
0

0


0


0

C
dx dx dx dx dx dx
3 4 6 7 2 3 1Dx5  1  x1  x3  x2  x7  x6  x4
C=(x1+x2+x3+x4)*(x4+x5+x6+x7)-x4*x4
cc=x1*m12+x2*m22+x3*m32+x4*m42+x5*m52+x6*m62+x7*m72
D=-c*cc.+s*(x1*x2*(x4+x5+x6+x7)+x5*x6*(x1+x2+x3+x4)+x1*x
4*x6+x2*x4*x5).+t*x3*x4*x7.+p12*(x1*x3*(x4+x5+x6+x7)+x3*x4*x5).+p22*(x2*x3*(x4+x5+x6+x7)+x3*x4*x6).+p32*(x5*x7*(x1+x2+x3+x4)+x1*x4*x7).+p42*(x6*x7*(x1+x2+x3+x4)+x2*x4*x7)

テストデータm12=1.0d0
m22=1.0d0
m32=1.0d0
m42=1.0d0
m52=1.0d0
m62=1.0d0
m72=1.0d0
p12=1.0d0
p22=1.0d0
p32=1.0d0
p42=1.0d0
s =1.0d0
t =1.0d0

解析近似解=0.1036407209893

aportai
1 1 x1 1 x1  x2 1 x1  x2  x3 1 x1  x2  x3  x4 1 x1  x2  x3  x4  x5

I  
0 0

 
0

0


0


0

C
dx dx dx dx dx dx
3 6 5 4 3 2 1Dx7  1  x1  x2  x3  x4  x5  x6
C=(x1+x2+x3+x4+x5)*(x1+x2+x3+x6+x7)-(x1+x2+x3)**2
cc=x1*m12+x2*m22+x3*m32+x4*m42+x5*m52+x6*m62+x7*m72
D=-c*cc.+s*(x1*x2*(x4+x5+x6+x7)+x1*x5*x6+x2*x4*x7x3*x4*x6).+t*x3*(-x4*x6+x5*x7).+p12*(x1*x3*(x4+x5+x6+x7)+x3*x4*(x6+x7)).+p22*(x2*x3*(x4+x5+x6+x7)+x3*x6*(x4+x5)).+p32*(x4*x5*(x1+x2+x3+x6+x7)+x4*x6*(x2+x3)+x1*x5*x7).+p42*(x6*x7*(x1+x2+x3+x4+x5)+x4*x6*(x1+x3)+x2*x5*x7)

テストデータm12=1.0d0
m22=1.0d0
m32=1.0d0
m42=1.0d0
m52=1.0d0
m62=1.0d0
m72=1.0d0
p12=1.0d0
p22=1.0d0
p32=1.0d0
p42=1.0d0
s =1.0d0
t =1.0d0
解析近似解=0.0853513981538

6次元積分実行結果一覧表
演算量(N=120)case
GFLOP
laporta g
312064laporat h243328laporta i
232896
実行時間一覧表(秒)cpu
SR16000
E5-2670
Phi5110P
HD5870プログラム 32core 64smp c
fortran cfortran 二面
laporta g 2223.744 1418.553 2385.352 2907.167 962.356 943.418 589.471
laporat h 1929.049 1256.823 2096.437 2717.488 925.742 959.321 568.483
laporta i 2212.909 1408.712 2467.079 4071.100 888.431 1063.900 606.803

演算量がこれだけ大きくなると,Phi5110P,
HD5870の性能がSR16000を大きく上回ります。
また,Phi5110P,HD5870のカタログ性能はほぼ同じですがHD5870の性能がPhi5110Pの倍近く良くなっています。これはPhi5110Pは60CPU,
240coreという構成によるものと考えられます。

laporta 計算結果一覧表
演算量(N=120)プログラム
LAPORTA DLAPORTA FLAPORTA GLAPORTA HLAPORTA I

GFLOP
2514
2638
312064243328232896

実行時間一覧表(秒)SR16000
プログラム
LAPORTA DLAPORTA FLAPORTA GLAPORTA HLAPORTA I

32core64smp

E5-2670

HOST

16smp

プログラム
LAPORTA DLAPORTA FLAPORTA GLAPORTA HLAPORTA I

cfortranPhi

240smp(60core)

プログラム
LAPORTA DLAPORTA FLAPORTA GLAPORTA HLAPORTA I

c

8.404
9.841
2223.744
1929.049
2212.909

5.426
6.373
1418.553
1256.823
1408.712

27.75628.2752385.352
2096.437
2467.079

27.62127.9212907.167
2717.488
4071.100

fortran
7.474
7.324
962.356925.742888.431

7.665
7.462
943.418
959.3211063.900

laporta

Phi5110P

実行時間(秒)一覧表

プログラム
laporta dlaporta flaporta glaporta hlaporta i

1smp /core 2smp /core 3smp /core 4smp /core
15.894091
8.930061
8.245382
7.477543
15.883854
9.015369
7.781932
7.302511
1567.222350 947.766614 956.547483 893.604456
1475.760124 890.712335 840.753795 794.503370
1682.613524 973.251626 936.788364 908.041320
グラフィックボード性能比較表
実行時間(秒)プログラム
laporta dlaporta flaporta glaporta hlaporta i

HD7970
W8000
HD79807.541418
7.764761
8.173080
7.600745
7.821768
8.249741
718.902868 738.792693 745.477537
684.978227 710.901497 745.486218
766.602135 788.137704 834.070693

SR16000/M1
LAPORTA (D)
MPI数
ノード数
128
256
256
512
LAPORTA(G)
MPI数
ノード数
128
256
256
512
LAPORTA (D)
(H)
MPI数
ノード数
128
256
256
512
LAPORTA (I)
MPI数
ノード数
128
256
256
512
S221

N=2048

MPI数

ノード数
128
256
256
512

GPGPU

複数ノード

実行時間
(秒）
4
3.231716
42.165075
8
1.617797
8
1.107155

実行時間
(秒）
4
684.465022
4
544.483556
8
345.817823
8277.552501

実行時間
(秒）
4
541.664086
4387.666776
8
270.937051
8
198.333612

実行時間
(秒）
4
728.709051
4
508.013126
8368.304897
8
258.708141

実行時間
(秒）
4
467.774509
4412.210191
8
233.916246
8
206.360301

LAPORTA
プログラム 実行時間（秒）D
5.85
G
589.553H
569.934I
611.884
S221
2166.834527

BG/Q
LAPORTA(D)
MPI数
ノード数
512
1024
2048

LAPORTA(G)
MPI数
ノード数
512
1024
2048

LAPORTA(H)
MPI数
ノード数
512
1024
2048

LAPORTA(I)
MPI数
ノード数
512
1024
2048

S221

N=2048

MPI数

ノード数

512
1024
2048

実行時間
(秒）
32 12.844877
32
6.63312732
3.967958

実行時間(秒）
32 977.217064
32 541.747376
32 439.703755

実行時間
(秒）
32 946.719043
32 514.122564
32 391.466375

実行時間
(秒）
32 954.046258
32 532.563084
32 430.696899

実行時間
(秒）
32 2144.31316
32 1161.5669
32 801.675789

理論最大性能(GFLOPs)
BG/Q 6553.6
SR16000 4node 3921.92
8node 7843.84
GPGPU 1088

4.6.5 3loop積分
今回扱った 3loop積分は
1 1 x1 1 x1  x 2 1 x1  x 2  x 3 1 x1  x 2  x 3  x 4 1 x1  x 2  x 3  x 4  x 5 1 x1  x 2  x 3  x 4  x 5  x 6

 
0

0

0





0



0



0

0

d  dx 7 dx 6 dx 5 dx 4 dx 3 dx 2 dx 1
でDが次の 3ケース ( N0, L01, L02)の場合です。
解析解は 20 5  20.7385551028673985...。

3 loop

実行結果一覧表N  25刻み幅N0h  4 . 734897

/ 24

L 01

h  4 . 696708

/ 24

L 02

h  4 . 705669/ 24

SR16000

64smp

プログラム 実行時間(秒)N0
93.6390L01
81.0711L02
67.2363

4倍精度

1
d
2
D

N0
x8=1.0-x1-x2-x3-x4-x5-x6-x7
x23=x2+x3
x24=x2+x4x247=x2+x4+x7
x1235=x1+x2+x3+x5
x2346=x2+x3+x4+x6
x5678=x5+x6+x7+x8

D = x24**2*x5**2-x2346*x247*x5**2-x2**2*x2346*x5678+2*x2*x23*x24*x5678-x1235*x24**2*x5678-x23**2*x247*x5678+x1235*x2346*x247*x5678+2*x2*x24*x5*x6-2*x23*x247*x5*x6+x2**2*x6**2-x1235*x247*x6**2+2*x2*x2346*x5*x7-2*x23*x24*x5*x7+2*x2*x23*x6*x7-2*x1235*x24*x6*x7+x23**2*x7**2-x1235*x2346*x7**2

L01
x8=1.0-x1-x2-x3-x4-x5-x6-x7
x123=x1+x2+x3
x167=x1+x6+x7
x348=x3+x8+x4
x2578=x2+x5+x7+x8

D = -x167*x2578*x3**2-x167*x2**2*x348-x1**2*x2578*x348+x123*x167*x2578*x348-2*x1*x2*x348*x7+x3**2*x7**2-x123*x348*x7**2-2*x167*x2*x3*x8-2*x1*x3*x7*x8+x1**2*x8**2-x123*x167*x8**2

L02
x8=1.0-x1-x2-x3-x4-x5-x6-x7
x123=x1+x2+x3
x1456=x1+x4+x5+x6
x678=x6+x7+x8
x247=x2+x4+x7

D = x2**2*x6**2-x123*x247*x6**2-x1456*x2**2*x678-x1**2*x247*x678+x123*x1456*x247*x678-2*x1*x2*x4*x678-x123*x4**2*x678-2*x1*x2*x6*x7-2*x123*x4*x6*x7+x1**2*x7**2-x123*x1456*x7**2

4.6.6 4loop積分

今回扱った 4loop積分は
1 1 x1 1 x1  x 2 1 x1  x 2  x 3 1 x1  x 2  x 3  x 4 1 x1  x 2  x 3  x 4  x 5 1 x1  x 2  x 3  x 4  x 5  x 6 1 x1  x 2  x 3  x 4  x 5  x 6  x 7

 
0

0

0


0


0


0


0

1
d  dx 8 dx 7 dx 6 dx 5 dx 4 dx 3dx 2 dx 1
CD
でC, Dが次の 2ケース (M 44 , M 45 )の場合です。
解析解は441 7
M 44は
 55.5852539156784...。
8
2
M 45は36 3  52.017868743610.....。
演算精度に関しては , 解析解と10進10桁以上合わせるにはM 44 では倍精度演算で事足 ります。
M 45 では , 拡張倍精度 ,4倍精度演算が必要。

0

M 44実行結果一覧表
条件 : 分点数25
刻み幅h  4.813324 / 24
t max  2.406662,  2  2  49.88652029 ≒ 9.6  10 16

E5-2670

16smp

ーＯ１ －ＯＰＥＮＭＰ

言語FORTRANCCFORTRANFORTRAN精度
実行時間(秒) 備考
倍精度
206.9346
倍精度
227.7403
拡張倍精度
491.31104倍精度
22487.9081
4倍精度
16751.9997 ーＯ２

Hpi5110P

240SMP

言語
FORTRAN
C
C

精度
実行時間(秒)
倍精度
1569.4017
倍精度
1455.5575
拡張倍精度
2334.1927

e5430

8smp

ーＯ１ －ＯＰＥＮＭＰ

言語
C

精度
倍精度

実行時間(秒)597.7637

ーＯ１ －ＯＰＥＮＭＰ

M 45 実行結果一覧表
条件 : 分点数 25
刻み幅 h  4 .78558154 / 24
t max  2 .39279077 ,  2  2  49.18797875 ≒ 1 .6  10 15
E5-2670

16smp

言語
C
FORTRANFORTRANFORTRAN精度
実行時間(秒) 備考
拡張倍精度
688.6699
4倍精度
39713.7063
4倍精度
30367.7276 ーＯ２
4倍精度
30745.1106 ーＯ３

SR16000

64smp

ーＯ３ －OMP

言語FORTRAN精度
4倍精度

実行時間(秒)
6387.4925

Phi5110P

240 smpーＯ１ －ＯＰＥＮＭＰ

言語
C

精度
実行時間(秒)
拡張倍精度
3379.0798

ーＯ１ －ＯＰＥＮＭＰ

性能的には拡張倍精度演算を使用するのが良い事を示しています。

M 44
x9 = 1.0 - x1 - x2 - x3 - x4 - x5 - x6 - x7 - x8
x123 = x1 + x2 + x3
x167 = x1 + x6 + x7
x278 = x2 + x7 + x8
x349 = x3 + x4 + x9
x589 = x5 + x8 + x9
C = x167 * x278 * x349 * x589 - x349 * x589 * x7 * *2 - x167 * x349 * x8 * *2
- x167 * x278 * x9 * *2 + x7 * *2 * x9 * *2
D = - (x167 * x278 * x3 * *2 * x589) - x167 * x2 * *2 * x349 * x589
- x1 * *2 * x278 * x349 * x589
+ x123 * x167 * x278 * x349 * x589 - 2 * x1 * x2 * x349 * x589 * x7
+ x3 * *2 * x589 * x7 * *2
- x123 * x349 * x589 * x7 * *2 + x167 * x3 * *2 * x8 * *2
+ x1 * *2 * x349 * x8 * *2
- x123 * x167 * x349 * x8 * *2 - 2 * x167 * x2 * x3 * x8 * x9
- 2 * x1 * x3 * x7 * x8 * x9
+ x167 * x2 * *2 * x9 * *2 + x1 * *2 * x278 * x9 * *2
- x123 * x167 * x278 * x9 * *2
+ 2 * x1 * x2 * x7 * x9 * *2 + x123 * x7 * *2 * x9 * *2

M 45
x9 = 1.0 - x1 - x2 - x3 - x4 - x5 - x6 - x7 - x8
x123 = x1 + x2 + x3
x1567 = x1 + x5 + x6 + x7
x279 = x2 + x7 + x9x349 = x3 + x4 + x9
x589 = x5 + x8 + x9
c = -x279 * x349 * x5 * *2 + x1567 * x279 * x349 * x589 - x349 * x589 * x7 * *2
- 2 * x349 * x5 * x7 * x9 - x1567 * x279 * x9 * *2 - x1567 * x349 * x9 * *2
+ x5 * *2 * x9 * *2 - x1567 * x589 * x9 * *2 + 2 * x5 * x7 * x9 * *2
+ x7 * *2 * x9 * *2 + 2 * x1567 * x9 * *3
d = x279 * x3 * *2 * x5 * *2 + x2 * *2 * x349 * x5 * *2 - x123 * x279 * x349 * x5 * *2
- x1567 * x279 * x3 * *2 * x589 - x1567 * x2 * *2 * x349 * x589
- x1 * *2 * x279 * x349 * x589 + x123 * x1567 * x279 * x349 * x589
- 2 * x1 * x2 * x349 * x589 * x7 + x3 * *2 * x589 * x7 * *2
- x123 * x349 * x589 * x7 * *2 + 2 * x1 * x279 * x3 * x5 * x9
- 2 * x1 * x2 * x349 * x5 * x9 + 2 * x2 * x3 * x5 * *2 * x9 - 2 * x1567 *x2 * x3 * x589 * x9
+ 2 * x2 * x3 * x5 * x7 * x9 + 2 * x3 * *2 * x5 * x7 * x9 - 2 * x123 * x349 *x5 * x7 * x9
- 2 * x1 * x3 * x589 * x7 * x9 + x1567 * x2 * *2 * x9 * *2 + x1 * *2 * x279 * x9 * *2
- x123 * x1567 * x279 * x9 * *2 + 2 * x1567 * x2 * x3 * x9 * *2
+ x1567 * x3 * *2 * x9 * *2 + x1 * *2 * x349 * x9 * *2 - x123 * x1567 * x349 * x9 * *2
+ 2 * x1 * x2 * x5 * x9 * *2 - 2 * x1 * x3 * x5 * x9 * *2 + x123 * x5 * *2 * x9 * *2
+ x1 * *2 * x589 * x9 * *2 - x123 * x1567 * x589 * x9 * *2 + 2 * x1 * x2 * x7 * x9 * *2
+ 2 * x1 * x3 * x7 * x9 * *2 + 2 * x123 * x5 * x7 * x9 * *2 + x123 * x7 * *2 * x9 * *2
- 2 * x1 * *2 * x9 * *3 + 2 * x123 * x1567 * x9 * *3

多倍長計算手法まとめ

平成２3年度-平成27年度
1. はじめに
2.多倍長精度演算
2.1 浮動小数点演算方式
2.2 整数演算方式
2.3 ルーチン作成前に決定すべき事項
2.4 ルーチン作成の要点
2.5 注意点詳細
2.6 整数演算方式の並列化に関して
3.実アプリケーション
3.1 ヒルベルト行列Hによる連立一次方程式Hx=bの求解
3.1.1 ヒルベルト行列の条件数
3.1.2 ヒルベルト行列実測結果
3.1.2.1 64倍精度演算までの結果
3.1.2.2 128-248倍精度演算での結果
3.1.2.3 308-488倍精度演算での結果
3.1.2.4 548倍精度演算での結果
3.2 対称疎行列Aによる連立一次方程式Ａx=bの求解
3.3 非対称疎行列Aによる連立一次方程式Ａx=bの求解
3.3.1 nx=ny=nz=65の場合
3.3.1.1 bcg法の収束状況
3.3.1.2 ｃｇｓ法の収束状況
3.3.2 サイズを拡大した場合
3.3.3 精度改善例

多倍長計算手法まとめ
3.4 ファインマンループ積分
3.4.1 数値積分法
3.4.2 massless計算
3.4.3 bsgamma
3.4.4 4次元積分
3.4.5 two loop vertex
3.4.6 3loop積分
3.4.7 4loop積分
3.4.8 3loop,4loop積分作業手順まとめ
3.5 量子モンテカルロ法による物性
スペクトル計算
3.5.1 指数部11ビットの限界について
3.5.2 必要演算精度見積り

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

c

implicit real*4(a-h,o-z)
s=0.0
a=1.0
do i=1,100000000
s=s+a
end dowrite(6,1) s
1 format(1h ,2x,f20.5)
stop
end
正しい結果は100000000ですが次ページ以降に示すように最適化オプションや並列化
で様子が異なります。

SR16000/xm1 f90 並列化をしない場合% f90 a.f -O3 -o a
f90: compile start : a.f*OFORT90 V03-02-/B entered.
*program name = MAIN
*end of compilation : MAIN
% ./a &%
16777216.00000% f90 a.f -Os -noloopexpand noparallel -o a
f90: compile start : a.f*OFORT90 V03-02-/B entered.
*program name = MAIN
*end of compilation : MAIN
% ./a16777216.00000

% f90 a.f -Os -noparallel -o a
f90: compile start : a.f*OFORT90 V03-02-/B entered.
*program name = MAIN
*end of compilation : MAIN
% ./a100000000.00000
SR16000/xm1 f90 並列化をした場合% f90 a.f -Os -o a
f90: compile start : a.f*OFORT90 V03-02-/B entered.
*program name = MAIN
*end of compilation : MAIN
*end of compilation :_parallel_func_1_MAIN% ./a100000000.00000%

次にスレッド数指定で異なる例を以下の例で示します。先のプログラムの100000000を1000000000に変更しただけのものです。

cimplicit real*4(a-h,o-z)s=0.0d0
a=1.0d0
C$OMP PARALLEL DO
C$OMP& reduction(+:s)
do i=1,loop
s=s+a
end doC$OMP END PARALLEL DO
write(6,1) s
1 format(1h ,2x,f20.5)
stop
end

E5  2670
16smp
268435456.00000  2 28  2 24  16  109
より結果は 2 28  2 24  16となる。
100smp
1000000000.00000  109  2 24  100
より正しい結果 となる。
Phi5110P
100smp 1000000000.00000
240smp 1000000000.00000
ここで 2 24 は単精度加算演算で正 確に計算できる個数となっています。

2. 多倍長精度演算演算方式には,浮動小数点演算方式と整数演算方式がある.

2.1 浮動小数点演算方式
多倍長精度変数を複数個の倍精度変数の和で表す。

表現できる数値範囲は,ieee754の倍精度と同じ.
a : 2n倍精度変数

a1 , a 2 ,...., a n : 倍精度変数a  a1  a 2  .....  a n
有効ビット数＝仮数部のビット数  1
2n倍精度(52  1)  n

ビット

n2

106ビットn3

159ビットn4212ビット

乗加算命令がある場合とない場合の処理に分かれる。
コンパイラの最適化方式に性能が大きく依存する。

次に示す基本演算例よりわかります。

139

基本の演算倍精度変数の 加減算 , 乗除算の結果を２つの 倍精度変数の和で表す 。
a , b, c, d倍精度変数 ; a  b, a  b, a  b  c  d
加算 : c  a  bt  cad  (a  (c  t ))  (b  t )減算 : c  a  bt  cad  (a  (c  t ))  (b  t )加算 , 減算ともに最適化オプ ションで括弧をはずし た演算順序にならない 様に注意が必要
乗算 (乗加算命令あり .)c  a*bd  a*b  c（最適化オプシ ョンにより d  0とならないように注意が必 要 )

乗算 (乗加算命令なし )r  134217729.0d 0(  2 27  1)c  a*bt1  a * r
a1  t1  ( t1  a )a 2  a  a1
t2  b * r
b1  t 2  ( t 2  b)b 2  b  b1d  ((a1 * b1  c)  a1 * b 2  a 2 * b1)  a 2 * b 2

演算量

加減算

n２に比例

乗算

n 3に比例n  3ではn個の変数の絶対値が大きい順にならべるソート処理が必要.
2．２ 整数演算方式
多倍長精度変数はieee754-2008の４倍精度の仮数部を変えたものであらわす。
表現できる数値範囲はieee754-2008の４倍精度と同じ.

a : p倍精度( p  3)変数

ｐ倍精度

符号部 1ビット

演算量

指数部 15ビット
仮数部 32  p  16ビット

加減算 pに比例
乗算

p 2に比例

有効ビット数p  4 113ビット
p  6 177ビット
p  8 241ビット

符号部
１ビット

指数部
15ビット

仮数部
32*p-16ビット(p=4: 112,p=6: 176,p=8: 240)

141

2.3

ルーチン作成前に決定すべき事項

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

2.4 ルーチン作成の要点
最初に加減算,乗算ルーチンを作成し平方根は既存の倍精度演算を使用して多倍長演算の初期値として反復法により求める。除算に関して
は,浮動小数点演算方式では既存の倍精度除算
を複数回使用する。6倍精度以上(使用する倍精度変数が3以上の場合,絶対値が大きい順になる様にする注意が必要。整数演算方式では逆数の初期値を既存の4倍精度除算を使用して反復法により逆数を求め,被除数との積をとり求める。

数学関数等の作成は上記の5演算の組み合わせで行う。具体的に立方根,対数関数,指数関数,正弦関数計算でその注意事項を記述する。

引数の取り得る範囲 (ここでは表現可能な数 値の絶対値
最大値を Mとする )は立方根 , 指数関数 , 正弦関数は | x | M.
対数関数は 0  x  M。ここで関数値 yのとる範囲は ,立方根 , 対数関数は | y | M, 指数関数は 0  y  M。
正弦関数は | y | 1。ここで , 指数関数は x  0でy  e x
1
より | x | log(M )となる.
ex
立方根に関しては x  y  z ( yの立方根は容易に計算 出来るより0  x  log(M ), e  x 

数値 , zをある標準区間 )とすると , 3 x  3 y  3 z より, x  y  zの分解によって誤差は発 生しにくい。立方根に 関してはこの関数が準備されていないシ ステムもあり , x  e
3

1
x log( x )3(べき乗 )で計算される事があるが x  0ではエラーとなるので 注意が必要である。
判れば簡単な事ではあ るが , 多くの人が同じ間違い を繰り返した問題であるので注意が 必要である。
指数関数及び正弦関数 はテーラー展開式がある が, 引数範囲すべてで計算 すると
性能, 精度の問題が発生する 。また対数関数では引 数範囲のテーラー展開はなく , x  y  zの様に分解してある近 似区間で zで計算して xを求めるが , e x  e y  e z なので指数関数は誤差 が発生しにくい。対数関 数と正弦関数を比較す ると, 結果の範囲
が | y | M, と | y | 1から , 正弦関数の計算の方が 誤差が出易い
事がわかる。正弦関数 の計算方法の一例をあ げる。





], t | x |  I  とし 0  t  ,0   t  より
2
2
2
2
2

sin( x )   sin( t ), sin(  t )で計算する。
2
ここで問題なのは , | x | が大きい場合 , xとtの持つ精度が異なるI  [| x | /ことである。（ tを求める際の桁落ちの ため ).また, | x | が小さくてもxがnに近い場合 , tは桁落ちが発生するが この場合は sin( x ) ≒ 0のため , 相対誤算は大きいが絶 対誤差は小さくて済む 。精度上の問題が発生した場合に は, 作成者により精度に関 する思想、扱いが異なるという意識を 持つことが必要となる 。

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

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

これを以下のRump’s例題を使用して検証しました。
Rump' s例題a  77617.0

b  33096.0

f  333.75b 6  a 2 (11a 2 b 2  b 6  121b 4  2)a
 5.5b 
2b54767
理論解   0.827396059946
66192[333.75b 6  a 2 (11a 2 b 2  b 6  121b 4  2)]  5.5b 88の演算で121ビットの桁落ちが発生 。
桁落ちが影響を受ける と0, 受けないと  2となる。
4倍精度演算では有効ビ ット数が113のため桁落ちの影響を受ける ため1.172603940053となり,有効ビット数が121より大きい演算(5倍超精度演算）では桁 落ちの影響を受けないので  0.827396059946となる。

Core内simdは適用できないのでcore間ｓｍｐで実行している。

多倍長rump' s例題一覧表
1,000,000回実行
実行時間一覧表(秒)E5-2670
演算精度
68倍精度
128倍精度188倍精度248倍精度308倍精度368倍精度428倍精度488倍精度548倍精度

Phi5110P
4smp/core
16smp
240smp11.9698
46.3840
42.3808 156.0239
89.0890 331.2860
152.0543 585.7436
248.2340 966.5418
355.4690 1378.1645
476.9385 1859.8102
617.0996 2410.4360
833.8722 3244.1159

Core間のスレッド並列の効果が出ている。

多倍長rump' s例題一覧表
1,000,000回実行

実行時間一覧表(秒)E5-2670
演算精度
68倍精度
128倍精度188倍精度248倍精度

16smp
8smp
4smp
2smp
1smp
11.9698
23.8582
47.858595.5189 191.5384
42.3808
84.7253 169.0772 338.4619 677.3167
89.0890 177.8675 355.6048 711.6836 1425.8922
152.0543 303.7099 607.2627 1214.3922 2438.3971

Phi5110P

4smp/core

演算精度
68倍精度
128倍精度188倍精度248倍精度

240smp180smp120smp60smp
46.3840
48.7955
50.8138
69.1542156.0239 178.3118 181.1152 244.4755
331.2860 378.7431 384.5338 525.7136
585.7436 669.6312 694.6440 938.0996

演算精度
68倍精度

4smp
1smp
868.7135 3463.0626

3.実アプリケーション
3.1 ヒルベルト行列Hによる連立一次方程式Hx=bの求解

3.1.1 ヒルベルト行列の条 件数
1(1  i, j  n )i  j 1(n  i  1)!(n  j  1)!
1
逆行列 H ji  (1) i  j(i  j  1)[(i  1)!( j  1)!] 2 (n  i)!(n  j)!
1
1
H ij のノルムは1   ...... nが大とすると ,2nln(n )     0.57721566490153286060 (オイラー数 )ヒルベルト行列 H ij 1

逆行列 H ji のノルムを | H ji

1

| の最大値とする。

また nは大として、スターリ ングの公式 n!~n

2 e nn

12を使用する。

逆行列の絶対値 Pは変形して ,(n  i)!(n  j)!i 2 j2P(n  i)(n  j)(i  j  1)[i! j!] 2 (n  i)!(n  j)!
i  an, j  bnとすると ,a 2 b 2 n 2 ((1  a )n )!((1  b)n )!
P
((a  b)n  1)(1  a )(1  b)[(an )!(bn )!] 2 ((1  a )n )!((1  b)n )!
1 a 1 b
1  a (1a ) n 1  b (1 b ) n 1  a 2 an 1  b 2 bn
1 a 1 b

(
)
(
)
(
) (
)
1 b
a
b
(1  a )(1  b)((a  b)n  1)4 2 1  a
ab

Pはaとbに対して対称より , a  bとすると、
a2
1  a 2 (1a ) n 1  a 4 an
P
(
)
(
)
2
2
a
(1  a )(2an  1)4 1  a
1  a (1a ) 1  a 2 a
1 a2
f (a )  ln[()(
) ]とすると f ' (a )  ln( 2 )
1 a
a
a

0.5  a  0.75では f (a )はa  0.5  0.7071で最大となる .
i  an, j  bnが整数となる事から , a  0.7とする.
0.0243 17 2.8 n 17 0.6 n 0.0243 1.530979 n
10( ) ( )
1 .4 n  1
3
1 .4 n  1 7より詳細にすると以下 の様になる。

すると Pmax a  0. 7

0.02433695097 17 0.6 n 17 2.8 n( ) ( )
73
1. 4 n  1
 101.530979068n 1.613733833log10 (1.4 n 1) (1)Pa  0 .5 

2
2

P

1
( 2n  1)4(2

2  2 ( 2 2 ) n( 2  1) 2 2n)2 2

 101.531102741n 1.596359739log10 ( 2n 1) (2)(2) /(1)  10 0.000123673n  0.017374096 log10 (1.4 n 1) log10 ( 2n 1)
10 m のm(1)( 2)(2) /(1)

条件数比率n  200 302.1364756 302.1741817 0.0377

1.091

n  400 608.0304816 608.0929301 0.0624

1.155

近似的には条件数  101.5 N となる.

条件数は Ax  bを(A  A )( x  x )  (b  b)で計算する場合の xの大きさに関する目安 となるものである。ヒルベル ト行列 H (ここでは Aとする )の係数
を8倍精度演算で求めると します.ここで解 x  (1,1,.....,1) Tとなる様に bを8倍精度演算で計算しま す.すなわち
A, b  0の場合 ,8倍精度演算で計算する と,40
DIMENSION SIZE 
int 8bai seido
gosa  5.453971473889753996929348189724412E - 0016
60
DIMENSION SIZE 
int 8bai seido
gosa  86.1703510314325363937588497909226となる.
ここで Aの係数を 4倍精度演算でもとめる と, A, b  0となり
200
DIMENSION SIZE 
int 8bai seido
gosa  3.305918377126923593570282216668663E - 0036
Aの係数を倍精度演算で 求め, 他の演算を4倍精度で行った場合で も200
DIMENSION SIZE GOSAMAX  7.568098984088555143981688648697012E - 0017となる。
この事は , 入力データは低精度 (もともとデータで 10進 30桁の精度がある事は稀である .)で演算のみを高精度で 行い , 計算結果を低精度でも出力すると いう場合に適用できる 。
多倍長精度演算のテス トでは , A, b  0の場合で実施。

3.1.2 ヒルベルト行列実測結果
3.1.2.1 64倍精度演算までの結果
ヒルベルト行列最大誤差一覧表
次元数

倍精度
４倍精度
８倍精度
20
51.616
1.44E-07 2.13E-47
40
182.10762.034 8.72E-16
60 1971.569
62.804251.13280642.074
300.734
252.653100
201.355
378.428
242.672次元数

16倍精度 ３２倍精度
20 1.47E-99 1.44E-276
40 4.09E-86 4.30E-243
60 2.66E-60 1.36E-210
80 3.65E-31 1.08E-178
100 2.58E+00 1.41E-150

ヒルベルト行列の精度
N

演算精度 最大誤差 条件数(ビット数表示)
170 32倍精度 4.11E-48
853
180 32倍精度 2.00E-32
904
190 32倍精度 7.11E-17
955
200 32倍精度 1.27E-02
1006
210 32倍精度 1.37E+02
1057
370 64倍精度 9.63E-35
1870
380 64倍精度 7.98E-30
1920
390 64倍精度 7.61E-15
1971
400 64倍精度 7.10E+01
2022

3.1.2.2 128-248倍精度演算結果
ヒルベルト行列の精度

N

最大誤差
400
500
600
700
800
1000
1100
1400

精度

1.014E-621 128倍精度
1.681E-467 128倍精度
1.321E-314 128倍精度
1.349E-161 128倍精度
7.603E-009 128倍精度
1.013E-280 188倍精度
4.211E-127 188倍精度
1.013E-245 248倍精度

最大誤差のビット数  有効ビット数  条件数のビット数となっています。
248倍精度演算ではN  1400が使用できるメモリ容量の上限となっています。

3.1.2.3 308-488倍精度演算での結果

N=100でのヒルベルト行列計算
条件数0.127E+152
精度

ビット数501

有効ビット数 誤差
実測
理論値
308
9841 4.73E-2814 2.80E-2812
368
11761 3.38E-3392 2.52E-3390
428
13681 1.61E-3971 2.66E-3968
488
15601 5.36E-4198 2.80E-4546 (注)488
15601 1.00E-4547 2.80E-4546(注)逆数計算で4倍精度の初期値から7回反復した場合の値
308,368,428倍精度は7回反復,488倍精度は8回反復が必要

3.1.2.4 548倍精度演算での結果
誤差は解x  (1,1,....,1) T と比較していますので、小さなNに対しては差が 0と正しく計算しているかの検証とは言えませ ん。そこでNを変えながら正しい事を検証しました。

x5570
最大誤差

実行時間
条件数
有効ビット-条件数(秒)(ビット数)
0.000D-0000 182.566
501 0.000D-0000
0.000D-0000 1304.570
1010 0.000D-0000
1.647D-4819 4368.843
1518 4.140D-4818
1.656D-466710544.205
2022 2.168D-4666
4.445D-451321630.467
2530 1.817D-4513

誤差が出るサイズ
条件数1138ビットとなるＮN

条件数(ビット数)
223
1127
224
1132
225
1137
226
1142
227
1147
228
1152

3.2 対称疎行列Aによる連立一次方程式
Ax=bの求解
扱った問題
ポアソン方程式u  f (), u  0 ()  [0,1]  [0,1]  [0,1]N  200  200  200
収束判定値 : 共役残差10 12
u  sin( x ) sin( y) sin( z)f  3   2  sin( x ) sin( y) sin( z)

初期値は x (i)  1.0
cgs, cgsilu, cgsmilの非対称疎行列用反復法では収束しな い。
この様な場合 , 初期値は
x (i) : (0,1)の一様乱数値としてい る。

倍精度演算実行での反復回数
反復回数
解法
sor
sor2
sor3
adi
cg
bcg
bicgs
cgs
cgsm

反復回数 備考
4293 陰解法
141799 陽解法
4275 陰解法+odd-even法
34195 陽解法
690
704
495
640 初期値x(i)=1.0では収束しない。
288 cgsiluの並列版

前処理付き
scg
cgsmilcgsilucgsmm
gpbicg693 前処理付きの最も簡単なケース
241 初期値x(i)=1.0では収束しない。
272 初期値x(i)=1.0では収束しない。
306 ｃｇｓｍｉｌの並列版
142

並列実行にすると , 初期値 x (i)  1.0でも収束するように改善さ れる場合がある。
反復回数の削減には ,初期値の変更 , 並列化などがある事になる。

3.3 非対称疎行列Aによる連立一次方程式
Ax=bの求解問題は 3次元ポアソン方程式
uu  R fx
領域[0,1]  [0,1]  [0,1]解析解 u ( x , y, z)  e xyz  sin( x )  sin( y)  sin( z)
非対称問題。
解法
bcg(biconjugate

gradient )法R  100, nx  ny  nz  65
R  100, nx  ny  nz  129
cgs(conjugate

gradient

square)法R  100, nx  ny  nz  65
R  100, nx  ny  nz  113
R  100, nx  ny  nz  129
収束判定値は共役残差 10 12
初期値 u ( x , y, z)  0
サイズは精度改善によ り拡大しています。

3.3.1 nx=ny=nz=65 の場合

3.3.1.1 BCG法の収束状況

bcg 法

反復回数
3500
3000

2911

倍精度演算では反復回数10000
でも収束せず

2500
2000

反復回数

1500

1000
500

1075

346

173

163

0
60 68 76 81 86 94 102113177241

ﾋﾞｯﾄ数

14
12
10
8
6
4
2
0
-2
-4
-6
-8
1
9
17
505
513
1001
1009
1017
2005
2013
2790
2798
2806

-log10(残差)

倍精度

反復回数

BCG法の収束状況
nx=ny=nz=65 収束判定値=1.0q-12-log10(残差)

15
10
5

56ﾋﾞｯﾄ
80ﾋﾞｯﾄ

-5

1
8
15
22
29
36
43
50
320
345
1091
1098

0

-10
反復回数

-log10(残差)

15
10
5

1
17
33
49
65
81
97
113
129
145
161
177

0
-5
-10
反復回数

4倍精度
5倍精度
6倍精度
7倍精度
8倍精度

CGS法の収束状況
nx=ny=nz=65 収束判定値=1.0q-12
倍精度演算では収束しない。

cgs 法

反復回数
3500
3000

2911

2500
2000

反復回数

1500
1000
587

500

364
151

128

0
97

99

102 106 145 177 212

ﾋﾞｯﾄ数

3.3.1.2

CGS法の収束状況-log10(残差)
４倍精度

15
10
5
4倍精度(n=2)
301
320
339
358
377
396
515
534
553
572
591

0
-5
-10
反復回数

-log10(残差)６倍精度,8倍精度

15
10
5
-5

51
60
69
78
87
96
105
114
123
132
141

0
-10
-15
-20
-25
-30

反復回数

6倍精度(n=3)8倍精度(n=4)

ｃｇｓ法反復回数
倍精度は10000回で収束せず。

精度
4倍精度
5倍精度
6倍精度
7倍精度
8倍精度

反復回数
ieee
ieee smp dd
392
448
670
226
155 none
130
131
144
128
128 none
128
128
128

q4sum,q5sum 使用
精度
4倍精度
5倍精度

反復回数
378
155

(1)ieee形式,dd形式ともに6倍精度演算が最も効率が良い事を示しています。
(2)q4sum,q5sumの様に総和演算を無限精度演算を使用すると,5倍精度の収束が良くなり,ieeeでのsmp実行と同じ反復回数となっています。

詳細な収束状況

dd形式cgs反復回数
ビット数

反復回数
95
4931
96
2155
97
3506
98
2811
99
1250
100
1103
101
1080
102
1027
103
815
104
702
105
630
106
670

3.3.2 サイズを拡大した場合

通常サイズを大きくすると収束に要する反復回数は増えますが、その状況を以下の条件で検証しました。
解法
bcg ( biconjugate

gradient ) 法

R  100, nx  ny  nz  129
cgs(conjugate

gradient

square) 法R  100, nx  ny  nz  113
収束判定値は共役残差 10 12
初期値 u ( x , y, z )  0の場合

反復回数10000まででは、
bcg法は倍精度演算,
ｃｇｓ法は4倍精度演算まで収束しません。

bcg 法反復回数一覧
演算精度 SR16000 x5570
E5-2670 Phi5110P
3倍精度
482
4倍精度
316
308
308
314
5倍精度
317
6倍精度
317
7倍精度
317
8倍精度
317

cgs 法反復回数一覧表
演算精度 SR16000 x5570
E5-2670 Phi5110P
5倍精度
373
368
368
6倍精度
281
237
240
240
7倍精度
217
217
217
8倍精度
217
217
217
217

(注 )
1.赤字はその計算機で最 も性能がよかったもの 。
2.6倍精度の有効ビット数 はSR16000は
159ビット, 他は177ビット
3.E5  2670, Phi 5110Pでは内積部分に無限精度演算を適 用しています。

bcg法反復回数
精度
56ビット
60ビット
64ビット
68ビット
72ビット
76ビット
3倍精度
84ビット
88ビット
92ビット
96ビット
100ビット
104ビット108ビット
4倍精度

反復回数 正規化
実測/正規化
50001 none
none
4315 node
none
1579
931
1.696
938
826
1.136
784
738
1.062
593
664
0.893
482
600
0.803
439
545
0.806
396
497
0.797
368
455
0.809
336
418
0.804
325
386
0.842
319
357
0.894
313
331
0.946
308
308
1

(注)50001回は50000回反復で収束しなかった事を示しています。
正規化
P有効ビット数=精度のビット数+1

正規化反復回数308  (

113
)2P

性能的には3倍精度が最も効率が良いと言えます。

3.3.3 精度改善例

以下の様に精度改善を行う事により,これまでのサイズおよびサイズをさらに拡大しても4倍精度演算で収束するようになりました。
前処理付き cgs法により, 並列処理で精度改善
cgs(conjugate

gradient

square)法R  100, nx  ny  nz  129
収束判定値は共役残差 10 12
初期値 u ( x , y, z)  0
4倍精度演算CGS (PARALLEL)SIZE NX, NY, NZ = 129
EPS = 1.00000000E - 0012
ITERATION =

142

129

129

3.4 ファインマンループ積分
ファインマンループ積分の一般式は以下の様になる。
N  n ( L 1)

2
1
N
1 nL
C
I()  (1) N ( ) 2 ( N  nL )   dx i (1  x 1  x 2  .....  x N )N  nL2 0 i1
4
(D  iC) 2

n : 時空次元 (n  4), N : ループ内の素粒子の数 , L : ループの多重度
 : 実定数 , C, D : x 1 , x 2 ,..., x N の多項式演算精度が問題となるのは,以下の例の様に仮想光子λが含まれる場合である。
(1) Infra vtx

1
2
I(s)    dydx (D   xys  ( x  y) 2 m e  (1  x  y)2 )0 0 D
1 1 x(2) Infra box

1
dzdydx2
0 0
0D(D   xys  tz (1  x  y  z)  ( x  y)2
1 1 x1 x  yI(s)    

 (1  x  y  z)(1  x  y)m e  z(1  x  y)m f )22

３．４．１数値積分法
数値積分法としては、(1)誤差を評価しながら領域を二分割して行く
Gauss-Kronrod積分法(2)端点特異点を持つ被積分関数の積分に強く,基本的には領域は分割せず, 標本点をふやして積分計算結果の精度をあげていく
二重指数関数型積分法がある。

これ以外にも一重指数関数型積分法もあるが,二重指数関数型積分法に比べ標本点数を大きくとる必要がある。
端点と領域内部で同時 に特異点をもつ場合の計算を , sin c求積法による Hadamard
有限部分計算で行う方 法もあるがファインマ ン
ループ積分への適用出 来る例は少ない。
数値積分で積分領域内 に特異点がある場合
微小量 を使用して D   D  iと有理化して   0で得た値を求める積分 値とする場合に加速法を使用するが主 にRicardson加速法 ,エイトケン加速法、 Wynnの  算法等を用いる。

Gauss  Kronrod積分法と二重指数関数 型積分法の選択に関しては以下 の差があります。
分点と重み係数のテー ブル作成方法ではGauss  Kronrod法では , まず Gauss  積分法で2n  1の分点と重み係数を求 めます。
 1  x 1  x 2  ....  x n  x n 1  0  x n  2  x n 3  ...  x 2 n 1  1
これに 2n  2の分点を追加して計 4n  3個の分点と重み係数を求めます。これは反 復法により求めるため , 積分計算の前に分点と重み係数の 計算が必要となります 。
これに対し , 二重指数関数型積分は 分点と重み係数は関数計算 (sinh, cosh, exp)で求まるため , 積分計算のプログラムに入れる事が可能です 。二重指数関数型積分 の誤差評価
は, N  2n  1で求めた値 w 1 , w 2 ,......., w 2 n 1で
2 n 1n 1

I1  h   w とI 2  2h   w 2i 1を比較して求められま すし ,i 1i

誤差  exp(c 2i 1N)で評価ができます。
log( N)

内部特異点がある場合 は領域分割によって計 算が可能です。
幾つかの例を次ページ に記しました。

例題
4倍精度演算 , 分点数 496, 変数変換区間 [ 1,1]で計算.分母 dでd  0となる点で領域を分割 し ,それぞれの領域で d   d  i( 0  1.2 30 ,...., 14  1.2  44 )
  算法で計算
1xdx  3.141592653589793
1
1 ( x )e
計算結果  3.141592653589793
2x0 (sin x  2 cos x ) 2 dx  0.7669479668299477
計算結果  0.7669479668299477
2

x sin 2 x
0 1  2  2.0 cos 2x  2.0 2 dx  0.1592257756149500
計算結果  0.1592257756149501
11
0 log | log(x ) |dx  0.154479641320
計算結果  0.154479641320(注 )この積分値は約 200年間 0と考えられていました 。

1I   f ( x )dx

1

f (x) 1

3
4

( 2  x )(1  x ) (1  x )14

h  0.2, N  40 10 13 の相対誤差
C 2  2.76
f (x) exp(C 2N)ln( N )

1
4

3
4

1
4

(1  x ) (1  x )1 (1  x ) (1  x )
( 2  x )(1  x )(1  x )
3
2x

3
4

3

1

3
1
1
1
4
4
4
4
 (1  x ) (1  x )  (1  x ) (1  x )
6
2

y

1
( x  1)2dx  2dy

1

x  1  2 y,1  x  2(1  y)

1
4

1

3
4

1

3

3

1

1

1



4 y (1  y)1
I   f ( x )dx   dy   y 4 (1  y) 4 dy   y 4 (1  y) 4 dy
3 0 3  2y
30
1
0
1
4

1



3
4

4 y (1  y)
1 1
3
dy


(
)

(
)
3 0 3  2 y
2 4
4
1

1
4

3
4

1

1

1

1

3

1




y (1  y)11J
dy   y 4 (1  y) 4 dy   y 4 (1  y) 4 dy
3  2y
20
40
0
1



3
4

3 y (1  y)
 40
3  2y

14

1 1 431
3dy  (  ) ( ) ( )
8 4
4
4
4

4
1 4 3 431
3
12I [  ( )]( )( )   1.949054259166747
2 3 8
4
4
4
3
実測値 ( 4倍精度演算 )  1.949054259166746

被積分関数が内部特異点を持つ場合,
D->D-iεで有理化して, ε->0として積分値を求めるのが一般的であるが収束が遅くなったり,異なる積分値を得る場合がある。
(1) inf ravertex

1 1 x1
0 0 Ddydx( D  sxy  ( x  y) 2 m 2  (1  x  y)2 )ではs  0の場合の解析近似解は
2

2

me
1
2
1
2
2 me
(
)[ln(
) ln(
)  ln (
)
]
2
s
s
2
s
6
me
でs  0のときは , s   s  iで解析接続して求める 。

s  500 2 , m  0.0005でGauss  Kronrod法
4倍精度演算で計算した 時の実行時間 (秒 )は
以下の様になっていま す。
D   D  ies   s  ie   xys   s  iで計算非常に良い精度の結果が得られている。

(>>2) inf rabox

1 1 x 1 x  y

I 
>>

0

0


0

1
dzdydxD2

2
D>>   xys  tz (1  x  y  x )  ( x  y)  (1  x  y  z )(1  x  y) m e

 z (1  x  y) m f

2

2>> s  0の場合の解析近似解
で2( t  m f ) 2
sIln( 2 ) ln
222 s(  t  m f )me mf
1

から s  0の場合 s   s  iと解析接続して計算し ます。
1 1 x 1 x  y

I 
0

0


0

1
dzdydx
D2で D   D  iで有理化して得た虚数 部の値は の初期値  0 の取り方により結果が 大きく異なる。
(数列としては十分な精 度で収束している。)(1)  0 >> m e
2

20.142647572794761143d - 09(2) m e >>  0 >> 2

0.371536892867592250d - 08(3) 2 >>  00.743073784900021472d - 08

収束の速さでは圧倒的に(1)が速い。
これに対し , s   s  iでは (1)と同じ収束の速さ
で (3)の結果が得られる。

ループ積分で注意すべ き点は , 積分値が log の関数を含んだ場合である。
例えば ,
1

1
1
M2
0 M 2 x  (1  x )2 dx  M 2  2 log( 2 )で
2  0とすると数学上では で, プログラム実行時にエラーが発生すると考 えられる。ところが 2  0としてもプログラムは正常に終 了し , それらしい結果が得ら れる。
数値積分法で広義積分 と結果が一致するのは 広義積分が有限値を持つ場合で ある。二重指数関数型 積分法で
１1
0 x dxのz  重み  被積分関数の値は
x 2
z  (1  x ) log(
)   2 で 0  xでもっとも小さい値を 1 x1とすると , この値は log( )となる。この値は大き い様に1見えるが , 倍精度   10 308 なら log( ) ≒ 709

1
4倍精度   10  4931 なら log( ) ≒ 11354で例題で M  100
なら , 倍精度演算で 0.0709,4倍精度演算で1.1354とエラーが発生する様な値にはな らない。

二重指数型積分法では , x  ( t )| f (( t ))' ( t ) | exp(c exp | t |), | t | の様に減衰するとき最適の公 式が得られる。
変数変換区間 [0,1]での DEでは ,

ey
e y
y
sinh( t ), x  y,1  x  yより
2e  e ye  e y1  cosh( t )x (1  x ) , w  ' ( t ) yy 2(e  e )(e y  e  y ) 2xx
 e 2 y ,2 y  log(
),
1 x
1 x
2
  cosh( t )    sinh( t ) 2  1    ( y) 2  1

 


4y 2   2
2

4y 2   2

x(log()) 2   2
1 xx
w  ' ( t )  x (1  x ) (log()) 2   2
1 x
このため被積分関数 f ( x )の[0,1]における積分でDEが効果があるかどうか はx  ,1  (0    1)でのx
f ( x ) x (1  x ) (log()) 2   2 が十分減衰しているか
1 x
どうかを見れば良い。

Beta関数
1

 1
 1
x
(
1

x
)
dx (  0,   0)では ,0x
I  f (( t ))' ( t )  (log()) 2   2 x  (1  x ) 
1 x
1
1
x
I  ( )     (0  ,   1)
2
2
1
2
2 

x   I  (log()      log( )1x  1   I    log( )
一般に , f ( x )  x  log x (0    1,0  x  1)
n log 10x  10  n とおくと , | f ( x ) |10 n
| f ( x ) | 10  m (mは得たい精度を表わす 数 )となる
m。
m  10,   0.01なら n  1000より
指数部は15ビットが必要になる事 がわかる。またnは n 

x  1  での から,   10 1000 が必要となり , 指数部
15ビットの変数を 2つつなげた演算が必要 になる
事が容易にわかる。

具体的計算例(1) inf ra

vtx

変数変換後の fはfx(s  500 2 , m e  0.0005)

2

 sx 2 y(1  y)  x 2 m e  (1  x )2

これから | f (( t 1 ), ( t 2 ))' ( t 1 )' ( t 2 ) | を求めると ,x  , y  ,1  
1 3 (log( )) 2

17
30
 27

で
≒
1
.
51

10


10
,



2

10
2
me
 s 3   2 m e  2
x  1  , y  ,1  
1 2 (log( )) 2
2
me
15
12

で
≒
1
.
53

10



10
2
s
 s  m e  2
これより, x  , y  1  に重点を置けば良い事 がわかる。
(2) inf rabox

変数変換後の fはx (1  x )D2
D  sx 2 y(1  y)  x2  (1  x ) 2 A(z)

f

2

2

A(z)  tz 2  ( t  m f  m e )z  m e

2(s  500 2 , t  150 2 , m f  150, m e  0.0005)y  ,1－で, | f (( t 1 ), ( t 2 ), ( t 3 ))' ( t 1 )' ( t 2 )' ( t 3 ) | は
x 2
z 2
1
x 2 (1  x ) 2 z(1  z) (log())   2 (log())   2  log( )
1 x
1 z

D2
D  sx 2   x2  (1  x ) 2 A(z)22z  のとき A (z)  m e , z  1  のとき A (z)  m f から
z  のときが最大値x1x 2 (1  x ) 2 (log()) 2   2  2 (log( )) 2
1 x

2
D
D  sx 2   x2  (1  x ) 2 m e となる。
x  ,1  ともに分子は同じで
1
 4 (log( )) 3

2
D
2
2
x
D  s 3  2  m e  m e
x  1 

2

D  s  2   2 m e  2

これより1
 4 (log( )) 3z  , x  1  の場合のを考えれば良い。
4  10  m ,   10  n   10 ( n  m ) とすると ,1 4 (log( )) 3

 10  4 n ((m  n ) log 10) 3
4
m  30, n  3だと 4.387  10 7n  2.9だと1.09  10 6
これは ,10進 6桁の精度を得るには 1  で  2 109 が必要と示している .これは 4倍精度で ieee754  2008形式では10進 6桁の精度の結果が得ら れ , DD形式では1に近い数の表現に注意しないと 10進 6桁の精度が得られず ,1に近い数の表現に注意すると , ieee754  2008形式よりも精度の良い 結果が得られる事と良く対 応している。

3.4.2 massless計算(log(a )) 2 2
a  1  (log(a )) x x  ..... (a  0, | x | )2!
ここで , | x | が小として一次の項ま で取ってx

a x  1  (log(a )) xとする.
1
a
a 1
log( ) 
a
c




1

1
log( )c1
1
1
log(a )1
1  (log(a ))
a 1log(a )  log( ) 



a
a
c
a
a
a


(1) vtx
1 1 x1
0 0  sxy  ( x  y) 2 m 2  (1  x  y)2 dydxm2
2s  0,s  m   , e,
 gと正規化して ,ss
1 1 x
1
I 
dydx
2
xy

(
x

y
)
e

(
1

x

y
)
g
0 0
2

2

1 1

p
dpdq
2
p
q
(
1

q
)

p
e

(
1

p
)
g
0 0

 

2

D  g 2  4g (q (1  q )  e)  0より
1

1
A
log( )dq
2A
g
0

I

A  q (1  q )  e11
A 1dq

2 0



1
1
log( )g

1 1 x

1 1

1

1
p
1
A 1dq
dydx    2 1 dpdq 
J 
1 
2

2 0
A
0 0 p
0 0 ( xy  ( x  y ) e)
よって I  J
11

1 1
1 log(A )dqdq IJ 
2 0 A
20 A

1log( 2 )の一次式


m 2  0, e  0とすると
１

I＝ q

1

(1  q )1

０

1 ((1  )) 2
1
B(, )  2
dq 2
 (1  2)1
log( 2 )の二次式

一般公式
1

0  t  1

1
1
)
log(
2

dx
0  x 2  x  tta  0,0  c  11

0

 ax (1  x )  c
)
2
c
1
2
c
]
)
[log (dx 3aa
 ax (1  x )  clog(s  500 2 , m  0.0005,   10  2500
解析近似解

1.272295747

計算解

1.272296056

s  500 2 , m  0.0005,   10  2500
解析近似解

実数 部

 1..272276008

虚数部 0.1448318783
計算解

実 数部

 1..272276295

虚数部 0.1448319591

s  500 2 , m  0.0005
の値が充分小さくない と結果の精度が良くな い
事があります。
実 測値と解析式の相対誤 差


実数 部

虚数部

10 15

6.23  10 3

1.38  10  2

10 150

6.30  10 5

1.49  10  4

10 1500

6.68  10  7

1.59  10  6

10  2500

2.26  10  7

5.58  10  7

10 15000

6.42  10 91.61  10 8m 2  0, e  0とすると
1 1 x  ( xy)
0

0

1

1

1
1 ((1  )) 2
1

dydx   x (1  x ) dx  2
0 (1  2)１

I＝ x 1 (1  x )  dxで
０

x 2
(log(
)   2 x  (1  x )1 
1 x
1

x
| f (( t ))' ( t ) |2
21
x  1   | f (( t ))' ( t ) | log( ) 1
  0.01
  10 12 | f (( t ))' ( t ) | 2.5  10 11| f (( t ))' ( t ) |x1
| f (( t ))' ( t ) | log( ) 
  0.01
  10 300
  10  2450

| f (( t ))' ( t ) | 0.691| f (( t ))' ( t ) | 2  10  21

計算結果
case1
η=0.01 λ=1.928749847963918E-022
解析解 9998.3788657943168703256823523950561
倍精度 9328.03614724728 N=1024
倍精度(改) 9998.37886579432 N=1024
拡張倍精度 9998.37886579431643 N=7384倍精度(ieee) 9998.37886579431687032567867936706 N=7724倍精度(dd) 9362.3766317791673695811894621156875 N=1184

case2
η=7.238241365054197E-003 λ=1.000000000000000E-030
解析解 19085.205538478496689369143315971833
倍精度 15919.3830485798
N=1024倍精度(改) 19085.2055384785 N=1024
拡張倍精度 19085.2055384784966e N=7384倍精度(ieee) 19085.2055384784966414971500589276 N=7724倍精度(dd) 16034.373970898940278247867121568078 N=1184

解析結果
1 1 x1
1 ((1  )) 2
0 0 ( xy)1 dydx   2 (1  2)
1 1 x

1

1
1
1

0 0 ( xy)1 dydx   0 x (1  x ) dx
f ( x )  x 1 (1  x )  とすると
x(log()) 2   2 x  (1  x )1 
1 x
1x  , | f (( t ))' ( t ) |   log( )となり, 表現可能な の
値に精度が依存する事 になる。

| f (( t ))' ( t ) |

1 1 x11((1  )) 21
0 0 ( xy) dydx    1 0 x (1  x ) dx  (1  2)(2  2)(1  2)

1 1 x

1 1 x1(1  2)(2  2)より  
dydx 
( xy)  dydx
1 
2



0 0 ( xy )0 0
1 1 x  ( xy) dydxはx  p(1  q), y  pqで変数変換して ,0

01 1 x1  ( xy) dydx   P
0

1

0



0

1 2 

1

dp   q  (1  q )  dq
0

1 2 
 P dpでは, | f (( t ))' ( t ) |
0

x
(log(
)) 2   2 x 2   (1  x )1 xx  ,1  で | f (( t ))' ( t ) | は十分小さくなる。
1

1


1
q
(
1

q
)
dq
では
,
x


,
1


で
|
f
(

(
t
))

'
(
t
)
|


log(
)
0

となり, 倍精度変数で の値は取れる。

一般化
1 1 x1 x1 ........ x n 1

  ........
0



0

0

1
dx n dx n 1 .....dx 1( x 1 x 2 .....x n )1

1 ((1  )) n n (1  n)((1  )) n
1  n
( m )  1
exp[  ( 1) m( n  n m ) m ]n
(1  n)m
(1  )m2((1  )) nは以下の様にして求め る事ができる。
(1  n)a m  (-1)

mn - nm(m)m

f 0  1, f 1  0, f 2  a 2 , f 3  a 3 , f 4  a 4 

1 2
a 2 , f5  a 5  a 2a 32

1 2 1 3a3  a22
6
1 2
f 7  a 7  a 2 a 5  a 3a 4  a 2 a 32
1 2 1 211
2
4
f 8  a 8  a 2 a 6  a 3a 5  a 4  a 2 a 4  a 2 a 3 
a2
2
2
2
24

f 6  a 6  a 2a 4 n 5
a 2  10( 2), a 3  40(3), a 4  155( 4)a 5  624(5)展開式
1  10( 2)2  40(3)3 ( 155( 4)  50(( 2)) 2 )4
 (624(5)  400( 2)(3))5n6
a 2  15( 2), a 3  70(3), a 4  645
 ( 4)
2a 5  1554(5), a 6  7775(6)展開式
1  15( 2)2  70(3)3645
225
 ( 4) (( 2)) 2 )422
 (1554(5)  1050( 2)(3))5

 ( ( 7775(6) 

9675
( 2)( 4)  2450((3)) 22

1125(( 2)) 3 )6
2

n 7
1197
 ( 4)
2a 5  3360(5), a 6  19607(6), a 7  117648(7)a 2  21( 2), a 3  112(3), a 4  
展開式
1  21( 2) 2  112(3) 3
1197
441
 (
 ( 4) (( 2)) 2 ) 422
 (3360(5)  2352( 2)(3)) 5
 ( 19607(6) 25137
( 2)( 4)  6272((3)) 2
23087(( 2)) 3 ) 6
2
 (117648(7)  70560( 2)(5)  67032(3)( 4) 24696(( 2)) 2 (3)) 7n 8
a 2  28( 2), a 3  168(3), a 4  1022( 4)a 5  6552(5), a 6  131068
(6), a 7  299592(7)3a 8  2097151(8)展開式
1  28( 2) 2  168(3) 3
 ( 1022( 4)  392(( 2)) 2 ) 4
 (6552(5)  4704( 2)(3)) 5
131068 ((6)  28616( 2)( 4)  14112((3)) 2
3
21952

(( 2)) 3 ) 6
6
 ( 299592(7)  183456( 2)(5)  171696(3)( 4) 65856(( 2)) 2 (3)) 7
3669904
( 2)(6)  110736(3)(5)3 522242(( 4)) 2  400624(( 2)) 2 ( 4)  395136( 2)((3)) 2

 ( 2097151(8) 



76382
(( 2)) 4 )8
3

2
4
6
8
,  ( 4) ,  ( 6) , (8) 
6
90
945
9450
実測結果   0.01
 ( 2) 

dim5  0.998402809124425953047924794752494
dim5  0.998402809640216913557719749246
dim6  0.997616251685321110907418255207708
dim6  0.997616251653673359393296190115
dim7  0.996679645316743857719951968389312
dim7  0.996679645317184684228686126118
dim8  0.995595435273981970066429062072135
dim8  0.995595435274204817137389673495
実測結果   0.001
dim5  0.99998359870899406828716309047912985
dim5  0.999983598708994592536696357810
dim6  0.99997541008787161963796582035430388
dim6  0.999975410087871616448538275141
dim7  0.99996559096266763917441568060487196
dim7  0.999965590962667639178603625416
dim8  0.99995414374373370183843190221152245
dim8  0.999954143743733701839646669370

積分値の精度は の値により表現可能な 数値範囲に依存するため , この影響を少なくする ために以下の公式を用いて計算す る。
1 1 x1 x1 ........ x n 1

  ........



0

0



1 ((1  )) n n (1  n)
1 1 x01 x1 ........ x n 1

   ........
n

0

0

1
dx n dx n 1 .....dx 1( x 1 x 2 .....x n )1
0

1dx n dx n 1 .....dx 1( x 1 x 2 .....x n )1
1 1 x1 x1 ........ x n 1 (1  n)(2  n)....(n  n)   ........
00 ( x x .....x ) dx dx
1

0



2

n

n

n 1

.....dx 1

一般ベータ関数
I n   x1

y1 1



x2

y 2 1.......x n

y n 1dJ n   x 1 x 2 .......x n d
y1

y2

yn

(1)
( 2)



n 1

  {x i  0 |  x i  1}, d  dx n dx n 1 ......dx 1
i 1( y1 y 2 ......y n )I n  (z1z 2 ........z n )J n:nzi  i   y jj1(1)で計算する場合 ,4倍精度または 4倍精度変数を
2つつなげた 8倍精度演算が必要。(DQ形式 )(2)で計算する場合は倍精 度または DD形式の 4倍精度演算で精度の良い 結果が得られる。

( 2) box
1



log(

1
)
2

1 1 x 1 x  y

I 
0

0

1

 D

2

dzdydx

0

D  xz  y(1  x  y  z)
では ,I

1 2 2  2 1 ((1  )) 2( 
2 )
1   
4
(1  2)((1  )) 2
1  2
( k )  1 kexp[( 1) k ( 2  2 k )() ]2(1  2)k
(1  )k 2

( k )  1 k) x より。
k
k 2
a
a
5
5
I  -22  1  a 0 とすると , a  2  4, a 1  4, a 0  4(  2  1)    2  4

12
3



log (1  x )   log(1  x )  (1   ) x   ( 1) k (これは以下の事からき ています。
1

1(1  q )1  q 1
1
1
I  [
dq
]

[
dp]
1 
1   0 p1 (1  p)1
(1  q )1
0 (1  2q )q1

1
1
1 2 ((1  )) 2
から[  1dp]
。
1   0 p (1  p)1
1    (1  2)
1(1  q )1  q 1
2  2 1
0 (1  2q)q1 (1  q)1 dq    4 2 となるがより詳細には次ぎの様になります。

(1  q )1  q 1
q(1  q ) 0 (1  2q)q1 (1  q)1 dq  0 [ (1  2q)q  (1  2q)(1  q) ]dq
1

1

1

1

q   (1  q ) 
2 2(1  t )   (1  t ) 
1 
 2
dq    2 dt
1  2q
 
t
0
0

2
(  1)(  2).......(  2k )
 21 2(1  )( 2k  1)!( 2k  1)
k 1






1
1
3
)


2
(
)  2 ( 2)224
k 1 ( 2k  1)
k 1 ( 2 k  1)( )の定数項は 2(1  24( )のの一次の項は, 負の数の和の絶対値で計算す る と,

分子の項

 (2k  1)!(2k  1) .
k 1

2k

1
 ( 2k )![log( 2k  1)   ]rr 1
  0.5772156649 (オイ ラー の定数)

分子の項  ( 2k )!

( 2k )![log( 2k  1)   ]log( 2k  1)
1

 

2
2
( 2k  1)!( 2k  1)
k 1
k 1 ( 2 k  1)
k 1 ( 2 k  1)
log( 2k  1) (2k  1)k 1

A  lim23
1
1
  2 [ log A  log( 2)  log( )  ]  0.4181157384
2
6881  2 2  33  .....  n n
2n n / 2  n / 2 1 / 12
log A  0.2487544703
n  

en2
4

 1.282427130


1
1
2


[

1
]


(
 1)  0.1348956184
2
2
8
k 1 ( 2 k  1)
k 1 ( 2 k  1)

これ から( )内のの一次の項の係数は 0.5530113568
全体では  21  2を かけて2.212 2 と なり これ が一次式まで取った場合の誤差と なる 。

実測例としては以下のものがあります。
解析近似式を使用して,   10 1000000001 ~ 10 1000000045
で計算しました。ただしgamma関数はサポートされていない機種もあり,展開式より計算しています。
a2= 4.00000000000000065094194454995095
a1= 4.00000000000001117756660514563709
a0= -12.4493406676118712066445306454111
変数変換区間 SR16000 [10 100 ,1  10 100 ], x 5570とx 2670では[10  200 ,1  2 113 ]の 4倍精度では
epsilon算法の収束具合は  10 150 で
値

誤差

x5570

0 .191143398445070834D + 07

0.140633312313127879D - 05

x2670

0.191143398444577885D + 07

0.210270925388344072D - 07

sr16000 0.191143398444582948Q + 07

0.156658222921909210Q - 05

  10 -1500 で
x 55700.190895950695064375D + 09

0.919500023320091177D - 01

x2670

0.190895950598091269D + 09

0.477814929853799862D - 01

sr16000 0.190895950686589134Q + 09

0.470812225711499829Q - 01

  10 -15000 で
x 55700.190871021209837876D + 11

0.771928148538333458D + 04

x2670

0.190871096344010433D + 11

0.183847060803184605D + 03

sr16000

0.190871111864473018Q + 11

0.190871121114962521Q + 11

(注 )SR16000では他に変化の小さい 要素がある。

となっています。

あとSR16000でのある の範囲での a 2 , a 1 , a 0は
以下の様になっていま す。

ramda=10^{-16}~10~[-30}a2= 4.00000344724677645844034046757416
a1= 3.99861407619387554668745765583265
a0= -12.2639453402406000000000000020020
ramda=10^{-31}~10~[-45}a2= 4.00000099577869056434484930528836
a1= 3.99939509233538395973479843855165
a0= -12.3269905127989999999999999759162
ramda=10^{-46}~10~[-60}a2= 4.00000041264201279800645140064705
a1= 3.99966395945103743579061852729121
a0= -12.3581962818229999999999999796046

ramda=10^{-61}~10~[-75}a2= 4.00000019676625921874246419452751
a1= 3.99979490092391388668793381657374
a0= -12.3781474808740000000000001848136
ramda=10^{-76}~10~[-90}a2= 4.00000014459869220092907590162562
a1= 3.99983240900972265003642249066011
a0= -12.3848690023040000000000001603800
rambda=10^{-91}~10^{-105}
a2= 4.00000003694782985017339116664484
a1= 3.99996494892422507975597908036244
a0= -12.4154587133217526000000000000001

の値が小さくなるほど5 2
a 0は 
 4に近付いています。
3

3.4.3 bsgamma
bs planar
2

2

2

2

2

2

D  c( x 1 ( p 1  m 1 )  x 2 ( p 2  m 2 )  x 3 m 3  x 6 m 6 )
2

2

2 2
4 1

2

2

2

2

2 2

2

2

 c1 (x 5 p 2  x p  x 4 x 5 (p1  p 2 ))
 c 2 (x 2 p 2  x 1 p1  x 1 x 2 (p1  p 2 ))
2

2

2

2

 2x 2 x 3 x 5 p 2  2x 1 x 3 x 4 p1  x 3 (x 2 x 4  x 1 x 5 )(p1  p 2 )x 6  1  x1  x 2  x 3  x 4  x 5
c1  x 1  x 2  x 3
c 2  1  x1  x 2
c  x 3 (1  x 1  x 2  x 3 )  (x 1  x 2 )(1  x 1  x 2 )  x 3 (c 2  x 3 )  (1  c 2 )c 2
2

2

2

2

2

2

2

p1  m b , p 2  m s , m1  m 2  1.5 2 , m 3  80.3477 2 ,2

2

2

m b  4.7 2 , m s  0.094 2 , m 6  2
1 1 x1 1 x1  x 2 1 x1  x 2  x 3 1 x1  x 2  x 3  x 4

I 
0 0

 
0

0


0

1
dx dx dx dx dx
2 5 4 3 2 1Dbs planar
5次元用の変数変換
x 1  xy(1  z )x 1  x 2  xy

x 2  xyz

x1  x 2  x 3  x

x 3  x (1  y)c1  xx 4  (1  x ) u (1  v)

c 2  1  xy

x 5  (1  x ) uv

x 4  x 5  (1  x ) u

x 6  (1  x )(1  u )c  x[(1  x )(1  y)  y(1  xy)]J  x 2 (1  x ) 2 yu
J
J
 02
2
D
D0

J 0  (1  x ) 2 yu

2

D 0  [(1  x )(1  y)  y(1  xy)][ xym1  x (1  y) m 3
2

 (1  x )(1  u ) m 6 ]2

2

 (1  x ) 2 u 2 [(1  v) p1  vp 2 ]2 x (1  x ) y(1  y)[(1  z ) p1  zp 2222

 2u ( zvp 2  (1  z )(1  v) p1 )22

 u ( z (1  v)  (1  z ) v)(p1  p 2 )]これを更に 4次元化する。

2

4次元化後の bs(1  x ) 2 yu 5
(1  x ) 2 yu 4I d  d
2(A  B)B
D
5
4
2

2

2

2

D  [1  x  xy(1  y)][xy(1  z)(p1  m1 )  xyz(p 2  m 2 )]2

 [1  x  xy(1  y)][x (1  y)m 3  (1  x )(1  u )2 ]
22

 (1  x ) 2 u 2 [(1  v)p1  vp 2 ]22

 (1  xy) xy 2 [(1  z)p1  zp 2 ]22 x (1  x ) y(1  y)u[(2  z  v)p1  (z  v)p 2 ]22A  (p1  p 2 )(1  x )u[(1  x )u  xy(1  y)]
2

2

2

2

B  [1  x  xy(1  y)][xy(1  z)(p1  m1 )  xyz(p 2  m 2 )]2

 [1  x  xy(1  y)][x (1  y)m 3  (1  x )(1  u )2 ]222

 (1  x ) 2 u 2 p1  x (1  x ) y(1  y)u[(2  z)p1  zp 2 ]22

 (1  xy) xy 2 [(1  z)p1  zp 2 ]
2

2

2

2

2

2

2

p1  m b , p 2  m s , m1  m 2  1.5 2 , m 3  80.3477 2 ,2

2

2

m b  4.7 2 , m s  0.094 2 , m 6  2

bsgamma

planar  10  nnvalue
16 0.013945911
17 0.014826698
18 0.015707767
19 0.016590158
20 0.017473234
21 0.018355518
22 0.019238562
23 0.020130180
24 0.021035437
25 0.021938463
26 0.022811082
27 0.023666177
28 0.024552533
29 0.025458382
30 0.026327504

value0.030000000
0.025000000
0.020000000
0.015000000

value

0.010000000

0.005000000
0.000000000
16 18 20 22 24 26 28 30

  10nbsgamma 一次差分結果
n N=728
N=1024
N=1456
16 0.000881676036303 0.000881676024318 0.000881676024308
17 0.000881675957819 0.000881676024256 0.000881676024308
18 0.000881676127056 0.000881676024379 0.000881676024308
19 0.000881676192593 0.000881676024096 0.000881676024308
20 0.000881675929067 0.000881676025850 0.000881676024307
21 0.000881675933258 0.000881676026869 0.000881676024308
22 0.000881675979610 0.000881676024874 0.000881676024305

3.4.5twoloop

vertexbs planar
2

2

2

2

を変形して , p1  p 2  m 1  m 2 2  m  として ,g 2を作成する。
2

2

D   C( x 3 m 3  x 6 m 6 ) [( x 1  x 2  x 3 )( x 4  x 5 ) 2  ( x 1  x 2 ) 2 (1  x 1  x 2 ) 2 x 3 ( x 1  x 2 )( x 4  x 5 )]m 

2807m   0.1057, m 3  m 6  808

m   0.1057, m 3  M z  91.19, m 6  809

m   0.1057, m 3  M Ｈ  125.5, m 6  

変数変換
x 1  xy(1  z )x 1  x 2  xy

x 2  xyz

x1  x 2  x 3  x

x 3  x (1  y)c1  xx 4  (1  x ) u (1  v)

c 2  1  xy

x 5  (1  x ) uv

x 4  x 5  (1  x ) u

x 6  (1  x )(1  u )c  x[(1  x )(1  y)  y(1  xy)]J  x 2 (1  x ) 2 yu
J
J
 02
2
D
D0

J 0  (1  x ) 2 yu

g2
2

2

D 0  [(1  x )(1  y)  y(1  xy)][ x (1  y) m 3  (1  x )(1  u ) m 6 ] [ xy 2 (1  xy)  (1  x ) 2 u 2  2 x (1  x ) y(1  y) u ]m 

22と3次元化して m   1と正規化して計算する 。
15の  10  n (16  n  30)に対して計算し ,807は二次差分を ,808,809
は一次差分をとり , log(1)の二次式 , 一次式となる事を検証 。
2

g2 808一次差分
N=2912で10進11桁(全桁)一致n

g2 808 一次差分
N=728
N=1456
N=2912
16 0.000038733241268 0.000038733248507 0.000038733248507
17 0.000038733234064 0.000038733248507 0.000038733248507
18 0.000038733266331 0.000038733248507 0.000038733248507
19 0.000038733326864 0.000038733248507 0.000038733248507
20 0.000038733174086 0.000038733248506 0.000038733248507
21 0.000038733353236 0.000038733248506 0.000038733248507
22 0.000038733544343 0.000038733248510 0.000038733248507
23 0.000038731850666 0.000038733248499 0.000038733248507
24 0.000038732625311 0.000038733248540 0.000038733248507
25 0.000038734471892 0.000038733248546 0.000038733248507
26 0.000038731896322 0.000038733248424 0.000038733248507
27 0.000038724804488 0.000038733248425 0.000038733248507
28 0.000038738845469 0.000038733249112 0.000038733248507
29 0.000038740941884 0.000038733249063 0.000038733248507

g2 809一次差分
N=2912で10進11桁(全桁)一致n

g2 809一次差分
N=728
N=1456
N=2912
16 0.000021493099355 0.000021493102821 0.000021493102821
17 0.000021493093286 0.000021493102821 0.000021493102821
18 0.000021493115638 0.000021493102821 0.000021493102821
19 0.000021493142483 0.000021493102821 0.000021493102821
20 0.000021493082925 0.000021493102821 0.000021493102821
21 0.000021493117428 0.000021493102821 0.000021493102821
22 0.000021493302447 0.000021493102823 0.000021493102821
23 0.000021492401183 0.000021493102817 0.000021493102821
24 0.000021492629664 0.000021493102839 0.000021493102821
25 0.000021493704654 0.000021493102844 0.000021493102821
26 0.000021492781172 0.000021493102769 0.000021493102821
27 0.000021488800060 0.000021493102792 0.000021493102821
28 0.000021495565878 0.000021493103126 0.000021493102821
29 0.000021496356431 0.000021493103198 0.000021493102821

3.4.4 4次元積分
4次元積分として , 以下の S 221を検討した。
S

221

1 1 x 1 x  y 1 x  y  z(s; m1 , m 2 , m 3 , m 4 , m 5 )   
2

2

2

2

2

0

0


0


0

1
dudzdydxDCC  ( x  y  z  u )(1  x  y  z  u )  ( x  y)(z  u )E  (1  x  y  z  u )( x  z)( y  u )  ( x  y)zu  (z  u ) xy
2

2

2

2

M 2  xm1  ym 2  zm 3  um 4  (1  x  y  z  u )m 5

2

D  sE  M 2 C
積分領域は で表す。
精度チェックのための 基本演算の解析解と計 算結果は下記のとうり 。
解析解
case1

S 221 (1;0,0,0,0,0)  6(3)  7.21234141895756568

case2

S 221 (1;0,0,0,0,0)  6(3)  7.21234141895756568

case3

S 221 (0;1,1,1,1,0)  1

case4

case5

2
 1.644934066848226436472
6
ln x (1  x )
 1)
1 x(
2
S 221 (0;1,1,1,1,1)    1  x  x 2
 0.78130241289648
1

x

x
0
S 221 (0;0,1,1,1,0)  (2) 

計算解
case1 result = - 7.21234141895747793
case2

result = 7.21234141895747793

case3 result = 0.999999999999997335
case4 result = 1.64493406684822352
case5 result = 0.781302412896484055

特殊な計算例は以下の とうりs  1, m 2  1

S 221 (s; m 2 ,0,0,0,0)

解析解  3(3)  3.6061707947878284

s  1, m 2  1
計算解（ε - 算法なし）= 3.60617070947862928
計算解(y = -s/m^2 ) = 0.107304279342169218 + 162 0.363 + 149
計算解(x = -m^2/s ) = 0.360617069930836021D + 01 0.257D - 10s  1, m 2  1
計算解（ε - 算法なし）= + * * * * * * * * * * * * * * * * * * * * * * * * *計算解(y = -s/m^2 ) = - 0.200213566025878417 + 162 0.261 + 162
計算解(x = -m^2/s ) = - 0.360611840086198487D + 01 0.161D - 04
尚以下の式は発散する 。
1ln(1  x )
dx
2
x
0

S (0;0,0,1,1,0)   
221

1

1

ln(1  x )1
S (0;0,0,0,1,0)   
dx

dy
2

y
x
0
0
221

さらに詳細に解析する と
1 1 x 1 x  y 1 x  y  z

S (s; m ,0,0,0,0)   
221

2

0 0





0

0

1
1
dudzdydx  
d
 CD
CD

D  sE  xm C s  1, m  1  D  E  xC
1
1
x


CD CE E(E  xC)
1
 CEd  6(3)  実測結果 result = 7.21234141895747793
x E(E  xC)d  3(3)  実測結果 result = 3.60617070947877494
1
221
 C 2 d  S (0;1,1,1,1,1)  実測結果 result = 0.781302412896484055
1
1
d


 C 6  実測結果 result = 0.1666666666666665747
1
11
d


(3)  実測結果 result = 1.46918065941728160
 E
9
S 221 (1;1,1,1,1,1) result = 0.680876150423024185
2

S 221 (1;1,1,1,1,1)

2

result = 0.923631826519863641

尚, 以下の場合は発散する 。
11ln(1  x )d xC 2
0 1  x dx

m 2  50と固定し , sを5,10,15,...,100と変化させたときの二 重指数関数型積分の積分結果は以 下のとうり.二重指数関数型積分の サイズは n  1024,   算法
での の初期値は1.2 15

s

積分値

積分誤差

5 0.142933799933564448D+00

0.864D-06

10 0.874770761633400645D-01

0.495D-04

15 0.580754664645131494D-01

0.123D-04

20 0.377928730070589761D-01

0.444D-06

25 0.218394449541923323D-01

0.113D-03

30 0.773426500218699423D-02

0.511D-05

35 -0.570584090394023972D-02

0.168D-04

40 -0.198070245154403426D-01

0.125D-06

45 -0.369777391702599822D-01

0.383D-04

50 -0.721357975080699487D-01

0.333D-05

55 -0.989714303728692485D-01

0.183D-05

60 -0.102840960786493671D+00

0.923D-06

65 -0.101914671715339883D+00

0.185D-06

70 -0.991426158779158406D-01

0.275D-05

75 -0.956234990743119462D-01

0.330D-06

80-0.917437470939333061D-01

0.166D-03

85 -0.877601377899803609D-01

0.324D-03

90 -0.848383517809659093D-01

0.422D-03

95 -0.814307958141493454D-01

0.575D-03

100 -0.772582671356500383D-01

0.244D-04

さらに複雑にした場合 のS 221 (s; m 2 , m 2 ,0,0, m 2 )計算

m 2  100と固定し , sを5,10,15,...,100と変化させたときの二 重指数関数型積分の積分結果は以 下のとうり.二重指数関数型積分の サイズは n  1024,   算法
での の初期値は1.2 15
s

積分値

積分誤差

5 0.300712966040258094D-01

0.109D-06

10 0.266732292799612622D-01

0.190D-08

15 0.247016613851884191D-01

0.650D-07

20 0.233122764814605241D-01

0.675D-08

25 0.222407169211559473D-01

0.115D-08

30 0.213693658017195323D-01

0.797D-08

35 0.206356969069118670D-01

0.595D-06

40 0.200022139331064649D-01

0.213D-08

45 0.194450748750508856D-01

0.746D-08

50 0.189476908636132708D-01

0.297D-06

55 0.184989268813915836D-01

0.590D-07

60 0.180892912432545142D-01

0.141D-06

65 0.177130501245448997D-01

0.180D-07

70 0.173657808967870794D-01

0.151D-06

75 0.170416613603775538D-01

0.238D-06

80 0.167388840687285462D-01

0.178D-06

85 0.164538254991338749D-01

0.124D-07

90 0.161835797358507970D-01

0.157D-05

95 0.159310205462771873D-01

0.299D-06

100 0.156908197395704066D-01

0.802D-06

3.4.6 3loop積分
今回扱った 3loop積分は
1 1 x1 1 x1  x 2 1 x1  x 2  x 3 1 x1  x 2  x 3  x 4 1 x1  x 2  x 3  x 4  x 5 1 x1  x 2  x 3  x 4  x 5  x 6

 
0

0

0


0


0


0

d  dx 7 dx 6 dx 5 dx 4 dx 3 dx 2 dx 1
でDが次の 3ケース ( N0, L01, L02)の場合です。
解析解は 20 5  20.7385551028673985...。


0

1
d
2
D

N0
x8=1.0-x1-x2-x3-x4-x5-x6-x7
x23=x2+x3
x24=x2+x4x247=x2+x4+x7
x1235=x1+x2+x3+x5
x2346=x2+x3+x4+x6
x5678=x5+x6+x7+x8

D = x24**2*x5**2-x2346*x247*x5**2-x2**2*x2346*x5678+2*x2*x23*x24*x5678-x1235*x24**2*x5678-x23**2*x247*x5678+x1235*x2346*x247*x5678+2*x2*x24*x5*x6-2*x23*x247*x5*x6+x2**2*x6**2-x1235*x247*x6**2+2*x2*x2346*x5*x7-2*x23*x24*x5*x7+2*x2*x23*x6*x7-2*x1235*x24*x6*x7+x23**2*x7**2-x1235*x2346*x7**2

L01
x8=1.0-x1-x2-x3-x4-x5-x6-x7
x123=x1+x2+x3
x167=x1+x6+x7
x348=x3+x8+x4
x2578=x2+x5+x7+x8

D = -x167*x2578*x3**2-x167*x2**2*x348-x1**2*x2578*x348+x123*x167*x2578*x348-2*x1*x2*x348*x7+x3**2*x7**2-x123*x348*x7**2-2*x167*x2*x3*x8-2*x1*x3*x7*x8+x1**2*x8**2-x123*x167*x8**2

L02
x8=1.0-x1-x2-x3-x4-x5-x6-x7
x123=x1+x2+x3
x1456=x1+x4+x5+x6
x678=x6+x7+x8
x247=x2+x4+x7

D = x2**2*x6**2-x123*x247*x6**2-x1456*x2**2*x678-x1**2*x247*x678+x123*x1456*x247*x678-2*x1*x2*x4*x678-x123*x4**2*x678-2*x1*x2*x6*x7-2*x123*x4*x6*x7+x1**2*x7**2-x123*x1456*x7**2

L02での計算手順例(1) h=0.2 -2.4<t<2.4
とてつもない結果
3.483422085800906E-008 0.999999965165779

n=25
h= 0.200000000000000
result= 1463432.72964967(2) hの目安をつける。
10進5桁まで一致
h=4.35/24
1.182374362483700E-006 0.999998817625638
n=25
h= 0.181250000000000
result= 20.7383465515643

(3)分点数を増やす
効果なし.

1.182374362483700E-006 0.999998817625638
n=41
h= 0.108750000000000
result= 20.7294941799961(4)解析値を超すhを求める。
h=4.36/24
1.102223275929152E-006 0.999998897776724
n=25
h= 0.181666666666667
result= 20.7866519155393(5) h=(4.351から4.359)/24で計算
結果が解析値のまわりを振動する。

演算精度を上げて(1)から(5)まで行い
振動しなければ二分法により精度を上げていくことが出来る。

3.4.7 4loop積分
今回扱った 4loop積分は
1 1 x1 1 x1  x 2 1 x1  x 2  x 3 1 x1  x 2  x 3  x 4 1 x1  x 2  x 3  x 4  x 5 1 x1  x 2  x 3  x 4  x 5  x 6 1 x1  x 2  x 3  x 4  x 5  x 6  x 7

 
0 0

0


0


0


0

1
d  dx 8 dx 7 dx 6 dx 5 dx 4 dx 3 dx 2 dx 1
CD
でC, Dが次の 2ケース (M 44 , M 45 )の場合です。
解析解は441 7
M 44 は 55.5852539156784...。
8
2
M 45 は36 3  52.017868743610.....。


0


0

M 44
x9 = 1.0 - x1 - x2 - x3 - x4 - x5 - x6 - x7 - x8
x123 = x1 + x2 + x3
x167 = x1 + x6 + x7
x278 = x2 + x7 + x8
x349 = x3 + x4 + x9
x589 = x5 + x8 + x9
C = x167 * x278 * x349 * x589 - x349 * x589 * x7 * *2 - x167 * x349 * x8 * *2
- x167 * x278 * x9 * *2 + x7 * *2 * x9 * *2
D = - (x167 * x278 * x3 * *2 * x589) - x167 * x2 * *2 * x349 * x589
- x1 * *2 * x278 * x349 * x589
+ x123 * x167 * x278 * x349 * x589 - 2 * x1 * x2 * x349 * x589 * x7
+ x3 * *2 * x589 * x7 * *2
- x123 * x349 * x589 * x7 * *2 + x167 * x3 * *2 * x8 * *2
+ x1 * *2 * x349 * x8 * *2
- x123 * x167 * x349 * x8 * *2 - 2 * x167 * x2 * x3 * x8 * x9
- 2 * x1 * x3 * x7 * x8 * x9
+ x167 * x2 * *2 * x9 * *2 + x1 * *2 * x278 * x9 * *2
- x123 * x167 * x278 * x9 * *2
+ 2 * x1 * x2 * x7 * x9 * *2 + x123 * x7 * *2 * x9 * *2

M 45
x9 = 1.0 - x1 - x2 - x3 - x4 - x5 - x6 - x7 - x8
x123 = x1 + x2 + x3
x1567 = x1 + x5 + x6 + x7
x279 = x2 + x7 + x9x349 = x3 + x4 + x9
x589 = x5 + x8 + x9
c = -x279 * x349 * x5 * *2 + x1567 * x279 * x349 * x589 - x349 * x589 * x7 * *2
- 2 * x349 * x5 * x7 * x9 - x1567 * x279 * x9 * *2 - x1567 * x349 * x9 * *2
+ x5 * *2 * x9 * *2 - x1567 * x589 * x9 * *2 + 2 * x5 * x7 * x9 * *2
+ x7 * *2 * x9 * *2 + 2 * x1567 * x9 * *3
d = x279 * x3 * *2 * x5 * *2 + x2 * *2 * x349 * x5 * *2 - x123 * x279 * x349 * x5 * *2
- x1567 * x279 * x3 * *2 * x589 - x1567 * x2 * *2 * x349 * x589
- x1 * *2 * x279 * x349 * x589 + x123 * x1567 * x279 * x349 * x589
- 2 * x1 * x2 * x349 * x589 * x7 + x3 * *2 * x589 * x7 * *2
- x123 * x349 * x589 * x7 * *2 + 2 * x1 * x279 * x3 * x5 * x9
- 2 * x1 * x2 * x349 * x5 * x9 + 2 * x2 * x3 * x5 * *2 * x9 - 2 * x1567 *x2 * x3 * x589 * x9
+ 2 * x2 * x3 * x5 * x7 * x9 + 2 * x3 * *2 * x5 * x7 * x9 - 2 * x123 * x349 *x5 * x7 * x9
- 2 * x1 * x3 * x589 * x7 * x9 + x1567 * x2 * *2 * x9 * *2 + x1 * *2 * x279 * x9 * *2
- x123 * x1567 * x279 * x9 * *2 + 2 * x1567 * x2 * x3 * x9 * *2
+ x1567 * x3 * *2 * x9 * *2 + x1 * *2 * x349 * x9 * *2 - x123 * x1567 * x349 * x9 * *2
+ 2 * x1 * x2 * x5 * x9 * *2 - 2 * x1 * x3 * x5 * x9 * *2 + x123 * x5 * *2 * x9 * *2
+ x1 * *2 * x589 * x9 * *2 - x123 * x1567 * x589 * x9 * *2 + 2 * x1 * x2 * x7 * x9 * *2
+ 2 * x1 * x3 * x7 * x9 * *2 + 2 * x123 * x5 * x7 * x9 * *2 + x123 * x7 * *2 * x9 * *2
- 2 * x1 * *2 * x9 * *3 + 2 * x123 * x1567 * x9 * *3

M 44 計算作業手順(1)DEの積分変数変換区間 [,1  ]でそれぞれの演算精度 に応じて   2 53 (倍精度 ),   2 65 (拡張倍精度 ),  2 113 (4倍精度 )で計算。全て nanまたはゼロ割りが発生 。
(2)  2 53 (4倍精度 )(h  0.25), 倍精度 t  2.4からt  2.4,
拡張倍精度 t  2.6からt  2.6(h  0.2)にして実行 .
結果の精度は倍精 度演算 , 拡張倍精度 ,4倍精度の順 .
すなわち計算結 果精度は ではなく hに大きく依存する .
(3) t  2.6からt  2.6倍精度演算で実行する と, nanが発生。
これは  2  2 53 の条件を満たしていな い。
(4)  2  m  l (m  仮数部のビット数 , l 

乗算の最高次数  1
 2)
2

より ,   2 50 , tの絶対値最大値 t max
1log( )
 , h  t /( N  1)
t max  sinh 1max
2
t max  2.408897648

N  2n  1

1t6(5)(7) dtよりt (7 ) 0 e  11t6(log(x )) 6
0 e t  1dt  0 1  x dxからDEの項が最大になるのは ,t  1.5付近。これにより , Nを固定し hをわずかずつ増やす
方法を取る事 にした。

(6) t  2.4からt  2.44までの最適値付近では tの絶対値の単調増加の値を示すので , ( Nを固定して )一桁ずつ計算。一回の 計算で200  9秒  30分で 5回実行すれば , 求める値にたどりつく .
（1回で一桁精度が向上し , t  2.4で10進 5桁の精度があったので 。)
E5  2670, Phi5110Pでの実行結果はすべて 10進11桁一致
E5  2670 FORTRAN

倍精度 ,4倍精度

C

倍精度 , 拡張倍精度

Phi5110P

FORTRAN

倍精度

C

倍精度 , 拡張倍精度

最適化オプションは倍 精度, 拡張倍精度
 O1  openmp (mmic)
4倍精度は  O1  openmp,O2  openmp

M 45 計算作業手順(1) tのとる範囲を M 44 の結果をもとに t  2.4からt  2.4までと想定した。
(2)3loopの３ケースと M 45 では倍精度演算では nan, inf inityが発生したため , 拡張倍精度 , DD形式 4倍精度演算を使用。
(3)時間を考え ,3loopはSR16000で 4倍精度で実行。
(4)M 45 はE5  2670 拡張倍精度で実行し , 最終結果を SR16000
で 4倍精度で実行。
(5)M 45 の結果は10進10桁一致 ,11桁一致 ,12桁一致する結果が得られた。(次ページ )(6)M 45 4倍精度は SR16000で実行中。
10進8桁一致する場合の実行 時間は 6374秒。
(7)オプション E5  2670
SR16000
サイズは全て N  25

 O1  openmp (16smp)
 O3

 omp

(64smp)拡張倍精度実行結果
E5 - 2670 16smp
拡張倍精度

N = 25h = 4.785582/24

10桁一致

x30[0] = 3.94884181984283738570e - 08
x30[n - 1] = 9.99999960511581801627e - 01
size =

25

h = 1.99399250000000000e - 01
eps = 0.00000000000000000e + 00
result = 5.20178687487030232e + 01
wallclock time = 688.694365 sec

h = 4.7855815/24 11桁一致
x30[0] = 3.94885893249829809237e - 08
x30[n - 1] = 9.99999960511410675059e - 01
size =
25
h = 1.99399229166666664e - 01
eps = 0.00000000000000000e + 00
result = 5.20178687432160053e + 01
wallclock time = 688.619698 sec
h = 4.78558154/24 12桁一致
x30[0] = 3.94885756348328792603e - 08
x30[n - 1] = 9.99999960511424365172e - 01
size =
25
h = 1.99399230833333330e - 01
eps = 0.00000000000000000e + 00
result = 5.20178687436344852e + 01
wallclock time = 688.669904 sec

3 . 4 .8

3loop,4loop積分作業手順まとめ

DEでの変数変換区間 [,1  ]として , それに対応するtの値をそれぞれ  t max , t max ( t max  0)とする。
t max  sinh 1 (

log

1
)


目安としては  ≒ 10 8 は t max  2.4,  ≒ 10 14 は t max  3.0
倍精度演算で t max  2.4でまともな値が出たの は , M 44 のみ.
M 45 , N 0, L0(case1), L0(case2)では nan , inf inity,10 6 という結果。
ここで t max の値を小さくしてある 程度妥当な値がでた場 合のt max の最大値は L0(case2)のt max  2.175.
ここで , L0(case2)で t max  2.175を少しずつ変えてより 良い精度の結果を求めよう とすると , 値が解析値のまわりを 振動する.M 44 の場合は解析値付近で は単調増加になってい てより良い精度が求めら れた。
L0(case2)の振動の原因を調べる ため ,1(log( x )) 4
dxで調査した。(この式だと倍精度演算
1 x0 (5) (5)  

でも , t max  3.0でも計算可能なため )t max  2.175の場合
i, f =

1 0.574525277298576

i, f =
i, f =

2 2.33890219974041
3 6.50618703452639

i, f =

4 13.1385801164068

i, f =

5

20.2148537888924

i, f =
i, f =

6
7

24.5889578882103
24.2698704409496

i, f =

8 19.7433827737906

i, f =

9 13.3148078153239

i, f =
i, f =

10 7.42443710006566
11 3.40083347814004

i, f =

12 1.32758205957475

i, f =

13 0.362594924949645

ここで , iとfは,t  ( 13  i)  h , と ( 26  i)  hの場合の被積分関数  重み係数の和。
i  13はt  0の場合の被積分関数  重み係数。
ここでわかるとうり , i  13の場合が最小値となっ て
いて , i  1の場合の方が大きい .これは , DEの持つ ,tの絶対値が大きいとこ ろでは , 急速に値が減少すると言う性質を表わ していない。 t max を変更せず
N  25をN  41に変えても L0(case2)の精度は改善されなかった事や , t max を少しかえても値が振 動するのはこれが原因による .
t max  2.4だと i, f =1 5.000558641136945E - 002

t max  3だと

i, f =

1 6.662369942894555E - 007

i, f =

2 2.085441754992190E - 004

i, f =

3 1.398608912042730E - 002

result =kai =24.8862661234180

24.8862661234409

i  13の値は t max の値 ( Nを固定した場合の h )には依存しない。
(log(2)) 4
2すなわち , t max の値を少しずつ変えて 精度の良い結果を得るには , t   t max , t max の場合の被積分関数  重み係数の和が十分小 さいかどうか( t  0での値は一定なのでこ の時の値に対して )が判断基準となります 。

以上の結果から次の手 順が考えられます。
(1)倍精度演算で   10 8 に対する , t max  2.4, N  25
, h  0.2で妥当と思われる値が 出るかどうか ?
(  10 8 は 2  10 16 から来ている。)( 2)妥当と思われる値が出 てれば , t max のとき被積分関数  重み係数の和が十分小 さい
値かどうか ?
この方法には以下 の3つの方法が考えられる 。
(ア )同じ t max で N  25から N  41で精度改善があるか調べる .( DEの性質を満たしていれ ば , 誤差1 / 4以下となる .)（演算量から言 うと 6次元積分 N  128に対し ,7次元積分 N  64,8次元積分  38で同じになり倍精度の 6次元積分の N  128での実行時間はわかっている ため。）
( イ )全次元 t  0の場合の被積分関数  重み係数の値と , 一次元だけ , t max, とt max のでの被積分関数  重み係数の値の和を比 較.
( n次元積分の場合 , それぞれ n個の値を比較して , 問題のある次元はない かチェック )(ウ )全次元 t  0の場合の被積分関数  重み係数の値と最小値 (被積分関数  重み係数の値 )及び全次元 t  0の場合の被積分関数  重み係数の値  閾値
以下となる被 積分関数  重み係数の個数でチェ ック
する .
(3)問題のある場合は拡張 倍精度 ,4倍精度で ,h  0.2, N  25, t max  2.4から始める。可能なら ば拡張倍精度
から始めるのが良 い.
今回の実行では 3loopは 4倍精度演算 , M 45 は拡張倍精度で結果が得られ , M 45 のE5  2670(16smp)拡張倍精度演算
は SR16000 / M1 (64smp) 4倍精度演算の10倍 (実行時間
1 / 10)となっている。

3.5 量子モンテカルロ法による物性スペクトル計算

量子モンテカルロ法で は, 安定した計算を行うには LはL  20   の関係が満たされている事が必 要になります。
U  10ではプログラムを実行 中に変数がとる値の最大値 は
e 200  となり, 指数部が11ビットだと  50.14,15ビットだと   802.978となります。
また粒子数 Nは100程度が必要となります 。
結果の精度からみると Nが 20でも100とはあまり差は出ていませ ん。

3.5.1 指数部が11ビットの限界について
量子モンテカルロ法  可変, L  20  , U  10
倍精度の最大値 1.7976931348623157E  0308
l=980
beta = 49.0000000000000000000000000000000

ispin =1

QDR = 1.785864139172194426883476026500380E + 0303
I + QDR = 1.785864139172194426883476026500380E + 0303

ispin =2
QDR = 1.016926695260681160826110539840909E + 0304
I + QDR = 1.016926695260681160826110539840909E + 0304
l=1000
beta = 50.0000000000000000000000000000000
spin =1
QDR = 2.862712242794756891557352959861312E + 0309
I + QDR = 2.862712242794756891557352959861312E + 0309

ispin =2
QDR = 1.651604970572476540194914368437567E + 0310
I + QDR = 1.651604970572476540194914368437567E + 0310

  49.6

OK,   49.7

NG

l =

992

beta =

49.6000000000000000000000000000000

ispin =1

QDR = 9.436165589130513054255944449201070E + 0306
I + QDR = 9.436165589130513054255944449201070E + 0306
ispin =2

QDR = 5.415677874136262454106581048120274E + 0307
I + QDR = 5.415677874136262454106581048120274E + 0307
l =994
beta = 49.7000000000000000000000000000000
ispin =QDR

1

= 3.938127827940264707461183253762617E + 0307

I + QDR = 3.938127827940264707461183253762617E + 0307
ispin =2
QDR = 2.263160356072587757351389081936904E + 0308
I + QDR = 2.263160356072587757351389081936904E + 0308
l =beta =
ispin =

996
49.8000000000000000000000000000000
1QDR = 1.643557573940390154190429173363875E + 0308
I + QDR = 1.643557573940390154190429173363875E + 0308
ispin =2
QDR = 9.457544390246214860747340892604545E + 0308
I + QDR = 9.457544390246214860747340892604545E + 0308
l =beta =

998
49.9000000000000000000000000000000

ispin =1
QDR = 6.859317195871621637905579118618770E + 0308
I + QDR = 6.859317195871621637905579118618770E + 0308
ispin =2
QDR = 3.952228149418581143027425919256247E + 0309
I + QDR = 3.952228149418581143027425919256247E + 0309

3.5.2 必要演算精度見積り
 と必要演算精度は以下 の様に見積もれます。

限界テスト

 / L  a (一定), U  10とする。(0  a  1)
  10, L  200, a  0.05
UL  U
 10 2  a

必要ビット数  ln 2 e 10 2 β

10 2  
 20.40278893  
ln 2

必要ビット数 必要精度
100
2041
180
3673
250
5101
300
6121
400
8162
500
10202
600
12242
700
14282
800
16323

L
68
128
188
248
308
368
428
488
548

2000
3600
5000
6000
8000
10000
12000
14080
16000

548倍精度演算では限界に 近い事もあり  800では他の演算精度の場 合と異なる結果となる。
n=
8
l=
16000
beta, dt= 800.000000000000000000000000000000
5.00E-02
u, t0= 10.0000000000000000000000000000000
1
fac= 0.707106781186547524400844362104849
For denominator: mstep=
1
For green func: mstep2=
1
ispin=1
QDR = 1.139496372849864562316010744931467E+4911
I+QDR= 1.139496372849864562316010744931467E+4911
ispin=2
QDR = 4.844939675477544024337394425713083E+4911
I+QDR= 4.844939675477544024337394425713083E+4911
0.0000000000D+00 -0.2500000000D+00
n=
8
l=
15744
beta, dt= 800.000000000000000000000000000000
5.08E-02
u, t0= 10.0000000000000000000000000000000
1
fac= 0.712832435640251297431604369230305
For denominator: mstep=
1
For green func: mstep2=
1
ispin=1
QDR = 3.370632870032681675889173302418504E+4929
I+QDR= 3.370632870032681675889173302418504E+4929
ispin=2
QDR = 1.056365855103840555758447648870414E+4931
I+QDR= 1.056365855103840555758447648870414E+4931
0.0000000000D+00 -0.1000000000D+01

各演算精度でまとめた結果は下記のとうり。
調査一覧
項番
β
N
1
100
2
180
3
250
4
300
5
400
6
500
7
600
8 10000/14
9
770
10
780
11
790
12
800
13
800

備考

U
100
100
100
100
20
20
20
20
20
20
20
20
20

最小値一覧
理論値
1 0.567D-307
2 0.120D-552
3 0.130D-767
4 0.370D-921
5 0.299D-1228
6 0.143D-1536
7 0.686D-1845
8 0.296D-2197
9 0.176D-2366
10 0.246D-2401
11 0.159D-2426
12 0.625D-2432
13 0.688D-2437

L
10
10
10
10
10
10
10
10
10
10
10
10
10

実測値
0.330D-307
0.512D-553
0.468D-768
0.121D-921
0.927D-1229
0.424D-1537
0.199D-1845
0.849D-2198
0.705D-2367
0.703D-2402
0.454D-2427
0.179D-2432
0.197D-2437

精度2000 68倍精度
3600 128倍精度
5000 188倍精度
6000 248倍精度
8000 308倍精度
10016 368倍精度
12032 428倍精度
14336 488倍精度
15424 548倍精度
15680 548倍精度
15808 548倍精度
15680 548倍精度
15744 548倍精度

論文リスト
外国雑誌
1.Acceleration of Feynman loop integrals in high-energy physics
on many core GPUs,F. Yuasa, T. Ishikawa, N. Hamaguchi,
T. Koike N. Nakasato. 2013. 7 pp.
J.Phys.Conf.Ser. 454 (2013) 012081
2.Numerical Computation of Two-loop Box Diagrams with Masses,
F. Yuasa , E. de Doncker, N. Hamaguchi, T. Ishikawa, K. Kato,
Y. Kurihara, J. Fujimoto, Y. Shimizu . Dec 2011. 20 pp.
Comput.Phys.Commun. 183 (2012) 2136-2144
国際会議プロシーディング
3.Numerical approach to multi-loop integrals
K. Kato, E. de Doncker, N. Hamaguchi, T. Ishikawa, T. Koike,
Y. Kurihara, Y. Shimizu, F. Yuasa. Jan 2012. 7 pp.
PoS QFTHEP2011 (2013) 029
4. The GRACE project: QCD, SUSY, multi-loop
J. Fujimoto et al.. 2011. 8 pp,
PoS RADCOR2011 (2011) 012
5. Numerical Approach to Calculation of Feynman Loop Integrals
F. Yuasa, T. Ishikawa, Y. Kurihara, J. Fujimoto, Y. Shimizu,
N. Hamaguchi , E. de Doncker , K. Kato . Sep 2011. 8 pp.
PoS CPP2010 (2010) 017
6. Recursive box and vertex integrations for the one-loop hexagon
reduction in the physical region
E. de Doncker, J. Fujimoto, N. Hamaguchi, T. Ishikawa,
Y. Kurihara, Y. Shimizu, F. Yuasa. 2010. 9 pp.
PoS ACAT2010 (2010) 073

7.Numerical Evaluation of Feynman Integrals by a Direct
Computation Method F. Yuasa, T. Ishikawa, J. Fujimoto,
N. Hamaguchi, E. de Doncker, Y. Shimizu. Apr 2009. 5 pp.
PoS ACAT08 (2008) 122
8. Status reports from the GRACE group
Y. Yasui, T. Ueda, E. de Doncker, J. Fujimoto, N. Hamaguchi,
T. Ishikawa, Y. Shimizu. F. Yuasa Oct 2007. 6 pp.
eConf C0705302 (2007) LOOP04 (LCWS-2007-LOOP04)
9. Precise Numerical Evaluation of the Scalar One-Loop
Integrals with the Infrared Divergence、
F. Yuasa E. de Doncker, J. Fujimoro, N. Hamaguchi, T. Ishikawa,Y. Shimizu Sep 2007. 10 pp. PoS ACAT2007 (2007) 087
国内研究会論文等
10. Feynman 4-loop積分のチューニングとアクセラレータを使用した高速化
濱口 信行、石川 正、湯浅 富久子
情報処理学会 研究報告ハイパフォーマンス
コンピューティング（HPC）,2015-HPC-152,(2015-12-16)
11. ファインマンループ積分によるアクセラレータの精度，
性能評価
濱口 信行 , 石川 正
情報処理学会研究報告ハイパフォーマンス
コンピューティング（HPC）,2014-HPC-147(8),1-8 (2014-12-02)
12. 量子モンテカルロ法による物性スペクトル計算への多倍長演算の適用
濱口信行 , 石川正 , 岩野薫
情報処理学会研究報告ハイパフォーマンス
コンピューティング（HPC）,2013-HPC-141(12),1-7 (2013-09-23）

13. SR16000/M1とBG/Qにおける性能，精度評価
濱口信行 , 石川正情報処理学会研究報告ハイパフォーマンス
コンピューティング（HPC）,2012-HPC-136(6),1-5 (2012-09-26)
14. 数値解析から見たファインマンループ積分の特徴と多倍長精度積分の適用
濱口 信行
京都大学数理解析研究所RIMS研究集会「科学技術計算における理論と応用の新展開」
数理解析研究所講究録 1791:131-140, 201215. SR11000での多倍長演算使用例
濱口 信行
九州大学情報基盤研究開発センター
全国共同利用システム広報. 2, (3), pp. 102-108,2009.
16. SR11000モデルK1における多倍長ルーチンによる数値積分
の高速化
濱口信行
情報処理学会情報処理学会研究報告ハイパフォーマンス
コンピューティング（HPC）,2008(125(2008-HPC-118)),
109-114 (2008-12-09)
17. 二重指数関数型積分法の素粒子物理学への応用
湯浅 富久子 , 濱口 信行
情報処理学会情報処理学会研究報告ハイパフォーマンス
コンピューティング（HPC）,2008(19(2008-HPC-114)),31-36 (2008-03-05)

18. Baileyアルゴリズムによる多倍長演算の性能評価
石川 正 , 濱口 信行
情報処理学会情報処理学会研究報告ハイパフォーマンス
コンピューティング（HPC）,2008(19(2008-HPC-114)),25-30 (2008-03-05)
19. 理論物理学計算への多倍長ライブラリの適用
濱口信行（日立）
情報処理学会情報処理学会研究報告ハイパフォーマンス
コンピューティング(HPC）,2007-HPC-113 (2007-12-07)