share/odepack/fortran/dstoka.f

   1 *DECK DSTOKA
   2       SUBROUTINE DSTOKA (NEQ, Y, YH, NYH, YH1, EWT, SAVF, SAVX, ACOR,
   3      1   WM, IWM, F, JAC, PSOL)
   4       EXTERNAL F, JAC, PSOL
   5       INTEGER NEQ, NYH, IWM
   6       DOUBLE PRECISION Y, YH, YH1, EWT, SAVF, SAVX, ACOR, WM
   7       DIMENSION NEQ(*), Y(*), YH(NYH,*), YH1(*), EWT(*), SAVF(*),
   8      1   SAVX(*), ACOR(*), WM(*), IWM(*)
   9       INTEGER IOWND, IALTH, IPUP, LMAX, MEO, NQNYH, NSLP,
  10      1   ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L,
  11      2   LYH, LEWT, LACOR, LSAVF, LWM, LIWM, METH, MITER,
  12      3   MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU
  13       INTEGER NEWT, NSFI, NSLJ, NJEV
  14       INTEGER JPRE, JACFLG, LOCWP, LOCIWP, LSAVX, KMP, MAXL, MNEWT,
  15      1   NNI, NLI, NPS, NCFN, NCFL
  16       DOUBLE PRECISION CONIT, CRATE, EL, ELCO, HOLD, RMAX, TESCO,
  17      2   CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND
  18       DOUBLE PRECISION STIFR
  19       DOUBLE PRECISION DELT, EPCON, SQRTN, RSQRTN
  20       COMMON /DLS001/ CONIT, CRATE, EL(13), ELCO(13,12),
  21      1   HOLD, RMAX, TESCO(3,12),
  22      2   CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND,
  23      3   IOWND(6), IALTH, IPUP, LMAX, MEO, NQNYH, NSLP,
  24      3   ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L,
  25      4   LYH, LEWT, LACOR, LSAVF, LWM, LIWM, METH, MITER,
  26      5   MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU
  27       COMMON /DLS002/ STIFR, NEWT, NSFI, NSLJ, NJEV
  28       COMMON /DLPK01/ DELT, EPCON, SQRTN, RSQRTN,
  29      1   JPRE, JACFLG, LOCWP, LOCIWP, LSAVX, KMP, MAXL, MNEWT,
  30      2   NNI, NLI, NPS, NCFN, NCFL
  31 C-----------------------------------------------------------------------
  32 C DSTOKA performs one step of the integration of an initial value
  33 C problem for a system of Ordinary Differential Equations.
  34 C
  35 C This routine was derived from Subroutine DSTODPK in the DLSODPK
  36 C package by the addition of automatic functional/Newton iteration
  37 C switching and logic for re-use of Jacobian data.
  38 C-----------------------------------------------------------------------
  39 C Note: DSTOKA is independent of the value of the iteration method
  40 C indicator MITER, when this is .ne. 0, and hence is independent
  41 C of the type of chord method used, or the Jacobian structure.
  42 C Communication with DSTOKA is done with the following variables:
  43 C
  44 C NEQ    = integer array containing problem size in NEQ(1), and
  45 C          passed as the NEQ argument in all calls to F and JAC.
  46 C Y      = an array of length .ge. N used as the Y argument in
  47 C          all calls to F and JAC.
  48 C YH     = an NYH by LMAX array containing the dependent variables
  49 C          and their approximate scaled derivatives, where
  50 C          LMAX = MAXORD + 1.  YH(i,j+1) contains the approximate
  51 C          j-th derivative of y(i), scaled by H**j/factorial(j)
  52 C          (j = 0,1,...,NQ).  On entry for the first step, the first
  53 C          two columns of YH must be set from the initial values.
  54 C NYH    = a constant integer .ge. N, the first dimension of YH.
  55 C YH1    = a one-dimensional array occupying the same space as YH.
  56 C EWT    = an array of length N containing multiplicative weights
  57 C          for local error measurements.  Local errors in y(i) are
  58 C          compared to 1.0/EWT(i) in various error tests.
  59 C SAVF   = an array of working storage, of length N.
  60 C          Also used for input of YH(*,MAXORD+2) when JSTART = -1
  61 C          and MAXORD .lt. the current order NQ.
  62 C SAVX   = an array of working storage, of length N.
  63 C ACOR   = a work array of length N, used for the accumulated
  64 C          corrections.  On a successful return, ACOR(i) contains
  65 C          the estimated one-step local error in y(i).
  66 C WM,IWM = real and integer work arrays associated with matrix
  67 C          operations in chord iteration (MITER .ne. 0).
  68 C CCMAX  = maximum relative change in H*EL0 before DSETPK is called.
  69 C H      = the step size to be attempted on the next step.
  70 C          H is altered by the error control algorithm during the
  71 C          problem.  H can be either positive or negative, but its
  72 C          sign must remain constant throughout the problem.
  73 C HMIN   = the minimum absolute value of the step size H to be used.
  74 C HMXI   = inverse of the maximum absolute value of H to be used.
  75 C          HMXI = 0.0 is allowed and corresponds to an infinite HMAX.
  76 C          HMIN and HMXI may be changed at any time, but will not
  77 C          take effect until the next change of H is considered.
  78 C TN     = the independent variable. TN is updated on each step taken.
  79 C JSTART = an integer used for input only, with the following
  80 C          values and meanings:
  81 C               0  perform the first step.
  82 C           .gt.0  take a new step continuing from the last.
  83 C              -1  take the next step with a new value of H, MAXORD,
  84 C                    N, METH, MITER, and/or matrix parameters.
  85 C              -2  take the next step with a new value of H,
  86 C                    but with other inputs unchanged.
  87 C          On return, JSTART is set to 1 to facilitate continuation.
  88 C KFLAG  = a completion code with the following meanings:
  89 C               0  the step was succesful.
  90 C              -1  the requested error could not be achieved.
  91 C              -2  corrector convergence could not be achieved.
  92 C              -3  fatal error in DSETPK or DSOLPK.
  93 C          A return with KFLAG = -1 or -2 means either
  94 C          ABS(H) = HMIN or 10 consecutive failures occurred.
  95 C          On a return with KFLAG negative, the values of TN and
  96 C          the YH array are as of the beginning of the last
  97 C          step, and H is the last step size attempted.
  98 C MAXORD = the maximum order of integration method to be allowed.
  99 C MAXCOR = the maximum number of corrector iterations allowed.
 100 C MSBP   = maximum number of steps between DSETPK calls (MITER .gt. 0).
 101 C MXNCF  = maximum number of convergence failures allowed.
 102 C METH/MITER = the method flags.  See description in driver.
 103 C N      = the number of first-order differential equations.
 104 C-----------------------------------------------------------------------
 105       INTEGER I, I1, IREDO, IRET, J, JB, JOK, M, NCF, NEWQ, NSLOW
 106       DOUBLE PRECISION DCON, DDN, DEL, DELP, DRC, DSM, DUP, EXDN, EXSM,
 107      1   EXUP, DFNORM, R, RH, RHDN, RHSM, RHUP, ROC, STIFF, TOLD, DVNORM
 108 C
 109       KFLAG = 0
 110       TOLD = TN
 111       NCF = 0
 112       IERPJ = 0
 113       IERSL = 0
 114       JCUR = 0
 115       ICF = 0
 116       DELP = 0.0D0
 117       IF (JSTART .GT. 0) GO TO 200
 118       IF (JSTART .EQ. -1) GO TO 100
 119       IF (JSTART .EQ. -2) GO TO 160
 120 C-----------------------------------------------------------------------
 121 C On the first call, the order is set to 1, and other variables are
 122 C initialized.  RMAX is the maximum ratio by which H can be increased
 123 C in a single step.  It is initially 1.E4 to compensate for the small
 124 C initial H, but then is normally equal to 10.  If a failure
 125 C occurs (in corrector convergence or error test), RMAX is set at 2
 126 C for the next increase.
 127 C-----------------------------------------------------------------------
 128       LMAX = MAXORD + 1
 129       NQ = 1
 130       L = 2
 131       IALTH = 2
 132       RMAX = 10000.0D0
 133       RC = 0.0D0
 134       EL0 = 1.0D0
 135       CRATE = 0.7D0
 136       HOLD = H
 137       MEO = METH
 138       NSLP = 0
 139       NSLJ = 0
 140       IPUP = 0
 141       IRET = 3
 142       NEWT = 0
 143       STIFR = 0.0D0
 144       GO TO 140
 145 C-----------------------------------------------------------------------
 146 C The following block handles preliminaries needed when JSTART = -1.
 147 C IPUP is set to MITER to force a matrix update.
 148 C If an order increase is about to be considered (IALTH = 1),
 149 C IALTH is reset to 2 to postpone consideration one more step.
 150 C If the caller has changed METH, DCFODE is called to reset
 151 C the coefficients of the method.
 152 C If the caller has changed MAXORD to a value less than the current
 153 C order NQ, NQ is reduced to MAXORD, and a new H chosen accordingly.
 154 C If H is to be changed, YH must be rescaled.
 155 C If H or METH is being changed, IALTH is reset to L = NQ + 1
 156 C to prevent further changes in H for that many steps.
 157 C-----------------------------------------------------------------------
 158  100  IPUP = MITER
 159       LMAX = MAXORD + 1
 160       IF (IALTH .EQ. 1) IALTH = 2
 161       IF (METH .EQ. MEO) GO TO 110
 162       CALL DCFODE (METH, ELCO, TESCO)
 163       MEO = METH
 164       IF (NQ .GT. MAXORD) GO TO 120
 165       IALTH = L
 166       IRET = 1
 167       GO TO 150
 168  110  IF (NQ .LE. MAXORD) GO TO 160
 169  120  NQ = MAXORD
 170       L = LMAX
 171       DO 125 I = 1,L
 172  125    EL(I) = ELCO(I,NQ)
 173       NQNYH = NQ*NYH
 174       RC = RC*EL(1)/EL0
 175       EL0 = EL(1)
 176       CONIT = 0.5D0/(NQ+2)
 177       EPCON = CONIT*TESCO(2,NQ)
 178       DDN = DVNORM (N, SAVF, EWT)/TESCO(1,L)
 179       EXDN = 1.0D0/L
 180       RHDN = 1.0D0/(1.3D0*DDN**EXDN + 0.0000013D0)
 181       RH = MIN(RHDN,1.0D0)
 182       IREDO = 3
 183       IF (H .EQ. HOLD) GO TO 170
 184       RH = MIN(RH,ABS(H/HOLD))
 185       H = HOLD
 186       GO TO 175
 187 C-----------------------------------------------------------------------
 188 C DCFODE is called to get all the integration coefficients for the
 189 C current METH.  Then the EL vector and related constants are reset
 190 C whenever the order NQ is changed, or at the start of the problem.
 191 C-----------------------------------------------------------------------
 192  140  CALL DCFODE (METH, ELCO, TESCO)
 193  150  DO 155 I = 1,L
 194  155    EL(I) = ELCO(I,NQ)
 195       NQNYH = NQ*NYH
 196       RC = RC*EL(1)/EL0
 197       EL0 = EL(1)
 198       CONIT = 0.5D0/(NQ+2)
 199       EPCON = CONIT*TESCO(2,NQ)
 200       GO TO (160, 170, 200), IRET
 201 C-----------------------------------------------------------------------
 202 C If H is being changed, the H ratio RH is checked against
 203 C RMAX, HMIN, and HMXI, and the YH array rescaled.  IALTH is set to
 204 C L = NQ + 1 to prevent a change of H for that many steps, unless
 205 C forced by a convergence or error test failure.
 206 C-----------------------------------------------------------------------
 207  160  IF (H .EQ. HOLD) GO TO 200
 208       RH = H/HOLD
 209       H = HOLD
 210       IREDO = 3
 211       GO TO 175
 212  170  RH = MAX(RH,HMIN/ABS(H))
 213  175  RH = MIN(RH,RMAX)
 214       RH = RH/MAX(1.0D0,ABS(H)*HMXI*RH)
 215       R = 1.0D0
 216       DO 180 J = 2,L
 217         R = R*RH
 218         DO 180 I = 1,N
 219  180      YH(I,J) = YH(I,J)*R
 220       H = H*RH
 221       RC = RC*RH
 222       IALTH = L
 223       IF (IREDO .EQ. 0) GO TO 690
 224 C-----------------------------------------------------------------------
 225 C This section computes the predicted values by effectively
 226 C multiplying the YH array by the Pascal triangle matrix.
 227 C The flag IPUP is set according to whether matrix data is involved
 228 C (NEWT .gt. 0 .and. JACFLG .ne. 0) or not, to trigger a call to DSETPK.
 229 C IPUP is set to MITER when RC differs from 1 by more than CCMAX,
 230 C and at least every MSBP steps, when JACFLG = 1.
 231 C RC is the ratio of new to old values of the coefficient  H*EL(1).
 232 C-----------------------------------------------------------------------
 233  200  IF (NEWT .EQ. 0 .OR. JACFLG .EQ. 0) THEN
 234         DRC = 0.0D0
 235         IPUP = 0
 236         CRATE = 0.7D0
 237       ELSE
 238         DRC = ABS(RC - 1.0D0)
 239         IF (DRC .GT. CCMAX) IPUP = MITER
 240         IF (NST .GE. NSLP+MSBP) IPUP = MITER
 241         ENDIF
 242       TN = TN + H
 243       I1 = NQNYH + 1
 244       DO 215 JB = 1,NQ
 245         I1 = I1 - NYH
 246 CDIR$ IVDEP
 247         DO 210 I = I1,NQNYH
 248  210      YH1(I) = YH1(I) + YH1(I+NYH)
 249  215    CONTINUE
 250 C-----------------------------------------------------------------------
 251 C Up to MAXCOR corrector iterations are taken.  A convergence test is
 252 C made on the RMS-norm of each correction, weighted by the error
 253 C weight vector EWT.  The sum of the corrections is accumulated in the
 254 C vector ACOR(i).  The YH array is not altered in the corrector loop.
 255 C Within the corrector loop, an estimated rate of convergence (ROC)
 256 C and a stiffness ratio estimate (STIFF) are kept.  Corresponding
 257 C global estimates are kept as CRATE and stifr.
 258 C-----------------------------------------------------------------------
 259  220  M = 0
 260       MNEWT = 0
 261       STIFF = 0.0D0
 262       ROC = 0.05D0
 263       NSLOW = 0
 264       DO 230 I = 1,N
 265  230    Y(I) = YH(I,1)
 266       CALL F (NEQ, TN, Y, SAVF)
 267       NFE = NFE + 1
 268       IF (NEWT .EQ. 0 .OR. IPUP .LE. 0) GO TO 250
 269 C-----------------------------------------------------------------------
 270 C If indicated, DSETPK is called to update any matrix data needed,
 271 C before starting the corrector iteration.
 272 C JOK is set to indicate if the matrix data need not be recomputed.
 273 C IPUP is set to 0 as an indicator that the matrix data is up to date.
 274 C-----------------------------------------------------------------------
 275       JOK = 1
 276       IF (NST .EQ. 0 .OR. NST .GT. NSLJ+50) JOK = -1
 277       IF (ICF .EQ. 1 .AND. DRC .LT. 0.2D0) JOK = -1
 278       IF (ICF .EQ. 2) JOK = -1
 279       IF (JOK .EQ. -1) THEN
 280         NSLJ = NST
 281         NJEV = NJEV + 1
 282         ENDIF
 283       CALL DSETPK (NEQ, Y, YH1, EWT, ACOR, SAVF, JOK, WM, IWM, F, JAC)
 284       IPUP = 0
 285       RC = 1.0D0
 286       DRC = 0.0D0
 287       NSLP = NST
 288       CRATE = 0.7D0
 289       IF (IERPJ .NE. 0) GO TO 430
 290  250  DO 260 I = 1,N
 291  260    ACOR(I) = 0.0D0
 292  270  IF (NEWT .NE. 0) GO TO 350
 293 C-----------------------------------------------------------------------
 294 C In the case of functional iteration, update Y directly from
 295 C the result of the last function evaluation, and STIFF is set to 1.0.
 296 C-----------------------------------------------------------------------
 297       DO 290 I = 1,N
 298         SAVF(I) = H*SAVF(I) - YH(I,2)
 299  290    Y(I) = SAVF(I) - ACOR(I)
 300       DEL = DVNORM (N, Y, EWT)
 301       DO 300 I = 1,N
 302         Y(I) = YH(I,1) + EL(1)*SAVF(I)
 303  300    ACOR(I) = SAVF(I)
 304       STIFF = 1.0D0
 305       GO TO 400
 306 C-----------------------------------------------------------------------
 307 C In the case of the chord method, compute the corrector error,
 308 C and solve the linear system with that as right-hand side and
 309 C P as coefficient matrix.  STIFF is set to the ratio of the norms
 310 C of the residual and the correction vector.
 311 C-----------------------------------------------------------------------
 312  350  DO 360 I = 1,N
 313  360    SAVX(I) = H*SAVF(I) - (YH(I,2) + ACOR(I))
 314       DFNORM = DVNORM (N, SAVX, EWT)
 315       CALL DSOLPK (NEQ, Y, SAVF, SAVX, EWT, WM, IWM, F, PSOL)
 316       IF (IERSL .LT. 0) GO TO 430
 317       IF (IERSL .GT. 0) GO TO 410
 318       DEL = DVNORM (N, SAVX, EWT)
 319       IF (DEL .GT. 1.0D-8) STIFF = MAX(STIFF, DFNORM/DEL)
 320       DO 380 I = 1,N
 321         ACOR(I) = ACOR(I) + SAVX(I)
 322  380    Y(I) = YH(I,1) + EL(1)*ACOR(I)
 323 C-----------------------------------------------------------------------
 324 C Test for convergence.  If M .gt. 0, an estimate of the convergence
 325 C rate constant is made for the iteration switch, and is also used
 326 C in the convergence test.   If the iteration seems to be diverging or
 327 C converging at a slow rate (.gt. 0.8 more than once), it is stopped.
 328 C-----------------------------------------------------------------------
 329  400  IF (M .NE. 0) THEN
 330         ROC = MAX(0.05D0, DEL/DELP)
 331         CRATE = MAX(0.2D0*CRATE,ROC)
 332         ENDIF
 333       DCON = DEL*MIN(1.0D0,1.5D0*CRATE)/EPCON
 334       IF (DCON .LE. 1.0D0) GO TO 450
 335       M = M + 1
 336       IF (M .EQ. MAXCOR) GO TO 410
 337       IF (M .GE. 2 .AND. DEL .GT. 2.0D0*DELP) GO TO 410
 338       IF (ROC .GT. 10.0D0) GO TO 410
 339       IF (ROC .GT. 0.8D0) NSLOW = NSLOW + 1
 340       IF (NSLOW .GE. 2) GO TO 410
 341       MNEWT = M
 342       DELP = DEL
 343       CALL F (NEQ, TN, Y, SAVF)
 344       NFE = NFE + 1
 345       GO TO 270
 346 C-----------------------------------------------------------------------
 347 C The corrector iteration failed to converge.
 348 C If functional iteration is being done (NEWT = 0) and MITER .gt. 0
 349 C (and this is not the first step), then switch to Newton
 350 C (NEWT = MITER), and retry the step.  (Setting STIFR = 1023 insures
 351 C that a switch back will not occur for 10 step attempts.)
 352 C If Newton iteration is being done, but using a preconditioner that
 353 C is out of date (JACFLG .ne. 0 .and. JCUR = 0), then signal for a
 354 C re-evalutation of the preconditioner, and retry the step.
 355 C In all other cases, the YH array is retracted to its values
 356 C before prediction, and H is reduced, if possible.  If H cannot be
 357 C reduced or MXNCF failures have occurred, exit with KFLAG = -2.
 358 C-----------------------------------------------------------------------
 359  410  ICF = 1
 360       IF (NEWT .EQ. 0) THEN
 361         IF (NST .EQ. 0) GO TO 430
 362         IF (MITER .EQ. 0) GO TO 430
 363         NEWT = MITER
 364         STIFR = 1023.0D0
 365         IPUP = MITER
 366         GO TO 220
 367         ENDIF
 368       IF (JCUR.EQ.1 .OR. JACFLG.EQ.0) GO TO 430
 369       IPUP = MITER
 370       GO TO 220
 371  430  ICF = 2
 372       NCF = NCF + 1
 373       NCFN = NCFN + 1
 374       RMAX = 2.0D0
 375       TN = TOLD
 376       I1 = NQNYH + 1
 377       DO 445 JB = 1,NQ
 378         I1 = I1 - NYH
 379 CDIR$ IVDEP
 380         DO 440 I = I1,NQNYH
 381  440      YH1(I) = YH1(I) - YH1(I+NYH)
 382  445    CONTINUE
 383       IF (IERPJ .LT. 0 .OR. IERSL .LT. 0) GO TO 680
 384       IF (ABS(H) .LE. HMIN*1.00001D0) GO TO 670
 385       IF (NCF .EQ. MXNCF) GO TO 670
 386       RH = 0.5D0
 387       IPUP = MITER
 388       IREDO = 1
 389       GO TO 170
 390 C-----------------------------------------------------------------------
 391 C The corrector has converged.  JCUR is set to 0 to signal that the
 392 C preconditioner involved may need updating later.
 393 C The stiffness ratio STIFR is updated using the latest STIFF value.
 394 C The local error test is made and control passes to statement 500
 395 C if it fails.
 396 C-----------------------------------------------------------------------
 397  450  JCUR = 0
 398       IF (NEWT .GT. 0) STIFR = 0.5D0*(STIFR + STIFF)
 399       IF (M .EQ. 0) DSM = DEL/TESCO(2,NQ)
 400       IF (M .GT. 0) DSM = DVNORM (N, ACOR, EWT)/TESCO(2,NQ)
 401       IF (DSM .GT. 1.0D0) GO TO 500
 402 C-----------------------------------------------------------------------
 403 C After a successful step, update the YH array.
 404 C If Newton iteration is being done and STIFR is less than 1.5,
 405 C then switch to functional iteration.
 406 C Consider changing H if IALTH = 1.  Otherwise decrease IALTH by 1.
 407 C If IALTH is then 1 and NQ .lt. MAXORD, then ACOR is saved for
 408 C use in a possible order increase on the next step.
 409 C If a change in H is considered, an increase or decrease in order
 410 C by one is considered also.  A change in H is made only if it is by a
 411 C factor of at least 1.1.  If not, IALTH is set to 3 to prevent
 412 C testing for that many steps.
 413 C-----------------------------------------------------------------------
 414       KFLAG = 0
 415       IREDO = 0
 416       NST = NST + 1
 417       IF (NEWT .EQ. 0) NSFI = NSFI + 1
 418       IF (NEWT .GT. 0 .AND. STIFR .LT. 1.5D0) NEWT = 0
 419       HU = H
 420       NQU = NQ
 421       DO 470 J = 1,L
 422         DO 470 I = 1,N
 423  470      YH(I,J) = YH(I,J) + EL(J)*ACOR(I)
 424       IALTH = IALTH - 1
 425       IF (IALTH .EQ. 0) GO TO 520
 426       IF (IALTH .GT. 1) GO TO 700
 427       IF (L .EQ. LMAX) GO TO 700
 428       DO 490 I = 1,N
 429  490    YH(I,LMAX) = ACOR(I)
 430       GO TO 700
 431 C-----------------------------------------------------------------------
 432 C The error test failed.  KFLAG keeps track of multiple failures.
 433 C Restore TN and the YH array to their previous values, and prepare
 434 C to try the step again.  Compute the optimum step size for this or
 435 C one lower order.  After 2 or more failures, H is forced to decrease
 436 C by a factor of 0.2 or less.
 437 C-----------------------------------------------------------------------
 438  500  KFLAG = KFLAG - 1
 439       TN = TOLD
 440       I1 = NQNYH + 1
 441       DO 515 JB = 1,NQ
 442         I1 = I1 - NYH
 443 CDIR$ IVDEP
 444         DO 510 I = I1,NQNYH
 445  510      YH1(I) = YH1(I) - YH1(I+NYH)
 446  515    CONTINUE
 447       RMAX = 2.0D0
 448       IF (ABS(H) .LE. HMIN*1.00001D0) GO TO 660
 449       IF (KFLAG .LE. -3) GO TO 640
 450       IREDO = 2
 451       RHUP = 0.0D0
 452       GO TO 540
 453 C-----------------------------------------------------------------------
 454 C Regardless of the success or failure of the step, factors
 455 C RHDN, RHSM, and RHUP are computed, by which H could be multiplied
 456 C at order NQ - 1, order NQ, or order NQ + 1, respectively.
 457 C in the case of failure, RHUP = 0.0 to avoid an order increase.
 458 C the largest of these is determined and the new order chosen
 459 C accordingly.  If the order is to be increased, we compute one
 460 C additional scaled derivative.
 461 C-----------------------------------------------------------------------
 462  520  RHUP = 0.0D0
 463       IF (L .EQ. LMAX) GO TO 540
 464       DO 530 I = 1,N
 465  530    SAVF(I) = ACOR(I) - YH(I,LMAX)
 466       DUP = DVNORM (N, SAVF, EWT)/TESCO(3,NQ)
 467       EXUP = 1.0D0/(L+1)
 468       RHUP = 1.0D0/(1.4D0*DUP**EXUP + 0.0000014D0)
 469  540  EXSM = 1.0D0/L
 470       RHSM = 1.0D0/(1.2D0*DSM**EXSM + 0.0000012D0)
 471       RHDN = 0.0D0
 472       IF (NQ .EQ. 1) GO TO 560
 473       DDN = DVNORM (N, YH(1,L), EWT)/TESCO(1,NQ)
 474       EXDN = 1.0D0/NQ
 475       RHDN = 1.0D0/(1.3D0*DDN**EXDN + 0.0000013D0)
 476  560  IF (RHSM .GE. RHUP) GO TO 570
 477       IF (RHUP .GT. RHDN) GO TO 590
 478       GO TO 580
 479  570  IF (RHSM .LT. RHDN) GO TO 580
 480       NEWQ = NQ
 481       RH = RHSM
 482       GO TO 620
 483  580  NEWQ = NQ - 1
 484       RH = RHDN
 485       IF (KFLAG .LT. 0 .AND. RH .GT. 1.0D0) RH = 1.0D0
 486       GO TO 620
 487  590  NEWQ = L
 488       RH = RHUP
 489       IF (RH .LT. 1.1D0) GO TO 610
 490       R = EL(L)/L
 491       DO 600 I = 1,N
 492  600    YH(I,NEWQ+1) = ACOR(I)*R
 493       GO TO 630
 494  610  IALTH = 3
 495       GO TO 700
 496  620  IF ((KFLAG .EQ. 0) .AND. (RH .LT. 1.1D0)) GO TO 610
 497       IF (KFLAG .LE. -2) RH = MIN(RH,0.2D0)
 498 C-----------------------------------------------------------------------
 499 C If there is a change of order, reset NQ, L, and the coefficients.
 500 C In any case H is reset according to RH and the YH array is rescaled.
 501 C Then exit from 690 if the step was OK, or redo the step otherwise.
 502 C-----------------------------------------------------------------------
 503       IF (NEWQ .EQ. NQ) GO TO 170
 504  630  NQ = NEWQ
 505       L = NQ + 1
 506       IRET = 2
 507       GO TO 150
 508 C-----------------------------------------------------------------------
 509 C Control reaches this section if 3 or more failures have occured.
 510 C If 10 failures have occurred, exit with KFLAG = -1.
 511 C It is assumed that the derivatives that have accumulated in the
 512 C YH array have errors of the wrong order.  Hence the first
 513 C derivative is recomputed, and the order is set to 1.  Then
 514 C H is reduced by a factor of 10, and the step is retried,
 515 C until it succeeds or H reaches HMIN.
 516 C-----------------------------------------------------------------------
 517  640  IF (KFLAG .EQ. -10) GO TO 660
 518       RH = 0.1D0
 519       RH = MAX(HMIN/ABS(H),RH)
 520       H = H*RH
 521       DO 645 I = 1,N
 522  645    Y(I) = YH(I,1)
 523       CALL F (NEQ, TN, Y, SAVF)
 524       NFE = NFE + 1
 525       DO 650 I = 1,N
 526  650    YH(I,2) = H*SAVF(I)
 527       IPUP = MITER
 528       IALTH = 5
 529       IF (NQ .EQ. 1) GO TO 200
 530       NQ = 1
 531       L = 2
 532       IRET = 3
 533       GO TO 150
 534 C-----------------------------------------------------------------------
 535 C All returns are made through this section.  H is saved in HOLD
 536 C to allow the caller to change H on the next step.
 537 C-----------------------------------------------------------------------
 538  660  KFLAG = -1
 539       GO TO 720
 540  670  KFLAG = -2
 541       GO TO 720
 542  680  KFLAG = -3
 543       GO TO 720
 544  690  RMAX = 10.0D0
 545  700  R = 1.0D0/TESCO(2,NQU)
 546       DO 710 I = 1,N
 547  710    ACOR(I) = ACOR(I)*R
 548  720  HOLD = H
 549       JSTART = 1
 550       RETURN
 551 C----------------------- End of Subroutine DSTOKA ----------------------
 552       END