なんとなーく CPU をぶん回してみようかなとしょうもないネタを。
この辺を少し書き直してみました。ほぼ、マニュアル https://t.co/SpnKeEDZqw を眺めながら取り敢えず動くやつ。
— 1円切手 (&o53) (@jp_yen) 2018年4月10日
4コアやから、マルチスレッドにせんと意味が無いねぇ。表示ルーチンも毎回モジュロ取るよりタイマー割り込みの方が軽いやろうし。F90 ぐらいの固定形式？ pic.twitter.com/6UiJLt6U40
MonteCarlo
! 楕円の面積をモンテカルロ法で求める
! 原点を中心とする横長の楕円を考える、第一象限の面積を求め 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) の正方形の中にランダムに点を打ち、楕円の中に
! 入っているか判定する
!
! 初期化
PROGRAM EllipseArea_MonteCarlo
IMPLICIT NONE
INTEGER I, Istep, Imax, J, K
INTEGER SeedSize
INTEGER,ALLOCATABLE,DIMENSION(:) :: Seed
REAL*16 F, Ra, Rb, L, Lab, Px, Py, RaRb4, AT_1, CORRECT
INTEGER Time_Start, Time_End, t_rate, t_max
REAL Time_Diff
REAL CpuTime_Start, CpuTime_End
! 実行時間を測る
CALL SYSTEM_CLOCK(Time_Start, t_rate, t_max)! 開始時刻
CALL CPU_TIME(CpuTime_Start)
! 変数初期化
Ra = 160
Rb = 90
L = 2 * Ra
J = 0
F = SQRT( Ra ** 2 - Rb ** 2 ) ! 焦点を求める
RaRb4 = Ra * Rb * 4.0 ! 4 * Ra * Rb 計算省略用
AT_1 = ATAN(1.0) ! arctan(1) 計算省略用
CORRECT = RaRb4 * AT_1 ! 面積の真値 計算省略用
Imax = HUGE(Imax) ! Int の最大値まで頑張る
Imax = 50000000 ! いや、頑張れませんでした。
Istep= 1000000 ! 途中経過を表示する単位
Imax = Istep * INT(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(:)) ! 種を与える
! まずは求める真値を表示
WRITE (*, '("Correct ",F13.7)') CORRECT
! モンテカルロ法開始
I = 0
K = 0
DO WHILE(K .LT. Imax)
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
! 計算経過の表示
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
ENDDO
! 計算終了
! 真値 (再掲)
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 EllipseArea_MonteCarlo_Solo.F90 && ./a.exe
SeedSize = 33
Correct 45238.9354706
Aeria 45227.4624000 diff: 0.0254% (Total: 1000000 Inside: 78.520% / 2.00% Done / esp 1.11 sec)
Aeria 45235.6704000 diff: 0.0072% (Total: 2000000 Inside: 78.534% / 4.00% Done / esp 2.20 sec)
Aeria 45239.8272000 diff:-0.0020% (Total: 3000000 Inside: 78.541% / 6.00% Done / esp 3.30 sec)
...snip...
Aeria 45242.3640000 diff:-0.0076% (Total: 48000000 Inside: 78.546% / 96.00% Done / esp 52.45 sec)
Aeria 45242.4125388 diff:-0.0077% (Total: 49000000 Inside: 78.546% / 98.00% Done / esp 53.55 sec)
Aeria 45242.4142080 diff:-0.0077% (Total: 50000000 Inside: 78.546% / 100.00% Done / esp 54.64 sec)
Correct 45238.9354706
elapsed time 54.641sec (←時計の経過時間)
elapsed cpu time 54.547sec (← CPUを使っていた時間)
5千万回で 0.0077% って結構良いですねぇ。
まぁ、でもシングルスレッドだからねぇ。OpenMP にしたらもっと早いのかなぁ？