3.1 JAXA LBM コードの性能測定および改善
宇宙航空研究開発機構 石田 崇
富士通株式会社 三吉 郁夫
3.1.1 はじめに航空宇宙分野，特に航空分野の数値流体力学（Computational Fluid Dynamics : CFD）
においては，巡航状態の航空機翼面上で発生する衝撃波を精度良く捉えるために，
Navier-Stokes 方程式をベースとした計算アルゴリズムが発展してきた．非構造格子法，近似リーマン解法，陰的時間積分，乱流モデルといった計算アルゴリズムの進展に加えて，大型計算機の性能向上によって，航空機の巡航状態における空力性能を精度良くかつ実用的な解析時間で予測できるようになってきたため，CFD は風洞試験と並んで航空機設計開発には欠かせないツールとしての地位を確立している．近年では CFD を航空機のオフデザイ
ン（バフェットなどの非定常状態）における空力性能予測に適用するニーズが増えてきている．しかしながら非定常状態に対する Navier-Stokes 方程式をベースとする計算アプロ
ーチでは，計算精度・計算コスト・物理モデル依存の観点でまだ課題が多く，効率よく解析することが難しい．このような状況を打破するアプローチとして，近年格子ボルツマン法（Lattice Boltzmann Method : LBM）が注目されている．LBM では支配方程式がボルツマ
ン輸送方程式であり，以下に示す特徴を持つ．
・支配方程式に非線形項が無い，・支配方程式は非圧縮 Navier-Stokes 方程式に漸近するが，音の伝搬を直接計算できる，・衝突（Collision）と並進（Stream）の 2 ステップで時間発展する，・時間積分は CFL 数が１に固定された陽解法である，・Navier-Stokes 方程式に比べて計算コストが小さい，・格子点毎に依存関係が無いので容易に並列化が出来る．
JAXA では航空機の離着陸形態における空力騒音問題や低速バフェットの解析に向けた大規模非定常流体解析コードとして LBM コードを開発・整備しており，本稿ではメニーコア
WG 内で行った FX100 を用いたコードの性能評価および改善について報告する．

3.1.2 JAXA LBM コードの概要
JAXA LBM コードでは，直交格子法の一つである Building-Cube Method (BCM)をフレームワークとして採用し，LBM と組み合わせて大規模非定常解析コードとして整備している．以降では，LBM および BCM の概要について説明する．
3.1.2.1

LBM のアルゴリズム概要

LBM における支配方程式は，
𝐷𝜕
𝑓(𝐱, 𝐞, 𝑡) = [ + 𝐞 ∙ 𝜵] 𝑓 = 𝛀
𝐷𝑡
𝜕𝑡

( 1 )

であり，格子上で離散化すると以下の離散式が得られる．
𝑓𝑖 (𝐱 + 𝐞𝑖 × 𝑑𝑡, 𝑡 + 𝑑𝑡) = 𝑓𝑖 (𝐱, 𝑡) + 𝑑𝑡 × 𝛀𝑖 ,𝑖 = 1, ⋯ , 𝑏( 2 )

ここで𝑓は状態分布関数，𝐱は位置ベクトル，𝐞は粒子の速度，𝛀は衝突項，𝑏は粒子の自由度の数であり，𝛀には以下に示す BGK モデルを用いた Single Relaxation Time (SRT)を用いるのが一般的である．

1
𝑒𝑞
𝛀𝐵𝐺𝐾
= − (𝑓𝑖 (𝐱, 𝑡) − 𝑓𝑖 (𝐱, 𝑡))
𝑖
𝜏

( 3 )

ここで𝑓 eq は平衡状態分布関数，𝜏は緩和係数である．時間積分は衝突（Collision）と並進（Stream）の 2 ステップからなり，衝突は( 4 )，並進は( 5 )のようにそれぞれ記述できる．
𝑓𝑖∗ (𝐱, 𝑡) = 𝑓𝑖 (𝐱, 𝑡) + 𝛀𝑖 ,𝑖 = 1, ⋯ , 𝑏𝑓𝑖 (𝐱 + 𝐞𝑖 × 𝑑𝑡, 𝑡 + 𝑑𝑡) = 𝑓𝑖∗ (𝐱, 𝑡),𝑖 = 1, ⋯ , 𝑏

( 4 )
( 5 )

SRT の特徴は分散誤差が少ない点であるが，航空機の空力性能予測で対象とするような高Re 数（低粘性）の流れでは計算が非常に不安定ですぐに破綻してしまう．SRT では全ての状態分布関数を同一の緩和係数で緩和するが，この過程に計算を不安定化させる高次のモードをダンプする仕組みが無いためである．この問題を解決するために様々な衝突項モデ
ルが世界的に研究されており，JAXA LBM コードでは Dr. Mertin Geier によって提案された
Cascaded LBM を衝突項に実装している．Cascaded LBM は状態分布関数を( 6 )に示す central
moment に変換し，central moment space で異なる緩和係数で緩和（特に計算を不安定化させる高次のモーメントを dissipative な緩和係数で緩和）させることにより高 Re 数流れでも安定に計算することが出来る．

̃ 𝑝,𝑞,𝑟 = ∑(𝑒𝑖𝑥 − 𝑢𝑥 )𝑝 ∙ (𝑒𝑖𝑦 − 𝑢𝑦 )𝑞 ∙ (𝑒𝑖𝑧 − 𝑢𝑧 )𝑟 ∙ 𝑓𝑖
𝜌𝐌
𝑖

( 6 )

=𝐂∙𝐟
状態分布関数から central moment への変換行列𝐂は，𝐂 = 𝐂(𝐮)とローカルな速度の関数になっているため，格子点毎に計算しなくてはならない．また，central moment space で緩和した後には状態分布関数に逆変換しなければならず，これらの行列演算にかかる計算負荷は大きい．参考に 2 次元 9 速度（D2Q9）モデルにおける変換行列𝐂およびその逆行列である𝐂 −1 を( 7 )および( 8 )にそれぞれ示す．

3.1.2.2

Building-Cube Method (BCM)

Building-Cube Method では，計算領域を立方体領域（Cube）に分割し，各 Cube 内に等間隔直交格子（Cell）を生成して解析するフレームワークである．流体解析のプログラム部
分は single block 用に作成し，Cube ループで全体の計算を回す間に適宜 Cube 間の情報交換を実施する．

図 6 Cube（左図）と Cell（右図）

( 7 )( 8 )3.1.3 コード改善
本章では，FX100（JSS2 MA システム）を用いたコード改善による性能評価について報告する．性能評価には C++で記述された 2 次元 9 速度（D2Q9）モデルのプログラムを用い，Lid
Driven Cavity を計算対象として格子点数128 × 128の single block の直交格子を用いた．
時間積分は 1000 ステップ分実施して性能を評価した．なお，計算は 1core で実施し，コン
パイルには以下オプションを使用した．
コンパイルオプション：-Kfast -Nlst=t -std=c++11以降で出てくるサンプルコード内の変数を以下に定義する．
Q：自由度の数（2 次元：Q=9，3 次元：Q=27）

ftmp：衝突後の状態分布関数

rho

：密度

CM

̃
：central moment 𝐌

omg：モード毎の緩和係数

M：変換行列𝐂

Minv：逆変換行列𝐂 −1

3.1.3.1

Collision の SIMD 化

central moment を central moment space でモード毎に緩和させた後に，元の状態分布関数に逆変換する操作において，衝突後の分布関数を格納する配列：ftmp の定義にダブルポ
インタを用いていたため，コンパイル最適化が十分に効かなかった．そこで SIMD 化を促進するため，以下の変更を実施した．その結果，表 に示す通り 1.74 倍程度高速化すること
が出来た．
変更前
for (int n = 0; n < Q; ++n) {ftmp[I][n] = 0;
for (int m = 0; m < Q; ++m) {ftmp[I][n] += rho[I] * (Minv[m][n] * (CM[m] * (1 - omg[m]) + CMeq[m] * omg[m]));}}
変更後
T ftmp_[Q];for (int n = 0; n < Q; ++n) {ftmp_[n] = 0;}
for (int m = 0; m < Q; ++m) {for (int n = 0; n < Q; ++n) {ftmp_[n] += rho[I] * (Minv[m][n] * (CM[m] * (1 - omg[m]) + CMeq[m] * omg[m]));}}for (int n = 0; n < Q; ++n) {ftmp[I][n] = ftmp_[n];
}
表 1
変更前

変更後

Elapsed(s)29.0696

16.7274

MFLOPS

588.2525

895.9001

MFLOPS/PEAK(%)

1.66552.5452

3.1.3.2

̃ = 𝐂 ∙ 𝐟を求める行列演算の書き下し
central moment 𝐌

central moment を求める行列は速度の関数になっており，格子点毎に実施する変換行列の定義およびそれらの演算はコストが大きい．そこで行列演算を書き下して性能が変化するか調査した．その結果，表 に示す通り更に 1.13 倍程度高速化することが出来，2 次元 9速度モデルでは行列演算の書き下しにある程度効果があることが分かった．
変更前//! Transform matrix from f to Central Moment
T M[Q][Q] = {1, 1, 1, 1, 1, 1, 1, 1, 1,-u, 1 - u,-u,-1 - u, -u, 1 - u, -1 - u, -1 - u, 1 - u,-v, -v,1 - v, -v, -1 - v,1 - v, 1 - v, -1 - v, -1 - v,
u*v, -v * (1 - u),-u * (1 - v), -v * (-1 - u),-u * (-1 - v),(1 - u)*(1 - v),(-1 - u)*(1 - v),(-1 - u)*(-1 - v),(1 - u)*(-1 - v),
u2 - v2, ua - v2,u2 - va,ub - v2,u2 - vb,ua - va, ub - va,ub - vb, ua - vb,
u2 + v2, ua + v2,u2 + va,ub + v2,u2 + vb,ua + va, ub + va,ub + vb, ua + vb,-u * v2, v2*(1 - u),-u * va,v2*(-1 - u),-u * vb,(1 - u)*va, (-1 - u)*va,(-1 - u)*vb, (1 - u)*vb,-u2 * v, -v * ua,u2*(1 - v),-v * ub,u2*(-1 - v),ua*(1 - v), ub*(1 - v),ub*(-1 - v), ua*(-1 - v),
u2*v2, v2*ua,u2*va,v2*ub,u2*vb,ua*va, ub*va,ub*vb, ua*vb};//! Calculation of Central Momentfor (int l = 0; l < Q; ++l) {CM[l] = 0;
for (int n = 0; n < Q; ++n) {CM[l] += M[l][n] * f[I][n];}CM[l] = CM[l] / ro;}

変更後
CM[0] = f_[0] + f_[1] + f_[2] + f_[3] + f_[4] + f_[5] + f_[6] + f_[7] + f_[8];CM[1] = (-u)*f_[0] + (1 - u)*f_[1] + (-u)*f_[2] + (-1 - u)*f_[3] + (-u) * f_[4] +(1 - u) * f_[5] + (-1 - u) * f_[6] + (-1 - u) * f_[7] + (1 - u) * f_[8];CM[2] = (-v)*f_[0] + (-v)*f_[1] + (1 - v)*f_[2] + (-v)*f_[3] + (-1 - v) * f_[4] +(1 - v) * f_[5] + (1 - v) * f_[6] + (-1 - v) * f_[7] + (-1 - v) * f_[8];CM[3] = (u * v)*f_[0] + (-v * (1 - u))*f_[1] + (-u * (1 - v))*f_[2] + (-v * (-1 - u))*f_[3] + (-u * (-1 - v))*f_[4] +((1 - u)*(1 - v))*f_[5] + ((-1 - u)*(1 - v))*f_[6] + ((-1 - u)*(-1 - v))*f_[7] + ((1 - u)*(-1 - v))*f_[8];CM[4] = (u2 - v2)*f_[0] + (ua - v2)*f_[1] + (u2 - va)*f_[2] + (ub - v2)*f_[3] + (u2 - vb)*f_[4] +(ua - va)*f_[5] + (ub - va)*f_[6] + (ub - vb)*f_[7] + (ua - vb)*f_[8];CM[5] = (u2 + v2)*f_[0] + (ua + v2)*f_[1] + (u2 + va)*f_[2] + (ub + v2)*f_[3] + (u2 + vb)*f_[4] +(ua + va)*f_[5] + (ub + va)*f_[6] + (ub + vb)*f_[7] + (ua + vb)*f_[8];CM[6] = (-u * v2)*f_[0] + (v2*(1 - u))*f_[1] + (-u * va)*f_[2] + (v2*(-1 - u))*f_[3] + (-u * vb)*f_[4] +((1 - u)*va)*f_[5] + ((-1 - u)*va)*f_[6] + ((-1 - u)*vb)*f_[7] + ((1 - u)*vb)*f_[8];CM[7] = (-u2 * v)*f_[0] + (-v * ua)*f_[1] + (u2*(1 - v))*f_[2] + (-v * ub)*f_[3] + (u2*(-1 - v))*f_[4] +(ua*(1 - v))*f_[5] + (ub*(1 - v))*f_[6] + (ub*(-1 - v))*f_[7] + (ua*(-1 - v))*f_[8];CM[8] = (u2 * v2)*f_[0] + (v2*ua)*f_[1] + (u2*va)*f_[2] + (v2*ub)*f_[3] + (u2*vb)*f_[4] +(ua*va)*f_[5] + (ub*va)*f_[6] + (ub*vb)*f_[7] + (ua*vb)*f_[8];

表 2
変更前

変更後

Elapsed(s)16.7274

14.7741

MFLOPS

895.9001

849.2416

MFLOPS/PEAK(%)

2.54522.4126

3.1.3.3

̃ を求める行列演算の書き下し
状態分布関数 𝐟 = 𝐂 −𝟏 ∙ 𝐌

2.6.3.2 で行列演算の書き下しによる効果が得られたため，central moment から状態分布関数を求める逆変換の行列演算も書き下して性能が変化するか調査した．その結果，表に示す通り 1.64 倍程度高速化することが出来，2 次元 9 速度モデルにおける逆行列演算に関する書き下しには大きな効果があることが分かった．
変更前
T ftmp_[Q];for (int n = 0; n < Q; ++n) {ftmp_[n] = 0;}
for (int m = 0; m < Q; ++m) {for (int n = 0; n < Q; ++n) {ftmp_[n] += rho[I] * (MiC2[m][n] * (CM[m] * (1 - omg[m]) + CMeq[m] * omg[m]));}}for (int n = 0; n < Q; ++n) {ftmp[I][n] = ftmp_[n];}変更後
T CMpost[Q];for (int n = 0; n < Q; ++n) {CMpost[n] = ro * (CM[n] * (1 - omg[n]) + CMeq[n] * omg[n]);}ftmp[I][0] = (u2 - 1) * (v2 - 1) * CMpost[0] + 0 * CMpost[1] + 0 * CMpost[2] + 4 * u * v * CMpost[3] + 0.5 * (v2 - u2)* CMpost[4] + 0.5 * (u2 + v2 - 2) * CMpost[5] + 2.*u * CMpost[6] + 2.*v * CMpost[7] + 1 * CMpost[8];ftmp[I][1] = -0.5 * u * (u + 1) * (v2 - 1) * CMpost[0] + 0.5 * CMpost[1] + 0 * CMpost[2] + -(2 * u + 1) * v *CMpost[3] + 0.25 * (u2 + u - v2 + 1) * CMpost[4] + 0.25 * (-u2 - u - v2 + 1) * CMpost[5] + -(u + 0.5) * CMpost[6] +-v * CMpost[7] + -0.5 * CMpost[8];ftmp[I][2] = -0.5 * v * (v + 1) * (u2 - 1) * CMpost[0] + 0 * CMpost[1] + 0.5 * CMpost[2] + -u * (2 * v + 1) *CMpost[3] + 0.25 * (u2 - v2 - v - 1) * CMpost[4] + 0.25 * (-u2 - v2 - v + 1) * CMpost[5] + -u * CMpost[6] + -(v + 0.5)
* CMpost[7] + -0.5 * CMpost[8];ftmp[I][3] = -0.5 * u * (u - 1) * (v2 - 1) * CMpost[0] + -0.5 * CMpost[1] + 0 * CMpost[2] + v * (1 - 2 * u) * CMpost[3]+ 0.25 * (u2 - u - v2 + 1) * CMpost[4] + 0.25 * (-u2 + u - v2 + 1) * CMpost[5] + (0.5 - u) * CMpost[6] + -v *CMpost[7] + -0.5 * CMpost[8];ftmp[I][4] = -0.5 * v * (v - 1) * (u2 - 1) * CMpost[0] + 0 * CMpost[1] + -0.5 * CMpost[2] + u * (1 - 2 * v) * CMpost[3]+ 0.25 * (u2 - v2 + v - 1) * CMpost[4] + 0.25 * (-u2 - v2 + v + 1) * CMpost[5] + -u * CMpost[6] + (0.5 - v) *CMpost[7] + -0.5 * CMpost[8];ftmp[I][5] = 0.25 * u * v * (u + 1) * (v + 1) * CMpost[0] + 0 * CMpost[1] + 0 * CMpost[2] + 0.25 * (2 * u + 1) * (2 * v+ 1) * CMpost[3] + 0.125 * (-u2 - u + v2 + v) * CMpost[4] + 0.125 * (u2 + u + v2 + v) * CMpost[5] + 0.25 * (2 * u +1)* CMpost[6] + 0.25 * (2 * v + 1) * CMpost[7] + 0.25 * CMpost[8];ftmp[I][6] = 0.25 * u * v * (u - 1) * (v + 1) * CMpost[0] + 0 * CMpost[1] + 0 * CMpost[2] + 0.25 * (2 * u - 1) * (2 * v +1) * CMpost[3] + 0.125 * (-u2 + u + v2 + v) * CMpost[4] + 0.125 * (u2 - u + v2 + v) * CMpost[5] + 0.25 * (2 * u - 1) *CMpost[6] + 0.25 * (2 * v + 1) * CMpost[7] + 0.25 * CMpost[8];ftmp[I][7] = 0.25 * u * v * (u - 1) * (v - 1) * CMpost[0] + 0 * CMpost[1] + 0 * CMpost[2] + 0.25 * (2 * u - 1) * (2 * v 1) * CMpost[3] + 0.125 * (-u2 + u + (v - 1) * v) * CMpost[4] + 0.125 * (u2 - u + (v - 1) * v) * CMpost[5] + 0.25 * (2 *u - 1) * CMpost[6] + 0.25 * (2 * v - 1) * CMpost[7] + 0.25 * CMpost[8];ftmp[I][8] = 0.25 * u * v * (u + 1) * (v - 1) * CMpost[0] + 0 * CMpost[1] + 0 * CMpost[2] + 0.25 * (2 * u + 1) * (2 * v- 1) * CMpost[3] + 0.125 * (-u2 - u + (v - 1) * v) * CMpost[4] + 0.125 * (u2 + u + (v - 1) * v) * CMpost[5] + 0.25 * (2 *u + 1) * CMpost[6] + 0.25 * (2 * v - 1) * CMpost[7] + 0.25 * CMpost[8];

表 3
変更前

変更後

Elapsed(s)14.77419.0030MFLOPS

849.2416

1457.6530

MFLOPS/PEAK(%)

2.41264.1411

3.1.3.4

Stream の方針変更

Stream では【自分⇒周囲】の順番で状態分布関数をコピーする操作をしていたが，キャッシュの効率を考慮して【周囲⇒自分】の順番に修正して性能が変化するか調査した．その結果，表 に示す通り更に 1.3 倍程度高速化することが出来た．
変更前
for (int I = 0; I < size; ++I) {ijk[0] = I % NDIM[0];ijk[1] = I / NDIM[0];for (int l = 0; l < Q; ++l) {indx[0] = ijk[0] - c[l][0];indx[1] = ijk[1] - c[l][1];indx[0] = (indx[0] >= 0) ? (indx[0] = indx[0]) : (indx[0] = 0);indx[1] = (indx[1] >= 0) ? (indx[1] = indx[1]) : (indx[1] = 0);indx[0] = (indx[0] < NDIM[0]) ? (indx[0] = indx[0]) : (indx[0] = NDIM[0] - 1);
indx[1] = (indx[1] < NDIM[1]) ? (indx[1] = indx[1]) : (indx[1] = NDIM[1] - 1);
J = indx[0] + NDIM[0] * indx[1];flag = !mask[I].isInner() * !mask[J].isInner();f[J][l] = ftmp[I][l] * (flag) + f[J][l] * (1 - flag);}

変更後
for (int I = 0; I < size; ++I) {ijk[0] = I % NDIM[0];ijk[1] = I / NDIM[0];for (int l = 0; l < Q; ++l) {indx[0] = ijk[0] - c[l][0];indx[1] = ijk[1] - c[l][1];indx[0] = (indx[0] >= 0) ? (indx[0] = indx[0]) : (indx[0] = 0);indx[1] = (indx[1] >= 0) ? (indx[1] = indx[1]) : (indx[1] = 0);indx[0] = (indx[0] < NDIM[0]) ? (indx[0] = indx[0]) : (indx[0] = NDIM[0] - 1);
indx[1] = (indx[1] < NDIM[1]) ? (indx[1] = indx[1]) : (indx[1] = NDIM[1] - 1);
J = indx[0] + NDIM[0] * indx[1];f_[l] = ftmp[J][l];}f[I][0] = f_[0];f[I][1] = f_[1];f[I][2] = f_[2];f[I][3] = f_[3];f[I][4] = f_[4];f[I][5] = f_[5];f[I][6] = f_[6];f[I][7] = f_[7];f[I][8] = f_[8];
}

表 4
変更前

変更後

Elapsed(s)9.00306.9164MFLOPS

1457.6530

1897.4156

MFLOPS/PEAK(%)

4.14115.3904

3.1.3.5

コード改善後の性能評価

3.1.3.2~3.1.3.4 の改善によって，検証ケース（格子点数128 × 128）では最終的におよ
そ 4.2 倍高速化することが出来た．そこで格子点数を変化させて SRT とコード改善した
Cascaded LBM の解析時間の比較を行い，その結果を図 7 に示す．コード改善によりコストの高い Cascaded LBM が SRT と比較して，元々5.5 倍であった解析コストが高々1.25 倍程度にまで削減出来ていることが分かった．また，解析コスト増分は格子点数によらずほぼ一定であった．コード改善無しの 3 次元コードの場合，Cascaded LBM の解析コストは SRT の約 13 倍程度となっており（図 8），今後は 2 次元におけるコード改善の結果を反映して 3次元における性能調査を実施する予定である．

図 7 SRT と Cascaded LBM の解析時間の比較（2 次元）

図 8 SRT と Cascaded LBM の解析時間の比較（3 次元）

3.1.4 スレッド並列性能
FX100 1node を用いて行ったスレッド並列性能調査の結果について以下に報告する．性能評価には C++で記述された 3 次元 27 速度（D3Q27）モデルのプログラムを用い，衝突項にはSRT を採用した．計算格子には64𝑏𝑙𝑜𝑐𝑘 × 323 𝑐𝑒𝑙𝑙のマルチブロック直交格子を用いて，スレッド数を 1~32 に変化させて実行時間の変化を調べ，その結果を図 9 に示す．またその際のプロファイラの結果を図 10 に示すが，並列計算時のホットスポットはブロック間で情報通信を実施する interfaceQ という関数であった．LBM は Navier-Stokes 方程式に比べて計算コストが小さいが使用変数が多く，コード改善によって演算に比べて通信がより際立つ結果となっている．今後は通信に関係する時間積分の箇所を分離する等，通信コストをうまく隠蔽できる実装に変更して性能調査を実施する予定である．

図 9 FX100 1node を用いたスレッド並列性能

図 10 プロファイラによるホットスポットの抽出

3.1.5 まとめ
メニーコア WG での JAXA LBM コードの性能調査の結果，以下の知見が得られた．


central moment に関わる行列演算は，2 次元の場合は SIMD 化を促進するより書き下した方が性能向上が期待できる．



SRT と比較して，コード改善した Cascaded LBM の解析コストは 2 次元で約 1.25 倍程度にまで削減することが出来た．3 次元の場合は27 × 27の行列を書き下すことになり，演算量や演算の複雑さが増すため，調査する必要がある．



LBM は Navier-Stokes 方程式ベースのアプローチより格子点当りの演算量は少ないが，格子点当りの変数が 2 次元では 9 変数（NS では 4 変数＋α），3 次元では 27 変数（NSでは 5 変数＋α）であり，ノード内・ノード間の情報通信量が相対的に増えるため，ブロック間の情報通信を隠ぺいする等の処置が必要である．



現状では変数配列をダブルポインタで管理しているが，アーキテクチャやコンパイラの観点でどのようなデータの持たせ方が良いか（SoA or AoS），検討する必要がある．

3.2 MUTSU コード高速化の検討
核融合科学研究所 三浦英昭
3.2.1 はじめに
MUTSU コードは、核融合・トーラスプラズマの不安定性研究や、これに関わる基礎研究のために使用する汎用流体コードである。MUTSU コードでは、基礎研究用の矩形形状（周期境界、非周期境界条件を含む；直角座標系を使用）や、核融合研究向けのトーラス構造を近似的に表現するための歪んだ矩形形状(この場合には非直交曲線座標系で数式を表現する)、
さらには直角座標や非直交曲線座標領域を組み合わせた座標において、拡張 MHD 方程式やNavier-Stokes 方程式を解くことを想定している。対象となる方程式によって使用するモジ
ュールを変える（コード名称も、方程式に応じて MUTSU/cXMHD3D, MUTSH/cNS3D などと名前を変える）が、空間微分や時間発展パートは共通となっている。
MUTSU コードは、入力パラメータで二次精度中心差分、三次精度風上差分(KK スキーム)、
四次精度中心差分、六次精度コンパクト差分、八次精度コンパクト、十次精度コンパクト
差分を切り替えるように設計されている。数値フィルターも、陽的なフィルターから十次精度コンパクトフィルターまでをサポートしている。また、ほぼ同じ構造のモジュールを使って、3 方向に周期的な場合のフーリエ・擬スペクトル法によるシミュレーションが可能になっている。
（こちらは、コード名末尾に”-T3”をつける。
）また、矩形格子形状で複雑物体周りの流れを扱うために、volume penalization 用のモジュールも用意されている。
このコードの性能は、方程式の選択や微分・フィルターに関わるスイッチの選択で性能が大きく変わるため、ここではコンパクト差分法のパートと、MUTSU-T3 で使う 3 次元 FFTのパートについて報告する。
3.2.2 MUTSU, MUTSU-T3 コードの開発・コンパイル
このコードの開発は、主に核融合科学研究所の”プラズマシミュレータ”(富士通株式会
社製 FX100)で行っている。また、一部の開発は、東京大学の Oakforest-PACS(富士通株式会社製)で行っている。FX100 における主要なコンパイルオプションは、以下の通りである。
FC= mpifrtpx

FFLAGS = -Cpp -Cfpp -Kfast -Kreduction -Kopenmp -Kocl -Kparallel -Koptmsg=2 -Qt
また、Oakforest-PACS での主要なコンパイルオプションは以下のとおりである。
FC= mpiifort

FFLAGS = –fpp –O2 -axMIC-AVX512 –qopt-report –qopenmp-mcmodel=medium –convert big_endian3.2.3 コンパクト差分パートの高速化
コンパクト差分では、
以下のような 3 重もしくは 5 重対角連立方程式を解く必要がある。

MUTSU コードにおいてこれらの微分評価を行うサブルーチンは、上記の方程式 3 セットの右辺の計算と、連立方程式の求解パートで構成される。右辺項の計算は、中心差分法などと同じ構造なので、比較的容易に、そこそこ高速な性能（FX100 でピーク性能比 6-8%前後）を達成可能である。連立方程式の解法としては、ここでは LU 分解法を採用した。（メニー
コア WG の会合では、問題規模から考えて、行列反転型の直接解法の方が速いとの指摘があった。
）
八次精度のコンパクト差分法を使用した場合（上の式で係数β=0）について、LU 分解パートのループ融合などの最適化を行った。具体的には、


ベクトル３成分の勾配（速度勾配テンソル）など、複数のコンパクト差分を一括して計算することによって流量を増加



LU 分解で使用する作業配列の次元を増やし、ループ融合
たとえば、入力データ b1, b2, b3 について
do k = 1, n2
do I = 1, n1dm(1,k,i) = b1(i,k,ijk)
dm(2,k,i) = b2(i,k,ijk)
dm(3,k,i) = b3(i,k,ijk)
end doend doなどの処理について、事前に入力データ b1,b2,b3 を一つの配列ｂにまとめ、
do k = 1, n2
do i = 1, n1
do j = 1, 3dm(j,k,i) = b(i,k,ijk,j)end doend doend doとする



融合配列の次元の入れ替え（最内側⇔再外側）を行う

などの工夫を施した。この結果、ピーク性能比で 4-9%を実現した。なお、配列の第２成分(x, y, z の直角座標では y 成分)に関する微分操作については、ストライドアクセスが発生することによる速度の低下のため、３方向についての微分の中では最も低速な演算(ピーク性比 4%)となった。
3.2.4 3 次元 FFT の通信時間隠蔽による周期境界用コード MUTSU-T3 コードの高速化
3 方向に対して周期境界条件が課される流体シミュレーションでは、中心差分法やコンパクト差分法に代えて、擬スペクトル法を使うことが多い。以下は、非圧縮性 Hall MHD シミ
ュレーション版(MUTSU /iHallMHD3D-T3)の場合についての報告である。
擬スペクトル計算の場合、3 次元 FFT がコードの性能を決定づける。3 次元 FFT では通信・データの転置が計算時間のかなりの部分を占めるため、スーパーコンピュータがその性能を十全に発揮するのは難しい。ここでは、3 次元高速フーリエ変換ライブラリ“FFTE”(http://www.ffte.jp)をベースに、そのパーツを分解し、複数のフーリエ変換を束ねることで計算による通信時間の隠ぺいを図った。なお、この隠蔽によるシミュレーションコードの高速化は JHPCN 平成 30 年度共同研究課題にもなっており、共同研究者・FFTEの開発者である高橋大介博士から FFTE Ver.6.0 および 6.2αの提供を受けて行った。
この研究における通信時間の隠ぺいは、流体コードが往々にして保存形式で記述されることを利用している。たとえば速度勾配テンソルの演算は速度場３成分それぞれの勾配を求めることになり、３つのフーリエ変換を３セット行うと考えることができる。１回のフーリエ変換は、一方向へのフーリエ変換、転置・通信を３回繰り返すことで構成されるため、一つの変数のフーリエ変換と他の変数の通信・転置を重ねることで、通信時間をある程度隠蔽できると考えられる。
3.2.5 3 次元 FFT における計算時間の隠ぺい結果
計算時間の隠ぺいを行う場合と行わない場合の MUTSU /iHallMHD-T3 における 1 ステップあたりの計算時間を、FX100 および Oakforest-PACS で測定した。この隠蔽の際の関数の呼
び出しは、もう一つの著名な 3 次元 FFT パッケージである P3DFFT (https://www.p3dfft.net)
において、p3dfft_ftran_r2c_many, p3dfft_ftran_c2r_many といった名前の末尾に”_MANY”がつく関数と対応がつくため、P3DFFT と FFTE の比較も行った。
この計測結果は下表（表 1, 表 2）のとおりである。表 1 は FX100、表 2 は Oakforest-PACSで計測した結果である。表中で「(通常版)」は、3 次元 FFT を個別にコールする場合を示している。P3DFFT の項目で「(_MANY)」は、p3dfft_ftran_r2c_many, p3dfft_ftran_c2r_many
を用いて複数の変数に対する一括フーリエ変換を行う場合を示す。FFTE の「（隠蔽版）」項
は、FFTE のコードを使用して複数のフーリエ変換（ここでは３変数）を行いつつ、通信について隠蔽を行った場合を示している。FX100 上での計測においては、アシスタントコアを使用していない。この点において、FX100 に特化した最適化によってさらに改善をもたらす余地は残っている。
この表からわかる通り、多くの場合において P3DFFT よりも FFTE の通常版の方が短時間で演算を終えていることがわかる。特に表 2 では、P3DFFT の_MANY 版の方が通常版よりも

7-30%高速であること、FFTE では通常版よりも隠蔽版が 11%程度高速であるが、最大規模の計算(N3=40963)では通常版の方が速い。最後の点についてはさらに調査が必要であり、また、P3D は比較的スレッド並列化を苦手とする（フラット MPI に近いほど、性能が顕著に高い）ためにスレッド並列数についての精査が必要であるなど、いくつかの課題は残る。しかし、
概 ね 所 期 の 目 標 を 達 成 す る こ と が で き た と 考 え ら れ る 。 他 方 、 FX100 に 比 べ て 、Oakforest-PACS での性能はやや低いものとなっている。この点から、Oakforest-PACS についての最適化をさらに進める必要がある。

表 1. 核融合科学研究所 F100”プラズマシミュレータ”における 3 次元 FFT の通信時
間の隠ぺい(MUTSU / iHallMHD3D-T3 の 1 ステップに要する時間；単位[秒])
N3

2048

3

10243
512

3

2563

Nodes

P3DFFT

FFTE( 通 常FFTE( 隠 蔽 あ(Processes)(_MANY)版)

り)512 (1024)

37.939―24.990128 (256)

14.78414.094

12.58432 (64)

5.000

4.846

3.742

8 (16)

2.276

3.096

1.751

表 2. Oakforest-PACS に お け る 3 次 元 FFT の 通 信 時 間 の 隠 ぺ い (MUTSU /iHallMHD3D-T3 の 1 ステップに要する時間；単位[秒])Nodes

P3DFFTP3DFFT

FFTE

FFTE(Processes)(通常版)(_MANY)

(通常版)

(隠蔽版)

N3=40963

1024(16384)

327.221306.214231.546247.666

N3=20483

512(16384)

42.55338.87930.66827.524

N3=10243

128(4096)

18.57317.84413.40812.867

N3=5123

32(1024)

8.662

8.228

5.398

4.906

N3=25638(256)

4.321

2.700

2.700

2.550

格子点数

コード開発の目的• 流体現象に関わる様々なシミュレーションを、矩形格子と高次差分で実行する• 直角座標、極座標等直交座標、
非直交曲線座標• 周期境界・非周期境界

コード開発の経緯と現在の方針• AMR コードとしてスタート• AMR 部分はpending• AMR構造を活かした柔軟な構造による多目的化 (たとえば C-grid, O-gridの変形など)を優先

N
W E
S

N
W E
S

N
W E
S

コンパクト差分の実装
ベクトル３成分それぞれのgradientの計算(多成分化による流量の増加)

ベクトル３成分それぞれのコンパクト差分右辺計算(中心差分と似た構造なので、この部分の効率は低くない）

3成分のコンパクト公式左辺のx,y,z一括LU分解
X-方向(Case1)

Y-方向(Case2)
Z-方向(Case3)

Case1:b( i ,j,k)↑この成分への入出力(外側 k で並列化可能)

Case2:b(i, j, k)↑この成分への入出力(ストライドアクセスの発生)

Case3:b(i, j, k)↑この成分への入出力(内側 j で並列化可能)コンパクト差分パート最適化• 最適化前の高コスト部の状態
高コストルーチン

SIMD化

SWP化 自動並列化

経過時間(コード全体200ス
テップ 898秒)

Banbks_2d_3vars

〇

△

△

74.91

Banbks_2d_u2

〇

△

△

61.91

Compact_der_xyz

〇

〇

〇

12.80

Compact_der_z1

〇

〇

△

3.54

• OpenMP化が速い場合とFX100自動並列化の方が速い場合が混
在 (ループ長に関する分岐判断のため）
全ループをOpenMP化すると低速化するので, 個別に切り分け• キャッシュミス等の分析と配列融合で最適化

配列融合について（1-1)
LU分解 banbks_*_3vars … 3変数一括変換
配列融合前
Real(dd),dimension(n1,n2,n3)::b1,b2,b3
配列融合後（内側）
Real(dd),dimension(3,n1,n2,n3)::b
配列融合後（外側）Real(dd),dimension(n1,n2,n3,3)::b(中略)配列融合について（1-2）
最内で融合

(中略)

(中略)

配列融合について（1-3）
最外で融合

(中略)
(中略)

配列融合について（1-4)
結果比較(1ノード4プロセス500回コール時間)
Banbks_2D_3varsの実行時間
Case 1Case 2

Case 3

計(秒)

改善率

修正前

3.24

1.44

1.90

6.59

-

最内

3.17
(4.09%)

1.53(9.09%)

1.74(6.52%)

6.43

2.30%

最外

3.63

1.45

1.90

6.99

-6.12%

Banbks_2D_u2 (Banbks_2D_3vars類似;変数長偶数限定)
実行時間
Case 1Case 2

Case 3

計(秒)

改善率

修正前

0.76

0.39

0.66

1.81

-

最内

0.94

0.72

0.97

2.63

-45.01%

最外

0.75

0.39

0.66

1.80

0.59%

Case 1-3の実行回数にも依存するが、最内がやや有利
非保存的な方程式などでは、高速化の保証がない
Case 1 がコンパクト差分最適化（ピーク性能比向上）のための問題点

3次元FFT“FFTE”を用いた計算時間の隠蔽• 平成30年度JHPCN 共同研究課題「電磁流体力学乱流の高精度・高並列LESシミュ
レーションコード開発研究」
共同研究者︓高橋大介先生（筑波大）他• 3次元FFT … FFTE ( http://www.ffte.jp ) ver 6.0• 基本方針
3次元FFT を分解
多重FFT=>転置⇒多重FFT=>転置⇒多重FFT
通信と計算・転置のオーバーラップ* p3dfft の _MANY 機能と同様

FX100性能測定詳細（ループ交換
等最適化前後比較）
１次キャッシュミス率(%)

Elapsed(secs)
変更前

変更後

改善率

変更前

変更後

改善率

FFTINVXL1 42.56

34.00

20.1%

29.97%

26.94%

10.1%

FFTINVYL1 31.14

26.74

14.1%

30.32%

26.93%

11.2%

FFTINVZL1 36.86

29.26

20.6%

30.20%

26.96%

10.7%

擬スペクトル計算における”_MANY”
𝑢1

𝑢2

𝑢1

𝐹

𝐹

𝐹

𝑢෤𝑘 , 𝑘 , 𝑘ଵଶଷ𝑣෤ 𝑘 , 𝑘 , 𝑘
ଵ

ଶ

ଷ

𝑤
෥𝑘 , 𝑘 , 𝑘
ଵ

ଶ

ଷ

𝐹 ିଵ

𝜕
𝑢
ଵ ଶ
ଷ
𝜕𝑥1 1
𝐹 ିଵ 𝜕
𝑖𝑘2 𝑢෤𝑘 , 𝑘 , 𝑘
𝑢
ଵ ଶ
ଷ
𝜕𝑥2 1
𝐹 ିଵ 𝜕
𝑖𝑘3 𝑢෤𝑘 , 𝑘 , 𝑘
𝑢
ଵ ଶ
ଷ
𝜕𝑥3 1
𝐹 ିଵ 𝜕
𝑖𝑘1 𝑣෤ 𝑘 , 𝑘 , 𝑘
𝑢
ଵ ଶ
ଷ
𝜕𝑥1 2
𝐹 ିଵ 𝜕
𝑖𝑘2 𝑣෤ 𝑘 , 𝑘 , 𝑘
𝑢
ଵ ଶ
ଷ
𝜕𝑥2 2
𝐹 ିଵ 𝜕
𝑖𝑘3 𝑣෤ 𝑘 , 𝑘 , 𝑘
𝑢
ଵ ଶ
ଷ
𝜕𝑥3 2
𝐹 ିଵ 𝜕
𝑖𝑘1 𝑤
෥𝑘 , 𝑘 , 𝑘
𝑢
ଵ ଶ
ଷ
𝜕𝑥1 3
𝐹 ିଵ 𝜕
𝑖𝑘2 𝑤
෥𝑘 , 𝑘 , 𝑘
𝑢
ଵ ଶ
ଷ
𝜕𝑥2 3
𝐹 ିଵ 𝜕
𝑖𝑘3 𝑤
෥𝑘 , 𝑘 , 𝑘
𝑢
ଵ ଶ
ଷ
𝜕𝑥3 3
𝑖𝑘1 𝑢෤𝑘 , 𝑘 , 𝑘

𝜕𝑢1
𝜕𝑢1
𝜕𝑢1
𝑢1
+ 𝑢2
+ 𝑢3
𝜕𝑥1
𝜕𝑥2
𝜕𝑥3

𝜕𝑢2
𝜕𝑢2
𝜕𝑢2
+ 𝑢2
+ 𝑢3
𝜕𝑥1
𝜕𝑥2
𝜕𝑥3
𝜕𝑢3
𝜕𝑢3
𝜕𝑢3
𝑢1
+ 𝑢2
+ 𝑢3
𝜕𝑥1
𝜕𝑥2
𝜕𝑥3
3次元FFT:
正変換*3
逆変換*9

𝑢1

•
•
•

一括処理（通信隠蔽のチャンス）
P3dfft に同様の機能⇒複雑で最適化に難FFTEを用いて実装

隠蔽効果の検証 (FX100)
測定データ 格子数︓2048x2048x2048
ステップ数︓10
スレッド数︓16

Time step loopffte(Xペンシル)P3dfft 2.7.5

ffte(Zペンシル)

通常版

many版

通常版

隠蔽版

通常版

隠蔽版

256

2283.016

2138.858

1577.800

1139.746

1092.112821.8535121036.223

986.4241160.650

973.610762.255549.5881024

482.646428.720511.279349.656

342.041246.9632048

198.744199.216201.124135.184132.365

93.8024096

88.05489.913140.027106.548

86.40661.901nproc• 川島康弘氏（富士通(株)）による計測• 注︓現在 FFTE 6.0⇒6.2α移行、隠蔽版は若干速度低下(最適化中)

隠蔽効果の検証(Oakforest-PACS)• ライブラリ
fftw … fftw 3.3.6-pl2
p3dfft … p3dfft3.7.7
ffte … ffte ver 6.0• コンパイル
mpiifort –fpp –O2-axMIC-AVX512 –qopt-report –qopenmp-mcmodel=medium –convert big_endian

隠蔽効果の検証(Oakforest-PACS)
Time for 1-time-step (secs)
P3dfftFFTE (z-pencil)

Node,procs

通常版

MANY版

通常版

MANY版

N3=40963 1024, 16384

327.221306.214231.546247.666

N3=20483 512, 16384

42.55338.87930.66827.524

N3=10243 128, 4096

18.57317.84413.40812.867

N3=5123

32,1024

8.662

8.228

5.398

4.906

N3=2563

8,256

4.321

2.700

2.700

2.550

FFTE N3=40963 で通常版とMANY版の計算時間の逆転については、要調査。

3.3 生体分子粗視化シミュレータ CafeMol の FX100 での性能測定とチューニング
理化学研究所 情報システム部
検崎博生
富士通株式会社 TC ソリューション事業本部
渡邉健太
3.3.1 はじめに
CafeMol はタンパク質や核酸のような生体分子の粗視化シミュレータである。CafeMol で用いる粗視化モデルでは、全原子モデルに比べて粒子数が 100 分の 1 程度まで減り、計算量が大幅に減る。このことにより長時間の時間発展の計算が可能になり、生体分子の大規
模な構造変化を取り扱うことができるのであるが、同時に並列性能を出しにくいことにもつながっている。具体的には、時間発展の 1 step が 1 ms 以下になることも多く、ノード間通信が問題になってくるので複数ノードでの並列計算を行うことは難しく、1 ノードでの性能向上を図ることが大事になってきている。
本 WG では、CafeMol の性能測定とチューニングについて検崎よる 2 回の発表と富士通渡邉による 1 回の発表があったので、ここではそれらをまとめることとする。最初に検崎による 1 回目の発表では、FX100 と skylake の比較と計算量が一番多い部分について配列の構造を Structure of Array(SoA)から Array of Structure(AoS)への変更を行うチューニング
を行った。次に、渡邉の発表では FX100 が Skylake に比べて性能が出ていない部分について性能測定とチューニングが行われた。最後に、検崎による 2 回目の発表では、渡邊により示されたタイマーのオーバーヘッドの確認と追加のチューニングを行った。
3.3.2 プログラム概要
分子動力学シミュレーションでは、分子の動きを時間積分することにより系を時間発展させるわけであるが、分子間に働く力の計算が最も計算時間が掛かる部分となっている。
CafeMol では力の計算にネイバリングリスト方式を用いて計算量を減らしており、ネイバリ
ングリストを分割する形で MPI と OpenMP によるハイブリッド並列化を行っている。
粗視化された粒子間に働く相互作用はさまざまなものがあり、相互作用毎に異なるチューニングが必要になってくる。ここではタンパク質/DNA 複合体であるヌクレオソームの10,000 step の時間発展の性能測定を行う。ヌクレオソームは粗視化粒子数は 2,000 程度の系で、計算量が一番多いのが静電相互作用の計算で粒子数の 2 乗程度の計算量であり、それ以外に粒子数の数倍程度の計算量であるさまざまな相互作用がある。
3.3.3 測定環境
検崎による性能測定には、理化学研究所の HOKUSAI システムを使用した。


FX100：SPARC64TMXIfx(1.975H)、富士通コンパイラ、-Kfast

Intel Xeon：Intel Xeon Gold 6148(3.1GHz)、インテルコンパイラ
富士通渡邉による性能測定には、以下のシステムを用いた。



FX100：SPARC64TMXIfx(1.5H)、富士通コンパイラIntel Xeon：Intel Xeon Gold 6148(3.1GHz)、インテルコンパイラ

3.3.4 検崎による 1 回目の発表のまとめ以前の性能測定で、FX100 で 1 ノード 32 コアを用いた時は、MPI による 2 プロセス並列と OpenMP による 16 スレッド並列の組合せが一番速くなることを確認していた。Skylake の
1 ノード 40 コアでは、4 プロセスと 10 スレッド並列の組合せが一番速くなった。FX100 とSkylake を比較すると、FX100 は一番計算量の大きい静電相互作用は SIMD 化により Skylake
よりも速くなった。一方、他の相互作用は SIMD 化が困難で Skylake に比べて数倍程度遅くなるものが多く、全体の計算時間も FX100 の方が遅くなった。
次に、skylake において静電相互作用のさらなるチューニングを行った。具体的には静電相互作用につかうネイバリングリストの配列の構造を Array of Structure(AoS)からStructure of Array(SoA)に変更した。結果として、静電相互作用の計算時間が 6.1 s から
4.6 s に高速化された。条件によってどちらの構造を使った方がいいかは変わるようだが、配列の構造を変えることがチューニングの 1 手法として有効であることを確認できた。
3.3.5 富士通渡邉による発表のまとめ
上記の検崎による測定結果で、多くの相互作用で FX100 の性能が Skylake に比べて性能が悪かったので、富士通渡邉により性能測定とチューニングが行われた。状況を簡単にするために、ボンド長、ボンド角、Go ポテンシャルの 3 つの相互作用について調べることとした。
まず、MPI_wtime により時間計測を行っていたが、オーバーヘッドが大きいことが示された。オーバーヘッドが少ないタイマーとして、FX100 では gettod を、Skylake では clockxを使うこととした。
チューニングとしてはボンド長について行われ、手動によるループアンスイッチングと”!ocl simd_redundant_vl()”を配列式の直前に挿入することが行われ、10%程度の性能向上が得られた。また、ループ最後のリダクション演算部にインダイレクトアクセスがあり、データの重なりがないようなレイアウトにすると大きな最適化の効果を得られることが推測された。
3.3.6 検崎による 2 回目の発表のまとめ
富士通渡邉による指摘を受け時間計測の方法についての確認とボンド長とボンド角の相互作用計算のチューニングを行った。
タイマーについては、MPI_wtime と gettod に加えて clock_time も試してみたところ、
MPI_wtime のオーバーヘッドが大きいことが再確認され、gettod と clock_time の差は小さかったので、移植性を考え clock_time を採用することとした。
ボンド長のチューニングについては、OpenMP のスレッド並列のチャンクサイズを 1 にすることによって行った。これは、ボンド長のネイバリングリストでは、リダクション部分
がスレッド並列毎に異なる領域を持ち、ネイバリングリストの作成方法から 4 つ以上離れていればリダクション部分が重ならいないようになっていることを利用している。このことから、チャンクサイズを 1 にしてスレッド並列数が 3 以上ならば、リダクション部分に

重なりがないことが保証されるのである。よって、”!ocl norecurrence()”を付けること
により、ボンド長の計算が SIMD 化できた。ただし、ボンド長ではチャンクサイズを 1 にすることによるオーバーヘッドが大きく、SIMD 化による性能向上は小さかったので、結果として高速化にはならなかった。
また、ボンド角についても同様のチューニングが行った。こちらは SIMD 化による性能向上の効果が比較的大きく、オーバーヘッドを込みで 30％程度高速化した。ただし、Skylakeに比べると性能はまだまだ不十分ではあり、さらなる高速化が求められる。
3.3.7 まとめ
CafeMol を FX100 と Skylake 上で性能測定を行い、いくつかのチューニングを試みた。最も計算量の多い静電相互作用については、skylake について AoS から SoA の変更により性能が向上した。それ以外の計算量が少ない相互作用について FX100 についてチューニングを行ったところ若干の性能向上に成功したが、今回のチューニングにより差が縮まったとは
いえ Skylake に比べて十分な性能を出すことが出来なかった。この理由は CPU やコンパイ
ラによる out of order などの命令の充実に差があるためと考えられ、ポスト京では一部改善が見込めると考えらえる。
また、MPI_Wtime については無視できないオーバーヘッドがあり、時間計測を行う際はsystem_clock やシステム依存なルーチンの利用などを行う必要がある場合がある。この指摘に対し、FX100 では経過時間を測定する際、NTP による時間自動調整結果が補正されないように clock_gettime システムコールを用いているためであると富士通より報告があった。なお、ポスト京で MPI_Wtime のオーバーヘッドは FX100 で用意されている gettod 相
当の処理に変更され、MPI_Wtime の性能は改善される見込みであることも併せて富士通から報告されている。

2017/10/20
SS研メニーコア時代のアプリ性能WG
第4回会合

生体分子粗視化シミュレータCafeMol
Skylakeでのチューニング

検崎博生
理研 情報システム部

概要• 粗視化モデルとCafeMolについて• 生体分子シミュレーションと粗視化• CGタンパク質/DNAモデル• 計算方法と並列化• Skylakeでのベンチマーク結果• 1ノードでのMPI並列数とスレッド並列数• Skylakeでのチューニング結果• 静電相互作用のSIMD化• データ構造変換

HOKUSAI System• GreatWave Massively Pallarell Computer (GWMPC)• 1,080 nodes, 32 cores/node, 34,560 cores• 1.092 PFLOPS (1.975 Hz x 16 FP x 32 cores x 1,080 CPUs)• CPU: SPARCTMXIfx (1.975GHz, 32 cores, 1 CPU/node)• Memory: 32GB/node• Memory BW: 480 GB/s/node• Interconnect: Tofu2 6D-Mesh/Torus (12.5GB/s, bidirectional)• Application Computing Server (ACS)• Large memory server: 2 nodes, 1.5 TB/node• GPU server: 30 nodes, 4 GPU/node, Tesla K20X
• BigWaterfall Massively Parallel Computer (BWMPC)• 840 nodes, 40 cores/node, 33,600 cores• 2.58 PFLOPS (2.4 GHz x 32 FP X 20 cores x 1680 CPUs)• CPU: Intel Xeon Gold 6148(2.4 GHz, 20 cores, 2 CPUs/node)• Memory: DDR4-2666 96GB/node• Memory BW: 255 GB/s/node• Interconnect: InfiniBand EDR (12.6 GB/s, bidirectional)CafeMol: Coarse-grained biomolecularsimulation software H.Kenzaki, et al, J. Chem. Theor. Chem. (2011)• CafeMol• Coarse-grained (CG) Protein/Nucleic
acid/lipid models• CafeMol 3.0 source code anddocumentation are released at
http://www.cafemol.org• Takada Lab (Kyoto Univ)CG protein model: Go-like model

C. Clementi, H. Nymeyer, and J.N. Onuchic, J. Mol. Biol. (2000)

Based on the energy landscape theory
: bond angle
: dihedral angle(0 means native state)

Vprotein  Vlocal  Vgo  Vex






 K  1  cos     K  1  cos 3   

Vlocal  K b  ri ,i1  r0 i ,i 1  K   i   0 i
2

i

2

i

1

3



i

0i

i



i

0i

i

  r 12  r 10 
0 ij
0 ij
Vgo  go  5   6   rij r  
i, j 

 ij  nativeVex   ex

nonnative


i, j

12

 
 
r 
 ij 

X. Yao, H. Kenzaki, S. Murakami,
and S. Takada, Nature Comm. (2010)

CG DNA model: 3SPN.1 force field(local, base pair and excluded volume interactions)
E.J. Sambrisiki, D.C. Schwartz, and J.J. de Pablo,Knotts, Biophys. J. (2009)

Vdna  Vlocal  Vstack  Vbp  Vex  Vqq  Vsolv
Vlocal  K  r  r 
 K  r  r
 K       K  1  cos   
2

b1

4

0 i ,i 1

i ,i 1b2i

0 i ,i 1

i ,i 1

i

2



i

i

0i



i

0i

i

  12    6 
0 ij 
0 ij 

Vstack  4 1  
 r  i , j  rij ij  Nst  r 12  r 10 
0 ij
0 ij
Vbp   4bpi 5   6  
r  
  rij i, j ij  Nbp 12   6 Vex  4 1   0    0     1 (if rij  dcut ), r  
i , j  rij  ij  Nex 0 (if rij  dcut ): bond angle
: dihedral angle(0 means B-type DNA)CG DNA model: 3SPN.1 force field(electrostatic and solvation interactions)  rij / D
qiqj

e
Vqq  
i , j  4 0 (T , C ) rij  (T , C )   (T )a (C )N

Debye-Huckel theory

=78

 (T )  249.4  0.788T / K  7.20  10  4 (T / k ) 2
a (C )  1.000  0.2551C / M
 5.15110  2 (C / M ) 2  6.889  10 3 (C / M ) 3
Nsolv



Vsolv    s 1  ei js   N AI  a ( rij  rs ) 2

sa-1 = 5.333A
rs = 13.38A
0  0.504982

eN  e0 (1  [1.40418  0.268231N nt ]1 )AI  0.474876(1  {0.148378  10.9553[ Na  ]}1 )Debye length  0RT  D  2 2 
 2 N Aeq I 計算規模と並列化• 計算規模• 数十から数万粒子数程度• ヌクレオソーム1個の系：1,854粒子• 並列化• 系の時間発展をMPIとOpenMPでハイブリッド並列化• 数十並列程度の並列化• 力の計算はネイバリングリスト方式• レプリカ交換によるMPI並列• 通信量が少ないので大規模並列化• 平衡量を求める計算に使える
ネイバリングリストの例 (2体の相互作用、i<jだけをリストアップ)
i

1

j

2

2
4

7

9

5

3
8

4

4
5

9

7

5
8

6

6
7

9

8

7
9

8

8
9

9

力場の計算はネイバリングリスト方式• 100stepに1回程度計算• ネイバリングリストはノード毎に分割• ローカルな相互作用以外は基本的に2体の相互作用の形になるので、あらかじめ各プロセスで計算する1つめの粒子iのリストを作っておく• 各プロセスは割り当てられた粒子iのリストについて力場を計算する。
• 最も大きな配列(粒子数x100〜10,000)• 相互作用の到達範囲毎にネイバリングリストを作成• ローカルな相互作用• ボンド長、ボンド角、二面角• 計算を始める前にネイバリングリストを作成しておく• 近距離相互作用(<20A)• 郷ポテンシャル、排除体積、塩基対• 中距離相互作用(20-50A)• 疎水相互作用• 多体のため複雑なネイバリングリストに• 遠距離相互作用(>50A)• 静電相互作用• イオン強度により相互作用距離が変化

Simulation of nucleosome
H.Kenzaki and S.Takada, PLoS. Comp. Biol. (2015)

Ion strength=300mM

Skylakeでの性能測定 nucleosome (105step, 100mM)
1 node (40 cores)利用で、MPI並列数xスレッド並列数=40を固定して、MPI並列数を変えて性能測定
コンパイラは、インテルコンパイラ17.0.4
最適化オプションは-fopenmp –O3 –xCORE-AVX512

1nodeでは4MPI並列x10OpenMP並列のときに一番高い性能

FX100との比較 nucleosome (105step, 100mM)
FX100(32 cores)

Skylake(40 cores)-Kopenmp -Kfast -Kparallel-fopenmp -O3 -xCORE-AVX512

2MPIx16OpenMP

4MPIx10OpenMP

force

36.7

18.1

_force(comm)

6.2

5.5

_force(local)

8.8

2.0

_force(go)

1.4

0.2

_force(pnl)

11.8

2.2

_force(ele)

5.8

7.9

random7.0

4.0

neighbor

12.7

1.9

update3.3

2.5

output

0.5

1.0

ope

59.9

21.9

comm

6.8

5.9

main_loop

65.7

27.8

Skylakeの方が2倍以上高速化
静電相互作用以外の相互作用とネイバリングリストで4-6倍

静電相互作用の計算方法
1. ネイバリングリストiele2mp(2, lele)を各プロセスで計算
i < jとなるペアだけを格納

2. simu_force.F90!$omp parallel private(tn)call simu_force_ele(force_mp_l(1, 1, tn))

スレッド毎にforce_mpの異なる領域を用意

3. simu_force_ele.F90!$omp do private(…)
do iele1 = 1, lele
imp1 = iele2mp(1, iele1)
imp2 = iele2mp(2, iele2)
v21(1:3) = xyz_mp(1:3, imp2) – xyz_mp(1:3, imp1)
dist2 = v21(1)**2 + v21(2)**2 + v21(3)**2
if(dist2 > cutoff2) cycle
dist1 = sqrt(dist2)
rdist1 = 1.0/dist1
dvdw_dr = coef(iele)*rdist1*rdist1*(rdist1+rcdist)*exp(-dist1*rcdist)
force_mp(1:3, imp1) = force_mp(1:3, imp1) – dvdw_dr*v21(1:3)
force_mp(1:3, imp2) = force_mp(1:3, imp2) + dvdw_dr*v21(1:3)
end do
SIMD化困難!$omp end do nowait!$omp atomicでは遅くなる

静電相互作用のSIMD化
1. ネイバリングリストiele2charge_k(ncharge, ncharge_mpi)を各プロセスで計算。
電荷をもっている粒子のリストicharge2mp(ncharge)と座標xyz_ele(3, ncharge)を作っておく。
タンパク質は4/20が、DNAは1/3の粒子が電荷をもっているため。

2. simu_force.F90!$omp parallel private(tn)call simu_force_ele2(force_mp_l(1, 1, tn))

3. simu_force_ele2.F90!$omp do private(…)
do icharge1 = 1, ncharge
imp1 = icharge2mp(icharge)
for(1:3) = 0.0
do iele = 1, lele_k(icharge)
電荷を持った粒子だけの座標
jcharge = iele2charge_k(iele, icharge)v21(1:3) = xyz_ele(1:3, icharge) – xyz_ele(1:3, jcharge)…
for(1:3) = for(1:3) + dvdw_dr*v21(1:3)
SIMD化OKend do
ただし計算量はforce_mp(1:3, imp1) = force_mp(1:3, imp1) + for(1:3)
本来必要な量の2倍
end do
計算時間は7.9s→6.1s!$omp end do nowait

静電相互作用のデータ構造変換• Array of Structure(AoS)からStructure of Array(SoA)に
データ構造を変換。
• xyz_ele(1:3, ncharge)をxyz_ele(ncharge, 1:3)に変換。
Skylake(40 cores)

Skylake(40 cores)
SIMD

Skylake(40 cores)
SIMD+DATA

force

18.1

16.1

14.8

_force(comm)

5.5

5.1

5.2

_force(local)

2.0

2.1

2.0

_force(go)

0.2

0.2

0.2

_force(pnl)

2.2

2.2

2.2

_force(ele)

7.9

6.1

4.6

random

4.0

4.0

3.9

neighbor

1.9

1.8

1.9

update

2.5

2.5

2.5

output

1.0

1.1

1.0

ope

21.9

20.6

19.0

comm

5.9

5.6

5.6

main_loop

27.8

26.2

24.6

計算時間は6.1s→4.6s

まとめ• Skylakeでのベンチマーク結果• 1ノードでは、4MPI並列x10スレッドの時に一番高い性能。
• FX100に対して2倍以上の高速化。
• 静電相互作用は余り変わらない。
• 静電相互作用以外やネイバリングリングリストの計算で4-6倍の差。

• Skylakeでの静電相互作用のチューニング
• SIMD化により1.3倍高速化。
• データ構造を変えることにより1.3倍高速化。

SS研 メニーコア時代のアプリ性能検討WG[A.I. 3]
CafemolのFX100とIntel Skylakeの性能差について
2018年7月20日
富士通株式会社
TCソリューション事業本部
渡邊 健太
0Copyright 2018 FUJITSU LIMITED

目次
 目的
 評価環境 現状分析
 チューニング まとめ

1Copyright 2018 FUJITSU LIMITED目的[A.I. 3] FX100おける cafemol の性能を分析する
・ Intel Skylakeとの性能差の原因を調査する
・ FX100で高速化の検討

剣崎さんから評価対象のサブルーチンを選定いただいています。
[メールの一部を抜粋]> 〇調べていただきたいサブルーチン> simu_force_bond.F90(タイマーの名称: __force(bond))> simu_force_angle.F90(タイマーの名称: __force(angle))> simu_force_nlocal_go.F90(タイマーの名称: __force(go))>> bondとangleは鎖で繋がった局所的な2体、3体の相互作用で、メモリアクセスが連続的> になっています。
> goはランダムアクセスを含むような2体の相互作用となっています。
> どれもskylakeにくらべて6倍程度遅くなっています。

2

Copyright 2018 FUJITSU LIMITED評価環境 FX100 と Intel Skylake
FX100

Intel Skylake

名称

SPARC64 XIfxIntel Xeon Gold 6148

動作周波数

1.50 GHz[Base] 2.40 GHz[MAX] Normal : 3.10 GHz
AVX512 : 2.20 GHz

コア数

32コア + アシスタントコア

20コア

キャッシュメモリ

L1I$ , L1D$ : 64KB/コア
L2$ : 24MB/CPU

L1I$ , L1D$ : 32KB/コア
L2$ : 1MBL3$ : 27.5MB

理論ピーク性能

768GFlops/CPU (倍精度)

1,408GFlops/CPU (AVX512,倍精度)

メモリ帯域(240+240)GB/s/CPU

128GB/s/CPU(DDR-2666, 6chnnels)

CPU数, コア数

1CPU/ノード,(32コア + アシスタントコア)/ノード

２CPU/ノード, 40コア/ノード

理論ピーク性能

768GFlops/CPU (倍精度)

2,816GFlops/CPU (AVX512,倍精度)

メモリ帯域(240+240)GB/s

256GB/s

コンパイラー

富士通コンパイラ

Intel Compiler v18.1

C
P
U

ノ
|
ド

3

Copyright 2018 FUJITSU LIMITED現状分析 (1/5)
 各計算機の性能を確認
 測定条件[コンパイルオプション]
FX100 : -Kopenmp,fast,parallel,optmsg=2 –Qa,m,p,t,x
SKL : -O3 –qopenmp -qopt-report=5

 測定結果
評価環境

測定者
計算機

タイマーの結果

並列数__force(bond)__force(angle)__force(go)

FX100
①
渡邊
FX100
(1.5GHz)2p x 16t
0.63162.06881.8975

Intel Skylake

②
剣崎さん
FX100
(1.975GHz)
2p x 16t
0.462
1.57371.4399#1③
渡邊
Skylake(2.40GHz)
4p x 10t0.11420.29700.2090④
剣崎さん
Skylake(2.40GHz)
4p x 10t
0.09120.28030.1935性能値の比較
②/④
5.1
5.6
7.4

#2

#1 クロック周波数の差と同じ性能差#2 BIOS設定SNC(Sub-NUMA Clustering)が違う?
③ : SNC disable④ : SNC enable
4

Copyright 2018 FUJITSU LIMITED現状分析 (2/5)
 コンパイラの最適化状況_force(bond)コンパイラ最適化_force(angle)_force(go)

FX100

SKL

FX100

SKL

FX100

SKL

SIMD

×

×

×

×

×

×

SWP

×

-

×

-

×

-

FULLUNROLLING*1

○

○

○

○

○

○

UNSWITCHING*2

○

-

○

-

×

-

PREFETCH

×

×

○

×

×

×

*1 : doループにではなく、ループ内部の配列式に対して適応*2 : ループ内で不変な不等式をループ外に出す最適化 実行効率の差
SIMD化されていないため、各サブルーチンの実行効率の差は下表となる
サブルーチン

① FX100(1.975GHz)② SKL (3.1GHz)*No-SIMD_force(bond)

0.462

0.0912_force(angle)

1.57370.2803_force(go)1.43990.1935③ No-SIMD時のピーク性能比
FX100 / SKL

実行効率の差①/②x③

FX100 / SKL = 0.51

2.6

FX100 : 252.8GFlops(1.975GHz x 32c x 4)
SKL : 496.0GFlops(3.100GHz x 40c x 4)

5

2.9
3.8
Copyright 2018 FUJITSU LIMITED現状分析 (3/5)
 PA情報(実行時間内訳)

force_bondforce_angle_force(bond)_force(angle)_force(go)4命令コミット

9.14%

9.68%

9.27%

2/3命令コミット

10.24%

11.42%

13.08%

１命令コミット

12.09%

12.42%

9.98%

浮動小数点演算待ち

26.53%

50.01%

35.12%

整数演算待ち

2.41%

4.77%

11.27%

浮動小数点ロードL1アクセス待ち

5.85%

2.26%

2.96%

浮動小数点ロードL2アクセス待ち

9.21%

2.52%

8.57%

整数ロード
L1アクセス待ち

5.97%

1.90%

2.65%

整数ロード
L2アクセス待ち

8.88%

1.81%

3.58%

force_go 命令スケジュールの改善で高速化できる部分(浮動小数点/整数演算待ち、浮動小数点/整数L1アクセス待ち)_force(bond)_force(angle)_force(go)40.76%

58.94%

52.00%

 各サブルーチンの××命令コミットの割合は30%を超えており、命令数もボト
ルネックになっている。
6

Copyright 2018 FUJITSU LIMITED現状分析 (4/5)
 PA情報(メモリ・キャッシュスループット情報)@FX100
 メモリ・キャッシュスループット情報
サブルーチン

L1ビジー率

メモリビジー率

L2スループット

メモリスループット_force(bond)

21%

7%

0%

32.08GB/s

0.00GB/s_force(angle)

19%

2%

0%

8.13GB/s

0.00GB/s_force(go)

23%

5%

0%

20.77GB/s

0.00GB/s

L2ビジー率 キャッシュミス情報
L2ミスdm率μDTLBミス率

mDTLBミス率

0.00%

70.27%

1.44%

0.00%

81.49%

0.00%

83.91%

0.33%

0.00%

87.43%

0.00%

85.24%

0.34%

0.00%

サブルーチン

L1Dミス率L1Dミスdm率_force(bond)

2.63%

91.05%_force(angle)

0.51%_force(go)

1.15%

L2ミス率

 メモリビジー率・メモリスループットから、メモリアクセスはほとんど無い。
 L1キャッシュミス率は、理論値(3.15%)以下である。

7

Copyright 2018 FUJITSU LIMITED現状分析 (5/5)
 基本プロファイル情報の採取[実行時のコマンド]
fipp -C –d fipp_dir -I hwm,cpu –P nouserfunc[プロファイル情報のテキスト出力]
Procedures profile*************************************************************************************
Application - procedures*************************************************************************************
Cost% Operation (S)Barrier% StartEnd------------------------------------------------------------------------------------19792 100.0000
1873.0138
1695 8.5641--- Application------------------------------------------------------------------------------------2551 12.8890
241.4149
1600 62.7205
406
515 simu_tintegral.get_random_number._OMP_1_
2311 11.6764
218.7002
00.0000--- __lll_lock_wait
2029 10.2516
192.0169
00.00005
143 simu_force_pnl_
2015 10.1809
190.6891
00.0000--- __pthread_mutex_unlock_usercnt
1990 10.0546
188.322200.000045
87 simu_neighbor_list._OMP_1_
1276 6.4470
120.753700.00004
88 simu_force_ele2_1023 5.1688
96.811200.00007424 simu_force_pnl2_
806
4.072476.277100.0000 510
549 mt_stream.mt_genrand_int32_
538
2.718350.913300.0000 1032 1080 gf2xe.mult_i32_496
2.506146.938800.000084
124 simu_neighbor_list_ele2._PRL_3_
443
2.238341.924000.0000--- __g_dsin4162.101939.369000.0000 708
729 mt_stream.mt_genrand_double1_
413
2.086739.083900.000014
113 simu_force_nlocal_go_378
1.909935.772500.0000 288
376 simu_force_fdih.calc_force_fdih_
231
1.167121.858200.00006329 mloop_flexible_local_
186
0.939817.602300.0000 166
236 simu_force_fdih.calc_phi_
1520.768014.384000.0000--- __GI___gettimeofday_internal
146
0.737713.8166
00.00004
191 simu_force_dih_
1390.702313.155600.0000--- __g_dscn2
131
0.661912.396800.000013770 simu_neighbor_assign_

8

Copyright 2018 FUJITSU LIMITEDコンパイルリスト (_force(bond))
55
56

57
58
59
60
61

62
63
64
65
66
67
68
69
70

71
72
73

74
75

76 2
!$omp do private(imp1,imp2,v21,dist,ddist,ddist2,for, &77 2 p 2v
else!$omp&
force,efull,iunit,junit,isys,istat)78 2! calc energy<<< Loop-information Start >>>79 2 p 2v
efull = (coef_bd(1, ibd) + coef_bd(2, ibd) * ddist2) * ddist2
 メインループの最適化情報<<< [OPTIMIZATION]
80 2 p 2v
iunit = imp2unit(imp1)<<< UNSWITCHING
SIMD/SWPはなし
81 2 p 2v
junit = imp2unit(imp2)<<< Loop-information End >>>82 2 p 2m
ene_unit(iunit, junit) = ene_unit(iunit, junit) + efull
1 p 2sdo ibd = ksta, kend
83 2
1
84 2 p 2v
isys = ibd2sysmbr_mgo(1, ibd)
1 p 2v
imp1 = ibd2mp(1, ibd)
 インダイレクトアクセス
85 3 p 2s
if(isys == 0) then If文
1 p 2v
imp2 = ibd2mp(2, ibd)<<< Loop-information Start >>>1<<< [OPTIMIZATION]<<< Loop-information Start >>><<< FULL UNROLLING<<< [OPTIMIZATION]<<< Loop-information End >>><<< FULL UNROLLING
86 3 p 6s
force_mp(1:3, imp1) = force_mp(1:3, imp1) - force(1:3)<<< Loop-information End >>>87 3 p 6s
force_mp(1:3, imp2) = force_mp(1:3, imp2) + force(1:3)
1 p 6v
v21(1:3) = xyz_mp_rep(1:3, imp2,irep) - xyz_mp_rep(1:3, imp1,irep)
88 3 p 2selse
1
89 3 p 2s
istat = ibd2sysmbr_mgo(2, ibd)
1 p 2v
dist = sqrt(v21(1)**2 + v21(2)**2 + v21(3)**2)<<< Loop-information Start >>>1 p 2v
ddist = dist - bd_nat(ibd)<<< [OPTIMIZATION]
1 p 2v
ddist2 = ddist**2<<< FULL UNROLLING
1<<< Loop-information End >>>1! calc force
90 3 p 6s
force_mp_mgo(1:3, imp1, istat, isys) = force_mp_mgo(1:3,
1 p 2v
for = (coef_bd(1, ibd) + 2.0e0_PREC * coef_bd(2, ibd) * ddist2) * &imp1, istat, isys) - force(1:3)
1(-2.0e0_PREC * ddist / dist)
91 3 p 6s
force_mp_mgo(1:3, imp2, istat, isys) = force_mp_mgo(1:3,<<< Loop-information Start >>>imp2, istat, isys) + force(1:3)<<< [OPTIMIZATION]
92 3 p 2
end if<<< FULL UNROLLING
93 2 p 2v
end if<<< Loop-information End >>>94 1 p 2v
end do1 p 6mforce(1:3) = for * v21(1:3)
95!$omp end do nowait
1
2 p 2vif(inmgo%i_multi_mgo == 0) then  If文 ループ中で不変<<< Loop-information Start >>><<< [OPTIMIZATION]<<< FULL UNROLLING<<< Loop-information End >>>2 p 6force_mp(1:3, imp1) = force_mp(1:3, imp1) - force(1:3)
2 p 6
force_mp(1:3, imp2) = force_mp(1:3, imp2) + force(1:3)

 サブルーチン内のメインループはIntel SkylakeでもSIMD化されていない
 FULLUNROLLINGはループ内の配列式に対して適応
9

Copyright 2018 FUJITSU LIMITEDコンパイルリスト (_force(angle))
52
53
54

!$omp do private(imp1,imp2,imp3,v21,v32,c11,c22,c21, &!$omp&
co_theta,dba,t3,for,force_21,force_32, &!$omp&
efull,iunit,junit,isys,istat)<<< Loop-information Start >>><<< [OPTIMIZATION]<<< UNSWITCHING メインループの最適化情報<<< PREFETCH:2
SIMD/SWPはなし<<<
ba_nat: 2<<< Loop-information End >>>s
do iba=ksta,kend

55 1 p56 1
57 1 p m
if (coef_ba(1, iba) < ZERO_JUDGE .and. coef_ba(2, iba) <ZERO_JUDGE) cycle
 cycle文
58 1
59 1 p s
imp1 = iba2mp(1, iba)
60 1 p s
imp2 = iba2mp(2, iba)
 インダイレクトアクセス
61 1 p s
imp3 = iba2mp(3, iba)
62 1<<< Loop-information Start >>><<< [OPTIMIZATION]<<< FULL UNROLLING<<< Loop-information End >>>63 1 p 3s
v21(1:3) = xyz_mp_rep(1:3, imp2, irep) - xyz_mp_rep(1:3, imp1,
irep)
64 1 p 3s
v32(1:3) = xyz_mp_rep(1:3, imp3, irep) - xyz_mp_rep(1:3, imp2,
irep)
65 1
66 1 p s
c11 = v21(1) * v21(1) + v21(2) * v21(2) + v21(3) * v21(3)67 1 p s
c22 = v32(1) * v32(1) + v32(2) * v32(2) + v32(3) * v32(3)68 1 p s
c21 = v32(1) * v21(1) + v32(2) * v21(2) + v32(3) * v21(3)69 1
70 1 p s
co_theta = - c21 / sqrt(c11 * c22)
71 1
72 2 p sif(co_theta > 1.0e0_PREC) then
 If文73 2 p m
co_theta = 1.0e0_PREC
74 2 p selse if(co_theta < -1.0e0_PREC) then75 2 p s
co_theta = -1.0e0_PREC
76 2 p s
end if
77 1
78 1 p s
dba = acos(co_theta) - ba_nat(iba)
79 1
80 1 p s
t3 = c11 * c22 - c21**2
81 2 p s
if(t3 <= 1.0e0_PREC) then If文82 2 p s
t3 = 1.0e0_PREC
83 2 p s
end if
84 1

85
86
87

1
1 p
1

! calc force
for = 2.0e0_PREC * coef_ba(1, iba) * dba / sqrt(t3) - &coef_ba(2, iba) / sqrt(c11 * c22)<<< Loop-information Start >>><<< [OPTIMIZATION]<<< FULL UNROLLING<<< Loop-information End >>>88 1 p 3s
force_21(1:3) = for * (v21(1:3) * (c21 / c11) - v32(1:3))
89 1 p 3s
force_32(1:3) = for * (v32(1:3) * (c21 / c22) - v21(1:3))
90 1
91 2 p s
if(inmgo%i_multi_mgo == 0) then  If文 ループ中で不変<<< Loop-information Start >>><<< [OPTIMIZATION]<<< FULL UNROLLING<<< Loop-information End >>>92 2 p 3
force_mp(1:3, imp1) = force_mp(1:3, imp1) - force_21(1:3)
93 2 p 3
force_mp(1:3, imp2) = force_mp(1:3, imp2) + force_21(1:3) force_32(1:3)
94 2 p 3
force_mp(1:3, imp3) = force_mp(1:3, imp3) + force_32(1:3)
95 2
96 2 p selse:
114 2 p s
end if
115 1
116 1 p v
end do117!$omp end do nowaits サブルーチン内のメインループはIntel SkylakeでもSIMD化されていない
 FULLUNROLLINGはループ内の配列式に対して適応 cycle文あり

10

Copyright 2018 FUJITSU LIMITEDソースコード (force_go)67!$omp do private(imp1,imp2,rcut_off2,v21,dist2,roverdist2,roverdist4, &68!$omp&
roverdist8,roverdist12,roverdist14,dgo_dr,for,imirror)
69 1 pdo icon=ksta,kend
70 1
71 1 pimp1 = icon2mp(1, icon) インダイレクトアクセス
72 1 pimp2 = icon2mp(2, icon)
73 2 pif (iclass_mp(imp1) == CLASS%RNA .AND. iclass_mp(imp2) ==CLASS%RNA) then
 If文74 2 prcut_off2 = rcut_off2_rna
75 2 pelse
76 2 prcut_off2 = rcut_off2_pro
77 2 pendif
78 1
79 2 pif(inperi%i_periodic == 0) then
 If文<<< Loop-information Start >>><<< [OPTIMIZATION]<<< FULL UNROLLING<<< Loop-information End >>>80 2 p 3
v21(1:3) = xyz_mp_rep(1:3, imp2, irep) - xyz_mp_rep(1:3, imp1,
irep)
81 2 pelse<<< Loop-information Start >>><<< [OPTIMIZATION]<<< FULL UNROLLING<<< Loop-information End >>>82 2 p 3
v21(1:3) = pxyz_mp_rep(1:3, imp2, irep) - pxyz_mp_rep(1:3, imp1,
irep)
83 2 pcall util_pbneighbor(v21, imirror)
84 2 pend if
85 1
86 1! v21(1:3) = xyz_mp_rep(1:3, imp2, irep) - xyz_mp_rep(1:3, imp1,
irep)
87 1
88 1 pdist2 = v21(1)*v21(1) + v21(2)*v21(2) + v21(3)*v21(3)89 1
90 1 proverdist2 = go_nat2(icon) / dist2

91 1 pif(roverdist2 < rcut_off2) cycle
 cycle文
92 1
93 1 proverdist4 = roverdist2 * roverdist2
94 1 proverdist8 = roverdist4 * roverdist4
95 1 proverdist12 = roverdist4 * roverdist8
96 1 proverdist14 = roverdist12 * roverdist2
97 1
98 1 pdgo_dr = 60.0e0_PREC * coef_go(icon) / go_nat2(icon) *(roverdist14 - roverdist12)
99 1
100 2 pif(dgo_dr > DE_MAX) then
 If文
101 2
!
write (*, *) "go", imp1, imp2, dgo_dr
102 2 pdgo_dr = DE_MAX
103 2 pend if104 1! if(dgo_dr > 5.0e0_PREC) dgo_dr = 5.0e0_PREC
105 1<<< Loop-information Start >>><<< [OPTIMIZATION]<<< FULL UNROLLING<<< Loop-information End >>>106 1 p 3
for(1:3) = dgo_dr * v21(1:3)
107 1 p 3
force_mp(1:3, imp2) = force_mp(1:3, imp2) + for(1:3)
108 1 p 3
force_mp(1:3, imp1) = force_mp(1:3, imp1) - for(1:3)
109 1
110 1 pend do111!$omp end do nowait

 サブルーチン内のメインループはIntel SkylakeでもSIMD化されていない
 UNSWITCHINGは適応されていない
 FULLUNROLLINGはループ内の配列式に対して適応 cycle文あり
11

Copyright 2018 FUJITSU LIMITED現状分析まとめ 3つのサブルーチンの1ノードのFX100とSKLの性能は、5~7倍SKLが良い
 FX100 (1.5GHz)とFX100 (1.975GHz)の性能差から、CPUの演算性能(周波数)の差が実測性能に直結
 PA情報からメモリアクセスはほとんど無い FX100/SKLでSIMD化されていない。
1ノードのFX100/SKLのピーク演算性能差が2倍程度異なるため、
FX100とSKLの実行効率の差は2.6~3.6倍になる

 PA情報(実行時間内訳)から、
 命令スケジューリングの改善で高速化できる部分は、40~58%ある。
 SWP化の検討
 命令数(XX命令コミット)の削減で高速化できる部分は、約30%ある。
 SIMD化の検討

12

Copyright 2018 FUJITSU LIMITEDプログラム全体に対する性能改善 1
 基本プロファイラから得られた “__lll_lock_wait”に対する対応
 実行時環境変数を指定
export XOS_MMM_L_ARENA_LOCK_TYPE=0

 「0」を指定することで、malloc要求を並列処理することができる。
 「1」がデフォルト値、逐次処理される。

 測定結果 (プログラム全体の経過時間)
指定なし

指定あり

改善率

130.75 (s)115.6313.1%# メインループおよび調査対象のサブルーチンの性能への影響はない。

13

Copyright 2018 FUJITSU LIMITEDプログラム全体に対する性能改善 2
 オーバーヘッドが大きい時間計測ルーチンを置き換える タイマールーチンの置換え
オリジナル

変更後

MPI_wtime

FX100

SKL

gettodclockx タイマー結果
FX100
サブルーチンSKL

FX100 (1.5GHz)

MPI_wtimeの性能差

タイマー置換え後の性能差

FX100 (1.975GHz)

MPI_wtime

gettodgettodMPI_wtimeclockxforce_bond

0.63160.45090.34240.11420.1005

5.1

3.4

force_angle

2.06881.89571.43980.29700.2809

5.6

5.1

force_go

1.89751.71191.30020.20900.19997.4

6.5

 経過時間が短いサブルーチンほど、タイマーのオーバーヘッドの影響が大きくなる14

Copyright 2018 FUJITSU LIMITEDforce_bondのチューニング
 チューニング方針
1. 手動でループアンスイッチングを適応
2. !ocl simd_redundant_vl()を配列式の直前に挿入 この指示行を指定しないとSWP化されない。コンパイラのメッセージは、「スケジューリング結果を得られなかったため、ソフトウェアパイプライニングを適用できません。」
simd_redundant_vl指示行
ループ回転数がSIMD長で割り切れない場合でも、マスク付きSIMD命令を利用し、SIMD実行するループを生成する。
57
58
59
60

1
1
1
1

61
62
63
64
65
66
67

2 p
2
2 p
2 p
2
2
2

if(inmgo%i_multi_mgo == 0) then 手動ループアンスイッチング!$omp do private(imp1,imp2,v21,dist,ddist,ddist2,for, &!$omp&
force,efull,iunit,junit,isys,istat)<<< Loop-information Start >>><<< [OPTIMIZATION]<<< SOFTWARE PIPELINING<<< Loop-information End >>>do ibd = ksta, kend
imp1 = ibd2mp(1, ibd)
imp2 = ibd2mp(2, ibd)

! write(*,*) "fj*** imp1=",imp1,"imp2=",imp2!ocl simd_redundant_vl(3)<<< Loop-information Start >>><<< [OPTIMIZATION]<<< SIMD(VL: 4)<<< FULL UNROLLING<<< Loop-information End >>>2 p 1v
v21(1:3) = xyz_mp_rep(1:3, imp2,irep) - xyz_mp_rep(1:3,

68
imp1,irep)
69 2
70 2 p71 2 p72 2 p
73 2

74
75
76
77

2
2 p
2
2

78
79
80
81
82
83
84
85
86

2
2
2
2
2
2
1
1
1

p
p
p
p

! calc force
for = (coef_bd(1, ibd) + 2.0e0_PREC * coef_bd(2, ibd) * ddist2) * &(-2.0e0_PREC * ddist / dist)!ocl simd_redundant_vl(3)<<< Loop-information Start >>><<< [OPTIMIZATION]<<< SIMD(VL: 4)<<< FULL UNROLLING<<< Loop-information End >>>1v
force(1:3) = for * v21(1:3)!ocl simd_redundant_vl(3)1vforce_mp(1:3, imp1) = force_mp(1:3, imp1) - force(1:3)!ocl simd_redundant_vl(3)1vforce_mp(1:3, imp2) = force_mp(1:3, imp2) + force(1:3)
enddo!$omp end do nowaitelse 手動ループアンスイッチング

dist = sqrt(v21(1)**2 + v21(2)**2 + v21(3)**2)
ddist = dist - bd_nat(ibd)ddist2 = ddist**2

15

Copyright 2018 FUJITSU LIMITEDforce_bondのチューニングの効果
 実行時間内訳の比較
 チューニング効果(#) タイマー(gettod)の経過時間から

force_bond

Original

Tune

向上率

0.4509秒

0.4120秒

9.4%①① 命令コミットのコストは少し改善
 ループ内の配列式をSIMD化したので命令数が
減少した。

②

② 命令スケジューリングで改善できる部分のSWPによる効果が小さい- SWPに必要な回転数 :44(コンパイラメッセージから)- 実際のループ回転数 : 約900
Original

Tune
16

Copyright 2018 FUJITSU LIMITEDforce_bondのSWPの効果
 SWPの効果が小さい = 命令を重ねられていない原因
 ループ最後のインダイレクトアクセスのあるリダクション演算部に対して、
コンパイラは依存関係を見切れず、逐次スケジューリングになっていると考えられる。
79
80
81
82

2
!ocl simd_redundant_vl(3)2 p 1v
force_mp(1:3, imp1) = force_mp(1:3, imp1) - force(1:3)
2!ocl simd_redundant_vl(3)2 p 1v
force_mp(1:3, imp2) = force_mp(1:3, imp2) + force(1:3)

 しかし、データには重なりがあるため、!OCL NORECURRENCEが使えない
fj*** imp1= 1 imp2= 2
fj*** imp1= 1 imp2= 3
fj*** imp1= 3 imp2= 4
fj*** imp1= 4 imp2= 5
fj*** imp1= 4 imp2= 6
fj*** imp1= 6 imp2= 7
fj*** imp1= 7 imp2= 8
fj*** imp1= 7 imp2= 9
fj*** imp1= 9 imp2= 10
fj*** imp1= 10 imp2= 11:

データに重なりが無いようなレイアウトにすることでコンパイラ最適化(SWP)の効果を最大限に引き出せると考える
17

Copyright 2018 FUJITSU LIMITEDまとめ FX100とSkylakeの5~7.4倍の性能差の分析
 タイマールーチンの置換えにより、性能差は3.4~6.5倍
 さらに、ピーク演算性能の差を考慮すると、その差は、1.7~3.3倍
 FX100では、PA情報の採取から、命令数と演算/データアクセス待ちがボトル
ネックとなっていた。

 高速化の検討
 手動ループアンスイッチングと!ocl simd_redundant_vlの指定により、SIMD化/SWP化を促進することで、約10%性能が改善した。
 インダイレクトアクセスのデータに対して、データの重なりがないようレイアウトを変更できれば、より最適化の効果を得ることができると考える。

 残り2ルーチンに対して
 _force(bond)との一番の違いはcycle文があること。
 cycle文直前でループ分割し、後半のループに対して、_force(bond)と同様のチューニングをすることが可能
18
Copyright 2018 FUJITSU LIMITED2019/2/25
SS研メニーコア時代のアプリ性能WG
第9回会合

生体分子粗視化シミュレータCafeMol
タイマーとローカルな力のSIMD化

検崎博生
理研 情報システム部

概要• 粗視化モデルとCafeMolについて• 生体分子シミュレーションと粗視化• CGタンパク質/DNAモデル• 計算方法と並列化• タイマーによるオーバーヘッド• mpi_wtime, system_clock, gettod• 1ノードでのMPI並列数とスレッド並列数• ローカルな相互作用のSIMD化• ボンド長とボンド角相互作用• OpenMPのスレッド並列のチャンクサイズを変更

Model
H.Kenzaki, et al,
J. Chem. Theor. Chem. (2011)• CafeMol• Protein• AICG2+ model W Li, et al PNAS (2012)• FLP model for disordered regions in X-ray structure
Knotts et al, J. Chem. Phys. (2007)• DNA• 3SPN.1 model

ローカルなポテンシャル

: bond angle
: dihedral angle(0 means native state) K  1  cos     K  1  cos 3   

Vlocal  K b  ri ,i1  r0 i ,i 1  K   i   0 i
2

i

2

i

1

3



i

i

0i



i

0i

i

ボンド長、ボンド角、二面角などの相互作用からなるが、モデルによって詳細はさまざまである。

今回の計算対象と計算方法• 計算対象• ヌクレオソーム1個の系：1,854粒子• 計算方法• ネイバリングリストを作っておく• ローカルな相互作用は最初に作っておく• それ以外は100stepに1回更新• MPIとOpenMPでネイバリングリストを分割して並列計算• プロセスとスレッド毎に力を保存して、最後に通信を行い足し合わせる• 速度や座標のアップデートは全プロセスで同じものを行う。
ネイバリングリストの例 (2体の相互作用)
1
i

1

j

2

2

3

4

5

6

2
4

7

9

5

7

8

9

3
8

4

10

11

4
5

9

7

12

13

14

5
8

6

15

16

6
7

9

8

17

18

7
9

8

19
8

9

9

タイマーの比較 nucleosome (105step, 100mM)• FX100 1 node: SPARC64 XIfx (1.975GHz, 32 cores, 1 CPU/node)• コンパイラは、Fujitsu Fortran Compiler 2.0.0-07• 最適化オプションは-Kopenmp,fast,parallel,ocl,optmsg=2 -Qa,m,p,t,x• 1ノード（32コア）で2MPI並列X16スレッド並列で実行
FX100(32 cores) mpi̲wtime FX100(32 cores) system̲clockFX100(32 cores) gettod
force

38.33

34.92

34.94

̲force(comm)

6.12

6.15

6.19

̲force(local)

10.42

8.72

8.68

̲̲force(bond)

0.48

0.42

0.38

̲̲force(bangle)

1.54

1.39

1.41

̲force(go)

1.43

1.31

1.30

̲force(pnl)

11.70

11.55

11.51

̲force(ele)

5.87

5.69

5.68

ope

60.52

53.03

52.75

comm

6.75

6.60

6.65

main̲loop

67.27

59.62

59.40

mpi_wtimeはオーバーヘッドが大き過ぎる。
system_clockとgettodは同程度。

FX100での性能測定 nucleosome (105step, 100mM)• FX100 1 node: SPARC64 XIfx (1.975GHz, 32 cores, 1 CPU/node)• コンパイラは、Fujitsu Fortran Compiler 2.0.0-07• 最適化オプションは-Kopenmp,fast,parallel,ocl,optmsg=2 -Qa,m,p,t,x• 1ノード（32コア）でMPI並列数X16スレッド並列数=32に固定して実行

1nodeでは2MPI並列x16OpenMP並列のときに一番高い性能

force(bond)の計算方法
1. ネイバリングリストiele2mp(2, lele)を各プロセスで計算
i < jとなるペアだけを格納

2. simu_force.F90!$omp parallel private(tn)call simu_force_bond(force_mp_l(1, 1, tn))

スレッド毎にforce_mpの異なる領域を用意

3. simu_force_bond.F90!$omp do private(....)
ibd=1 imp1=1 imp2=2
ibd=2 imp1=1 imp3=3
do ibd = ksta, kend
ibd=3 imp1=3 imp2=4
imp1 = ibd2mp(1, ibd)
ibd=4 imp1=4 imp2=5
imp2 = ibd2mp(2, ibd)
ibd=5 imp1=4 imp2=6
v21(1:3) = xyz_mp_rep(1:3, imp2,irep) - xyz_mp_rep(1:3, imp1,irep)
ibd=6 imp1=6 imp2=7
dist = sqrt(v21(1)**2 + v21(2)**2 + v21(3)**2)
ddist = dist - bd_nat(ibd)ddist2 = ddist**2
for = (coef_bd(1, ibd) + 2.0e0_PREC * coef_bd(2, ibd) * ddist2) * (-2.0e0_PREC * ddist / dist)
force(1:3) = for * v21(1:3)
force_mp(1:3, imp1) = force_mp(1:3, imp1) - force(1:3)
force_mp(1:3, imp2) = force_mp(1:3, imp2) + force(1:3)
end doforce_mpへの足しこみが!$omp end do nowait

原因でSIMD化できないforce(bond)のSIMD化• 方針• force_mpへの足し込みをループ内で被らないようする。
• force_mpがスレッド毎に別なので、各スレッド毎に被らないようする。
• ボンド長の場合、ネイバリングリストで3つ以上離れていれば被らない。
• よって、チャンクサイズを1にし、スレッド並列数が3以上であれば、SIMD化が可能なはず。

• tune1• 今回のモデルに関係ない部分を削除しておく。
• チェック用のifなど、if文は全て削除する。

• tune2• !omp do schedule(static,1)としてチャンクサイズを1に設定する。
• この段階では、メモリアクセスが遅くなると予想できる。

• tune3• !ocl norecurrence()を付けてSIMD化を行う。

• ボンド各force(bangle)についても同様にSIMD化を行うが、FLPモデルについては、中で関数を読んでいるのでSIMD化に失敗。
• 他の相互作用でもif文を外せないなどSIMD化できないことはよくある。

force(bond)とforce(bangle)のSIMD化• tune1: 今回の計算に使わない部分を削除しソースを単純化• tune2: OpenMPのスレッド並列の順序を変更• tune3: !ocl norecurrence()を付けてSIMD化• コンパイルオプションは-Kopenmp,fast,parallel,ocl,optmsg=2 -Qa,m,p,t,x• 1ノード（32コア）で2MPI並列X16スレッド並列で実行
FX100(32 cores) tune1

FX100(32 cores) tune2

FX100(32 cores) tune3

force

34.90

35.44

35.00

̲force(comm)

6.18

6.21

6.22

̲force(local)

8.71

8.99

8.79

̲̲force(bond)

0.37

0.62

0.51

̲̲force(bangle)

1.42

1.44

1.06

̲force(go)

1.27

1.29

1.29

̲force(pnl)

11.52

11.65

11.47

̲force(ele)

5.68

5.67

5.68

ope

52.80

53.24

52.80

comm

6.65

6.66

6.65

main̲loop

59.44

59.90

59.44

bondはSIMD化の副作用の方が大きく遅くなった。
bangleはSIMD化により30%程度速くなった。

まとめ• Timerの比較• MPI_WTIMEはオーバーヘッドが大きい。
• system_clockとgettodは同程度。

• MPI並列とOpenMP並列の兼ね合い。
• 1ノードでは、2MPI並列x16スレッドの時に一番高い性能。

• FX100でのローカルな相互作用のSIMD化• スレッド並列の割当を変えることにより、SIMD化に成功。
• 複雑な力場を使う場合はif文や関数呼び出しを含むので難しい。

• 性能の向上は限定的。
• スレッド並列の割当を変えるオーバーヘッドとの兼ね合い。

3.4

ADVENTURE における多倍長精度演算の利用に向けた検討
名古屋大学情報基盤センター 荻野 正雄

3.4.1

はじめに

ADVENTURE と は 、設 計 用 大 規 模 計 算 力 学 シ ス テ ム 開 発 プ ロ ジ ェ ク ト (通 称 ADVENTURE プ ロ ジ ェ ク ト 、 https://adventure.sys.t.u-tokyo.ac.jp/)に お い て 開 発
さ れ て い る オ ー プ ン ソ ー ス CAE シ ス テ ム で あ る 。 大 規 模 メ ッ シ ュ を 用 い て 自 然
物や人工物を丸ごと詳細にモデル化し、多様な並列分散計算機環境のもとで固体の変形や熱・流体の流れ等の力学解析から可視化、設計最適化までを行える
特徴がある。無料公開されていることもあり、教育・研究から産業応用まで幅
広く利用されているソフトウェアである。ここでは、有限要素法による 3 次元
電 磁 場 解 析 モ ジ ュ ー ル ADVENTURE_Magnetic を 取 り 上 げ る 。 電 磁 場 に は 静 磁 場 、
低周波電磁場、高周波電磁場など様々な種類があり、対象とする問題の特性に応 じ て Maxwell 方 程 式 か ら 偏 微 分 方 程 式 を 導 出 す る 。 そ れ ら に 辺 要 素 有 限 要 素法を適用することで、様々な線形方程式が得られる。例えば、時間調和渦電流問題や高周波電磁場問題では複素対称線形方程式を解くことになる。しかし、係 数 行 列 が 悪 条 件 に な り や す く 、 共 役 直 交 共 役 勾 配 (COCG)法 な ど の 反 復 解 法 で
高 い 精 度 の 収 束 解 を 得 る こ と が 困 難 な 場 合 が あ る 。そ こ で 現 在 、IEEE 754-2008
が 規 定 す る 四 倍 精 度 浮 動 小 数 点 数 や D.H. Bailey や K. Briggs な ど が 提 案 し た
倍精度数 2 つを用いる疑似四倍精度数などの多倍長精度演算に着目し、反復計算における丸め誤差を低減させることで、電磁場解析を効率化することを目指し て い る 。 今 回 は 、 富 士 通 FX100 を 含 む い く つ か の 環 境 に お い て 、 任 意 精 度 ・
多倍長精度ライブラリの性能評価を行った。また、ライブラリ利用における知識 や 技 術 の 蓄 積 と 共 有 を 行 い 、 さ ら に 富 士 通 TCS 環 境 の 課 題 を 明 ら か に す る 。
3.4.2

任意精度・多倍長精度演算ライブラリの概要

表 3.4.1 に 今 回 用 い た 任 意 精 度 ・ 多 倍 長 精 度 演 算 ラ イ ブ ラ リ 、 実 数 及 び 複 素
数の浮動小数点数型を示す。
表 2.9.1

主な任意精度・多倍長精度浮動小数点数

Library

Precision

Real data type

Complex data type

C data typeextendedlong double

long double _Complex

Fortran data type

quadruple

real(kind=16)

complex(kind=16)

libquadmath

quadruple__float128__complex128

Intel's _Quadquadruple_Quadcomplex<_Quad>

QD

pseudo-quad

dd_realcomplex<dd_real>

Fujitsu's fast_dd

pseudo-quad

dd_realcomplex<dd_real>

ARPREC

arbitrarymp_realmp_complex

exflibarbitrary

exfloatcomplex<exfloat>

mpfr_t

-

-

mpc_t

MPFR/GMP
MPC

arbitrary

四 倍 精 度 と し て は 、 Fortran デ ー タ 型 、 GCC 拡 張 で あ る C 言 語 の libquadmath
ラ イ ブ ラ リ (https://gcc.gnu.org/onlinedocs/libquadmath/)、 Intel コ ン パ イ
ラ 拡 張 で あ る C 言 語 の _Quad を 用 い る 。libquadmath ラ イ ブ ラ リ は GCC 4.6.0 以
降 か ら 利 用 可 能 で あ る 。 QD ラ イ ブ ラ リ (http://crd-legacy.lbl.gov/~dhbailey/mpdist/)は 倍 精 度 浮 動 小 数 点 数 を 2 つ 用 い る double-double 形 式 の 疑 似 四 倍 精
度 演 算 ラ イ ブ ラ リ で あ る 。 富 士 通 fast_dd は 、 富 士 通 TCS に 含 ま れ る 高 速 4 倍
精 度 基 本 演 算 ラ イ ブ ラ リ で あ り 、 QD と 同 じ く double-double 形 式 の 疑 似 四 倍 精
度 で あ る 。比 較 用 と し て 、C デ ー タ 型 の 拡 張 倍 精 度 、任 意 精 度 演 算 ラ イ ブ ラ リ で
あ る ARPREC(http://crd-legacy.lbl.gov/~dhbailey/mpdist/)、 exflib(http://www-an.acs.i.kyoto-u.ac.jp/~fujiwara/exflib/)、 並 び に MPFR(http://www.
mpfr.org/)と そ の 複 素 数 版 で あ る MPC(http://www.multiprecision.org/)を 用いる。
こ こ で 、 Intel _Quad、 QD、 fast_dd、 並 び に exflib で は 複 素 数 の 基 本 演 算 機
能 が 提 供 さ れ て い な い た め 、C++複 素 数 テ ン プ レ ー ト ク ラ ス と 組 み 合 わ せ て 用 い
た。
また、任意精度・多倍長精度演算ライブラリは利用マニュアルが十分でないことが多く、サンプルコードからプログラミング方法を学ぶ必要がある。今回
作 成 し た い く つ か の プ ロ グ ラ ム を ス ラ イ ド 1～ 8 に 示 す 。
3.4.3 測 定 環 境
性 能 測 定 に は 、 汎 用 PC、 名 古 屋 大 学 の FX100、 及 び 名 古 屋 大 学 の CX400 を 使
用 し た 。そ れ ぞ れ の 環 境 に お け る CPU、OS、コ ン パ イ ラ を 表 2.9.2 に 示 す 。ま た 、
そ れ ぞ れ の 環 境 に お い て 性 能 測 定 を 行 う こ と が で き た ラ イ ブ ラ リ を 表 2.9.3 に示す。性能測定できなかった理由は次の通りである。
・libquadmath は 、 FX100(fcc)は 富 士 通 コ ン パ イ ラ の た め 、 FX100(gcc)は gcc6.3.0 だ が 当 該 ラ イ ブ ラ リ 未 提 供 の た め 、 CX400 は gcc-4.4.7 の た め・Intel's _Quad は 、 PC で は ラ イ セ ン ス 未 導 入 の た め 、 FX100 は SPARC64 プ ロ
セッサのため・

fast_dd は 、 FX100 の み 利 用 可 能 で あ る た め・exflib は x64 バ イ ナ リ 版 の み を 利 用 し た た め

表 3.4.2
Name
PC

実験に用いた計算機とソフトウェア
CPUOS

Compiler

SiteIntel Core i7-8650U (Kabyla

Ubuntu 18.04

gcc-7.3.

-

ke-R)

TLS

0

FX100(fcc)

Fujitsu SPARC64 XIfx

FX100(gc

XTC OS 2.1.1

fcc-2.0.

Nagoy

0

a

gcc-6.3.

Nagoy

0

a

icc-14.

Nagoy

0.3

a

c)
CX400

Intel Xeon E5-2697v3 (Haswe
ll-EP)
表 3.4.3

Library

RHEL

使用ライブラリと計算環境における利用可否
PC

FX100(fcc)

FX100(gcc)

CX400

〇

-

-

-

icc-14.0.3

-

-

-

〇

QD

2.3.20〇〇〇〇

fast_dd

2.3.2

-

〇

〇

-

ARPREC

2.2.19〇〇〇〇exflibx64-bin-20180620

〇

-

-

〇

GMP

4.0.1

MPFR

6.1.2

〇

〇

〇

〇

MPC

1.1.0

libquadmath

Version(gcc-4.6.0 or later)Intel's_Quad

3.4.4 計 算 精 度 に 関 す る 数 値 実 験a) 問 題 設 定
実 数 を 係 数 に も つ 2 次 方 程 式 𝑎𝑥 2 + 𝑏𝑥 + 𝑐 = 0を 考 え る 。こ の 式 は 解 の 公 式 に より 、𝑥1 = (−𝑏 − √𝑏2 − 4𝑎𝑐)/(2𝑎)と 𝑥2 = (−𝑏 + √𝑏2 − 4𝑎𝑐)/(2𝑎)で 得 ら れ る 。こ の と き 、
桁 落 ち し や す い 係 数 𝑎 = 1.01、 𝑏 = 2718281、 𝑐 = 0.01を 用 い て 𝑥2 を 計 算 し た と き の精度を評価する。この問題では加減算、乗算、除算、平方根の計算を行う。
真 の 解 は 𝑥2 = −3.678795532912165323920499 × 10−9で あ る 。 比 較 と し て 、 exflib ラ イ ブ ラ リ で 100 桁 精 度 で 計 算 し た 結 果 は 𝑥̃2 = −3.
6787955329121653239204992803347 × 10−9で あ っ た 。
b) 実 験 結 果
ス ラ イ ド 9 に 実 験 結 果 を 示 す 。任 意 精 度 演 算 ラ イ ブ ラ リ に つ い て は 、4 倍 精 度相当となるように計算精度を設定している。表の結果に書かれている数値は、太字・下線のところまで正しく計算できていることを表している。表より、普通に倍精度数で計算したのでは小数点以下 1 桁までしか正しく計算できていな
い こ と が 分 か る 。 こ れ に 対 し 、 疑 似 四 倍 精 度 の QD と fast_dd は 16～ 17 桁 ま で
正 し く 計 算 で き て お り 、Fortran や C の libquadmath と い っ た 四 倍 精 度 の 結 果 (18～ 19 桁 )に 近 い 計 算 精 度 が 得 ら れ て い る 。 こ れ よ り 、 名 大 FX100 環 境 な ど に おいて、多倍長精度計算ライブラリを正しく利用できていることが示された。
こ こ で 、 x86 系 で あ る PC の 環 境 で は long double 型 は 小 数 点 以 下 3 桁 程 度 の精度で拡張倍精度として倍精度より高い精度で計算しているようであるがこの
問 題 に と っ て は 精 度 不 足 で あ る 。一 方 、FX100 の long double 型 は 四 倍 精 度 相 当の 結 果 が 得 ら れ て い る 。long double 型 変 数 が 占 め る サ イ ズ を 調 べ た と こ ろ 128
ビ ッ ト で あ り 、き ち ん と 128 ビ ッ ト 分 の 精 度 が 出 て い る こ と に な る 。富 士 通 の C
言語処理系でとりあえず四倍精度数を使いたい場合に有用である。
ま た 、Intel's _Quad は 倍 精 度 と 同 じ 程 度 の 結 果 と な っ た 。こ れ は 、平 方 根 を
含 め 数 学 関 数 の 四 倍 精 度 版 が 提 供 さ れ て い な か っ た た め で あ る 。 実 験 当 時 の In
tel コ ン パ イ ラ の 説 明 と し て _Quad は 実 験 的 な も の と 述 べ ら れ て い た 。
な お 、 任 意 精 度 演 算 ラ イ ブ ラ リ の ARPREC と exflib は 想 定 以 上 の 精 度 が 出 ているが理由は不明である。
3.4.5 計 算 時 間 に 関 す る 数 値 実 験a) 問 題 設 定
以下のロジスティック写像を考える。
𝑥𝑛+1 = 𝑎𝑥𝑛 (1 − 𝑥𝑛 )𝑎 = 4と す る と 、𝑥𝑛 は 区 間 [0,1]全 体 で 非 周 期 に 変 動 す る カ オ ス と な る こ と が 知 ら
れ て い る 。 こ の と き 、 初 期 値 𝑥0 = 0.7501と し 、 𝑛が 106に な る ま で 計 算 し た と き の
計 算 時 間 を 評 価 す る 。ま た 、𝑥𝑛 を 実 数 と し た と き 、複 素 数 と し た と き の 2 ケ ー ス
を評価する。この問題では、乗算と減算のみを行う。
b) 実 験 結 果
ス ラ イ ド 10 と 11 に 実 数 と 複 素 数 の と き の 実 験 結 果 を そ れ ぞ れ 示 す 。 表 に お
い て 、 QD(C++)と QD(C)は そ れ ぞ れ プ ロ グ ラ ム 記 述 言 語 が C++と C で あ る こ と を
表 し て い る 。 QD ラ イ ブ ラ リ 自 体 は C++で 記 述 さ れ て お り 、 C++プ ロ グ ラ ム か ら 利
用 す る 場 合 向 け に 四 則 演 算 等 の 演 算 子 オ ー バ ー ロ ー ド を し て い る 一 方 で 、C プ ロ
グ ラ ム 向 け に 関 数 コ ー ル で の 利 用 も 可 能 と な っ て い る 。C++で あ れ ば プ ロ グ ラ ム
は 書 き や す い が 、計 算 時 間 の 観 点 で 最 適 な コ ー ド が 生 成 さ れ な い 可 能 性 が あ る 。
そこで比較のために 2 つの実装を行った。
ス ラ イ ド 10 よ り 、 PC で は QD(C++)が 最 速 で あ っ た 。 こ れ は C++の 最 適 化 が 十
分 に 行 え て い る の だ と 推 測 さ れ る 。 ま た 、 QD(C++)に よ る 疑 似 四 倍 精 度 演 算 は 倍
精度演算に比べて 3 倍遅い程度であった。また、疑似四倍精度演算は四倍精度演算よりも 3 倍程度速く、前の実験で計算精度が同程度であったことから、実用 性 が 高 い と 言 え る 。FX100(fcc)で は 、同 じ く QD(C++)が 最 速 で あ っ た が 、倍 精

度 演 算 よ り も 7 倍 程 度 の 計 算 時 間 で あ り 、C++コ ー ド の さ ら な る 最 適 化 が 望 ま れる 。 CX400 で は Intel's _Quad が 最 速 で あ り 、 QD(C++)も ほ ぼ 同 程 度 で あ っ た 。
次 に 、 タ ー ゲ ッ ト と す る 複 素 数 演 算 の 性 能 を 、 ス ラ イ ド 11 で 見 て み る 。 実 数の と き と 異 な り 、 全 て の 環 境 に お い て QD(C)が 最 速 と な っ た 。 よ っ て 、 QD(C)の利 用 が 良 い よ う に 思 え る が 、プ ロ グ ラ ム 例 の ス ラ イ ド 4 を 見 て わ か る よ う に 、Q
D ライブラリを C 言語から複素数で利用する場合は複素数の四則演算などの関数を自作する必要があり、さらにプログラミングコストが高くなり、可搬性も損な わ れ る 。 C++の 演 算 子 オ ー バ ー ロ ー ド を 利 用 す れ ば 、 疑 似 四 倍 精 度 数 で あ っ ても 、複 素 数 で あ っ て も 、z=x+y の よ う に 記 述 で き る 。し か し 、実 際 に は 、実 部 と虚部それぞれの加算に展開され、加算は疑似四倍精度数の加算になり、その加
算 1 回 あ た り で は Knuth に よ る 加 算 の 無 誤 差 変 換 に 基 づ い て 倍 精 度 数 の 加 減 算が 11 回 行 わ れ る 。 つ ま り 、 z=x+y と い う 記 述 は 22 回 の 演 算 に 展 開 さ れ る 。 z=x*y は 100 回 以 上 、 z=x/y は 250 回 以 上 の 演 算 に も な る 。 実 際 の プ ロ グ ラ ム で はも っ と 長 い 計 算 式 が 書 か れ 、for ル ー プ ブ ロ ッ ク 内 で あ れ ば ル ー プ ボ デ ィ が 非 常に 長 く な る こ と を 意 味 す る 。 ス ラ イ ド 14 と 15 に 、 QD(C++)を 用 い た 疑 似 四 倍 精
度 の 密 行 列 ベ ク ト ル 積 コ ー ド を FX100 の 富 士 通 C/C++コ ン パ イ ラ で コ ン パ イ ル した と き の 最 適 化 メ ッ セ ー ジ を 示 す 。 こ れ よ り 、 SIMD 化 も ソ フ ト ウ ェ ア パ イ プ ラ
イン化も行われていないことが分かる。コンパイラによる最適化が望まれる。
なお、インライン展開に関するメッセージが大量にでてくるため、メッセージを調査しづらい問題があった。
ま た 、 名 大 FX100 に お い て 富 士 通 C/C++コ ン パ イ ラ の 32 レ ジ ス タ 制 限 版 を 用
い た 実 験 も 行 っ た 。 そ の 結 果 を ス ラ イ ド 12 と 13 に 示 す 。 こ の ケ ー ス で は 通 常
環境と性能差はなかった。しかし、上述したように現状でもコンパイラにはもう 少 し 最 適 化 を 頑 張 っ て ほ し い た め 、 32 レ ジ ス タ 制 限 環 境 の 影 響 を 受 け な い かどうかは判断できないと言える。
3.4.6 複 素 対 称 線 形 方 程 式 の 数 値 実 験複素対称線形方程式における多倍長精度計算の有効性を評価する数値実験を
PC 環 境 で 行 っ た 。 疎 行 列 デ ー タ ベ ー ス UF Sparse Matrix Collection で 公 開 され て い る 2 つ の 小 規 模 行 列 を 用 い て 、IC 前 処 理 付 き COCG 法 で 解 い た 実 験 結 果 をス ラ イ ド 16～ 18 に 示 す 。 こ れ よ り 、 疑 似 四 倍 精 度 数 を 用 い る と 倍 精 度 数 の 場 合と比べて反復回数を削減できることが分かる。このケースでは総計算時間は倍精度の場合が高速であったが、反復回数の削減効果は問題依存のため、その効
果が大きい場合であれば、倍精度数よりも疑似四倍精度数を使う方が高速になる可能性があると言える。
ス ラ イ ド 19 と 20 に 、 高 周 波 電 磁 場 の 有 限 要 素 解 析 で 得 ら れ た 複 素 対 称 行 列による数値実験結果を示す。このケースでは、わずかにであるが疑似四倍精度が高速であった。

3.4.7 富 士 通 fast_dd に つ い て
富 士 通 TCS に 含 ま れ る 富 士 通 fast_dd は QD と 同 様 に double-double 形 式 の 疑
似 四 倍 精 度 演 算 機 能 を 提 供 す る ラ イ ブ ラ リ で あ る が 、ス ラ イ ド 10～ 13 に 示 す よう に 計 算 時 間 の 面 で は QD ラ イ ブ ラ リ よ り も 低 い 性 能 と な っ て い る 。 し か し 、 富士 通 fast_dd は 2 要 素 ず つ 同 時 に 計 算 す る マ ル チ 関 数 や ベ ク ト ル 関 数 を 提 供 しており、それらが利用できるケースでの性能改善が目的の 1 つにある。前述の計算は単純な漸化式であったことからマルチ関数等の利用は行っていない。そこ で 、 第 9 回 WG に お い て 富 士 通 側 報 告 と し て 行 わ れ た マ ル チ 関 数 ・ ベ ク ト ル 関
数 の 性 能 評 価 例 を ス ラ イ ド 21～ 23 に 示 す 。 こ れ よ り 、 複 数 要 素 の 同 時 計 算 が 行
え る 場 合 で は fast_dd は QD よ り も 高 い 性 能 が 期 待 で き る こ と が 明 ら か に な っ た 。
3.4.8 HOWTO
今回の数値実験を通して得られた多倍長精度演算ライブラリの利用に関する
知 識 に つ い て 、HOWTO と し て ま と め 、知 識 共 有 を 行 う 。得 ら れ た 知 識 は 以 下 の 通りである。
① 任 意 精 度・多 倍 長 精 度 演 算 ラ イ ブ ラ リ に は 、計 算 順 序 の 入 れ 替 え を 行 う と 計
算 精 度 が 低 下 す る も の が あ る 。 よ っ て 、 例 え ば 富 士 通 C/C++コ ン パ イ ラ で は"-Keval"オ プ シ ョ ン は 指 定 す べ き で は な い 。"-Keval"は "-Kfast"に 含 ま れ て
い る の で "-Kfast"の 利 用 も 避 け る べ き で あ る 。し か し 、環 境 に よ っ て は 、コン パ イ ラ の 標 準 オ プ シ ョ ン と し て "-Kfast"が 指 定 さ れ て い る 場 合 が あ る ので 注 意 が 必 要 で あ り 、 "-Knoeval"を 使 う の が 望 ま し い 。
② 富 士 通 C/C++処 理 系 で 任 意 精 度 ・ 多 倍 長 精 度 演 算 ラ イ ブ ラ リ を 利 用 す る と き
は 、 コ ン パ イ ラ の イ ン ラ イ ン 展 開 オ プ シ ョ ン "-x-"を 指 定 し た 方 が 高 速 と なる。
③ 富 士 通 C/C++コ ン パ イ ラ で QD ラ イ ブ ラ リ を ビ ル ド す る と き は 、 コ ン パ イ ラに "-Kfp_contract"オ プ シ ョ ン を 指 定 す れ ば FMA 命 令 を き ち ん と 使 っ て く れる 。 ① 、 ② と あ わ せ て 、 "-O3 -Kfp_contract-Knoeval -x-"が 基 本 と な る 。
④ 名 大 FX100 の gcc-6.3.0 環 境 で は FMA 命 令 を 使 う よ う 指 定 し て QD ラ イ ブ ラ
リ を ビ ル ド す る と 、計 算 結 果 が 誤 っ た も の に な る 。SPARC64 XIfx の 正 し い F
MA 命 令 が 出 さ れ な い の だ と 思 わ れ る 。ま た 、libquadmath は ビ ル ド さ れ て いない。
⑤ 富 士 通 fast_dd は SSL II の 一 部 と し て 提 供 さ れ て い る が 、GCC な ど 富 士 通 TCS 以 外 の コ ン パ イ ラ で SSL II ラ イ ブ ラ リ を 使 う と き は "-SSL2"で は リ ン ク

で き な い 。 例 え ば 今 回 指 定 し た の は 、 イ ン ク ル ー ド パ ス 指 定 は "-I/opt/FJS
Vmxlang/include"、 リ ン ク 指 定 は "-L/opt/FJSVmxlang/lib64 -lssl2mtfour
simd_prexi -lssl2mtfoursimd_com -lssl2mtfoursimd_postxi -lfj90i -l
fj90fmt-lfj90f -lfj90rt -lfjrtcl -ltrt -L/usr/lib64 -lm -lelf"で あ る 。Intel MKL Link Line Advisor の よ う な リ ン ク 方 法 を 調 べ る 手 助 け が あ る の が望ましい。
⑥ 富 士 通 C/C++処 理 系 で は long double 型 は 128 ビ ッ ト の 4 倍 精 度 浮 動 小 数 点数である。
⑦ 任 意 精 度・多 倍 長 精 度 の 浮 動 小 数 点 数 に 定 数 を 代 入 す る と き は 、定 数 の 精 度に 気 を 付 け る 必 要 が あ る 。 例 え ば QD ラ イ ブ ラ リ で は 文 字 列 定 数 を 使 え ば 、
疑似四倍精度数に変換した上で代入してくれる。
⑧ 富 士 通 fast_dd は 疑 似 四 倍 精 度 定 数 の 指 定 方 法 が 不 明 だ が 、富 士 通 C/C++処
理 系 の と き は 4 倍 精 度 で あ る long double 型 定 数 で 代 用 で き る 。
⑨ 任 意 精 度・多 倍 長 精 度 演 算 ラ イ ブ ラ リ を 使 う と き は 、出 力 関 数 が 対 応 し て い
るかも気にする必要がある。
3.4.9 ま と め
今 回 は 将 来 的 に ADVENTURE へ 多 倍 長 精 度 演 算 を 組 み 込 む こ と を 目 的 に 、 任 意
精 度 ・ 多 倍 長 精 度 演 算 ラ イ ブ ラ リ を 名 大 FX100 環 境 に 移 植 し 、 数 値 実 験 に よ って性能評価を行った。当初は、想定する計算結果が得られないことが多く、丸め 回 数 が 少 な い FMA 命 令 を 使 え て い る か 、 コ ン パ イ ラ に よ る 計 算 順 序 入 れ 替 えは起こっていないか、高精度な定数はどう指定するか、高精度な数学関数はあ
るのか、あたりを確認していけば正しく移植できることが見えてきた。疑似四倍精度数かつ複素数となると 1 回の四則演算あたりの計算量が多く、それを含む DO ル ー プ ／ for ル ー プ の 最 適 化 は 難 し い と 思 わ れ る が 、 ポ ス ト 京 の コ ン パ イ
ラが効率良いコードを出してくれることに期待する。また、多倍長精度計算の専門家との協調作業が必要と思われるが、そのためには大量に出力されるコンパイラの最適化メッセージの改善も望まれる。

複素数版プログラム例 (FORTRAN)
program quadformula_complex_real16
implicit none
4倍精度複素数の宣言
complex(kind=16) :: a, b, c, x2
a = 1.01Q0
4倍精度定数の代入
b = 2718281.0Q0
c = 0.01Q0x2 = (-b+sqrt(b*b-4.0*a*c))/(2.0*a)
print *, real(x2)end program quadformula_complex_real161

複素数版プログラム例 (libquadmath, C)#include <stdio.h>#include <stdlib.h>#include <quadmath.h>int main(void){__complex128 a, b, c, x2;
char str[1024];

4倍精度複素数の宣言

4倍精度定数の代入
a = 1.01Q;
b = 2718281.0Q;
4倍精度数学関数の呼び出し
c = 0.01Q;
x2 = (-b+csqrtq(b*b-4.0*a*c))/(2.0*a);quadmath_snprintf(str, 1024, "%31.30Qe , %31.30Qe",__real__ x2, __imag__ x2);
4倍精度専用の出力関数
printf("x2= %s¥n", str);}2

複素数版プログラム例 (QD, C++)
#include <iostream>#include <iomanip>#include <complex>#include <stdlib.h>#include <qd/dd_real.h>using namespace std;
int main(void){
疑似4倍精度変数の宣言
complex<dd_real> x2;
unsigned int oldcw;
80ビット拡張倍精度を使うIntel x86向け設定
fpu_fix_start(&oldcw);

4倍精度定数の代入
complex<dd_real> a("1.01", "0.0");
complex<dd_real> b("2718281.0", "0.0");
complex<dd_real> c("0.01", "0.0");x2 = (-b+sqrt(b*bcomplex<dd_real>(4.0,0.0)*a*c))/(complex<dd_real>(2.0,0.0)*a);cout << "x2= " << scientific << setprecision(31) << x2 <<
"¥n";
}

3

複素数版プログラム例 (QD, C)...
#include <qd/dd_real.h>
マクロ関数の定義#include <qd/c_dd.h>#define TO_DOUBLE_PTR(a, ptr) ptr[0] = a.x[0]; ptr[1] = a.x[1];/* (a_r + b_r) + (a_i + b_i) I */void c_zdd_add(const double *a_r, const double *a_i, const double *b_r,
const double *b_i, double *c_r, double *c_i) {dd_real cr, ci;
cr = dd_real(a_r) + dd_real(b_r);
倍々精度複素数のci = dd_real(a_i) + dd_real(b_i);
四則演算関数は自作
TO_DOUBLE_PTR(cr, c_r);
TO_DOUBLE_PTR(ci, c_i);}...
変数1個につき, 実部と虚部それぞれ大きさ2個の倍精度配列
int main(void){
double a_r[2], a_i[2], b_r[2], b_i[2], c_r[2], c_i[2], x2_r[2], x2_i[2];...
c_dd_copy_d(1.01, a_r);
四則演算は関数呼び出し...
c_zdd_mul(b_r, b_i, b_r, b_i, t1_r, t1_i); /* b*b */c_zdd_mul(a_r, a_i, c_r, c_i, t2_r, t2_i); /* a*c */...
c_dd_write(x2_r);

4

複素数版プログラム例 (Fujitsu fastdd, C++)
#include <iostream>#include <iomanip>#include <stdlib.h>#include <complex>#include <fast_dd.h>using namespace std;
int main(void)
疑似4倍精度複素数の宣言{
complex<dd_real> a, b, c, x1, x2;
long double lx;
a = 1.01L;
b = 2718281.0L;
c = 0.01L;

4倍精度定数の作り方が用意されていないので，
long double定数を代入

x2 = (-b+sqrt(b*b-4.0*a*c))/(2.0*a);lx = to_long_double(real(x2));
cout << "x2= " << scientific << setprecision(31) << lx <<"¥n";
4倍精度の出力機能は用意されていないのでlong}

doubleとかに変換してから出力する必要あり

5

複素数版プログラム例 (ARPREC, C++)
#include <iostream>#include <iomanip>#include <stdlib.h>#include <arprec/mp_complex.h>using namespace std;
int nr_digits = 32;

精度の設定

int main(void){
mp::mp_init(nr_digits);
任意精度変数の宣言
mp_complex x2;
mp_complex a("1.01E+0","0E+0");
任意精度定数の代入
mp_complex b("2718281.0E+0","0E+0");
mp_complex c("0.01E+0","0E+0");x2 = (-b+sqrt(b*b-4*a*c))/(2*a);cout << "x2= " << scientific << setprecision(31) << x2.real << "" << x2.imag << "¥n";
mp::mp_finalize();}6

複素数版プログラム例 (exflib, C++)#define PRECISION 32#include <iostream>#include <iomanip>#include <complex>#include <stdlib.h>#include <exfloat.h>

精度の設定

using namespace std;
int main(void){
complex<exfloat> a, b, c, x2;
a = "1.01";
b = "2718281.0";
c = "0.01";

任意精度複素数の宣言

任意精度定数の代入

x2 = (-b+sqrt(b*bcomplex<exfloat>(4,0)*a*c))/(complex<exfloat>(2,0)*a);cout << "x2= " << scientific << setprecision(31) << x2 <<
"¥n";
}

7

複素数版プログラム例 (MPC, C++)
#include <iostream>#include <iomanip>#include <stdlib.h>#include <gmp.h>#include <mpfr.h>#include <mpc.h>

mpc_set_str(a, "(1.01 0.0)", 10, MPC_RNDNN);
mpc_set_str(b, "(2718281.0 0.0)", 10, MPC_RNDNN);
mpc_set_str(c, "(0.01 0.0)", 10, MPC_RNDNN);
mpc_set_str(_four, "(4.0 0.0)", 10, MPC_RNDNN);
mpc_set_str(_two, "(2.0 0.0)", 10, MPC_RNDNN);

精度の設定

mpc_mul(t1, b, b, MPC_RNDNN);
mpc_mul(t2, a, c, MPC_RNDNN);
mpc_mul(t3, t2, _four, MPC_RNDNN);
mpc_sub(t2, t1, t3, MPC_RNDNN);
mpc_sqrt(t1, t2, MPC_RNDNN);
mpc_sub(t2, t1, b, MPC_RNDNN);
mpc_mul(t3, a, _two, MPC_RNDNN);
mpc_div(x2, t2, t3, MPC_RNDNN);

乗算#define PRECISION 106
using namespace std;
int main(void){
mpc_t a, b, c, _four, _two, x2, t1, t2, t3;

任意精度変数の宣言

mpc_init2(a,PRECISION);
mpc_init2(b,PRECISION);
mpc_init2(c,PRECISION);
mpc_init2(x2,PRECISION);
mpc_init2(t1,PRECISION);
mpc_init2(t2,PRECISION);
mpc_init2(t3,PRECISION);
mpc_init2(_four,PRECISION);
mpc_init2(_two,PRECISION);

変数の初期化や任意精度定数の代入

cout << "x2= ";mpc_out_str(stdout, 10, 0, x2, MPC_RNDNN);
cout << "¥n";mpc_clear(a);mpc_clear(b);mpc_clear(c);mpc_clear(x2);mpc_clear(t1);mpc_clear(t2);mpc_clear(t3);mpc_clear(_four);
mpc_clear(_two);

オブジェクトの解放}

8計算精度に関する数値実験結果• 最適化オプション “-O3”, MPFRの丸めモードは最近接偶数丸め• 任意精度ライブラリは4倍精度相当で計算
言語

精度

結果

doubleC倍-3.688406236100904979616609403840e-09

long doubleC

拡張倍-3.678838531936214554367889299446e-09

long double@FX

C

128ビット-3.678795532912165323960311671534e-09

real(kind=16)

F90

4倍

-3.678795532912165323960311671534e-09

libquadmath

C

4倍-3.6787955329121653237603627688244e-09

_Quad

C++

4倍

-3.6884062361009050123322488651663e-09

QD

C++

疑似4倍 -3.678795532912165322960567157986e-09

fast_ddC++

疑似4倍 -3.678795532912165394583236142940e-09

ARPRECC++

32桁-3.6787955329121653239204992803347e-09

exflibC++

32桁-3.678795532912165323920499457919e-09

MPFR

C

106ビット-3.6787955329121653309585232663659e-09
9計算時間[s]に関する実験結果 (実数)言語

PC (gcc)

FX100 (fcc)

FX100 (gcc)

CX400 (icc)

doubleC

0.01050.01160.00546

0.0060long double

C

0.01100.78820.781

0.0095real(kind=16)F

0.10500.00000.781

0.0746libquadmath

C0.0921_Quad

C++

QD (C++)

C++

0.03040.07320.03330.0798QD (C)

C

0.04090.28840.02210.1364

fast_ddC++

0.1441

0.144

ARPRECC++

0.89686.6298

7.493

exflibC++

0.0856

MPFR

C

0.3329

N/A

N/A

N/A

N/A

N/A

N/A

N/A
1.3545

N/A
0.597

N/A
0.0714N/A
0.69340.10600.3939赤字は四倍精度相当の計算結果になっているもので最速
青字は計算結果の傾向が異なるもの
10計算時間[s]に関する実験結果 (複素数)言語

PC / gcc

FX100 / fcc

FX100 / gcc

CX400 / icc

doubleC

0.02270.02340.02280.0087long double

C

0.04571.46691.288

0.0155complex(kind=16)F

0.32500.00001.453

0.2133libquadmath

C0.2141

N/A

N/A

N/A

_Quad

C++

N/A

N/A

N/A

0.2004

QD (C++)

C++

0.13490.32320.342

0.1552QD (C)

C

0.13340.28620.241

0.1472

fast_ddC++

N/A

0.74640.803

N/A

ARPRECC++

1.556513.6597

5.515

1.3135

exflibC++

0.2591

N/A

MPC

C

0.19280.5868N/A
0.422

0.28640.2342

赤字は四倍精度相当の計算結果になっているもので最速
青字は計算結果の傾向が異なるもの
11計算時間[s]に関する実験結果2 (実数)言語

PC / gcc

FX100 / fcc

FX100 /
fcc(32reg)

FX100 / gcc

doubleC

0.00679

0.005480.00546

0.00546

long double

C

0.00645

0.788

0.785

0.781

real(kind=16)F900.0730

0.313

0.324

0.781

libquadmath

C

0.0636

QD (C++)

C++

0.02850.03330.03220.0692QD (C)

C

0.03230.02210.217

0.0851

fast_dd

C++

N/A

0.144

0.143

0.155

ARPREC

C++

0.408

7.493

7.247

2.90

exflibC++

0.0727

MPFR

C

0.166

N/A

N/A
0.597

N/A

N/A
0.574

N/A

N/A
0.469

赤字は四倍精度相当の計算結果になっているもので最速
12
青字は計算結果の傾向が異なるもの

計算時間[s]に関する実験結果2 (複素数)言語

PC / gcc

FX100 / fcc

FX100 /
fcc(32reg)

FX100 / gcc

doubleC

0.005420.02050.02050.0228long double

C

0.0248

1.455

1.452

1.288

complex(kind=16)

F90

0.186

0.787

0.584

1.453

libquadmath

C

0.138

QD (C++)

C++

0.0579

0.173

0.128

0.342

QD (C)

C

0.0569

0.112

0.113

0.241

fast_dd

C++

N/A

0.743

0.733

0.803

ARPRECC++

0.704

13.54913.801

5.515

exflib

C++

0.205

N/A

N/A

MPC

C

0.132

N/A

0.528

N/A

0.487

N/A

N/A
0.422

赤字は四倍精度相当の計算結果になっているもので最速
13
青字は計算結果の傾向が異なるもの

FX100における疑似四倍精度実行列ベクト
ル積演算(QD(C++))の最適化メッセージ例
38
39
40
41
42
43
44
...

i

i

2
2
2

for (i = 0; i < n; i++) {q[i] = 0.0;
#pragma loop noaliasfor (j = 0; j < n; j++) {q[i] += A[i*n+j] * p[j];}}・SIMD化されず
・SWPLされず・ループ展開2倍

jwd6101s-i "logisticmap_real_qd_quad.cpp", line 41: SIMD conversion is notapplied because a statement that prevents SIMD conversion exists.
jwd8670o-i "logisticmap_real_qd_quad.cpp", line 41: This loop is not
software pipelined because the loop contains a branch instruction which is
not for loop iteration.
jwd8202o-i "logisticmap_real_qd_quad.cpp", line 41: Loop unrolling
expanding 2 times is applied to this loop.
jwd8101o-i "logisticmap_real_qd_quad.cpp", line 42: Inline expansion isapplied to the user defined function '_ZmlRK7dd_realS1_'.
jwd8101o-i "logisticmap_real_qd_quad.cpp",line 42: Inline expansion is
インライン展開のメッセージが大量にでる
applied to the user defined function '_ZN7dd_realpLERKS_'.

14

FX100における疑似四倍精度複素行列ベク
トル積演算(QD(C++))の最適化メッセージ例
40
41
42
43
44
45
46
...

i

i

for (i = 0; i < n; i++) {q[i] = c_zero;#pragma loop noaliasfor (j = 0; j < n; j++) {q[i] += A[i*n+j] * p[j];}}・SIMD化されず
・SWPLされず・ループ展開なし

jwd6101s-i "logisticmap_complex_qd_quad.cpp", line 43: SIMD conversion is
not applied because a statement that prevents SIMD conversion exists.
jwd8670o-i "logisticmap_complex_qd_quad.cpp", line 43: This loop is not
software pipelined because the loop contains a branch instruction which is
not for loop iteration.
jwd8101o-i "/opt/FJSVmxlang/bin/../include/c++/std/stl/_string_io.c", line
44: Inline expansion is applied to the user defined function'_ZNKSt8ios_base5flagsEv'.
jwd8101o-i "logisticmap_complex_qd_quad.cpp", line 44: Inline expansion isapplied to the user defined function '_ZStmlI7dd_realESt7complexIT_ERKS3_S5_'.
jwd8101o-i "logisticmap_complex_qd_quad.cpp", line 44: Inline expansion isapplied to the user defined function '_ZNSt7complexI7dd_realEpLERKS1_'.

15

Numerical experimentson the convergence of the iterative method• Solving complex symmetric linear systems
𝐴𝐱 = 𝐟• coefficient matrices (from the SuiteSparse Matrix Collection)• dwg961b: n= 961, num. of nonzero= 19,591• qc2534: n= 2,534, num. of nonzero= 463,360
an electromagnetic field problem and a complex symmetricmatrix• 𝐟 = 𝐴 × 1 + 𝑖, 1 + 𝑖, … , 1 + 𝑖 ்• Iterative methods• COCG (Conjugate Orthogonal Conjugate Gradient) method w/
incomplete Cholesky preconditioning• Multiple-precision libraries• libquadmath (GCC 5.4.0), QD 2.3.17, ARPREC 2.2.18• Computers• Intel Core i7-6500U / Ubuntu 16.04 TLS / gcc-5.4.0
16

Comparison of number of iterations withdifferent precision arithmetic
Convergence criterion: 10-7

double100digits

100digits
1000digits

doublepseudo-quad
quadruple

Convergence history for dwg961b

pseudo-quad
quadruple

Convergence history for qc2534• “100digits” means 100 decimal digits of precision calculated by ARPREC asreference.
• In solving complex symmetric linear systems,• multiple-precision got faster convergence compared with double• (pseudo-)quadruple is enough precision
17

Comparison of calculation time withdifferent precision arithmetic
Convergence criterion: 10-7
Matrixdwg961bqc2534

Precision name

Decimaldigitsnon-usedouble

15

1,058

0.480

QD

pseudo-quadruple

32

479

1.216

quadmath

quadruple

34

457

4.132

ARPREC

pseudo-quadruple

32

383

32.503

non-usedouble

15

572

9.437

QD

pseudo-quadruple

32

248

20.841

quadmath

quadruple

34

249

73.272

ARPREC

pseudo-quadruple

32

240

748.447

Library# of Iter. Time [s]• QD is faster than GCC’s libquadmath and ARPREC in total.
• QD takes at most 3 times longer than double precision in small-scale cases.
18

TEAM Workshop Problem 29• Benchmark model for hyperthermia simulation• Problem size: 134,573• Frequency: 300 (MHz)• Relative permittivity of Phantom: 80.0
• Electrical conductivity of Phantom: 0.52 (S/m)• Linear solver• Iterative method: COCG w/ Symmetric SOR (𝜔 = 1.05) preconditioning• Convergence tolerance: relative residual norm < 10-7

• Computers• Intel Core i7-6500U / Ubuntu 16.04 TLS / gcc-5.4.0
sensing coil

Phantomfeeding coil

Model

Finite element mesh

19Computational performance of 300 (MHz)doublequadruple

mixed

Precision# of iter.

Time [hr.]double24,5105

5.78

pseudo-quadruple (QD)

43,598

5.69

quadruple (libquadmath)

38,64314.09
20

4倍精度演算 IEEE 754-2008のbinary128形式
 ソフトでエミュレーションするため低速 double-double形式
 倍精度を2つ並べた形式• 倍精度演算器を使って計算できるためIEEE形式よりも高速 加減乗算の演算方法は広く知られている
 doube-double形式の実装例• フリーソフト : QDなど
• 富士通 : 高速4倍精度基本演算ライブラリ(以降fast_dd)• 他ベンダー• IBMでは4倍精度としてdouble-double形式をサポートするシステム有り• NEC : ASLQUAD• この他に、ユーザ自身が実装する場合もあり
21
FUJITSU CONFIDENTIAL

21

Copyright 2019 FUJITSU LIMITEDQD と fast_dd
 QDの実装 C++版はヘッダファイル内に演算が書かれており、インライン展開される• 翻訳時に最適化が適用されることで高速化を実現

 Fortran版は、wrapperのみでCの関数を呼び出している
 注意点• C++版は、インライン展開されないと高速化されず、また、演算順序を変更する最適化が適用されると結果が正しくなくなる• Fortranは1要素ずつ関数を呼び出すため高速化されない fast_ddの実装
 ユーザの指定する最適化オプションによって性能や結果に影響が出ることを避けるため、インライン展開方式は採用しない• 演算は基本的に1要素ずつ関数を呼び出すことで実現

 高速化のために、マルチ関数（2要素ずつ計算）・ベクトル関数（入出力が1次元配列）を用意マルチ関数・ベクトル関数を使用して頂くことを期待• 7/20ご報告の評価で使われたロジスティック写像は漸化式でベクトル関数が使えな
いためQDのC++版（インライン）が高速
22
FUJITSU CONFIDENTIAL

22

Copyright 2019 FUJITSU LIMITED

実行性能
 加算の性能
FX100 1.5Ghz
N=1024 (L1に乗るサイズ)

use fast̲dd
INTEGER::N
TYPE(dd̲real)::X(N),Y(N),Z(N)DO I=1,N
Z(I)=X(I)+Y(I)END DO※ マルチ関数、ベクトル関数は関数呼出し化した版を使用

性能 mflops(注1)QD
fast_dd 単体関数
fast_dd マルチ関数
fast_ddベクトル関数

備考

C++

68.2

インライン、SWPLのみ

Fortran

21.9

1要素ずつ関数呼出し

C++

13.1

1要素ずつ関数呼出し

Fortran

22.2

1要素ずつ関数呼出し

C++

43.0

2要素ずつ関数呼出し

Fortran114.0

2要素ずつ関数呼出し

C++

805.0

SIMD+SWPL

Fortran874.0

SIMD+SWPL

注1) double-double形式の加算1回を1flopとしたときの性能
23
FUJITSU CONFIDENTIAL

23

Copyright 2019 FUJITSU LIMITED