— 20 —SENAC Vol. 48, No. 3（2015. 7）［大規模科学計算システム］高速化推進研究活動報告第６号より転載

ベクトルコンピュータにおける高速化
小林広明 1 江川隆輔 1 小松一彦 1 岡部公起 1 大泉健治 2 小野 敏 2 山下 毅 2
佐々木大輔 2 森谷友映 2 齋藤敦子 2 撫佐昭裕 3 松岡浩司 3 渡部修 3 曽我 隆 4 山口健太 4
1

スーパーコンピューティング研究部
情報部情報基盤課
3
日本電気株式会社
4
NEC ソリューションイノベータ株式会社
2

本章では，本センターにおいてこれまで培われてきた SX-9 のベクトル性能を向上するための手法や事例を紹介する．これらの手法や事例は今後導入されるベクトルスーパーコンピュータ SX-ACE においても有用である．
1.

ベクトルコンピュータの特徴
近年，高性能計算機（HPC）システムを活用したシミュレーションの応用範囲は益々広がっており，同時にシミュレーションの規模も拡大し続けている．シミュレーションを行う研究者からの要求は計算規模を拡大すると同時に，計算結果を得るまでの時間は短縮したいというものである．このような研究者の要求を満たすため HPC システムの性能も年々向上し続けている．
HPC システムに採用されるプロセッサあるいはコア(以下コア)には大きく分類して，スカラ型とベクトル型がある．スカラ型コアは命令レベル並列性に基づきデータを一つずつ実行するのに対して，ベクトル型
コアはデータレベルの並列性に着目して複数の演算器(パイプライン)を用いて一度に複数の処理を同時に実行する．図 1-1 にスカラ命令とベクトル命令の動作を示す．この図は A(I)=B(I)+C(I)と D(I)=E(I)+F(I)の 2 つの計算式を 100 回繰り返す DO ループの処理を表している．スカラ命令では DO 変数が 1 の時のA(1)と D(1)の計算結果を求め，次に DO 変数が 2 の時の A(2)と D(2)の計算結果を求めるという処理を，
DO 変数が 100 になるまで繰り返す．それに対してベクトル命令は A(I)=B(I)+C(I)の計算を DO 変数 1 から
100 まで一度に行い，次に D(I)=E(I)+F(I)の計算を DO 変数 1 から 100 までを一度に行うことにより実行時間の短縮を図る．ベクトル型コアを搭載する HPC システムは，このベクトル命令の特性を活かすため，高いメモリ帯域を具備してメモリからコアへの大量のデータ供給を可能にしている．

図 1-1 スカラ命令とベクトル命令の様子

— 21 —

ベクトルコンピュータにおける高速化

表 1-1 に近年の HPC システムが搭載しているプロセッサの性能と内包するコア数，プロセッサ当たりのメ
モリバンド幅を示す．SX-9 以外のプロセッサでは複数の演算コアを搭載するマルチコア構成により演算性能の向上を図っているため，コア当たりのメモリバンド幅は小さくなる．しかしスカラプロセッサではメモリのバンド幅を最大限に使うためには複数のコアを使用する必要があるが，SX-ACE では 1 つコアで最大のメモリバンド幅を使うことが可能となっている．
このようにベクトルコンピュータは，高い理論演算性能とメモリ帯域を実装することで，実アプリケーショ
ンの実行において高い実効性能の達成を可能にしている．
表 1-1 HPC システムのプロセッサ性能と B/F 値
HPC システム名
SX-9
LX406Re-2
SR16000
FX10「京」
SX-ACE

名称―
Intel Xeon E5-269v2
POWER7SPARC64 IXfx
SPARC64 VIIIfx―

プロセッサ
理論性能(GFlop/s)
102.4
230.4
245.1
236.5
128.0
256.0

コア数
1
12
8
16
8
4

メモリバンド幅(GB/s)
256
59.7
128
85.3
64
256

— 22 —

2.

SENAC Vol. 48, No. 3（2015. 7）

高速化の手法
SX システムの超高速性を十分に引き出すために重要なことは，プログラムの中でベクトル命令によって処理される部分の比率，すなわちベクトル化率と，生成されるベクトル命令の効率を可能な限り高めることである．そのためには，プログラムの最適化状況の把握・分析を行い，主要な性能指標に基づき，高速化を行う必要がある．ここでは，そのためのツール類及び，それらの性能指標と各指標に沿った高速化の取り組みかたについて述べる．
2.1 性能解析ツール(1) 編集リスト
編集リストはベクトル化および自動並列化に関する情報をソースプログラムの左側に表示したリストであり，どのループがベクトル化されたかを確認することができる．編集リストはコンパイラオプ
ション「-R5」を指定することにより"入力ファイル名.L"という名前で出力される．
!編集リストFILE NAME:t5.f
PROGRAM NAME: sub
FORMAT LIST
LINE
LOOP
FORTRAN STATEMENT
1:
subroutine sub(a,b,c,d,z,ix)
2:real,dimension(100)::a,b,c,d,x,y,z
3:
integer ix(100)
4: V------>
do I = 1,100
5: |I
call sub2(x,a,b,I)
6: |y(I) = c(I) + d(I)7: |Sz(ix(I)) = z(ix(I)) + x(I) + y(I)8: V-----enddo
9:end

図 2-1 編集リストの例

主なループ情報，スカラ情報，手続インライン情報など編集リストの出力イメージを次に示す．
① ループ全体がベクトル化される場合
ベクトル化されたループに“V”が表示される．
V------>|
V------

do I=1,100a(I)=b(I)+c(I)enddo

図 2-2 ベクトル化された DO ループの編集リスト②

ベクトル化されない場合
ベクトル化されないループには“+”が表示される．

+------>|+------

do I=1,100
print *,a(I)enddo

図 2-3 ベクトル化されない DO ループの編集リスト

ベクトルコンピュータにおける高速化③— 23 —

部分ベクトル化の場合
ベクトル化不可の処理がある行には“S”が表示される．
※部分ベクトルとは，ループにベクトル化を阻害する部分が含まれている場合，その前後で自動的にループを分割し，ベクトル化可能な部分だけをベクトル化する拡張機能である．

V------>||
S
V------

do I =1,100a(I)=b(I)+c(I)print *,a(I)enddo

図 2-4 部分ベクトル化の編集リスト④

配列式をベクトル化した場合(1)ループの先頭と最後の行が同じである場合，ループの構造は“=”で表示される．
V=======

real a(90),b(90),ca(1:90)=b(1:90)+c

図 2-5 ベクトル化された配列式の編集リスト⑤

配列式をベクトル化した場合(２)ループ融合している場合は，その範囲について“V”で表示される．
※ループ融合とは，繰返し数が等しい DO ループまたは配列式が複数個，連続して存在している場合一つのループに融合することである．ただし，各ループで使われているデー
タの定義・参照関係に，融合によって矛盾が生じる場合は行わない．

real a(90),b(90),c
integer d(90),e(90)a(1:90)=b(1:90)+c
V------>|
d(1:90)=int(a(1:90))
V-----e(1:90)=d(1:90)+1

図 2-6 ループ融合とベクトル化された配列式の編集リスト⑥

手続呼出しがインライン展開された場合
インライン展開された手続がある行には“I”が表示される．
I

call sub2(x,a,b,c,I)

図 2-7 インライン展開の編集リスト⑦

多重ループが一重化された場合一重化されたループの外側ループに“W”，内側ループに“＊”が表示される．

W------>|＊---->|||＊---W------

do J =1,100
do I =1,100a(I,J)=b(I,J)+c(I,J)
enddo
enddo

図 2-8 多重 DO ループの一重化の編集リスト⑧

ループの入れ換えが行われた場合入れ換えた結果ベクトル化されるループに“X”が表示され，ベクトル化されなくなるループには“+”が表示される．

— 24 —

SENAC Vol. 48, No. 3（2015. 7）

X------>|+----->|||+----X------

do J =1,1000
do I =1,10
a(I,J)=b(I,J)+c(I,J)
enddo
enddo

図 2-9 DO ループの入れ換えの編集リスト⑨

ADB にデータがバッファリングされる場合
ADB を使用したベクトルロード/ストアがある行に文字“A”が表示される．

+------>|V----->||
A|V----+------

do J =1,100
do I =1,100a(I,J)=a(I,J+1)+b(I,J)
enddo
enddo

図 2-10 ADB の使用の編集リスト⑩

ベクトル化と並列化される場合
ベクトル化と並列化が行われたループに“Y”が表示される．

Y------>|
Y------

do i=1, 10000
c(i) = c(i) + a(i) * b(i)end do

図 2-11 ベクトル化と並列化された DO ループの編集リスト⑪

並列化される場合
並列化が行われたループに“P”が表示される．

P------>|V----->|||V----P------

do j = 1, 100
do i = 1, 100d(i,j) = 0.0d0
end doend do

図 2-12 並列化された DO ループの編集リスト⑫

複数の情報がある場合
一行に対して複数の情報がある場合，“M”が表示される．
※配列 a はベクトル化され，配列 b はループ長が短いためループが展開されるので，この一行には複数の情報がある．

M======a(1:10000)=0.0d0 ; b(1:3)=0.0d0

図 2-13 複数の情報がある場合の編集リスト(2) 簡易性能解析機能(FTRACE)
簡易性能解析情報はコンパイラオプションに「-ftrace」を指定することで，手続（サブルーチンや関数）ごとの性能解析情報を採取することができ，プログラム実行後に解析情報が標準出力ファイルに出力される．特に，チューニング対象とする手続を選択する際に活用すると良い．図2-14 に FTRACE 機能による解析リストの例を示す．FTRACE では実効性能と平均ベクトル長などの情報を手続（サブルーチンや関数）単位で取得することができる．

— 25 —

ベクトルコンピュータにおける高速化①
PROC.NAME②
FREQUENCY③
EXCLUSIVE
TIME[sec](④
AVER.TIME% )[msec]

⑤
MOPS

⑥
⑦
⑧
MFLOPS V.OP AVER.
RATIO V.LEN⑨⑩⑪⑫
VECTOR I-CACHE O-CACHEBANK CONFLICT
TIME
MISS
MISS CPU PORT NETWORK

subc
100
7.594( 42.2)
75.936 79426.5 52676.9 99.48 256.0
7.594 0.0000 0.0000
0.032
0.000
subb
100
5.227( 29.1)
52.272 77121.3 38261.9 99.22 256.0
5.227 0.0000 0.0000
0.000
0.000
suba
100
5.173( 28.7)
51.728 58600.1 1933.29 98.97 256.0
5.173 0.0000 0.0000
0.000
0.000
test
1
0.000( 0.0)
0.518
772.0
0.0 0.00
0.0
0.000 0.0000 0.0000
0.000
0.000
------------------------------------------------------------------------------------------------------------------------total
301
17.994(100.0)
59.780 7276.83 38902.9 99.28 256.0
17.993 0.0001 0.0001
0.032
0.000

図 2-14 解析リストの出力例
表示される各項目の意味は以下のとおりである．
①PROC.NAME：手続名②FREQUENCY：手続の呼び出し回数③EXCLUSIVE TIME

：手続の実行に要した占有の CPU 時間（秒)と，手続全体の実行に要した CPU 時間に対する比率④AVER.TIME：手続の 1 回の実行に要した平均 CPU 時間（ミリ秒）⑤MOPS：1 秒間に実行された演算数を 100 万単位で示した値⑥MFLOPS：1 秒間に実行された浮動小数点演算数を 100 万単位で示した値

⑦V.OP RATIO：ベクトル演算率⑧AVER V.LEN：平均ベクトル長⑨VECTOR TIME

：ベクトル命令実行時間（秒）⑩I-CACHE MISS

：命令キャッシュミスにより発生した合計時間（秒）⑪O-CACHE MISS

：オペランドキャッシュミスにより発生した合計時間（秒）⑫BANK CONFLICT：バンクコンフリクト時間（秒）CPU PORT：CPU ポート競合時間（秒）NETWORK：メモリネットワーク競合時間（秒）2.2. ベクトル化率の向上(1) ベクトル化率
ベクトル型スーパーコンピュータではプログラム中のベクトル処理可能な部分を高速に実行している．図 2-15 に同一のプログラムをスカラ処理する場合とベクトル処理する場合の実行時間に関する概念図を示す．

— 26 —

SENAC Vol. 48, No. 3（2015. 7）

T1（スカラ処理で実行した時の総実行時間）
T1V（T1のうちベクトル処理で実行可能な時間）

スカラ処理の場合

ベクトル処理不可部分

ベクトル処理の場合

ベクトル処理不可部分

ベクトル処理可能部分（DOループ部分）

ベクトル処理T2（ベクトル処理可能部分をベクトル処理で実行した時の総実行時間）

図 2-15 ベクトル処理による実行時間短縮のイメージ
ここで，あるプログラムをスカラ処理で実行する場合の総実行時間を T1，そのプログラムでベク
トル処理が可能な部分の実行時間を T1V とすると，ベクトル化率αは以下の式で定義される．
ベクトル化率 α =T1V
T1

また，そのプログラムをベクトル処理する場合の性能向上比 P は，スカラ処理性能とベクトル処
理性能の比をβとすると，
1

性能向上比 P =(1 - α) +αβと表される．
この式から，ベクトル化率と性能向上比の関係は図 2-16 のようになる（アムダールの法則）．この図からベクトル化率 80％程度では大きな性能向上は見られず，ベクトル化率 90％を超えたあた
りから急激に性能が向上していることがわかる．ベクトル型スーパーコンピュータで高い実効性能を得るためにはベクトル化率を 100%にできる限り近付ける必要がある．
β

性能向上比
（
倍）

P

0
0

20

40

60

80

100

ベクトル化率α （％）

図 2-16 ベクトル化率と性能向上比（アムダールの法則）

ベクトルコンピュータにおける高速化— 27 —

ベクトル化率を求めるため，プログラムをスカラ実行する場合とベクトル実行する場合のそれぞれの実行時間が必要である．しかし，これらの実行時間は同時に得られないのでベクトル化率の算出は困難である．そのため SX-9 ではベクトル化率の代わりにベクトル演算率をベクトル化の指標としている．FTRACE などの性能解析機能に出力される情報はベクトル演算率であり，ベクトル演算率は，プログラムで処理される全演算要素数に対するベクトル演算命令で処理される演算要
素数の割合で算出し，ほぼベクトル化率とみなせる指標であるため，後述のベクトル化率はベクト
ル演算率のこととする．
(2) ベクトル化を阻害する要因
ベクトル化を阻害する要因を以下に示す．
①

依存関係の有無が判断できない※前の繰り返しで定義した値を参照する場合や間接参照②

I/O 処理(OPEN/CLOSE/READ/WRITE 文)③

ユーザ関数(SUBROUTINE，FUNCTION)の呼び出し④

ループの繰り返し数がループ本体の実行前に決定しない⑤

ループの入口および出口が一つではない⑥

ベクトル化対象外の精度を使用している※2 バイト整数型，4 倍精度実数型，4 倍精度複素数型，1 バイト論理型，文字型,構造型
はベクトル化対象外の精度

2.3 平均ベクトル長の拡大(1) 平均ベクトル長
ベクトル処理を効率的に行う上で重要な指標にベクトル長がある．ベクトル長はベクトル化対象の DO ループの繰り返し回数のことであり，効率良くベクトル処理を行うためには，できるだけループ長を長くする必要がある．
一般に命令を発行してから結果が返ってくるまでに遅延時間が存在する．この遅延時間を立ち
上がり時間と呼ぶ．図 2-17 にベクトル処理における立ち上がり時間および演算時間の概念図を示す．網掛け部分は演算器が稼働する時間を示しており，演算を行うレジスタのデータを読み込
み，演算結果をレジスタに格納するまでの時間となる．図 2-18 にベクトル長と立ち上がり時間の関係を示す．立ち上がり時間はベクトル長が異なる場合でも一定である．よって総演算量が等しい場合，短ベクトル長処理を繰り返すよりも長ベクトル長処理を一括して行う方が実行時間を短くすることができる．このように高速化を図るためには，DO ループのベクトル長を十分に確保し，演算時間に対する立ち上がり時間の占める割合を低くすることが重要である．図 2-19 にベクトル長と立ち上がり時間の関係を示す．ベクトル処理の立ち上り時間が発生するため，ベクトル長が交差
ループ長より短い(4 以下)場合はベクトル化を行わない方が良い場合がある．交差ベクトル長とは，ベクトル化した場合とベクトル化しない場合とで実行時間が等しくなるループ長のことである．

— 28 —

SENAC Vol. 48, No. 3（2015. 7）

図 2-17 立ち上がり時間

図 2-18 ベクトル長と立ち上がり時間

図 2-19 ベクトル長と立ち上がり時間の関係

ベクトルコンピュータにおける高速化— 29 —

2.4 メモリアクセスの効率化(1) メモリ競合
SX-9 ではマルチメモリバンクシステムを採用しており，メモリを 32 のバンクグループに分割している．各バンクグループには 2 要素(16 バイト)ずつ配列データが格納される．
図 2-20 は連続アクセスとなる DO ループの例，図 2-21 は格納される配列データ(二次元の配列サイズが(128,64)のとき)の概念図を示しており，CPU からメモリにアクセスする際，CPU 内の全ポートを使用して効率的にメモリアクセスを行う．
!連続アクセス
do j=1,ny
do i=1,nx:
Rs(i,j)=Rs(i,j)+Cf(i,j)
end doend do

図 2-20 連続アクセスの DO ループの例

図 2-21 メモリバンクに格納される配列データの概念図(連続アクセス)

一方，図 2-22 はストライドアクセスとなる DO ループの例，図 2-23 は格納される配列データ(二次元の配列サイズが(128,64)のとき)の概念図を示しており，CPU からメモリをアクセスする際，アクセスするバンクグループが限定され CPU 内の特定のポート(00，16)だけを使用してメモリアクセスを行うことになる．このように CPU 内における同一のポートにロード・ストアが集中した時に発生する CPU ポート競合や，同一のメモリバンクへのアクセスなどで発生するメモリネットワーク競合
はメモリアクセスが遅延するため高速化の妨げとなる．

— 30 —

SENAC Vol. 48, No. 3（2015. 7）!ストライドアクセス
do j=1,ny
do i=1,nx,32:
Rs(i,j)=Rs(i,j)+Cf(i,j)
end doend do

図 2-22 ストライドアクセスの DO ループの例

図 2-23 メモリバンクに格納される配列データのイメージ(ストライドアクセス)

以降にメモリアクセス性能を改善するための指針として，メモリアクセスパターンの改善と ADBの活用について述べる．
(2) メモリアクセスパターンの改善
図 2-24 にメモリへアクセスパターンの種類を示す．アクセスパターンの種類としては，連続アクセス，ストライドアクセス，間接アクセスの 3 種類に分けることができる．

ベクトルコンピュータにおける高速化— 31 —

図 2-24 メモリのアクセスパターンの種類
図 2-21 に示すように，連続アクセスの場合が最も効率的なメモリアクセスを行うことができる．
ストライドアクセスの場合には，図 2-23 に示すように，使用できるバンク数が限られてくる．この場
合，ループ入れ換えなどにより，なるべく最初の次元でメモリアクセスを行うようにする．その際，ループインデックスが繰り返しごとに 1 ずつ増加あるいは減少するようにできでば連続アクセスとなり，より効率的なメモリアクセスが可能となる．間接アクセスの場合，間接アクセスとなる配列(図 2-24の間接アクセスにおける配列 A)が何度も参照されるならば，最初に，この配列のデータを作業用の配列にコピーしておき，以降のループでは作業用の配列をアクセスするようにすることにより，連続アクセスにすることができる．
(3) ADB の活用
SX-9 では CPU-主記憶装置間に ADB(Assinable Data Buffer)と呼ばれるベクトルデータを選択的にバッファする機能を有している．この ADB を活用することによって，メモリアクセス性能を改善できる．
コンパイラは，再利用性があると判断した配列を自動的に ADB にバッファリングするが，コンパイラが判断できなかった配列については，ON_ADB 指示行を用いることによって，明示的に該当配列を ADB にバッファリングすることを指示することができる．例えば，間接アクセスにおいて，間接アクセスとなる配列(図 2-24 の間接アクセスにおける配列 A)に対し ON_ADB 指示行を用いてADB にバッファリングしておくことにより，この配列のロードを高速に行えるようになる．
なお，ADB の容量は限られているため，ADB にバッファリングするデータのサイズに注意し，再利用性のあるデータを ADB にバッファリングするように指定することが重要である．

— 32 —

3.

SENAC Vol. 48, No. 3（2015. 7）

ベクトル化率向上の事例

3.1. ループ内の WRITE 文の括り出し(1) チューニング方針
SX-9 では，DO ループの一部にベクトル化対象外の文が存在すると，ベクトル化されない．もしくは，部分ベクトル化となり，完全にはベクトル化されない． 図 3-1 にチューニング前のコードを示す．このコードでは，DO ループ内にエラー検出時に実行される WRITE 文が含まれている．入出力文はベクトル化対象外の処理であるため，DO ループ全体がベクトル化されていない．図3-2 のチューニング前の FTRACE 情報を確認すると，ベクトル化率が 0.0%であり，全くベクトル化
されていないことが分かる．そこで，ベクトル化阻害要因である WRITE 文を DO ループの外に移動させることにより，DO ループをベクトル化可能にする．
!チューニング前
12: +------>
13: |+----->:
18: ||
19: ||
20: ||
21: ||
22: ||
23: ||
24: ||
25: ||
28: ||
29: ||
30: ||
31: ||
32: ||
34: ||
36: ||
37: |+----38: +------

DO J=2, JF
DO I=2, IF
IF (HZ(I,J) .GE. -30.0) THENZZ = Z(I,J,1) - RX*(M(I,J,1)-M(I-1,J,1))&- RY*(N(I,J,1)-N(I,J-1,1))
IF (ABS(ZZ) .LT. GX) ZZ = 0.0
DD = ZZ + HZ(I,J)
IF (DD .LT. GX) DD = 0.0
DZ(I,J,2) = DD
Z(I,J,2) = DD - HZ(I,J)
IF(ABS(Z(I,J,2)).GT.100.0) THEN
NC=1
WRITE(*,'(A24,3I6)')'Over flow Z at (K,I,J) :', KK,I,J
WRITE(*,*)'Within Region :', NREG
WRITE(*,*)'Computation is unstable.'
END IF
END IF
END DOEND DO

図 3-1 WRITE 文括り出し前のコード
FREQUENCY

EXCLUSIVE
AVER.TIME
TIME[sec]( % )[msec]
40
6.863( 3.1)
171.564

MOPS
558.3

MFLOPS V.OP AVER.
RATIO V.LEN
195.4 0.00
0.0

VECTOR I-CACHE O-CACHEBANK CONFLICT
PROC.NAMETIMEMISSMISS CPU PORT NETWORK
0.000
0.000
2.239
0.000
0.000 【チューニング前】

図 3-2 WRITE 文括り出し前の FTRACE 情報(2) チューニング内容
図 3-3 にチューニング後のコードを示す．チューニング前の DO ループではエラー検出時にWRITE 文によるエラー情報の出力を行っているが，本修正においては，DO ループ内ではエラー検出時にフラグ検出用配列にフラグをセットするようにした．そして，WRITE 文によるエラー情報の出力処理は DO ループ外に移動させ，出力が必要な場合のみ処理を行うようにした．これにより，
DO ループ内にはベクトル化対象外の文がなくなるため，ベクトル化が可能となる．

— 33 —

ベクトルコンピュータにおける高速化!チューニング後! フラグの初期化
23:
IFLG1 = IF * JF + 1
24: W*=====
IDX = IFLG1:! メイン処理!
エラーの検出のみ実施して，エラー出力処理をループの外に出す
27: +------>
DO J=2, JF
DO I=2, IF
28: |V----->:
33: ||IF (HZ(I,J) .GE. -30.0) THEN
34: ||ZZ = Z(I,J,1) - RX*(M(I,J,1)-M(I-1,J,1))
35: ||&- RY*(N(I,J,1)-N(I,J-1,1))
36: ||IF (ABS(ZZ) .LT. GX) ZZ = 0.0
37: ||
DD = ZZ + HZ(I,J)
38: ||
39: ||IF (DD .LT. GX) DD = 0.0
40: ||
DZ(I,J,2) = DD
41: ||
Z(I,J,2) = DD - HZ(I,J):
45: ||
IF(ABS(Z(I,J,2)).GT.100.0) THEN
IDX(I,J) = I + (J - 1) * IF
53: ||
54: ||END IF56: ||END IF
57: |V----END DO
58: +-----END DO:! エラー出力処理
61: +------>
DO J=2, JF
62: |V----->
DO I=2, IF
63: ||IFLG1 = MIN(IDX(I,J), IFLG1)
64: |V----ENDDO
65: +-----ENDDO
66:
IF ( IFLG1 .LT. IF * JF + 1 ) THEN67:
INDX_J = (IFLG1 - 1) / IF + 1
68:
INDX_I = IFLG1 - (INDX_J - 1) * IF
69:
NC=1
70:
71:
WRITE(*,'(A24,3I6)')'Over flow Z at (K,I,J) :',73:
WRITE(*,*)'Within Region :', NREG
74:
WRITE(*,*)'Computation is unstable.'
75:
ENDIF

図 3-3

WRITE 文括り出し後のコード(3) 性能分析
図 3-4 にチューニング前後の性能値を示す．
FREQUENCY4040

EXCLUSIVE
TIME[sec](

AVER.TIME[msec]MOPS% )MFLOPS V.OP AVER.
RATIO V.LEN

6.863( 3.1)
0.079( 0.0)

171.564
1.966

558.3
63568.9

195.4 0.00
0.018871.6 99.54 238.9

VECTOR I-CACHE O-CACHEBANK CONFLICT
PROC.NAMETIMEMISSMISS CPU PORT NETWORK
0.000
0.078

0.000
0.000

2.239
0.000

0.000
0.004

0.000 【チューニング前】
0.012 【チューニング後】

図 3-4 WRITE 文括り出し前後 FTRACE 情報
DO ループ全体がベクトル化されたことにより，ベクトル化率が 0.0%から 99.5%に向上し，実行時間は 6.86 秒から 0.08 秒に短縮することができた．

— 34 —

SENAC Vol. 48, No. 3（2015. 7）

3.2. 部分ベクトル化の回避(1) チューニング方針
図 3-5 にチューニング前のコードを示す．39 行目の変数 vmax は前の繰り返しで定義された値を参照しているためベクトル化の阻害要因となるが，コンパイラは最大値を求める演算をマクロ演
算に置き換えてベクトル化を行おうとする．最大値もしくは最小値を求める演算がある場合，ベクト
ルパイプライン毎にベクトル演算を行い，最後に各演算結果のマスクを取るマクロ演算を適用してベクトル化を行う．しかし，39 行目で求める vmax の値を 40 行目で参照する必要があるため，ベク
トル化を阻害する依存関係となり，39 行目と 40 行目をスカラ処理する部分ベクトル化となる．図3-6 に示すチューニング前の FTRACE 情報を確認するとベクトル化率が 67.5%で十分にベクトル化がされていないことが分かる．ベクトル化を阻害する依存関係を排除し，マクロ演算を適用してループ全体をベクトル化する．
!チューニング前
23: +------>
24: |V----->:
37: ||
38: ||
S
39: ||
40: ||
S
41: |V----42: +------

do iy = 1,ny
do ix = 1,nx
absv =max(absvxi,absvxj,absvyi,absvyj)
cf= sqrt(cs2(ix,iy)+ca2(ix,iy))
vmax = max(vmax,absv+cf)
dtmin = min(dtmin,dlmin/vmax)
enddo
enddo

図 3-5 部分ベクトル化回避前のコード
FREQUENCY821

EXCLUSIVE
TIME[sec](

AVER.TIME[msec]MOPS% )MFLOPS V.OP AVER.
RATIO V.LEN

VECTOR I-CACHE O-CACHEBANK CONFLICT
PROC.NAMETIMEMISSMISS CPU PORT NETWORK

26.020(7.2)

31.693

1664.6496.5 67.53 255.9

2.017

0.001

1.503

0.003

1.113 【チューニング前】

図 3-6 部分ベクトル化回避前の FTRACE 情報(2) チューニング内容
チューニング前では vmax の最大値判定と更新，dtmin の最小値判定と更新をループ回数分行っている．dtmin の判定式は不変値 dlmin を vmax で除算しているため，dtmin は vmax が更新されたときのみ dtmin が更新される．そのため，DO ループで vmax の最大値を求め，DO ループ終了後に一回 dtmin の判定を行うことで同じ結果が得られる．図 3-7 にチューニング後のコードを示す．前述したとおり，MAX 関数で求める vmax を MIN 関数で参照することが依存関係であるため，MIN 関数の演算を DO ループの外に移動することで依存関係が解消され，vmax を求める処
理はマクロ演算が適用されてベクトル化が行われる．

— 35 —

ベクトルコンピュータにおける高速化!チューニング後
23: +------>
24: |V----->:
37: ||
A
38: ||
39: ||
A
40: |V----41: +-----42:

do iy = 1,ny
do ix = 1,nx
absv = max(absvxi,absvxj,absvyi,absvyj)
cf= sqrt(cs2(ix,iy)+ca2(ix,iy))
vmax = max(vmax,absv+cf))
enddo
enddo
dtmin = min(dtmin,dlmin/vmax)

図 3-7 部分ベクトル化回避後のコード(3) 性能分析
図 3-8 にチューニング前後の性能値を示す．
FREQUENCY

821
821

EXCLUSIVE
TIME[sec](

AVER.TIME[msec]MOPS% )MFLOPS V.OP AVER.
RATIO V.LEN

VECTOR I-CACHE O-CACHEBANK CONFLICT
PROC.NAMETIMEMISSMISS CPU PORT NETWORK

26.020( 7.2)
2.273( 0.8)

31.6932.769

1664.614592.6

496.5 67.53 255.9
4215.6 99.19 255.9

2.017
1.971

0.001
0.000

1.503
0.124

0.003
0.019

1.113 【チューニング前】
1.029 【チューニング後】

図 3-8 部分ベクトル化回避前後の FTRACE 情報
部分ベクトルが解消しベクトル化率が 67.5%から 99.2%に改善したことで実行時間は 26.0 秒から
2.3 秒に短縮することができた．

4.

平均ベクトル長の拡大の事例

4.1. ループに関する情報の追加による最適化の促進(1) チューニング方針
多重ループの最内にあるループが自動ベクトル化の対象となるため，最内ループのベクトル長が短い場合，十分な性能を発揮できない．図 4-1 のチューニング前のコードではループの中に配列式や配列関数を使用しており，その部分が最内ループとなりベクトル化の対象となる．図4-2 の FTRACE 情報では平均ベクトル長が 67.0 であり十分なベクトル長でないことが分かる．ここでは，コンパイラがループ長を判断してベクトル長を伸ばす方法で最適化を行う．ループの演算
処理量が多いなどの問題からこのような最適化が行われない場合もあるが，ループに関する情報をコンパイラに明示することでベクトル長を伸ばす最適化が行われる場合がある．このような高速化の事例を以下に示す．

— 36 —

SENAC Vol. 48, No. 3（2015. 7）!チューニング前
5:real*8:: y(Nch), xp(Ndim,Nch,Np)6:
real*8:: loglhTP(Np), dy(Nch):
15: V------> do ip=1,Np
dy = y - xp(1,:,ip) &16: |V===== A
17: |- xp(Mu+1,:,ip) - xp(Mu+Ms,:,ip)
18: |loglhTP(ip) = - 0.5d0 * ( &19: |V===== A
20: |dble(Nch) * log( 2.0d0 * pi )21: |+ log( product( xp(Ndim,:,ip) ))22: |+ sum( dy * dy / xp(Ndim,:,ip) )23: |)24: V-----enddo&&&

図 4-1 ループに関する情報を追加する前のコード
FREQUENCY60

EXCLUSIVE
TIME[sec](

AVER.TIME[msec]

MOPS% )73.899( 73.1)

1231.647

2404.9MFLOPS V.OP AVER.
RATIO V.LEN
576.8 89.76

67.0

VECTOR I-CACHE O-CACHEBANK CONFLICT
PROC.NAMETIMEMISSMISS CPU PORT NETWORK
63.853

0.000

0.000

0.001

41.961 【チューニング前】

図 4-2 ループに関する情報を追加する前の FTRACE 情報(2) チューニング内容
ベクトル化対象となる配列式と配列関数のループ長(変数 Nch)は 3，外側 ip の DO ループの長さは 1,280,000 であり，この情報を追加することでコンパイラはより効率の良い最適化を実施する．
図 4-3 にチューニング後コードを示す．PARAMETER 文を追加することでコンパイラは最内 DO
ループのループ長が 3 であることを認識するため，ループの展開を行い，外側のループでベクトル化を行う．ただし，このチューニング内容は変数 Nch，Np の値が不変であることを前提としている場合のみ使用することが出来るので注意が必要となる．
!チューニング後
6:parameter(Nch=3,Np=1280000)
7:
real*8:: y(Nch), xp(Ndim,Nch,Np)8:
real*8:: loglhTP(Np), dy(Nch)
15: V------> do ip=1,Np
dy = y - xp(1,:,ip)
16: |*=====
17: |- xp(Mu+1,:,ip) - xp(Mu+Ms,:,ip)
18: |loglhTP(ip) = - 0.5d0 * ( &19: |*=====
20: |dble(Nch) * log( 2.0d0 * pi )21: |+ log( product( xp(Ndim,:,ip) ))22: |+ sum( dy * dy / xp(Ndim,:,ip) )23: |)24: V------ enddo&&&

図 4-3 ループに関する情報を追加後のコード

— 37 —

ベクトルコンピュータにおける高速化(3) 性能分析
図 4-4 にチューニング前後の性能値を示す．
FREQUENCY60
60

EXCLUSIVE
TIME[sec](

AVER.TIME[msec]MOPS% )MFLOPS V.OP AVER.
RATIO V.LEN

VECTOR I-CACHE O-CACHEBANK CONFLICT
PROC.NAMETIMEMISSMISS CPU PORT NETWORK

73.899( 73.1)
0.289( 1.0)

1231.647
4.813

2404.923763.3

576.8 89.76 67.0
12232.7 99.60 256.0

63.853
0.267

0.000
0.000

0.000
0.000

0.001
0.017

41.961 【チューニング前】0.111 【チューニング後】

図 4-4 ループに関する情報を追加する前後の FTRACE 情報
外側のループをベクトル化対象とすることで平均ベクトル長が 67.0 から 256.0 となり，十分な性能を発揮することで，実行時間が 73.9 秒から 0.3 秒に高速化された．
4.2. 部分配列を DO ループに書き替え，ループ交換を実施(1) チューニング方針
チューニング前コードおよびチューニング前の FTRACE 情報は前項と同じ図 4-1 と図 4-2 である．前項では，変数 Nch，Np の値が不変であるという条件のもと PARAMETER 文を追加することでコンパイラの最適化を促進させた．ここでは変数 Nch，Np の値が不定の場合のチューニング
方法を示す．
(2) チューニング内容
変数 Nch，Np の値が不定の場合のチューニング方法を示すが，ここで Nch の値が小さく，Npの値が十分に大きいことが分かっているものとする．図 4-5 にチューニング後コードを示す．配列
式と配列関数を DO ループの形に変形し，ループ分割を行っている．このとき作業配列を新たに用意し使用する．また，Nch＜Np であるため，Np の DO ループがベクトル化の対象となるように最内にループの入れ換えを行っている．
!チューニング後
5:real*8:: y(Nch), xp(Ndim,Nch,Np)6:
real*8:: loglhTP(Np), dy(Nch)8:real*8:: wrk_dy,wrk
9:real*8:: wrk_sum(Np),wrk_prd(Np)
24: V------> wrk_prd=1.d0
25: V-----wrk_sum=0.d0
26: +------> do j=1,Nch
27: |V-----> do ip=1,Np
28: ||
A
wrk_dy=y(j)-xp(1,j,ip)
29: ||-xp(Mu+1,j,ip)-xp(Mu+Ms,j,ip)
30: ||
A
wrk_prd(ip)=wrk_prd(ip)31: ||* xp(Ndim,j,ip)
32: ||
A
wrk_sum(ip)=wrk_sum(ip)33: ||+(wrk_dy**2/xp(Ndim,j,ip))34: |V----- end do
35: +-----end do
36:
wrk=dble(Nch)*log(2.0d0*pi)
37: V------> do ip=1,Np
38: |A
loglhTP(ip)=-0.5d0*(wrk
39: |+log(wrk_prd(ip))+wrk_sum(ip))
40: V------ enddo

&
&
&

&

図 4-5 部分配列を書き換え後のコード

— 38 —

SENAC Vol. 48, No. 3（2015. 7）(3) 性能分析
図 4-6 にチューニング前後の性能値を示す．
FREQUENCY60
60

EXCLUSIVE
TIME[sec](

AVER.TIME[msec]MOPS% )MFLOPS V.OP AVER.
RATIO V.LEN

VECTOR I-CACHE O-CACHEBANK CONFLICT
PROC.NAMETIMEMISSMISS CPU PORT NETWORK

73.899( 73.1)
0.314( 1.1)

1231.647
5.232

2404.925804.5

576.8 89.76 67.0
11742.5 99.54 256.0

63.853
0.287

0.000
0.000

0.000
0.000

0.001
0.026

41.961 【チューニング前】
0.101 【チューニング後】

図 4-6 部分配列を書き換え前後の FTRACE 情報
外側のループをベクトル化対象とすることで平均ベクトル長が 67.0 から 256.0 となり，実行時間
が 73.9 秒から 0.3 秒に短縮された．
5.

メモリ競合の回避の事例

5.1. ADB によるメモリアクセス低減(1) チューニング方針
0.(3)で述べているように，ADB を利用することにより，より効率的なメモリアクセスを行うことができる．その事例として，図 5-1 を挙げる．
図 5-1 において最内側の DO ループがベクトル化されているが，配列 b はこのベクトル化されたループに入る度にメモリからのロードを行うため，メモリ競合が発生する．図 5-2 の FTRACE 情報が示すように，チューニング前のコードでは，実行時間 26.6 秒に対してバンクコンフリクト時間
が 18.1 秒を占めている．この配列 b に着目してみると，ベクトル化されている DO ループの外側のDO ループ（DO 変数 l のループ）による繰り返しによって配列の参照位置が変わることはなく，最内 DO ループにおいては，同じ添え字 j のデータを毎回メモリからロードすることが分かる．従って，
配列 b は最内 DO ループにおいて再利用性があることが分かる．
通常，コンパイラはループに対する最適化処理の中で配列の再利用性解析を行う．しかし，本事例のように，ループからの飛び出しがあるような場合にはコンパイラはループの最適化処理を行わない．そのため，配列の再利用性解析も行われず，再利用性がある配列に対してもコンパイ
ラは自動で ADB に載せるための処理を行わない．このような場合，コンパイラ指示行により明示的に配列を ADB に載せることを指示し，ADB 経由でのメモリアクセスを行うことにより，メモリに直接アクセスする頻度を削減し，バンクコンフリクト時間を短縮する．
!チューニング前
335: ||+---->
336: |||V--->
337: ||||
338: |||V--339: |||340: |||
341: ||+----

do l=1,lmax
do i=is,imax(j)+1
if(a(l).lt.b(i,j)) exit
end dois=iiwk(l)=i
end do

図 5-1 ADB によるメモリアクセス低減前のコード
FREQUENCY126

EXCLUSIVE
TIME[sec](

AVER.TIME[msec]MOPS% )MFLOPS V.OP AVER.
RATIO V.LEN

VECTOR I-CACHE O-CACHEBANK CONFLICT
PROC.NAMETIMEMISSMISS CPU PORT NETWORK

26.574(100.0)

210.906

1166.8342.1 88.72 114.4

22.488

0.000

0.000

0.012

図 5-2 ADB によるメモリアクセス低減前の FTRACE 情報

18.072 【チューニング前】— 39 —

ベクトルコンピュータにおける高速化(2) チューニング内容
ON_ADB 指示行を用いて，配列 b を ADB に載せることを指示する．チューニング後のコードを
図 5-3 に示す．これにより 1 回目のロード命令はメモリから行われるが，2 回目以降のロード命令
は ADB から行われることになる．この結果，メモリへのアクセスを削減し，バンクコンフリクトによる実行時間の増加を解消する．
!チューニング後
335: ||+---->
336: |||
337: |||V--->
338: ||||
A
339: |||V--340: |||341: |||
342: ||+----

do l=1,lmax!cdir on_adb(b)do i=is,imax(j)+1
if(a(l).lt.b(i,j)) exit
end dois=iiwk(l)=i
end do

図 5-3 ADB によるメモリアクセス低減後のコード(3) 性能分析
図 5-4 にチューニング前後の性能値を示す．
FREQUENCY

126
126

EXCLUSIVE
TIME[sec](

AVER.TIME[msec]MOPS% )MFLOPS V.OP AVER.
RATIO V.LEN

VECTOR I-CACHE O-CACHEBANK CONFLICT
PROC.NAMETIMEMISSMISS CPU PORT NETWORK

26.574(100.0)
9.547( 99.9)

210.906
75.7711166.83247.7342.1 88.72 114.4
952.3 88.72 114.4

22.488
5.455

0.000
0.000

0.000
0.000

0.012
0.012

18.072 【チューニング前】
1.036 【チューニング後】

図 5-4 ADB によるメモリアクセス低減前後の FTRACE 情報
メモリへのアクセスが減少することにより，バンクコンフリクト時間が短縮され，実行時間が 26.6
秒から 9.5 秒に短縮することができた．
5.2. ループ融合によるベクトルロード・ストアの削減(1) チューニング方針
通常，コンパイラによる最適化では，同じループ構造の DO ループが連続していると自動でループ融合を行う．ループ融合を行うことで同一の配列データに対するメモリアクセス(ロード，ストア)回数を減らすことができる．図 5-5 にチューニング前のコードを示す．チューニング前コードでは同じループ構造の DO ループが連続しているが OpenMP 指示行を含んでいるためループ融合が行われない．コンパイラは PARALLEL リージョン単位で最適化を行うため，ループ融合が行われず連続する複数の DO ループで同じ配列データのメモリアクセスが発生し，配列 sm3dh はメモリ
ロード 2 回，メモリストア 2 回の命令が発生している．図 5-6 の FTRACE 情報を確認すると実行
時間 13.1 秒に対してバンクコンフリクト時間が 7.4 秒を占めているため，ループ融合によりメモリア
クセス回数を減らしバンクコンフリクト時間が短縮される．

— 40 —

SENAC Vol. 48, No. 3（2015. 7）!チューニング前
7414:
:
7419:
7420: P------>
7421: |+----->
7422: ||V---->:
7437: |||
A
7438: |||:
7443: ||V---7444: |+----7445: P-----7446:
7447:
:
7454:
:
7457:
:
7462:
7463: P------>
7464: |+----->
7465: ||V---->
7466: |||
A
7467: |||:
7473: ||V---7474: |+----7475: P-----7476:
7477:
:
7498:
:
7502:
7503: P------>
7504: |+----->
7505: ||V---->
7506: |||
A
7510: |||
A:
7519: ||V---7520: |+----7521: P-----7522:
7523:!$OMP PARALLEL SHARED (sm3dh,sn3dh,sl3dh,sr3dh,se3dh!$OMP DOdo iz=1,lz
do iy=1,ly
do ix=3,lx-2&

sm3dh(ix,iy,iz)=sm3d(ix,iy,iz)+sdltxinv*( du(ix-1,iy,iz)-du(ix,iy,iz) )

enddo
enddo
enddo
!$OMP END DO NOWAIT!$OMP END PARALLEL
sdltthalfgx=0.5d0*adltt*sgravx!$OMP PARALLEL SHARED (sm3dh,sn3dh,sl3dh,sr3dh,se3dh!$OMP DOdo iz=1,lz
do iy=1,ly
do ix=3,lx-2
sm3dh(ix,iy,iz)=sm3dh(ix,iy,iz)&+sdltthalfgx*sr3d(ix,iy,iz)
enddo
enddo
enddo
!$OMP END DO NOWAIT!$OMP END PARALLEL!$OMP PARALLEL SHARED (sm3dh,sn3dh,sl3dh,sr3dh,se3dh!$OMP DOdo iz=1,lz
do iy=1,ly
do ix=3,lx-2
sr3dhinv=1.0d0/sr3dh(ix,iy,iz)
su3dh(ix,iy,iz)=sr3dhinv*sm3dh(ix,iy,iz)
enddo
enddo
enddo
!$OMP END DO NOWAIT!$OMP END PARALLEL

図 5-5 ループ融合によるメモリアクセス低減前のコード
FREQUENCY200

EXCLUSIVE
TIME[sec](

AVER.TIME[msec]MOPS% )MFLOPS V.OP AVER.
RATIO V.LEN

VECTOR I-CACHE O-CACHEBANK CONFLICT
PROC.NAMETIMEMISSMISS CPU PORT NETWORK

13.111(

2.6)

65.55736092.3

16377.9 99.26 129.0

13.102

0.002

0.005

0.802

6.589 【チューニング前】

図 5-6 ループ融合によるメモリアクセス低減前の FTRACE 情報(2) チューニング内容
コンパイラによる自動ループ融合が行われない場合，ユーザチューニングによるループ融合が効果的である．チューニング後コードを図 5-7 に示す．連続する複数の DO ループは同じループ構造であり，一つの DO ループに修正することでコンパイラが最適化を行う．ループ融合後の配列

— 41 —

ベクトルコンピュータにおける高速化

sm3dh に注目すると，演算で求めた結果をレジスタに格納しその後の演算ではレジスタ内のデー
タを参照することになるため，メモリアクセス回数はメモリストアの 1 回になり，メモリアクセス回数を削減することができる．表 5-1 にチューニング前後のメモリアクセス回数の一覧を示す．
!チューニング後
7858:
7859:
7860:
7861: P------>
7862: |+----->
7863: ||V---->:
7868: |||
A
7869: |||:
7875: |||
7876: |||:
7881: |||
7882: |||:
7891: ||V---7892: |+----7893: P-----7894:
7895:

sdltthalfgy=0.5d0*adltt*sgravy!$OMP PARALLEL PRIVATE (ix,iy,iz,sr3dhinv)!$OMP DOdo iz=1,lz
do iy=1,ly
do ix=3,lx-2&

sm3dh(ix,iy,iz)=sm3d(ix,iy,iz)+sdltxinv*( du(ix-1,iy,iz)-du(ix,iy,iz) )&sm3dh(ix,iy,iz)=sm3dh(ix,iy,iz)+sdltthalfgx*sr3d(ix,iy,iz)
sr3dhinv=1.0d0/sr3dh(ix,iy,iz)
su3dh(ix,iy,iz)=sr3dhinv*sm3dh(ix,iy,iz)

enddo
enddo
enddo
!$OMP END DO NOWAIT!$OMP END PARALLEL

図 5-7 ループ融合によるメモリアクセス低減後のコード
表 5-1 ループ融合によるメモリアクセス低減前後のメモリアクセス命令数
チューニング前
36
17

メモリロードメモリストアチューニング後169(3) 性能分析
図 5-8 にチューニング前後の性能値を示す．
FREQUENCY200
200

EXCLUSIVE
TIME[sec](

AVER.TIME[msec]MOPS% )MFLOPS V.OP AVER.
RATIO V.LEN

VECTOR I-CACHE O-CACHEBANK CONFLICT
PROC.NAMETIMEMISSMISS CPU PORT NETWORK

13.111(
8.791(

2.6)
1.8)

65.55743.95336092.3
52951.2

16377.9 99.26 129.0
24427.6 99.47 129.0

13.102
8.770

0.002
0.002

0.005
0.003

0.802
0.769

6.589 【チューニング前】
3.745 【チューニング後】

図 5-8 ループ融合によるメモリアクセス低減前後の FTRACE 情報
メモリアクセス回数が減少することでバンクコンフリクト時間が削減され，実行時間が 13.1 秒から
8.8 秒に短縮することができた．
5.3. 内側ループの展開によるメモリアクセスの削減(1) チューニング方針
図 5-9 にチューニング前のコードを示す．最内側 n の DO ループの長さ(変数 nmax)はPARAMETER 文で 5 と宣言されており，ベクトル化の対象ループとなるが 121,122 行目の変数
vcx1，vcy2 はベクトル化の阻害要因となるため，コンパイラは最内側 n の DO ループでベクトル化を行わない．一つ外側 ic の DO ループの長さは 4 のため，コンパイラは更に外側 i の DO ループ

— 42 —

SENAC Vol. 48, No. 3（2015. 7）でベクトル化を行う．i の DO ループのベクトル化に伴い，ループ分割とループ入れ換えを行う．ループ分割では，i の DO ループのみに依存する演算と i，ic，n の DO ループに依存する演算に分割し，演算結果を作業配列に格納する．次に i，ic，n の DO ループの演算では，i の DO ループを内側に移動し事前に求めた作業配列の値を使用してベクトル処理する．このように外側の DO ル
ープをベクトル化することで作業配列を用いたループ分割とループ入れ換えを行うためメモリアク
セスの回数が増大する．図 5-10 の FTRACE 情報では，実行時間 267 秒の内バンクコンフリクト
時間が 196 秒を占めるため，メモリアクセス回数を減らしてバンクコンフリクト時間を短縮することを検討する．
!チューニング前
92: ||X---->
93: |||
94: |||&
95: |||
96: |||
97: |||
98: |||
99: |||
100: |||
101: |||
102: |||&
103: |||&
104: |||&
105: |||
106: |||
108: |||
109: |||
110: |||
111: |||
112: |||
113: |||
114: |||+--->
115: ||||
メモリストア116: ||||
117: ||||
118: ||||
119: ||||+-->
121: |||||
122: |||||
123: |||||
124: |||||
126: |||||
127: |||||
128: |||||
129: |||||
130: |||||
131: |||||
132: |||||
133: |||||
134: |||||
135: |||||
136: ||||+-:

do i = istart, iend+1
sde = sqrt( 0.25_DP*(dsx(i ,j-1,k) + dsx(i,j,k)+dsy(i-1,j ,k) + dsy(i,j,k)) )sdei = 1.0_DP/sde
vxe = 0.5_DP*( vxsx(i ,j-1,k) + vxsx(i,j,k) )vye = 0.5_DP*( vysy(i-1,j ,k) + vysy(i,j,k) )vaxe = 0.5_DP*( bx(i ,j-1,k) + bx(i,j,k) )*sdei
vaye = 0.5_DP*( by(i-1,j ,k) + by(i,j,k) )*sdei
cs_tmp = 0.25_DP*(cs2(i-1,j-1,k)+cs2(i ,j-1,k)+cs2(i-1,j ,k)+cs2(i ,j ,k))cs_tmp = sqrt(cs_tmp)
vcx0 = abs( vxe - vaxe )vcy0 = abs( vye + vaye )= 0.0_DPwc
vcx1c = 0.0_DP
vcy2c = 0.0_DP
do ic = 1,4
sign1 = sc(ic,1)
sign2 = sc(ic,2)
メモリロードvcx1 = vcx0*sign1
vcy2 = vcy0*sign2
do n = 1, nmax
vcx1old = vcx1
vcy2old = vcy2
if (vcx1 .eq. 0.0_DP) vcx1old = cs_tmp
if (vcy2 .eq. 0.0_DP) vcy2old = cs_tmp= i - nint( 0.5_DP + 0.5_DP*sign1 )i1delx = 0.5_DP*( dx – abs(vcx1*dt) )*sign1*fslp
vy2 = vysy(i1,j,k) + dvysydx(i1,j,k)*delxby2 =by(i1,j,k) +dbydx(i1,j,k)*delx
vcy2 = vy2 + by2*sdei
j2= j - nint( 0.5_DP + 0.5_DP*sign2 )dely = 0.5_DP*( dy – abs(vcy2*dt) )*sign2*fslp
vx1 = vxsx(i,j2,k) + dvxsxdy(i,j2,k)*delybx1 =bx(i,j2,k) +dbxdy(i,j2,k)*dely
vcx1 = vx1 - bx1*sdei
enddo

図 5-9 内側ループの展開によるメモリアクセス低減前のコード

— 43 —

ベクトルコンピュータにおける高速化

FREQUENCY50

EXCLUSIVE
TIME[sec](

AVER.TIME[msec]MOPS% )MFLOPS V.OP AVER.
RATIO V.LEN

VECTOR I-CACHE O-CACHEBANK CONFLICT
PROC.NAMETIMEMISSMISS CPU PORT NETWORK

266.912( 17.9)

5338.239

22400.9

6851.1 99.59 213.8

265.674

0.002

0.002

17.916178.408 【チューニング前】

図 5-10 内側ループの展開によるメモリアクセス低減前の FTRACE 情報(2) チューニング内容
内側の DO ループの長さがコンパイル時に判明しているため，EXPAND，UNROLL 指示行を挿入しループの展開を行う．内側ループを展開することで単純なベクトルループとなるため，ループ分割とループ入れ換えのための作業配列を使用する必要が無く，レジスタのデータで演算を行うことが可能になる．作業配列を使用しないためメモリアクセスの回数を削減することができる．図 5-11 にチューニング後コード，表 5-2 にチューニング前後のメモリアクセス回数の一覧を示す．
!チューニング後
do i = istart, iend+1
93: ||V---->
94: |||
sde = sqrt( 0.25_DP*(dsx(i ,j-1,k) + dsx(i,j,k)
95: |||&+
dsy(i-1,j ,k) + dsy(i,j,k)) )96: |||
sdei = 1.0_DP/sde
97: |||
vxe = 0.5_DP*( vxsx(i ,j-1,k) + vxsx(i,j,k) )98: |||
vye = 0.5_DP*( vysy(i-1,j ,k) + vysy(i,j,k) )99: |||
vaxe = 0.5_DP*( bx(i ,j-1,k) + bx(i,j,k) )*sdei
100: |||
vaye = 0.5_DP*( by(i-1,j ,k) + by(i,j,k) )*sdei
101: |||
102: |||
cs_tmp = 0.25_DP*(cs2(i-1,j-1,k)103: |||&+cs2(i ,j-1,k)104: |||&+cs2(i-1,j ,k)105: |||&+cs2(i ,j ,k))106: |||
cs_tmp = sqrt(cs_tmp)
107: |||
109: |||
vcx0 = abs( vxe - vaxe )110: |||
vcy0 = abs( vye + vaye )111: |||
wc= 0.0_DP
112: |||
vcx1c = 0.0_DP
113: |||
vcy2c = 0.0_DP
114: |||!cdir unroll=4115: |||
do ic = 1,4
116: |||*--->
117: ||||
sign1 = sc(ic,1)
118: ||||
sign2 = sc(ic,2)
119: ||||
vcx1 = vcx0*sign1
120: ||||
vcy2 = vcy0*sign2!cdir expand=nmax
121: ||||
122: ||||*-->
do n = 1, nmax
124: |||||
vcx1old = vcx1
125: |||||
vcy2old = vcy2
126: |||||
if (vcx1 .eq. 0.0_DP) vcx1old = cs_tmp
127: |||||
if (vcy2 .eq. 0.0_DP) vcy2old = cs_tmp
129: |||||
i1= i - nint( 0.5_DP + 0.5_DP*sign1 )130: |||||
delx = 0.5_DP*( dx – abs(vcx1*dt) )*sign1*fslp
131: |||||
vy2 = vysy(i1,j,k) + dvysydx(i1,j,k)*delx
132: |||||by2 =by(i1,j,k) +dbydx(i1,j,k)*delx
133: |||||
vcy2 = vy2 + by2*sdei
134: |||||
j2= j - nint( 0.5_DP + 0.5_DP*sign2 )
135: |||||
dely = 0.5_DP*( dy – abs(vcy2*dt) )*sign2*fslp
136: |||||
vx1 = vxsx(i,j2,k) + dvxsxdy(i,j2,k)*dely

— 44 —

SENAC Vol. 48, No. 3（2015. 7）

137: |||||
138: |||||
139: ||||*-:bx1 =bx(i,j2,k) +vcx1 = vx1 - bx1*sdei
enddo

dbxdy(i,j2,k)*dely

図 5-11 内側ループの展開によるメモリアクセス低減後のコード

表 5-2 内側ループの展開によるメモリアクセス低減前後のメモリアクセス命令数
チューニング前
471
329
160

メモリロードメモリストア
間接参照

チューニング後
332
40(3) 性能分析
図 5-12 にチューニング前後の性能値を示す．
FREQUENCY5050

EXCLUSIVE
TIME[sec](% )AVER.TIME[msec]MOPS

MFLOPS V.OP AVER.
RATIO V.LEN

266.912( 17.9) 5338.239 22400.9
6851.1 99.59 213.8
82.086( 8.5) 1641.720 46808.5 22759.7 99.46 214.4

VECTOR I-CACHE O-CACHEBANK CONFLICT
PROC.NAMETIMEMISSMISS CPU PORT NETWORK
265.674
70.789

0.002
3.235

0.002
0.002

17.9160.086

178.408 【チューニング前】
14.105 【チューニング後】

図 5-12 内側ループの展開によるメモリアクセス低減前後の FTRACE 情報
メモリアクセス回数が減少することでバンクコンフリクト時間が短縮され，実行時間が 266.9 秒から 82.1 秒に高速化することができた．
5.4. 複素数のバンク競合回避(1) チューニング方針
図 5-13 に FTRACE 情報を示す．FTRACE 情報では，実行時間に対してバンクコンフリクトの占める割合が大きいことが分かる．該当のサブルーチンでは複素数型配列を利用しているために飛びアクセスが発生していると考えられる．このため，複素数を実数部と虚数部に分け実数型配列に保存し，メモリ上の連続する領域へのアクセスとすることにより，速度低下やバンクコンフリクトを招く要因となりうる飛びアクセスの削減が期待できる．
FREQUENCY

44000

EXCLUSIVE
TIME[sec](

AVER.TIME[msec]

MOPS% )1453.723( 41.5)

33.03931198.3

MFLOPS V.OP AVER.
RATIO V.LEN

VECTOR I-CACHE O-CACHEBANK CONFLICT
PROC.NAMETIMEMISSMISS CPU PORT NETWORK

13481.2 99.62 255.9 1451.882

0.380

0.709

255.084

図 5-13 複素数のバンク競合回避前の FTRACE 情報

821.097 【チューニング前】ベクトルコンピュータにおける高速化!チューニング前
951: V------>
952: V-----:
955:
956: V------>
957: |
958: V-----:
991: V------>
992: |
993: V-----:
1009: V------>
1010: |
1011: V-----:
1028:
1029: V------>
1030: |
1031: V------

101— 45 —

DO 101 JG=1,NXYZ
RHO2(JG)=(0.D0,0.D0)*VDIR NODEP(RHO1)
DO 100 IG=1,NXYZ
JG=J2G(IG)100
RHO1(JG)=P(IG)do ig=1,nxyz
VG(ig)=VGG(ig)+Vloc(ig)
enddo

300

DO 300 I=1,NXYZ
fac=dt*dreal( vg(i) )RHO2(I)=dcmplx( dcos(fac),-dsin(fac) )*RHO1(I)*VDIR NODEP(RHO2)
DO 110 IG=1,NXYZ
JG=J2G(IG)110
P(IG)=RHO2(JG)

図 5-14 複素数のバンク競合回避前の前コード(2) チューニング内容
図 5-14 にチューニング前の複素数型配列を利用したコード，図 5-15 に最適化後の複素数を実数部と虚数部に分け実数型配列にしたチューニング後コードを示す．チューニング前コード
では倍精度複素数型を使用し配列 RHO1，RHO2 に演算を行っている．この該当箇所では倍精度実数型を使用し実数部は RHO1_R，RHO2_R とし，虚数部は RHO1_I ,RHO2_I にそれぞれ分け演算を行い処理するようにする．
!チューニング後
953: V------>
954: |
955: |
956: V-----:
968:
969: V------>
970: |
971: |
972: |
973: V-----:
984: V------>
985: |
986: |
987: |
988: |
989: V-----:
995:
996: V------>
997: |
998: |
999: |
1000: V------

DO

JG=1,NXYZ
RHO2_R(JG)=real(0.D0)
RHO2_I(JG)=aimag((0.D0,0.D0))
ENDDO*CDIR NODEP
DO

IG=1,NXYZ
JG=J2G(IG)
RHO1_R(JG)=dble(P(IG))
RHO1_I(JG)=aimag(P(IG))
ENDDO

do I=1,nxyz
VG_R(I)=VGG(I)+Vloc(I)fac=dt*VG_R(I)RHO2_R(I)=(dcos(fac)*RHO1_R(I))-((-dsin(fac))*RHO1_I(I))
RHO2_I(I)=(dcos(fac)*RHO1_I(I))-((dsin(fac))*RHO1_R(I))
enddo*CDIR NODEP
DO

IG=1,NXYZ
JG=J2G(IG)
P(IG)=dcmplx(RHO2_R(JG),RHO2_I(JG))
VG(IG)=dcmplx(VG_R(IG),0.0d0)
ENDDO

図 5-15 複素数のバンク競合回避後のコード

— 46 —

SENAC Vol. 48, No. 3（2015. 7）(3) 性能分析
図 5-16 は，チューニング前後の性能値を示す．
FREQUENCY EXCLUSIVE
AVER.TIME
TIME[sec]( % )[msec]
33.03944000
1453.723( 41.5)
44000
996.058( 34.9)
22.638MOPS
31198.345680.0

MFLOPS V.OP AVER.
VECTOR I-CACHE O-CACHEBANK CONFLICT
PROC.NAMERATIO V.LEN
TIME
MISS
MISS CPU PORT NETWORK
13481.2 99.62 255.9 1451.882
0.380
0.709 255.084 821.097 【チューニング前】
19675.5 99.63 255.9 993.704
0.676
1.039 234.856 466.970 【チューニング後】

図 5-16 複素数のバンク競合回避前後の FTRACE 情報
メモリ上の連続する領域へのアクセスすることにより，不連続なメモリアクセスが軽減され，実行
時間が 1,453 秒から 996 秒に短縮することができた．
6.

ベクトル演算の効率化の事例

6.1. 総和演算の効率化(1) チューニング方針
ベクトルループ内の総和演算は一回前の繰り返しの結果を参照するためベクトル化の阻害要因となる．しかし，コンパイラは総和演算を特別なパターンであることを認識し，専用のベクトル命
令(総和型マクロ演算)を用いることによりベクトル化を行う．図 6-1 はチューニング前のコードである．配列 Wn は i の次元を持たないため i の DO ループをベクトル化すると，コンパイラは総和型マクロ演算を用いてベクトル化を行う．この総和型マクロ演算により実行を高速化できるが，総和演
算を含む多重ループの場合，ループを入れ替え，総和型マクロ演算をベクトル演算にすることによりさらなる高速化が可能となる場合がある．このような事例を以下に示す．
!チューニング前
11: +------>
12: |V----->
13: ||
14: ||
15: ||
16: |V----17: +------

do j=1,ny
do i=1,nx
Wn(j) = Wn(j) + Wm(i,j)
end doend do

図 6-1 総和演算の効率化前のコード(2) チューニング内容
図 6-2 にループ入れ換え後のコードを示す．ループ入れ換えにより j の DO ループでベクトル化することで配列 Wm のメモリアクセスがストライドになることに注意されたい．
!チューニング後
12: +------>
13: |V----->
14: ||
15: ||
16: ||
17: |V----18: +------

do i=1,nx
do j=1,ny
Wn(j) = Wn(j) + Wm(i,j)
end doend do

図 6-2 総和演算の効率化後のコード

— 47 —

ベクトルコンピュータにおける高速化(3) 性能分析
図 6-3 にチューニング前後の性能値を示す．
FREQUENCY

1
1

EXCLUSIVE
TIME[sec](% )AVER.TIME[msec]MOPS

MFLOPS V.OP AVER.
RATIO V.LEN

232.086( 71.8)232086.438 17813.0
8620.3 99.30 173.6
91.148( 28.2) 91148.138 36054.1 20772.5 99.11 191.0

VECTOR I-CACHE O-CACHEBANK CONFLICT
PROC.NAMETIMEMISSMISS CPU PORT NETWORK
232.083
91.147

0.001
0.000

0.002
0.001

0.000
27.53632.327 【チューニング前】
4.160 【チューニング後】

図 6-3 総和演算の効率化前後の FTRACE 情報
チューニングを行うことで実効性能が向上し，実行時間が 232.1 秒から 91.1 秒に高速化することができた．
6.2. IF 文の括り出しによる演算数の削減(1) チューニング方針
ベクトルループ内に IF 文がある場合，SX-9 では IF 文の真偽に関わらず全ての演算を行い最後に IF 文の真にあたる結果のみを採用する方法でベクトル処理を行う． 図 6-4 にチューニング
前のコードを示す．
!チューニング前
12: +------>
13: |V----->
14: ||
15: ||
16: ||
17: |V----18: +------

do j=1,my
do i=1,mx
if(itbl(i)==1)thenWn(i,j) = Wn(i,j) + Wm(i,j)
end ifend doend do

図 6-4

IF 文の括り出しによる演算数削減前のコード(2) チューニング内容
図 6-5 にチューニング後のコードを示す．IF 文の判定式である配列 itbl はベクトル化対象である i の DO ループの次元しか持たないため，IF 文を含めたループ入れ換えを行うことでベクトルループ内に IF 文の処理が無くなる．ベクトル計算機では，ベクトルループ内に IF 文がある場合，条件式の真偽に関わらず両方の演算を行い，最後に条件式が真となる場合の結果を選択する．両方の演算を行うため演算数は増加するが，ベクトル演算を行うことにより性能が向上する．ループ
入れ換えを行うことにより，ベクトル演算の効率を低下させることなく演算数を削減し，実行時間を短縮することができる場合がある．このような高速化の事例を以下に示す．なお，ループの入れ換えにより，配列 Wn，Wm のメモリアクセスがストライドとなるため，メモリアクセスの観点から性能が低下する可能性があり，実際に高速化できるかについては確認が必要である．
!チューニング後
12: +------>
13: |
14: |V----->
15: ||
16: |V----17: |18: +------

do i=1,mx
if(itbl(i)==1)then
do j=1,myWn(i,j) = Wn(i,j) + Wm(i,j)
end doend ifend do

図 6-5

IF 文の括り出しによる演算数削減後のコード

— 48 —

SENAC Vol. 48, No. 3（2015. 7）(3) 性能分析
図 6-6 にチューニング前後の性能値を示す．
FREQUENCY

1
1

EXCLUSIVE
TIME[sec](

AVER.TIME[msec]MOPS

MFLOPS V.OP AVER.
RATIO V.LEN

VECTOR I-CACHE O-CACHEBANK CONFLICT
PROC.NAMETIMEMISSMISS CPU PORT NETWORK

98.349( 68.2) 98349.331
45.849( 31.8) 45848.683

36691.736026.1

20319.8 98.99 191.0
20756.3 99.11 191.0

98.34845.847

% )

図 6-6

0.000
0.000

0.001
0.001

22.75013.6834.486 【チューニング前】
2.121 【チューニング後】IF 文の括り出しによる演算数削減前後の FTRACE 情報

ベクトル演算量(実行時間×実効性能)が約 1/2 となるため，実行時間が 98.3 秒から 45.9 秒に短縮することができた．
6.3. ライブラリの置き換えによる高速化(1) チューニング方針
SX-9 では，HPC 用に高度に最適化された数学ライブラリ集の Mathkeisan が利用できる．
BLAS/LAPACK も最適化されており，プログラムの高速化が容易に行える．
図 6-7 にチューニング前のコードを示す．このコードでは，LAPACK のサブルーチン zheev が呼び出されている．図 6-8 のチューニング前の性能情報を確認すると，ベクトル化率が 98.17%となっている．zheev は QR 法を用いてエルミート行列の固有値・固有ベクトルを求めるが，分割統治
法でこれらを求めるライブラリ zheevd に置き換えることでベクトル化率の向上と計算時間の短縮を図る．
なお，分割統治法は QR 法より高速と知られているが，QR 法の演算量が行列の次数の 3 乗のオーダーであるのに対し，分割統治法の演算量は 2 乗から 3 乗のオーダーとばらつきがあることから，場合によっては，有意な性能差がみられないケースもありうる．また，必要となるメモリ領域の大きさについては，QR 法は次数の 1 乗のオーダーであるのに対し，分割統治法は 2 乗のオーダーとなる点も考慮する必要がある．
!チューニング前
983:
call zheev('V','U',nb*2,v2,nb*2,eig2,work2,2*(nb*2)-1,rwork2, info)

図 6-7 ライブラリ置き換え前のコード
FREQUENCY4EXCLUSIVE
TIME[sec](

AVER.TIME[msec]MOPS% )MFLOPS V.OP AVER.
RATIO V.LEN

VECTOR I-CACHE O-CACHEBANK CONFLICT
PROC.NAMETIMEMISSMISS CPU PORT NETWORK

78.111( 5.0)

19527.870

7546.64231.4 98.17 210.4

27.7660.022

23.0093.657

18.205【チューニング前】

図 6-8 ライブラリ置き換え前の FTRACE 情報(2) チューニング内容
図 6-9 にチューニング後のコードを示す．チューニング前のコードで呼び出されていた
LAPACK のサブルーチン zheev を zheevd に置き換えている．zheev と zheevd は呼び出しの際の引数が異なり，作業用の配列も 2 つ増える．作業配列のサイズは計算前に事前に一度 zheevd を呼び出すことで，計算に最適なサイズを取得することができる．取得したサイズの作業配列を確保
し，zheevd の呼び出す．

— 49 —

ベクトルコンピュータにおける高速化!チューニング後184:
lwork2= -1
185:
lrwork2= -1
186:
liwork2= -1
187:
call zheevd('V','U',nb*2,v2,nb*2,eig2,work2,lwork2,rwork2,
188:&
lrwork2,iwork2,liwork2,info)
189:
lwork2=work2(1)190:lrwork2=rwork2(1)191:
liwork2=iwork2(1)192:
deallocate (work2,rwork2,iwork2)
193:
allocate (work2(lwork2),rwork2(lrwork2),iwork2(liwork2)):
986:
call zheevd('V','U',nb*2,v2,nb*2,eig2,work2,lwork2,rwork2,
987:&
lrwork2,iwork2,liwork2,info)

図 6-9 ライブラリ置き換え後のコード(3) 性能分析
図 6-10 はチューニング後の性能情報である．適切なライブラリを選択することでベクトル化率
が 98.17％から 99.38％へ向上し，78.11 秒から 6.954 秒に短縮することができた．また，zheev からzheevd に置換することによる計算結果への影響もみられなかった．
FREQUENCY

44EXCLUSIVE
TIME[sec](

AVER.TIME[msec]MOPS% )MFLOPS V.OP AVER.
RATIO V.LEN

VECTOR I-CACHE O-CACHEBANK CONFLICT
PROC.NAMETIMEMISSMISS CPU PORT NETWORK

78.111(5.0)
6.954(0.5)

19527.870
1738.465

7546.653869

4231.4 98.17 210.4
31038.0 99.38 221.5

27.766
6.729

0.022
0.014

23.009
0.059

3.657
1.026

18.205【チューニング前】
2.812【チューニング後】

図 6-10 ライブラリ置き換え前後の FTRACE 情報
ここではライブラリルーチンの置き換えによる高速化の事例を挙げたが，置換する際には，置換によるメリット・デメリットを把握した上で置換を行い，置換後の計算結果の妥当性も必ず確認する必要がある．
6.4.

ファイルアクセスの高速化(1) チューニング方針
図 6-11 に FTRACE 情報を示す．FTRACE 情報では，ベクトル化率が 1.2％となっている．詳
細な解析により，図 6-12 に DO ループの中でファイルからデータを読み込んだ直後，読み込んだデータに対して演算が実行されていることが分かった．そこで，ファイルの読み込みと演算を別々に処理するよう最適化を行う．
FREQUENCY

1

EXCLUSIVE
TIME[sec](% )AVER.TIME[msec]

219.190( 12.5)219189.905

469.7

MOPS

4.1

MFLOPS V.OP AVER.
RATIO V.LEN
1.20 248.1

0.082

VECTOR I-CACHE O-CACHEBANK CONFLICT
PROC.NAMETIMEMISSMISS CPU PORT NETWORK
34.416

5.627

0.395

図 6-11 ファイルアクセス高速化前の FTRACE 情報

0.057 【チューニング前】— 50 —

!チューニング前
2776: +------>:|
2783: |
2784: |+----->
2785: ||
2786: ||+---->
2787: |||
2788: |||V--->
2789: ||||
2790: |||V--2791: ||+---2792: ||2793: ||+---->
2794: |||V--->
2795: ||||
2796: |||V--2797: ||+---2798: ||2799: ||
2800: ||
2801: ||
2802: ||
2803: |+----2804: |2805: |
2806: |
2807: |
2808: |+----->
2809: ||V---->
2810: |||
2811: ||V---2812: |+----2813: |2814: +------

SENAC Vol. 48, No. 3（2015. 7）

DO 450 IK = 1, NUMK
if ( my_rank.eq.0 ) thendo icpu=0,ncpu
nbleng=nend(icpu)-nbegin(icpu)+1
do 451 ib=1,nblengread(22)( rho1(ig),rho2(ig),ig=1,nxyz )do ig=1,nxyz
coef0(ig,ib,ik)=dcmplx( rho1(ig),rho2(ig) )enddo
451 continue
if ( icpu.eq.0 ) thendo ib=1,nbleng
do ig=1,nxyz
coef(ig,ib,ik)=coef0(ig,ib,ik)
enddo
enddo
else
call MPI_Send(coef0(1,1,ik),nxyz*nbleng,& MPI_DOUBLE_COMPLEX,icpu,tag,MPI_COMM_WORLD,ierr)endif
enddo ! end of icpu loop
else
nbleng=nend(my_rank)-nbegin(my_rank)+1
call MPI_Recv(coef0(1,1,ik),nxyz*nbleng& ,MPI_DOUBLE_COMPLEX,0,tag,MPI_COMM_WORLD,status,ierr)do ib=1,nbleng
do ig=1,nxyz
coef(ig,ib,ik)=coef0(ig,ib,ik)
enddo
enddo
endif
450 CONTINUE! ik roop

図 6-12 ファイルアクセスの高速化前のコード(2) チューニング内容
図 6-12 にチューニング前コード，図 6-13 にチューニング後コードを示す．チューニング前コードでは read(22) が DO ループの中にあり，その中でさらに演算を行っている．この該当箇所に対して read(22) を DO ループの外へ移動し，演算に関連する部分の配列を 1 次元配列から 2 次
元配列に変更しファイルの読み込みと演算を別々に処理するようにする．

— 51 —

ベクトルコンピュータにおける高速化!チューニング後
2776: +------>:|
2783: |
2784: |+----->
2785: ||
2786: |+----2787: |2788: |+----->
2789: ||
2790: ||+---->
2791: |||V--->
2792: ||||
2793: |||V--2794: ||+---2795: ||2796: ||
2797: ||+---->
2798: |||V--->
2799: ||||
2800: |||V--2801: ||+---2802: ||2803: ||
2804: ||
2805: ||
2806: ||
2807: ||
2808: |+----:|
2811: |
2812: |
2813: |
2814: |
2815: |+----->
2816: ||V---->
2817: |||
2818: ||V---2819: |+----2820: |
2821: +------DO 450 IK = 1, NUMK
if ( my_rank.eq.0 ) thendo ib=1,nend(ncpu)
read(22)(tmp_22(ig,ib),ig=1,nxyz)
enddo
icnt=0do icpu=0,ncpu
nbleng=nend(icpu)-nbegin(icpu)+1
do ib=1,nbleng
do ig=1,nxyz
coef0(ig,ib,ik)=tmp_22(ig,ib+icnt)
enddo
enddo
icnt=icnt+nbleng
if ( icpu.eq.0 ) thendo ib=1,nbleng
do ig=1,nxyz
coef(ig,ib,ik)=coef0(ig,ib,ik)
enddo
enddo
else
call MPI_Send(coef0(1,1,ik),nxyz*nbleng,
MPI_DOUBLE_COMPLEX,icpu,tag,MPI_COMM_WORLD,ierr)

&

endif
enddo ! end of icpu loop
else
nbleng=nend(my_rank)-nbegin(my_rank)+1
call MPI_Recv(coef0(1,1,ik),nxyz*nbleng,&
MPI_DOUBLE_COMPLEX,0,tag,MPI_COMM_WORLD,status,ierr)do ib=1,nbleng
do ig=1,nxyz
coef(ig,ib,ik)=coef0(ig,ib,ik)
enddo
enddo
endif
450 CONTINUE! ik roop

図 6-13 ファイルアクセスの高速化後のコード(3) 性能分析
図 6-14 にチューニング前後の性能値を示す．
FREQUENCY

1
1

EXCLUSIVE
TIME[sec](

219.190( 12.5)
0.095( 0.0)

% )AVER.TIME[msec]MOPS

MFLOPS V.OP AVER.
RATIO V.LEN

219189.905
94.761469.713163.1

4.1 1.20 248.1
0.0 98.99 248.0

VECTOR I-CACHE O-CACHEBANK CONFLICT
PROC.NAMETIMEMISSMISS CPU PORT NETWORK
0.082
0.087

34.416
0.001

5.627
0.001

0.395
0.017

0.057 【チューニング前】
0.061 【チューニング後】

図 6-14 ファイルアクセスの高速化前後の FTRACE 情報
ファイルの読み込みと演算を別々に処理することにより，ベクトル化率が 1.2％から 98.99％に改善したことで実行時間が 219.190 秒から 0.095 秒に短縮することができた．