モンテカルロ法で楕円の面積を求めてみる (OpenMP 化)

# マルチスレッド化する

[前回](https://qiita.com/jp_yen/items/5d561ada97dd4fc6d558 "モンテカルロ法で楕円の面積を求めてみる")、モンテカルロ法で楕円の面積を求めたが 1 Core しか使ってなくて、マルチスレッドにしたら早くなるんじゃないかと OpenMP でのマルチスレッドを試してみた。
4Core CPU (Intel i5-4690) で、マルチスレッドは効果があるのか？

## プログラム

```fortran:EllipseArea_MonteCarlo_Multi.F90
! 楕円の面積をモンテカルロ法で求める

! 原点を中心とする横長の楕円を考える、第一象限の面積を求め 4倍する
! 長半径, 短半径 Ra, Rb
! 焦点の座標 (F,0), (-F,0)
!
! Ra, Rb ; 長半径, 短半径
! F      ; 中心から楕円の焦点までの距離 Fa(-F, 0), Fb(F, 0)
! Px, Py ; ランダムなポイント P(Px, Py)
! I      ; カウンタ
! Imax   ; 最大試行回数
! J      ; 楕円の範囲内であった数
! K      ; 試行回数
! L      ; 焦点 Fa, Fb から円周上点P までの距離の合計 2 * Ra
! Lab    ; 焦点 Fa, Fb からランダムな点P までの距離の合計
!
! (0,0) - (Ra,Rb) の正方形の中にランダムに点を打ち、楕円の中に
! 入っているか判定する
!
! export OMP_NUM_THREADS=スレッド数
!
      PROGRAM EllipseArea_MonteCarlo
      !$ use omp_lib
      IMPLICIT NONE
      INTEGER I, Istep, Imax, J, K
      INTEGER SeedSize
      INTEGER,ALLOCATABLE,DIMENSION(:) :: Seed
      REAL*16 Lab, Px, Py

      ! 定数初期化
      REAL*16,PARAMETER :: Ra = 160
      REAL*16,PARAMETER :: Rb =  90
      REAL*16,PARAMETER :: L = 2 * Ra
      REAL*16,PARAMETER :: F = SQRT( Ra ** 2 - Rb ** 2 )     ! 焦点を求める
      REAL*16,PARAMETER :: RaRb4   = Ra * Rb * 4.0           ! 4 * Ra * Rb 計算省略用
      REAL*16,PARAMETER :: AT_1    = ATAN(1.0)               ! arctan(1)   計算省略用
      REAL*16,PARAMETER :: CORRECT = RaRb4 * AT_1            ! 面積の真値  計算省略用
      ! 時間計測用
      INTEGER Time_Start, Time_End, t_rate, t_max
      REAL Time_Diff
      REAL CpuTime_Start, CpuTime_End

      ! 変数初期化
      J = 0
      K = 0

      ! 開始時刻を記録
      CALL SYSTEM_CLOCK(Time_Start, t_rate, t_max)! 開始時刻
      CALL CPU_TIME(CpuTime_Start)

      Imax = 50000000                   ! 最大値までは頑張れない
      Istep= 1000000                    ! 途中経過を表示する単位
      Imax = Imax - MOD(Imax,Istep)     ! Imax を突き抜けないように計算しなおす

      ! 乱数初期化
      !
      CALL RANDOM_SEED(SIZE = SeedSize) ! seed のサイズを求める
      WRITE (*,*) 'SeedSize = ', SeedSize
      ALLOCATE(Seed(SeedSize))          ! 配列の割り当て
      DO I = 1, SeedSize
        CALL SYSTEM_CLOCK(count = Seed(i))! 現在時間を乱数の種に
      END DO
      CALL RANDOM_SEED(put = Seed(:))   ! 種を与える
      DEALLOCATE(Seed)

      ! まずは求める真値を表示
      WRITE (*, '("Correct ",F13.7)') CORRECT
      ! モンテカルロ法開始

!$omp parallel private(I, Px, Py, Lab)
!$omp single
      print *, "threads no ", omp_get_num_threads()
!$omp end single
      print *, "IN SINGLE: My Thread Number is ", omp_get_thread_num()

      DO WHILE(K .LT. Imax)
!$omp do reduction(+:J, K)
!ループを出た後に必要となるのは J, K のみ。
!こいつらは、ループを出た後に足し合わせる必要あり。
        DO I = 1, Istep
          ! (0,0) - (Ra, Rb) の中のランダムな点P
          CALL RANDOM_NUMBER(Px)
          CALL RANDOM_NUMBER(Py)
          Px = Px * Ra
          Py = Py * Rb

          ! 各焦点から、点P までの距離の合計
          Lab = SQRT((-F + Px) ** 2 + Py ** 2) + &
                SQRT(( F + Px) ** 2 + Py ** 2)
          IF (Lab .LE. L) THEN
              J = J + 1         ! 楕円の中ならカウントアップ
          END IF
          K = K + 1             ! 試行回数をカウント
        END DO
!$omp end do
!$omp single
!経過表示をマルチスレッドで行うとおかしくなるので、
!ここで待ち合わせをして 1スレッドで表示する。
        ! 計算経過の表示
        CALL SYSTEM_CLOCK(Time_End)
        WRITE (*,100) RaRb4 * J / K, 100 * (1 - J/(K * AT_1)), K, &
                100.0*J/K, 100.0*K/Imax, DBLE(Time_End - Time_Start)/t_rate
!$omp end single
      END DO
!$omp end parallel

! 計算終了
! 真値 (再掲)
      WRITE (*, '("Correct ",F13.7)') CORRECT

      ! 終了時刻 実時間 & CPU 時間
      CALL SYSTEM_CLOCK(Time_End)
      CALL CPU_TIME(CpuTime_End)
      IF ( Time_End < Time_Start ) THEN
        Time_Diff = (t_max - Time_Start) + Time_End + 1
      ELSE
        Time_Diff = Time_End - Time_Start
      ENDIF

      Time_Diff = DBLE(Time_Diff/t_rate)        ! 単位を秒へ変換
      WRITE (*,'("elapsed time ", F0.3, "sec")') Time_Diff
      WRITE (*,'("elapsed cpu time ", F0.3, "sec")') &
        (CpuTime_End - CpuTime_Start)

! 面積 F13.7 (誤差、試行回数 I10、楕円の内側の割合、進行割合 F6.2%
100   FORMAT('  Aeria ', F13.7, ' diff:', F7.4, '% (Total:',I10 , &
        ' Inside:', F7.3,'% / ', F6.2'% Done / esp ', F0.2, ' sec)')

      END
```
## 実行！

```shell:tcsh
% gfortran -fopenmp EllipseArea_MonteCarlo_Multi.F90 && ./a.exe
 SeedSize =           33
Correct 45238.9354706
 threads no            4
 IN SINGLE: My Thread Number is            2
 IN SINGLE: My Thread Number is            0
 IN SINGLE: My Thread Number is            3
 IN SINGLE: My Thread Number is            1
  Aeria 45251.5968000 diff:-0.0280% (Total:   1000000 Inside: 78.562% /   2.00% Done / esp .31 sec)
  Aeria 45247.7376000 diff:-0.0195% (Total:   2000000 Inside: 78.555% /   4.00% Done / esp .61 sec)
  Aeria 45231.8016000 diff: 0.0158% (Total:   3000000 Inside: 78.527% /   6.00% Done / esp .92 sec)
...snip...
  Aeria 45234.7308000 diff: 0.0093% (Total:  48000000 Inside: 78.533% /  96.00% Done / esp 14.00 sec)
  Aeria 45234.2027755 diff: 0.0105% (Total:  49000000 Inside: 78.532% /  98.00% Done / esp 14.28 sec)
  Aeria 45234.5114880 diff: 0.0098% (Total:  50000000 Inside: 78.532% / 100.00% Done / esp 14.58 sec)
Correct 45238.9354706
elapsed time 14.578sec (← 経過時間)
elapsed cpu time 55.531sec (← CPUをぶん回した時間延べ時間 4スレッド分)
```

おぉ、すごい。なら、強制的に 8 スレッドにしてみたら？`$omp parallel` に `num_threads(n)` を付けて...

```
!$omp parallel private(I, Px, Py, Lab), num_threads(8)
```

```shell:tcsh
% setenv OMP_NUM_THREADS 8
% ./a.exe
 SeedSize =           33
Correct 45238.9354706
 threads no            8
 IN SINGLE: My Thread Number is            1
 IN SINGLE: My Thread Number is            3
 IN SINGLE: My Thread Number is            2
 IN SINGLE: My Thread Number is            5
 IN SINGLE: My Thread Number is            6
 IN SINGLE: My Thread Number is            4
 IN SINGLE: My Thread Number is            0
 IN SINGLE: My Thread Number is            7
  Aeria 45265.9968000 diff:-0.0598% (Total:   1000000 Inside: 78.587% /   2.00% Done / esp .30 sec)
  Aeria 45237.3408000 diff: 0.0035% (Total:   2000000 Inside: 78.537% /   4.00% Done / esp .58 sec)
  Aeria 45241.4016000 diff:-0.0055% (Total:   3000000 Inside: 78.544% /   6.00% Done / esp .88 sec)
...snip...
  Aeria 45240.3141120 diff:-0.0030% (Total:  50000000 Inside: 78.542% / 100.00% Done / esp 15.05 sec)
Correct 45238.9354706
elapsed time 15.047sec (← 経過時間)
elapsed cpu time 55.500sec (← CPUをぶん回した時間延べ時間 8スレッド分)
```

むりか...まぁ、4スレッド の時点で CPU を 95% 以上使ってたからねぇ。
そして、Ryzen のノート PC (HUAWEI MateBook D 15) を入手したのでこちらもテスト。

## タイムアタック
それぞれの環境での実行時間。
[google スプレッドシート での結果](https://docs.google.com/spreadsheets/d/10He09b08ZezuzDAIxgUyKm4UkkQaj_bQlCIkZPqS5Yk/edit?usp=sharing)

|CPU|シングルスレッド|2スレッド<br>(1スレッド比)| 4スレッド<br>(1スレッド比) | 8スレッド<br>(1スレッド比)|16スレッド|
|---|--:|--:|--:|--:|--:|
|Core 2 Duo T9500<br>実時間|101.060 秒|55.229 秒<br>(-45%)|-|-||
|Core 2 Duo T9500<br>CPU時間|100.891 秒|99.844 秒<br>(-1%)|-|-||
|i5-4690<br>実時間|54.641 秒|-|14.578 秒<br>(-73%)|15.047 秒<br>(-72%)|-|
|i5-4690<br>CPU時間|54.547 秒|-|55.531 秒<br>(+2%)|55.500 秒<br>(+2%)|-|
|E5-2623 v4×2<br>実時間|68.601 秒|35.662 秒<br>(-48%)|23.974 秒<br>(-65%)|18.080 秒<br>(-74%)|7.872 秒<br>(-89%)|
|E5-2623 v4×2<br>CPU時間|68.600 秒|70.177 秒<br>(+2%)|81.360 秒<br>(+19%)|109.550 秒<br>(+60%)|123.634 秒<br>(+80%)|
|i5-8200Y<br>実時間|70.436 秒|-|35.039 秒|-|-|
|i5-8200Y<br>CPU時間|70.436 秒|-|139.482 秒|-|-|
|Ryzen 5 3500U<br>実時間|50.747 秒|-|16.242 秒<br>(-68%)|10.018 秒<br>(-80%)|-|
|Ryzen 5 3500U<br>CPU時間|50.688 秒|-|63.422 秒<br>(+25%)|78.594 秒<br>(+55%)||
|i5-10210U<br>実時間|47.352 秒|24.308 秒<br>(-49%)|15.365 秒<br>(-68%)|12.498 秒<br>(-74%)|-|
|i5-10210U<br>CPU時間|47.312 秒|48.453 秒<br>(+2%)|59.500 秒<br>(+26%)|96.500 秒<br>(+104%)|-|

CPUの諸元

|CPU|発売開始|コア数|キャッシュ|バススピード|クロック|TDP|プロセスルール|内蔵グラフィック|
|:--|---|---|---|---|---|---|---|---|
|Intel Core 2 Duo T9500|2009年1月|2コア|6MB L2 Cache|800MHz|2.60GHz<br>(ターボブーストなし)|35W|45nm (Penryn)|なし|
|Intel Core i5-4690|2014年5月11日|4コア|6MB Intel® Smart Cache|4 GT/s|ベース 3.5GHz<br>ブースト 3.9GHz|84W|22nm (Haswell Refresh)<br>Core i シリーズ 第4世代|インテル® HD グラフィックス 4600|
|Intel® Xeon® E5-2623 v4 x 2個|2016 6月20日|4C 8T|10 MB Intel® Smart Cache|8 GT/s|ベース 2.60 GHz<br>ブースト 3.20 GHz|85W|14nm (Broadwell)|－| 
|Core™ i5-8200Y|2018年08月|2C 4T|L1:128 KB, L2:512 KB, L3:4MB|5GT/s|ベース 1.30 GHz<br>ブースト 3.90GHz|5W|14 nm (Amber Lake Y)|インテル® UHD グラフィックス 615|
|AMD Ryzen 5 3500U|2019年1月7日|4C 8T|L1:384KB, L2:2MB, L3:4MB|?|ベース 2.1GHz<br>ブースト 3.7GHz|15W|12nm (Picasso)<br>Ryzen Mobile APU 第2世代|Radeon™ Vega 8|
|Intel® Core™ i5-10210U|2019年8月|4C 8T|L1:64K/Core, L2:256K/Core, L3 6MB/Total|4 GT/s|ベース 1.60 GHz<br>ブースト 4.20 GHz|15W|14nm (Comet Lake)|インテル® Core™ UHD グラフィックス|

[Intel® Core™ 2 Duo T9500](https://ark.intel.com/content/www/jp/ja/ark/products/33918/intel-core-2-duo-processor-t9500-6m-cache-2-60-ghz-800-mhz-fsb.html) : Lenovo G550 2958GCJ (CPU換装, Mem 8GB, SSD 1TB)
[Intel® Core™ i5-4690](https://ark.intel.com/content/www/jp/ja/ark/products/80810/intel-core-i5-4690-processor-6m-cache-up-to-3-90-ghz.html) : ASUS H97-Pro (Mem 32GB, SSD CT1000MX500SSD1/JP, HDD...)
[インテル® Core™ i5-8200Y](https://ark.intel.com/content/www/jp/ja/ark/products/185280/intel-core-i5-8200y-processor-4m-cache-up-to-3-90-ghz.html) : NEC VersaPro VKT13/H-5 (Mem 8GB, SSD 256GB)
[Intel® Xeon® E5-2623 v4](https://ark.intel.com/content/www/jp/ja/ark/products/92980/intel-xeon-processor-e52623-v4-10m-cache-2-60-ghz.html) : PowerEdge DL380 G9 (MEM 128GB)
[AMD Ryzen™ 5 3500U](https://www.amd.com/ja/products/apu/amd-ryzen-5-3500u) : HUAWEI MateBook D 15 BOH-WAQ9RP (Mem 16GB, SSD 256GB, HDD 1TB)
[Intel® Core™ i5-10210U](https://ark.intel.com/content/www/jp/ja/ark/products/185280/intel-core-i5-8200y-processor-4m-cache-up-to-3-90-ghz.html) : Panasonic CF-SV9HDLVS (Mem 8GB, SSD 256GB)

## 環境
それぞれ、Windows の WSL 上の Fedora Remix for WSL です。 

```shell:Core_2_Doo_T9500
% uname -a
Linux Yudai 4.4.0-18362-Microsoft #476-Microsoft Fri Nov 01 16:53:00 PST 2019 x86_64 GNU/Linux
% gfortran --version
GNU Fortran (Debian 6.3.0-18+deb9u1) 6.3.0 20170516
Copyright (C) 2016 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.
```

```shell:i5-8200Y
% uname -a
Linux WIN 4.19.128-microsoft-standard #1 SMP Tue Jun 23 12:58:10 UTC 2020 x86_64 x86_64 x86_64 GNU/Linux
% cat /etc/redhat-release
Generic release 32 (Generic)
% gfortran --version
GNU Fortran (GCC) 10.2.1 20201016 (Red Hat 10.2.1-6)
Copyright (C) 2020 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.
```

```shell:Ryzen_5_3500U
% uname -a
Linux HINAKO 4.4.0-19041-Microsoft #1-Microsoft Fri Dec 06 14:06:00 PST 2019 x86_64 GNU/Linux
% gfortran --version
GNU Fortran (Debian 8.3.0-6) 8.3.0
Copyright (C) 2018 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.
```

```shell:i5-10210U
% uname -a
Linux Win 4.4.0-19041-Microsoft #1202-Microsoft Wed Aug 25 18:06:00 PST 2021 x86_64 x86_64 x86_64 GNU/Linux
% gfortran --version
GNU Fortran (GCC) 11.2.1 20210728 (Red Hat 11.2.1-1)
Copyright (C) 2021 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.

% cat /etc/os-release
NAME=Generic
VERSION="34 (Generic)"
ID=generic
ID_LIKE=fedora
VERSION_ID=34
PRETTY_NAME="Generic 34 (Generic)"
ANSI_COLOR="0;34"
LOGO=generic-logo-icon
CPE_NAME="cpe:/o:generic:generic:34"
HOME_URL="http://www.zombo.com/"
SUPPORT_URL="https://en.wikipedia.org/wiki/Help!_(album)"
BUG_REPORT_URL="https://youtu.be/CSemARaqGqE"
REDHAT_BUGZILLA_PRODUCT="Generic"
REDHAT_BUGZILLA_PRODUCT_VERSION=%{bug_version}
REDHAT_SUPPORT_PRODUCT="Generic"
REDHAT_SUPPORT_PRODUCT_VERSION=%{bug_version}
PRIVACY_POLICY_URL="http://nsa.gov"
```

```shell:Intel® Xeon® E5-2623 v4 x 2個
% uname -a
FreeBSD FreeBSD 13.2-RC1 FreeBSD 13.2-RC1 releng/13.2-n254563-13264ea9a370 GENERIC amd64
% gfortran12 --version
GNU Fortran (FreeBSD Ports Collection) 12.2.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.


```

当然だけど Intel 4Core CPU に対して 4スレッドはほぼ 4倍になりましたねぇ。そして、8スレッドにしても効率はあまり落ちないんですね。これは、16Core とか 32Core とかでも試してみたいな。ネカフェでオンラインゲーム用席を借りたら試せるかな？
Core 2 Duo はシングルスレッドで大体 2倍の時間がかかる。全体で 4～5倍。10年というのはそれぐらいと。(ハイエンドと、ミドルレンジの比較ですが)
Ryzen 5 電力を考えるとワッパ圧倒的ですねぇ。8コアで、5倍ってマルチスレッドの効率...いや、4コアか。4コアで5倍なら凄いな。