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

###　はじめに

https://qiita.com/soybean/items/f732cb71ff978a2352f8

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


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

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

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


```math 
\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コア数の半分くらい速くなったら良いなと期待していました．この期待は裏切られることになりました．

### 方法

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

```Fortran
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の指示文を入れます．

```Fortran
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で高速化して，オチを付けます．