2.2 3 次元非圧縮性熱流体計算プログラム COSMOS の測定評価（株）豊田中央研究所

2.2.1

堀之内 成明はじめに自動車の設計・開発における熱流体 CAE (Computer-Aided Engineering) は，日常的なものとして様々な製品・部品を対象に用いられてきているが，実際には，非定常な乱流に加えて，熱，騒音，運動（移動境界），反応などが連成するような現象も多く，これらに対して信頼性の高い結果を得るためには，計算格子幅や，時間ステップ間隔を十分に小さくする必要が出てくるため，その結果として必然的に計算規模が大きくなる．
このような大規模熱流体シミュレーションは，従来ではベクトル演算方式の計算機で行われることが多く，プログラムもこの上で高速計算が可能となるようにチューニング，ならびにアルゴリズム開発がなされてきた．本稿では，このようなプログラムを，ベクトル計算機での性能も維持しつつ，「京」
をはじめとするマルチコアのスカラー型計算機においても高い実効性能を満たすことを目的として，高速化に向けた改良を試みる．

2.2.2

プログラム概要

ここでは，（株）豊田中央研究所で開発した非定常非圧縮性乱流計算プログラム COSMOS を対象にチューニングを行う．このプログラムは，空力や空力騒音といった問題に適用するため，自動車まわりの高 Reynolds 数非定常流れを精度よく計算することを目的に開発された [1]．

2.2.2.1 プログラムの特徴
本プログラムは，以下の特徴を有する．

• Large Eddy Simulation (LES) による非定常乱流の計算に関して，以下の独自技術を開発し，実問題に適用し得る計算精度と実用性を両立させている．

– Grid-scale (GS) 成分計算の離散化スキームとして，質量，運動量のみならず，運動エネル
ギーまでをも高い精度で保存させ，数値粘性を用いなくても安定な計算を可能とする「改良型コロケーション格子法 [2]」– Subgrid-scale (SGS) 成分の乱流モデルとして，様々な問題に対してモデル係数のチューニングが不要で，かつダンピング関数を必要としない「混合時間スケール SGS モデル [3]」• 計算格子には，計算精度に優れているというメリットから「境界適合型構造格子」を用いる．これは，碁盤目状の格子を，対象となる物体の表面に沿わせ，且つ表面近くに格子点を集中させたものである．また，三角形（2 次元の場合）やテトラ（3 次元の場合）を用いる非構造格子に比べて，プログラミングに際して間接参照（リストアクセス）の必要性が軽減でき，計算効率の観点からも有利である．但し，車両などの複雑形状に対しては，格子生成に多くの工数やノウ
ハウを必要とするというデメリットもあり，そのような計算領域全体を一つの格子で覆うこと自体が現実的には困難であったり，計算格子の歪みが大きくなり質が低下するといった課題もある．これを解決するために，局所的な形状に適合させた高品質な部分的構造格子（以下，部分格子）をパッチワーク的に重ね合わせる「重合格子法」を採用している．車両形状に対する重合格子の例を図 1 に示す．

• 非圧縮性流れの基礎方程式は，連続の式と Navier-Stokes 方程式（運動方程式）であるが，SMAC
法に準じて，これらから圧力に関する Poisson 方程式を導いて解く．時間方向の離散化には，Crank-Nicolson 法に基づく陰解法を採用しており，これから速度場に関する連立一次方程式が導かれる．一方，圧力の Poisson 方程式も，連立一次方程式を解くことに帰着される．本プログラムの計算時間の 8 割程度は，これらの連立一次方程式を解くことに費やされるため，これらをいかに効率的に計算できるようにするかがキーとなる．この時，上述した重合格子法との兼ね合いにより，単一格子での行列方程式計算よりも，アルゴリズムや配列の取り方などが，やや複雑になり，この点も含めたチューニングが求められる．

図 1 自動車まわりの重合格子の例

2.2.2.2 連立一次方程式の数値解法
重合格子法の計算アルゴリズムは，領域分割法として提案されている Schwarz 法の中でも，乗法的Schwarz 法 (Multiplicative Schwarz algorithm) の考え方に等しい（付録 p1 参照）．その反復計算自体も数値計算上の研究課題となり得るが，今回は，各部分格子の中で，離散化された連立一次方程式を，
係数行列がスパースな場合に用いられる手法の一つである SOR (Successive Over Relaxation) 法で計算する際の効率化を対象とする．
連立一次方程式の係数行列を A，解ベクトルを x，右辺ベクトルを b とする．すなわち，Ax = b である．また，反復法の第 k ステップ目での近似解ベクトルを x(k) と表すこととする．SOR 法のベースとなる Gauss-Seidel 法は，元の係数行列 A を，A = L + D + U (L: A の対角要素を除く下三角行列，D:A の対角行列，U: A の対角要素を除く上三角行列) と分解した (L + D + U)x = b から，(L + D) x = −U x + b(1)(L + D) x(k+1) = −U x(k) + bx(k+1) = − (L + D)−1 U x(k) + (L + D)−1 b

(2)
(3)

とし，これに基づき，として近似解ベクトルを更新（緩和）する手法である．式 (2) から，近似解ベクトルを要素ごとに逐次計算していく際に，下三角部分 L にかかる要素は，常に最新の値が用いられていることがわかる．このことから，プログラミング上は回帰参照が発生することになる．なお，実際のプログラムでは，式(3) に基づいて x(k+1) を計算するのではなく，
r(k) ≡ b − Lx(k+1) − (D + U) x(k)
とおいた残差ベクトル r(k) を用いて，

x(k+1) = x(k) + D−1 r(k)(4)

のように計算する．さらに，この式 (4) の右辺で，前の近似解ベクトル x(k) に第 2 項を足し込むときに，加速係数（緩和係数とも言う）ω を掛け，x(k+1) = x(k) + ωD−1 r(k)(5)

とするのが SOR 法である．ω は 0 < ω < 2 の範囲の値をとり，ω = 1 の場合が Gauss-Seidel 法ということになる．

Navier-Stokes 方程式，ならびに圧力の Poisson 方程式を，それぞれ 2 次中心差分に相当するスキー
ムで空間的に離散化することで，連立一次方程式を得るわけだが，図 1 のように格子線の曲りや不等間隔も許すため，計算空間（格子線座標 i, j, k の空間）で考えると，計算点から見て，斜め方向の点も参照することとなる．Navier-Stokes 方程式に関しては，対流項の 3 次風上差分（QUICK 法）の可能性を考慮し，i, j, k の各方向に 5 点差分と，斜めの点を合わせて，全部で 27 点の差分，圧力 Poisson 方
程式に関しては全部で 19 点の差分となる．図 2 にその関係を示す．そして，この参照関係から，係数行列の構造は図 3 のようになる．

図 2 格子点の参照関係

図 3 係数行列の形（イメージ）

さて，上記のような依存関係を含む手続きを正しく計算するために，ベクトル計算機向けのプログ
ラミングと同様に，マルチコアで並列計算においても回帰演算を回避する必要があり，本プログラム
ではマルチカラーオーダリングで実現している．オリジナル版では，ベクトル機で用いたものと同じ

1 次元化された配列を最小の色数となる 7 色でオーダリングする．圧力 Poisson 方程式に対するソー
スプログラムを付録 p2 に示す．

2.2.3

測定環境

今回は，円柱まわりの非圧縮性非定常流れの計算を性能評価用の例題とし，付録 p3 にある断面を示
すように，円柱のまわりのみと計算領域全体を覆う２種類の構造格子を重合させた計算格子を用いる．

表 1 利用計算機の概要
機種CPU スペック

CPU 数 / ノード

コア数 / CPU

4

8

RX

Intel Xeon X7560(PRIMERGY RX600S5)(Nehalem-EX), 2.26GHzFX1SPARC64 VII, 2.5GHz

1

4

京

SPARC64 VIII, 2.0GHz

1

8

FX10

SPARC64 IX, 1.848GHz

116

測定に利用した計算機の一覧を表 1 に示す．
また，それぞれの計算機で測定する際のコンパイルオプションは次の通りである．

• RX: -Am -Kfast,parallel,ocl,vppocl,tl trt,prefetch -X9 -w -KXEON55 -Qt -Ectfiu -xdir=inline dir-Karray subscript,array subscript rank=2,array subscript elementlast=25-Kprefetch iteration=15,preex• FX1: -Am -Kimpact,ocl,vppocl,tl trt -X9 -w -Qt -Ectfiu -xdir=$(MKDIR)/inline dir-Kprefetch model=FX1-Karray subscript,array subscript rank=2,array subscript elementlast=25-Kprefetch iteration=15,preex• 京,FX10: -Kfast,parallel,ocl,vppocl -Ntl trt -Nrt tune -X9-Qt -xdir=inline dir -W0,-zmpb=Swpl:message=detail-Karray subscript,array subscript rank=2,array subscript elementlast=25-Kprefetch iteration=15,preex
但し，FX10 については，後述するチューニング実行時 (2.2.4.2 節) には，これに加えてオプション(-Kprefetch stride,prefetch sequential=soft) が追加されている．
2.2.4

測定結果

2.2.4.1 オーダリングの違いに関する性能評価
オリジナル版のままでは，計算点に対するストライドアクセスが発生するために，キャッシュ利用
効率が低下することが予想された．そこで，計算空間座標 (i, j, k) に対して，各 2 色の組み合わせとする計 8 色でのオーダリングを行い（以下，改良版１），オリジナル版と比較した（但し，圧力 Poisson方程式から得られる行列方程式のみ）．後者のソースプログラムを付録 p4 に，FX1 の 1 ノード 4Coreでの比較結果，PA 情報を p5,p6 に示す．これによると，予想通り，改良版１の方がキャッシュアクセ
ス待ちを 7 割程度減らすことができ，プログラム全体で 8% 程度の実効速度向上につながった．
しかしながら，改良版１であっても，実効性能のピーク性能比は 4% 強にしか達しておらず，さらなる性能向上が求められるため，次に，付録 p7 に示す改良を試みた（以下，改良版２）．これは，計算点から遠い参照点を別格子とし，すなわち k 方向に交互にグルーピングする扱いとすることで，さらにキャッシュ利用効率を上げることを狙ったものである．ソースプログラムを付録 p8 に，社内計算機での評価結果を付録 p9 に示す．コア数が少ない場合は，若干の効果は見られたが，残念ながら期待したほどの結果とはならなかった．また，FX1，京，FX10 でのプログラム全体の計算時間計測結果を付録 p10 に示す．

2.2.4.2 コンパイルオプションによるチューニングの結果
前節でより良い性能を示した改良版においても，プロファイラの結果を見ていくと，浮動小数点演
算待ちや，ロード・ストア待ちが多く発生していることがわかる（付録 p11）
．これらの改良のために，
以下のチューニングを行った．

1. 圧力 Poisson 方程式から得られる連立一次方程式の求解ルーチン sor8p b h（改良版２）に対して，ソフトウェアパイプラインによる命令スケジューリング（具体的な詳細は付録 p12〜p13）

2. Navier-Stokes 方程式から得られる連立一次方程式の求解ルーチン sor7v b（オリジナル版）に対して，配列 clp に対するソフトウェアプリフェッチと，ストライドアクセスオプションの指定（具体的な詳細は付録 p16〜p17）
これらのチューニングの結果，付録 p14〜p15, p18 に示すように，それぞれの連立一次方程式の SOR法による求解部分において，10 ％〜11 ％の性能向上（該当サブルーチンのみでの実行時間の短縮）を得ることができた．

2.2.5

まとめ

重合格子を用いた非定常乱流計算プログラムにおける連立一次方程式の反復ソルバーを対象に，マルチコア計算機に対する性能評価，ならびに計算効率化のためのチューニングを行った．しかしながら，ベクトル計算機との併用を前提にした配列構造やオーダリング（演算順序）の範囲内での性能向上には限界があり，必ずしも，マルチコア本来の並列計算性能を十分に引き出せてはいないと考えられる．また，SOR 法を用いたこの反復ソルバーのアルゴリズムをそのまま書き下すと，一つの計算点に対して周辺の 19 点のデータを参照しながらの反復計算が必要となり，このままでは，どうしても

Byte/Flops（演算数に対するデータのロード・ストア数の割合）が大きくなってしまうことも，性能向上が容易でない要因の一つと考えられる．
以上，今回の性能評価とチューニングから得られた知見から，さらに性能を向上させるためには，• 反復解法レベルでの計算アルゴリズムの見直し• マルチコア計算機に特化した配列構造への書き換えなどの，より深い取組みが必要になってくると考えられる．今後の課題としたい．
参考文献[1] 堀之内成明, 稲垣昌英, 加藤由博, 自動車空力騒音のシミュレーション技術, 応用数理, 15 (2005),
60–64.
[2] 稲垣昌英, 安倍賢一, コロケーション格子を用いた LES の計算精度改善に関する一考察, 日本機械学会論文集, B 編, 64, 623 (1998), 1981–1988.

[3] 稲垣昌英, 近藤継男, 長野靖尚, 実用的な LES のための混合時間スケール SGS モデル, 日本機械学会論文集, B 編, 68, 673 (2002), 2572–2579.

3次元非圧縮性熱流体計算プログラムCOSMOSの測定評価 （補足資料）
株式会社 豊田中央研究所
堀之内 成明1

連立一次方程式の求解について

計算時間の大半（6割～8割）は連立一次方程式を解くことに費やされる．
次元数（=格子点数）: 数10万～数100万z 規則非対称帯行列z 各行の非零要素数: ～ 25z

解法: 定常反復法 (SOR, Gauss-Seidel)，クリロフ部分空間法 (Bi-CGSTAB, ...)z

重合格子の場合は，このような行列を格子ブロックごとに解くが， 反復の途中で境界領域の格子点に対する補間を行う．

(1),(2)の変化量が小さくなるまで
繰り返し
Gauss-Seidel like

格子 m=1m=2(1)格子
補間点をディリクレ境界とし，
r ... r格子m=(1)
A補間点をディリクレ境界とし，xΩ m Ωrm = bΩrm を解く
格子
m=NVOL(1)
A補間点をディリクレ境界とし，Ω m xΩrm = bΩrm を解く(2) 他の格子の最新の値から補間(1)rA
補間点をディリクレ境界とし，Ω xΩ =rbΩr を解く(2) 他の格子の最新の値から補間xΓAm m ←xrmΙ x=Ω ′bm を解くrΩ Ω r Ω
常に最新の値(2) 他の格子の最新の値から補間xΓ mm ← mΙ xΩ ′ mrr(2) 他の格子の内部の値から補間x ←Ι x[ ][ ][ Ωr′ ]
ΓrxΓ ← Ι [xΩ ′ ]mm

他の格子の内部領域

Multiplicative（乗法的） Schwarz Method

2

7色オーダリングによるノード内並列化（オリジナル版）
 SOR法を採用（実際には緩和係数を1としているので，Gauss-Seidel法）
 各格子ブロックごとにマルチカラーオーダリングによる自動並列化で対応
z ベクトル機でのオーダリングをそのまま流用 ⇒ 7色飛び（1次元的）配列のとり方も，それに合せて1次元並びになっている．
real*8 :: p0(n1max), bp(n1max)
real*8 :: clp(n1max,19)

real*8 :: p0(n1max), bp(n1max)
real*8 :: fj_clp(19, n1max/7+1, 0:6)

do m=1,nvol
do ic=0,6
do n1=n1min(m)+ic,n1max(m),7res=p0(n1-iip1-1)*clp(n1,1)$+p0(n1-ijp1-1)*clp(n1,2)$+p0(n11)*clp(n1,3)
$
+ ...
$
+p0(n1-ijp1)*clp(n1,9)$+p0(n1)$+p0(n1+ijp1)*clp(n1,11)
$
+ ...
$
+p0(n1+iip1+1)*clp(n1,19)&-bp(n1)p0(n1)=p0(n1)-res*omegap
enddo
enddo
オリジナルのベクトル機用コード
enddo

do m=1,nvol
連続アクセスとした
do ic=0,6
i=1
do n1=n1rs(m)+ic,n1re(m),7res=p0(n1-iip1-1)*fj_clp(1,i,ic)$+p0(n1-ijp1-1)*fj_clp(2,i,ic)$+p0(n1- 1)*fj_clp(3,i,ic)
$
+ ...
$
+p0(n1-ijp1)*fj_clp(9,i,ic)$+p0(n1)$+p0(n1+ijp1)*fj_clp(11,i,ic)
$
+ ...
$
+p0(n1+iip1+1)*fj_clp(19,i,ic)&-bp(n1)i=i+1
p0(n1)=p0(n1)-res*omegap
enddo
enddo
少し修正したコード（これで計測）
enddo

ｃｌｐのみ色ごとに3

性能評価計算の概要
 LESを用いた円柱まわりの非定常流れ場の計算
z （実用計算を踏まえ）重合格子を使用

円柱まわり:
161x56x81=757,188点

背景:121x61x81=620,248点

4

8色オーダリング（改良版１）の概要
 I,j,k各方向に2色，計8色のオーダリングに変更．但し，配列の構造（宣言）は基本的に変えていない．
do m=1,nvol

real*8 :: p0(n1max), bp(n1max)
real*8 :: fj_clp0(19, m2max/7+1, nvol),&fj_clp1(19, m2max/7+1, nvol),&...
&fj_clp7(19, m2max/7+1, nvol)

色ごとのループ

do i=1,ico0max(m)kk=(i-1)/ijp8aa
ij=-kk*ijp8aa+(i-1)
jj=ij/iip8aii=-jj*iip8a+ij
n1=(kk*2)*ijp1+(jj*2)*iip1+(ii*2)+n1rstr
res=p0(n1-iip1-1)*fj_clp0(1,i,m)+p0(n1-ijp1-1)*fj_clp0(2,i,m)+p0(n11)*fj_clp0(3,i,m)+ ...
+p0(n1-ijp1)*fj_clp0(9,i,m)+p0(n1)+p0(n1+ijp1)*fj_clp0(11,i,m)+ ...
+p0(n1+iip1+1)*fj_clp0(19,i,m)-bp(n1)p0(n1)=p0(n1)-res*omegap
enddo

$
$
$
$
$
$
$
$
&

do i=1,ico1max(m)
...
enddo
...

（圧力Poissonのみに適用）

enddo

5

評価結果の概要


JAXA殿FX1（1ノード4Core）での計測結果
スレッド並列：自動並列（VISIMPACT)z ｺﾝﾊﾟｲﾙｵﾌﾟｼｮﾝ: -Am -Kimpact,ocl,vppocl,tl_trt -X9 -w -Qt -Ectfiu -xdir=inline_dirz

7色オーダリング
ルーチン名
sor7p_b_
sor7v_b_sgs_
corvcn_sor7p_
全体

プロファイラ(count数)
31,883
10,455
7,608
4,036
2,995
64,9668色オーダリング

CPU時間
[sec]
[%]
318.8349.1%
104.55
16.1%
76.08
11.7%
40.36
6.2% 4174.03 Mips
29.95
4.6%
649.66100.0% 1608.46 Mflops

ルーチン名
sor8p_b_
sor7v_b_sgs_
corvcn_sgsbnd_全体

プロファイラ(count数)
28,331
10,482
7,618
4,018
1,153
59,442

csms_sor7(全体） - 時間 (sec)整数ロー ドメモリアクセス待ち
整数ロー ドキ ャッシュアクセス待ち
分岐命令待ち1命令コミット
4命令コミット

浮動小数点ロードメモリアクセス待ち
浮動小数点ロードキ ャッシュアクセス待ち
命令フェッチ待ち
整数レ ジスタ書き込み制約

ストア待ち
整数演算待ちuOPコミット
浮動小数点レ ジスタ書き込み制約

CPU時間
[sec]
[%]
283.3147.7%
104.82
17.6%
76.18
12.8%
40.18
6.8%
11.53
1.9%
594.42100.0%

5516.60 Mips
1742.24 Mflops

csms_sor8（全体） - 時間 (sec)
I/Oアクセス待ち
浮動小数点演算待ち
その他の待ち2/3命令コミット(その他)

時間(sec)

整数ロー ドメモリアクセス待ち
整数ロー ドキャッシュアクセス待ち
分岐命令待ち1命令コミット
4命令コミット

浮動小数点ロードメモリアクセス待ち
浮動小数点ロードキ ャッシュアクセス待ち
命令フェッチ待ち
整数レ ジスタ書き込み制約

ストア待ち
整数演算待ちuOPコミット
浮動小数点レジスタ書き込み制約

I/Oアクセス待ち
浮動小数点演算待ち
その他の待ち2/3命令コミット(その他)

時間(sec)7.0E+02

FP演算待ち: 10%

7.0E+02

FP演算待ち: 12% (+9sec)
6.0E+02

6.0E+02

5.0E+02

5.0E+02

4.0E+02

4.0E+02

FPﾛｰﾄﾞｷｬｯｼｭｱｸｾｽ待ち: 23%
3.0E+02

3.0E+02

FPﾛｰﾄﾞｷｬｯｼｭｱｸｾｽ待ち: 13% (-74sec)

2.0E+02

2.0E+02

FPﾛｰﾄﾞﾒﾓﾘｱｸｾｽ待ち: 41%

1.0E+02

1.0E+02

0.0E+00

FPﾛｰﾄﾞﾒﾓﾘｱｸｾｽ待ち: 38% (-40sec)

0.0E+00
Process 0 Thread 0

Process 0 Thread 1

Process 0 Thread 2

Process 0 Thread 3

Process 0 Thread 0

Process 0 Thread 1

Process 0 Thread 2

Process 0 Thread 3

6

コスト上位ルーチン別のPA情報

圧力Poisson方程式のSOR計算
7色オーダリング

3.5E+02

8色オーダリング

時間(sec)3.0E+02

3.0E+02

2.5E+02

2.5E+02

2.0E+02

FPﾛｰﾄﾞｷｬｯｼｭｱｸｾｽ待ち: 34%

2.0E+02

1.5E+02
1.5E+02

1.0E+02
1.0E+02

FPﾛｰﾄﾞﾒﾓﾘｱｸｾｽ待ち: 46%

FPﾛｰﾄﾞﾒﾓﾘｱｸｾｽ待ち: 43%
5.0E+01

5.0E+01

0.0E+00

0.0E+00
Process 0 Thread 0

Process 0 Thread 1

Process 0 Thread 2

Process 0 Thread 0

Process 0 Thread 3

Process 0 Thread 1

Process 0 Thread 2

 運動方程式のSOR計算

 SGS項算出

1.2E+02

9.0E+01

Process 0 Thread 3

8.0E+01
1.0E+02
7.0E+01

8.0E+01

6.0E+01

5.0E+01
6.0E+01
4.0E+01

FP演算待ち: 33%
4.0E+01

3.0E+01

FPﾛｰﾄﾞﾒﾓﾘｱｸｾｽ待ち: 70%

2.0E+01
2.0E+01
1.0E+01

0.0E+00

0.0E+00
Process 0 Thread 0

Process 0 Thread 1

Process 0 Thread 2

Process 0 Thread 3

Process 0 Thread 0

Process 0 Thread 1

Process 0 Thread 2

Process 0 Thread 3

7

更なる性能向上に向けた試み（改良版２）
今回
前回

6
8

全てを
一つの配列
p0に持つ

4

5

2

4

1

3
8

7

6

8

8

7

6
8

8

3

2

4

p1

8

p2

4
5

7

7

6

k方向を交互に二つの配列
p1, p2に持つ8を計算するときは同じ配列をみる
はを計算するときに参照される

数字は色分けの番号

キャッシュに乗りやすい近くの点は別の配列をみる
キャッシュに乗りにくい遠くの点

8

ソースのイメージ
real*8 :: p1(m2max,0:1,nvol)
real*8 :: p0(n1max), bp(n1max)
real*8 :: fj_clp0(19, m2max/7+1, nvol),&fj_clp1(19, m2max/7+1, nvol),&...
&fj_clp7(19, m2max/7+1, nvol)

do m=1,nvol
do i=1,ico0max(m)kk=(i-1)/ijp8aa
ij=-kk*ijp8aa+(i-1)
jj=ij/iip8aii=-jj*iip8a+ij
nn1=kk*ijp1+(jj*2)*iip1+(ii*2)+nn1rst
n1 =nn1+kk*ijp1+n1rsg
nn2=nn1-ijp1
$
$
$
$
$
$
$
$
$
$
$
&

ここでは p1,p2 をp1(:,0:1,:) として一つの配列にしてある

res=p1(nn1-iip1-1,0,m)*fj_clp0(1,i,m)+p1(nn2-1,1,m)*fj_clp0(2,i,m)+p1(nn11,0,m)*fj_clp0(3,i,m)+ …
+p1(nn1-iip1 ,1,m)*fj_clp0(8,i,m)+p1(nn2,1,m)*fj_clp0(9,i,m)+p1(nn1,0,m)+p1(nn1,1,m)*fj_clp0(11,i,m)+ …
+p1(nn1+1,0,m)*fj_clp0(17,i,m)+p1(nn1+1,1,m)*fj_clp0(18,i,m)+p1(nn1+iip1+1,0,m)*fj_clp0(19,i,m)-bp(n1)p1(nn1,0,m)=p1(nn1,0,m)-res*omegap

enddo
do i=1,ico1max(m)
...
enddo
...
enddo

9

結果

圧力Poisson部のみの計算時間[秒] 結果的にはあまり速度向上はしなかった．
120
100
80
前回

60

今回

40
20
0
2core

4core

8core

32core但し，以下の社内計算機で評価
PRIMERGY RX600S5 x 11node
Intel Xeon X7560 (Nehalem-EX)
8core/chip, 4chip/node
2.26GHz, 24MB Chache

測定結果
 CSMS_SOR8最新版（改良版２）のFX1、京、FX10での測定結果
 イタレーションは、step=50で測定■FX1、京、FX10 実行時間比較 (timeコマンドのreal)
FX1
京
実行時間
実行時間
スケール Thread
スケール
Thread(秒)(秒)
1th 931.36
1.00
1th 997.14
1.00
csms_
sor8
4th 405.91
2.31
8th 144.59
6.90

FX10
実行時間
スケール
Thread(秒)
1th 1187.14
1.00
16th
90.64
13.10

■FX1比
csms_
sor8

Thread
1th
4th

FX1
1.00
1.00

(目安：FX1とのF演算ピーク性能比)
Thread
FX1
1th
1.00
コア
4th
1.00
ノード

Thread
1th
8th

京
0.93
2.81

Thread
1th
16th

FX10
0.78
4.48

Thread
1th
8th

京
1.59
3.17

Thread
1th
16th

FX10
1.47
5.87

10

Copyright 2013 FUJITSU LIMITED

高コストのPA情報（FX10）
 PA情報スレッド０
step=50

メモリスループットは少し余裕がある。
L2ミス(dm)率が高い。
演算待ちが多い。

FX10 (16スレッド)

浮動小数点待ちが多く
命令スケジューリング
の余地があるか？

メモリスループットは
少し余裕がある。
L2ミス(dm)率が高い。

sor7v_b T0
sgs T0sor8p_b_h T0③②①
Performance

SIMD
実行時間 (sec)浮動小数点演算
ピーク比

11.11
12.71
36.80

2.94%
4.80%
3.55%

MFLOPS
434
710
525

458
1390
1321

SIMD演算命令率 SIMDロード命令率 SIMDストア命令率
SIMD命令率(/SIMD対象命令 (/SIMD対象演算 (/SIMD対象ロード (/SIMD対象ストア
命令数)
命令数)
命令数)
数）
25.21%
31.41%
99.99%
0.00%
0.00%
0.00%
0.00%
0.00%
0.00%
0.00%
10.65%
21.57%
99.98%0.00%
0.00%

SIMD命令率(/有効総命令数)

浮動小数点演算数

MIPS

4.82E+09
9.02E+09
1.93E+10

Cache

sor7v_b T0
sgs T0sor8p_b_h T0

L1I ミス率(/有効総命令数)

L1D ミス率(/ロード・ストア数)

0.05%
0.00%
0.02%

25.99%
4.12%
6.24%

L2 ミス率(/ロード・ストア数)
sor7v_b T0
sgs T0sor8p_b_h T0

9.62%
0.34%
5.46%

L2 ミス数
2.71E+08
2.20E+07
1.03E+09

ロード・ストア数
2.81E+09
6.45E+09
1.89E+10
L2 ミス dm 率(/L2 ミス数)
23.17%
27.26%
7.47%

L1D ミス数

L1D ミス dm率(/L1D ミス数)

L1D ミス hwpf
率(/L1D ミス数)

L1D ミス swpf 率(/L1D ミス数)

52.05%
70.90%
41.06%

26.33%
3.43%
58.93%

21.62%
25.67%
0.00%

7.32E+08
2.66E+08
1.18E+09
L2 ミス pf 率(/L2 ミス数)

メモリスループット(GB/sec)

L2 スループット(GB/sec)

3.61
0.24
3.92

8.43
2.68
4.09

76.83%
72.74%
92.53%11µDTLB ミス率
mDTLB ミス率(/ロード・ストア(/ロード・ストア数)
数)
0.00292%0.00002%0.00083%0.00001%
0.00424%0.00001%高コスト③sor8p_b_hチューニングについて
 Loop1チューニング
 「浮動小数点演算待ち」と「整数演算待ち」が多い。
ソフトウェアパイプライン(SWP)による命令スケ
ジューリングがされてない。
-KSWPオプションは有効であったが、コンパイラがSWPしても効果が無いと判断して、SWPしていなかった。
そのため、最適化制御行（OCL）で強制的にSWPしたところ、効果が見られた。

12

チューニングソース
 チューニング後のソース
 DOループの前の66行目に !ocl swp を追加
62
63
64
65
66

67
68
69
70
71
72
73
74
75

1 pp
1 p
1 p
1 p
1 p
1 p
1 p
1 p
1 p

*vocl loop,novrec(p1,div)*voption indep(p1,div)!cdir nodep(p1),nosync
追加!cdir on_adb(p1)!ocl swp<<< Loop-information Start >>><<< [PARALLELIZATION]<<< Standard iteration count: 69<<< [OPTIMIZATION]<<< SIMD<<< SOFTWARE PIPELINING<<< Loop-information End >>>v
do i=1,ico0max(m)v
kk=(i-1)/ijp8aav
ij=-kk*ijp8aa+(i-1)
v
jj=ij/iip8av
ii=-jj*iip8a+ijv
k1 =kk*ijp1v
nn1=k1+(jj*2)*iip1+(ii*2)+nn1rstv
n1 =k1+nn1+n1rsgv
nn2=nn1-ijp1

76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98

13

1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1

p

v
$
$
$
$
$
$
$
$
$
$
$
$
$
$
$
$
$
$
&

p
p
p

v
v
v

res=p1(nn1-iip1-1,0,m)*fj_clp0(1,i,m)+p1(nn2-1,1,m)*fj_clp0(2,i,m)+p1(nn11,0,m)*fj_clp0(3,i,m)+p1(nn1-1,1,m)*fj_clp0(4,i,m)+p1(nn1+iip1-1,0,m)*fj_clp0(5,i,m)+p1(nn2-iip1 ,1,m)*fj_clp0(6,i,m)+p1(nn1-iip1 ,0,m)*fj_clp0(7,i,m)+p1(nn1-iip1 ,1,m)*fj_clp0(8,i,m)+p1(nn2,1,m)*fj_clp0(9,i,m)+p1(nn1,0,m)+p1(nn1,1,m)*fj_clp0(11,i,m)+p1(nn2+iip1 ,1,m)*fj_clp0(12,i,m)+p1(nn1+iip1 ,0,m)*fj_clp0(13,i,m)+p1(nn1+iip1 ,1,m)*fj_clp0(14,i,m)+p1(nn1-iip1+1,0,m)*fj_clp0(15,i,m)+p1(nn2+1,1,m)*fj_clp0(16,i,m)+p1(nn1+1,0,m)*fj_clp0(17,i,m)+p1(nn1+1,1,m)*fj_clp0(18,i,m)+p1(nn1+iip1+1,0,m)*fj_clp0(19,i,m)-bp(n1)p1(nn1,0,m)=p1(nn1,0,m)-res*omegap
div(n1)=abs(res)enddo高コスト③sor8p_b_hチューニング結果(Loop1)
 Loop1チューニング
チューニングの結果、「浮動小数点演算待ち」と「整数演算待ち」のコストが減った。
その結果として、
メモリスループットが約68GB/secになった。

約11%
性能向上

性能は約11%向上した。

Performance
実行時間 (sec)浮動小数点演算
ピーク比

0.49
0.43

3.39%
3.79%

オリジナル
SWPL

MFLOPS

MIPS

501
561

浮動小数点演算
数

1272
1537

2.43E+08
2.43E+08

メモリスループットが3.79 ⇒ 4.24 GB/secに改善。
16スレッドで、約68GB/sとなっており、
メモリスループットネックとなっている。

Cache
L1I ミス率L1D ミス率
ロード・ストア数(/有効総命令数) (/ロード・ストア数)
オリジナル
SWPL

0.06%
0.06%

6.18%
5.09%

2.39E+08
2.94E+08

L1D ミス数
1.48E+07
1.50E+07

L1D ミス dm率(/L1D ミス数)

L1D ミス hwpf
率(/L1D ミス数)

42.23%
38.24%

57.75%
61.73%

L1D ミス swpf 率(/L1D ミス数)
0.03%
0.03%

L2 ミス率(/ロード・ストア数)
5.42%
4.38%

L2 ミス数
1.30E+07
1.29E+07

L2 ミス dm 率(/L2 ミス数)

L2 ミス pf 率(/L2 ミス数)

8.05%
9.23%

メモリスループット(GB/sec)

L2 スループット(GB/sec)

3.79
4.24

3.90
4.41

91.95%
90.77%

14

高コスト③sor8p_b_hチューニング結果(Loop1-8)
 Loop1-8チューニング結果（16スレッド並列のスレッド０）

Loop1-8は、ほぼ同じ処理のループです。
Loop1-8の合計でもLoop1と同様に、約11%の性能改善とメモリースループット
ネックになっています。

実行時間
(sec)
Loop1
Loop2
Loop3
Loop4
Loop5
Loop6
Loop7

15

Loop8

オリジナル
SWPL
オリジナル
SWPL
オリジナル
SWPL
オリジナル
SWPL
オリジナル
SWPL
オリジナル
SWPL
オリジナル
SWPL
オリジナル
SWPL

0.49
0.43
0.48
0.43
0.49
0.43
0.48
0.42
0.48
0.44
0.47
0.43
0.48
0.43
0.47
0.43

メモリ
スループット(GB/sec)
3.79
4.24
3.83
4.33
3.82
4.28
3.83
4.30
3.90
4.28
3.96
4.33
3.87
4.35
3.88
4.34

高コスト①のチューニングについて 高コスト①のsor7v_b
 -Karray_subscript, array_subscript_rank=2, array_subscript_elementlast=25

が既に指定されていて、 配列clp は次元移動により連続アクセス化がされていました。
 L2ミス(dm)率が高いため、配列clpのプリフェッチ(HWPF)が効果的でないのではないか？と言う観点で、HWPFを止めてSWPで検証しました。
また、コンパイラからは配列clp は 38 2 pp v
do 100 n1=n1rs(m)+ic,n1re(m),7
coef=omegav/clp(n1,10)
ストライドアクセスに見えるため、 404243 222 pp v $ resu=u(n1-iip1-1)*clp(n1,1) +u(n1-ijp1-1)*clp(n1,2)+u(n1- 1)*clp(n1,3) +u(n1+ijp1-1)*clp(n1,4)
44
2
$
+u(n1
+iip1-1)*clp(n1,5) +u(n1-ijp1-iip1 )*clp(n1,6)
ストライドアクセスオプションも
45
2$+u(n1-iip1)*clp(n1,7) +u(n1+ijp1-iip1 )*clp(n1,8)46
2$+u(n1-ijp1)*clp(n1,9) +u(n1)*clp(n1,10)
指定しました。
47
2$+u(n1+ijp1)*clp(n1,11)+u(n1-ijp1+iip1 )*clp(n1,12)48
2$+u(n1+iip1)*clp(n1,13)+u(n1+ijp1+iip1 )*clp(n1,14)
49
2
$
+u(n1
-iip1+1)*clp(n1,15)+u(n1-ijp1+1)*clp(n1,16)
オプションは以下です。
50
2$+u(n1+ 1)*clp(n1,17)+u(n1+ijp1+1)*clp(n1,18)
51
2
$
+u(n1
+iip1+1)*clp(n1,19)-Kprefetch_stride56
2$-bu(n1)58
2pv
resv=v(n1-iip1-1)*clp(n1,1) +v(n1-ijp1-1)*clp(n1,2)-Kprefetch_sequential=soft
59
2$+v(n1- 1)*clp(n1,3) +v(n1+ijp1-1)*clp(n1,4)
60
2
$
+v(n1
+iip1-1)*clp(n1,5) +v(n1-ijp1-iip1 )*clp(n1,6)-Kprefetch_stride (使用手引書より抜粋)
ループ内で使用されるキャッシュのラインサイズよりも大きなストライドでアクセスされる配列データに対して、
prefetch命令を使用したオブジェクトを生成します。
プリフェッチするアドレスが翻訳時に確定しないループを含みます。

61
62
63
64
65
66
67
72
：
：
103
106

2
2
2
2
2
2
2
2

$
$
$
$
$
$
$
$

2
1

p

100
110

+v(n1-iip1)*clp(n1,7) +v(n1+ijp1-iip1 )*clp(n1,8)+v(n1-ijp1)*clp(n1,9) +v(n1)*clp(n1,10)+v(n1+ijp1)*clp(n1,11)+v(n1-ijp1+iip1 )*clp(n1,12)+v(n1+iip1)*clp(n1,13)+v(n1+ijp1+iip1 )*clp(n1,14)+v(n1-iip1+1)*clp(n1,15)+v(n1-ijp1+1)*clp(n1,16)+v(n1+ 1)*clp(n1,17)+v(n1+ijp1+1)*clp(n1,18)+v(n1+iip1+1)*clp(n1,19)-bv(n1)

continue
continue

16

sor7v_b ①
251

do 110 ic=0,6<<< Loop-information Start >>><<< [PARALLELIZATION]<<<
Standard iteration count: 37<<< [OPTIMIZATION]
<<<
SIMD
<<<
SOFTWARE PIPELINING<<<
PREFETCH: 90<<<
v: 54, bw2: 6, bw: 6, w: 18, bv: 6<<< Loop-information End >>>v
do 100 n1=n1rs(m)+ic,n1re(m),7
coef=omegav/clp(n1,10)v
resu=u(n1-iip1-1)*clp(n1,1) +u(n1-ijp1-1)*clp(n1,2)$+u(n11)*clp(n1,3) +u(n1+ijp1-1)*clp(n1,4)$+u(n1+iip1-1)*clp(n1,5) +u(n1-ijp1-iip1 )*clp(n1,6)$+u(n1-iip1)*clp(n1,7) +u(n1+ijp1-iip1 )*clp(n1,8)$+u(n1-ijp1)*clp(n1,9) +u(n1)*clp(n1,10)$+u(n1+ijp1)*clp(n1,11)+u(n1-ijp1+iip1 )*clp(n1,12)$+u(n1+iip1)*clp(n1,13)+u(n1+ijp1+iip1 )*clp(n1,14)$+u(n1-iip1+1)*clp(n1,15)+u(n1-ijp1+1)*clp(n1,16)$+u(n1+
1)*clp(n1,17)+u(n1+ijp1+1)*clp(n1,18)
26
1
$
+u(n1
+iip1+1)*clp(n1,19)$-bu(n1)：100
continue
110
continue

プリフェッチオプションはデフォルト

38
40
42
43
44
45
46
47
48
49
50
51
56

2 pp
2 p
2 p
2
2
2
2
2
2
2
2
2
2

103 2
106 1

p

-Kprefetch_stride-Kprefetch_sequential=soft 指定-Kprefetch_stride,-Kprefetch_sequential=soft
により、配列 clp に対して、
prefetch命令が生成されて
いる。

39
41
43
44
45
46
47
48
49
50
51
52
57

2 pp
2 p
2 p
2
2
2
2
2
2
2
2
2
2

104
107

2
1

17

p

do 110 ic=0,6<<< Loop-information Start >>><<< [PARALLELIZATION]<<<
Standard iteration count: 37<<< [OPTIMIZATION]
<<<
SIMD
<<<
SOFTWARE PIPELINING<<<
PREFETCH: 222<<<
u: 54, v: 54, bu: 6, w: 54, bw2: 6<<<
bw: 6, clp: 36, bv: 6<<< Loop-information End >>>v
do 100 n1=n1rs(m)+ic,n1re(m),7
coef=omegav/clp(n1,10)v
resu=u(n1-iip1-1)*clp(n1,1) +u(n1-ijp1-1)*clp(n1,2)$+u(n11)*clp(n1,3) +u(n1+ijp1-1)*clp(n1,4)$+u(n1+iip1-1)*clp(n1,5) +u(n1-ijp1-iip1 )*clp(n1,6)$+u(n1-iip1)*clp(n1,7) +u(n1+ijp1-iip1 )*clp(n1,8)$+u(n1-ijp1)*clp(n1,9) +u(n1)*clp(n1,10)$+u(n1+ijp1)*clp(n1,11)+u(n1-ijp1+iip1 )*clp(n1,12)$+u(n1+iip1)*clp(n1,13)+u(n1+ijp1+iip1 )*clp(n1,14)$+u(n1-iip1+1)*clp(n1,15)+u(n1-ijp1+1)*clp(n1,16)$+u(n1+
1)*clp(n1,17)+u(n1+ijp1+1)*clp(n1,18)$+u(n1+iip1+1)*clp(n1,19)$-bu(n1)：100
continue
110
continuePA情報 [ sor7v_b ] sor7v_b ①
プリフェッチオプション
デフォルト
[秒]
1.2E+00(PA情報はstep=5)-Kprefetch_stride-Kprefetch_sequential=soft

L2ミスdm率が高く、「浮動小数点ロードメモリアクセス待ち」が発生↓約10%の性能向上

1.0E+00

プリフェッチ生成により約10%の性能向上

8.0E-01
6.0E-01

浮動小数点ロードキャッシュアクセス待ち・オプションによるプリフェッチの生成をした。

4.0E-01
2.0E-01

浮動小数点ロード
メモリアクセス待ち⇒ L2のdmミス削減

0.0E+00

Before T0

OPT Tune T0

Performance

SIMD
実行時間 (sec)浮動小数点演算
ピーク比

1.11
1.00

2.94%
3.25%

Before T0
OPT Tune T0

MFLOPS
434
480

浮動小数点演算
数

MIPS
458
632

4.82E+08 Before T0
4.82E+08 OPT Tune T0

SIMD命令率 SIMD演算命令率 SIMDロード命令率 SIMDストア命令率
SIMD命令率(/SIMD対象命令 (/SIMD対象演算 (/SIMD対象ロード (/SIMD対象ストア(/有効総命令数)
命令数)
命令数)
命令数)
数）
25.20%
31.41%
99.99%
0.00%
0.00%
20.20%
31.41%
99.99%
0.00%
0.00%
16スレッドで、約70.4GB/sとなっており、
メモリスループットネックとなっている。

Cache
L1I ミス率L1D ミス率
ロード・ストア数(/有効総命令数) (/ロード・ストア数)
Before T0
OPT Tune T0

0.05%
0.05%

25.99%
26.91%

2.82E+08
3.02E+08

L1D ミス数

L1D ミス dm率(/L1D ミス数)

L1D ミス hwpf
率(/L1D ミス数)

7.32E+07
8.14E+07

52.05%
54.89%

26.33%
0.01%

L1D ミス swpf 率(/L1D ミス数)
21.62%
45.11%
18

L2 ミス率L2 ミス数(/ロード・ストア数)
9.62%
9.97%

2.71E+07
3.02E+07

L2 ミス dm 率(/L2 ミス数)
23.15%
3.38%

L2 ミス pf 率 メモリスループット L2 スループット(/L2 ミス数) (GB/sec)(GB/sec)
76.85%
96.62%

3.63
4.40

8.42
10.37