是否有进一步优化stochastic_rk Fortran 90代码的空间？

Question

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

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

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

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

请注意，我使用了我的随机数模块，而高斯是我的随机数函数，这部分并不重要。

我只是好奇，

对于代码是否可以进一步成为optimized?

• Does，任何人都知道什么是最好的/最快的subroutine子程序，有人能给出一些建议吗？或者哪种算法是最好的？

非常感谢!

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变量。