stochastic_rk Fortran 90 代码还有进一步优化的余地吗？

Question

我需要使用 Fortran 代码来求解随机微分方程 (SDE)。我看了 Burkardt 著名的 Fortran 代码网站，

https://people.math.sc.edu/Burkardt/f_src/stochastic_rk/stochastic_rk.html

我特别查看了 stochastic_rk.f90 代码中的 rk4_ti_step 子例程，

https://people.math.sc.edu/Burkardt/f_src/stochastic_rk/stochastic_rk.f90

下面是我优化后的版本，

subroutine rk4_ti_step_mod ( x, t, h, q, fi, gi, seed, xstar )
use random  
implicit none
real ( kind = 8 ), external :: fi
real ( kind = 8 ), external :: gi
real ( kind = 8 ) h
real ( kind = 8 ) k1
real ( kind = 8 ) k2
real ( kind = 8 ) k3
real ( kind = 8 ) k4
real ( kind = 8 ) q
real ( kind = 8 ) r8_normal_01
integer ( kind = 4 ) seed
real ( kind = 8 ) t
real ( kind = 8 ) t1
real ( kind = 8 ) t2
real ( kind = 8 ) t3
real ( kind = 8 ) t4
real ( kind = 8 ) w1
real ( kind = 8 ) w2
real ( kind = 8 ) w3
real ( kind = 8 ) w4
real ( kind = 8 ) x
real ( kind = 8 ) x1
real ( kind = 8 ) x2
real ( kind = 8 ) x3
real ( kind = 8 ) x4
real ( kind = 8 ) xstar   
real ( kind = 8 ) :: qoh
real ( kind = 8 ) :: normal(4)
real ( kind = 8 ), parameter :: a21 = 2.71644396264860D+00 &
,a31 = - 6.95653259006152D+00 &
,a32 =   0.78313689457981D+00 &
,a41 =   0.0D+00 &
,a42 =   0.48257353309214D+00 &
,a43 =   0.26171080165848D+00 &
,a51 =   0.47012396888046D+00 &
,a52 =   0.36597075368373D+00 &
,a53 =   0.08906615686702D+00 &
,a54 =   0.07483912056879D+00 &
,q1 =   2.12709852335625D+00 &
,q2 =   2.73245878238737D+00 &
,q3 =  11.22760917474960D+00 &
,q4 =  13.36199560336697D+00
real ( kind = 8 ), parameter, dimension(4) :: qarray = [ 2.12709852335625D+00 &
    ,2.73245878238737D+00 &
    ,11.22760917474960D+00 &
    ,13.36199560336697D+00 ]
real ( kind = 8 ) :: warray(4)
integer (kind = 4) :: i
qoh = q / h
normal = gaussian(4) 
do i =1,4
    warray(i) = normal(i)*sqrt(qarray(i)*qoh)
enddo
t1 = t
x1 = x
k1 = h * ( fi ( x1 ) + gi ( x1 ) * warray(1) ) 
t2 = t1 + a21 * h
x2 = x1 + a21 * k1
k2 = h * ( fi ( x2 ) + gi ( x2 ) * warray(2) )
t3 = t1 + ( a31  + a32 )* h
x3 = x1 + a31 * k1 + a32 * k2
k3 = h * ( fi ( x3 ) + gi ( x3 ) * warray(3) )
t4 = t1 + ( a41  + a42 + a43 ) * h
x4 = x1 + a41 * k1 + a42 * k2
k4 = h * ( fi ( x4 ) + gi ( x4 ) * warray(4) )
xstar = x1 + a51 * k1 + a52 * k2 + a53 * k3 + a54 * k4
return
end

注意我用的是我的随机数模块，gaussian是我的随机数函数，这部分无所谓

我只是想知道，

任何人都可以就代码进一步优化提出一些建议吗？
有谁知道 best/fastest SDE Fortran 子程序是什么？或者什么算法最好？

非常感谢！

Answer 1

x 和 c 的相互依赖意味着你不能像我最初想的那样把太多的东西变成线性代数，但我仍然希望通过将所有东西分组到适当的数组中来加快速度:

subroutine rk4_ti_step_mod ( x, t, h, q, fi, gi, seed, xstar )
  use random  
  implicit none
  
  integer, parameter :: dp = selected_real_kind(15,307)
  integer, parameter :: ip = selected_int_kind(9)
  
  real(dp), intent(in) :: x
  real(dp), intent(in) :: t  
  real(dp), intent(in) :: h
  real(dp), intent(in) :: q
  real(dp), external :: fi
  real(dp), external :: gi
  integer(ip), intent(in) :: seed
  real(dp), intent(out) :: xstar

  real(dp), parameter :: as(4,5) = reshape([ &
     &  0.0_dp,              0.0_dp,              0.0_dp,              0.0_dp, &
     &  2.71644396264860_dp, 0.0_dp,              0.0_dp,              0.0_dp, &
     & -6.95653259006152_dp, 0.78313689457981_dp, 0.0_dp,              0.0_dp, &
     &  0.0_dp,              0.48257353309214_dp, 0.26171080165848_dp, 0.0_dp, &
     &  0.47012396888046_dp, 0.36597075368373_dp, 0.08906615686702_dp, 0.07483912056879_dp &
     & ], [4,5])
  real(dp), parameter :: qs(4) = [ &
     &  2.12709852335625_dp, &
     &  2.73245878238737_dp, &
     & 11.22760917474960_dp, &
     & 13.36199560336697_dp ]

  real(dp) :: ks(4)
  real(dp) :: r8_normal_01
  real(dp) :: ts(4)
  real(dp) :: ws(4)
  real(dp) :: xs(4)
  real(dp) :: normal(4)
  real(dp) :: warray(4)
  
  normal = gaussian(4) 
  warray = normal*sqrt(qs)*sqrt(q/h)
  
  do i=1,4
    ts(i) = t + sum(as(:i-1,i)) * h
    xs(i) = x + dot_product(as(:i-1,i), ks(:i-1))
    ks(i) = h * (fi(xs(i)) + gi(xs(i))*warray(i))
  enddo
 
  xstar = x + dot_product(as(:,5), ks)
end subroutine

虽然在不了解 fi 和 gi 的情况下很难分辨。

另请注意，您似乎没有使用 t1 到 t4 变量。

stochastic_rk Fortran 90 代码还有进一步优化的余地吗？

Is there room to further optimize the stochastic_rk Fortran 90 code?

algorithm

optimization

fortran

stochastic

differential-equations