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,
Thank you very much!
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,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
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.