はじめに
本記事は2019年5月4日開催のモダンFortran勉強会.f01での質問発表用に準備しました．筆者がたまに使用する数値計算法の一つであるスペクトル法をモダンFortranで実装しようとする際に，少し問題となることをまとめました．
追記
モダンFortran勉強会.f01に参加してアドバイスをもらい，その後の暗黙さんの記事でもアドバイスを諸先輩方からいただきましたので，解決法の一つを試してみました．
スペクトル法とは？
スペクトル法は重み付き残差法の1種で，有限要素法のガラーキン法とは親戚みたいなものである．流体解析でのスペクトル法を本格的に学びたい方は，Spectral Methods in Fluid Dynamics1を一読されることをお勧めする．
高速フーリエ変換
スペクトル法では高速フーリエ変換(Fast Fourier Transform, FFT)を利用する．計算時間の大半をFFTが占めるので，FFTさえ早く計算できさえすればよい．
計算格子は有限個数でしか確保できないから，1次元の座標を
$$
x_j = \frac{2\pi j}{N} \qquad j=0,\dots ,N-1
$$
とすると，速度 $u(x_j)$ の離散フーリエ変換は以下のように定義される．
$$
\tilde{u}_k = \frac{1}{N}\sum_{j=0}^{N-1}u(x_j)e^{-ikx_j}
\qquad -N/2 \le k \le N/2-1
$$
ここで，$N$ は格子点数，$i$ は虚数単位，$k$ は波数である．逆変換は以下のようになる．
$$
u(x_j) = \sum_{k=-N/2}^{N/2-1}\tilde{u}_ke^{ikx_j}
\qquad j = 1,\dots ,N-1
$$
これらの離散フーリエ変換の順変換・逆変換を繰り返して偏微分方程式を解くことがスペクトル法でのメインの計算である．
微分演算
$N$ 個の格子点上での微分演算子を $\mathcal{D}_N$ と表すと，速度 $u$ の微分は
以下のように表せる．
$$
(\mathcal{D}_Nu)_l = \sum_{k=-N/2}^{N/2-1}a_ke^{2ikl\pi/N}
$$
ここで，
$$
a_k = \frac{ik}{N}\sum_{j=0}^{N-1}u_je^{-2ikj\pi/N}
$$
であるから，
のように波数空間で $ik$ を掛け合わせることで簡単に微分演算ができる．
これを利用すると，ラプラス演算子で構成されるポアソン方程式は，微分演算を2回行えばよいので，
$$
\Delta \phi_j = Q_j
$$
を解くために，
- スカラー$\phi_j$と生成項$Q_j$をフーリエ変換し，$\phi_k$と$Q_k$を得る．
- $(ik)\times(ik)=-k^2$で$Q_k$を割り，それを$\phi_k$に代入する．
- $\phi_k$を逆フーリエ変換する．
という計算をすればよい．これがスペクトル法が有限差分法や有限体積法に対する優位点の一つであり，一般に流体解析で最も計算時間がかかるポアソン方程式の反復計算を単純な代数演算に置き換えることができる．
ナビエ・ストークス方程式の解法
上述の1階微分，2階微分ができれば，対流項(非線形項)，粘性項は計算でき，ポアソン方程式を解くことによって得られる圧力から圧力勾配項も求まり，時間積分に関しては他の有限差分法などと同様な手法を用いることができる．非線形項の計算時にはエイリアジング誤差(aliasing error)の除去が問題となるが，本記事では省略する．
FFTライブラリによる実フーリエ変換
微分演算の説明で述べたように，スペクトル法では $\mathrm{FFT}$ または $\mathrm{FFT}^{-1}$ によって物理空間と波数空間を行ったり来たりするが，それぞれの場の変数は物理空間と波数空間の両方のメモリを確保しているわけではなく，共通のメモリ領域に割り付ける．格子点数が $N$ の場合，物理空間の配列大きさは $N$ であるが，波数空間では複素数で表現するために，実部・虚部を合わせて $2N$ の大きさになってしまうが，スペクトル法では複素共役，すなわち
$$
\tilde{u}_k^* = \tilde{u}_{-k}
$$
の性質があるため，例えば，クラシックFORTRANで物理空間で
REAL U(N)
と実数型宣言された変数は，その後で
COMPLEX UC(N/2)
EQUIVALENCE (U,UC)
として複素数型として宣言し，EQUIVALENCEを使うことで同じメモリ領域に割り付けることができる．
このような共役複素数の性質を用いた実装の仕方はスペクトル法ではよく用いられており，科学技術計算ライブラリに含まれるFFTライブラリでもそのような使い方を前提とした仕様となっていることが多い．例えば，筆者がNECのSXシリーズスパコンで利用したことがあるASLライブラリでは1次元だけでなく3次元の実フーリエ変換ライブラリDFR3FB(倍精度実数版)，RFR3FB(単精度実数版)が提供されており，使用方法及びライブラリコール時の引数は以下のようになっている．
CALLL DFR3FB(NX,NY,NZ,R,LX,LY,LZ,ISW,IFAX,TRIGS,WK,IERR)
CALL RFR3FB(NX,NY,NZ,R,LX,LY,LZ,ISW,IFAX,TRIGS,WK,IERR)
| 項番 | 引数名 | 型 | 大きさ | 入出力 | 内容 |
|---|---|---|---|---|---|
| 1 | NX | I | 1 | 入力 | 1次元目のデータ数 $n_x$ |
| 2 | NY | I | 1 | 入力 | 2次元目のデータ数 $n_y$ |
| 3 | NZ | I | 1 | 入力 | 3次元目のデータ数 $n_z$ |
| 4 | R | D/R | LX,LY,LZ | 入力 | 入力データ$r_{kx,ky,kz}$(順変換),または$c_{jx,jy,jz}$(逆変換) |
| 出力 | 出力データ$c_{jx,jy,jz}$(順変換),または$r_{kx,ky,kz}$(逆変換) | ||||
| 5 | LX | I | 1 | 入力 | 配列Rの整合寸法 |
| 6 | LY | I | 1 | 入力 | 配列Rの第2寸法 |
| 7 | LZ | I | 1 | 入力 | 配列Rの第3寸法 |
| 8 | ISW | I | 1 | 入力 | 処理スイッチ(ISW = 0:初期化のみ, ISW = 1:初期化を含む順変換, ISW =−1:初期化を含む逆変換) |
| 9 | IFAX | I | 60 | 出力 | 基数分け情報 |
| 10 | TRIGS | D/R | NX+2×(NY+NZ) | 出力 | 三角関数テーブル |
| 11 | WK | D/R | LX×LY×LZ | ワーク | 作業領域 |
| 12 | IERR | I | 1 | 出力 | エラーインディケータ |
実際のDFR3FBライブラリの使用サンプルコードは以下のとおりである．
コード0: DFFT_code0.f
PROGRAM BFR3BF
C *** EXAMPLE OF DFR3FB AND DFR3BF ***
PARAMETER (NX=6,NY=4,NZ=3,LX=10,LY=5,LZ=3,NT=2)
DOUBLE PRECISION R(LX,LY,LZ)
DOUBLE PRECISION TRIGS(NX+2*(NY+NZ)),WK(LX*LY*LZ)
COMPLEX*16 C(LX/2,LY,LZ)
INTEGER IFAX(60)
EQUIVALENCE(R,C)
C**** INPUT ****
DO K=1,NZ
DO J=1,NY
DO I=1,NX
R(I,J,K)=DBLE(I*J*K)/DBLE(NX*NY*NZ)
ENDDO
ENDDO
ENDDO
WRITE(6,1000)
WRITE(6,1010) NX,NY,NZ,LX,LY,LZ,NT
DO K=1,NZ
WRITE(6,1020) K,((R(I,J,K),J=1,NY),I=1,NX)
ENDDO
C**** OUTPUT ****
WRITE(6,1030)
C**** INITIALIZATION + FORWARD TRANSFORM ****
ISW= 1
CALL DFR3FB(NX,NY,NZ,R,LX,LY,LZ,ISW,IFAX,TRIGS,WK,IERR)
C**** NORMALIZATION ****
DO K=1,NZ
DO J=1,NY
DO I=1,(NX+2)/2
C(I,J,K)=C(I,J,K)/DBLE(NX*NY*NZ)
ENDDO
ENDDO
ENDDO
WRITE(6,1040) IERR
DO K=1,NZ
WRITE(6,1050) K
DO I=1,(NX+2)/2
WRITE(6,1060) (C(I,J,K),J=1,NY)
ENDDO
ENDDO
C**** BACKWARD TRANSFORM ****
ISW=-1
CALL DFR3BF(NX,NY,NZ,R,LX,LY,LZ,ISW,IFAX,TRIGS,WK,IERR)
WRITE(6,1070) IERR
DO K=1,NZ
WRITE(6,1020) K,((R(I,J,K),J=1,NY),I=1,NX+2)
ENDDO
STOP
C**** FORMAT ****
1000 FORMAT(1X,'*** DFR3FB AND DFR3BF **'//
& 1X,' ** INPUT **'/)
1010 FORMAT(1X,' NX =',I3,' NY =',I3,' NZ =',I3/
& 1X,' LX =',I3,' LY =',I3,' LZ =',I3/
& 1X,' NT =',I3/)
1020 FORMAT(/3X,'R(IX,IY,',I2,')'
& /4(4X,F7.4))
1030 FORMAT(/1X,' ** OUTPUT **')
1040 FORMAT(/1X,' ( FORWARD TRANSFORM )'//4X,'IERR =',I6)
1050 FORMAT(/3X,'C(IX,IY,',I2,') ')
1060 FORMAT(3X,4(' (',F6.3,',',F6.3,')'))
1070 FORMAT(/1X,' ( BACKWARD TRANSFORM )'//4X,'IERR =',I6)
END
実際にモダンFortranで実装するには??
モダンFortranでスペクトル法のカーネル部分を実装する際の問題点は，
-
EQUIVALENCEは廃止事項 - 多くのFFTライブラリの内部はクラシックFORTRANなので引数の制約がある
- 複素共役の性質を用いたメモリの節約
などがある．
また，超大規模並列計算を最新のスパコン環境で実施するためには既存のFFTライブラリでは不十分な場合があり，FFTライブラリのチューニング(オープンソースの場合)または自前ライブラリの用意が必要となる．
EQUIVALENCEを使わない解決方法
いろいろな方からの情報を超簡単にまとめると，EQUIVALENCEを使わないためには，
- Cのポインタを使う
- Cray pointerを使う
-
transferを使う - 野良サブルーチンを介する
という方法があることがわかりました． 1.のCのポインタを実際に使ってみました．
固定形式から自由形式への変更など
EQUIVALENCEを取り替える前に，前述のコード0のみてくれをモダンFortranぽくしてみました．
コード1: DFFT_code1.f90
program bfr3bf
use,intrinsic :: iso_fortran_env
implicit none
! *** example of dfr3fb and dfr3bf ***
integer(int32), parameter :: nx=6,ny=4,nz=3,lx=10,ly=5,lz=3,nt=2
real(real64) :: r(lx,ly,lz)
real(real64) :: trigs(nx+2*(ny+nz)),wk(lx*ly*lz)
complex(real64) :: c(lx/2,ly,lz)
integer(int32) :: ifax(60)
integer(int32) :: i,j,k,isw,ierr
equivalence(r,c)
!**** input ****
do k = 1,nz
do j = 1,ny
do i = 1,nx
r(i,j,k) = dble(i*j*k)/dble(nx*ny*nz)
enddo
enddo
enddo
print 1000
print 1010,nx,ny,nz,lx,ly,lz,nt
do k = 1,nz
print 1020,k,((r(i,j,k),j=1,ny),i=1,nx)
enddo
!**** output ****
print 1030
!**** initialization + forward transform ****
isw = 1
call dfr3fb(nx,ny,nz,r,lx,ly,lz,isw,ifax,trigs,wk,ierr)
!**** normalization ****
do k = 1,nz
do j = 1,ny
do i = 1,(nx+2)/2
c(i,j,k) = c(i,j,k)/dble(nx*ny*nz)
enddo
enddo
enddo
print 1040,ierr
do k = 1,nz
print 1050, k
do i = 1,(nx+2)/2
print 1060, (c(i,j,k),j=1,ny)
enddo
enddo
!**** backward transform ****
isw = -1
call dfr3bf(nx,ny,nz,r,lx,ly,lz,isw,ifax,trigs,wk,ierr)
print 1070,ierr
do k = 1,nz
print 1020,k,((r(i,j,k),j=1,ny),i=1,nx+2)
enddo
stop
!**** format ****
1000 format(1x,'*** dfr3fb and dfr3bf **'// &
1x,' ** input **'/)
1010 format(1x,' nx =',i3,' ny =',i3,' nz =',i3/ &
1x,' lx =',i3,' ly =',i3,' lz =',i3/ &
1x,' nt =',i3/)
1020 format(/3x,'r(ix,iy,',i2,')' &
/4(4x,f7.4))
1030 format(/1x,' ** output **')
1040 format(/1x,' ( forward transform )'//4x,'ierr =',i6)
1050 format(/3x,'c(ix,iy,',i2,') ')
1060 format(3x,4(' (',f6.3,',',f6.3,')'))
1070 format(/1x,' ( backward transform )'//4x,'ierr =',i6)
end program bfr3bf
format文の置き換えは面倒なのでしませんでした．
Cのポインタを導入
次にCのポインタを使ってequivalenceを削除します．
コード2: DFFT_code2.f90
program bfr3bf
use,intrinsic :: iso_fortran_env
use,intrinsic :: iso_c_binding
implicit none
! *** example of dfr3fb and dfr3bf ***
integer(int32), parameter :: nx=6,ny=4,nz=3,lx=10,ly=5,lz=3,nt=2
! real(real64) :: r(lx,ly,lz)
real(real64),target :: r(lx,ly,lz)
real(real64) :: trigs(nx+2*(ny+nz)),wk(lx*ly*lz)
! complex(real64) :: c(lx/2,ly,lz)
complex(real64),dimension(:,:,:),contiguous,pointer :: c
integer(int32) :: ifax(60)
integer(int32) :: i,j,k,isw,ierr
! equivalence(r,c)
call c_f_pointer(c_loc(r),c,[lx/2,ly,lz])
!**** input ****
do k = 1,nz
do j = 1,ny
do i = 1,nx
r(i,j,k) = dble(i*j*k)/dble(nx*ny*nz)
enddo
enddo
enddo
print 1000
print 1010,nx,ny,nz,lx,ly,lz,nt
do k = 1,nz
print 1020,k,((r(i,j,k),j=1,ny),i=1,nx)
enddo
!**** output ****
print 1030
!**** initialization + forward transform ****
isw = 1
call dfr3fb(nx,ny,nz,r,lx,ly,lz,isw,ifax,trigs,wk,ierr)
!**** normalization ****
do k = 1,nz
do j = 1,ny
do i = 1,(nx+2)/2
c(i,j,k) = c(i,j,k)/dble(nx*ny*nz)
enddo
enddo
enddo
print 1040,ierr
do k = 1,nz
print 1050, k
do i = 1,(nx+2)/2
print 1060, (c(i,j,k),j=1,ny)
enddo
enddo
!**** backward transform ****
isw = -1
call dfr3bf(nx,ny,nz,r,lx,ly,lz,isw,ifax,trigs,wk,ierr)
print 1070,ierr
do k = 1,nz
print 1020,k,((r(i,j,k),j=1,ny),i=1,nx+2)
enddo
stop
!**** format ****
1000 format(1x,'*** dfr3fb and dfr3bf **'// &
1x,' ** input **'/)
1010 format(1x,' nx =',i3,' ny =',i3,' nz =',i3/ &
1x,' lx =',i3,' ly =',i3,' lz =',i3/ &
1x,' nt =',i3/)
1020 format(/3x,'r(ix,iy,',i2,')' &
/4(4x,f7.4))
1030 format(/1x,' ** output **')
1040 format(/1x,' ( forward transform )'//4x,'ierr =',i6)
1050 format(/3x,'c(ix,iy,',i2,') ')
1060 format(3x,4(' (',f6.3,',',f6.3,')'))
1070 format(/1x,' ( backward transform )'//4x,'ierr =',i6)
end program bfr3bf
コードの変更箇所は，
-
iso_c_bindingのモジュールを使用する． - FFTに用いる実数・複素数の配列宣言の記述を変更する．
- サブルーチン
c_f_pointerを使用する．
です．実際にこのコードをコンパイルして実行すると，正しい結果が得られました．
% sxf03 -lasl DFFT_code2.f90
% ./a.out
*** dfr3fb and dfr3bf **
** input **
nx = 6 ny = 4 nz = 3
lx = 10 ly = 5 lz = 3
nt = 2
r(ix,iy, 1)
0.0139 0.0278 0.0417 0.0556
0.0278 0.0556 0.0833 0.1111
0.0417 0.0833 0.1250 0.1667
0.0556 0.1111 0.1667 0.2222
0.0694 0.1389 0.2083 0.2778
0.0833 0.1667 0.2500 0.3333
r(ix,iy, 2)
0.0278 0.0556 0.0833 0.1111
0.0556 0.1111 0.1667 0.2222
0.0833 0.1667 0.2500 0.3333
0.1111 0.2222 0.3333 0.4444
0.1389 0.2778 0.4167 0.5556
0.1667 0.3333 0.5000 0.6667
r(ix,iy, 3)
0.0417 0.0833 0.1250 0.1667
0.0833 0.1667 0.2500 0.3333
0.1250 0.2500 0.3750 0.5000
0.1667 0.3333 0.5000 0.6667
0.2083 0.4167 0.6250 0.8333
0.2500 0.5000 0.7500 1.0000
** output **
( forward transform )
ierr = 0
c(ix,iy, 1)
( 0.243, 0.000) (-0.049, 0.049) (-0.049, 0.000) (-0.049,-0.049)
(-0.035, 0.060) (-0.005,-0.019) ( 0.007,-0.012) ( 0.019,-0.005)
(-0.035, 0.020) ( 0.003,-0.011) ( 0.007,-0.004) ( 0.011, 0.003)
(-0.035, 0.000) ( 0.007,-0.007) ( 0.007, 0.000) ( 0.007, 0.007)
c(ix,iy, 2)
(-0.061, 0.035) ( 0.005,-0.019) ( 0.012,-0.007) ( 0.019, 0.005)
( 0.000,-0.020) ( 0.004, 0.004) (-0.000, 0.004) (-0.004, 0.004)
( 0.006,-0.010) ( 0.001, 0.003) (-0.001, 0.002) (-0.003, 0.001)
( 0.009,-0.005) (-0.001, 0.003) (-0.002, 0.001) (-0.003,-0.001)
c(ix,iy, 3)
(-0.061,-0.035) ( 0.019,-0.005) ( 0.012, 0.007) ( 0.005, 0.019)
( 0.017,-0.010) (-0.001, 0.005) (-0.003, 0.002) (-0.005,-0.001)
( 0.012, 0.000) (-0.002, 0.002) (-0.002, 0.000) (-0.002,-0.002)
( 0.009, 0.005) (-0.003, 0.001) (-0.002,-0.001) (-0.001,-0.003)
( backward transform )
ierr = 0
r(ix,iy, 1)
0.0139 0.0278 0.0417 0.0556
0.0278 0.0556 0.0833 0.1111
0.0417 0.0833 0.1250 0.1667
0.0556 0.1111 0.1667 0.2222
0.0694 0.1389 0.2083 0.2778
0.0833 0.1667 0.2500 0.3333
0.0000 0.0000 0.0000 0.0000
0.0000 0.0000 0.0000 0.0000
r(ix,iy, 2)
0.0278 0.0556 0.0833 0.1111
0.0556 0.1111 0.1667 0.2222
0.0833 0.1667 0.2500 0.3333
0.1111 0.2222 0.3333 0.4444
0.1389 0.2778 0.4167 0.5556
0.1667 0.3333 0.5000 0.6667
0.0000 0.0000 0.0000 0.0000
0.0000 0.0000 0.0000 0.0000
r(ix,iy, 3)
0.0417 0.0833 0.1250 0.1667
0.0833 0.1667 0.2500 0.3333
0.1250 0.2500 0.3750 0.5000
0.1667 0.3333 0.5000 0.6667
0.2083 0.4167 0.6250 0.8333
0.2500 0.5000 0.7500 1.0000
0.0000 0.0000 0.0000 0.0000
0.0000 0.0000 0.0000 0.0000
出力を見ると，順変換・逆変換をそれぞれ1回行い，元の値に戻っていることがわかります．このプログラムを実行時に，
* 253 Invalid operation FUNC=zwcmf3 ELN=247(4000b0e58)
というエラーが出たのですが，原因はよくわからず，特に悪さをしているわけでもないので，とりあえずはスルーしておきます．
大きい配列計算での性能評価
上記のプログラムでは，配列の大きさが $L_x\times L_y\times L_z = 10\times 5\times 3 = 150$ と小さすぎるので，スペクトル法で使用するときの性能については予測できません．そこで3次元配列の大きさを $N = 256^3$，$512^3$，$1024^3$にして計算してみました．コードはほとんど変わりませんが，以下のようになります．
コード3: DFFT_code3.f90 ($N=256^3$，Cポインタ使用版)
program bfr3bf
use,intrinsic :: iso_fortran_env
use,intrinsic :: iso_c_binding
implicit none
! *** example of dfr3fb and dfr3bf ***
! integer(int32), parameter :: nx=6,ny=4,nz=3,lx=10,ly=5,lz=3,nt=2
integer(int32), parameter :: nx=256,ny=256,nz=256
integer(int32), parameter :: lx=258,ly=256,lz=256,nt=2
real(real64),target :: r(lx,ly,lz)
real(real64) :: trigs(nx+2*(ny+nz)),wk(lx*ly*lz)
complex(real64),dimension(:,:,:),contiguous,pointer :: c
integer(int32) :: ifax(60)
integer(int32) :: i,j,k,isw,ierr
call c_f_pointer(c_loc(r),c,[lx/2,ly,lz])
!**** input ****
do k = 1,nz
do j = 1,ny
do i = 1,nx
r(i,j,k) = dble(i*j*k)/dble(nx*ny*nz)
enddo
enddo
enddo
print 1000
print 1010,nx,ny,nz,lx,ly,lz,nt
! do k = 1,nz
! print 1020,k,((r(i,j,k),j=1,ny),i=1,nx)
! enddo
!**** output ****
print 1030
!**** initialization + forward transform ****
isw = 1
call dfr3fb(nx,ny,nz,r,lx,ly,lz,isw,ifax,trigs,wk,ierr)
!**** normalization ****
do k = 1,nz
do j = 1,ny
do i = 1,(nx+2)/2
c(i,j,k) = c(i,j,k)/dble(nx*ny*nz)
enddo
enddo
enddo
print 1040,ierr
! do k = 1,nz
! print 1050, k
! do i = 1,(nx+2)/2
! print 1060, (c(i,j,k),j=1,ny)
! enddo
! enddo
!**** backward transform ****
isw = -1
call dfr3bf(nx,ny,nz,r,lx,ly,lz,isw,ifax,trigs,wk,ierr)
print 1070,ierr
! do k = 1,nz
! print 1020,k,((r(i,j,k),j=1,ny),i=1,nx+2)
! enddo
stop
!**** format ****
1000 format(1x,'*** dfr3fb and dfr3bf **'// &
1x,' ** input **'/)
1010 format(1x,' nx =',i3,' ny =',i3,' nz =',i3/ &
1x,' lx =',i3,' ly =',i3,' lz =',i3/ &
1x,' nt =',i3/)
1020 format(/3x,'r(ix,iy,',i2,')' &
/4(4x,f7.4))
1030 format(/1x,' ** output **')
1040 format(/1x,' ( forward transform )'//4x,'ierr =',i6)
1050 format(/3x,'c(ix,iy,',i2,') ')
1060 format(3x,4(' (',f6.3,',',f6.3,')'))
1070 format(/1x,' ( backward transform )'//4x,'ierr =',i6)
end program bfr3bf
計算結果確認用の標準出力(print文)は配列の中身を全て出力して大変なことになるので，コメントアウトしました．また，ここには掲載しませんが，equivalence使用版のプログラムも用意して計算時間等を比較しました．性能評価ツールはいろいろあると思いますが，今回使用している計算機システムではPROGINFとftraceを利用しました．コンパイル，実行は以下のとおりです．
% setenv F_PROGINF DETAIL
% sxf03 -lasl -ftrace DFFT_code3.f90
% ./a.out > log
% sxftrace > ftrace.log
$N = 256^3$のPROGINFの結果を以下に示します．
equivalence 版の結果
****** Program Information ******
Real Time (sec) : 0.100999
User Time (sec) : 0.093046
Sys Time (sec) : 0.005181
Vector Time (sec) : 0.092248
Inst. Count : 39755068.
V. Inst. Count : 13127702.
V. Element Count : 3273784614.
V. Load Element Count : 506582455.
FLOP Count : 1388180405.
MOPS : 35470.756185
MFLOPS : 14919.291587
A. V. Length : 249.379870
V. Op. Ratio (%) : 99.193211
Memory Size (MB) : 448.031250
MIPS : 427.262515
I-Cache (sec) : 0.000099
O-Cache (sec) : 0.000291
Bank Conflict Time
CPU Port Conf. (sec) : 0.002544
Memory Network Conf. (sec) : 0.037090
ADB Hit Element Ratio (%) : 22.639856
Cポインタ版の結果
****** Program Information ******
Real Time (sec) : 0.100042
User Time (sec) : 0.091974
Sys Time (sec) : 0.005187
Vector Time (sec) : 0.091267
Inst. Count : 40013091.
V. Inst. Count : 13265437.
V. Element Count : 3260941350.
V. Load Element Count : 493991607.
FLOP Count : 1388180405.
MOPS : 35745.852132
MFLOPS : 15093.182910
A. V. Length : 245.822384
V. Op. Ratio (%) : 99.186430
Memory Size (MB) : 448.031250
MIPS : 435.047850
I-Cache (sec) : 0.000094
O-Cache (sec) : 0.000241
Bank Conflict Time
CPU Port Conf. (sec) : 0.002595
Memory Network Conf. (sec) : 0.037046
ADB Hit Element Ratio (%) : 18.103757
メモリサイズは同一であり，きちんとメモリの割り付けができていることがわかります．ftraceはほぼ同じ情報を出力するので，全てのケースについてftraceの結果をまとめたものを以下に示します．
*----------------------*
FTRACE ANALYSIS LIST
*----------------------*
*** N = 256^3 ***
PROC.NAME FREQUENCY EXCLUSIVE AVER.TIME MOPS MFLOPS V.OP AVER. VECTOR I-CACHE O-CACHE BANK CONFLICT ADB HIT
TIME[sec]( % ) [msec] RATIO V.LEN TIME MISS MISS CPU PORT NETWORK ELEM.%
(equivalence)
BFR3BF 1 0.093(100.0) 92.897 35527.2 14943.3 99.19 249.4 0.092 0.000 0.000 0.003 0.037 22.64
----------------------------------------------------------------------------------------------------------------------------------
total 1 0.093(100.0) 92.897 35527.2 14943.3 99.19 249.4 0.092 0.000 0.000 0.003 0.037 22.64
(C-pointer)
BFR3BF 1 0.092(100.0) 91.824 35803.7 15117.9 99.19 245.8 0.091 0.000 0.000 0.003 0.037 18.10
----------------------------------------------------------------------------------------------------------------------------------
total 1 0.092(100.0) 91.824 35803.7 15117.9 99.19 245.8 0.091 0.000 0.000 0.003 0.037 18.10
*** N = 512^3 ***
PROC.NAME FREQUENCY EXCLUSIVE AVER.TIME MOPS MFLOPS V.OP AVER. VECTOR I-CACHE O-CACHE BANK CONFLICT ADB HIT
TIME[sec]( % ) [msec] RATIO V.LEN TIME MISS MISS CPU PORT NETWORK ELEM.%
(equivalence)
BFR3BF 1 0.780(100.0) 779.561 37593.6 16154.6 99.35 252.5 0.777 0.000 0.001 0.019 0.301 15.38
----------------------------------------------------------------------------------------------------------------------------------
total 1 0.780(100.0) 779.561 37593.6 16154.6 99.35 252.5 0.777 0.000 0.001 0.019 0.301 15.38
(C-pointer)
BFR3BF 1 0.776(100.0) 776.106 37870.8 16226.6 99.34 249.1 0.775 0.000 0.000 0.019 0.302 11.80
----------------------------------------------------------------------------------------------------------------------------------
total 1 0.776(100.0) 776.106 37870.8 16226.6 99.34 249.1 0.775 0.000 0.000 0.019 0.302 11.80
*** N = 1024^3 ***
PROC.NAME FREQUENCY EXCLUSIVE AVER.TIME MOPS MFLOPS V.OP AVER. VECTOR I-CACHE O-CACHE BANK CONFLICT ADB HIT
TIME[sec]( % ) [msec] RATIO V.LEN TIME MISS MISS CPU PORT NETWORK ELEM.%
(equivalence)
BFR3BF 1 6.395(100.0) 6395.043 38721.5 18112.3 99.46 253.8 6.384 0.000 0.005 0.127 2.566 17.37
----------------------------------------------------------------------------------------------------------------------------------
total 1 6.395(100.0) 6395.043 38721.5 18112.3 99.46 253.8 6.384 0.000 0.005 0.127 2.566 17.37
(C-pointer)
BFR3BF 1 6.288(100.0) 6288.001 39519.9 18420.6 99.46 252.2 6.281 0.000 0.002 0.134 2.527 14.20
----------------------------------------------------------------------------------------------------------------------------------
total 1 6.288(100.0) 6288.001 39519.9 18420.6 99.46 252.2 6.281 0.000 0.002 0.134 2.527 14.20
結果を比較すると，どの配列大きさでもCポインタの方がFLOPSが少し高く，ADB(Assignable Data Buffer)のヒット率はequivalenceの方が高くなっています．ですが，差は微少でCポインタにすることによって速度が遅くなったわけではないので，equivalenceを使わないで実フーリエ変換を利用し，スペクトル法をモダンFortranで実装できそうです．
まとめ
EQUIVALENCEを使わない方法の1つとしてCポインタを利用する方法を試してみて，特に問題なく代替できることがわかりました．教えていただいた他の方法を試していませんが，Cポインタを利用する方法は(私としては)特に違和感なく利用できるので今回は良しとします．違う計算機システムや違うコンパイラ，違うFFTライブラリでも機会があったら検証してみたいと思います．
参考文献
-
C. Canuto, M. Y. Hussaini, A. Quarteroni, T. A. Zand, Specral Method in Fluid Dynamics, Springer-Verlag (1988). ↩