Go to Qiita Advent Calendar 2024 Top
LoginSignup
8
7

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

More than 5 years have passed since last update.

@yotakagi77(yotakagi77)

スペクトル法をモダンFortranで実装する際の課題

Last updated at Posted at 2019-06-29

はじめに

本記事は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}
$$
であるから，

diff

のように波数空間で $ik$ を掛け合わせることで簡単に微分演算ができる．

これを利用すると，ラプラス演算子で構成されるポアソン方程式は，微分演算を2回行えばよいので，
$$
\Delta \phi_j = Q_j
$$
を解くために，

  1. スカラー$\phi_j$と生成項$Q_j$をフーリエ変換し，$\phi_k$と$Q_k$を得る．
  2. $(ik)\times(ik)=-k^2$で$Q_k$を割り，それを$\phi_k$に代入する．
  3. $\phi_k$を逆フーリエ変換する．
diff

という計算をすればよい．これがスペクトル法が有限差分法や有限体積法に対する優位点の一つであり，一般に流体解析で最も計算時間がかかるポアソン方程式の反復計算を単純な代数演算に置き換えることができる．

ナビエ・ストークス方程式の解法

上述の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でスペクトル法のカーネル部分を実装する際の問題点は，

  1. EQUIVALENCEは廃止事項
  2. 多くのFFTライブラリの内部はクラシックFORTRANなので引数の制約がある
  3. 複素共役の性質を用いたメモリの節約

などがある．

また，超大規模並列計算を最新のスパコン環境で実施するためには既存のFFTライブラリでは不十分な場合があり，FFTライブラリのチューニング(オープンソースの場合)または自前ライブラリの用意が必要となる．

EQUIVALENCEを使わない解決方法

いろいろな方からの情報を超簡単にまとめると，EQUIVALENCEを使わないためには，

  1. Cのポインタを使う
  2. Cray pointerを使う
  3. transferを使う
  4. 野良サブルーチンを介する

という方法があることがわかりました． 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ライブラリでも機会があったら検証してみたいと思います．

参考文献

  1. C. Canuto, M. Y. Hussaini, A. Quarteroni, T. A. Zand, Specral Method in Fluid Dynamics, Springer-Verlag (1988). ↩

8
7
2

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
8
7

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?