Intel® Fortran Compiler
Build applications that can scale for the future with optimized code designed for Intel® Xeon® and compatible processors.
27303 Discussions

Why heap array make my code 6 times slower?

CRquantum
New Contributor I
213 Views

Dear all,

 

I am using Intel Fortran OneAPI. 

A similar and detailed topic I put it here, https://fortran-lang.discourse.group/t/why-heap-array-make-my-code-6-times-slower/1858

 

But my intension is not to duplicate the topic. Since I am basically always using Intel Fortran, I would like to hear Intel's experts suggestions too. 


I have a code will show you below, which is based on John Burkardt’s code. John’s version is scaler version can solve 1 SDE. I just made it can slove 2 SDEs. There is plenty room to improve the code but it is not the topic here. 

See the code below. You can directly compile and run it.

Thing is, if I enable ‘heap array’, then the code is slow, like 8-9 seconds.
in windows, flag is
/heap-arrays0
linux is,
-heap-arrays

 

If I do not heap arrays, it only tool 1.5s or so. 

 

Do you why heap-arrays can make the code so slow?

 

Thank you very much in advance!

 

    module random2
       implicit none
       integerprivateparameter :: i4=selected_int_kind(9)
       integerprivateparameter :: i8=selected_int_kind(15)
       integerprivateparameter :: r8=selected_real_kind(15,9)
 
    contains
       subroutine ran1(rn,irn)
    ! multiplicative congruential with additive constant
       integer(kind=i8),  parameter :: mask24 = ishft(1_8,24)-1
       integer(kind=i8),  parameter :: mask48 = ishft(1_8,48_8)-1_8
       real(kind=r8),  parameter :: twom48=2.0d0**(-48)
       integer(kind=i8),  parameter :: mult1 = 44485709377909_8
       integer(kind=i8),  parameter :: m11 = iand(mult1,mask24)
       integer(kind=i8),  parameter :: m12 = iand(ishft(mult1,-24),mask24)
       integer(kind=i8),  parameter :: iadd1 = 96309754297_8
       integer(kind=i8) :: irn
       real(kind=r8) :: rn
       integer(kind=i8) :: is1,is2
       is2 = iand(ishft(irn,-24),mask24)
       is1 = iand(irn,mask24)
       irn = iand(ishft(iand(is1*m12+is2*m11,mask24),24)+is1*m11+iadd1,mask48)
       rn = ior(irn,1_8)*twom48
       return
       end subroutine ran1
 
       subroutine ran2(rn,irn)
    ! multiplicative congruential with additive constant
       integer(kind=i8),  parameter :: mask24 = ishft(1_8,24)-1
       integer(kind=i8), parameter :: mask48 = ishft(1_8,48)-1
       real(kind=r8), parameter :: twom48=2.0d0**(-48)
       integer(kind=i8), parameter :: mult2 = 34522712143931_8
       integer(kind=i8), parameter :: m21 = iand(mult2,mask24)
       integer(kind=i8), parameter :: m22 = iand(ishft(mult2,-24),mask24)
       integer(kind=i8), parameter :: iadd2 = 55789347517_8
       integer(kind=i8) :: irn
       real(kind=r8) :: rn
       integer(kind=i8) :: is1,is2
       is2 = iand(ishft(irn,-24),mask24)
       is1 = iand(irn,mask24)
       irn = iand(ishft(iand(is1*m22+is2*m21,mask24),24)+is1*m21+iadd2,mask48)
       rn = ior(irn,1_8)*twom48
       return
       end subroutine ran2
   
    end module random2
 
    module random
       implicit none
       integerprivateparameter :: i4=selected_int_kind(9)
       integerprivateparameter :: i8=selected_int_kind(15)
       integerprivateparameter :: r8=selected_real_kind(15,9)
       integer(kind=i8), privateparameter :: mask48 = ishft(1_8,48)-1
       integer(kind=i8), privatesave :: irn = 1_8,irnsave
 
    contains   
       function randn(n)
    ! return an array of random variates (0,1)
       use random2
       integer(kind=i8) :: n
       real(kind=r8), dimension(n) :: randn
       integer(kind=i8) :: i
       do i=1,n
	      call ran1(randn(i),irn)
       enddo
       return
       end function randn
 
       function gaussian(n)
    ! return an array of gaussian random variates with variance 1
       integer(kind=i8) :: n
       real(kind=r8), dimension(n) :: gaussian
       real(kind=r8), dimension(2*((n+1)/2)) :: rn
       real(kind=r8) :: rn1,rn2,arg,x1,x2,x3
       real(kind=r8), parameter :: pi=4.0d0*atan(1.0_r8)
       integer(kind=i8) :: i
       rn=randn(size(rn,kind=i8))
       do i=1,n/2
	      rn1=rn(2*i-1)
	      rn2=rn(2*i)
	      arg=2.0_r8*pi*rn1
	      x1=sin(arg)
	      x2=cos(arg)
	      x3=sqrt(-2.0_r8*log(rn2))
	      gaussian(2*i-1)=x1*x3
	      gaussian(2*i)=x2*x3
       enddo
       if (mod(n,2).ne.0) then
	      rn1=rn(n)
	      rn2=rn(n+1)
	      arg=2.0_r8*pi*rn1
	      x2=cos(arg)
	      x3=sqrt(-2.0_r8*log(rn2))
	      gaussian(n)=x2*x3
	      endif
       return
       end function gaussian
 
       subroutine setrn(irnin)
    ! set the seed
       integer(kind=i8) :: irnin
       integer(kind=i8) :: n
       integer(kind=i8), allocatable :: iseed(:)
       irn=iand(irnin,mask48)
       call savern(irn)
       call RANDOM_SEED(SIZE=n)
       allocate(iseed(n))
       iseed=int(irnin)
       call RANDOM_SEED(PUT=iseed) ! set the internal ran function's seed.
       return
       end subroutine setrn
 
       subroutine savern(irnin)
    ! save the seed
       integer(kind=i8) :: irnin
       irnsave=irnin
       return
       end subroutine savern
   
       subroutine showirn(irnout)
    ! show the seed
       integer(kind=i8) :: irnout
       irnout=irnsave
       return
       end subroutine showirn   
   
       function randomn(n) ! the default ramdom_number subroutine.
    ! return an array of random variates (0,1)
       integer(kind=i8) :: n
       real(kind=r8), dimension(n) :: randomn
       call RANDOM_NUMBER(randomn)
       return
       end function randomn
   
       subroutine random1(x) ! generate random number x using the default ramdom_number subroutine.
       real(kind=r8) :: x
       call RANDOM_NUMBER(x)
       return
       end subroutine random1
 
    end module random
    module constants
	    implicit none
	    integerpublicparameter :: i4=selected_int_kind(9)
	    integerpublicparameter :: i8=selected_int_kind(15)
	    integerpublicparameter :: r8=selected_real_kind(15,9)
	    real(kind=r8), publicparameter :: zero=0.0_r8,one=1.0_r8,two=2.0_r8,three=3.0_r8,four=4.0_r8 &
		    ,five=5.0_r8,six=6.0_r8,seven=7.0_r8,eight=8.0_r8,nine=9.0_r8 &
		    ,ten=10.0_r8,tenth=.1_r8,half=.5_r8,third=1.0_r8/3.0_r8,sixth=1.0_r8/6.0_r8 &
		    ,pi=4.0_r8*atan(1.0_r8) &
		    ,normpdf_factor=1.0_r8/sqrt(8.0_r8*atan(1.0_r8)) ! 1/sqrt(2 pi)
	    complex(kind=r8), publicparameter :: czero=(0.0_r8,0.0_r8),ci=(0.0_r8,1.0_r8),cone=(1.0_r8,0.0_r8)  
    end module constants
    module stochastic_rk
    use constants  
    implicit none
    contains    
    subroutine rk4_ti_step_mod ( x, t, h, q, fi, gi, seed, xstar )
    use random
 
    !*****************************************************************************80
    !
    !! RK4_TI_STEP takes one step of a stochastic Runge Kutta scheme.
    !
    !  Discussion:
    !
    !    The Runge-Kutta scheme is fourth-order, and suitable for time-invariant
    !    systems in which F and G do not depend explicitly on time.
    !
    !    d/dx X(t,xsi) = F ( X(t,xsi) ) + G ( X(t,xsi) ) * w(t,xsi)
    !
    !  Licensing:
    !
    !    This code is distributed under the GNU LGPL license.
    !
    !  Modified:
    !
    !    20 June 2010
    !
    !  Author:
    !
    !    John Burkardt
    !
    !  Reference:
    !
    !    Jeremy Kasdin,
    !    Runge-Kutta algorithm for the numerical integration of
    !    stochastic differential equations,
    !    Journal of Guidance, Control, and Dynamics,
    !    Volume 18, Number 1, January-February 1995, pages 114-120.
    !
    !    Jeremy Kasdin,
    !    Discrete Simulation of Colored Noise and Stochastic Processes
    !    and 1/f^a Power Law Noise Generation,
    !    Proceedings of the IEEE,
    !    Volume 83, Number 5, 1995, pages 802-827.
    !
    !  Parameters:
    !
    !    Input, real ( kind = 8 ) X, the value at the current time.
    !
    !    Input, real ( kind = 8 ) T, the current time.
    !
    !    Input, real ( kind = 8 ) H, the time step.
    !
    !    Input, real ( kind = 8 ) Q, the spectral density of the input white noise.
    !
    !    Input, external real ( kind = 8 ) FI, the name of the deterministic
    !    right hand side function.
    !
    !    Input, external real ( kind = 8 ) GI, the name of the stochastic
    !    right hand side function.
    !
    !    Input/output, integer ( kind = 4 ) SEED, a seed for the random 
    !    number generator.
    !
    !    Output, real ( kind = 8 ) XSTAR, the value at time T+H.
    !
      implicit none
 
      real ( kind = r8 ), external :: fi
      real ( kind = r8 ), external :: gi
      real ( kind = r8 ) h
      real ( kind = r8 ) k1
      real ( kind = r8 ) k2
      real ( kind = r8 ) k3
      real ( kind = r8 ) k4
      real ( kind = r8 ) q
      real ( kind = r8 ) r8_normal_01
      integer ( kind = i8 ) seed
      real ( kind = r8 ) t
      real ( kind = r8 ) t1
      real ( kind = r8 ) t2
      real ( kind = r8 ) t3
      real ( kind = r8 ) t4
      real ( kind = r8 ) w1
      real ( kind = r8 ) w2
      real ( kind = r8 ) w3
      real ( kind = r8 ) w4
      real ( kind = r8 ) x
      real ( kind = r8 ) x1
      real ( kind = r8 ) x2
      real ( kind = r8 ) x3
      real ( kind = r8 ) x4
      real ( kind = r8 ) xstar
  
      real ( kind = r8 ) :: qoh
      real ( kind = r8 ) :: normal(4)
  
      real ( kind = r8 ), parameter :: a21 = 2.71644396264860_r8 &
      ,a31 = - 6.95653259006152_r8 &
      ,a32 =   0.78313689457981_r8 &
      ,a41 =   0.0_r8 &
      ,a42 =   0.48257353309214_r8 &
      ,a43 =   0.26171080165848_r8 &
      ,a51 =   0.47012396888046_r8 &
      ,a52 =   0.36597075368373_r8 &
      ,a53 =   0.08906615686702_r8 &
      ,a54 =   0.07483912056879_r8
  
      real ( kind = r8 ), parameterdimension(4) :: qarray = [ 2.12709852335625_r8 &
	      ,2.73245878238737_r8 &
	      ,11.22760917474960_r8 &
	      ,13.36199560336697_r8 ]
  
      real ( kind = r8 ) :: 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 + a43 * k3 ! the a43 * k3 seems need to be added.
      k4 = h * ( fi ( x4 ) + gi ( x4 ) * warray(4) )
 
      xstar = x1 + a51 * k1 + a52 * k2 + a53 * k3 + a54 * k4
 
    return
    end subroutine rk4_ti_step_mod
 
    subroutine rk4_ti_step_vec_mod ( x, t, h, q, fi, gi, n, xstar )
	    use random
	    implicit none
	    interface				   !! Declare function        
		    function fi ( x ) result (f)
		    use constants
		    real ( kind = r8 ) :: f(2)
		    real ( kind = r8 ) :: x(2)
		    end function
		    function gi ( x ) result (g)
		    use constants
		    real ( kind = r8 ) :: g(2)
		    real ( kind = r8 ) :: x(2)
		    end function      
	    end interface
	    integer ( kind = i8 ) n
	    real ( kind = r8 ) h
	    real ( kind = r8 ) k1(n)
	    real ( kind = r8 ) k2(n)
	    real ( kind = r8 ) k3(n)
	    real ( kind = r8 ) k4(n)
	    real ( kind = r8 ) q
	    real ( kind = r8 ) r8_normal_01
	
	    real ( kind = r8 ) t
	    real ( kind = r8 ) t1
	    real ( kind = r8 ) t2
	    real ( kind = r8 ) t3
	    real ( kind = r8 ) t4
	    real ( kind = r8 ) w1
	    real ( kind = r8 ) w2
	    real ( kind = r8 ) w3
	    real ( kind = r8 ) w4
	    real ( kind = r8 ) x(n)
	    real ( kind = r8 ) x1(n)
	    real ( kind = r8 ) x2(n)
	    real ( kind = r8 ) x3(n)
	    real ( kind = r8 ) x4(n)
	    real ( kind = r8 ) xstar(n)
   
	    real ( kind = r8 ), parameter :: a21 = 2.71644396264860_r8 &
	    ,a31 = - 6.95653259006152_r8 &
	    ,a32 =   0.78313689457981_r8 &
	    ,a41 =   0.0_r8 &
	    ,a42 =   0.48257353309214_r8 &
	    ,a43 =   0.26171080165848_r8 &
	    ,a51 =   0.47012396888046_r8 &
	    ,a52 =   0.36597075368373_r8 &
	    ,a53 =   0.08906615686702_r8 &
	    ,a54 =   0.07483912056879_r8
  
	    real ( kind = r8 ), parameterdimension(4) :: qs = [ 2.12709852335625_r8 &
		    ,2.73245878238737_r8 &
		    ,11.22760917474960_r8 &
		    ,13.36199560336697_r8 ]
	    real ( kind = r8 ), parameterdimension(2) :: aaa = [one,one]
	
	    real ( kind = r8 ) :: qoh
	    real ( kind = r8 ) :: warray(n,4)
	    integer (kind = i8) :: i
  
	    qoh = q / h
  
	    do i =1,4
		    warray(:,i) = gaussian(n)*sqrt(qs(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 + a43 * k3 ! the a43 * k3 seems need to be added.
	    k4 = h * ( fi ( x4 ) + gi ( x4 ) * warray(:,4) )
 
	    xstar = x1 + a51 * k1 + a52 * k2 + a53 * k3 + a54 * k4
 
    return
    end subroutine rk4_ti_step_vec_mod
 
 
	    subroutine rk4_ti_step_vec_mod1 ( x, t, h, q, fi_vec, gi_vec, n, xstar )
		    use random
		    implicit none
		    interface				   !! Declare function        
			    function fi_vec ( x ) result (f)
			    use constants
			    real ( kind = r8 ) :: f(2)
			    real ( kind = r8 ) :: x(2)
			    end function
			    function gi_vec ( x ) result (g)
			    use constants
			    real ( kind = r8 ) :: g(2)
			    real ( kind = r8 ) :: x(2)
			    end function      
		    end interface       
		    integer(kind = i8), intent(in) :: n
		    real(kind=r8), intent(in) :: x(n)
		    real(kind = r8), intent(in) :: t  
		    real(kind = r8), intent(in) :: h
		    real(kind = r8), intent(in) :: q
		    !real(kind = r8), external :: fi_vec
		    !real(kind = r8), external :: gi_vec
		    real(kind = r8), intent(out) :: xstar(n)
		    real ( kind = r8 ), parameter :: alphas(4) = [ 0.47012396888046_r8 &
		    ,0.36597075368373_r8 &
		    ,0.08906615686702_r8 &
		    ,0.07483912056879_r8 ]
		    real(kind = r8), parameter :: as(4,4) = reshape([ &
			    &  0.0_r8,              0.0_r8,              0.0_r8,              0.0_r8, &
			    &  2.71644396264860_r8, 0.0_r8,              0.0_r8,              0.0_r8, &
			    & -6.95653259006152_r8, 0.78313689457981_r8, 0.0_r8,              0.0_r8, &
			    &  0.0_r8,              0.48257353309214_r8, 0.26171080165848_r8, 0.0_r8 ], [4,4])
		    real(kind = r8), parameter :: qs(4) = [  2.12709852335625_r8, &
			    &  2.73245878238737_r8, &
			    & 11.22760917474960_r8, &
			    & 13.36199560336697_r8 ]
		    real(kind = r8) :: ks(n,4)
		    real(kind = r8) :: ts(n,4)
		    real(kind = r8) :: ws(n,4)
		    real(kind = r8) :: xs(n,4)
		    real(kind = r8) :: warray(n,4)
		    real(kind = r8) :: ak(n)
		    integer(kind = i8) :: j,i
		    real(kind = r8) :: qoh
		    qoh = q / h
		    xstar = x    
		    do j=1,4
			    warray(:,j) = gaussian(n)*sqrt(qs(j)*qoh)
			    !ts(j) = t + sum(as(:j-1,j)) * h      
			    do concurrent (i=1:n)
				    ak(i)=dot_product(as(:j-1,j), ks(i,:j-1))  
			    enddo        
			    !ak=zero
			    !do i=1,j-1
			    !    ak = ak + as(i,j)*ks(:,i)
			    !enddo
			    xs(:,j) = x + ak  
			    ks(:,j) = h * (  fi_vec(xs(:,j)) + gi_vec(xs(:,j))*warray(:,j)  )
			    xstar = xstar + alphas(j)*ks(:,j)     
		    enddo
	    return
	    end subroutine rk4_ti_step_vec_mod1
 
    subroutine rk4_ti_step_mod1 ( x, t, h, q, fi, gi, seed, xstar ) 
    ! mostly suggested by veryreverie from stackoverflow, 
    ! https://stackoverflow.com/questions/69147944/is-there-room-to-further-optimize-the-stochastic-rk-fortran-90-code?noredirect=1#comment122226757_69147944
	    use random  
	    implicit none
	    real(kind=r8), intent(in) :: x
	    real(kind = r8), intent(in) :: t  
	    real(kind = r8), intent(in) :: h
	    real(kind = r8), intent(in) :: q
	    real(kind = r8), external :: fi
	    real(kind = r8), external :: gi
	    integer(kind = i8), intent(in) :: seed
	    real(kind = r8), intent(out) :: xstar
	    real ( kind = r8 ), parameter :: alphas(4) = [ 0.47012396888046_r8 &
	    ,0.36597075368373_r8 &
	    ,0.08906615686702_r8 &
	    ,0.07483912056879_r8 ]
	    real(kind = r8), parameter :: as(4,4) = reshape([ &
		    &  0.0_r8,              0.0_r8,              0.0_r8,              0.0_r8, &
		    &  2.71644396264860_r8, 0.0_r8,              0.0_r8,              0.0_r8, &
		    & -6.95653259006152_r8, 0.78313689457981_r8, 0.0_r8,              0.0_r8, &
		    &  0.0_r8,              0.48257353309214_r8, 0.26171080165848_r8, 0.0_r8 ], [4,4])
	    real(kind = r8), parameter :: qs(4) = [ &
		    &  2.12709852335625_r8, &
		    &  2.73245878238737_r8, &
		    & 11.22760917474960_r8, &
		    & 13.36199560336697_r8 ]
	    real(kind = r8) :: ks(4)
	    real(kind = r8) :: r8_normal_01
	    real(kind = r8) :: ts(4)
	    real(kind = r8) :: ws(4)
	    real(kind = r8) :: xs(4)
	    real(kind = r8) :: normal(4)
	    real(kind = r8) :: warray(4)
	    integer(kind = i8) :: i
	    real(kind = r8) :: qoh
	    qoh = q/h
	    normal = gaussian(4)  
	    xstar = x
	    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))
	    warray(i) = normal(i)*sqrt(qs(i)*qoh)
	    ks(i) = h * (fi(xs(i)) + gi(xs(i))*warray(i))
	    xstar = xstar + alphas(i)*ks(i)
	    enddo
    return
    end subroutine rk4_ti_step_mod1
 
 
    subroutine rk4_tv_step ( x, t, h, q, fv, gv, seed, xstar )
 
    !*****************************************************************************80
    !
    !! RK4_TV_STEP takes one step of a stochastic Runge Kutta scheme.
    !
    !  Discussion:
    !
    !    The Runge-Kutta scheme is fourth-order, and suitable for time-varying
    !    systems.
    !
    !    d/dx X(t,xsi) = F ( X(t,xsi), t ) + G ( X(t,xsi), t ) * w(t,xsi)
    !
    !  Licensing:
    !
    !    This code is distributed under the GNU LGPL license.
    !
    !  Modified:
    !
    !    08 June 2010
    !
    !  Author:
    !
    !    John Burkardt
    !
    !  Reference:
    !
    !    Jeremy Kasdin,
    !    Runge-Kutta algorithm for the numerical integration of
    !    stochastic differential equations,
    !    Journal of Guidance, Control, and Dynamics,
    !    Volume 18, Number 1, January-February 1995, pages 114-120.
    !
    !    Jeremy Kasdin,
    !    Discrete Simulation of Colored Noise and Stochastic Processes
    !    and 1/f^a Power Law Noise Generation,
    !    Proceedings of the IEEE,
    !    Volume 83, Number 5, 1995, pages 802-827.
    !
    !  Parameters:
    !
    !    Input, real ( kind = 8 ) X, the value at the current time.
    !
    !    Input, real ( kind = 8 ) T, the current time.
    !
    !    Input, real ( kind = 8 ) H, the time step.
    !
    !    Input, real ( kind = 8 ) Q, the spectral density of the input white noise.
    !
    !    Input, external real ( kind = 8 ) FV, the name of the deterministic
    !    right hand side function.
    !
    !    Input, external real ( kind = 8 ) GV, the name of the stochastic
    !    right hand side function.
    !
    !    Input/output, integer ( kind = 4 ) SEED, a seed for the random 
    !    number generator.
    !
    !    Output, real ( kind = 8 ) XSTAR, the value at time T+H.
    !
      implicit none
 
      real ( kind = r8 ) a21
      real ( kind = r8 ) a31
      real ( kind = r8 ) a32
      real ( kind = r8 ) a41
      real ( kind = r8 ) a42
      real ( kind = r8 ) a43
      real ( kind = r8 ) a51
      real ( kind = r8 ) a52
      real ( kind = r8 ) a53
      real ( kind = r8 ) a54
      real ( kind = r8 ), external :: fv
      real ( kind = r8 ), external :: gv
      real ( kind = r8 ) h
      real ( kind = r8 ) k1
      real ( kind = r8 ) k2
      real ( kind = r8 ) k3
      real ( kind = r8 ) k4
      real ( kind = r8 ) q
      real ( kind = r8 ) q1
      real ( kind = r8 ) q2
      real ( kind = r8 ) q3
      real ( kind = r8 ) q4
      real ( kind = r8 ) r8_normal_01
      integer ( kind = i8 ) seed
      real ( kind = r8 ) t
      real ( kind = r8 ) t1
      real ( kind = r8 ) t2
      real ( kind = r8 ) t3
      real ( kind = r8 ) t4
      real ( kind = r8 ) w1
      real ( kind = r8 ) w2
      real ( kind = r8 ) w3
      real ( kind = r8 ) w4
      real ( kind = r8 ) x
      real ( kind = r8 ) x1
      real ( kind = r8 ) x2
      real ( kind = r8 ) x3
      real ( kind = r8 ) x4
      real ( kind = r8 ) xstar
 
      a21 =   0.66667754298442_r8
      a31 =   0.63493935027993_r8
      a32 =   0.00342761715422_r8
      a41 = - 2.32428921184321_r8
      a42 =   2.69723745129487_r8
      a43 =   0.29093673271592_r8
      a51 =   0.25001351164789_r8
      a52 =   0.67428574806272_r8
      a53 = - 0.00831795169360_r8
      a54 =   0.08401868181222_r8
 
      q1 = 3.99956364361748_r8
      q2 = 1.64524970733585_r8
      q3 = 1.59330355118722_r8
      q4 = 0.26330006501868_r8
 
      t1 = t
      x1 = x
      !w1 = r8_normal_01 ( seed ) * sqrt ( q1 * q / h )
      k1 = h * fv ( t1, x1 ) + h * gv ( t1, x1 ) * w1 
 
      t2 = t1 + a21 * h
      x2 = x1 + a21 * k1
      !w2 = r8_normal_01 ( seed ) * sqrt ( q2 * q / h )
      k2 = h * fv ( t2, x2 ) + h * gv ( t2, x2 ) * w2
 
      t3 = t1 + a31 * h  + a32 * h
      x3 = x1 + a31 * k1 + a32 * k2
      !w3 = r8_normal_01 ( seed ) * sqrt ( q3 * q / h )
      k3 = h * fv ( t3, x3 ) + h * gv ( t3, x3 ) * w3
 
      t4 = t1 + a41 * h  + a42 * h  + a43 * h
      x4 = x1 + a41 * k1 + a42 * k2 + a43 * k3
      !w4 = r8_normal_01 ( seed ) * sqrt ( q4 * q / h )
      k4 = h * fv ( t4, x4 ) + h * gv ( t4, x4 ) * w4
 
      xstar = x1 + a51 * k1 + a52 * k2 + a53 * k3 + a54 * k4
 
      return
    end
    subroutine timestamp ( )
 
    !*****************************************************************************80
    !
    !! TIMESTAMP prints the current YMDHMS date as a time stamp.
    !
    !  Example:
    !
    !    31 May 2001   9:45:54.872 AM
    !
    !  Licensing:
    !
    !    This code is distributed under the GNU LGPL license.
    !
    !  Modified:
    !
    !    18 May 2013
    !
    !  Author:
    !
    !    John Burkardt
    !
    !  Parameters:
    !
    !    None
    !
      implicit none
 
      character ( len = 8 ) ampm
      integer ( kind = 4 ) d
      integer ( kind = 4 ) h
      integer ( kind = 4 ) m
      integer ( kind = 4 ) mm
      character ( len = 9 ), parameterdimension(12) :: month = (/ &
	    'January  ''February ''March    ''April    ', &
	    'May      ''June     ''July     ''August   ', &
	    'September''October  ''November ''December ' /)
      integer ( kind = 4 ) n
      integer ( kind = 4 ) s
      integer ( kind = 4 ) values(8)
      integer ( kind = 4 ) y
 
      call date_and_time ( values = values )
 
      y = values(1)
      m = values(2)
      d = values(3)
      h = values(5)
      n = values(6)
      s = values(7)
      mm = values(8)
 
      if ( h < 12 ) then
	    ampm = 'AM'
      else if ( h == 12 ) then
	    if ( n == 0 .and. s == 0 ) then
	      ampm = 'Noon'
	    else
	      ampm = 'PM'
	    end if
      else
	    h = h - 12
	    if ( h < 12 ) then
	      ampm = 'PM'
	    else if ( h == 12 ) then
	      if ( n == 0 .and. s == 0 ) then
		    ampm = 'Midnight'
	      else
		    ampm = 'AM'
	      end if
	    end if
      end if
 
      write ( *, '(i2.2,1x,a,1x,i4,2x,i2,a1,i2.2,a1,i2.2,a1,i3.3,1x,a)' ) &
	    d, trim ( month(m) ), y, h, ':', n, ':', s, '.', mm, trim ( ampm )
 
      return
    end
    end module stochastic_rk
    module stats
    use constants 
    implicit none
    contains
    subroutine stat2(x,xave,xvar)
    use constants
    real (kind=r8) :: x(:,:),val(5),val2(5),xvar(5),xave(5)
    integer (kind =i8) :: i,j,xsize,np
    val=zero
    val2=zero
    np = size(x,2)
    do i=1,np
	    do j=1,5
		    val(j) = val(j) + x(j,i)
		    val2(j) = val2(j) + x(j,i)**2
	    enddo
    enddo
    xave=val/np
    xvar=abs( val2/np - xave**2 ) 
    return
    end subroutine stat2 
 
    subroutine stat(x,xave,xvar)
    use constants
    real (kind=r8) :: x(:),val,val2,xvar,xave
    integer (kind =i8) :: i,xsize,np
    val=0.0
    val2=0.0
    np = size(x)
    do i=1,np
	    val = val + x(i)
	    val2= val2 + x(i)**2
	    !write (6,*) i, x(i)
    enddo
    xave=val/np
    xvar=abs(  val2/np - xave**2  )   
    end subroutine stat 
	
    end module stats
 
	
	    program main
	    use random
	    use random2
	    use constants
	    use stochastic_rk
	    use stats
 
	    !*****************************************************************************80
	    !
	    !! MAIN is the main program for STOCHASTIC_RK_TEST.
	    !
	    !  Discussion:
	    !
	    !    STOCHASTIC_RK_TEST tests the STOCHASTIC_RK library.
	    !
	    !  Licensing:
	    !
	    !    This code is distributed under the GNU LGPL license.
	    !
	    !  Modified:
	    !
	    !    20 June 2010
	    !
	    !  Author:
	    !
	    !    John Burkardt
	    !
	    implicit none
  
	    integer(kind=i8) :: i,np,seed,k
	    integer(kind=i8) :: day,hour,minute
	    real (kind=r8) :: t0,t1,t2
	    real (kind=r8), allocatable :: x(:,:)
	    real (kind=r8) :: xvar(5),xave(5)
 
	    call timestamp ( )
	    write ( *, '(a)' ) ' '
	    write ( *, '(a)' ) 'STOCHASTIC_RK_TEST'
	    write ( *, '(a)' ) '  FORTRAN90 version'
	    write ( *, '(a)' ) '  Test the STOCHASTIC_RK library.'
 
  
	    seed = 123456789
	    call setrn(int(seed,8))
  
	    np = 10000
	    allocate(x(5,np))
  
	    call CPU_TIME(t1) 
	    do i=1,np
	    call test02 ( x(:,i) )   
	    !write (6,*) i, x(i) 
	    enddo
	    call CPU_TIME(t2)  
	    write(6,'(''Time cost = '',t20,g12.5,'' sec'')') t2-t1  
	    ! statistics
	    call stat2(x,xave,xvar)
	    write(6,'(''step size : '',t20,g12.5)''1/1000'
	    write(6,'(''Np: '',t20,g12.5)') np
	    write(6,'('' '')')
	    write(6,'(''Simulated Mean at t=0.2:1.0:0.2 is: '',t50,5(g12.5,1x))') (xave(k),k=1,5)
	    write(6,'(''Theoretical Mean at t=0.2:1.0:0.2 is: '',t50,5(g12.5,1x))') ( 20.0*exp(-k*0.2),k=1,5 )
	    write(6,'('' '')')
	    write(6,'(''Simulated Variance at t=0.2:1.0:0.2 is: '',t50,5(g12.5,1x))') (xvar(k),k=1,5)
	    write(6,'(''Theoretical Variance at t=0.2:1.0:0.2 is: '',t50,5(g12.5,1x))') ( 0.5*(1.0-exp(-2*k*0.2)),k=1,5 )
  
	    !
	    !  Terminate.
	    !
	    call timestamp ( )
	    stop ('Program end normally.')
 
	    end program main 
	
 
	    subroutine test02 ( xout )
	    use constants
	    use stochastic_rk
	    implicit none
	    interface				   !! Declare function        
	    function fi_vec2 ( x ) result (f)
	    use constants
	    real ( kind = r8 ) :: f(2)
	    real ( kind = r8 ) :: x(2)
	    end function
	    function gi_vec2 ( x ) result (g)
	    use constants
	    real ( kind = r8 ) :: g(2)
	    real ( kind = r8 ) :: x(2)
	    end function      
	    end interface 
	    integer ( kind = i8 ), parameter :: n = 1000
	    real ( kind = r8 ) h
	    integer ( kind = i8 ) :: i, itot, istart, istep
	    real ( kind = r8 ) q
	    integer ( kind = i8 ) seed
	    real ( kind = r8 ) t
	    real ( kind = r8 ), parameter :: t0 = 0.0
	    real ( kind = r8 ), parameter :: tn = 1.0
	    real ( kind = r8 ) x(2,0:n)
	    real (kind = r8) :: xout(5)
	    integer ( kind = i8 ) :: nd ! dimension, =2
 
	    h = ( tn - t0 ) / real(n)
	    q = 1.0
	    seed = 123456789
 
	    i = 0
	    t = t0
	    x(1,i) = 20.0_r8
	    x(2,i) = one
	    nd = 2
 
	    !write ( *, '(a)' ) ' '
	    !write ( *, '(a)' ) '         I           T             X'
	    !write ( *, '(a)' ) ' '
	    !write ( *, '(2x,i8,2x,f14.6,2x,g14.6)' ) i, t, x(i)
 
	    do i = 1, n
	    t = t0 + (real(i)*(tn-t0))/real(n)
	    call rk4_ti_step_vec_mod1 ( x(:,i-1), t, h, q, fi_vec2, gi_vec2, nd, x(:,i) )
	    !write ( *, '(2x,i8,2x,f14.6,2x,g14.6)' ) i, t, x(i)
	    end do
  
	    itot = i-1
	    istart = itot/5
	    istep = istart
	    xout(1:5:1) = x(1,istart:itot:istep)
	    !write ( *, '(2x,i8,2x,f14.6,2x,g14.6)' ) i-1, t, xout
	    return
	    end subroutine test02   
	
	    function fi_vec2 ( x ) result (f)
	    use constants
	    implicit none
	    real ( kind = r8 ) :: f(2)
	    real ( kind = r8 ) :: x(2)
	      f(1) = -x(2)*x(1) 
	      f(2) = -x(2) + one
	    return
	    end function fi_vec2
	    function gi_vec2 ( x ) result(g)
	    use constants
	    implicit none
	    real ( kind = r8 ), parameter :: gp(2) = [one,zero] 
	    real ( kind = r8 ) :: x(2),g(2)
	    g=gp
	    return
	    end function gi_vec2    
	
	    subroutine test01 ( xout )
	    use constants
	    use stochastic_rk
	    implicit none
	    integer ( kind = i8 ), parameter :: n = 1000
	    real ( kind = r8 ), external :: fi
	    real ( kind = r8 ), external :: gi
	    real ( kind = r8 ) h
	    integer ( kind = i8 ) :: i, itot, istart, istep
	    real ( kind = r8 ) q
	    integer ( kind = i8 ) seed
	    real ( kind = r8 ) t
	    real ( kind = r8 ), parameter :: t0 = 0.0
	    real ( kind = r8 ), parameter :: tn = 1.0
	    real ( kind = r8 ) x(0:n)
	    real (kind = r8) :: xout(5)
 
	    h = ( tn - t0 ) / real(n)
	    q = 1.0
	    seed = 123456789
 
	    i = 0
	    t = t0
	    x(i) = 20.0
 
	    !write ( *, '(a)' ) ' '
	    !write ( *, '(a)' ) '         I           T             X'
	    !write ( *, '(a)' ) ' '
	    !write ( *, '(2x,i8,2x,f14.6,2x,g14.6)' ) i, t, x(i)
 
	    do i = 1, n
	    t = t0 + (real(i)*(tn-t0))/real(n)
	    call rk4_ti_step_mod ( x(i-1), t, h, q, fi, gi, seed, x(i) )
	    !write ( *, '(2x,i8,2x,f14.6,2x,g14.6)' ) i, t, x(i)
	    end do
  
	    itot = i-1
	    istart = itot/5
	    istep = istart
	    xout(1:5:1) = x(istart:itot:istep)
	    !write ( *, '(2x,i8,2x,f14.6,2x,g14.6)' ) i-1, t, xout
	    return
	    end subroutine test01
	
	
	    function fi ( x )
	    use constants
 
	    !*****************************************************************************80
	    !
	    !! FI is a time invariant deterministic right hand side.
	    !
	    !  Licensing:
	    !
	    !    This code is distributed under the GNU LGPL license.
	    !
	    !  Modified:
	    !
	    !    20 June 2010
	    !
	    !  Author:
	    !
	    !    John Burkardt
	    !
	    !  Parameters:
	    !
	    !    Input, real ( kind = 8 ) X, the argument.
	    !
	    !    Output, real ( kind = 8 ) FI, the value.
	    !
	      implicit none
 
	      real ( kind = r8 ) fi
	      real ( kind = r8 ) x
 
	      fi = -x  !1.0_r8
 
	      return
	    end function fi
	    function gi ( x )
	    use constants
 
	    !*****************************************************************************80
	    !
	    !! GI is a time invariant stochastic right hand side.
	    !
	    !  Licensing:
	    !
	    !    This code is distributed under the GNU LGPL license.
	    !
	    !  Modified:
	    !
	    !    20 June 2010
	    !
	    !  Author:
	    !
	    !    John Burkardt
	    !
	    !  Parameters:
	    !
	    !    Input, real ( kind = 8 ) X, the argument.
	    !
	    !    Output, real ( kind = 8 ) GI, the value.
	    !
	      implicit none
 
	      real ( kind = r8 ) gi
	      real ( kind = r8 ) x
 
	      gi = 1.0
 
	      return
	    end function gi
Labels (1)
0 Kudos
0 Replies
Reply