概要
SX-Aurora TSUBASA (SXAT)で効率的に動作するプログラムの作り方を勉強するために，簡単な問題（ベクトル和）を計算してみました．
環境
- Vector Engine Type 10C
- nfort 3.5.1
ベクトル和
ベクトル（プログラム上では配列）$A, B$の和$C = A + B$を計算してみます．特に難しいことは何もなく，これをやっても新しいことは判りません．ですが，既存のベンチマーク問題をやって性能がこうなりましたと言われても，「左様でございますか」としか言えないので，自身でプログラムを作る中で現れる個々の計算（ベクトル和，総和など）がどう扱われるのか（最適化されるのか）を勉強する意味はあるかと思います．
ベクトル和の計算には，Fortranにおける実行時間の測定で作成したプログラムを修正して（ループのi, jの順序を入れ替えて）利用しました．実行時間の測定には，組込手続cpu_timeを用いました．
program main
use, intrinsic :: iso_fortran_env
!$ use omp_lib
implicit none
integer(int32), parameter :: N = 2**13
real(real64), allocatable :: a(:, :), b(:, :), c(:, :)
integer(int32) :: i, j
real(real32) :: time_begin_s, time_end_s
allocate (a(N, N))
allocate (b(N, N))
allocate (c(N, N))
!$omp parallel
! 配列の初期化
!$omp do
do i = 1, N
do j = 1, N
a(i, j) = 1._real64
b(i, j) = 2._real64
c(i, j) = 0._real64
end do
end do
!$omp end do
!$omp master
call cpu_time(time_begin_s)
!$omp end master
!$omp do
do j = 1, N
do i = 1, N
c(i, j) = a(i, j) + b(i, j)
end do
end do
!$omp end do
!$omp master
call cpu_time(time_end_s)
print *, time_end_s - time_begin_s, "sec", sum(c)/N/N
!$omp end master
!$omp end parallel
deallocate (a)
deallocate (b)
deallocate (c)
end program main
OpenMP並列化
最初は，OpenMPを利用して並列化し，スレッド数を替えて実行時間を測定しました．
$ nfort -fopenmp -report-format vectoradd.f90
診断結果を見ると，外側ループが並列化され，内側ループがベクトル化されていることを確認できます．
31: !$omp do
32: P------> do j = 1, N
33: |V-----> do i = 1, N
34: || c(i, j) = a(i, j) + b(i, j)
35: |V----- end do
36: P------ end do
37: !$omp end do
実行結果は下記のようになりました．
| スレッド数 | 実行時間 [ms] |
|---|---|
| 1 | 4.66 |
| 2 | 2.71 |
| 4 | 2.24 |
| 8 | 2.29 |
stream ベンチマークの時と同様に，4スレッドあたりで実行時間が最短になり，8スレッドになると若干遅くなりました．
最適化オプション-O3や-O4を付与しても，実行時間には影響がほとんどありませんでした．
$ nfort -O4 -fopenmp -report-format vectoradd.f90
一方で，配列Cの平均の計算は，最適化オプションを付けることで，ループのベクトル化だけでなく，多重ループの1重化も行われるようになりました．平均の計算が高速化されたことで，プログラムが早く終了するようになったことを，体感でも感じることができます．
$ nfort -fopenmp -report-format vectoradd.f90
$ cat vectoradd.L
中略
40: V======> print *,time_end_s - time_begin_s,"sec",sum(c)/N/N
以下略
$ nfort -O4 -fopenmp -report-format vectoradd.f90
$ cat vectoradd.L
中略
40: W======> print *,time_end_s - time_begin_s,"sec",sum(c)/N/N
以下略
自動並列化
OpenMPによる並列化ではなく，nfortの自動並列化を利用し，実行時間を測定しました．
$ nfort -mparallel -report-format vectoradd.f90
スレッド数を環境変数VE_OMP_NUM_THREADSで変更して実行しました．結果はOpenMPを使った場合と同じ傾向を示しました．
| スレッド数 | 実行時間 [ms] |
|---|---|
| 1 | 4.66 |
| 2 | 2.70 |
| 4 | 2.23 |
| 8 | 2.27 |
ベクトル和の書き方による違い
ここまでは，同じ書き方でコンパイルオプションを変えていましたが，この節ではベクトル和の書き方を変えてみます．
do concurrent
do concurrentは，Fortran 2008で追加された構文です．do concurrentを用いて繰り返しを記述すると，構文内の処理に依存性がなく，並列に実行できることが明示されるので，コンパイラがdo concurrentを自動で並列化することが期待されます．
OpenMPのディレクティブを全て削除し，doをdo concurrentに書き換えました．
call cpu_time(time_begin_s)
do concurrent(j=1:N, i=1:N)
c(i, j) = a(i, j) + b(i, j)
end do
call cpu_time(time_end_s)
自動並列化オプションを利用してコンパイルしました．
$ nfort -mparallel -report-format vectoradd.f90
診断結果を確認してはみましたが，初めて見る表記になっており，何が行われているのかよくわかりません．Yはループが並列化・ベクトル化されたことを示していますが，+はベクトル化されなかったことを示しているので，つまり・・・なにが行われたのでしょうか･･･？
24: Y------> do concurrent(j=1:N, i=1:N)
25: || c(i, j) = a(i, j) + b(i, j)
26: |+----- end do
環境変数VE_OMP_NUM_THREADSでスレッド数を変更して実行すると，do文を用いているときと同じ結果が得られます．診断結果はよくわかりませんでしたが，無事に並列化されたようです．
| スレッド数 | 実行時間 [ms] |
|---|---|
| 1 | 4.67 |
| 2 | 2.69 |
| 4 | 2.23 |
| 8 | 2.28 |
配列式
do文を用いず，単純に配列式としてベクトル和を記述できます．
c(:, :) = a(:, :) + b(:, :)
また，配列式を記述する際，全範囲を表すコロン(:, :)を記述しないようにもできます．
c = a + b
どちらの書き方をしても，全要素に対して要素同士の和を計算してくれます．配列がallocatable属性を持つ場合には，後者（範囲を記述を記述しない）の書き方をすると，左辺の配列が右辺の配列サイズに合わせて自動的に再割付けされます．
コンパイルは，これまでと同じように自動的並列化オプションを付与して行いました．
$ nfort -mparallel -report-format vectoradd.f90
どちらの書き方をしても，並列化＋ベクトル化が行われます．
24: Y======> c(:, :) = a(:, :) + b(:, :)
24: Y======> c = a + b
ただし，実行結果は大きく異なります．範囲を記述したときの実行結果は，これまでの結果と同じです．しかし，範囲を記述しないと，倍ほど時間がかかっています．自動再割付けの影響なのでしょうか？配列サイズが同じ場合は再割付けされないはずですが，自動再割付けが行われる可能性があることで最適化が阻害されているのかもしれません．-O4オプションを付けても，時間は短縮されませんでした．コンパイラのマニュアルを読んでも，この状況に関する記述はありませんでした．
| スレッド数 | 実行時間 [ms]c(:,:)=a(:,:)+b(:,:)
|
実行時間 [ms]c=a+b
|
実行時間 [ms]c=a+b with -O4
|
|---|---|---|---|
| 1 | 4.66 | 8.26 | 8.84 |
| 2 | 2.72 | 4.66 | 4.94 |
| 4 | 2.23 | 3.79 | 3.80 |
| 8 | 2.28 | 3.84 | 3.89 |
配列式の書き方は非常に長い間議論の対象になっているようで，範囲を書かない方が優勢だったようです．かつては，(:)をつけると配列が一度テンポラリにコピーされるため遅くなったとか，コンパイラの最適化の問題があったとか，色々と教えていただきました．
https://qiita.com/implicit_none/items/9164c4d35d84e6e77f29#comment-110a5430fd8a9ebf4534
https://qiita.com/implicit_none/items/9164c4d35d84e6e77f29#comment-ed185dea986705d2b43c
しかし，NEC Fortranでは範囲を書かない方が有意に高速であるため，自動再割付けを利用する意図がない場合は範囲を書くべきでしょう．
SXユーザの方から，NEC Fortranで範囲を書かないと遅くなるのは昔からだよと教えていただきました．
まとめ
SX-Aurora TSUBASA (SXAT)で動作するプログラムの特徴を調べるために，色々な書き方でベクトル和を計算してみました．
- OpenMPを使っても，コンパイラの自動並列化を使っても，実行速度には影響しませんでした．
-
do concurrentを使ってみたところ，きちんと並列化されました． - 配列式を書く場合は，添字の範囲を表すコロン
(:, :)を省略すると実行時間が倍ほどかかるようになりました．