NEC FortranとSX-Aurora TSUBASAを使ってみた（その4）

# 概要
SX-Aurora TSUBASA (SXAT)を使える機会を得たので，姫野ベンチマークをやってみました．また，検索すると東北大学とNECが開催した（と思われる）講習会の資料[^sscc_lecture]が出てきたので，資料を参考にプロファイリングをやってみました．

[^sscc_lecture]: https://www.ss.cc.tohoku.ac.jp/lectures/ の講習会資料→令和3年度開催みｇSX-Aurora_TSUBASAにおける性能分析と高速化

- [その1](https://qiita.com/implicit_none/items/8573afae2ac87532b7e9)
- [その2 - fpm, stdlibを使ってみた](https://qiita.com/implicit_none/items/abe5104e980ccc7d92b1)
- [その3 - stream ベンチマークやってみた](https://qiita.com/implicit_none/items/cebda1fbb26920326a41)
- [その4 - 姫野ベンチマークやってみた](https://qiita.com/implicit_none/items/ed280b88083d46044d1e)
- [その5 - ベクトル和](https://qiita.com/implicit_none/items/c0cf96b582bf62ba4d4e)

# 環境
- Vector Engine Type 10C
- nfort 3.5.1
- ifort 2021.5.0

# 姫野ベンチマーク
[姫野ベンチマーク](https://i.riken.jp/supercom/documents/himenobmt/)は，理化学研究所の姫野龍太郎博士が考案したベンチマークです．公式ページでは，下記のように解説されています．

> 姫野ベンチマークは情報基盤センター・センター長の姫野龍太郎氏が非圧縮流体解析コードの性能評価のために考えたものでポアッソン方程式解法をヤコビの反復法で解く場合に主要なループの処理速度を計るものです。

姫野ベンチマークの背景にある理論については，[解説](https://qiita.com/implicit_none/items/679972853029134a5ab5)を参考にしてください．

姫野ベンチマークには，4通りの問題サイズが用意されています．問題サイズはそれぞれ下記の通りです．

|問題サイズ|格子点数|
|:--|:--|
|S|128x64x64|
|M|256x128x128|
|L|512x256x256|
|XL|1024x512x512|

実行すると，MFLOPSと，その値がPentium III比で何倍かを表示します．SXATの演算性能は2.15 TFLOPSなので，MFLOPSで表示すると判りにくいと考え，GFLOPSに変換して表示するようにしました．また，ソースは動作に影響のない範囲でモダナイズしています．

オプションには`-O4 -mparallel`を使いました．環境変数`VE_OMP_NUM_THREADS`は`8`に設定しています．
```console
$ nfort -O4 -mparallel himenoBMT.f90
```

問題サイズS, M, L, XLの順に実行したときのコンソール出力を示します．

```console
$ ./a.out 
  mimax= 129  mjmax= 65  mkmax= 65
   imax= 128   jmax= 64   kmax= 64
Time measurement accuracy :  0.10000E-05
  Start rehearsal measurement process.
  Measure the performance in 3 times.
   GFLOPS:  89.3538437   time(s):   5.5289268493652344E-04   3.2964125E-03
 Now, start the actual measurement process.
 The loop will be excuted in 325560  times.
 This will take about one minute.
 Wait for a while.
  Loop executed for  325560  times
  Error :   7.2469045E-11
  GFLOPS:  92.0694275   time(s):  58.2302198410034180
  Score based on Pentium III 600MHz :   1.1114126E+03
```

```console
$ ./a.out 
  mimax= 257  mjmax= 129  mkmax= 129
   imax= 256   jmax= 128   kmax= 128
Time measurement accuracy :  0.10000E-05
  Start rehearsal measurement process.
  Measure the performance in 3 times.
   GFLOPS:   1.3701706E+02   time(s):   3.0019283294677734E-03   1.6942389E-03
 Now, start the actual measurement process.
 The loop will be excuted in 59961  times.
 This will take about one minute.
 Wait for a while.
  Loop executed for  59961  times
  Error :   6.0708083E-10
  GFLOPS:   1.3662019E+02   time(s):  60.1738350391387939
  Score based on Pentium III 600MHz :   1.6492056E+03
```

```console
$ ./a.out 
  mimax= 513  mjmax= 257  mkmax= 257
   imax= 512   jmax= 256   kmax= 256
Time measurement accuracy :  0.10000E-05
  Start rehearsal measurement process.
  Measure the performance in 3 times.
   GFLOPS:   1.6466167E+02   time(s):   2.0381927490234375E-02   8.6125336E-04
 Now, start the actual measurement process.
 The loop will be excuted in 8831  times.
 This will take about one minute.
 Wait for a while.
  Loop executed for  8831  times
  Error :   2.2421236E-04
  GFLOPS:   1.6490797E+02   time(s):  59.9079880714416504
  Score based on Pentium III 600MHz :   1.9906805E+03
```

```console
$ ./a.out 
  mimax= 1025  mjmax= 513  mkmax= 513
   imax= 1024   jmax= 512   kmax= 512
Time measurement accuracy :  0.10000E-05
  Start rehearsal measurement process.
  Measure the performance in 3 times.
   GFLOPS:   1.6935811E+02   time(s):   0.1600978374481201   4.3863800E-04
 Now, start the actual measurement process.
 The loop will be excuted in 1124  times.
 This will take about one minute.
 Wait for a while.
  Loop executed for  1124  times
  Error :   3.6194670E-04
  GFLOPS:   1.6997499E+02   time(s):  59.7656250000000000
  Score based on Pentium III 600MHz :   2.0518467E+03
```

|問題サイズ|格子点数|演算性能 [GFLOPS]|
|:--|:--|:--|
|S|128x64x64|92|
|M|256x128x128|137|
|L|512x256x256|165|
|XL|1024x512x512|170|

問題サイズが大きくなるにつれて，GFLOPSの値が上昇していきます．XLでは170 GFLOPSでした．演算性能のカタログスペックとは一桁ほど離れた値ですが，公式には

> 姫野ベンチマークテストは主に計算機のメモリバンド幅などの性能による効果が大きいベンチマークです。

と書いてあるので，演算性能を引き出せる問題ではないようです．二昔ほど前のCPU（Intel Xeon W-2125 @ 4 GHz）で，Intel Fortranを用いてビルドして実行したところ，性能はおおよそ6 GFLOPSでした．これと比較すると，SXATが遅いということはなく，同じコードで30倍ほど高速化しています．異なるプロセッサで，かつ最適化も何もしていないので公平な比較ではないかもしれませんが，バンド幅が影響するような問題に対するSXATの優位性は示されていると判断しています．

# 診断とプロファイリング
講習会の資料には，コンパイラの診断メッセージの取得方法やプロファイリングの方法が書かれていました．

## 診断
診断メッセージを出力するには，`-fdiag-vector=2`や`-report-diagnostics`オプションを付与します．また，`-report-format`を付与すると，ループがベクトル化されたかなど，コンパイラによる最適化結果が出力されます．

`-report-format`オプションを付与すると，結果が`ソースファイル名.L`に出力されます．
```console
$ nfort -O4 -mparallel -report-format himenoBTM.f90 
$ cat himenoBTM.L
中略
   235: +------>          do loop=1,numItr
   236: |                     error= 0.0
   237: |P----->              do k=2,kmax-1
   238: ||U---->                  do j=2,jmax-1
   239: |||V--->                      do i=2,imax-1
   240: ||||    F                         pnew =  a(i,j,k,x )*  p(i+1,j  ,k  ) &
   241: ||||                                     +a(i,j,k,y )*  p(i  ,j+1,k  ) &
   242: ||||                                     +a(i,j,k,z )*  p(i  ,j  ,k+1) &
   243: ||||                                     +b(i,j,k,xy)*( p(i+1,j+1,k  )-p(i+1,j-1,k  )  &
   244: ||||                                                   -p(i-1,j+1,k  )+p(i-1,j-1,k  )) &
   245: ||||                                     +b(i,j,k,yz)*( p(i  ,j+1,k+1)-p(i  ,j-1,k+1)  &
   246: ||||                                                   -p(i  ,j+1,k-1)+p(i  ,j-1,k-1)) &
   247: ||||                                     +b(i,j,k,zx)*( p(i+1,j  ,k+1)-p(i-1,j  ,k+1)  &
   248: ||||                                                   -p(i+1,j  ,k-1)+p(i-1,j  ,k-1)) &
   249: ||||                                     +c(i,j,k,x )*  p(i-1,j  ,k  ) &
   250: ||||                                     +c(i,j,k,y )*  p(i  ,j-1,k  ) &
   251: ||||                                     +c(i,j,k,z )*  p(i  ,j  ,k-1) &
   252: ||||                                     +src(i,j,k)
   253: ||||      
   254: ||||    F                         dp = (pnew*a(i,j,k,ctr) - p(i,j,k))*bnd(i,j,k)
   255: ||||    F                         error = error + dp*dp
   256: ||||    F                         wrk(i,j,k) = p(i,j,k) + rlx*dp
   257: |||V---                       end do
   258: ||U----                   end do
   259: |P-----               end do
   260: |Y=====>              p(2:imax-1,2:jmax-1,2:kmax-1) = wrk(2:imax-1,2:jmax-1,2:kmax-1)
   261: +------           end do
   262:                   !End of iteration
   263:               end subroutine jacobi 
```

この結果はかなり見やすくて，非常によいと思います．`P-P`はループが並列化されたこと，`V-V`はループがベクトル化されたこと，`U-U`はループがアンロールされたことを示しています．`Y=`は，配列式に対して，並列化とベクトル化が行われたことを示しています．
行の先頭にある`F`は，ベクトル積和命令が生成されたことを意味しています．

## プロファイリング
プログラム実行時に，プログラム全体の性能情報を取得できます．

環境変数`VE_PROGINF`を`YES`または`DETAIL`に設定してプログラムを実行すると，取得された情報がプログラム終了後に表示されます．
```console
            ********  Program  Information  ********
  Real Time (sec)                         :             0.327568
  User Time (sec)                         :             2.268602
  Vector Time (sec)                       :             1.932322
  Inst. Count                             :           1395926797
  V. Inst. Count                          :            289407129
  V. Element Count                        :          71896752500
  V. Load Element Count                   :          30120927213
  FLOP Count                              :          36219453027
  MOPS                                    :         27516.939766
  MOPS (Real)                             :        257763.897543
  MFLOPS                                  :         11803.782932
  MFLOPS (Real)                           :        110571.492335
  A. V. Length                            :           248.427718
  V. Op. Ratio (%)                        :            98.689496
  L1 Cache Miss (sec)                     :             0.049389
  VLD LLC Hit Element Ratio (%)           :            33.564993
  Memory Size Used (MB)                   :         16198.000000
  Non Swappable Memory Size Used (MB)     :           100.000000
```

## 簡易性能解析
`-ftrace`オプションを付与してコンパイルし，プログラムを実行すると，性能情報が記録されたファイル`ftrace.out`が出力されます．当該ファイルがあるディレクトリで`ftrace`コマンドを実行すると，性能解析結果が一覧で表示されます．

```console
$ nfort -O4 -mparallel -ftrace himenoBTM.f90
$ ./a.out
$ ls
a.out  ftrace.out  himenoBTM.L  himenoBTM.f90
$ ftrace
*----------------------*
  FTRACE ANALYSIS LIST
*----------------------*

Execution Date : Thu Dec 15 19:46:28 2022 JST
Total CPU Time : 0:00'01"588 (1.588 sec.)


FREQUENCY  EXCLUSIVE       AVER.TIME     MOPS   MFLOPS  V.OP  AVER.    VECTOR L1CACHE CPU PORT VLD LLC PROC.NAME  
           TIME[sec](  % )    [msec]                    RATIO V.LEN      TIME    MISS     CONF HIT E.%            

       32     1.595( 78.7)    49.849  42428.3  22705.5  99.29 255.2     1.503   0.044    0.000   34.79 HIMENOBMTXP_F90::JACOBI$1
        4     0.197(  9.7)    49.366  43011.9  23017.8  99.29 255.2     0.190   0.004    0.000   34.81  -thread0
        4     0.200(  9.9)    50.073  42404.4  22692.7  99.29 255.2     0.189   0.005    0.000   34.73  -thread1
        4     0.201(  9.9)    50.127  42358.8  22668.3  99.29 255.2     0.189   0.005    0.000   34.81  -thread2
        4     0.201(  9.9)    50.137  42350.5  22663.9  99.29 255.2     0.189   0.006    0.000   34.87  -thread3
        4     0.200(  9.9)    50.117  42367.6  22673.1  99.29 255.2     0.188   0.006    0.000   34.76  -thread4
        4     0.200(  9.9)    50.063  42413.2  22697.4  99.29 255.2     0.188   0.006    0.000   34.82  -thread5
        4     0.201(  9.9)    50.159  42331.9  22653.9  99.29 255.2     0.188   0.006    0.000   34.79  -thread6
        4     0.195(  9.6)    48.754  42190.8  22578.4  99.29 255.2     0.183   0.006    0.000   34.74  -thread7

以下略
```

講習会の資料によると，このプロファイリングを利用することで，プログラムがメモリバンド幅依存なのか，LLC（Last Level Cache）のバンド幅依存なのか，演算性能依存なのかを判定することができるようです．

判定には，Byte/FLOPという指標を用います．Byte/FLOPは，1回の浮動小数点演算を行うのに必要なメモリの量（メモリに何 Byte読み書きするか）を表しています．1回の演算をするのに必要なメモリ読み込み量が多いのであれば，いくら演算性能が高くてもメモリバンド幅から読み込む時間より実行時間を短くすることはできません．

ftraceでByte/FLOPを取得するには，環境変数`VE_PERF_MODE`を`VECTOR-MEM`に設定します．再度プログラムを実行してftraceで性能解析結果を確認すると，いくつか出力される項目が増えます．

```console
REQ. REQ.ST REQ.LD ACT.VLD    FLOP COUNT     FMA ELEM. PROC.NAME  
 B/F    B/F    B/F     B/F                                        
3.58   0.27   3.31    1.77   36219451384   11431397040 HIMENOBMTXP_F90::JACOBI$1
3.58   0.27   3.31    1.77    4545186056    1434528256  -thread0
3.58   0.27   3.31    1.77    4545186056    1434528256  -thread1
3.58   0.27   3.31    1.77    4545186056    1434528256  -thread2
3.58   0.27   3.31    1.77    4545186056    1434528256  -thread3
3.58   0.27   3.31    1.77    4545186056    1434528256  -thread4
3.58   0.27   3.31    1.77    4545186056    1434528256  -thread5
3.58   0.27   3.31    1.77    4545186056    1434528256  -thread6
3.58   0.27   3.31    1.77    4403148992    1389699248  -thread7
```

講習会の資料によると，横軸にRequired Byte/FLOP (Req. B/F)，縦軸にActual Byte/FLOP(Act. B/F)をとったとき，Act. B/F=0.51 Req. B/Fの直線上が演算依存，その直線より上がメモリバンド幅依存，下がLLCバンド幅依存だということです．また，Act. B/Fは，ftraceの結果のうちREQ.ST B/F + ACT.VLD B/Fを用いているとのことです．

測定されたAct. B/FはREQ.ST B/F + ACT.VLD B/F=2.04となり，直線0.51 Req. B/F=1.8258の上にあることから，姫野ベンチマークは演算性能依存に近いものの，メモリバンド幅依存であることが判ります．

# まとめ
SX-Aurora TSUBASAで姫野ベンチマークをやってみました．特にソースをいじることなく，コンパイルオプションを付けるだけでかなりの高速化が達成されました．

また，既に公開されている講習会資料を参考に，コンパイル時の診断結果出力，プロファイリングや性能解析を試してみました．著者はベクトル計算機を使ったことがないのですが，経験豊富な方は，これらの情報からチューニングの方針を立てられるのだろうと思います．公開されている資料を探して，試行錯誤してみようと思います．