LoginSignup
1
1
@jp_yen(切手 一円)

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

Last updated at Posted at 2018-11-07

マルチスレッド化する

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

プログラム

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

実行！

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)
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 スプレッドシート での結果

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

CPUの諸元

CPU 発売開始 コア数 キャッシュ バススピード クロック TDP プロセスルール 内蔵グラフィック
Intel Core 2 Duo T9500 2009年1月 2コア 6MB L2 Cache 800MHz 2.60GHz
(ターボブーストなし)		 35W 45nm (Penryn) なし
Intel Core i5-4690 2014年5月11日 4コア 6MB Intel® Smart Cache 4 GT/s ベース 3.5GHz
ブースト 3.9GHz		 84W 22nm (Haswell Refresh)
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
ブースト 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
ブースト 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
ブースト 3.7GHz		 15W 12nm (Picasso)
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
ブースト 4.20 GHz		 15W 14nm (Comet Lake) インテル® Core™ UHD グラフィックス

Intel® Core™ 2 Duo T9500 : Lenovo G550 2958GCJ (CPU換装, Mem 8GB, SSD 1TB)
Intel® Core™ i5-4690 : ASUS H97-Pro (Mem 32GB, SSD CT1000MX500SSD1/JP, HDD...)
インテル® Core™ i5-8200Y : NEC VersaPro VKT13/H-5 (Mem 8GB, SSD 256GB)
Intel® Xeon® E5-2623 v4 : PowerEdge DL380 G9 (MEM 128GB)
AMD Ryzen™ 5 3500U : HUAWEI MateBook D 15 BOH-WAQ9RP (Mem 16GB, SSD 256GB, HDD 1TB)
Intel® Core™ i5-10210U : Panasonic CF-SV9HDLVS (Mem 8GB, SSD 256GB)

環境

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

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.
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.
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.
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"
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倍なら凄いな。

1
1
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
1
1