I need to use a Fortran code to solve stochastic differential equation (SDE). I looked at the famous Fortran code website by Burkardt,

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

I particular looked at the rk4_ti_step subroutine in stochastic_rk.f90 code,

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

My optimized version is below,

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

Note that I use my module of random number, and gaussian is my random number function, this part does not matter.

I just wonder,

1. Can anyone give some suggestions as to can the code be further optimized?
2. Does anyone know what is the best/fastest SDE Fortran subroutine? Or what algorithm is the best?

Thank you very much!

CodePudding user response：

The interdependence of x and c means you can't turn as much into linear algebra as I first thought, but I'd still expect some speedup by grouping everything into appropriate arrays as:

subroutine rk4_ti_step_mod ( x, t, h, q, fi, gi, seed, xstar )
  use random  
  implicit none
  
  integer, parameter :: dp = selected_real_kind(15,300)
  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

although it's difficult to tell without knowing anything about fi and gi.

Also note you don't seem to be using the t1 to t4 variables.