OpenMPを使って昨日のスペクトル法のコードを高速化してみた by Fortran

###　はじめに

こちらの記事の続きです．

逆離散フーリエ変換のループが重い問題

スペクトル法により一次線形移流方程式を実装する場合，以下のような解法をとることとなります．

(1) 初期値を離散フーリエ変換して，フーリエスペクトルを求める．
(2) 時間発展方程式

\hat{u}_k(t) = \hat{u}_k(0) e^{-ikt}

により未来のフーリエスペクトルを求める．
(3) 逆離散フーリエ変換により，u(x,t)を求める．

ところが，この(1)および(3)は，波数Nについて約O(N^2)オーダーの計算量となり，そのままでは他の解法と比較して計算コストが大きく，辛いです．

そこで，通常は高速フーリエ変換(Fast Fourier Transform, FFT)とよばれるアルゴリズムを用いることで，O(NlogN)オーダーまで計算量を削減します．

以下では，このうち，特に計算量の大きい(3)に着目し，敢えてFFTを用いず，OpenMPによる並列化を軽率に用いた場合，どの程度高速化されるのかを実験してみました．

OpenMPとは

https://openmp.org
によれば，OpenMPとは，マルチプラットフォームな共有メモリ型並列計算プログラミングをサポートするAPIであり，C/C++/Fortranで用いることができます．

Fortranでは，!$OMPで始まる指示文をコード中に入れることで，指示分の内容に応じて並列処理が実行されます．

今回の(3)は，共有メモリを参照し，かつ衝突のない問題となっていますので，CPUコア数の半分くらい速くなったら良いなと期待していました．この期待は裏切られることになりました．

方法

以下のコードを題材とします．

program main

use plantfem
implicit none

integer(int32) ,parameter:: N = 2**15
complex(real64),allocatable :: x(:),uk0(:),uk(:),ut(:)
integer(int32) :: i,k,j
real(real64) :: T


type(MPI_) :: mpid
type(Math_) :: Math
type(IO_)   :: f
type(Console_) :: console


call mpid%start()


! 初期時刻を0として，何秒後の解がほしいか．
!call console%log("Input Time: T")
!call console%read()
!T = freal(console%line)
T = 1.0d0

! (1) 離散フーリエ変換により，初期値のフーリエスペクトルを求める．
allocate(uk0(-N:N) )

do k=-N,N
    if(k==0)then
        uk0(k) = 1.0d0/Math%PI*exp(-Math%i*dble(k)*Math%PI )*Math%PI/5.0d0
        cycle  
    endif
    uk0(k) =1.0d0/Math%PI/dble(k)*exp(-Math%i*dble(k)*Math%PI )*sin(Math%PI/5.0d0*dble(k) ) 
enddo


! (2) ある時刻T後のフーリエスペクトルを求める．
allocate(uk(-N:N) )
do k=-N,N
    uk(k) = uk0(k)*exp(-Math%i*dble(k)*T )
enddo

! (3) 逆フーリエ変換(by DFT)
x = linspace([0.0d0,2.0d0*Math%PI],2*N+1)
ut = zeros(2*N+1)

!>> ここをターゲットにします．
!>> ここをターゲットにします．

do i=1,size(x)
    do k=-N,N
        ut(i) = ut(i) + uk(k)*exp(Math%i*dble(k)*x(i) )
    enddo
enddo
!<< ここをターゲットにします．
!<< ここをターゲットにします．


call mpid%end()

end program

この，(3)の処理の箇所に，OpenMPの指示文を入れます．

program main

use plantfem
implicit none

integer(int32) ,parameter:: N = 2**15
complex(real64),allocatable :: x(:),uk0(:),uk(:),ut(:)
integer(int32) :: i,k,j
real(real64) :: T


type(MPI_) :: mpid
type(Math_) :: Math
type(IO_)   :: f
type(Console_) :: console


call mpid%start()


! 初期時刻を0として，何秒後の解がほしいか．
!call console%log("Input Time: T")
!call console%read()
!T = freal(console%line)
T = 1.0d0

! (1) 離散フーリエ変換により，初期値のフーリエスペクトルを求める．
allocate(uk0(-N:N) )

do k=-N,N
    if(k==0)then
        uk0(k) = 1.0d0/Math%PI*exp(-Math%i*dble(k)*Math%PI )*Math%PI/5.0d0
        cycle  
    endif
    uk0(k) =1.0d0/Math%PI/dble(k)*exp(-Math%i*dble(k)*Math%PI )*sin(Math%PI/5.0d0*dble(k) ) 
enddo


! (2) ある時刻T後のフーリエスペクトルを求める．
allocate(uk(-N:N) )
do k=-N,N
    uk(k) = uk0(k)*exp(-Math%i*dble(k)*T )
enddo

! (3) 逆フーリエ変換(by DFT)
x = linspace([0.0d0,2.0d0*Math%PI],2*N+1)
ut = zeros(2*N+1)

!>> ここをターゲットにします．
!>> ここをターゲットにします．
!$OMP parallel do private(i)
do i=1,size(x)
    do k=-N,N
        ut(i) = ut(i) + uk(k)*exp(Math%i*dble(k)*x(i) )
    enddo
enddo
!$OMP end parallel do
!<< ここをターゲットにします．
!<< ここをターゲットにします．


call mpid%end()

end program

ここに，!$OMP parallel do private(i)は，イテレータ(i)をそれぞれのプロセス固有の値として保持して，ループ内を並列化することを指示しています．16
計測には，以下の4種類の環境/条件を用いました．

• (1) Intel(R) Core(TM) i7-10510U CPU @ 1.80GHz (非並列化)
• (2) Intel(R) Core(TM) i9-9900 CPU @ 3.10GHz 4C8T　(非並列化)
• (3) Intel(R) Core(TM) i7-10510U CPU @ 1.80GHz (並列化)
• (4) Intel(R) Core(TM) i9-9900 CPU @ 3.10GHz　8C16T(並列化)

計測結果

N	(1) sec.	(2) sec.	(3) sec.	(4) sec.
2^10	0.66		0.48
2^11	1.34		0.77
2^12	4.08		1.75
2^13	11.23		5.09
2^14	40.45	34.22	12.61	3.76
2^15		137.14		14.48
2^16		547.30		58.79

4C8Tで5倍くらい速くなりました．わーい！
8C16Tで10倍くらい速くなりました．わーい！

まとめと補足

OpenMPでもちゃんとコア数分早くなるもんなんだなと感動しました．
また，do concurrent構文を使ってみたりもしましたが，そちらはむしろ遅くなりました．残念．
次回はFFTで高速化して，オチを付けます．