sstode.f

Go to the documentation of this file.

      SUBROUTINE sstode (NEQ, Y, YH, NYH, YH1, EWT, SAVF, ACOR,
     1   WM, IWM, F, JAC, PJAC, SLVS)
C***BEGIN PROLOGUE  SSTODE
C***SUBSIDIARY
C***PURPOSE  Performs one step of an ODEPACK integration.
C***TYPE      SINGLE PRECISION (SSTODE-S, DSTODE-D)
C***AUTHOR  Hindmarsh, Alan C., (LLNL)
C***DESCRIPTION
C
C  SSTODE performs one step of the integration of an initial value
C  problem for a system of ordinary differential equations.
C  Note:  SSTODE is independent of the value of the iteration method
C  indicator MITER, when this is .ne. 0, and hence is independent
C  of the type of chord method used, or the Jacobian structure.
C  Communication with SSTODE is done with the following variables:
C
C  NEQ    = integer array containing problem size in NEQ(1), and
C           passed as the NEQ argument in all calls to F and JAC.
C  Y      = an array of length .ge. N used as the Y argument in
C           all calls to F and JAC.
C  YH     = an NYH by LMAX array containing the dependent variables
C           and their approximate scaled derivatives, where
C           LMAX = MAXORD + 1.  YH(i,j+1) contains the approximate
C           j-th derivative of y(i), scaled by h**j/factorial(j)
C           (j = 0,1,...,NQ).  on entry for the first step, the first
C           two columns of YH must be set from the initial values.
C  NYH    = a constant integer .ge. N, the first dimension of YH.
C  YH1    = a one-dimensional array occupying the same space as YH.
C  EWT    = an array of length N containing multiplicative weights
C           for local error measurements.  Local errors in Y(i) are
C           compared to 1.0/EWT(i) in various error tests.
C  SAVF   = an array of working storage, of length N.
C           Also used for input of YH(*,MAXORD+2) when JSTART = -1
C           and MAXORD .lt. the current order NQ.
C  ACOR   = a work array of length N, used for the accumulated
C           corrections.  On a successful return, ACOR(i) contains
C           the estimated one-step local error in Y(i).
C  WM,IWM = real and integer work arrays associated with matrix
C           operations in chord iteration (MITER .ne. 0).
C  PJAC   = name of routine to evaluate and preprocess Jacobian matrix
C           and P = I - h*el0*JAC, if a chord method is being used.
C  SLVS   = name of routine to solve linear system in chord iteration.
C  CCMAX  = maximum relative change in h*el0 before PJAC is called.
C  H      = the step size to be attempted on the next step.
C           H is altered by the error control algorithm during the
C           problem.  H can be either positive or negative, but its
C           sign must remain constant throughout the problem.
C  HMIN   = the minimum absolute value of the step size h to be used.
C  HMXI   = inverse of the maximum absolute value of h to be used.
C           HMXI = 0.0 is allowed and corresponds to an infinite hmax.
C           HMIN and HMXI may be changed at any time, but will not
C           take effect until the next change of h is considered.
C  TN     = the independent variable. TN is updated on each step taken.
C  JSTART = an integer used for input only, with the following
C           values and meanings:
C                0  perform the first step.
C            .gt.0  take a new step continuing from the last.
C               -1  take the next step with a new value of H, MAXORD,
C                     N, METH, MITER, and/or matrix parameters.
C               -2  take the next step with a new value of H,
C                     but with other inputs unchanged.
C           On return, JSTART is set to 1 to facilitate continuation.
C  KFLAG  = a completion code with the following meanings:
C                0  the step was succesful.
C               -1  the requested error could not be achieved.
C               -2  corrector convergence could not be achieved.
C               -3  fatal error in PJAC or SLVS.
C           A return with KFLAG = -1 or -2 means either
C           abs(H) = HMIN or 10 consecutive failures occurred.
C           On a return with KFLAG negative, the values of TN and
C           the YH array are as of the beginning of the last
C           step, and H is the last step size attempted.
C  MAXORD = the maximum order of integration method to be allowed.
C  MAXCOR = the maximum number of corrector iterations allowed.
C  MSBP   = maximum number of steps between PJAC calls (MITER .gt. 0).
C  MXNCF  = maximum number of convergence failures allowed.
C  METH/MITER = the method flags.  See description in driver.
C  N      = the number of first-order differential equations.
C  The values of CCMAX, H, HMIN, HMXI, TN, JSTART, KFLAG, MAXORD,
C  MAXCOR, MSBP, MXNCF, METH, MITER, and N are communicated via COMMON.
C
C***SEE ALSO  SLSODE
C***ROUTINES CALLED  SCFODE, SVNORM
C***COMMON BLOCKS    SLS001
C***REVISION HISTORY  (YYMMDD)
C   791129  DATE WRITTEN
C   890501  Modified prologue to SLATEC/LDOC format.  (FNF)
C   890503  Minor cosmetic changes.  (FNF)
C   930809  Renamed to allow single/double precision versions. (ACH)
C   010413  Reduced size of Common block /SLS001/. (ACH)
C   031105  Restored 'own' variables to Common block /SLS001/, to
C           enable interrupt/restart feature. (ACH)
C***END PROLOGUE  SSTODE
C**End
      EXTERNAL f, jac, pjac, slvs
      INTEGER NEQ, NYH, IWM
      REAL Y, YH, YH1, EWT, SAVF, ACOR, WM
      dimension neq(*), y(*), yh(nyh,*), yh1(*), ewt(*), savf(*),
     1   acor(*), wm(*), iwm(*)
      INTEGER INIT, MXSTEP, MXHNIL, NHNIL, NSLAST, CNYH,
     1   ialth, ipup, lmax, meo, nqnyh, nslp,
     1   icf, ierpj, iersl, jcur, jstart, kflag, l,
     2   lyh, lewt, lacor, lsavf, lwm, liwm, meth, miter,
     3   maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
      INTEGER I, I1, IREDO, IRET, J, JB, M, NCF, NEWQ
      REAL CONIT, CRATE, EL, ELCO, HOLD, RMAX, TESCO,
     1   ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround
      REAL DCON, DDN, DEL, DELP, DSM, DUP, EXDN, EXSM, EXUP,
     1   r, rh, rhdn, rhsm, rhup, told, svnorm
      COMMON /sls001/ conit, crate, el(13), elco(13,12),
     1   hold, rmax, tesco(3,12),
     1   ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround,
     2   init, mxstep, mxhnil, nhnil, nslast, cnyh,
     3   ialth, ipup, lmax, meo, nqnyh, nslp,
     3   icf, ierpj, iersl, jcur, jstart, kflag, l,
     4   lyh, lewt, lacor, lsavf, lwm, liwm, meth, miter,
     5   maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
C
C***FIRST EXECUTABLE STATEMENT  SSTODE
      kflag = 0
      told = tn
      ncf = 0
      ierpj = 0
      iersl = 0
      jcur = 0
      icf = 0
      delp = 0.0e0
      IF (jstart .GT. 0) GO TO 200
      IF (jstart .EQ. -1) GO TO 100
      IF (jstart .EQ. -2) GO TO 160
C-----------------------------------------------------------------------
C On the first call, the order is set to 1, and other variables are
C initialized.  RMAX is the maximum ratio by which H can be increased
C in a single step.  It is initially 1.E4 to compensate for the small
C initial H, but then is normally equal to 10.  If a failure
C occurs (in corrector convergence or error test), RMAX is set to 2
C for the next increase.
C-----------------------------------------------------------------------
      lmax = maxord + 1
      nq = 1
      l = 2
      ialth = 2
      rmax = 10000.0e0
      rc = 0.0e0
      el0 = 1.0e0
      crate = 0.7e0
      hold = h
      meo = meth
      nslp = 0
      ipup = miter
      iret = 3
      GO TO 140
C-----------------------------------------------------------------------
C The following block handles preliminaries needed when JSTART = -1.
C IPUP is set to MITER to force a matrix update.
C If an order increase is about to be considered (IALTH = 1),
C IALTH is reset to 2 to postpone consideration one more step.
C If the caller has changed METH, SCFODE is called to reset
C the coefficients of the method.
C If the caller has changed MAXORD to a value less than the current
C order NQ, NQ is reduced to MAXORD, and a new H chosen accordingly.
C If H is to be changed, YH must be rescaled.
C If H or METH is being changed, IALTH is reset to L = NQ + 1
C to prevent further changes in H for that many steps.
C-----------------------------------------------------------------------
 100  ipup = miter
      lmax = maxord + 1
      IF (ialth .EQ. 1) ialth = 2
      IF (meth .EQ. meo) GO TO 110
      CALL scfode (meth, elco, tesco)
      meo = meth
      IF (nq .GT. maxord) GO TO 120
      ialth = l
      iret = 1
      GO TO 150
 110  IF (nq .LE. maxord) GO TO 160
 120  nq = maxord
      l = lmax
      DO 125 i = 1,l
 125    el(i) = elco(i,nq)
      nqnyh = nq*nyh
      rc = rc*el(1)/el0
      el0 = el(1)
      conit = 0.5e0/(nq+2)
      ddn = svnorm(n, savf, ewt)/tesco(1,l)
      exdn = 1.0e0/l
      rhdn = 1.0e0/(1.3e0*ddn**exdn + 0.0000013e0)
      rh = min(rhdn,1.0e0)
      iredo = 3
      IF (h .EQ. hold) GO TO 170
      rh = min(rh,abs(h/hold))
      h = hold
      GO TO 175
C-----------------------------------------------------------------------
C SCFODE is called to get all the integration coefficients for the
C current METH.  Then the EL vector and related constants are reset
C whenever the order NQ is changed, or at the start of the problem.
C-----------------------------------------------------------------------
 140  CALL scfode (meth, elco, tesco)
 150  DO 155 i = 1,l
 155    el(i) = elco(i,nq)
      nqnyh = nq*nyh
      rc = rc*el(1)/el0
      el0 = el(1)
      conit = 0.5e0/(nq+2)
      GO TO (160, 170, 200), iret
C-----------------------------------------------------------------------
C If H is being changed, the H ratio RH is checked against
C RMAX, HMIN, and HMXI, and the YH array rescaled.  IALTH is set to
C L = NQ + 1 to prevent a change of H for that many steps, unless
C forced by a convergence or error test failure.
C-----------------------------------------------------------------------
 160  IF (h .EQ. hold) GO TO 200
      rh = h/hold
      h = hold
      iredo = 3
      GO TO 175
 170  rh = max(rh,hmin/abs(h))
 175  rh = min(rh,rmax)
      rh = rh/max(1.0e0,abs(h)*hmxi*rh)
      r = 1.0e0
      DO 180 j = 2,l
        r = r*rh
        DO 180 i = 1,n
 180      yh(i,j) = yh(i,j)*r
      h = h*rh
      rc = rc*rh
      ialth = l
      IF (iredo .EQ. 0) GO TO 690
C-----------------------------------------------------------------------
C This section computes the predicted values by effectively
C multiplying the YH array by the Pascal Triangle matrix.
C RC is the ratio of new to old values of the coefficient  H*EL(1).
C When RC differs from 1 by more than CCMAX, IPUP is set to MITER
C to force PJAC to be called, if a Jacobian is involved.
C In any case, PJAC is called at least every MSBP steps.
C-----------------------------------------------------------------------
 200  IF (abs(rc-1.0e0) .GT. ccmax) ipup = miter
      IF (nst .GE. nslp+msbp) ipup = miter
      tn = tn + h
      i1 = nqnyh + 1
      DO 215 jb = 1,nq
        i1 = i1 - nyh
Cdir$ ivdep
        DO 210 i = i1,nqnyh
 210      yh1(i) = yh1(i) + yh1(i+nyh)
 215    CONTINUE
C-----------------------------------------------------------------------
C Up to MAXCOR corrector iterations are taken.  A convergence test is
C made on the R.M.S. norm of each correction, weighted by the error
C weight vector EWT.  The sum of the corrections is accumulated in the
C vector ACOR(i).  The YH array is not altered in the corrector loop.
C-----------------------------------------------------------------------
 220  m = 0
      DO 230 i = 1,n
 230    y(i) = yh(i,1)
      CALL f (neq, tn, y, savf)
      nfe = nfe + 1
      IF (ipup .LE. 0) GO TO 250
C-----------------------------------------------------------------------
C If indicated, the matrix P = I - h*el(1)*J is reevaluated and
C preprocessed before starting the corrector iteration.  IPUP is set
C to 0 as an indicator that this has been done.
C-----------------------------------------------------------------------
      CALL pjac (neq, y, yh, nyh, ewt, acor, savf, wm, iwm, f, jac)
      ipup = 0
      rc = 1.0e0
      nslp = nst
      crate = 0.7e0
      IF (ierpj .NE. 0) GO TO 430
 250  DO 260 i = 1,n
 260    acor(i) = 0.0e0
 270  IF (miter .NE. 0) GO TO 350
C-----------------------------------------------------------------------
C In the case of functional iteration, update Y directly from
C the result of the last function evaluation.
C-----------------------------------------------------------------------
      DO 290 i = 1,n
        savf(i) = h*savf(i) - yh(i,2)
 290    y(i) = savf(i) - acor(i)
      del = svnorm(n, y, ewt)
      DO 300 i = 1,n
        y(i) = yh(i,1) + el(1)*savf(i)
 300    acor(i) = savf(i)
      GO TO 400
C-----------------------------------------------------------------------
C In the case of the chord method, compute the corrector error,
C and solve the linear system with that as right-hand side and
C P as coefficient matrix.
C-----------------------------------------------------------------------
 350  DO 360 i = 1,n
 360    y(i) = h*savf(i) - (yh(i,2) + acor(i))
      CALL slvs (wm, iwm, y, savf)
      IF (iersl .LT. 0) GO TO 430
      IF (iersl .GT. 0) GO TO 410
      del = svnorm(n, y, ewt)
      DO 380 i = 1,n
        acor(i) = acor(i) + y(i)
 380    y(i) = yh(i,1) + el(1)*acor(i)
C-----------------------------------------------------------------------
C Test for convergence.  If M.gt.0, an estimate of the convergence
C rate constant is stored in CRATE, and this is used in the test.
C-----------------------------------------------------------------------
 400  IF (m .NE. 0) crate = max(0.2e0*crate,del/delp)
      dcon = del*min(1.0e0,1.5e0*crate)/(tesco(2,nq)*conit)
      IF (dcon .LE. 1.0e0) GO TO 450
      m = m + 1
      IF (m .EQ. maxcor) GO TO 410
      IF (m .GE. 2 .AND. del .GT. 2.0e0*delp) GO TO 410
      delp = del
      CALL f (neq, tn, y, savf)
      nfe = nfe + 1
      GO TO 270
C-----------------------------------------------------------------------
C The corrector iteration failed to converge.
C If MITER .ne. 0 and the Jacobian is out of date, PJAC is called for
C the next try.  Otherwise the YH array is retracted to its values
C before prediction, and H is reduced, if possible.  If H cannot be
C reduced or MXNCF failures have occurred, exit with KFLAG = -2.
C-----------------------------------------------------------------------
 410  IF (miter .EQ. 0 .OR. jcur .EQ. 1) GO TO 430
      icf = 1
      ipup = miter
      GO TO 220
 430  icf = 2
      ncf = ncf + 1
      rmax = 2.0e0
      tn = told
      i1 = nqnyh + 1
      DO 445 jb = 1,nq
        i1 = i1 - nyh
Cdir$ ivdep
        DO 440 i = i1,nqnyh
 440      yh1(i) = yh1(i) - yh1(i+nyh)
 445    CONTINUE
      IF (ierpj .LT. 0 .OR. iersl .LT. 0) GO TO 680
      IF (abs(h) .LE. hmin*1.00001e0) GO TO 670
      IF (ncf .EQ. mxncf) GO TO 670
      rh = 0.25e0
      ipup = miter
      iredo = 1
      GO TO 170
C-----------------------------------------------------------------------
C The corrector has converged.  JCUR is set to 0
C to signal that the Jacobian involved may need updating later.
C The local error test is made and control passes to statement 500
C if it fails.
C-----------------------------------------------------------------------
 450  jcur = 0
      IF (m .EQ. 0) dsm = del/tesco(2,nq)
      IF (m .GT. 0) dsm = svnorm(n, acor, ewt)/tesco(2,nq)
      IF (dsm .GT. 1.0e0) GO TO 500
C-----------------------------------------------------------------------
C After a successful step, update the YH array.
C Consider changing H if IALTH = 1.  Otherwise decrease IALTH by 1.
C If IALTH is then 1 and NQ .lt. MAXORD, then ACOR is saved for
C use in a possible order increase on the next step.
C If a change in H is considered, an increase or decrease in order
C by one is considered also.  A change in H is made only if it is by a
C factor of at least 1.1.  If not, IALTH is set to 3 to prevent
C testing for that many steps.
C-----------------------------------------------------------------------
      kflag = 0
      iredo = 0
      nst = nst + 1
      hu = h
      nqu = nq
      DO 470 j = 1,l
        DO 470 i = 1,n
 470      yh(i,j) = yh(i,j) + el(j)*acor(i)
      ialth = ialth - 1
      IF (ialth .EQ. 0) GO TO 520
      IF (ialth .GT. 1) GO TO 700
      IF (l .EQ. lmax) GO TO 700
      DO 490 i = 1,n
 490    yh(i,lmax) = acor(i)
      GO TO 700
C-----------------------------------------------------------------------
C The error test failed.  KFLAG keeps track of multiple failures.
C Restore TN and the YH array to their previous values, and prepare
C to try the step again.  Compute the optimum step size for this or
C one lower order.  After 2 or more failures, H is forced to decrease
C by a factor of 0.2 or less.
C-----------------------------------------------------------------------
 500  kflag = kflag - 1
      tn = told
      i1 = nqnyh + 1
      DO 515 jb = 1,nq
        i1 = i1 - nyh
Cdir$ ivdep
        DO 510 i = i1,nqnyh
 510      yh1(i) = yh1(i) - yh1(i+nyh)
 515    CONTINUE
      rmax = 2.0e0
      IF (abs(h) .LE. hmin*1.00001e0) GO TO 660
      IF (kflag .LE. -3) GO TO 640
      iredo = 2
      rhup = 0.0e0
      GO TO 540
C-----------------------------------------------------------------------
C Regardless of the success or failure of the step, factors
C RHDN, RHSM, and RHUP are computed, by which H could be multiplied
C at order NQ - 1, order NQ, or order NQ + 1, respectively.
C In the case of failure, RHUP = 0.0 to avoid an order increase.
C The largest of these is determined and the new order chosen
C accordingly.  If the order is to be increased, we compute one
C additional scaled derivative.
C-----------------------------------------------------------------------
 520  rhup = 0.0e0
      IF (l .EQ. lmax) GO TO 540
      DO 530 i = 1,n
 530    savf(i) = acor(i) - yh(i,lmax)
      dup = svnorm(n, savf, ewt)/tesco(3,nq)
      exup = 1.0e0/(l+1)
      rhup = 1.0e0/(1.4e0*dup**exup + 0.0000014e0)
 540  exsm = 1.0e0/l
      rhsm = 1.0e0/(1.2e0*dsm**exsm + 0.0000012e0)
      rhdn = 0.0e0
      IF (nq .EQ. 1) GO TO 560
      ddn = svnorm(n, yh(1,l), ewt)/tesco(1,nq)
      exdn = 1.0e0/nq
      rhdn = 1.0e0/(1.3e0*ddn**exdn + 0.0000013e0)
 560  IF (rhsm .GE. rhup) GO TO 570
      IF (rhup .GT. rhdn) GO TO 590
      GO TO 580
 570  IF (rhsm .LT. rhdn) GO TO 580
      newq = nq
      rh = rhsm
      GO TO 620
 580  newq = nq - 1
      rh = rhdn
      IF (kflag .LT. 0 .AND. rh .GT. 1.0e0) rh = 1.0e0
      GO TO 620
 590  newq = l
      rh = rhup
      IF (rh .LT. 1.1e0) GO TO 610
      r = el(l)/l
      DO 600 i = 1,n
 600    yh(i,newq+1) = acor(i)*r
      GO TO 630
 610  ialth = 3
      GO TO 700
 620  IF ((kflag .EQ. 0) .AND. (rh .LT. 1.1e0)) GO TO 610
      IF (kflag .LE. -2) rh = min(rh,0.2e0)
C-----------------------------------------------------------------------
C If there is a change of order, reset NQ, l, and the coefficients.
C In any case H is reset according to RH and the YH array is rescaled.
C Then exit from 690 if the step was OK, or redo the step otherwise.
C-----------------------------------------------------------------------
      IF (newq .EQ. nq) GO TO 170
 630  nq = newq
      l = nq + 1
      iret = 2
      GO TO 150
C-----------------------------------------------------------------------
C Control reaches this section if 3 or more failures have occurred.
C If 10 failures have occurred, exit with KFLAG = -1.
C It is assumed that the derivatives that have accumulated in the
C YH array have errors of the wrong order.  Hence the first
C derivative is recomputed, and the order is set to 1.  Then
C H is reduced by a factor of 10, and the step is retried,
C until it succeeds or H reaches HMIN.
C-----------------------------------------------------------------------
 640  IF (kflag .EQ. -10) GO TO 660
      rh = 0.1e0
      rh = max(hmin/abs(h),rh)
      h = h*rh
      DO 645 i = 1,n
 645    y(i) = yh(i,1)
      CALL f (neq, tn, y, savf)
      nfe = nfe + 1
      DO 650 i = 1,n
 650    yh(i,2) = h*savf(i)
      ipup = miter
      ialth = 5
      IF (nq .EQ. 1) GO TO 200
      nq = 1
      l = 2
      iret = 3
      GO TO 150
C-----------------------------------------------------------------------
C All returns are made through this section.  H is saved in HOLD
C to allow the caller to change H on the next step.
C-----------------------------------------------------------------------
 660  kflag = -1
      GO TO 720
 670  kflag = -2
      GO TO 720
 680  kflag = -3
      GO TO 720
 690  rmax = 10.0e0
 700  r = 1.0e0/tesco(2,nqu)
      DO 710 i = 1,n
 710    acor(i) = acor(i)*r
 720  hold = h
      jstart = 1
      RETURN
C----------------------- END OF SUBROUTINE SSTODE ----------------------
      END