— 25 —スーパーコンピュータ SX-9 の高速化

［大規模科学計算システム］高速化推進研究活動報告第 5 号より転載

スー�ー����ー� SX-9 の���
スーパーコンピューティング研究部

江川隆輔

岡部公起

情報部情報基盤課

伊藤英一

小野敏

山下毅

日本電気株式会社

撫佐昭裕

神山典

小久保達信

吉村健二

遠藤清隆

小沢実希

坂本英顕

金野浩伸

坂口祐太

NEC システムテクノロジー株式会社

曽我隆

4.1 SX-9 の特徴
SX-9 は高速ベクトルプロセッサと大規模主記憶装置を有するベクトル型スーパーコンピュータである．東北大学サイバーサイエンスセンターでは 18 ノードの SX-9 を導入しており，そのうち 16 ノードがノード間接続装置（以下，IXS）に接続されたマルチノードシステムである．表 4.1-1 に SX-9 の主要諸元を，表 4.1-2に SX-9 マルチノードシステムの主要諸元を示す．
表 4.1-1 SX-9 主要諸元
項目
最大ベクトル演算性能
CPU 数

CPU

ベクトルレジスタ
ベクトルマスクレジスタ
レジスタ
スカラレジスタ
固定小数点
浮動小数点
データ形式
*128bit はスカラ命令のみサポート
論理
ベクトル演算パイプライン
ADB
スカラ演算パイプライン
キャッシュ

主記憶装置

容量
最大データ転送能力

ノード数
最大ベクトル演算性能
CPU 数
主記憶装置容量IXS 転送性能（片方向）諸元
1.6TFLOPS
16
144KB
256bit×16
64bit×128
32/64bit
32/64/128bit
IEEE
64bit
5 種類×8 セット
256KB
1 セット
命令:32KB
オペランド:32KB
1TB
4TB/s

表 4.1.2 SX-9 マルチノードシステム主要諸元
項目
諸元
16
26.2TFLOPS
256
16TB
128GB/s (ノード当たり)
2TB/s (システム全体)(1) CPU の特徴
SX-9 に搭載されている CPU は，単一チップで 102.4GFLOPS のベクトル演算性能を有するベクトル

— 26 —

SENAC Vol. 45, No. 2（2012.4）プロセッサである．図 4.1-1 に SX-9 のプロセッサ構成を示す．CPU 内部はベクトルユニットとスカラユ
ニットに分かれており，ベクトルユニットは 8 つのベクトルパイプラインセットで構成されている．各ベクトルパイプラインセットはマスク演算器，論理演算器，乗算器×2，加算器×2，除算器/平方根演算器を有している．また，SX-9 から MAX/MIN 関数がハードウェア命令化され，最大値・最小値の計算時間が短縮されている．

図 4.1-1 SX-9 プロセッサの構成図(2) メモリシステムの特徴
SX-9 は，CPU と主記憶装置（MMU）間の距離が論理的に等しい共有メモリ型計算機である．図 4.1-2 に SX-9 のメモリ構成を示す．SX-9 は CPU と主記憶装置間にネットワークルータ（以下，RTR）という
内部スイッチング機構を配することで，広いメモリバンド幅と大規模共有メモリを実現している．また，主記憶装置はノードあたり 32,768 個のメモリバンクで構成されている．バンクのアクセスに関しては SX-7
C までの SX シリーズでは一度に 1 演算要素（8Byte）を処理していたが，SX-9 から 2 演算要素（16Byte）を処理するように強化し，ロード・ストアの処理時間の短縮を図っている．

図 4.1-2 SX-9 のメモリ構成図

— 27 —

スーパーコンピュータ SX-9 の高速化

SX-9 では CPU－主記憶装置間に ADB（Assignable Data Buffer）と呼ばれるベクトルデータを選択的にバッファする機構を有している．図 4.1-3 に ADB の概念図を示す．ADB は CPU－主記憶装置間よりメモリバンド幅が広く，メモリレイテンシ（遅延）が短いという特徴を持っている．表 4.1-3 に SX-9 の主記憶装置と ADB のデータ供給性能を示す．データ供給性能は，演算性能あたりのデータ供給量を表す指標としてベクトル演算性能あたりのメモリバンド幅（Byte/FLOP）である．

図 4.1-3 ADB の概念図

表 4.1-3 SX-9 の主記憶装置と ADB のデータ供給性能
CPU-主記憶装置間
バンド幅256.0GB/s
ベクトル演算性能あたりのバンド幅
2.5Byte/FLOP

CPU-ADB 間
409.6GB/s
4.0Byte/FLOP(3) ソフトウェアの特徴①

コンパイラ
SX-9 用コンパイラとして FORTRAN90/SX と C++/SX が利用可能である．FORTRAN90/SX は
ISO/IEC 1539-1:1997，ISO/IEC1539:1991 に準拠している．C++/SX は ISO/IEC 9899:1999，

ISO/IEC 14882:1998 に準拠している．
FORTRAN90/SX，C++/SX は自動ベクトル化機能，自動並列化機能を有しており，ユーザは意識することなく，ベクトル化，並列化することができる．
②

MPI ライブラリ
SX-9 用 MPI ライブラリは MPI-1.2 および MPI-2 に準拠した機能を提供している．

③

科学技術計算ライブラリ
SX-9 用ライブラリとして，科学技術計算ライブラリ ASL，数値計算ライブラリ MathKeisan が利用可能である．

4.2 ベクトル処理による高速化(1) ベクトル処理の概要
科学技術計算プログラムの多くは特定の DO ループに実行の大部分が集中し，DO ループ内の配— 28 —

SENAC Vol. 45, No. 2（2012.4）

列データ（ベクトル）に対し繰り返し演算が行われている．この規則性に着目して高速化を図った方式
がベクトル処理である．ベクトル処理は DO ループによる配列データへの繰り返し演算を，ベクトルユニ
ットを用いて一括に実行するものである．一方，スカラ処理は繰り返し演算を逐次的に実行するものである．
(2) ベクトル化率
ベクトル型スーパーコンピュータではプログラム中のベクトル処理可能な部分を高速に実行している．
図 4.2-1 に同一のプログラムをスカラ処理する場合とベクトル処理する場合の実行時間に関する概念
図を示す．
T1（スカラ処理で実行した時の総実行時間）
T1V（T1のうちベクトル処理で実行可能な時間）

スカラ処理の場合

ベクトル処理不可部分

ベクトル処理の場合

ベクトル処理不可部分

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

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

図 4.2-1 ベクトル処理による実行時間短縮のイメージ

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

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

P=

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

— 29 —

スーパーコンピュータ SX-9 の高速化

β

性能向上比
（
倍）

P

0
0

20

40

�0

�0

�00

ベクトル化率α （％）

図 4.2-2 ベクトル化率と性能向上比（アムダールの法則）
ベクトル化率を求めるため，プログラムをスカラ実行する場合とベクトル実行する場合のそれぞれの実行時間が必要である．しかし，これらの実行時間は同時に得られないのでベクトル化率の算出は困難である．そのため SX-9 ではベクトル化率の代わりにベクトル演算率をベクトル化の指標としている．
後述の PROGINF や FTRACE などの性能解析機能に出力される情報はベクトル演算率である．ベクト
ル演算率は，プログラムで処理される全演算要素数に対するベクトル演算命令で処理される演算要素数の割合で算出し，ほぼベクトル化率とみなせる指標である．
(3) ベクトル長
ベクトル処理を効率的に行う上で重要な指標にベクトル長がある．ベクトル長はベクトル化対象のDO ループの繰り返し回数のことであり，効率良くベクトル処理を行うためには，できるだけループ長を長くする必要がある．
一般に命令を発行してから結果が返ってくるまでに遅延時間が存在する．この遅延時間を立ち上がり時間と呼ぶ．図 4.2-3 にベクトル処理における立ち上がり時間および演算時間の概念図を示す．図4.2-4 にベクトル長と立ち上がり時間の関係を示す．立ち上がり時間はベクトル長が異なる場合でも一定である．よって総演算量が等しい場合，短ベクトル長処理を繰り返すよりも長ベクトル長処理を一括して行う方が実行時間を短くすることができる．このように高速化を図るためには，DO ループのベクト
ル長を十分に確保し，演算時間に対する立ち上がり時間の占める割合を低くすることが重要である．

— 30 —

SENAC Vol. 45, No. 2（2012.4）【短ベクトル長】
C

C

C

C

演算器

C：演算時間（ベクトル長に比例）

L

L

L

L

【長ベクトル長】
4×C
演算器

L：立ち上がり時間
L

図 4.2-3 立ち上がり時間    

図 4.2-4 ベクトル長と立ち上がり時間(4) メモリアクセス
ベクトルユニットはベクトル処理によって複数の演算要素を一括で処理することが可能である．この処理を高速に行うためには，主記憶装置からベクトルユニットに対して，データを滞りなく供給することが必要である．
データ供給の遅延の要因にバンクコンフリクトがある．バンクコンフリクトは CPU ポート競合とメモリネットワーク競合の 2 つに分類される．CPU ポート競合は CPU 内における同一のポートにロード・ストアが集中した時に発生する競合を指す．メモリネットワーク競合は同一のメモリバンクへのアクセスで発生する競合や CPU と主記憶装置間の経路上で発生する競合を指す．バンクコンフリクトが発生するとデー
タ供給能力が低下し，高速化の妨げの要因になる．バンクコンフリクトを回避するためには，メモリアク
セスの間隔が偶数となることを避け，可能ならば連続アクセスとすることや，ループのアンロールを行うなどしてメモリアクセス回数を減らすことが重要である．
また，メモリアクセスの効率化では，前述の ADB の活用があげられる．ADB を利用することで主記憶装置からのデータ供給性能を補うことや，メモリレイテンシを短くすることが可能となり性能向上を図れる場合がある．後述のチューニング例 4.3.8(4)を参考に ADB の利用を試みて頂きたい．

4.3 高速化技法
4.3.1 多重 DO ループの一重化による高速化(1) 多重 DO ループの一重化の概要
多重 DO ループの一重化とは，図 4.3-1 のような多重構造の DO ループを，図 4.3-2 のように一重の DO ループに書き換えることをいう．図 4.3-1 では，内側の DO ループはベクトル長 10 でベクトル処理され，外側の DO ループはこのベクトル処理が 20 回実行されるので，ベクトル処理の立ち上がり時
間は 20 回分必要となる（4.2.(3) 図 4.2-4 参照）．一方，図 4.3-2 では，ベクトル長 200 のベクトル処理を１回だけ実行すればよいので，立ち上がり時間は１回分だけですみ，高速化することができる．ここで，配列 a(ij,1)，b(ij,1)，c(ij,1)は一重化された DO ループで，元の配列 a(i,j)，b(i,j),c(i,j)を連続してアクセスするために書きかえたものである．
多重 DO ループの一重化は，ここで示した例のように，内側の DO ループの長さが短いときに特に効果がある．

— 31 —

スーパーコンピュータ SX-9 の高速化

!多重 DO ループ
real,dimension(10,20)::a,b,c
do j=1,20
do i=1,10a(i,j)=b(i,j)+c(i,j)
enddo
enddo

図 4.3-1 多重 DO ループの例!図 4.3-1 の多重 DO ループの一重化
do ij = 1, 10*20a(ij,1) = b(ij,1) + c(ij,1)
enddo

図 4.3-2 多重 DO ループの一重化の例(2) コンパイラによる DO ループの一重化
多重 DO ループは，コンパイラによって自動的に一重化される場合がある．4.3.1.(1)であげた図
4.3-1 から図 4.3-2 への変形は，コンパイラが自動的に行うものである．コンパイラが変換可能なのは，多重ループが以下の条件を満たす場合である．
①

密な多重 DO ループ（注）である②

指標変数を含む添字式が，その配列の次元の下限から上限までの値をとる③

各 DO ループの指標変数は，ループ中に現れる配列要素の添字中に左から連続して同じ形式
かつ同じ順序で現れている④

DO ループ中に現れる配列要素の各添字式は異なるループの指標変数を含まない⑤

一重化しても定義・引用の順序関係が正しく保たれる⑥

ループ中に現れるすべての配列はポインタでも形状引継ぎ配列でもない（注）密な多重 DO ループとは，直接の入れ子関係をなす DO ループのうち，外側ループの DO 文と内側ループの DO 文の間および内側ループの ENDDO 文と外側ループの ENDDO 文の間に，実行文が現れないもの（およびそれと等価な構造の多重ループ）．
(3) 自動的に一重化されない多重 DO ループのチューニング
DO ループ中にポインタや形状引継ぎ配列がある場合は上述⑥の条件を満たさないため，コンパイ
ラによる自動的なループの一重化は行われない．このような場合，ループはコンパイル時にオプション「-pvctl collapse」（注 1）を指定することにより一重化することができる．
この他，自動的に一重化されない例として，図 4.3-3 のように DO ループの範囲の外に冗長な領域を持つ配列が含まれる場合をあげられる．ループの指標変数が配列 work の下限から上限までの値をとらず，上述②の条件を満たさないため，ループの一重化は行われない．図 4.3-4 のように配列 workの冗長な領域をなくすことで，コンパイラにより自動的に一重化されるようになる．なお，図 4.3-3 と図
4.3-4 は，コンパイラによる編集リスト（注 2）の形式で記載している．

— 32 —

SENAC Vol. 45, No. 2（2012.4）!オリジナル
42:
43:
44:
45:
46:
47:
48:
49:
50:
51: +------>
52: |+----->
53: ||V---->
54: |||
55: |||
56: |||
57: ||V---58: |+----59: +------

subroutine sample_f(a,b,work,p)implicit none
include 'parameter.h'real*8,intent(in) :: a(imax,jmax,kmax,3)
real*8,intent(in) :: b(imax,jmax,kmax,3)
real*8,intent(in) :: work(imax+1,jmax+1,kmax+1,3)
real*8,intent(out) :: p(imax,jmax,kmax)
integer :: i,j,k
do k=1,kmax
do j=1,jmax
do i=1,imax
p(i,j,k)=work(i,j,k,1)*(a(i,j,k,1)+b(i,j,k,1)) &+work(i,j,k,2)*(a(i,j,k,2)+b(i,j,k,2)) &+work(i,j,k,3)*(a(i,j,k,3)+b(i,j,k,3))
enddo
enddo
enddo

図 4.3-3 多重 DO ループが一重化されないコードの例!チューニング
42:
43:
44:
45:
46:
47:
48:
49:
50:
51: W------>
52: |*----->
53: ||*---->
54: |||
55: |||
56: |||
57: ||*---58: |*----59: W------

subroutine sample_f(a,b,work,p)implicit none
include 'parameter.h'real*8,intent(in) :: a(imax,jmax,kmax,3)
real*8,intent(in) :: b(imax,jmax,kmax,3)
real*8,intent(in) :: work(imax,jmax,kmax,3)
real*8,intent(out) :: p(imax,jmax,kmax)
integer :: i,j,k
do k=1,kmax
do j=1,jmax
do i=1,imax
p(i,j,k)=work(i,j,k,1)*(a(i,j,k,1)+b(i,j,k,1)) &+work(i,j,k,2)*(a(i,j,k,2)+b(i,j,k,2)) &+work(i,j,k,3)*(a(i,j,k,3)+b(i,j,k,3))
enddo
enddo
enddo

図 4.3-4 多重 DO ループが一重化されないコードのチューニング例（注 1）「-pvctl collapse」は，次元数が 2 以上の配列式，または多重 DO ループ中に現われる形状引き継ぎ配列，ポインタ配列の全要素がメモリ上の連続領域にあると仮定し，一重化を行うことを指示するコンパイルオプションである．
（注 2）編集リストとは，ソースコードと共に，ループと配列式のベクトル化情報，手続きのインライン展
開情報，ループと配列式の並列化情報を出力する機能である．以下に編集リストに用いられる記号の意味を示す．
V：ベクトル化されたループ・配列式
W：一重化されてベクトル化されたループ・配列式
P：並列化されたループ
Y：並列化され，かつベクトル化されたループ
I：インライン展開された手続き+：最適化されなかったループ

— 33 —

スーパーコンピュータ SX-9 の高速化

(4) 性能評価(3)の図 4.3-3，図 4.3-4 で示したチューニング前後のコードについて，問題サイズを「imax= 50,
jmax=200, kmax=1000」として，SX-9 1CPU を用いて性能評価を行う．SX-9 では，性能情報として
PROGINF（注 1）と FTRACE（注 2）がある．PROGINF からはプログラム全体の性能情報が取得でき，
FTRACE からは手続（サブルーチンや関数）ごとの性能情報を取得することができる．
①

PROGINF
図 4.3-5，図 4.3-6 はチューニング前後の PROGINF である．実効性能を比較するため，
MFLOPS（1 秒間に実行された浮動小数点演算数を 100 万単位で示した値）の項目を見る．チューニング前は 8,578MFLOPS であるのに対し，チューニング後には 20,063MFLOPS となっており，実効性能が 2.33 倍向上していることが分かる．また，A.V.Length（平均ベクトル長（注 3））の項目を見ると，チューニング前は 49.9 であるのに対し，チューニング後は 252.3 となっている．チューニ
ングにより多重 DO ループが一重化されることで，コード全体の平均ベクトル長が長くなっていることを確認できる．
****** Program Information ******
Real Time (sec):User Time (sec):
Sys Time (sec):
Vector Time (sec):Inst. Count:
V. Inst. Count:
V. Element Count:FLOP Count
:
MOPS
:
:
MFLOPS:
�. V. Length
V. Op. Ratio (%):
Memory Size (MB):MIPS:
I-Cache (sec):
O-Cache (sec):Bank Conflict TimeCPU Port Conf. (sec):Memory Network Conf. (sec) :4.699699
4.672104
0.023863
4.667600
3615133579.
1806406030.
90320127249.
40080201114.
19718.922095
8578.619208
49.999904
98.036742
960.031250
773.769929
0.000384
0.000478
0.053537
3.597606

図 4.3-5 チューニング前の PROGINF****** Program Information ******
Real Time (sec):User Time (sec):
Sys Time (sec):
Vector Time (sec):Inst. Count:
V. Inst. Count:
V. Element Count:FLOP Count
:
MOPS
:
:
MFLOPS:
�. V. Length
V. Op. Ratio (%):
Memory Size (MB):MIPS:
I-Cache (sec):
O-Cache (sec):Bank Conflict TimeCPU Port Conf. (sec):Memory Network Conf. (sec) :2.028015
1.997679
0.025952
1.993263
695335689.
357973030.
90320127249.
40080201166.
45381.410080
20063.384140
252.309866
99.627871
960.031250
348.071782
0.000325
0.000416
0.010343
0.779060

図 4.3-6 チューニング後の PROGINF

— 34 —②

SENAC Vol. 45, No. 2（2012.4）

FTRACE
図 4.3-7，図 4.3-8 はチューニング前後の FTRACE の情報である．①PROGINF と同様に，実効性能と平均ベクトル長を確認することができるが，FTRACE ではこれらの情報を手続（サブルーチンや関数）単位で取得することができる．
PROC.NAMEFREQUENCY

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

sample_f
500
4.638( 99.5)
9.276 19794.6
8624.8 98.04 50.0
4.637
0.000
0.000
0.051
3.577
tp_sample
1
0.022( 0.5)
22.060 14752.3
3635.5 98.33 50.0
0.022
0.000
0.000
0.000
0.015
------------------------------------------------------------------------------------------------------------------------total
501
4.660(100.0)
9.301 19770.7
8601.2 98.04 50.0
4.659
0.000
0.000
0.051
3.593

図 4.3-7 チューニング前の FTRACE 情報
PROC.NAMEFREQUENCY

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

sample_f
500
1.959( 98.9)
3.918 46114.9 20420.2 99.63 256.0
1.959
0.000
0.000
0.011
0.715
tp_sample
1
0.022( 1.1)
22.036 14749.9
3639.4 98.45 50.0
0.022
0.000
0.000
0.000
0.015
------------------------------------------------------------------------------------------------------------------------total
501
1.981(100.0)
3.954 45765.9 20233.5 99.63 252.3
1.980
0.000
0.000
0.011
0.730

図 4.3-8 チューニング後の FTRACE 情報（注 1）プログラム実行解析情報：プログラム全体の性能情報が標準エラー出力ファイルに出力される．以下に PROGINF 機能で取得可能な情報を示す．
Real Time：経過時間（プログラム実行に要した時間）

User Time：ユーザ時間（プログラム実行に要した CPU 時間の内，ユーザルーチンが実行に要した時間）

Sys Time：システム時間（プログラム実行に要した CPU 時間の内，システムルーチンが実行に要した時間）

Vector Time：ベクトル命令実行時間

Inst. Count：全命令実行数

V. Inst. Count：ベクトル命令実行数

V. Element Count：ベクトル命令で実行された演算要素の個数

FLOP Count：実行された浮動小数点演算の個数

MOPS：1 秒間に実行された演算数を 100 万単位で示した値

MFLOPS：1 秒間に実行された浮動小数点演算数を 100 万単位
で示した値

A. V. Length：平均ベクトル長（注 3 参照）

V. Op. Ratio：全実行命令の個数の内，ベクトル命令の占める割合

Memory Size

：メモリ使用量

MIPS：1 秒間に実行された命令数を 100 万単位で示した値

I-Cache：命令キャッシュミスにより発生した合計時間

O-Cache：オペランドキャッシュミスにより発生した合計時間

Bank Conflict Time：バンクコンフリクト時間

CPU Port Conf.

：CPU ポート競合時間

Memory Network Conf.

：メモリネットワーク競合時間

— 35 —

スーパーコンピュータ SX-9 の高速化

（注 2）性能解析情報：コンパイル時にオプション「-ftrace」を指定することで，手続（サブルーチンや関数）ごとの性能解析情報を採取することができ，プログラム実行後に解析情報が標準出力
ファイルに出力される．特に，チューニング対象とする手続を選択する際に活用すると良い．
以下に FTRACE 機能で取得可能な情報を示す．
PROC.NAME：手続名

FREQUENCY：手続の呼び出し回数

EXCLUSIVE TIME

：手続の実行に要した占有の CPU 時間と，手続全体の実行に要した CPU 時間に対する比率

AVER.TIME：手続の 1 回の実行に要した平均 CPU 時間

MOPS：1 秒間に実行された演算数を 100 万単位で示した値

MFLOPS：1 秒間に実行された浮動小数点演算数を 100 万単位
で示した値

V.OP RATIO：ベクトル演算率

AVER V.LEN：平均ベクトル長（注 3 参照）

VECTOR TIME

：ベクトル命令実行時間

I-CACHE MISS

：命令キャッシュミスにより発生した合計時間

O-CACHE MISS

：オペランドキャッシュミスにより発生した合計時間

BANK CONFLICT：バンクコンフリクト時間

CPU PORT： CPU ポート競合時間

NETWORK： メモリネットワーク競合時間（注 3）SX-9 の性能情報の一つに，平均ベクトル長がある．SX-9 は一つのベクトル命令で最大 256組の演算要素を処理することができるので，例えば DO ループの長さが 1000 の場合，256，
256，256，232 という 4 つの処理単位に DO ループを分割して処理する．この場合，平均ベクトル長は 250 になる．

4.3.2 多重 DO ループのアンローリングによる高速化(1) 多重 DO ループのアンローリングの概要
DO ループのアンローリングとは，ループ本体の演算を n 倍にし，繰り返し数を 1/n にする処理である．このとき n をアンローリングの段数という．
図 4.3-9 の外側の j のループについて，2 段のアンローリングをしたのが図 4.3-10 である．外側の jのループは増分 2 とし，内側の i のループ内に j+1 の演算を追加して，内側ループの演算を 2 倍にしている．これにより，図 4.3-10 の①の x(i)のストアと②の x(i)のロードが不要になる．つまり，図 4.3-10は図 4.3-11 のように変形できる．よって，多重 DO ループ全体で考えると，x(i)のロード回数，ストア回
数がアンローリング前の半分になる．
多重 DO ループのアンローリングでは，一般に，ベクトル化されない外側ループについてアンローリ
ングを行うことで性能向上が期待できる．特にロード・ストア回数が減少するような場合には，大きな効果が得られることが多い．

— 36 —

SENAC Vol. 45, No. 2（2012.4）!多重 DO ループ
do j=1,100
do i=1,nx(i)=x(i)+a(i,j)*b(j)
enddo
enddo

図 4.3-9 多重 DO ループの例!図 4.3-9 の多重 DO ループを，外側ループについて 2 段でアンローリング
do j=1,99,2
do i=1,nx(i)=x(i)+a(i,j)*b(j)・・・①x(i)=x(i)+a(i,j+1)*b(j+1) ・・・②
enddo
enddo

図 4.3-10 外側 DO ループのアンローリング後!図 4.3-10 の x(i)のロード・ストアを削減した結果
do j=1,99,2
do i=1,nx(i)=x(i)+a(i,j)*b(j)+a(i,j+1)*b(j+1)
enddo
enddo

図 4.3-11 x(i)のロード・ストア削減の結果(2) 差分コードのアンローリング
図 4.3-12 は，差分コードのチューニング例である．多重 DO ループにおいてコンパイラが外側ループのアンローリングを行うように，コンパイラ指示行 OUTERUNROLL（注 1）を挿入している．
!チューニング：
47: +------>
48: |
49: |+----->
50: ||V---->
51: |||
52: |||
53: |||
54: |||
55: |||
56: |||
57: |||
58: |||
59: ||V---60: |+----61: +-----：

do k=2,kmax-1!�D�� O�������O��=4
do j=2,jmax-1
do i=2,imax-1
s0=c(i,j,k,1)*(a(i+1,j+1,k)-a(i+1,j-1,k) &-a(i-1,j+1,k)+a(i-1,j-1,k)) &+c(i,j,k,2)*(a(i,j+1,k+1)-a(i,j-1,k+1) &-a(i,j+1,k-1)+a(i,j-1,k-1)) &+c(i,j,k,3)*(a(i+1,j,k+1)-a(i-1,j,k+1) &-a(i+1,j,k-1)+a(i-1,j,k-1))
ss=s0*b(i,j,k)p(i,j,k)=a(i,j,k)+omega*ss
enddo
enddo
enddo

図 4.3-12 差分コードのチューニング例
図 4.3-13 は図 4.3-12 のチューニング例について，アンローリング部分の変形結果を変形リスト（注2）より抜粋したものである．

スーパーコンピュータ SX-9 の高速化

— 37 —!図 4.3-12 のアンローリング部分の変形リスト
2 回参照される配列の例
1
do j = 3, 254, 4
2 !cdir
nodep
3
do i = 1, 2544
s01 = c(1+i,j+1,k,1)*(a(2+i,j+2,k) -a(2+i,j,k)-a(i,j+2,k)+a(i,j,k))5+ c(1+i,j+1,k,2)*(a(1+i,j+2,k+1)-a(1+i,j,k+1) -a(1+i,j+2,k-1)+a(1+i,j,k-1))
6+ c(1+i,j+1,k,3)*(a(2+i,j+1,k+1)-a(i,j+1,k+1) -a(2+i,j+1,k-1)+a(i,j+1,k-1))
7
s02 = c(1+i,j+2,k,1)*(a(2+i,j+3,k) -a(2+i,j+1,k) - a(i,j+3,k)+a(i,j+1,k))8+ c(1+i,j+2,k,2)*(a(1+i,j+3,k+1)-a(1+i,j+1,k+1)-a(1+i,j+3,k-1)+a(1+i,j+1,k-1))
9+ c(1+i,j+2,k,3)*(a(2+i,j+2,k+1)-a(i,j+2,k+1) -a(2+i,j+2,k-1)+a(i,j+2,k-1))
10
s03 = c(1+i,j+3,k,1)*(a(2+i,j+4,k) -a(2+i,j+2,k) -a(i,j+4,k)+a(i,j+2,k))11+ c(1+i,j+3,k,2)*(a(1+i,j+4,k+1)-a(1+i,j+2,k+1)-a(1+i,j+4,k-1)+ a(1+i,j+2,k-1))
12+ c(1+i,j+3,k,3)*(a(2+i,j+3,k+1)-a(i,j+3,k+1) -a(2+i,j+3,k-1)+a(i,j+3,k-1))
13
s04 = c(1+i,j+4,k,1)*(a(2+i,j+5,k) -a(2+i,j+3,k) -a(i,j+5,k)+ a(i,j+3,k) )14+ c(1+i,j+4,k,2)*(a(1+i,j+5,k+1)-a(1+i,j+3,k+1)-a(1+i,j+5,k-1)+a(1+i,j+3,k-1))
15+ c(1+i,j+4,k,3)*(a(2+i,j+4,k+1)-a(i,j+4,k+1) -a(2+i,j+4,k-1)+a(i,j+4,k-1))
16p(1+i,j+1,k) = a(1+i,j+1,k) + 8.00000000000000e-001*s01*b(1+i,j+1,k)17p(1+i,j+2,k) = a(1+i,j+2,k) + 8.00000000000000e-001*s02*b(1+i,j+2,k)18
p(1+i,j+3,k) = a(1+i,j+3,k) + 8.00000000000000e-001*s03*b(1+i,j+3,k)19
p(1+i,j+4,k) = a(1+i,j+4,k) + 8.00000000000000e-001*s04*b(1+i,j+4,k)
20
enddo
21
enddo

図 4.3-13 アンローリング後の変形リスト
この 4 段のアンローリングの例では，配列 a, b, c は全部で 68 個の要素がロードされるが，そのうち
配列 a は 8 個の要素が 2 回参照される（注 3）．重複する部分は再度ロードされることが無いため，DO
ループ内のロード回数が 68 回から 60 回に減る．同様に， 8 段でアンローリングした場合は，配列 a, b,
c は全部で 136 個の要素がロードされるが，そのうち配列 a は 24 個の要素が 2 回参照されるため，DO
ループ内のロード回数が 136 回から 112 回に減る．
（注 1）直後の外側 DO ループに対し，強制的に n 段のアンローリングを行うことを指定する指示行である．これを指定した場合，コンパイラはループ内のデータの依存関係をチェックせずにアンローリングを行う．なお OUTERUNROLL の段数は，指定した値以下の 2 のべき乗の値となる．
（注 2）変形リスト：コンパイラが行った最適化の結果，実際に生成されたオブジェクトコードの構造を
ソースコードのイメージで表示したリストである．
（注 3）図 4.3-13 では，変形リストを抜き出して行番号を追記し，2 回参照される配列がわかりやすい
よう，配列ごとに形式を変えて強調表記している．例えば，行番号 4 と 10 の矢印で示した
a(i,j+2,k)の配列は a(i,j+2,k)と表記している．
(3) 性能評価(2)の図 4.3-12，図 4.3-13 で示しているチューニング前後のコードについて，問題サイズを「imax=256, jmax=256, kmax=256」とし，SX-9 1CPU を用いて性能評価を行う．表 4.3-1 が性能測定結果である．4 段でアンローリングする場合には，チューニング前に比べて 1.16 倍，8 段でアンローリングする場合には，1.30 倍性能向上していることが分かる．
表 4.3-1 アンローリングの性能比較
性能
チューニング前
16.29GFLOPS
アンローリング後（4 段）
18.91GFLOPS
アンローリング後（8 段）
20.21GFLOPS

— 38 —

SENAC Vol. 45, No. 2（2012.4）(4) アンローリング段数による性能の違い(3)での評価結果からもわかるとおり，アンローリング段数によって性能に差が生じている．図 4.3-14
は，(2)の差分コードについて，段数を 2, 4, 8, 16, 32 と変えて性能測定している結果である．

25

GFLPOS

20
15
10
5
0
無し

2

4
8
16
アンローリング段数

3264

図 4.3-14 アンローリング段数による性能差
一般に，アンローリング段数を増やすと性能が向上する．しかし，段数が増えるに伴いループ中の中間変数が増加する．そのため段数を増やしすぎると途中で中間変数を保存するレジスタが足りなくなり，性能が劣化する．図 4.3-14 のとおり，(2)の差分コードの例では，16 段でアンローリングすると最も性能が高くなる．最適なアンローリング段数はコードにより異なる．
(5) アンローリングの自動/手動に関する補足
外側の DO ループをアンローリングすることで，内側の DO ループのベクトル演算パイプラインを有効に利用できる場合や，ロード・ストア回数を減らす効果がある場合，コンパイラは自動的にアンローリ
ングを行う．コンパイラによって自動的にアンローリングされる場合でも，デフォルトの段数（通常 4 段）
以外で実行すると更に性能向上が得られる場合がある．多重ループでは，ロード・ストア回数が最も削減される次元でアンロールするのが有効である．

4.3.3 IF 文を含む DO ループのチューニング
(1) IF 文を含む DO ループのベクトル処理の概要
図 4.3-15 のように IF 文を含む DO ループは，つぎのように処理される．
①

条件式「flag(n).eq.1」の真偽に応じて，マスクベクトルを作成②

条件式の真偽にかかわらず，b(n)+c(n) を演算③

マスクベクトルが真である要素のみ，b(n)+c(n)の演算結果を a(n)に代入

このように，IF 文を含まない DO ループに比べ，条件分岐のための処理が加わるので，分岐の数が多いと性能劣化の要因となる．

— 39 —

スーパーコンピュータ SX-9 の高速化

real,dimension(nsp)::a,b,c
integer::flag(nsp)do n=1,nsp
if(flag(n).eq.1)thena(n)=b(n)+c(n)
endif
enddo

図 4.3-15 IF 文を含んだ DO ループの例
なお，IF 文を含んでいても，図 4.3-16 のようにループの繰り返しの間，条件式の結果が不変の場合，
コンパイラは IF 文を DO ループの外に出してベクトル化するので，IF 文のオーバヘッドは発生しない．
do i=1,100a(i)=b(i)*c(i)if(isw.eq.1)thenc(i)=d(i)+e(i)
endif
enddo

コンパイラが自動的に最適化する：
do i=1,100a(i)=b(i)*c(i)enddo
if(isw.eq.1)then
do i=1,100c(i)=d(i)+e(i)
enddo
endif

図 4.3-16 IF 文が DO ループの外に展開される変形イメージ(2) マスクベクトルを用いてベクトル化しているループの性能
図 4.3-17 の DO ループは，マスクベクトルを用いてベクトル化するので TRUE，FALSE のどちらの場合も演算が実行される．演算結果は，マスクベクトルを用いて TRUE に対応する要素のみ代入される．
このため，flag_v の真の割合（以下，真率）に関係なく，DO ループの処理時間はほぼ一定になる．
do n=1,nsp
if(flag_v(n) .eq. 1) then
tmpx3(n)=(tmpx1(n)-tmpx2(n))*2.0d0
tmpy3(n)=(tmpy1(n)-tmpy2(n))*2.0d0
tmpz3(n)=(tmpz1(n)-tmpz2(n))*2.0d0
endif
enddo

図 4.3-17 一般的な IF 文を含む DO ループ
図 4.3-18 は，以下の 3 パターンで図 4.3-17 の DO ループを SX-9 1CPU で実行するときの性能値である．それぞれの実行時間がほとんど変わらないことがわかる．


Rfv_vt：flag_v の値が全て TRUE の場合

Rfv_vf：flag_v の値が全て FALSE の場合

Rfv_50：flag_v の値がランダムに TRUE/FALSE となる場合

— 40 —

PROC.NAMESENAC Vol. 45, No. 2（2012.4）

FREQUENCY

Rfv_vtRfv_vfRfv_50EXCLUSIVE
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

0.757(
0.755(
0.758(

0.3)
0.3)
0.3)

25.23725.17425.26248061.848182.2
48015.4

15954.8 99.59 256.0
15994.8 99.59 256.0
15939.4 99.59 256.0

0.757
0.755
0.758

30
30
30

0.000
0.000
0.000

0.000
0.000
0.000

0.006
0.000
0.014

0.327
0.327
0.329

図 4.3-18 FTRACE 情報の一覧(3) 事例１：IF 文をループの外に移動するチューニング
図 4.3-19 の多重ループは，TRUE 側の演算量が FALSE 側の演算量よりも大きい IF 文の例である．
マスクベクトルを用いるベクトル化は，条件式の結果に関わらず，どちらの演算も実行しているため，本来必要な演算数より多くなる．以下に，ループの入れ替えによって IF 文をループの外に移動するチューニングの例を示す．
!オリジナル
153: +------>
154: |+----->
155: ||V---->
156: |||
157: |||
158: |||
159: |||
160: |||
161: |||
162: |||
163: |||
164: |||
165: ||V---166: |+----167: +------

do k=1,nz
do j=1,ny
do i=1,nx
if(lc(i,k).eq.0) thenp(1,i,j,k) = ff*va(1,i,j,k)*(f1*vb(i,j,k)+gv10) + d1*c(i,j,k)p(2,i,j,k) = ff*va(2,i,j,k)*(f2*vb(i,j,k)+gv10) + d2*c(i,j,k)p(3,i,j,k) = ff*va(3,i,j,k)*(f3*vb(i,j,k)+gv10) + d3*c(i,j,k)else
p(1,i,j,k) = d1*c(i,j,k)p(2,i,j,k) = d2*c(i,j,k)p(3,i,j,k) = d3*c(i,j,k)
endif
enddo
enddo
enddo

図 4.3-19

TRUE 側の演算量が FALSE 側の演算量よりも大きい IF 文のコード

この IF 文の条件式「lc(i,k).eq.0」はループ変数 j に依存しておらず，ループ内の配列に依存関係がない．このため DO ループ j を IF 文の内側に移動できる．
チューニング後のコードを図 4.3-20 に示す．IF 文は DO ループ j よりも外側にあり，マスクベクトル
処理とならない．このチューニングで IF 文の条件式により内側の DO ループが選択され，真の時だけ実行される．
!チューニング
153: +------>
154: |+----->
155: ||
156: ||V---->
157: |||
158: |||
159: |||
160: ||V---161: ||162: ||V---->
163: |||
164: |||
165: |||
166: ||V---167: ||168: |+----169: +------

do k=1,nz
do i=1,nx
if(lc(i,k).eq.0) thendo j=1,nyp(1,i,j,k) = ff*va(1,i,j,k)*(f1*vb(i,j,k)+gv10) + d1*c(i,j,k)p(2,i,j,k) = ff*va(2,i,j,k)*(f2*vb(i,j,k)+gv10) + d2*c(i,j,k)p(3,i,j,k) = ff*va(3,i,j,k)*(f3*vb(i,j,k)+gv10) + d3*c(i,j,k)
enddo
else
do j=1,nyp(1,i,j,k) = d1*c(i,j,k)p(2,i,j,k) = d2*c(i,j,k)p(3,i,j,k) = d3*c(i,j,k)
enddo
endif
enddo
enddo

図 4.3-20 IF 文をループの外に移動したコード

— 41 —

スーパーコンピュータ SX-9 の高速化

チューニング前後の各ケースについて SX-9 1CPU を用いて性能評価を行う．問題サイズは「nx=256,
ny=256, nz=128」としている．IF 文の条件式の真率による評価を行うため，チューニング前後で真率をそれぞれ以下に設定する．


Res_vt： IF 文の条件式が全て TRUE の場合

Res_vf： IF 文の条件式が全て FALSE の場合



Res_25： IF 文の条件式が真率 25%となる場合



Res_50： IF 文の条件式が真率 50%となる場合



Res_75： IF 文の条件式が真率 75%となる場合!オリジナル
PROC.NAME FREQUENCY

Res_vtRes_vfRes_25
Res_50
Res_75

30
30
30
30
30

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

0.145( 11.8)
0.145( 11.8)
0.163( 13.2)
0.172( 13.9)
0.163( 13.2)

4.836
4.827
5.421
5.719
5.427

67836.667955.3
60511.5
57358.6
60447.5

36428.7 99.74 256.0
36492.4 99.74 256.0
32495.0 99.74 256.0
30801.9 99.74 256.0
32460.7 99.74 256.0

0.145
0.145
0.163
0.172
0.163

0.000
0.000
0.000
0.000
0.000

0.000
0.000
0.000
0.000
0.000

0.053
0.052
0.067
0.072
0.068

0.035
0.037
0.054
0.062
0.049

図 4.3-21 チューニング前の FTRACE 情報!チューニング
PROC.NAME FREQUENCY

Res_vtRes_vfRes_25
Res_50
Res_75

30
30
30
30
30

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

0.123( 13.9)
0.095( 10.6)
0.105( 11.8)
0.115( 13.0)
0.120( 13.4)

4.116
3.151
3.514
3.844
3.983

61321.9
18804.930567.040563.351256.7

36687.8 99.71 256.0
7987.7 99.11 256.0
16091.8 99.46 256.0
22932.7 99.60 256.0
30021.8 99.67 256.0

0.123
0.094
0.105
0.115
0.119

0.000
0.000
0.000
0.000
0.000

0.000
0.000
0.000
0.000
0.000

0.064
0.045
0.050
0.054
0.059

0.044
0.060
0.057
0.055
0.049

図 4.3-22 チューニング後の FTRACE 情報
チューニング前は，実行時間にわずかな差が出ているが，全ての場合で全浮動小数点演算数（MFLOPS×時間）の変化がない．このことから，IF 文の条件式の真率に関係なく TRUE/FALSE 両方の演算がなされていることがわかる．
チューニング後は，全浮動小数点演算数は IF 文の条件式の真率よって増減している．演算密度の違いで MFLOPS 値が低下している場合もあるが，どの場合でも処理時間の短縮が確認できる．
(4) 事例２：IF 文を他の処理に置き換えるチューニング
図 4.3-23 の例は，粒子コードに見られる周期的境界条件で位置情報の補正を行うコードである．
DO ループ中に多数の IF 文があるため，多くのマスク付きベクトル演算を行い，演算時間が長くなっている．この根本的な解決として IF 文に代わる処理への置き換えを実施する．
!チューニング前
54: V------>
55: |
56: |
57: |
58: |
59: |
60: |
61: V------

do n=1,nsp
if(pa(n,1) < 0.0d0) pa(n,1)=pa(n,1)+xdim
if(pa(n,1) >= xdim) pa(n,1)=pa(n,1)-xdim
if(pa(n,2) < 0.0d0) pa(n,2)=pa(n,2)+ydim
if(pa(n,2) >= ydim) pa(n,2)=pa(n,2)-ydim
if(pa(n,3) < 0.0d0) pa(n,3)=pa(n,3)+zdim
if(pa(n,3) >= zdim) pa(n,3)=pa(n,3)-zdim
enddo

図 4.3-23 IF 文による周期的境界条件のコード
図 4.3-23 のコードは，位置情報が各方向の最小値から最大値の範囲に収まっているかを，IF 文の— 42 —

SENAC Vol. 45, No. 2（2012.4）

条件式で判定している．収まっていない場合は，位置情報に各辺の長さを加減して周期的境界条件を満たす補正がなされる．
図 4.3-24 にチューニング後のコードを示す．ここでは，mod 関数を使用して周期的境界条件を満たす同等な補正を行っている．全要素一律に各辺の長さを加算し，その長さで割った余りを求めることで周期的境界条件が満たされ，IF 文による場合分けが不要となる．ただし，周期内の条件にある位置は，
mod 関数による補正計算が行われるため，チューニング前のコードと比べて若干の誤差が生じる．どちらのコードもループ長は同じである．
!チューニング後
54: V------>
55: |
56: |
57: |
58: V------

do n=1,nsp
pa(n,1)=dmod(pa(n,1)+xdim,xdim)
pa(n,2)=dmod(pa(n,2)+ydim,ydim)pa(n,3)=dmod(pa(n,3)+zdim,zdim)
enddo

図 4.3-24 mod 関数による周期的境界条件のコード
図 4.3-23 と図 4.3-24 の各場合について SX-9 1CPU を用いて性能評価を行う．問題サイズは「nx=256, ny=128, nz=64, ns=32」としている．要素数 nsp は nx, ny, nz, ns の積となり，各方向の最大値を実数型変数 xdim, ydim, zdim で与えている．
PROC.NAME

sample_f

FREQUENCY30

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

9.557( 85.8)

318.582

6327.7

1895.8 99.87 256.0

9.557

0.000

0.000

0.021

8.658

図 4.3-25 チューニング前の FTRACE 情報
PROC.NAME

sample_f

FREQUENCY30

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

1.720( 51.5)

57.34552707.6

17554.0 99.91 256.0

1.720

0.000

0.000

0.004

0.029

図 4.3-26 チューニング後の FTRACE 情報
図 4.3-25 と図 4.3-26 がそれぞれチューニング前後の性能値で，ベクトル演算率やベクトル長に大きな差異は見られない．チューニング前は配列 pa のロード・ストア回数が多く，条件分岐のオーバヘッ
ドで待ち時間も発生し，「BANK CONFLICT/NETWORK」が大きくなっている．チューニング後は mod関数を用い IF 文の条件分岐が取り除かれたことで，ベクトル演算の密度が上がり MOPS およびMFLOPS の値が向上している．

4.3.4 DO ループ内の作業配列の削除による高速化(1) 作業配列の削除の概要
一般にプログラムは入力データと出力結果の他にも，途中経過を作業配列に保存する場合がある．
しかし，不必要な作業配列の利用はメモリ負荷を増加させ，演算性能を低下させることになる．
図 4.3-27 の例は，ループ内でのみ使用する途中の演算結果 a(n)+b(n)を作業配列 work(n)に保存している．この途中の演算結果をその後の処理で利用することがなければ，作業配列 work(n)を使わず
スカラ変数に書き換え，作業配列への保存を回避できる．このようにすることで，途中の演算結果は作業レジスタに保存され，メモリ負荷が下がり演算性能が向上する．

— 43 —

スーパーコンピュータ SX-9 の高速化

subroutine sample_f(a,b,p,nsize)
real*8 :: work(nsize)：
do n=1,nsize
work(n)=a(n)+b(n)p(n)=work(n)*0.5
enddo：
end subroutine sample_f

ループ内でのみ使用する値は，配列に保存せずスカラ変数を使う．

図 4.3-27 作業配列からスカラ変数に書き換えることができる例(2) 作業配列を削減するチューニングの事例
図 4.3-28 のサブルーチン sample_f1, sample_f2 は，プログラム中の複数の箇所で同じ演算となる部分を括りだし，共通部品としてサブルーチンにしたものである．このプログラムは，3 次元配列のデータ
を 1 次元作業配列にコピーしサブルーチン sample_f1 へ渡し，サブルーチン sample_f1, sample_f2 で計算し，1 次元作業配列の結果を 3 次元配列に書き戻している．
!オリジナル
34: |+----->
35: ||+---->
36: |||
37: |||V--->
38: ||||
39: ||||
40: ||||
41: ||||
42: |||V--43: |||44: |||
45: |||
46: |||
47: |||V--->
48: ||||
49: |||V--50: |||51: ||+---52: |+-----

do k=1,kmax
do j=1,jmax
do i=1,imax
wi(1,i)=a(1,i,j,k)wi(2,i)=a(2,i,j,k)wi(3,i)=a(3,i,j,k)
wv(i) =b(i,j,k)
enddo
call sample_f1(imax,wi,c ,wt)
call sample_f2(imax,wv,wt,wo)
do i=1,imax
p(i,j,k)=wo(i)
enddo
enddo
enddo

85:
：
97: V------>
98: |
99: |
100: |
101: V------

subroutine sample_f1(lmax,wi,wc,wt)

106:：
118: V------>
119: |
120: V------

subroutine sample_f2(lmax,wv,wt,wo)

do l=2,lmax-1
wt(l)=wc(1)*(wi(1,l+1)-wi(1,l-1)) &+wc(2)*(wi(2,l+1)-wi(2,l-1)) &+wc(3)*(wi(3,l+1)-wi(3,l-1))
enddo

do l=2,lmax-1
wo(l)=wt(l)*wv(l)enddo

図 4.3-28 作業配列を削減前のコード
図 4.3-29 のチューニングは，作業配列 wi,wv と元の入力配列 a,b の 1 列が同じ値を参照していることに着目する．この場合，作業配列 wi,wv を使わないで配列 a,b の 1 列を直接引き渡すことができる．
同様に引数の作業配列 wo と出力配列 p の 1 列が同じ値を更新しているので，出力データも作業配列
wo を使わないで配列 p へ直接出力できる．このようにサブルーチン呼び出し前後に作業配列 wi, wv,
wo を使用しないでデータを引き渡すことができるため，作業配列へのメモリアクセス回数が削減できる．

— 44 —

SENAC Vol. 45, No. 2（2012.4）!チューニング①
58: |+----->
59: ||+---->
60: |||
61: |||
62: |||
63: |||
64: ||+---65: |+-----

do k=1,kmax
do j=1,jmax
call sample_f1(imax,a(1,1,j,k),c ,wt)
call sample_f2(imax,b(1,j,k) ,wt,p(1,j,k))
enddo
enddo

85:
：
97: V------>
98: |
99: |
100: |
101: V------

subroutine sample_f1(lmax,wi,wc,wt)

106:：
118: V------>
119: |
120: V------

subroutine sample_f2(lmax,wv,wt,wo)

do l=2,lmax-1
wt(l)=wc(1)*(wi(1,l+1)-wi(1,l-1)) &+wc(2)*(wi(2,l+1)-wi(2,l-1)) &+wc(3)*(wi(3,l+1)-wi(3,l-1))
enddo

do l=2,lmax-1
wo(l)=wt(l)*wv(l)enddo

図 4.3-29 ［第 1 段階］入力データと出力データを引数にする修正
図 4.3-30 の更なるチューニングは，サブルーチン sample_f1, sample_f2 を 1 つのサブルーチン
sample_f にまとめ，これらのサブルーチンそれぞれにあった DO ループを融合している．このループ融合によって，sample_f1 から sample_f2 に受け渡している途中の結果を，作業配列 wt を使わずに演算処理できる．
!チューニング②
71: |+----->
72: ||+---->
73: |||
74: |||
75: |||
76: ||+---77: |+----127:：
141: V------>
142: |
143: |
144: |
145: |
146: V------

do k=1,kmax
do j=1,jmax
call sample_f(imax,a(1,1,j,k),b(1,j,k),c,p(1,j,k))
enddo
enddo
subroutine sample_f(lmax,wi,wv,wc,wo)
do l=2,lmax-1
s=wc(1)*(wi(1,l+1)-wi(1,l-1)) &+wc(2)*(wi(2,l+1)-wi(2,l-1)) &+wc(3)*(wi(3,l+1)-wi(3,l-1))
wo(l)=s*wv(l)enddo

図 4.3-30 ［第 2 段階］サブルーチンにまたがる 2 つの DO ループを融合する修正(3) 性能評価
チューニング前後の各場合について SX-9 1CPU を用いて性能評価を行う．(2)で示したコードを用
い，問題サイズは「imax=1000, jmax=50, kmax=50」としている．図 4.3-31 は，チューニング前後の性能値である．
case_000：図 4.3-28 チューニング前の結果



case_001：図 4.3-29 チューニング［第 1 段階］での結果

— 45 —

スーパーコンピュータ SX-9 の高速化

case_002：図 4.3-30 チューニング［第 2 段階］での結果

PROC.NAMEFREQUENCY

case_000
case_001
case_002

1000
1000
1000

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

5.841(100.0)
1.821(100.0)
1.064(100.0)

5.841
1.821
1.064

12542.9
26312.640150.03844.5 98.83 249.7
12331.1 98.93 249.5
21097.7 99.26 249.5

5.840
1.820
1.064

0.000
0.000
0.000

0.000
0.000
0.000

1.862
0.520
0.516

4.740
1.145
0.423

図 4.3-31 FTRACE 情報一覧
case_000→case_001→case_002 とチューニングの段階を経て，「BANK CONFLICT/ NETWORK」の時間が短縮されている．これは作業配列へのデータ転送が減っているためである．データ転送が削減されメモリ負荷が下がったことで，結果的に演算器の待ち時間が減少し，MFLOPS 値が向上している．

4.3.5 組み込み関数 MAX/MIN の利用による高速化(1) 組み込み関数 MAX/MIN の利用による高速化の概要
最大値・最小値を求める計算はマクロ演算を利用してベクトル化され，図 4.3-32 に示すようにコンパ
イル診断メッセージに「Macro operation Max/Min」と表示される．また，4.1 節で述べたとおり，SX-9 から MAX/MIN 関数がハードウェア命令化され，最大値・最小値を求める計算時間の短縮が行われる．
本節では，最大値を求めるコードでベクトル化される場合とベクトル化されない場合を示し，さらにハードウェア命令が使われた場合と，使われない場合の性能の違いを示す．
61

vec

(

26): Macro operation Max/Min.

：（※中略）
57: W------>
58: |*----->
59: ||*---->
60: |||
61: |||
62: ||*---63: |*----64: W------

行番号 61 番で最大値・最小値のマクロ演算を適用するメッセージ

do k=1,kmax
do j=1,jmax
do i=1,imax
tmp=r1(i,j,k)*r2(i,j,k)s=max(s,tmp)
enddo
enddo
enddo

図 4.3-32 MAX/MIN 関数利用の例とコンパイル診断メッセージ(2) 最大値を求めるコードの事例
次に示す図は，3 種類の最大値を求めるコードで，以下の特徴がある．


コード①（図 4.3-33）は，組み込み関数 MAX を利用して最大値を求める

コード②（図 4.3-34）は，大小関係を IF 文で判断し，大きい値で変数 s を更新し最大値を求める

コード③（図 4.3-35）は，ループの 1 回目で初期値を設定し，2 回目以降のループで最大値を求める

コード①は，組み込み関数 MAX を利用しているので SX-9 のハードウェア命令が使われベクトル化
される．コード②はコンパイラが最大値を求める処理と認識できるのでベクトル化される．コード③は条件分岐が複雑であり，コンパイラが最大値を求める処理と認識できないのでベクトル化されない．

— 46 —

！コード①
131:
：
143:
144: V------>
145: |
146: |
147: V------

SENAC Vol. 45, No. 2（2012.4）

subroutine sample_b1(r1,r2,s)s=-999d0
do i=1,imax*jmax*kmax
tmp=r1(i)*r2(i)s=max(s,tmp)
enddo

図 4.3-33 最大値を求めるコード①
！コード②
154:
：
166:
167: V------>
168: |
169: |
170: V------

subroutine sample_b2(r1,r2,s)s=-999d0
do i=1,imax*jmax*kmax
tmp=r1(i)*r2(i)if(tmp .gt. s) s=tmp
enddo

図 4.3-34 最大値を求めるコード②
！コード③
177:
：
190:
191: +------>
192: |
193: |
194: |
195: |
196: |
197: |
198: |
199: +------

subroutine sample_b3(r1,r2,s)lflag=.true.
do i=1,imax*jmax*kmax
tmp=r1(i)*r2(i)if(lflag) thens=tmp
lflag=.false.
else
s=max(s,tmp)
endif
enddo

図 4.3-35 最大値を求めるコード③(3) 性能評価
各コードを SX-9 1CPU を用いて性能評価を行う．問題サイズは「imax*jmax*kmax=2563」で十分に長いループ長としている．

コード①
コード②
コード③

表 4.3-2 3 種類コードの性能
性能値
20,999MFLOPS
11,373MFLOPS
93MFLOPS

コード①は MAX/MIN 関数でベクトル化され SX-9 で強化されたハードウェア命令が使われるため，高い性能を示している．コード②はベクトル化されるが SX-9 のハードウェア命令が使われてないため，
コード①と比較して半分程度の性能である．コード③はベクトル化されていないので性能が出ていない．
SX-9 で最大値・最小値を求める場合，高い性能を得るには MAX/MIN 関数を使用してベクトル化することが必要である．

— 47 —

スーパーコンピュータ SX-9 の高速化

4.3.6 科学技術計算ライブラリＡＳＬの利用による高速化(1) 科学技術計算ライブラリ利用の概要
数値計算などで利用している連立 1 次方程式，固有値・固有ベクトル，高速フーリエ変換(FFT)，特殊関数，乱数などの機能は計算ライブラリとして準備されている．このような計算ライブラリをユーザが作成することはできるが，ハードウェアの性能を引き出すために複雑な最適化が必要であり，高い性能を得るのは容易ではない．また一般の計算ライブラリは，特定のハードウェアに向けての最適化は通常行われておらず，SX-9 の最大性能を得ることは難しい．SX-9 には科学技術計算ライブラリ ASL があり，
これは広範な分野の数値シミュレーションプログラムの作成を強力に支援する計算ライブラリである．
ASL はハードウェアの性能を最大限に引き出す最適化が行われており，ユーザはこれらのライブラリを使うことにより，SX-9 の性能を引き出すことが容易に行える．
これから科学技術計算ライブラリ ASL の FFT とフリーソフトウェアの FFT ライブラリとの性能を比較す
る．比較対象のフリーソフトウェアは，NETLIB の FFTPACK，および筑波大学高橋大介氏作成のFFTE 4.1 としている．それぞれの FFT ライブラリは，以下の URL からダウンロードできる．
FFTPACK

http://www.netlib.org/fftpack/

FFTE

http://www.ffte.jp/(2) ASL と FFTPACK の性能評価
ここでは１次元複素数フーリエ変換を ASL のサブルーチン JFCOSI/JFCOSF と FFTPACK のサブ
ルーチン ZFFTI/ZFFTF/ZFFTB を用いて比較する．それぞれの評価コードを図 4.3-36 と図 4.3-37に示す．これらサブルーチンは周期２πの１周期を n 等分した n 個の複素数データが与えられたとき，
n 個のデータに対する離散フーリエ変換（順変換と逆変換）を行う．
! 1 次元複素数 FFT（ASL）の評価コード
integer, parameter :: ld=2**26
complex(kind=8) :: c(ld),trigs((3*ld)/2),wc(ld)integer :: ifax(128)real(kind=8) :: t1,t2,tm,perf(-1:1)
do ii=11,26
n=2**ii
do isw=1,-1,-2
do i=1,nc(i) = dcmplx(cos(sqrt(2.d0)*i),sin(sqrt(2.d0)*i))
enddo
call �F��S�(n,ifax,trigs,ierr)call clock(t1)call �F��SF(n,c,ld,isw,ifax,trigs,wc,ierr)call clock(t2)tm=(t2-t1)
perf(isw)=5.d0*n*log(dble(n))/log(2.d0)/tm*1.d-9
enddo
write(6,*) n,perf(1),perf(-1)
enddo
end

図 4.3-36 1 次元複素数 FFT（ASL）の評価コード

ASL の FFT サブルーチン

— 48 —

SENAC Vol. 45, No. 2（2012.4）! 1 次元複素数 FFT（FFTPACK）の評価コード
integer, parameter :: ld=2**26
complex(kind=8) :: w(ld*2+15),cx(ld)real(kind=8) :: t1,t2,tm,perf(-1:1)
do ii=11,26
n=2**ii
do isw=1,-1,-2
do i=1,ncx(i) = dcmplx(cos(sqrt(2.d0)*i),sin(sqrt(2.d0)*i))
enddo
call �FFT� (n,w)
call clock(t1)if(isw.eq.1)thencall �FFTF (n,cx,w)
else
call �FFT� (n,cx,w)
endif
call clock(t2)tm=(t2-t1)
perf(isw)=5.d0*n*log(dble(n))/log(2.d0)/tm*1.d-9

FFTPACK のサブルーチン

enddo
write(6,*) n,perf(1),perf(-1)
enddo
end

図 4.3-37 1 次元複素数 FFT（FFTPACK）の評価コード
ここで性能測定プログラムに使う評価コードについて説明する．この評価コードでは，必要なデータ領域を確保し，FFT の入力データおよびワーク領域の初期化を行って１次元複素数フーリエ変換のサブルーチンを呼び出している．実行時間の計測は，フーリエ変換のサブルーチンを呼び出す前後で，タイマールーチン（CLOCK）を挟んで行っている．
つぎに性能の算出方法について説明する．FFT の代表的な計算方法である Cooley-Tukey のアル
ゴリズムからデータ数 n の１次元複素数 FFT の演算数を数え上げると，次式で表わされる．
演算数＝5n×log2(n)この Cooley-Tukey のアルゴリズムの演算数は基数が 2 の場合の結果であるが，基数が 2 以外でも大きく違わないため，演算数を求めるには良い近似値となる．この方法で FFT の演算数を求め，実行時間で除する事で性能(GFLOPS 値)を算出している．
図 4.3-38 は，データ数を 211 (=2,048) から 226(=67,108,864)まで 2 のべき乗で増加させた時の FFTの性能である．図中，fwd は順変換を bwd は逆変換を意味している．

— 49 —

スーパーコンピュータ SX-9 の高速化

１次元FFT
45
40
35
GFLOPS

30

FFTPACK(fwd)
FFTPACK(bwd)ASL(fwd)
ASL(bwd)

25
20
15
10
5
0

211 212 213 214 215 216 217 218 219 220 221 2 22 2 23 2 24 2 25 2 26
データ数

図 4.3-38 1 次元複素数 FFT の性能比較
FFTPACK は，データ数を増加しても大きな性能向上が得られない． ASL は，データ数を増加すると 217（=131,072）で 36GFLOPS と性能が向上し，最大性能は 224(=16,777,216）で 41GFLOPS となる．
ASL と FFTPACK の性能を比較すると，データ数が大きいところで約 5 倍 ASL が高速である．
(3) ASL と FFTE の性能評価
ここでは２次元複素数フーリエ変換を ASL のサブルーチン ZFC2FB/ZFC2BF と，FFTE はベクトル用のサブルーチン ZFFT2D を用いて行っている（注 1）．それぞれの評価コードを，図 4.3-39 と図
4.3-40 に示す．これらのサブルーチンは，２次元の各次元のデータで離散フーリエ変換（順変換と逆変換）を行う．

— 50 —

SENAC Vol. 45, No. 2（2012.4）! ２次元複素数 FFT（ASL）の評価コード
integer, parameter :: ld=10241
complex(kind=8) :: c(ld*ld),trigs(2*ld),wc(ld*ld)
integer :: ifax(40)real(kind=8):: t1,t2,tm,perf(-1:1)
do ii=1024,10240,1024
do isw=1,-1,-2
nx=ii
ny=ii
lx=nx+1
ly=ny+1
call init(c,lx,ly,nx,ny)
call �F�2F�(nx,ny,c,lx,ly,0,ifax,trigs,wc,ierr)call clock(t1)call �F�2�F(nx,ny,c,lx,ly,isw,ifax,trigs,wc,ierr)call clock(t2)tm=(t2-t1)
perf(isw)=5.d0*nx*ny*(log(dble(nx))+log(dble(ny)))/log(2.d0)/tm*1.d-9

ASL の FFT サブルーチン

enddo
write(6,*) ii,perf(1),perf(-1)
enddo
end
subroutine init(c,lx,ly,nx,ny)
complex(kind=8) :: c(lx,ly)
do j=1,ny
do i=1,nxc(i,j) = dcmplx(cos(sqrt(2.d0)*i),sin(sqrt(2.d0)*j))
enddo
enddo
end

図 4.3-39 2 次元複素数 FFT（ASL）の評価コード! ２次元複素数 FFT（FFTE）の評価コード
integer, parameter :: ld=10240
complex(kind=8) :: a(ld*ld)
real(kind=8) :: t1,t2,tm,perf(-1:1)
do ii=1024,10240,1024
if(mod(ii,7).eq.0) cycle
do isw=1,-1,-2
nx=ii
ny=ii
call init(a,nx,nx,nx,ny)
call �FFT2�(a,nx,ny,0)
FFTE のサブルーチン
call clock(t1)call �FFT2�(a,nx,ny,-isw)
call clock(t2)tm=(t2-t1)
perf(isw)=5.d0*nx*ny*(log(dble(nx))+log(dble(ny)))/log(2.d0)/tm*1.d-9
enddo
write(6,*) ii,perf(1),perf(-1)
enddo
end
subroutine init(c,lx,ly,nx,ny)
complex(kind=8) :: c(lx,ly)
do j=1,ny
do i=1,nxc(i,j) = dcmplx(cos(sqrt(2.d0)*i),sin(sqrt(2.d0)*j))
enddo
enddo
end

図 4.3-40 2 次元複素数 FFT（FFTE）の評価コード

— 51 —

スーパーコンピュータ SX-9 の高速化

図 4.3-41 は 2 次元データの各次元とも 1,024 の増分で 1,024×1,024 から 10,240×10,240 までデー
タ数を増加させた時の FFT の性能である．ただし， FFTE はデータ数に 2, 3, 5 以外の約数が入ると計算できないため，約数に 7 が入るデータ数は除外している（注 2）．実行時間の計測は，(2)と同様，
フーリエ変換のサブルーチンを呼び出す前後で行い，タイマールーチンを呼ぶ前に FFT の入力デー
タとワーク領域の初期化を行っている．演算数は１次元の FFT での算出方法を２次元の FFT に拡張し，性能を求めている．

２次元FFT60

GFLOPS

50
40
30
20
10
0
1,024 2,048 3,072

4,096 5,120

6,144

7,168

8,192

9,216 10,240

データ数
FFTE(fwd)

FFTE(bwd)

ASL(fwd)

ASL(bwd)

図 4.3-41 2 次元複素数 FFT の性能比較
FFTE はデータ数により 32～42GFLOPS の性能である．ASL はデータ数が 8,192×8,192 のとき，最大性能の 54GFLOPS となる．ASL と FFTE の性能を比較すると 1.1～1.4 倍 ASL が高速である．
（注 1）FFTE は，ベクトル用およびスカラ用に最適化された FFT のサブルーチンが用意されている．
それぞれのプラットフォームに適しているサブルーチンを利用すること．
（注 2）FFTE は，データ数に 2,3,5 以外の約数が入っていると FFT の計算ができない仕様である．
例えば 7,168 は約数に 7 が入り計算できない．一方 ASL の FFT は，データ数に約数の制
約はなく任意のデータ数で FFT の計算ができる．
(4) 性能評価のまとめ
科学技術計算ライブラリ ASL の FFT サブルーチンとフリーソフトウェアの FFT ライブラリである
FFTPACK，FFTE との間で性能を比較している．その結果，ASL の FFT サブルーチンが高い性能を示している．FFTE にはベクトル用に最適化されている FFT サブルーチンが用意されており，実行効率
が 30%～40%と比較的よい性能が得られているが，ASL の FFT サブルーチンは実行効率が最大で 50%を超え，より高速である．
SX で FFT の計算を行う場合は，SX のハードウェア向けに高度な最適化が施されている科学技術計

— 52 —

SENAC Vol. 45, No. 2（2012.4）

算ライブラリ ASL を使うことで，フリーソフトウェアの FFT ライブラリを上回る性能が引き出せることがわかる．
4.3.7 リストベクトルを含む DO ループの最適化(1) リストベクトルの概要
リストベクトルとは，図 4.3-42 に示すように間接参照される配列 IX にしたがって，配列 A の参照また
は更新をベクトル処理で行うものである．図 4.3-43 に図 4.3-42 に示したコードのメモリアクセスパター
ンを示す．リストベクトルを用いることで，連続でないランダムなメモリアクセスをベクトル処理として行うことができる．リストベクトルは DO ループ中の IF 文の除去や，ベクトル化できないデータ参照関係を回避する平面法などに活用されている．平面法については，「高速化活動報告第 4 号」を参照のこと．
!参照
do i=1,nT(i)=A(IX(i))
enddo
!更新
do i=1,nA(IX(i))=T(i)enddo

図 4.3-42 リストベクトルのコード参照する要素のアドレスを計算変更する要素のアドレスを計算

IX

A

レジスタ

T

IX

T

レジスタ

A

3

a1

a3

a3

1

t1

t1

t1

8

a2

a8

a8

8

t2

t2

t6

2

a3

a2

a2

5

t3

t3

t4

4

a4

a4

3

t4

t4

a4

10

a5

a4
a10

a10

10

t5

t5

t3

5

a6

a5

a5

2

t6

t6

a6

a7

a7

a8

t2

a9

a9

a10

t5
T(i)=A(IX(i))

A(IX(i))=T(i)図 4.3-43 リストベクトルのメモリアクセス(2) 事例１：リスト参照のメモリレイテンシの隠ぺい
リストベクトルによる参照や更新が同一 DO ループ内に複数回現れるプログラムの高速化の手法について説明する．
図 4.3-44 にリストベクトルによる参照や更新が同一 DO ループ内に複数回現れるコード（注）を示す．

スーパーコンピュータ SX-9 の高速化

— 53 —

間接参照 n2 による配列 p の参照は，間接参照 n1 による同配列 p への更新の完了を待つためメモリレ
イテンシ（遅延）が表面化し，演算性能を発揮できない．そこで配列 p のリストベクトルによる参照と更新のタイミングをずらすため，一度別の作業配列に連続データとして保存する．
121: +------>
122: |
123: |V----->
124: ||
125: ||：
141: ||
142: ||
143: ||
144: ||：
161: ||
162: ||：
191: ||
192: ||：
209: |V----210: +------

do icol=1,ncol!cdir nodep
do edge=esep(icol),esep(icol+1)-1
n1=e2n(edge,1)
n2=e2n(edge,2)x1=u(n1,1)y1=u(n1,2)x2=u(n2,1)y2=u(n2,2)!userfunc はインライン化によりベクトル化されるユーザ定義関数
v=userfunc(x1,y1,z1,x2)
p(n1)=p(n1)+v
p(n2)=p(n2)-v
enddo
enddo

図 4.3-44 リストベクトルが同一 DO ループに複数回現れるコード
図 4.3-45 に作業配列を使って最適化するコードを示す．先ず配列 p を間接参照 n1,n2 により参照
し，それぞれの作業配列 pn1,pn2 に連続データとして保存する．次に作業配列 pn1,pn2 を使って連続アクセスのデータで演算する．最後に演算結果を，作業配列から元の配列 p へ書き戻している．このよ
うにすることで，配列 p の間接参照による参照を更新の前にまとめて行い，更新完了の待ちによる間接
参照時のメモリレイテンシの影響を少なくできる．

— 54 —

SENAC Vol. 45, No. 2（2012.4）

121: +------>
122: |
123: |V----->
124: ||
125: ||：
141: ||
142: ||
143: ||
144: ||：
161: ||
162: ||：
179: |V----180: |181: |V----->
182: ||：
211: ||
212: ||：
229: |V----230: |
231: |
232: |V----->
233: ||
234: ||：
250: ||
251: ||：
268: |V----269: +------

do icol=1,ncol! リストベクトルによる参照処理の集約
do edge=egsep(icol),egsep(icol+1)-1
n1=e2n(edge,1)
n2=e2n(edge,2)x1_a(edge)=u(n1,1)y1_a(edge)=u(n1,2)x2_a(edge)=u(n2,1)y2_a(edge)=u(n2,2)pn1(edge)=p(n1)pn2(edge)=p(n2)enddo! 連続ベクトルによる演算
do edge=egsep(icol),egsep(icol+1)-1
v=userfunc(x1_a(edge),y1_a(edge),x2_a(edge),y2_a(edge))pn1(edge)=pn1(edge)+v
pn2(edge)=pn2(edge)-v
enddo! リストベクトルによる更新処理の集約!cdir nodep
do edge=egsep(icol),egsep(icol+1)-1
n1=e2n(edge,1)
n2=e2n(edge,2)p(n1)=pn1(edge)p(n2)=pn2(edge)
enddo
enddo

図 4.3-45 作業配列を使った演算コード
チューニング前後のコードを使って，SX-9 1CPU で性能評価を行う．図 4.3-46 に FTRACE の情報を示す．org がチューニング前のコード，tune がチューニング後のコードを表している．配列 p の間接参照を更新より先に一括して行っていることで，メモリレイテンシを隠ぺいでき，図 4.3-46 で示されるようにメモリネットワーク競合時間（NETWORK）が半減し，CPU 時間比にして 1.95 倍に性能が向上している．
PROC.NAME

org
tune

FREQUENCY100
100

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

24.849( 12.8)
12.707( 6.5)

248.493127.0746508.816339.5

1817.5 99.95 250.1
3554.2 99.80 250.1

24.84312.640

0.001
0.016

0.001
0.012

2.039
2.301

18.9728.792

図 4.3-46 チューニング前後の FTRACE 情報（注）本コードは，間接参照 n1, n2 がそれぞれ最内の DO ループの中で重複する値を取らないことと，間接参照 n1 と n2 の間にも値の重複がないことを想定している．いずれも重複がある場合，ベクトル化によりデータの不整合が起こる．コンパイラはデータ依存性を判断できないため，
図 4.3-44 の 122 行目および図 4.3-45 の 231 行目に指示行「!cdir nodep」を追加し，コンパイラに重複がないことを指示することでベクトル化を行っている．

— 55 —

スーパーコンピュータ SX-9 の高速化

(3) 事例２：配列の次元の入れ替えによるリストベクトルの回避
リストベクトル処理はメモリアクセスのレイテンシを隠ぺいすることが難しく，メモリを連続にアクセスする処理と比べて，処理時間を多く必要とし実効性能が低下する．
図 4.3-47 にリストベクトルのコードを示す．最内の i の DO ループの配列 ain はリスト il2, il1, ir1, ir2
により間接参照され，リストベクトルとなっている．
68: +------>
69: |+----->
70: ||V---->
71: |||
72: |||
73: |||
74: |||
75: |||
76: |||
77: |||
78: |||
79: |||
80: |||
81: |||

do k=1,kmax
do j=1,jmax
do i=2,imax
il2=ilir(i,1)
il1=ilir(i,2)
ir1=ilir(i,3)
ir2=ilir(i,4)
stbc=wstbc(i)dqr0=ain(i,j,k)
dqr1=ain(ir1,j,k)dqr2=ain(ir2,j,k)dql0=ain(i-1,j,k)dql1=ain(il1,j,k)dql2=ain(il2,j,k)

図 4.3-47 リストベクトル処理のコード
配列 ain に対して，i,j,k の次元から j,k,i の次元へ入れ替えを行うとともに，i の DO ループを最外に移動させる．図 4.3-48 に配列の次元を入れ替えたコードを示す．この変更により，配列 ain の参照は
j,k の DO ループに対し連続ベクトルとなり，リストベクトルが回避される．加えてこの変更により，添え字
j が配列の下限から上限までの値を取ることで j,k の DO ループが一重化でき，ベクトル長が長くなる（注）．
68: V----->
69: |
70: |
71: |
72: |
73: |
74: |W---->
75: ||*--->
76: |||
77: |||
78: |||
79: |||
80: |||
81: |||

do i=2,imax
il2=ilir(i,1)
il1=ilir(i,2)
ir1=ilir(i,3)
ir2=ilir(i,4)
stbc=wstbc(i)do k=1,kmax
do j=1,jmax
dqr0=ain(j,k,i)
dqr1=ain(j,k,ir1)
dqr2=ain(j,k,ir2)
dql0=ain(j,k,i-1)
dql1=ain(j,k,il1)
dql2=ain(j,k,il2)

図 4.3-48 配列の次元を入れ替えたコード
チューニング前後のコードを使って，SX-9 1CPU で性能評価を行う．図 4.3-49 に FTRACE の情報を示す． org がリストベクトル処理主体のコード，tune が配列の次元を入れ替えたコードを表している．
リストベクトルの回避とベクトル長が長くなることで，メモリネットワーク競合時間が大幅に削減され，CPU
時間比にして 11.41 倍に性能が向上している．
PROC.NAME

org
tune

FREQUENCY

1
1

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

33.964( 30.1) 33963.691
2.977( 2.6) 2976.662

7920.346796.1

1748.9 99.55 99.0
19955.2 99.50 250.0

33.964
2.976

% )

0.000
0.000

図 4.3-49 チューニング前後の FTRACE 情報

0.000
0.000

1.875
0.010

29.9671.418

— 56 —

SENAC Vol. 45, No. 2（2012.4）（注）図 4.3-48 のコードにおける j,k の DO ループは，ループ長 jmax×kmax の一重ループとして処理される．例えば ain(j,k,i)への参照は，1～jmax×kmax の値を取るループ変数 jk を用いることで ain(jk,1,i)への参照と同様に処理される．
(4) 事例３：IF 文の変更によるリストベクトルの回避
図 4.3-50 のコードは，22 行目の IF 文で決定される変数 j を用いて配列 c を間接参照しているため，リストベクトルとなる（注）．この場合，IF 文の処理内容を変更することでリストベクトルを回避できる．
21: V------>
22: |
23: |
24: |
25: |
26: |
27: |
28: |
29: V------

do i=1,nif ( x(i) .ge. 0.0 ) then
j=2
else
j=1
endif
z1(i)=z1(i)+y1(i)*c(1,j)z2(i)=z2(i)+y2(i)*c(2,j)
enddo

図 4.3-50 IF 文で決定した変数を用いて間接参照を行うコード
図 4.3-51 にリストベクトルの回避コードを示す．22 行目以降の IF 文を配列 c の値を直接参照し，新たに導入した作業変数 c1,c2 に格納する．これらの作業変数 c1,c2 を用いた 29 行目以降の演算は連続ベクトルとなるため，リストベクトルが回避される．加えて，IF 文の処理は配列 c の値を定数としてレジスタ上で利用できるため，メモリアクセスも大幅に低減される．
21: V------>
22: |
23: |
24: |
25: |
26: |
27: |
28: |
29: |
30: |
31: V------

do i=1,nif ( x(i) .ge. 0.0 ) thenc1=c(1,2)c2=c(2,2)else
c1=c(1,1)c2=c(2,1)
endif
z1(i)=z1(i)+y1(i)*c1z2(i)=z2(i)+y2(i)*c2
enddo

図 4.3-51 IF 文内で配列を直接参照するコード
チューニング前後のコードを使って，SX-9 1CPU で性能評価を行う．図 4.3-52 に FTRACE の情報を示す．org がリストベクトルのコード，tune がリストベクトル回避コードを表している．リストベクトルの回避とメモリアクセスの低減によりメモリネットワーク競合時間が大幅に削減され，CPU 時間比にして 9.90
倍に性能が向上している．
PROC.NAME

org
tune

FREQUENCY

1
1

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

22.664( 27.6) 22663.970
2.287( 2.8) 2286.695

8405.870192.4

1764.9 99.73 255.8
17492.5 99.68 255.8

22.664
2.287

% )

0.000
0.000

0.000
0.000

1.275
0.023

20.0020.693

図 4.3-52 チューニング前後の FTRACE 情報（注）ループ変数 i の値それぞれに応じた変数 j の値は，ベクトルレジスタもしくはコンパイラが自動生成する作業配列に，一度ベクトルデータとして格納される．そのため j による配列 c の間接参照は，この格納されたベクトルデータを元にリストベクトルとして処理される．

— 57 —

スーパーコンピュータ SX-9 の高速化

4.3.8 SOR 法のベクトル化（Red Black 法）(1) SOR 法の概要
SOR 法（Successive Over-Relaxation method）とは，n 元連立方程式 Ax=b を反復法で解く手法の一つである．図 4.3-53 に示すとおり，ある要素 p(i,j,k)の計算には近傍の要素 p(i-1,j,k)，p(i+1,j,k)，p(i,j-1,k)，p(i,j+1,k)，p(i,j,k-1)，p(i,j,k+1)が必要な計算を行う．

図 4.3-53 SOR 法計算の関係図
図 4.3-54 に SOR 法のコードを示す．p(i,j,k)の値を計算する際に，p(i-1,j,k)の要素を参照する必要がある．ベクトル化の対象となる最内の DO ループでは，p(i-1,j,k)の値は 1 回前で定義されており（注1），この定義・参照関係がベクトル化を阻害する依存関係になる．そのためベクトル化できず（注 2），
SX-9 で高い実効性能を得ることができない．
（注 1）p(i+1,j,k)は参照された後で値が更新されるので，ベクトル化を阻害する要因にはならない．
（注 2）図 4.3-54 の例のように，コンパイラはベクトル化可能な部分のみベクトル化を行う．これは部分ベクトル化と呼ばれ編集リストでは「S」で示される．
115
116
117
:
141
142
143
144
145
:
148
149
150
151

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

do k = n_bnd + 1, n_bnd + n_cell
do j = n_bnd + 1, n_bnd + n_cell
do i = n_bnd + 1, n_bnd + n_cell
S
S
S
S
S

|||
S
||V--|+---+-----

.
.
.
.

p_a = c0
* ( p(i-1,j ,k ) * xm + p(i+1,j ,k ) * xp
+ p(i ,j-1,k ) * ym + p(i ,j+1,k ) * yp
+ p(i ,j ,k-1) * zm + p(i ,j ,k+1) * zp
- rhs * ds * ds )p(i,j,k) = p(i,j,k) + cnst_p * ( p_a - p0 )end doend doend do

図 4.3-54 SOR 法のコード(2) Red Black 法
定義・参照の依存関係が存在することでベクトル化ができない DO ループに対して，依存関係が生じない順序で配列要素をアクセスすることによりベクトル化を可能にする手法の一つに Red Black 法が

— 58 —

SENAC Vol. 45, No. 2（2012.4）ある．図 4.3-55 に 3 次元の Red Black 法の概念図を示す．配列要素を赤・黒と順番に色付けを行い，赤の要素の計算と黒の要素の計算を行うループを分けることにより，ベクトル化を阻害する依存関係を排除できる．

図 4.3-55 Red Black 法の概念図
3 次元の場合，配列のアクセスは以下のように行われる（各次元のサイズが 8 の場合）．
赤のループ：p(1,1,1),p(3,1,1),…,p(7,1,1),p(2,2,1),p(4,2,1),….
黒のループ：p(2,1,1),p(4,1,1),…,p(8,1,1),p(1,2,1),p(3,2,1),….
図 4.3-56 に Red Black 法のコードを示す．最内 DO ループは上記に示した配列のアクセスのとおり
増分 2 で繰り返される．先頭要素から始まるループが赤のループ，後続の先頭要素+1 から始まるルー
プが黒のループに該当する．j，k の DO ループも増分 2 のアクセスになるように IF 文で処理を分岐している．
113 +----->
114 |+---->
115 ||
116 ||
117 ||V--->:133 |||
134 |||
135 |||
136 |||
137 |||
138 |||
139 ||V--140 ||141 ||
142 ||V--->:do k = n_bnd + 1, n_bnd + n_cell
do j = n_bnd + 1, n_bnd + n_cell
if(mod(k,2).eq.0.and.mod(j,2).eq.0.or.
+mod(k,2).eq.1.and.mod(j,2).eq.1) thendo i = n_bnd + 1, n_bnd + n_cell,2

.
.
.
.

+

p_a = c0
* ( p(i-1,j ,k ) * xm + p(i+1,j ,k ) * xp
+ p(i ,j-1,k ) * ym + p(i ,j+1,k ) * yp
+ p(i ,j ,k-1) * zm + p(i ,j ,k+1) * zp
- rhs * ds * ds )p(i,j,k) = p(i,j,k) + cnst_p * ( p_a - p0 )end doelse if(mod(k,2).eq.1.and.mod(j,2).eq.0.or.
mod(k,2).eq.0.and.mod(j,2).eq.1) thendo i = n_bnd + 2, n_bnd + n_cell,2

赤のループの実行

黒のループの実行

図 4.3-56 Red Black 法のコード(3) DO ループの一重化
図 4.3-56 に示したコードにおいて，最内の DO ループがベクトル化の対象になるため，ベクトル長は

— 59 —

スーパーコンピュータ SX-9 の高速化

i の DO ループの半分のループ長になる（増分 2 でアクセスする）ため，十分なベクトル長が得られない．
そこで，2 つのループを含めて 3 重ループを一重化することで，ベクトル長の拡大を図る．この例では，各次元に計算を行わない袖要素を含んでいるため，単純に DO ループの一重化ができない．そこで予め計算を行う要素の場合を「１」，行わない要素の場合を「０」とするマスクテーブルを用意しておき，IF文の判定で真になる場合のみ計算を実行するようにする．図 4.3-57 に DO ループの一重化を行うコードを示す．コンパイラが配列 p の依存関係を解析できないので，コンパイラ指示行 NODEP を指定して，
配列 p のアクセスに重なりがないことを明示する．
277 +----->
278 |
279 |V---->
280 ||
281 ||:
297 ||
298 ||
299 ||
300 ||
301 ||:
304 ||
305 ||
306 |V---307 +-----

do icolor = 1, 2
!cdir nodep(p)do ijk=icolor,max_cell*max_cell*max_cell,2if(matbl(ijk,1,1,icolor).eq.1.or.
+
matbl(ijk,1,1,icolor+2).eq.1 ) then

.
.
.
.

マスクテーブル配列 matbl の値が真の場合，ループの処理が実行される

p_a = c0
* ( p(ijk-1,1,1) * xm + p(ijk+1,
1,1) * xp+ p(ijk-max_cell,1,1) * ym + p(ijk+max_cell,1,1) * yp+ p(ijk-max_cell*max_cell,1,1) * zm + p(ijk+max_cell*max_cell,1,1) * zp
- rhs * ds * ds )p(ijk,1,1) = p(ijk,1,1) + cnst_p * ( p_a - p0 )endif
end doend do

図 4.3-57 マスクテーブルを使用するコード(4) ADB の利用
SOR 法の計算は，4.3.8(1)で示しているように，ある要素の計算に近傍の 6 つの要素を使用するステンシル計算である．各要素は異なる繰り返しにおける計算で再び参照される．ON_ADB 指示行で，これらの要素を格納する配列 p を ADB に乗せることを指示する．これにより，1 回目のロード命令はメモリから行われるが，2 回目以降のロード命令は ADB から行われることになる．この結果，メモリのアクセス
回数を削減することができ，処理効率の向上が期待される．図 4.3-58 に ON_ADB 指示行を指定した
コードを示す．
277 +----->
278 |
279 |
280 |V---->
281 ||
282 ||:
298 ||
299 ||
300 ||
301 ||
302 ||:
305 ||
306 ||
307 |V---308 +-----

do icolor = 1, 2
!cdir nodep(p)!cdir on_adb(p)do ijk=icolor,max_cell*max_cell*max_cell,2if(matbl(ijk,1,1,icolor).eq.1.or.
+
matbl(ijk,1,1,icolor+2).eq.1 ) then

.
.
.
.

p_a = c0
* ( p(ijk-1,1,1) * xm + p(ijk+1,
1,1) * xp+ p(ijk-max_cell,1,1) * ym + p(ijk+max_cell,1,1) * yp+ p(ijk-max_cell*max_cell,1,1) * zm + p(ijk+max_cell*max_cell,1,1) * zp
- rhs * ds * ds )p(ijk,1,1) = p(ijk,1,1) + cnst_p * ( p_a - p0 )endif
end doend do

図 4.3-58 ON_ADB 指示行の挿入

— 60 —

SENAC Vol. 45, No. 2（2012.4）(5) 性能評価
SX-9 1CPU を用いて性能評価を行う．評価を行うプログラムは以下のとおりである．
①

SOR 法を最適化行わず記述したプログラム②

Red Black 法を用いてベクトル化を行ったプログラム③②に対してマスクテーブルを用いてループの一重化を行ったプログラム④③に対して再利用性の高い配列 p を ADB に乗せることを指示したプログラム

表 4.3-3 に，各次元のループ長が 64 の場合の性能特性を示す．

プログラム
①
②
③
④

表 4.3-3 性能特性
ベクトル演算率(%)
平均ベクトル長
51.86
64.0
90.11
54.4
99.59
255.9
99.59
255.9

実効性能(MFLOPS)
189.5
970.1
8,247.3
18,940.6①のプログラムは，ベクトル化できない DO ループ構造であるため，コンパイラは部分的にしかベクトル化を行わず，ベクトル演算率は 51.9%となり 200MFLOPS に満たない性能値となる．②の 3 次元ループのままの Red Black 法のプログラムでは，ベクトル化は行えるものの，平均ベクトル長が短いことと，ベクトルループを外側で何度も繰り返すため，高い実効性能が得られない．③の 3 次元ループに対して
マスクテーブルを用いて一重化したプログラムでは，平均ベクトル長が 255.9 となり，8.3GFLOPS の性
能値となる．さらに④の ADB に再利用性の高い配列要素を乗せるプログラムでは，メモリへのアクセス
回数を削減することができ，③と比較して 2 倍近い性能向上が得られる．
図 4.3-59 に， それぞれのプログラムで次元サイズを変更した場合の性能グラフを示す．すべてのケースにおいて，④が最も高い性能であることが分かる．ベクトル型プロセッサの SX-9 では高いベクトル演算率，十分な平均ベクトル長だけでなく，メモリ負荷を軽減する最適化が有効であることが示される．特に再利用性の高い要素を ADB に載せることで実効性能の向上が期待される．
Red Black 法は，元の SOR 法と計算順序が異なるため，結果に差異が生じる．そのため SOR 法のよ
うな反復計算において，収束までの繰り返し回数が元の SOR 法の計算より増える場合がある．極端に収束回数が増加する場合は，Red Black 法による高速化の効果が相殺される可能性があるので注意されたい．

図 4.3-59 配列サイズと性能値