Fortran+OpenMPによる自動SIMD化で高速化できるらしい

:::note info
この記事は以下のサイトの受け売りで, `Hello, World!` 的なソースコードで検証をしようとしています.
詳しく知りたい場合は以下のサイトを読む方が良いでしょう.
:::

https://www.cc.u-tokyo.ac.jp/kashiwa_open/files/25simd.pdf

https://www.intel.com/content/www/us/en/developer/articles/technical/explicit-vector-programming-in-fortran.html


# やりたいこと
プログラムをもっと高速化したい!
並列化程度では満足できない!

# SIMDってなに?
Single Instruction Multiple Data の略です. 一回の操作で複数のデータを処理できます.
128ビットや256ビットのデータを一括で処理できるらしく, 例えば, 32ビット整数の場合は4つや8つを一気に計算できるため, 理想的には4, 8倍の速度が出ます.
夢のような話ですが, 手動で SIMD 化するには, [intel のサイト](https://www.intel.com/content/www/us/en/docs/intrinsics-guide/index.html) を見て命令を調べたり, `__mm256i` みたいな型をアレコレしたり, アドレス配置を意識したりする必要があります.
また, Fortran では手動で SIMD 化できません.
コンパイラが自動で SIMD 化してくれればよいのに...ということで, OpenMP の **SIMDディレクティブ** の出番です.

# Fortran
## プログラムの速度
プログラミング言語 Fortran は主に科学技術計算で(現在も)用いられている言語です.
初出は1950年代なので, 古いと思われがちですが, 今もなおアップデートされており, 1977, 1990, 1995, 2003, 2008, 2018, 2023年に規格が新しくなっています.
速さはC, C++, Rust 並に速いです. **速いのです**.
速い上に, ポインタみたいな低レベルな場所も触らなくて良いので **簡単で速い** です.
しかし, Fortran は低レベルな場所はあんまり触れないため, SIMD での高速化は **(コンパイラ)神頼み** になると思います.
やはり, 数値計算は速さがあればあるほど, たくさんデータが取れたり, スーパーコンピュータの使用ポイントを削減したり, シミュレーションを行える数が増えたりするので, **速くしたい** というのが人の情というものでしょう.
しかし, コンパイラが勝手に SIMD 化できないっぽいです. おそらく, 勝手にやるのは(データの独立性などの問題で)難しいのでしょう.
そこで登場するのが, OpenMP のバージョン4.5以降の **自動SIMD化** です.
コード書く人が, 処理を SIMD 化できるというヒントをコンパイラに教えることで高速化させることができます.

## [OpenMP](https://www.openmp.org/)
OpenMP は並列化するやつであり, Fortran の言語仕様に含まれているらしいです.
つまり, どの Fortran コンパイラでも OpenMP に対応しているということでしょう.
OpenMP 4.5 以降から自動SIMD化できるようです.
OpenMP は複数のコアで並列化を試みますが, OpenMP の SIMD 化はシングルコアでも動くと思います.

# 自動SIMD化してみる

## ソースコード
```test_simd.f90
module auto_simd_m
  use, intrinsic :: iso_fortran_env
  implicit none
contains
  subroutine vectorize_add(a, b, c)
    !$omp declare simd linear(ref(a, b, c))
    real(real32), intent(in) :: a, b
    real(real32), intent(inout) :: c
    c = c + a + b
  end subroutine vectorize_add
end module auto_simd_m

program test_simd
  use, intrinsic :: iso_fortran_env
  use auto_simd_m
  implicit none
  integer(int32), parameter :: n = 10 ** 8, nloop = 100
  real(real32), allocatable :: a(:), b(:), c(:)
  real(real32) :: summ
  integer(int32) :: i, j
  allocate(a(n), b(n), c(n))
  a(:) = 1
  b(:) = 2
  c(:) = 0
  do j = 1, nloop
      !$omp simd
      do i = 1, n
        call vectorize_add(a(i), b(i), c(i))
      end do
  end do
  write(*, *) c(n / 2), c(1) * n
  summ = 0
  !$omp simd reduction(+:summ)
  do i = 1, n
      summ = summ + c(i)
  end do
  write(*, *) summ
end program test_simd
```

`Hello, World` な感じなソースコードを書いてみました.
サブルーチンにして, `!$omp declare simd` みたいな感じにすれば良いらしいです.

## 性能評価
### 環境
```bash
$ lsb_release -a
No LSB modules are available.
Distributor ID:	Ubuntu
Description:	Ubuntu 22.04.5 LTS
Release:	22.04
Codename:	jammy
$ gfortran --version
GNU Fortran (Ubuntu 12.3.0-1ubuntu1~22.04) 12.3.0
Copyright (C) 2022 Free Software Foundation, Inc.
This is free software; see the source for copying conditions.  There is NO
warranty; not even for MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
$ ifx --version
ifx (IFX) 2025.0.1 20241113
Copyright (C) 1985-2024 Intel Corporation. All rights reserved.
$ grep -e 'avx' /proc/cpuinfo
(略) avx (略)
```
このCPUにはAVX2, AVX512は含まれていません.

実行結果
```bash
$ export OMP_NUM_THREADS=1
$ gfortran -O -fno-openmp-simd test_simd.f90 -fopt-info && time ./a.out
   300.000000       3.00000010E+10
   8.58993459E+09

real	0m13.577s
user	0m13.138s
sys	0m0.428s
$ gfortran -O -fopenmp-simd test_simd.f90 -fopt-info && time ./a.out
test_simd.f90:35:24: optimized: loop vectorized using 16 byte vectors
test_simd.f90:28:44: optimized: loop vectorized using 16 byte vectors
test_simd.f90:13:17: optimized: loop with 3 iterations completely unrolled (header execution count 1206603)
test_simd.f90:13:17: optimized: loop with 3 iterations completely unrolled (header execution count 1355734)
   300.000000       3.00000010E+10
   3.43597384E+10

real	0m11.568s
user	0m11.118s
sys	0m0.441s
```

一応速くなっていますが, 倍にもなっていないです.
gfortran なら `-fopt-info` で高速化の情報が, `-fopenmp-simd` でOpenMPのSIMDが有効になります.
また, `-O` オプションで手加減していますが, 処理が単純なので, `-O3` オプションでベクトル化(SIMD化)されます.
```bash
$ gfortran -O3 test_simd.f90 -fopt-info && time ./a.out
test_simd.f90:28:44: optimized:  Inlining vectorize_add/0 into test_simd/1.
test_simd.f90:25:8: optimized: applying unroll and jam with factor 2
test_simd.f90:24:10: optimized: Loop 3 distributed: split to 0 loops and 1 library calls.
test_simd.f90:34:8: optimized: loop vectorized using 16 byte vectors
test_simd.f90:27:12: optimized: loop vectorized using 16 byte vectors
test_simd.f90:23:10: optimized: loop vectorized using 16 byte vectors
test_simd.f90:22:10: optimized: loop vectorized using 16 byte vectors
test_simd.f90:37:18: optimized: basic block part vectorized using 16 byte vectors
test_simd.f90:37:18: optimized: basic block part vectorized using 16 byte vectors
   300.000000       3.00000010E+10
   8.58993459E+09

real	0m6.251s
user	0m5.684s
sys	0m0.562s
```

- 単純な場合は, コンパイラに任せた方が良いのかもしれません.
- `real(real32)` では精度が足りないので `3E+10` に一致しません.
- 和を `$omp reduction (+:summ)` として, `-fopenmp-simd` した場合は精度が多少良いです.
  + 足し算の順番が変わって誤差が累積しづらくなるためだと思います.

```bash
$ ifx -O3 -qno-openmp-simd -qopt-report -qopt-report-phase=all -qopt-report-file=no_openmpsimd.txt test_simd.f90 && time ./a.out
   300.0000      3.0000001E+10
  2.8862222E+10

real	0m23.602s
user	0m23.123s
sys	0m0.474s
$ ifx -O3 -qopenmp-simd -qopt-report -qopt-report-phase=all -qopt-report-file=openmpsimd.txt test_simd.f90 && time ./a.out
   300.0000      3.0000001E+10
  2.8862222E+10

real	0m11.290s
user	0m10.821s
sys	0m0.444s
```

- 2倍程度速くなっています.
- このソースコードの場合, intel fortran compiler の方が gfortran よりも遅いです. (不思議?)
- OpenMP 6.0 から `!$omp declare_simd(...)` になったっ[ぽい](https://www.openmp.org/wp-content/uploads/OpenMP-RefGuide-6.0-OMP60SC24-web.pdf) です.
- `-qopenmp-simd` で OpenMP の SIMD が有効になります.
- `-qopt-report -qopt-report-phase=all -qopt-report-file=(filename)` で　`(filename)` に高速化の情報が書き込まれます.
  + (filename) には `stdout` や `stderr` も指定できます.

```diff
$ diff -u no_openmpsimd.txt openmpsimd.txt 
--- no_openmpsimd.txt	2024-12-10 23:54:06.548175633 +0900
+++ openmpsimd.txt	2024-12-10 23:54:13.034202222 +0900
@@ -42,18 +42,13 @@
     remark #15541: loop was not vectorized: outer loop is not an auto-vectorization candidate.
 
     LOOP BEGIN at test_simd.f90 (27, 7)
-        remark #15344: Loop was not vectorized: vector dependence prevents vectorization
-        remark #15346: vector dependence: assumed FLOW dependence between c [ /tmp/test_simd.f90 (9, 5) ] and a [ /tmp/test_simd.f90 (28, 9) ] 
-        remark #15346: vector dependence: assumed FLOW dependence between c [ /tmp/test_simd.f90 (9, 5) ] and a [ /tmp/test_simd.f90 (28, 9) ] 
-        remark #15346: vector dependence: assumed FLOW dependence between c [ /tmp/test_simd.f90 (9, 5) ] and b [ /tmp/test_simd.f90 (28, 9) ] 
-        remark #15346: vector dependence: assumed FLOW dependence between c [ /tmp/test_simd.f90 (9, 5) ] and b [ /tmp/test_simd.f90 (28, 9) ] 
-        remark #15346: vector dependence: assumed FLOW dependence between c [ /tmp/test_simd.f90 (9, 5) ] and c [ /tmp/test_simd.f90 (28, 9) ] 
-        remark #15346: vector dependence: assumed FLOW dependence between c [ /tmp/test_simd.f90 (9, 5) ] and c [ /tmp/test_simd.f90 (28, 9) ] 
+        remark #15301: SIMD LOOP WAS VECTORIZED
+        remark #15305: vectorization support: vector length 16
     LOOP END
 LOOP END
 
 LOOP BEGIN at test_simd.f90 (34, 3)
-    remark #15300: LOOP WAS VECTORIZED
+    remark #15301: SIMD LOOP WAS VECTORIZED
     remark #15305: vectorization support: vector length 4
 LOOP END
 =================================================================
@@ -111,9 +106,9 @@
 
 Hardware registers used
  Reserved : 4 [ mxcsr rsp ssp rip  ]
- Available : 23 [ rbp r9 r10 r11 r12 r13 r14 r15 zmm5 zmm6 zmm7 zmm8 zmm9 zmm10 zmm11 zmm12 zmm13 zmm14 zmm15 k0_k1 k2_k3 k4_k5 k6_k7  ]
+ Available : 18 [ rbp r11 r12 r13 r14 r15 zmm8 zmm9 zmm10 zmm11 zmm12 zmm13 zmm14 zmm15 k0_k1 k2_k3 k4_k5 k6_k7  ]
  Callee Saved : 1 [ rbx  ]
- Assigned : 13 [ eflags rax rbx rcx rdi rdx rsi r8 zmm0 zmm1 zmm2 zmm3 zmm4  ]
+ Assigned : 18 [ eflags rax rbx rcx rdi rdx rsi r8 r9 r10 zmm0 zmm1 zmm2 zmm3 zmm4 zmm5 zmm6 zmm7  ]
 
 Routine stack 
  Variables : 104 Bytes
```

- OpenMP SIMD でベクトル化されると `remark #15301: SIMD LOOP WAS VECTORIZED` となるっぽいです.
- ベクトル化できないと `remark #15344: Loop was not vectorized: (理由)` となるようです.
- 単純なループを `-O3` でベクトル化できないのは意外です. (subroutine の書き方が悪いのかもしれません.)


# まとめ
`Hello, World!` 感覚なソースコードを高速化しようと試みました.
しかし, あんまり高速化の実感がありません.
おそらく, 
- AVX2ではない
- コードの規模が小さい
- 計算が重くない
などが原因ではないかと思います.
次は実践的なコードで試してみたいと思いました.

# 参考
https://www.cc.u-tokyo.ac.jp/kashiwa_open/files/25simd.pdf

https://www.intel.com/content/www/us/en/developer/articles/technical/explicit-vector-programming-in-fortran.html

## SIMD
https://www.intel.com/content/www/us/en/docs/intrinsics-guide/index.html

## OpenMP
https://www.openmp.org/

https://www.openmp.org/wp-content/uploads/OpenMP-RefGuide-6.0-OMP60SC24-web.pdf