2.10 非構造格子 CFD ソルバにおけるタイリング
宇宙航空研究開発機構

情報・計算工学センター

坂下 雅秀

2.10.1 はじめに864筐体，82,944CPU（663,552コア）という巨大なシステムで論理ピーク演算性能11.3ペタフロップスを実現した京速コンピュータ「京」（「京」）の共同利用が昨年（2012年）秋に開始された[1]．
さらにその先を見据えて，エクサフロップス時代を実現するための様々な検討も，既にはじまっている[2][3]．
そのような現状において，CFDの分野においても10,000並列を超えるような大規模並列計算への対応が求められている．
一方で，スーパーコンピュータ（スパコン）を構成する要素に目を向けてみると，CPUの演算性能向上にメインメモリのアクセススピード（メモリスループット）が追いつかない状況になりつつある．そのため，メモリに比べて少ない容量ではあるが，十分高速にアクセス可能なキャッシュメモリ（キャッシュ）と呼ばれる機構がCPU内に設けられている．演算に必要なデータのメモリからの供給が十分でなくなりつつある現状では，いかにしてこのキャッシュを有効に利用するかということも，大規模並列計算への対応と並んで，スーパーコンピューティングにおける重要な課題のひとつとなっている．
ここでは，後者のキャッシュ有効利用について，その代表的な手法のひとつであるタイリングを，非構造格子CFDソルバに適用する方法について考える．

2.10.2 演算性能とメモリスループット
まずはじめに，「京」を中心として演算性能とメモリスループットの現状について概観する．
「京」は，SPARC64 VIIIfxと呼ばれるCPUを搭載している．そこで，同じSPARC64系のCPUを搭載したハイエンドテクニカルコンピューティングサーバ FX1（FX1）とPRIMEHPC FX10（FX10）を比較対照として傾向を見ていくこととする．ここに，FX1およびFX10は富士通株式会社（富士通）のHPC（High Performance Computing）用途向け製品である．FX1のCPUは，SPARC64 VIIIfxの前世代にあたるSPARC64 VIIであり，FX10は次の世代となるSPARC64 IXfxである．また，いずれのCPUもopen SPARCアーキテクチャに基づいて富士通が開発したCPUである．表1に，各CPUの主な性能諸元を示す．参考として，2012年8月に発表されたCPUであるSPARC64 Xについても同時に示した．表中のBF比（Byte per FLOP）とは，メモリスループットを論理演算性能で除した値のことであり，演算性能に対してメモリスループットが十分であるか否かを判断するひとつの指標である．

CPU
コ ア 数
ピーク性能[GFLOPS]
メモリスループット[GB/s]
BF 比[－]
キャッシュ[MB]出 荷 開 始 時 期

表1 主な性能諸元[4][5][6] [7]FX1
FX10「京」SPARC64 VII SPARC64 VIIIfx SPARC64 IXfx
481640128236.5
40648510.50.3666122008年
2010年
2012年

―
SPARC64 X
163821020.2724―

表1よりBF比が年々減少していることがわかる．また，図1に1980年の性能を1とした場合のCPUとメモリの性能の推移を示す．この図から，CPUの性能は年率60%程度の性能向上率であるのに対して，メモリの性能向上率は7%程度にしかなっておらず，メモリスループットの問題が，決して今に始まったことではないことがわかる．

CPU

memory

図1 CPUとメモリの性能の推移[8]そこで，以下では圧縮性Navier-Stokes方程式を解く非構造格子CFDソルバ[9]にとって，現状のシステムのBF比がどの程度のものであるか評価することを試みる．
非圧縮性Navier-Stokes方程式を解く場合，計算に必要なデータ量に比べて計算量がもっとも多いのは，数値流束を計算する部分である．したがって，この部分を計算するのに必要なBF比について考える．表2に，FX1において実測した数値流束を計算する際に必要な浮動小数点演算量（FLOP）を示す．表中，HLLEおよびRoeは，数値流束を計算する場合の代表的なスキームである．実測は，1ノード（1cpu,4process）で，格子点数2,416,496のデータを用い，時間積分200ステップで行い，計算に要した時間（CPU time）および実効性能（MFLOPS)を実測し，これより演算量（GFLOP）を求めた．

HLLERoe

表2 数値流束の計算に必要な演算量
CPU time[s]MFLOPS
30.50
4198.5938.70
4852.63

GFLOP
128.07187.80

この結果より，HLLEスキームによる数値流束の計算に必要な1格子点あたりの演算量は，128  109 FLOP  2.4  106 grid  200step  300FLOPとなる．一方で，この計算には，最低限，入力として5次元の保存量ベクトル2本と，出力として，やはり5次元の数値流束ベクトル1本が必要となるので，これをバイト数に換算すれば，

5  2  5word 8byte/word  120byteとなる．ただし，数値流束は倍精度実数型で計算するものとした．この結果，HLLEスキームによる数値流束の計算に必要なBF比は，120byte  300FLOP  0.4と計算できる．この結果からは，FX1のBF比1は数値流束の計算に十分であり，FX10のBF比では不十分ということが言える．
同様に，数値流束をRoeスキームにより計算することにした場合に必要なBF比は，およそ0.27程度であると計算できる．
ここで注意しなければならないのは，CFDソルバにおける数値流束の計算以外の部分は，演算量が相対的に少なくなるため，より大きなBF比を要求することである．このことを考え合わせると，
FX1はともかく，「京」のBF比0.5でもCFDソルバの計算には不十分である．このBF比の不足を補うための機構がキャッシュであり，これをより有効に利用することが求められる時期に来ていることがわかる．

2.10.3 タイリング
キャッシュを有効に利用する手法のひとつにタイリング（Tiling)と呼ばれる手法がある．ここでは，行列・行列積の計算を例に，そのタイリングとはどのようなものかについて見ることにする．
簡単のため，正方行列の場合について考える．今，N行N列の正方行列A=[aij],B=[bij]について，その積C=[cij]C  ABを計算するには，その要素ごとにcij  cij  aik  bkj

i, j , k  1 ~ N(1)を計算することになる．この計算では，入力データとしてaik，bkjおよびcijの3語を必要とし，出力としてcijの1語が必要であるから，合計4語を必要とする．倍精度浮動小数点演算を行うこととすれば，1語は8byteであるから，計算に必要なデータは合計32byteである．一方で，この計算は和（+）と積（×）を含むので，浮動小数点演算数は2である．計算で要求されるBF比は，BF比[Byte/FLOP]=データ量[Byte]÷演算数で求めることができるので，行列・行列積の計算で要求されるBF比は16となる．
ところで，行列・行列積の計算全体に目を向けた場合，ひとつのN行N列の行列に含まれるデータはN2語であり，行列Cは入力と出力の両方に必要であるから，計算に必要なデータ量は4N2語である．
一方，行列・行列積全体を計算するためには，(1)式の計算をi，jおよびkのそれぞれついてN回づつ繰り返す必要があるので，(1)式を計算するのべ回数はN3となる．したがって，行列・行列積全体の計算で要求されるBF比は，



8  4  N 2 2  N 3  16 Nとなる．たとえば，N=1000とすれば，この計算に必要となるBF比は，わずか0.016でよいことになる．
このことは，行列・行列積の計算を以下のように要素の式として書き下してみると容易に理解できる．

c11  a11b11  a12b21   a1N bN 1
c21  a21b11  a22b21   a2 N bN 1

c N 1  a N 1b11  a N 2b21   a NN bN 1(2)

(2)式において，行列Bの要素b11，b12，…はN回繰り返して現れており，データが再利用されていることがわかる．タイリングでは，このことを利用する．
もし，Nが十分小さくて，行列A，BおよびCの全てがキャッシュに乗る程度の大きさであるならば，BF比は16/Nでよいので，特に問題はない．もし，Nが十分に大きくてキャッシュに乗らない場
合には，行列Aをp(=N/n)行q(=N/m)列の小行列Aij=[aijkl]に，行列Bをq(=N/m)行p(=N/n)列の小行列
Bij=[bijkl]に，それぞれ分割する．ただし，nおよびmは適当な整数である．このとき行列・行列積は，
C11 C12…C1n
C21 C22…C2n: ::Cn1 Cn2…Cnn=

A11 A12…A1m
A21 A22…A2m: ::An1 An2…Anm

B11 B12…B1m
B21 B22…B2m: ::Bn1 Bn2…Bnmとなる．ここで，行列Cij=[cijkl]は行列Cのp行p列の小行列であり，pおよびqはキャッシュに乗る程度

の大きさになるように選ぶ．こうすることにより，小行列どうしの計算，
mCij  Cij   A ik B kj

i, j  1 ~ nk 1の一部であるAikBkj，k=kの計算，

cij kl  cij kl  a ij ko  bij ol

k  1 ~ p, o, l  1 ~ q

は，キャッシュ上のデータによって計算可能となり，そのときに必要となるBF比は，



8  2  p 2  pq 2 pq 2で与えられる．簡単のためp=qとなるように小行列を選ぶことにすれば，BF比は16/p（p >1）となり，より狭いメモリバンド幅でキャッシュにデータを供給することが可能となる．
このような計算手法は，行列を小行列（タイル）に分割することから，タイリングと呼ばれる．
ここで，行列・行列積というひとつの計算であっても，(1)式を計算するのに必要なBF比は大きく，演算量は少ないという結果が得られるのに対して，行列・行列積の計算全体を見ると必要とするBF比は小さく，演算量が多いという正反対の結果が得られたことは極めて重要である．このことから，アプリケーションにおけるBF比を評価する場合には，どういう観点から評価するのかに十分な注意を払わなければならないことがわかる．
以下に，行列・行列積にタイリングを適用したことによる効果の一例を示す．リスト1は，一般的な行列・行列積のコーディングの例であり，それにタイリングを適用した場合のコーディング例がリスト2である．これらを実際にFX1で逐次（1コア）実行して経過時間と実効性能を測定した結果を表3に示す．測定は，リスト1および2に示されるようにN=4096，n=m=256，p=q=16として行い，実効性能は，2×N3 / tとして計算した．ただし，tは計算に要した経過時間の実測値である．
! N=p*n=NT*NB=256*16=4096
integer, parameter :: NN = 4096 ! NN=Ninteger, parameter :: NT =256 ! NT=n=m

real*8,allocatable :: A(:,:), B(:,:), C(:,:)

integer, parameter :: NB =16 ! NB=p=q

real*8,allocatable :: A(:,:,:,:),B(:,:,:,:),C(:,:,:,:)
allocate( A(NN,NN), B(NN,NN), C(NN,NN) )allocate( A(NB,NB,NT,NT), B(NB,NB,NT,NT),&
C(NB,NB,NT,NT) )(snip)(snip)

call gettod( t_start ) !経過時間計測開始

call gettod( t_start ) !経過時間計測開始
do n=1,NT
do m=1,NT
do l=1,NT

do j=1,NN

do j=1,NB

do k=1,NN

do k=1,NB

do i=1,NN

do i=1,NBC(i,j) = C(i,j) + A(i,k)*B(k,j)

C(i,j,l,n) = C(i,j,l,n) + A(i,k,l,m)*B(k,j,m,n)

end doend doend doend doend doend doend doend doend focall gettod( t_end ) !経過時間計測終了

call gettod( t_end ) !経過時間計測終了(snip)(sinp)

リスト1 一般的な行列・行列積のコーディング

リスト2 タイリングによる行列・行列積のコーディング例

表3 行列・行列積の測定結果
一般的なコーディング（リスト1）
127.891.07―

経過時間[秒]
実効性能[GFLOPS]
倍率[－]

タイリング適用（リスト2）
60.43
2.27
2.13

表3より，タイリングを適用することにより2倍以上の性能向上が得られたことがわかる．pを変化させた場合の測定結果については，「2.10.6 APPENDX」に示した．
ところで，リスト2のコーディング例では，行列を保存する配列A,BおよびCを2次元配列から小行列に対応する4次元配列に変更している．実際のアプリケーションでは，このような配列の変更には，膨大な工数を必要とするため実現困難な場合も多い．そのような場合には，リスト3に示すように作業用配列を用意し，そこにデータをコピーすることでタイリングを行う簡易的な方法も考えられる．
! N=p*n=NT*NB=256*16= 4096
integer, parameter :: NT = 256
integer, parameter :: NB =16

real*8,allocatable :: A(:,:), B(:,:), C(:,:)
real*8 :: AW(NB,NB), BW(NB,NB), CW(NB,NB)! 作業用配列の確保

allocate( A(N,N), B(N,N), C(N,N) )(snip)
call gettod( t_start ) !経過時間計測開始
do l=1,NT*NT*NT
is = mod(l-1,NT)*NB + 1ie = is + NB - 1
ks = mod(( l -1 )/ NT,NT)*NB + 1
ke = ks + NB - 1
js =( l - 1 )/(NT*NT)*NB + 1

je = js + NB - 1
do k=ks,ke ; do i=is,ieAW(i-is+1,k-ks+1) = A(i,k) ! 行列Aを作業用配列にコピー
end do ; end do
do j=js,je ; do k=ks,keBW(k-ks+1,j-js+1) = B(k,j) ! 行列Bを作業用配列にコピー
end do ;end dodo j=1,NB ; do i=1,NB
CW(i,j) = 0.d0
end do ;! 行列Cの作業用配列を初期化

end dodo j=1,NB ; do k=1,NB ;do i=1,NB

CW(i,j) = CW(i,j) + AW(i,k)*BW(k,j)
end do ;end do ; end do

do j=js,je ;do i=is,ie

C(i,j) = CW(i-is+1,j-js+1) ! 計算結果の書き戻し
end do ;end do!(続き）
end docall gettod( t_end ) !経過時間計測終了(snip)

リスト3 作業用配列を用いてタイリングを適用した行列・行列積のコーディング例

このコーディング例を，FX1で逐次実行したときの経過時間は78.51秒，したがって，実効性能は1.75GFLOPSであった．作業用配列へのコピーのオーバーヘッドがあるため，リスト2のコーディン
グ例に比べれば性能は若干低下するが，タイリングによる一定の効果は得らる．
なお，ここで示したコーディング例（リスト2）は，あくまでも例であって，最適なコーディングではない．小行列のサイズは，pとqを等しくとるよりも，p > qとした方がよい可能性が考えられる
し，その場合には，行列Bは転置行列として配列上に確保した方がよいかも知れない．さらに，当然のことではあるが，最適なpとqの組み合わせは，実行するハードウエアによって異なる．
また，タイリングはあくまでもメモリスループットの能力不足を補い，メインメモリからCPU内のキャッシュにデータを効率的に供給するための手法である．よって，キャッシュから演算器へのデータ供給を含めた演算器そのもの利用については，別途効率化を図る必要がある．

2.10.4 非構造格子CFDソルバへのタイリングの適用
CFDソルバの開発に係る人の多くは，CFDに関する専門知識を有しているが，スパコンのハードウエアやコンパイラによる最適化手法に関する知識を持っているとは限らない．ここでは，そのような人でも，なるべく追加の知識を必要とせずに，比較的容易にキャッシュの有効利用を図る方法について考える．
キャッシュを有効利用するための基本的な考え方は，キャッシュ上にあるデータを使ってなるべく多くの計算を行うようにすることである．一方で，CFD計算を行うためのプラットホームとしては，スパコンかPCクラスタが一般的であり，その結果MPIによりプロセス並列化されていることが多い．並列化の概略は以下の通りである（図2）．
(1) 解析領域をいくつかのリージョンに分割する．
(2) 隣接格子点の物理量を参照可能にするため「袖」と呼ばれる領域を確保する．
(3) 適宜，「袖」にデータを転送しながら計算を進める．
このリージョンは，一種のタイルと考えることができる．そこで，実際にプロセス数を順次増加させ，領域をより多くのリージョンに細分化したときにキャッシュミス率がどうなるかを実測した
例を図3に示す．これより，プロセス数の増加とともにキャッシュミス率は単調に減少していることがわかる．これは，リージョンに含まれる格子点数が減少することにより，一旦キャッシュにロードされたデータが他のデータをロードする必要からメモリに書き戻される可能性が減少したことによるものであり，タイリングによりキャッシュがより有効に利用されていると考えることができる．
しかし，領域の細分化には問題もある．図3の実測を行ったのと同じ条件のもと，プロセス数の増加とともに計算時間がどのように変化するかを測定し並列化効率を求めた結果を図4に示す．この結果からは，領域の細分化には性能的に限界があることがわかる．これは，プロセス間の負荷分散のアンバランスやプロセス間通信のオーバーヘッドなどの要因による性能低下が，キャッシュの有効利用による性能向上を上回ってしまうことが原因である．このことから，プロセス数を増大させずに処理を細分化する方法を考えなければならないことがわかる．
プロセス数を増大させずに処理を細分化するひとつの方法として，領域分割により生成された各リージョンに対して，再帰的に領域分割を適用し，サブリージョンを生成する方法が考えられる．
すなわち，図2において解析領域をリージョン，リージョンをサブリージョンと読み替えればよい．
適切なサブリージョンの大きさはハードウエアの構成に依存するが，1コア当たり1MBのキャッシュがある場合には，1,000格子点程度が目安であると考えられる．
ただし，データの転送（プロセス間通信）に関しては工夫を要する．今，「勾配の計算」―「勾配制限関数の計算」―「数値流束の計算」という非構造格子CFDソルバにおける典型的な計算の流れを例に考える．図5に一連の計算の流れの概略を示す．通常の領域分割による並列化では，まず計算に必要となる「袖」の部分のデータを転送した後（step1），「勾配の計算」を行う（step2）．この計算結果が次の「勾配制限関数の計算」で使用されるため，全てのリージョンで「勾配の計算」が行われた後，計算結果を「袖」の部分に転送する（step3）．転送が終わったのちに「勾配制限関数の計算」を始めることができる（step4）．「数値流束の計算」（step5）では，「勾配の計算」の計算結果とともに「勾配制限関数の計算」の計算結果も必要となるので，やはり計算を始める前に全てのリージョンで「勾配制限関数の計算」が終了し「袖」の部分へのデータの転送（step6）が終わっていなければならない．
しかし，この方法をそのままサブリージョンにおける計算に応用すると，あるサブリージョンで「勾配の計算」に使用するためにキャッシュに乗せられたデータが，別のサブリージョンの計算をするためにキャッシュ上から消去されてしまい，効率の良いキャッシュ利用法とは言えない．そこ
で，計算と「袖」へのデータ転送を交互に繰り返すのではなく，「袖」を余分に確保して，最初にそこへ必要十分なデータを転送しておくことにすれば，そのサブリージョンに関する計算だけを連続して行うことが可能となる．今の場合，「袖」を1格子分ではなく3格子分用意し，「勾配の計算」
では「袖」の部分を2格子分，「勾配制限関数の計算」では「袖」の部分を1格子分，余分に計算することにすれば良い（図6）．もし，n段階の計算を行いたい場合には，n格子分の「袖」を用意すればよい．

解析領域(1) 領域分割

キャッシュミス率[%]

0.40.30.2

0.10.0

リージョン10

リージョン210203040506070

プロセス数(2) 「袖」の確保

図3 プロセス数とキャッシュミス率の関係(3) データの転送

並列化効率[-]

201510500

図2 MPI並列化の概要

2040
プロセス数

60

図4 プロセス数と並列化効率の関係

80step1 勾配の計算に必要なデータの転送

step1 勾配の計算に必要なデータを3格子分転送

step2 勾配の計算
step2 勾配の計算（「袖」の部分も2格子分計算）

リージョン1

リージョン2

step3 勾配の計算結果の転送

サブリージョン1

サブリージョン2

step3 勾配制限関数の計算（「袖」の部分も1格子分計算）

step4 勾配制限関数の計算
step4 数値流束の計算（通常通り計算）

リージョン1

リージョン2

step5 勾配制限関数の計算結果の転送

図6 サブリージョンにおける計算の流れの概略

step6 数値流束の計算

リージョン1
リージョン2
図5 非構造格子CFDソルバの計算の流れの概略
この（「袖」の部分まで計算に含める）方法は，プロセス間通信を最適化する場合によく用いられる方法である．計算量は増えるが通信回数は減っており，通信量が同じであれば通信のオーバーヘッドが削減される．このことから，計算時間の増加に比べて，通信のオーバーヘッドの削減効果が大きければ全体の計算時間は減少する．非構造格子CFDソルバの場合には，勾配の計算結果を「袖」に転送するデータ量に比べ勾配の計算に必要となる「袖」のデータ量の方が少ないので，転送データ量の削減にもなっている．そして，メモリスループットが十分でないハードウエアの場合には，キャッシュを有効に利用することにより計算時間が大幅に削減されることが期待できる．
実装に際しては，サブリージョンiに含まれる格子点の番号をns,ii 1

ns ,i   m j  1j 0m0  0i  1,  n

から始め，サブリージョン内で連番になるように付け替えるとよい．ただし，mjはサブリージョンjに含まれる格子点の総数であり，nはサブリージョンの数である．この番号付け替えによりサブリージョンを生成したことによる既存プログラムのデータ構造の変更が不要になる．
そして，サブリージョン内のデータにアクセスするためには，各サブリージョンにおける格子点番号の最大値ne,i=ns,i+1－1に加えてne,0=0を保存したテーブルを用意しておけばよい．このとき，サブリージョンiに含まれる格子点の番号は，ne,i-1+1からne,iである．
以上の方法は，MPIによるプロセス並列化の手法を再帰的に適用するものであり，実際にCFDソルバの並列化を行ったことのある人であれば，新たな知識を必要とせずに実装可能であるところに最大の特徴がある．
2.10.5 まとめ「京」を中心とした現状のスパコンに対して，非構造格子CFDソルバの計算で必要となるBF比が確保されているか概算した結果，十分とは言えない状況になりつつあることが確認された．
そのため，BF比の不足を補うための機構であるキャッシュの有効利用法について考察した．CFDソルバがMPIにより領域分割の手法で並列化されていることを前提に，同様の方法を分割されたあとの領域に再帰的に適用することで，キャッシュに乗る程度の細かい領域を作ることが可能となる．
このとき，「袖」を十分な量確保し，必要に応じて袖の部分も計算することにより，他の領域の計算結果を参照することなく計算を進めることが可能となる．これにより，キャッシュ上のデータが他の領域の計算のためにメモリに書き戻されるのを避けることができる．
この方法には，MPIによる領域分割の手法で並列化することができるならば，他の知識を必要とせずにキャッシュの有効利用が計れるという利点もある．
2.10.6 APPENDIX
N行N列（N=4096）の正方行列についての行列・行列積にタイリングを適用（リスト2）し，小行列（タイル）の数n（=m）を適宜変更して，FX1において逐次実行したときの経過時間および実効性
能を表A-1に示す．また，そのときのL1データキャッシュミス率（L1$op miss）およびL2キャッシュミス率（L2$ miss）を表A-2に示す．また，図A-1にL1データキャッシュミス率と経過時間の関係を，
図A-2にL2キャッシュミス率の関係をそれぞれ示す．
図A-2よりn=8（p=512）のときに比べn=16（p=256）のときにL2キャッシュミス率が100分の1以下に減少しており，L2キャッシュが有効に利用されるようになったことがわかる．しかし，この段階では実効性能の向上は見られない．
一方で，図A-1よりn=32（p=128）のときに比べてn=64（p=64）のときにL1$op missが10分の1以下に減少し，L1データキャッシュが有効に利用されるようになったことがわかる．同時に，実効性能も明らかに向上しており，リスト2で示した性能向上の効果は，タイリングによってL1データキャッシュが有効に利用されるようになった結果であることがわかる．
また，L1データキャッシュとL2キャッシュのいずれにおいても，小行列のサイズをさらに小さくしていくと，再びミス率は増加に転じており，キャッシュを効果的に利用するには，タイルの大きさを適切に選択しなければならないこともわかる．
ただし，リスト2で示した小行列のサイズn=256（p=16）の場合には，特異的によい実効性能を示しており，単純にキャッシュの利用効率だけで小行列のサイズを決定することはできないこともわ
かる．

表A-1 行列・行列積（リスト2）の測定結果表

n(=m)

p(=q)

経過時間[秒]

12481632641282565124096
204810245122561286432168128.10125.61125.68123.91
123.49124.57
83.26
95.48
60.43
120.69

表A-2 行列・行列積（リスト2）のキャッシュミス率

実効性能[GFLOPS]
1.07
1.09
1.09
1.11
1.11
1.10
1.65
1.44
2.27
1.14

150

3.0

150

100

2.0

100n(=m)p(=q)

L1$op miss[%]

L2$miss[%]

12
4
8
16
32
64
128
256
512

4096
2048
1024
512
256
128
64
32
168

2.79952.78642.78512.77732.76312.71720.17910.24950.72590.96230.69470.69390.69340.00410.00560.00860.01250.02480.14220.3530

1.0

0.5
50

01

2481632641.0500.00128 256 512

小行列（タイル）の数 n(=m)

0.01248163264128256

小行列（タイル）の数 n(=m)

L1キャッシュミス率（右軸：％）

L2キャッシュミス率（右軸：％）

経過時間

経過時間（左軸：秒）

図A-1 L1キャッシュミス率と経過時間の関係

512（左軸：秒）

図A-2 L2キャッシュミス率と経過時間の関係

また，本稿では配列A，BおよびCについてallocate文による動的確保の例（リスト1）を示したが，これを静的に確保した場合には，FX1のコンパイラ1は，自動的にBLASの行列・行列積ライブラリ
MATMULを呼び出す．しかし，このライブラリの実効性能は，配列のサイズを4096×4096として確保すると，少なくともN=512およびN=4096の場合には，およそ430MFLOPS程度にとどまる．確保す
る配列のサイズを5000×5000とすると，N=512の場合の実効性能は8,525MFLOPSとなるので，配列サ
イズ4096×4096の場合に，MUTMALライブラリの性能が特異的に低下するものと推定される．高度に最適化されたライブラリの利用が可能な場合には，それを利用すべきであることは言うまでもな
いが，ライブラリの性能のこのような不安定性は，その利用を阻害するのに十分なものであり，改善が求められる．
また，配列を静的に確保すると自動的に呼び出されるライブラリが，動的に確保すると呼び出されなくなってしまうことに関しても改善の余地があるものと考える．

1

使用したコンパイラは，Fujitsu Fortran Version 8.1．

参考文献
[1] JHPCI
https://www.hpci-office.jp/[2] 文部科学省『「今後のハイパフォーマンス・コンピューティング技術の研究開発について」の報告書のとりまとめ』 平成23年7月15日
http://www.mext.go.jp/b_menu/houdou/23/07/1308508.htm[3] 東北大学サイバーサイエンスセンター他主催 「戦略的高性能計算システム開発」
http://www.open-supercomputer.org/workshop/[4] ハイエンドテクニカルコンピューティングサーバFX1
http://jp.fujitsu.com/solutions/hpc/products/fx1.htm[5] 世界最速スーパーコンピュータ「京」http://www.nsc.riken.jp/shirutsudoi/8-houkoku/p_01.pdf[6] PRIMEHPC FX10
http://jp.fujitsu.com/solutions/hpc/products/primehpc/[7] SPARC64X:Fujitsu's New Generation 16 Core Processor for the next generation UNIX server
http://jp.fujitsu.com/platform/server/sparcenterprise/event/12/hotchips24/pdf/HotChips24_Fujitsu_presentation.pdf[8] Mucci, P., "Memory Bandwidth and the Performance of Scientific Applications: A Studyof the AMD Opteron Processor" AMD Technical Whitepaper, June 2004.
http://web.eecs.utk.edu/~mucci/latest/mucci_pubs.php[9] Hashimoto, A., et al., "Toward the Fastest Unstructured CFD Code 'FaSTAR'" AIAA paper 2012-1075