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

Why heap array make my code 6 times slower?

CRquantum
New Contributor I
‎09-16-2021 09:14 PM
491 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
       integer, private, parameter :: i4=selected_int_kind(9)
       integer, private, parameter :: i8=selected_int_kind(15)
       integer, private, parameter :: 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
       integer, private, parameter :: i4=selected_int_kind(9)
       integer, private, parameter :: i8=selected_int_kind(15)
       integer, private, parameter :: r8=selected_real_kind(15,9)
       integer(kind=i8), private, parameter :: mask48 = ishft(1_8,48)-1
       integer(kind=i8), private, save :: 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
	    integer, public, parameter :: i4=selected_int_kind(9)
	    integer, public, parameter :: i8=selected_int_kind(15)
	    integer, public, parameter :: r8=selected_real_kind(15,9)
	    real(kind=r8), public, parameter :: 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), public, parameter :: 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 ), parameter, dimension(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 ), parameter, dimension(4) :: qs = [ 2.12709852335625_r8 &
		    ,2.73245878238737_r8 &
		    ,11.22760917474960_r8 &
		    ,13.36199560336697_r8 ]
	    real ( kind = r8 ), parameter, dimension(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 ), parameter, dimension(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)
Labels:
0 Kudos

Link Copied

0 Replies
Reply

Community support is provided during standard business hours (Monday to Friday 7AM - 5PM PST). Other contact methods are available here.

Intel does not verify all solutions, including but not limited to any file transfers that may appear in this community. Accordingly, Intel disclaims all express and implied warranties, including without limitation, the implied warranties of merchantability, fitness for a particular purpose, and non-infringement, as well as any warranty arising from course of performance, course of dealing, or usage in trade.

For more complete information about compiler optimizations, see our Optimization Notice.