2.2 3 次元非圧縮性熱流体計算プログラム COSMOS の測定評価

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) + b
x(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.26GHz
FX1
SPARC64 VII, 2.5GHz
1
4
京
SPARC64 VIII, 2.0GHz
1
8
FX10
SPARC64 IX, 1.848GHz
1
16
測定に利用した計算機の一覧を表 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割）は連立一次方程式を解くことに費やされる．
z
z
z

次元数（=格子点数）: 数10万～数100万
規則非対称帯行列
各行の非零要素数: ～ 25
解法: 定常反復法 (SOR, Gauss-Seidel)，
クリロフ部分空間法 (Bi-CGSTAB, ...)
z
重合格子の場合は，このような行列を格子ブロックごとに解くが，反復の途中
で境界領域の格子点に対する補間を行う．
(1),(2)の変化量が
小さくなるまで
繰り返し
Gauss-Seidel like
格子 m=1
m=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Ω ′ m
r
r
(2) 他の格子の内部の値から補間
x ←Ι x
Γm
r
xΓ m
[ ]
[ ]
[ Ωr′ ]
← Ι [xΩ ′ ]
他の格子の内部領域
Multiplicative（乗法的） Schwarz Method
2
7色オーダリングによるノード内並列化（オリジナル版）
SOR法を採用（実際には緩和係数を1としているので，Gauss-Seidel法）
各格子ブロックごとにマルチカラーオーダリングによる自動並
列化で対応

z
ベクトル機でのオーダリングをそのまま流用 ⇒ 7色飛び（1次元的）
配列のとり方も，それに合せて1次元並びになっている．
real*8 :: p0(n1max), bp(n1max)
real*8 :: clp(n1max,19)
do m=1,nvol
do ic=0,6
do n1=n1min(m)+ic,n1max(m),7
res=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
real*8 :: p0(n1max), bp(n1max)
real*8 :: fj_clp(19, n1max/7+1, 0:6)
ｃｌｐのみ色ごとに
do m=1,nvol
連続アクセスとした
do ic=0,6
i=1
do n1=n1rs(m)+ic,n1re(m),7
res=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/iip8a
ii=-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）での計測結果
z
z
スレッド並列：自動並列（VISIMPACT)
ｺﾝﾊﾟｲﾙｵﾌﾟｼｮﾝ: -Am -Kimpact,ocl,vppocl,tl_trt -X9 -w -Qt -Ectfiu -xdir=inline_dir
7色オーダリング
ルーチン名
sor7p_b_
sor7v_b_
sgs_
corvcn_
sor7p_
全体
プロファイラ
(count数)
31,883
10,455
7,608
4,036
2,995
64,966
8色オーダリング
CPU時間
[sec]
[%]
318.83
49.1%
104.55
16.1%
76.08
11.7%
40.36
6.2%
29.95
4.6%
649.66
100.0%
ルーチン名
4174.03 Mips
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.31
47.7%
104.82
17.6%
76.18
12.8%
40.18
6.8%
11.53
1.9%
594.42
100.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
運動方程式のSOR計算

Process 0 Thread 1
Process 0 Thread 2
Process 0 Thread 3
SGS項算出
9.0E+01
1.2E+02
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
4
5
7
7
6
8
k方向を交互に
二つの配列
p1, p2に持つ
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/iip8a
ii=-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)率が高い。
SIMD
実行時間 (sec)
sor7v_b T0
sgs T0
sor8p_b_h T0
③
②
①
Performance
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
L1I ミス率
(/有効総命令数)
sor7v_b T0
sgs T0
sor8p_b_h T0
0.05%
0.00%
0.02%
L2 ミス率
(/ロード・ストア数)
sor7v_b T0
sgs T0
sor8p_b_h T0
9.62%
0.34%
5.46%
L1D ミス率
(/ロード・ストア数)
25.99%
4.12%
6.24%
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 ミス数)
7.32E+08
2.66E+08
1.18E+09
L2 ミス pf 率
(/L2 ミス数)
52.05%
70.90%
41.06%
メモリスループット
(GB/sec)
76.83%
72.74%
92.53%
3.61
0.24
3.92
11
L1D ミス hwpf
率
(/L1D ミス数)
26.33%
3.43%
58.93%
L2 スループット
(GB/sec)
8.43
2.68
4.09
L1D ミス swpf 率
(/L1D ミス数)
21.62%
25.67%
0.00%
µ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
1
1
1
1
1
1
1
1
pp
p
p
p
p
p
p
p
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)/ijp8aa
v
ij=-kk*ijp8aa+(i-1)
v
jj=ij/iip8a
v
ii=-jj*iip8a+ij
v
k1 =kk*ijp1
v
nn1=k1+(jj*2)*iip1+(ii*2)+nn1rst
v
n1 =k1+nn1+n1rsg
v
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
オリジナル
SWPL
MFLOPS
3.39%
3.79%
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 ミス数)
42.23%
38.24%
L1D ミス hwpf
率
(/L1D ミス数)
L1D ミス swpf 率
(/L1D ミス数)
57.75%
61.73%
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)
91.95%
90.77%
3.79
4.24
L2 スループット
(GB/sec)
3.90
4.41
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_stride
56
2
$
-bu(n1)
58
2
p
v
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 ①
25
1
<<<
<<<
<<<
<<<
<<<
<<<
<<<
<<<
<<<
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
v
v
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 >>>
do 100 n1=n1rs(m)+ic,n1re(m),7
coef=omegav/clp(n1,10)
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)
-Kprefetch_stride,
-Kprefetch_sequential=soft
により、配列 clp に対して、
prefetch命令が生成されて
いる。
<<<
<<<
<<<
<<<
<<<
<<<
<<<
<<<
<<<
<<<
：
103 2
106 1
p
100
110
continue
continue
-Kprefetch_stride
-Kprefetch_sequential=soft 指定
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
v
v
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 >>>
do 100 n1=n1rs(m)+ic,n1re(m),7
coef=omegav/clp(n1,10)
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)
：
104
107
2
1
17
p
100
110
continue
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)
Before T0
OPT Tune T0
浮動小数点演算
ピーク比
1.11
1.00
2.94%
3.25%
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 ミス数
7.32E+07
8.14E+07
L1D ミス dm率
(/L1D ミス数)
L1D ミス hwpf
率
(/L1D ミス数)
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

Download Report