C.hr eulers.f
C....*...1.........2.........3.........4.........5.........6.........7.*
C    
C     eulers.f
C
C     Copyright (C) 1996.
C
C     A. Ronald Gallant, P.O. Box 659, Chapel Hill NC 27514-0659, USA
C
C     Permission to use, copy, modify, and distribute this software and 
C     its documentation for any purpose and without fee is hereby 
C     granted, provided that the above copyright notice appear in all
C     copies and that both that copyright notice and this permission 
C     notice appear in supporting documentation.
C     
C     This software is provided "as is" without any expressed or implied 
C     warranty.
C
C....*...1.........2.........3.........4.........5.........6.........7.*
C.hr eulers
C@
*....*...1.........2.........3.........4.........5.........6.........7.*
*
*     eulers     8/10/96
*
*     purpose
*     compute a realization of length n from an Euler-Maruyama scheme
*     for the nonautonomous Ito stochastic differential equation
*
*       dX = A(t,X)dt + B(t,X)dW
*
*     where X has dimension d and W has dimension m.
*
*     If iterations exceed a tolerance, then the process is stopped and 
*     iterations are put to last acceptable value of X thereafter.
*
*     usage
*     call eulers(drift,difuse,t,X,A,B,d,m,t0,X0,dt,iseed,n,Y,Xtst,Xtol)
*
*     libf calls
*     unsk
*
*     arguments
*     drift  - subroutine with usage call drift(t,X,d,A) that is
*              declared external in the calling program.
*              input, external
*     difuse - subroutine with usage call difuse(t,X,d,m,B) that is
*              declared external in the calling program.
*              input, external
*     t      - input to drift and difuse, on return t=t0+n*dt.
*              workspace, real*8
*     X      - input to drift and difuse, a vector of length d.  on
*              return X(i)=Y(i,n) for i=1,...,d.
*              workspace, real*8
*     A      - output of drift, a vector of length d.
*              workspace, real*8
*     B      - output of difuse, a matrix of order d by m.
*              workspace, real*8
*     d      - dimension of X and A, row dimension of B and Y.  d must
*              be less than parameter md which is currently set to 15.
*              input, integer*4
*     m      - column dimension of B.  
*              input, integer*4
*     t0     - time of initial condition.
*              input, real*8
*     X0     - initial condition at time t=t0, a vector of length d.
*              input, real*8
*     dt     - time increment for the simulations.
*              input, real*8
*     iseed  - seed for the simulations.
*              input and output, integer*4
*     n      - number of simulations, column dimension of Y.
*              input, integer*4
*     Y      - computed simulations, a d by n matrix.
*              output, real*8
*     Xtst   - value of the state variable to be used with the 
*              tolerence check, a vector of length d.
*              input, real*8
*     Xtol   - if |(X-Xtst)/Xtst|>Xtol, then interations are put to X
*              thereafer.  
*              input, real*8
*
*     remarks
*     To reduce memory requirements for a computation such as
*     (1/T)(integral from 0 to T of func[X(t)]) where T = n*dt,
*     dimension as Y(d,n0) and use
*       do k=1,n,n0
*         call eulers
*     &     (drift,difuse,t0,X0,A,B,d,m,t0,X0,dt,iseed,n0,Y,Xtst,Xtol)
*         do it=1,n0
*           do i=1,d
*             ave(i) = ave(i) + func(Y(i,it))/n
*           end do
*         end do
*       end do
*     instead of dimensioning as Y(d,n) and using
*       call eulers
*    &    (drift,difuse,t,X,A,B,d,m,t0,X0,dt,iseed,n,Y,Xtst,Xtol)
*       do it=1,n
*         do i=1,d
*           ave(i) = ave(i) + func(Y(i,it))/n
*         end do
*       end do
*     
*     reference
*     Kloeden, Peter E., and Eckhard Platen (1992), Numerical Solution
*     of Stochastic Differential Equations, Berlin, Springer-Verlag,
*     p. 341. 
*
      subroutine 
     &  eulers(drift,difuse,t,X,A,B,d,m,t0,X0,dt,iseed,n,Y,Xtst,Xtol)
      implicit none

*     arguments

      external drift, difuse

      real*8 X, A, B, X0, Y, Xtst
      real*8 t, t0, dt, Xtol

      integer*4 d, m, iseed, n

      dimension X(d), A(d), B(d,m), X0(d), Y(d,n), Xtst(d)

*     parameters

      integer*4 md
      integer*4 nout
      parameter (md=15)
      parameter (nout=3)

*     local variables

      intrinsic DSQRT

      real*4 z,unsk

      real*8 dW
      real*8 sum

      real*8 rdt
      real*8 ratio, rlow, rhigh

      logical toobig

      integer*4 i,j,it

      dimension sum(md)

*     error traps

      if (d.gt.md) then
        write(nout,'(a,i5)') 'Error, eulers, md < d =', d
        stop
      end if

*     Euler-Maruyama scheme

      rdt  = DSQRT(dt)

      do i=1,d
        X(i)=X0(i)
      end do

      t = t0

      rlow  = 1.d0 - Xtol
      rhigh = 1.d0 + Xtol

      do it=1,n

        toobig = .FALSE.

        do i=1,d
          ratio=X(i)/Xtst(i)
          if ((ratio.gt.rhigh).or.(ratio.lt.rlow)) toobig = .TRUE.
        end do

        if (toobig) then
          do j=it,n
            do i=1,d
               Y(i,j)=X(i)
            end do
	    t = t + dt
          end do
          return
        end if

        call drift(t,X,d,A)
        call difuse(t,X,d,m,B)

        do i=1,d
          sum(i) = 0.d0
        end do

        do j=1,m
          z = unsk(iseed)
          dW = rdt*z
          do i=1,d
            sum(i) = sum(i) + B(i,j)*dW
          end do
        end do

        do i=1,d
          Y(i,it) = X(i) + A(i)*dt + sum(i)
        end do

*       For a stiff system the above loop is replaced by an equation
*       solver to solve Y = X + A(it,Y)dt + sum for Y.  See p. 407 of
*       Kloeden and Platen (1992).  Subroutine drift should be revised
*       to return first derivatives with respect to its argument X to
*       allow use of a solver that uses analytic first derivatives.
 
        do i=1,d
          X(i) = Y(i,it)
        end do

        t = t + dt

      end do
 
      return
      end