relibc_openlibm/slatec/dstod.f

*DECK DSTOD
      SUBROUTINE DSTOD (NEQ, Y, YH, NYH, YH1, EWT, SAVF, ACOR, WM, IWM,
     +   DF, DJAC, RPAR, IPAR)
C***BEGIN PROLOGUE  DSTOD
C***SUBSIDIARY
C***PURPOSE  Subsidiary to DDEBDF
C***LIBRARY   SLATEC
C***TYPE      DOUBLE PRECISION (STOD-S, DSTOD-D)
C***AUTHOR  Watts, H. A., (SNLA)
C***DESCRIPTION
C
C   DSTOD integrates a system of first order odes over one step in the
C   integrator package DDEBDF.
C ----------------------------------------------------------------------
C DSTOD  performs one step of the integration of an initial value
C problem for a system of ordinary differential equations.
C Note.. DSTOD  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 DSTOD  is done with the following variables..
C
C Y      = An array of length .GE. N used as the Y argument in
C          all calls to DF and DJAC.
C NEQ    = Integer array containing problem size in NEQ(1), and
C          passed as the NEQ argument in all calls to DF and DJAC.
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 N elements with which the estimated local
C          errors in YH are compared.
C SAVF   = An array of working storage, of length N.
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 = DOUBLE PRECISION and INTEGER work arrays associated with
C          matrix operations in chord iteration (MITER .NE. 0).
C DPJAC   = Name of routine to evaluate and preprocess Jacobian matrix
C          if a chord method is being used.
C DSLVS   = Name of routine to solve linear system in chord iteration.
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 successful.
C              -1  The requested error could not be achieved.
C              -2  Corrector convergence could not be achieved.
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 METH/MITER = The method flags.  See description in driver.
C N      = The number of first-order differential equations.
C ----------------------------------------------------------------------
C
C***SEE ALSO  DDEBDF
C***ROUTINES CALLED  DCFOD, DPJAC, DSLVS, DVNRMS
C***COMMON BLOCKS    DDEBD1
C***REVISION HISTORY  (YYMMDD)
C   820301  DATE WRITTEN
C   890531  Changed all specific intrinsics to generic.  (WRB)
C   890911  Removed unnecessary intrinsics.  (WRB)
C   891214  Prologue converted to Version 4.0 format.  (BAB)
C   900328  Added TYPE section.  (WRB)
C   910722  Updated AUTHOR section.  (ALS)
C   920422  Changed DIMENSION statement.  (WRB)
C***END PROLOGUE  DSTOD
C
      INTEGER I, I1, IALTH, IER, IOD, IOWND, IPAR, IPUP, IREDO, IRET,
     1      IWM, J, JB, JSTART, KFLAG, KSTEPS, L, LMAX, M, MAXORD,
     2      MEO, METH, MITER, N, NCF, NEQ, NEWQ, NFE, NJE, NQ, NQNYH,
     3      NQU, NST, NSTEPJ, NYH
      DOUBLE PRECISION ACOR, CONIT, CRATE, DCON, DDN,
     1      DEL, DELP, DSM, DUP, DVNRMS, EL, EL0, ELCO,
     2      EWT, EXDN, EXSM, EXUP, H, HMIN, HMXI, HOLD, HU, R, RC,
     3      RH, RHDN, RHSM, RHUP, RMAX, ROWND, RPAR, SAVF, TESCO,
     4      TN, TOLD, UROUND, WM, Y, YH, YH1
      EXTERNAL DF, DJAC
C
      DIMENSION Y(*),YH(NYH,*),YH1(*),EWT(*),SAVF(*),ACOR(*),WM(*),
     1          IWM(*),RPAR(*),IPAR(*)
      COMMON /DDEBD1/ ROWND,CONIT,CRATE,EL(13),ELCO(13,12),HOLD,RC,RMAX,
     1                TESCO(3,12),EL0,H,HMIN,HMXI,HU,TN,UROUND,IOWND(7),
     2                KSTEPS,IOD(6),IALTH,IPUP,LMAX,MEO,NQNYH,NSTEPJ,
     3                IER,JSTART,KFLAG,L,METH,MITER,MAXORD,N,NQ,NST,NFE,
     4                NJE,NQU
C
C
C     BEGIN BLOCK PERMITTING ...EXITS TO 690
C        BEGIN BLOCK PERMITTING ...EXITS TO 60
C***FIRST EXECUTABLE STATEMENT  DSTOD
            KFLAG = 0
            TOLD = TN
            NCF = 0
            IF (JSTART .GT. 0) GO TO 160
            IF (JSTART .EQ. -1) GO TO 10
               IF (JSTART .EQ. -2) GO TO 90
C              ---------------------------------------------------------
C               ON THE FIRST CALL, THE ORDER IS SET TO 1, AND OTHER
C               VARIABLES ARE INITIALIZED.  RMAX IS THE MAXIMUM RATIO BY
C               WHICH H CAN BE INCREASED IN A SINGLE STEP.  IT IS
C               INITIALLY 1.E4 TO COMPENSATE FOR THE SMALL INITIAL H,
C               BUT THEN IS NORMALLY EQUAL TO 10.  IF A FAILURE OCCURS
C               (IN CORRECTOR CONVERGENCE OR ERROR TEST), RMAX IS SET AT
C               2 FOR THE NEXT INCREASE.
C              ---------------------------------------------------------
               LMAX = MAXORD + 1
               NQ = 1
               L = 2
               IALTH = 2
               RMAX = 10000.0D0
               RC = 0.0D0
               EL0 = 1.0D0
               CRATE = 0.7D0
               DELP = 0.0D0
               HOLD = H
               MEO = METH
               NSTEPJ = 0
               IRET = 3
            GO TO 50
   10       CONTINUE
C              BEGIN BLOCK PERMITTING ...EXITS TO 30
C                 ------------------------------------------------------
C                  THE FOLLOWING BLOCK HANDLES PRELIMINARIES NEEDED WHEN
C                  JSTART = -1.  IPUP IS SET TO MITER TO FORCE A MATRIX
C                  UPDATE.  IF AN ORDER INCREASE IS ABOUT TO BE
C                  CONSIDERED (IALTH = 1), IALTH IS RESET TO 2 TO
C                  POSTPONE CONSIDERATION ONE MORE STEP.  IF THE CALLER
C                  HAS CHANGED METH, DCFOD  IS CALLED TO RESET THE
C                  COEFFICIENTS OF THE METHOD.  IF THE CALLER HAS
C                  CHANGED MAXORD TO A VALUE LESS THAN THE CURRENT
C                  ORDER NQ, NQ IS REDUCED TO MAXORD, AND A NEW H CHOSEN
C                  ACCORDINGLY.  IF H IS TO BE CHANGED, YH MUST BE
C                  RESCALED.  IF H OR METH IS BEING CHANGED, IALTH IS
C                  RESET TO L = NQ + 1 TO PREVENT FURTHER CHANGES IN H
C                  FOR THAT MANY STEPS.
C                 ------------------------------------------------------
                  IPUP = MITER
                  LMAX = MAXORD + 1
                  IF (IALTH .EQ. 1) IALTH = 2
                  IF (METH .EQ. MEO) GO TO 20
                     CALL DCFOD(METH,ELCO,TESCO)
                     MEO = METH
C              ......EXIT
                     IF (NQ .GT. MAXORD) GO TO 30
                     IALTH = L
                     IRET = 1
C        ............EXIT
                     GO TO 60
   20             CONTINUE
                  IF (NQ .LE. MAXORD) GO TO 90
   30          CONTINUE
               NQ = MAXORD
               L = LMAX
               DO 40 I = 1, L
                  EL(I) = ELCO(I,NQ)
   40          CONTINUE
               NQNYH = NQ*NYH
               RC = RC*EL(1)/EL0
               EL0 = EL(1)
               CONIT = 0.5D0/(NQ+2)
               DDN = DVNRMS(N,SAVF,EWT)/TESCO(1,L)
               EXDN = 1.0D0/L
               RHDN = 1.0D0/(1.3D0*DDN**EXDN + 0.0000013D0)
               RH = MIN(RHDN,1.0D0)
               IREDO = 3
               IF (H .EQ. HOLD) GO TO 660
               RH = MIN(RH,ABS(H/HOLD))
               H = HOLD
               GO TO 100
   50       CONTINUE
C           ------------------------------------------------------------
C            DCFOD  IS CALLED TO GET ALL THE INTEGRATION COEFFICIENTS
C            FOR THE CURRENT METH.  THEN THE EL VECTOR AND RELATED
C            CONSTANTS ARE RESET WHENEVER THE ORDER NQ IS CHANGED, OR AT
C            THE START OF THE PROBLEM.
C           ------------------------------------------------------------
            CALL DCFOD(METH,ELCO,TESCO)
   60    CONTINUE
   70    CONTINUE
C           BEGIN BLOCK PERMITTING ...EXITS TO 680
               DO 80 I = 1, L
                  EL(I) = ELCO(I,NQ)
   80          CONTINUE
               NQNYH = NQ*NYH
               RC = RC*EL(1)/EL0
               EL0 = EL(1)
               CONIT = 0.5D0/(NQ+2)
               GO TO (90,660,160), 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
C               IS SET TO L = NQ + 1 TO PREVENT A CHANGE OF H FOR THAT
C               MANY STEPS, UNLESS FORCED BY A CONVERGENCE OR ERROR TEST
C               FAILURE.
C              ---------------------------------------------------------
   90          CONTINUE
               IF (H .EQ. HOLD) GO TO 160
               RH = H/HOLD
               H = HOLD
               IREDO = 3
  100          CONTINUE
  110          CONTINUE
                  RH = MIN(RH,RMAX)
                  RH = RH/MAX(1.0D0,ABS(H)*HMXI*RH)
                  R = 1.0D0
                  DO 130 J = 2, L
                     R = R*RH
                     DO 120 I = 1, N
                        YH(I,J) = YH(I,J)*R
  120                CONTINUE
  130             CONTINUE
                  H = H*RH
                  RC = RC*RH
                  IALTH = L
                  IF (IREDO .NE. 0) GO TO 150
                     RMAX = 10.0D0
                     R = 1.0D0/TESCO(2,NQU)
                     DO 140 I = 1, N
                        ACOR(I) = ACOR(I)*R
  140                CONTINUE
C     ...............EXIT
                     GO TO 690
  150             CONTINUE
C                 ------------------------------------------------------
C                  THIS SECTION COMPUTES THE PREDICTED VALUES BY
C                  EFFECTIVELY MULTIPLYING THE YH ARRAY BY THE PASCAL
C                  TRIANGLE MATRIX.  RC IS THE RATIO OF NEW TO OLD
C                  VALUES OF THE COEFFICIENT  H*EL(1).  WHEN RC DIFFERS
C                  FROM 1 BY MORE THAN 30 PERCENT, IPUP IS SET TO MITER
C                  TO FORCE DPJAC TO BE CALLED, IF A JACOBIAN IS
C                  INVOLVED.  IN ANY CASE, DPJAC IS CALLED AT LEAST
C                  EVERY 20-TH STEP.
C                 ------------------------------------------------------
  160             CONTINUE
  170             CONTINUE
C                    BEGIN BLOCK PERMITTING ...EXITS TO 610
C                       BEGIN BLOCK PERMITTING ...EXITS TO 490
                           IF (ABS(RC-1.0D0) .GT. 0.3D0) IPUP = MITER
                           IF (NST .GE. NSTEPJ + 20) IPUP = MITER
                           TN = TN + H
                           I1 = NQNYH + 1
                           DO 190 JB = 1, NQ
                              I1 = I1 - NYH
                              DO 180 I = I1, NQNYH
                                 YH1(I) = YH1(I) + YH1(I+NYH)
  180                         CONTINUE
  190                      CONTINUE
                           KSTEPS = KSTEPS + 1
C                          ---------------------------------------------
C                           UP TO 3 CORRECTOR ITERATIONS ARE TAKEN.  A
C                           CONVERGENCE TEST IS MADE ON THE R.M.S. NORM
C                           OF EACH CORRECTION, WEIGHTED BY THE ERROR
C                           WEIGHT VECTOR EWT.  THE SUM OF THE
C                           CORRECTIONS IS ACCUMULATED IN THE VECTOR
C                           ACOR(I).  THE YH ARRAY IS NOT ALTERED IN THE
C                           CORRECTOR LOOP.
C                          ---------------------------------------------
  200                      CONTINUE
                              M = 0
                              DO 210 I = 1, N
                                 Y(I) = YH(I,1)
  210                         CONTINUE
                              CALL DF(TN,Y,SAVF,RPAR,IPAR)
                              NFE = NFE + 1
                              IF (IPUP .LE. 0) GO TO 220
C                                ---------------------------------------
C                                 IF INDICATED, THE MATRIX P = I -
C                                 H*EL(1)*J IS REEVALUATED AND
C                                 PREPROCESSED BEFORE STARTING THE
C                                 CORRECTOR ITERATION.  IPUP IS SET TO 0
C                                 AS AN INDICATOR THAT THIS HAS BEEN
C                                 DONE.
C                                ---------------------------------------
                                 IPUP = 0
                                 RC = 1.0D0
                                 NSTEPJ = NST
                                 CRATE = 0.7D0
                                 CALL DPJAC(NEQ,Y,YH,NYH,EWT,ACOR,SAVF,
     1                                      WM,IWM,DF,DJAC,RPAR,IPAR)
C                          ......EXIT
                                 IF (IER .NE. 0) GO TO 440
  220                         CONTINUE
                              DO 230 I = 1, N
                                 ACOR(I) = 0.0D0
  230                         CONTINUE
  240                         CONTINUE
                                 IF (MITER .NE. 0) GO TO 270
C                                   ------------------------------------
C                                    IN THE CASE OF FUNCTIONAL
C                                    ITERATION, UPDATE Y DIRECTLY FROM
C                                    THE RESULT OF THE LAST FUNCTION
C                                    EVALUATION.
C                                   ------------------------------------
                                    DO 250 I = 1, N
                                       SAVF(I) = H*SAVF(I) - YH(I,2)
                                       Y(I) = SAVF(I) - ACOR(I)
  250                               CONTINUE
                                    DEL = DVNRMS(N,Y,EWT)
                                    DO 260 I = 1, N
                                       Y(I) = YH(I,1) + EL(1)*SAVF(I)
                                       ACOR(I) = SAVF(I)
  260                               CONTINUE
                                 GO TO 300
  270                            CONTINUE
C                                   ------------------------------------
C                                    IN THE CASE OF THE CHORD METHOD,
C                                    COMPUTE THE CORRECTOR ERROR, AND
C                                    SOLVE THE LINEAR SYSTEM WITH THAT
C                                    AS RIGHT-HAND SIDE AND P AS
C                                    COEFFICIENT MATRIX.
C                                   ------------------------------------
                                    DO 280 I = 1, N
                                       Y(I) = H*SAVF(I)
     1                                        - (YH(I,2) + ACOR(I))
  280                               CONTINUE
                                    CALL DSLVS(WM,IWM,Y,SAVF)
C                             ......EXIT
                                    IF (IER .NE. 0) GO TO 430
                                    DEL = DVNRMS(N,Y,EWT)
                                    DO 290 I = 1, N
                                       ACOR(I) = ACOR(I) + Y(I)
                                       Y(I) = YH(I,1) + EL(1)*ACOR(I)
  290                               CONTINUE
  300                            CONTINUE
C                                ---------------------------------------
C                                 TEST FOR CONVERGENCE.  IF M.GT.0, AN
C                                 ESTIMATE OF THE CONVERGENCE RATE
C                                 CONSTANT IS STORED IN CRATE, AND THIS
C                                 IS USED IN THE TEST.
C                                ---------------------------------------
                                 IF (M .NE. 0)
     1                              CRATE = MAX(0.2D0*CRATE,DEL/DELP)
                                 DCON = DEL*MIN(1.0D0,1.5D0*CRATE)
     1                                  /(TESCO(2,NQ)*CONIT)
                                 IF (DCON .GT. 1.0D0) GO TO 420
C                                   ------------------------------------
C                                    THE CORRECTOR HAS CONVERGED.  IPUP
C                                    IS SET TO -1 IF MITER .NE. 0, TO
C                                    SIGNAL THAT THE JACOBIAN INVOLVED
C                                    MAY NEED UPDATING LATER.  THE LOCAL
C                                    ERROR TEST IS MADE AND CONTROL
C                                    PASSES TO STATEMENT 500 IF IT
C                                    FAILS.
C                                   ------------------------------------
                                    IF (MITER .NE. 0) IPUP = -1
                                    IF (M .EQ. 0) DSM = DEL/TESCO(2,NQ)
                                    IF (M .GT. 0)
     1                                 DSM = DVNRMS(N,ACOR,EWT)
     2                                       /TESCO(2,NQ)
                                    IF (DSM .GT. 1.0D0) GO TO 380
C                                      BEGIN BLOCK
C                                      PERMITTING ...EXITS TO 360
C                                         ------------------------------
C                                          AFTER A SUCCESSFUL STEP,
C                                          UPDATE THE YH ARRAY.
C                                          CONSIDER CHANGING H IF IALTH
C                                          = 1.  OTHERWISE DECREASE
C                                          IALTH BY 1.  IF IALTH IS THEN
C                                          1 AND NQ .LT. MAXORD, THEN
C                                          ACOR IS SAVED FOR USE IN A
C                                          POSSIBLE ORDER INCREASE ON
C                                          THE NEXT STEP.  IF A CHANGE
C                                          IN H IS CONSIDERED, AN
C                                          INCREASE OR DECREASE IN ORDER
C                                          BY ONE IS CONSIDERED ALSO.  A
C                                          CHANGE IN H IS MADE ONLY IF
C                                          IT IS BY A FACTOR OF AT LEAST
C                                          1.1.  IF NOT, IALTH IS SET TO
C                                          3 TO PREVENT TESTING FOR THAT
C                                          MANY STEPS.
C                                         ------------------------------
                                          KFLAG = 0
                                          IREDO = 0
                                          NST = NST + 1
                                          HU = H
                                          NQU = NQ
                                          DO 320 J = 1, L
                                             DO 310 I = 1, N
                                                YH(I,J) = YH(I,J)
     1                                                    + EL(J)
     2                                                      *ACOR(I)
  310                                        CONTINUE
  320                                     CONTINUE
                                          IALTH = IALTH - 1
                                          IF (IALTH .NE. 0) GO TO 340
C                                            ---------------------------
C                                             REGARDLESS OF THE SUCCESS
C                                             OR FAILURE OF THE STEP,
C                                             FACTORS RHDN, RHSM, AND
C                                             RHUP ARE COMPUTED, BY
C                                             WHICH H COULD BE
C                                             MULTIPLIED AT ORDER NQ -
C                                             1, ORDER NQ, OR ORDER NQ +
C                                             1, RESPECTIVELY.  IN THE
C                                             CASE OF FAILURE, RHUP =
C                                             0.0 TO AVOID AN ORDER
C                                             INCREASE.  THE LARGEST OF
C                                             THESE IS DETERMINED AND
C                                             THE NEW ORDER CHOSEN
C                                             ACCORDINGLY.  IF THE ORDER
C                                             IS TO BE INCREASED, WE
C                                             COMPUTE ONE ADDITIONAL
C                                             SCALED DERIVATIVE.
C                                            ---------------------------
                                             RHUP = 0.0D0
C                       .....................EXIT
                                             IF (L .EQ. LMAX) GO TO 490
                                             DO 330 I = 1, N
                                                SAVF(I) = ACOR(I)
     1                                                    - YH(I,LMAX)
  330                                        CONTINUE
                                             DUP = DVNRMS(N,SAVF,EWT)
     1                                             /TESCO(3,NQ)
                                             EXUP = 1.0D0/(L+1)
                                             RHUP = 1.0D0
     1                                              /(1.4D0*DUP**EXUP
     2                                                + 0.0000014D0)
C                       .....................EXIT
                                             GO TO 490
  340                                     CONTINUE
C                                      ...EXIT
                                          IF (IALTH .GT. 1) GO TO 360
C                                      ...EXIT
                                          IF (L .EQ. LMAX) GO TO 360
                                          DO 350 I = 1, N
                                             YH(I,LMAX) = ACOR(I)
  350                                     CONTINUE
  360                                  CONTINUE
                                       R = 1.0D0/TESCO(2,NQU)
                                       DO 370 I = 1, N
                                          ACOR(I) = ACOR(I)*R
  370                                  CONTINUE
C     .................................EXIT
                                       GO TO 690
  380                               CONTINUE
C                                   ------------------------------------
C                                    THE ERROR TEST FAILED.  KFLAG KEEPS
C                                    TRACK OF MULTIPLE FAILURES.
C                                    RESTORE TN AND THE YH ARRAY TO
C                                    THEIR PREVIOUS VALUES, AND PREPARE
C                                    TO TRY THE STEP AGAIN.  COMPUTE THE
C                                    OPTIMUM STEP SIZE FOR THIS OR ONE
C                                    LOWER ORDER.  AFTER 2 OR MORE
C                                    FAILURES, H IS FORCED TO DECREASE
C                                    BY A FACTOR OF 0.2 OR LESS.
C                                   ------------------------------------
                                    KFLAG = KFLAG - 1
                                    TN = TOLD
                                    I1 = NQNYH + 1
                                    DO 400 JB = 1, NQ
                                       I1 = I1 - NYH
                                       DO 390 I = I1, NQNYH
                                          YH1(I) = YH1(I) - YH1(I+NYH)
  390                                  CONTINUE
  400                               CONTINUE
                                    RMAX = 2.0D0
                                    IF (ABS(H) .GT. HMIN*1.00001D0)
     1                                 GO TO 410
C                                      ---------------------------------
C                                       ALL RETURNS ARE MADE THROUGH
C                                       THIS SECTION.  H IS SAVED IN
C                                       HOLD TO ALLOW THE CALLER TO
C                                       CHANGE H ON THE NEXT STEP.
C                                      ---------------------------------
                                       KFLAG = -1
C     .................................EXIT
                                       GO TO 690
  410                               CONTINUE
C                    ...............EXIT
                                    IF (KFLAG .LE. -3) GO TO 610
                                    IREDO = 2
                                    RHUP = 0.0D0
C                       ............EXIT
                                    GO TO 490
  420                            CONTINUE
                                 M = M + 1
C                             ...EXIT
                                 IF (M .EQ. 3) GO TO 430
C                             ...EXIT
                                 IF (M .GE. 2 .AND. DEL .GT. 2.0D0*DELP)
     1                              GO TO 430
                                 DELP = DEL
                                 CALL DF(TN,Y,SAVF,RPAR,IPAR)
                                 NFE = NFE + 1
                              GO TO 240
  430                         CONTINUE
C                             ------------------------------------------
C                              THE CORRECTOR ITERATION FAILED TO
C                              CONVERGE IN 3 TRIES.  IF MITER .NE. 0 AND
C                              THE JACOBIAN IS OUT OF DATE, DPJAC IS
C                              CALLED FOR THE NEXT TRY.  OTHERWISE THE
C                              YH ARRAY IS RETRACTED TO ITS VALUES
C                              BEFORE PREDICTION, AND H IS REDUCED, IF
C                              POSSIBLE.  IF H CANNOT BE REDUCED OR 10
C                              FAILURES HAVE OCCURRED, EXIT WITH KFLAG =
C                              -2.
C                             ------------------------------------------
C                          ...EXIT
                              IF (IPUP .EQ. 0) GO TO 440
                              IPUP = MITER
                           GO TO 200
  440                      CONTINUE
                           TN = TOLD
                           NCF = NCF + 1
                           RMAX = 2.0D0
                           I1 = NQNYH + 1
                           DO 460 JB = 1, NQ
                              I1 = I1 - NYH
                              DO 450 I = I1, NQNYH
                                 YH1(I) = YH1(I) - YH1(I+NYH)
  450                         CONTINUE
  460                      CONTINUE
                           IF (ABS(H) .GT. HMIN*1.00001D0) GO TO 470
                              KFLAG = -2
C     ........................EXIT
                              GO TO 690
  470                      CONTINUE
                           IF (NCF .NE. 10) GO TO 480
                              KFLAG = -2
C     ........................EXIT
                              GO TO 690
  480                      CONTINUE
                           RH = 0.25D0
                           IPUP = MITER
                           IREDO = 1
C                 .........EXIT
                           GO TO 650
  490                   CONTINUE
                        EXSM = 1.0D0/L
                        RHSM = 1.0D0/(1.2D0*DSM**EXSM + 0.0000012D0)
                        RHDN = 0.0D0
                        IF (NQ .EQ. 1) GO TO 500
                           DDN = DVNRMS(N,YH(1,L),EWT)/TESCO(1,NQ)
                           EXDN = 1.0D0/NQ
                           RHDN = 1.0D0/(1.3D0*DDN**EXDN + 0.0000013D0)
  500                   CONTINUE
                        IF (RHSM .GE. RHUP) GO TO 550
                           IF (RHUP .LE. RHDN) GO TO 540
                              NEWQ = L
                              RH = RHUP
                              IF (RH .GE. 1.1D0) GO TO 520
                                 IALTH = 3
                                 R = 1.0D0/TESCO(2,NQU)
                                 DO 510 I = 1, N
                                    ACOR(I) = ACOR(I)*R
  510                            CONTINUE
C     ...........................EXIT
                                 GO TO 690
  520                         CONTINUE
                              R = EL(L)/L
                              DO 530 I = 1, N
                                 YH(I,NEWQ+1) = ACOR(I)*R
  530                         CONTINUE
                              NQ = NEWQ
                              L = NQ + 1
                              IRET = 2
C           ..................EXIT
                              GO TO 680
  540                      CONTINUE
                        GO TO 580
  550                   CONTINUE
                        IF (RHSM .LT. RHDN) GO TO 580
                           NEWQ = NQ
                           RH = RHSM
                           IF (KFLAG .EQ. 0 .AND. RH .LT. 1.1D0)
     1                        GO TO 560
                              IF (KFLAG .LE. -2) RH = MIN(RH,0.2D0)
C                             ------------------------------------------
C                              IF THERE IS A CHANGE OF ORDER, RESET NQ,
C                              L, AND THE COEFFICIENTS.  IN ANY CASE H
C                              IS RESET ACCORDING TO RH AND THE YH ARRAY
C                              IS RESCALED.  THEN EXIT FROM 680 IF THE
C                              STEP WAS OK, OR REDO THE STEP OTHERWISE.
C                             ------------------------------------------
C                 ............EXIT
                              IF (NEWQ .EQ. NQ) GO TO 650
                              NQ = NEWQ
                              L = NQ + 1
                              IRET = 2
C           ..................EXIT
                              GO TO 680
  560                      CONTINUE
                           IALTH = 3
                           R = 1.0D0/TESCO(2,NQU)
                           DO 570 I = 1, N
                              ACOR(I) = ACOR(I)*R
  570                      CONTINUE
C     .....................EXIT
                           GO TO 690
  580                   CONTINUE
                        NEWQ = NQ - 1
                        RH = RHDN
                        IF (KFLAG .LT. 0 .AND. RH .GT. 1.0D0) RH = 1.0D0
                        IF (KFLAG .EQ. 0 .AND. RH .LT. 1.1D0) GO TO 590
                           IF (KFLAG .LE. -2) RH = MIN(RH,0.2D0)
C                          ---------------------------------------------
C                           IF THERE IS A CHANGE OF ORDER, RESET NQ, L,
C                           AND THE COEFFICIENTS.  IN ANY CASE H IS
C                           RESET ACCORDING TO RH AND THE YH ARRAY IS
C                           RESCALED.  THEN EXIT FROM 680 IF THE STEP
C                           WAS OK, OR REDO THE STEP OTHERWISE.
C                          ---------------------------------------------
C                 .........EXIT
                           IF (NEWQ .EQ. NQ) GO TO 650
                           NQ = NEWQ
                           L = NQ + 1
                           IRET = 2
C           ...............EXIT
                           GO TO 680
  590                   CONTINUE
                        IALTH = 3
                        R = 1.0D0/TESCO(2,NQU)
                        DO 600 I = 1, N
                           ACOR(I) = ACOR(I)*R
  600                   CONTINUE
C     ..................EXIT
                        GO TO 690
  610                CONTINUE
C                    ---------------------------------------------------
C                     CONTROL REACHES THIS SECTION IF 3 OR MORE FAILURES
C                     HAVE OCCURRED.  IF 10 FAILURES HAVE OCCURRED, EXIT
C                     WITH KFLAG = -1.  IT IS ASSUMED THAT THE
C                     DERIVATIVES THAT HAVE ACCUMULATED IN THE YH ARRAY
C                     HAVE ERRORS OF THE WRONG ORDER.  HENCE THE FIRST
C                     DERIVATIVE IS RECOMPUTED, AND THE ORDER IS SET TO
C                     1.  THEN H IS REDUCED BY A FACTOR OF 10, AND THE
C                     STEP IS RETRIED, UNTIL IT SUCCEEDS OR H REACHES
C                     HMIN.
C                    ---------------------------------------------------
                     IF (KFLAG .NE. -10) GO TO 620
C                       ------------------------------------------------
C                        ALL RETURNS ARE MADE THROUGH THIS SECTION.  H
C                        IS SAVED IN HOLD TO ALLOW THE CALLER TO CHANGE
C                        H ON THE NEXT STEP.
C                       ------------------------------------------------
                        KFLAG = -1
C     ..................EXIT
                        GO TO 690
  620                CONTINUE
                     RH = 0.1D0
                     RH = MAX(HMIN/ABS(H),RH)
                     H = H*RH
                     DO 630 I = 1, N
                        Y(I) = YH(I,1)
  630                CONTINUE
                     CALL DF(TN,Y,SAVF,RPAR,IPAR)
                     NFE = NFE + 1
                     DO 640 I = 1, N
                        YH(I,2) = H*SAVF(I)
  640                CONTINUE
                     IPUP = MITER
                     IALTH = 5
C              ......EXIT
                     IF (NQ .NE. 1) GO TO 670
                  GO TO 170
  650             CONTINUE
  660             CONTINUE
                  RH = MAX(RH,HMIN/ABS(H))
               GO TO 110
  670          CONTINUE
               NQ = 1
               L = 2
               IRET = 3
  680       CONTINUE
         GO TO 70
  690 CONTINUE
      HOLD = H
      JSTART = 1
      RETURN
C     ----------------------- END OF SUBROUTINE DSTOD
C     -----------------------
      END