share/odepack/fortran/dstoda.f

   1 *DECK DSTODA
   2       SUBROUTINE DSTODA (NEQ, Y, YH, NYH, YH1, EWT, SAVF, ACOR,
   3      1   WM, IWM, F, JAC, PJAC, SLVS)
   4       EXTERNAL F, JAC, PJAC, SLVS
   5       INTEGER NEQ, NYH, IWM
   6       DOUBLE PRECISION Y, YH, YH1, EWT, SAVF, ACOR, WM
   7       DIMENSION NEQ(*), Y(*), YH(NYH,*), YH1(*), EWT(*), SAVF(*),
   8      1   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 IOWND2, ICOUNT, IRFLAG, JTYP, MUSED, MXORDN, MXORDS
  14       DOUBLE PRECISION CONIT, CRATE, EL, ELCO, HOLD, RMAX, TESCO,
  15      2   CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND
  16       DOUBLE PRECISION ROWND2, CM1, CM2, PDEST, PDLAST, RATIO,
  17      1   PDNORM
  18       COMMON /DLS001/ CONIT, CRATE, EL(13), ELCO(13,12),
  19      1   HOLD, RMAX, TESCO(3,12),
  20      2   CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND,
  21      3   IOWND(6), IALTH, IPUP, LMAX, MEO, NQNYH, NSLP,
  22      4   ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L,
  23      5   LYH, LEWT, LACOR, LSAVF, LWM, LIWM, METH, MITER,
  24      6   MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU
  25       COMMON /DLSA01/ ROWND2, CM1(12), CM2(5), PDEST, PDLAST, RATIO,
  26      1   PDNORM,
  27      2   IOWND2(3), ICOUNT, IRFLAG, JTYP, MUSED, MXORDN, MXORDS
  28       INTEGER I, I1, IREDO, IRET, J, JB, M, NCF, NEWQ
  29       INTEGER LM1, LM1P1, LM2, LM2P1, NQM1, NQM2
  30       DOUBLE PRECISION DCON, DDN, DEL, DELP, DSM, DUP, EXDN, EXSM, EXUP,
  31      1   R, RH, RHDN, RHSM, RHUP, TOLD, DMNORM
  32       DOUBLE PRECISION ALPHA, DM1,DM2, EXM1,EXM2,
  33      1   PDH, PNORM, RATE, RH1, RH1IT, RH2, RM, SM1(12)
  34       SAVE SM1
  35       DATA SM1/0.5D0, 0.575D0, 0.55D0, 0.45D0, 0.35D0, 0.25D0,
  36      1   0.20D0, 0.15D0, 0.10D0, 0.075D0, 0.050D0, 0.025D0/
  37 C-----------------------------------------------------------------------
  38 C DSTODA performs one step of the integration of an initial value
  39 C problem for a system of ordinary differential equations.
  40 C Note: DSTODA is independent of the value of the iteration method
  41 C indicator MITER, when this is .ne. 0, and hence is independent
  42 C of the type of chord method used, or the Jacobian structure.
  43 C Communication with DSTODA is done with the following variables:
  44 C
  45 C Y      = an array of length .ge. N used as the Y argument in
  46 C          all calls to F and JAC.
  47 C NEQ    = integer array containing problem size in NEQ(1), and
  48 C          passed as the NEQ argument in all calls to F and JAC.
  49 C YH     = an NYH by LMAX array containing the dependent variables
  50 C          and their approximate scaled derivatives, where
  51 C          LMAX = MAXORD + 1.  YH(i,j+1) contains the approximate
  52 C          j-th derivative of y(i), scaled by H**j/factorial(j)
  53 C          (j = 0,1,...,NQ).  On entry for the first step, the first
  54 C          two columns of YH must be set from the initial values.
  55 C NYH    = a constant integer .ge. N, the first dimension of YH.
  56 C YH1    = a one-dimensional array occupying the same space as YH.
  57 C EWT    = an array of length N containing multiplicative weights
  58 C          for local error measurements.  Local errors in y(i) are
  59 C          compared to 1.0/EWT(i) in various error tests.
  60 C SAVF   = an array of working storage, of length N.
  61 C ACOR   = a work array of length N, used for the accumulated
  62 C          corrections.  On a successful return, ACOR(i) contains
  63 C          the estimated one-step local error in y(i).
  64 C WM,IWM = real and integer work arrays associated with matrix
  65 C          operations in chord iteration (MITER .ne. 0).
  66 C PJAC   = name of routine to evaluate and preprocess Jacobian matrix
  67 C          and P = I - H*EL0*Jac, if a chord method is being used.
  68 C          It also returns an estimate of norm(Jac) in PDNORM.
  69 C SLVS   = name of routine to solve linear system in chord iteration.
  70 C CCMAX  = maximum relative change in H*EL0 before PJAC is called.
  71 C H      = the step size to be attempted on the next step.
  72 C          H is altered by the error control algorithm during the
  73 C          problem.  H can be either positive or negative, but its
  74 C          sign must remain constant throughout the problem.
  75 C HMIN   = the minimum absolute value of the step size H to be used.
  76 C HMXI   = inverse of the maximum absolute value of H to be used.
  77 C          HMXI = 0.0 is allowed and corresponds to an infinite HMAX.
  78 C          HMIN and HMXI may be changed at any time, but will not
  79 C          take effect until the next change of H is considered.
  80 C TN     = the independent variable. TN is updated on each step taken.
  81 C JSTART = an integer used for input only, with the following
  82 C          values and meanings:
  83 C               0  perform the first step.
  84 C           .gt.0  take a new step continuing from the last.
  85 C              -1  take the next step with a new value of H,
  86 C                    N, METH, MITER, and/or matrix parameters.
  87 C              -2  take the next step with a new value of H,
  88 C                    but with other inputs unchanged.
  89 C          On return, JSTART is set to 1 to facilitate continuation.
  90 C KFLAG  = a completion code with the following meanings:
  91 C               0  the step was succesful.
  92 C              -1  the requested error could not be achieved.
  93 C              -2  corrector convergence could not be achieved.
  94 C              -3  fatal error in PJAC or SLVS.
  95 C          A return with KFLAG = -1 or -2 means either
  96 C          ABS(H) = HMIN or 10 consecutive failures occurred.
  97 C          On a return with KFLAG negative, the values of TN and
  98 C          the YH array are as of the beginning of the last
  99 C          step, and H is the last step size attempted.
 100 C MAXORD = the maximum order of integration method to be allowed.
 101 C MAXCOR = the maximum number of corrector iterations allowed.
 102 C MSBP   = maximum number of steps between PJAC calls (MITER .gt. 0).
 103 C MXNCF  = maximum number of convergence failures allowed.
 104 C METH   = current method.
 105 C          METH = 1 means Adams method (nonstiff)
 106 C          METH = 2 means BDF method (stiff)
 107 C          METH may be reset by DSTODA.
 108 C MITER  = corrector iteration method.
 109 C          MITER = 0 means functional iteration.
 110 C          MITER = JT .gt. 0 means a chord iteration corresponding
 111 C          to Jacobian type JT.  (The DLSODA/DLSODAR argument JT is
 112 C          communicated here as JTYP, but is not used in DSTODA
 113 C          except to load MITER following a method switch.)
 114 C          MITER may be reset by DSTODA.
 115 C N      = the number of first-order differential equations.
 116 C-----------------------------------------------------------------------
 117       KFLAG = 0
 118       TOLD = TN
 119       NCF = 0
 120       IERPJ = 0
 121       IERSL = 0
 122       JCUR = 0
 123       ICF = 0
 124       DELP = 0.0D0
 125       IF (JSTART .GT. 0) GO TO 200
 126       IF (JSTART .EQ. -1) GO TO 100
 127       IF (JSTART .EQ. -2) GO TO 160
 128 C-----------------------------------------------------------------------
 129 C On the first call, the order is set to 1, and other variables are
 130 C initialized.  RMAX is the maximum ratio by which H can be increased
 131 C in a single step.  It is initially 1.E4 to compensate for the small
 132 C initial H, but then is normally equal to 10.  If a failure
 133 C occurs (in corrector convergence or error test), RMAX is set at 2
 134 C for the next increase.
 135 C DCFODE is called to get the needed coefficients for both methods.
 136 C-----------------------------------------------------------------------
 137       LMAX = MAXORD + 1
 138       NQ = 1
 139       L = 2
 140       IALTH = 2
 141       RMAX = 10000.0D0
 142       RC = 0.0D0
 143       EL0 = 1.0D0
 144       CRATE = 0.7D0
 145       HOLD = H
 146       NSLP = 0
 147       IPUP = MITER
 148       IRET = 3
 149 C Initialize switching parameters.  METH = 1 is assumed initially. -----
 150       ICOUNT = 20
 151       IRFLAG = 0
 152       PDEST = 0.0D0
 153       PDLAST = 0.0D0
 154       RATIO = 5.0D0
 155       CALL DCFODE (2, ELCO, TESCO)
 156       DO 10 I = 1,5
 157  10     CM2(I) = TESCO(2,I)*ELCO(I+1,I)
 158       CALL DCFODE (1, ELCO, TESCO)
 159       DO 20 I = 1,12
 160  20     CM1(I) = TESCO(2,I)*ELCO(I+1,I)
 161       GO TO 150
 162 C-----------------------------------------------------------------------
 163 C The following block handles preliminaries needed when JSTART = -1.
 164 C IPUP is set to MITER to force a matrix update.
 165 C If an order increase is about to be considered (IALTH = 1),
 166 C IALTH is reset to 2 to postpone consideration one more step.
 167 C If the caller has changed METH, DCFODE is called to reset
 168 C the coefficients of the method.
 169 C If H is to be changed, YH must be rescaled.
 170 C If H or METH is being changed, IALTH is reset to L = NQ + 1
 171 C to prevent further changes in H for that many steps.
 172 C-----------------------------------------------------------------------
 173  100  IPUP = MITER
 174       LMAX = MAXORD + 1
 175       IF (IALTH .EQ. 1) IALTH = 2
 176       IF (METH .EQ. MUSED) GO TO 160
 177       CALL DCFODE (METH, ELCO, TESCO)
 178       IALTH = L
 179       IRET = 1
 180 C-----------------------------------------------------------------------
 181 C The el vector and related constants are reset
 182 C whenever the order NQ is changed, or at the start of the problem.
 183 C-----------------------------------------------------------------------
 184  150  DO 155 I = 1,L
 185  155    EL(I) = ELCO(I,NQ)
 186       NQNYH = NQ*NYH
 187       RC = RC*EL(1)/EL0
 188       EL0 = EL(1)
 189       CONIT = 0.5D0/(NQ+2)
 190       GO TO (160, 170, 200), IRET
 191 C-----------------------------------------------------------------------
 192 C If H is being changed, the H ratio RH is checked against
 193 C RMAX, HMIN, and HMXI, and the YH array rescaled.  IALTH is set to
 194 C L = NQ + 1 to prevent a change of H for that many steps, unless
 195 C forced by a convergence or error test failure.
 196 C-----------------------------------------------------------------------
 197  160  IF (H .EQ. HOLD) GO TO 200
 198       RH = H/HOLD
 199       H = HOLD
 200       IREDO = 3
 201       GO TO 175
 202  170  RH = MAX(RH,HMIN/ABS(H))
 203  175  RH = MIN(RH,RMAX)
 204       RH = RH/MAX(1.0D0,ABS(H)*HMXI*RH)
 205 C-----------------------------------------------------------------------
 206 C If METH = 1, also restrict the new step size by the stability region.
 207 C If this reduces H, set IRFLAG to 1 so that if there are roundoff
 208 C problems later, we can assume that is the cause of the trouble.
 209 C-----------------------------------------------------------------------
 210       IF (METH .EQ. 2) GO TO 178
 211       IRFLAG = 0
 212       PDH = MAX(ABS(H)*PDLAST,0.000001D0)
 213       IF (RH*PDH*1.00001D0 .LT. SM1(NQ)) GO TO 178
 214       RH = SM1(NQ)/PDH
 215       IRFLAG = 1
 216  178  CONTINUE
 217       R = 1.0D0
 218       DO 180 J = 2,L
 219         R = R*RH
 220         DO 180 I = 1,N
 221  180      YH(I,J) = YH(I,J)*R
 222       H = H*RH
 223       RC = RC*RH
 224       IALTH = L
 225       IF (IREDO .EQ. 0) GO TO 690
 226 C-----------------------------------------------------------------------
 227 C This section computes the predicted values by effectively
 228 C multiplying the YH array by the Pascal triangle matrix.
 229 C RC is the ratio of new to old values of the coefficient  H*EL(1).
 230 C When RC differs from 1 by more than CCMAX, IPUP is set to MITER
 231 C to force PJAC to be called, if a Jacobian is involved.
 232 C In any case, PJAC is called at least every MSBP steps.
 233 C-----------------------------------------------------------------------
 234  200  IF (ABS(RC-1.0D0) .GT. CCMAX) IPUP = MITER
 235       IF (NST .GE. NSLP+MSBP) IPUP = MITER
 236       TN = TN + H
 237       I1 = NQNYH + 1
 238       DO 215 JB = 1,NQ
 239         I1 = I1 - NYH
 240 CDIR$ IVDEP
 241         DO 210 I = I1,NQNYH
 242  210      YH1(I) = YH1(I) + YH1(I+NYH)
 243  215    CONTINUE
 244       PNORM = DMNORM (N, YH1, EWT)
 245 C-----------------------------------------------------------------------
 246 C Up to MAXCOR corrector iterations are taken.  A convergence test is
 247 C made on the RMS-norm of each correction, weighted by the error
 248 C weight vector EWT.  The sum of the corrections is accumulated in the
 249 C vector ACOR(i).  The YH array is not altered in the corrector loop.
 250 C-----------------------------------------------------------------------
 251  220  M = 0
 252       RATE = 0.0D0
 253       DEL = 0.0D0
 254       DO 230 I = 1,N
 255  230    Y(I) = YH(I,1)
 256       CALL F (NEQ, TN, Y, SAVF)
 257       NFE = NFE + 1
 258       IF (IPUP .LE. 0) GO TO 250
 259 C-----------------------------------------------------------------------
 260 C If indicated, the matrix P = I - H*EL(1)*J is reevaluated and
 261 C preprocessed before starting the corrector iteration.  IPUP is set
 262 C to 0 as an indicator that this has been done.
 263 C-----------------------------------------------------------------------
 264       CALL PJAC (NEQ, Y, YH, NYH, EWT, ACOR, SAVF, WM, IWM, F, JAC)
 265       IPUP = 0
 266       RC = 1.0D0
 267       NSLP = NST
 268       CRATE = 0.7D0
 269       IF (IERPJ .NE. 0) GO TO 430
 270  250  DO 260 I = 1,N
 271  260    ACOR(I) = 0.0D0
 272  270  IF (MITER .NE. 0) GO TO 350
 273 C-----------------------------------------------------------------------
 274 C In the case of functional iteration, update Y directly from
 275 C the result of the last function evaluation.
 276 C-----------------------------------------------------------------------
 277       DO 290 I = 1,N
 278         SAVF(I) = H*SAVF(I) - YH(I,2)
 279  290    Y(I) = SAVF(I) - ACOR(I)
 280       DEL = DMNORM (N, Y, EWT)
 281       DO 300 I = 1,N
 282         Y(I) = YH(I,1) + EL(1)*SAVF(I)
 283  300    ACOR(I) = SAVF(I)
 284       GO TO 400
 285 C-----------------------------------------------------------------------
 286 C In the case of the chord method, compute the corrector error,
 287 C and solve the linear system with that as right-hand side and
 288 C P as coefficient matrix.
 289 C-----------------------------------------------------------------------
 290  350  DO 360 I = 1,N
 291  360    Y(I) = H*SAVF(I) - (YH(I,2) + ACOR(I))
 292       CALL SLVS (WM, IWM, Y, SAVF)
 293       IF (IERSL .LT. 0) GO TO 430
 294       IF (IERSL .GT. 0) GO TO 410
 295       DEL = DMNORM (N, Y, EWT)
 296       DO 380 I = 1,N
 297         ACOR(I) = ACOR(I) + Y(I)
 298  380    Y(I) = YH(I,1) + EL(1)*ACOR(I)
 299 C-----------------------------------------------------------------------
 300 C Test for convergence.  If M .gt. 0, an estimate of the convergence
 301 C rate constant is stored in CRATE, and this is used in the test.
 302 C
 303 C We first check for a change of iterates that is the size of
 304 C roundoff error.  If this occurs, the iteration has converged, and a
 305 C new rate estimate is not formed.
 306 C In all other cases, force at least two iterations to estimate a
 307 C local Lipschitz constant estimate for Adams methods.
 308 C On convergence, form PDEST = local maximum Lipschitz constant
 309 C estimate.  PDLAST is the most recent nonzero estimate.
 310 C-----------------------------------------------------------------------
 311  400  CONTINUE
 312       IF (DEL .LE. 100.0D0*PNORM*UROUND) GO TO 450
 313       IF (M .EQ. 0 .AND. METH .EQ. 1) GO TO 405
 314       IF (M .EQ. 0) GO TO 402
 315       RM = 1024.0D0
 316       IF (DEL .LE. 1024.0D0*DELP) RM = DEL/DELP
 317       RATE = MAX(RATE,RM)
 318       CRATE = MAX(0.2D0*CRATE,RM)
 319  402  DCON = DEL*MIN(1.0D0,1.5D0*CRATE)/(TESCO(2,NQ)*CONIT)
 320       IF (DCON .GT. 1.0D0) GO TO 405
 321       PDEST = MAX(PDEST,RATE/ABS(H*EL(1)))
 322       IF (PDEST .NE. 0.0D0) PDLAST = PDEST
 323       GO TO 450
 324  405  CONTINUE
 325       M = M + 1
 326       IF (M .EQ. MAXCOR) GO TO 410
 327       IF (M .GE. 2 .AND. DEL .GT. 2.0D0*DELP) GO TO 410
 328       DELP = DEL
 329       CALL F (NEQ, TN, Y, SAVF)
 330       NFE = NFE + 1
 331       GO TO 270
 332 C-----------------------------------------------------------------------
 333 C The corrector iteration failed to converge.
 334 C If MITER .ne. 0 and the Jacobian is out of date, PJAC is called for
 335 C the next try.  Otherwise the YH array is retracted to its values
 336 C before prediction, and H is reduced, if possible.  If H cannot be
 337 C reduced or MXNCF failures have occurred, exit with KFLAG = -2.
 338 C-----------------------------------------------------------------------
 339  410  IF (MITER .EQ. 0 .OR. JCUR .EQ. 1) GO TO 430
 340       ICF = 1
 341       IPUP = MITER
 342       GO TO 220
 343  430  ICF = 2
 344       NCF = NCF + 1
 345       RMAX = 2.0D0
 346       TN = TOLD
 347       I1 = NQNYH + 1
 348       DO 445 JB = 1,NQ
 349         I1 = I1 - NYH
 350 CDIR$ IVDEP
 351         DO 440 I = I1,NQNYH
 352  440      YH1(I) = YH1(I) - YH1(I+NYH)
 353  445    CONTINUE
 354       IF (IERPJ .LT. 0 .OR. IERSL .LT. 0) GO TO 680
 355       IF (ABS(H) .LE. HMIN*1.00001D0) GO TO 670
 356       IF (NCF .EQ. MXNCF) GO TO 670
 357       RH = 0.25D0
 358       IPUP = MITER
 359       IREDO = 1
 360       GO TO 170
 361 C-----------------------------------------------------------------------
 362 C The corrector has converged.  JCUR is set to 0
 363 C to signal that the Jacobian involved may need updating later.
 364 C The local error test is made and control passes to statement 500
 365 C if it fails.
 366 C-----------------------------------------------------------------------
 367  450  JCUR = 0
 368       IF (M .EQ. 0) DSM = DEL/TESCO(2,NQ)
 369       IF (M .GT. 0) DSM = DMNORM (N, ACOR, EWT)/TESCO(2,NQ)
 370       IF (DSM .GT. 1.0D0) GO TO 500
 371 C-----------------------------------------------------------------------
 372 C After a successful step, update the YH array.
 373 C Decrease ICOUNT by 1, and if it is -1, consider switching methods.
 374 C If a method switch is made, reset various parameters,
 375 C rescale the YH array, and exit.  If there is no switch,
 376 C consider changing H if IALTH = 1.  Otherwise decrease IALTH by 1.
 377 C If IALTH is then 1 and NQ .lt. MAXORD, then ACOR is saved for
 378 C use in a possible order increase on the next step.
 379 C If a change in H is considered, an increase or decrease in order
 380 C by one is considered also.  A change in H is made only if it is by a
 381 C factor of at least 1.1.  If not, IALTH is set to 3 to prevent
 382 C testing for that many steps.
 383 C-----------------------------------------------------------------------
 384       KFLAG = 0
 385       IREDO = 0
 386       NST = NST + 1
 387       HU = H
 388       NQU = NQ
 389       MUSED = METH
 390       DO 460 J = 1,L
 391         DO 460 I = 1,N
 392  460      YH(I,J) = YH(I,J) + EL(J)*ACOR(I)
 393       ICOUNT = ICOUNT - 1
 394       IF (ICOUNT .GE. 0) GO TO 488
 395       IF (METH .EQ. 2) GO TO 480
 396 C-----------------------------------------------------------------------
 397 C We are currently using an Adams method.  Consider switching to BDF.
 398 C If the current order is greater than 5, assume the problem is
 399 C not stiff, and skip this section.
 400 C If the Lipschitz constant and error estimate are not polluted
 401 C by roundoff, go to 470 and perform the usual test.
 402 C Otherwise, switch to the BDF methods if the last step was
 403 C restricted to insure stability (irflag = 1), and stay with Adams
 404 C method if not.  When switching to BDF with polluted error estimates,
 405 C in the absence of other information, double the step size.
 406 C
 407 C When the estimates are OK, we make the usual test by computing
 408 C the step size we could have (ideally) used on this step,
 409 C with the current (Adams) method, and also that for the BDF.
 410 C If NQ .gt. MXORDS, we consider changing to order MXORDS on switching.
 411 C Compare the two step sizes to decide whether to switch.
 412 C The step size advantage must be at least RATIO = 5 to switch.
 413 C-----------------------------------------------------------------------
 414       IF (NQ .GT. 5) GO TO 488
 415       IF (DSM .GT. 100.0D0*PNORM*UROUND .AND. PDEST .NE. 0.0D0)
 416      1   GO TO 470
 417       IF (IRFLAG .EQ. 0) GO TO 488
 418       RH2 = 2.0D0
 419       NQM2 = MIN(NQ,MXORDS)
 420       GO TO 478
 421  470  CONTINUE
 422       EXSM = 1.0D0/L
 423       RH1 = 1.0D0/(1.2D0*DSM**EXSM + 0.0000012D0)
 424       RH1IT = 2.0D0*RH1
 425       PDH = PDLAST*ABS(H)
 426       IF (PDH*RH1 .GT. 0.00001D0) RH1IT = SM1(NQ)/PDH
 427       RH1 = MIN(RH1,RH1IT)
 428       IF (NQ .LE. MXORDS) GO TO 474
 429          NQM2 = MXORDS
 430          LM2 = MXORDS + 1
 431          EXM2 = 1.0D0/LM2
 432          LM2P1 = LM2 + 1
 433          DM2 = DMNORM (N, YH(1,LM2P1), EWT)/CM2(MXORDS)
 434          RH2 = 1.0D0/(1.2D0*DM2**EXM2 + 0.0000012D0)
 435          GO TO 476
 436  474  DM2 = DSM*(CM1(NQ)/CM2(NQ))
 437       RH2 = 1.0D0/(1.2D0*DM2**EXSM + 0.0000012D0)
 438       NQM2 = NQ
 439  476  CONTINUE
 440       IF (RH2 .LT. RATIO*RH1) GO TO 488
 441 C THE SWITCH TEST PASSED.  RESET RELEVANT QUANTITIES FOR BDF. ----------
 442  478  RH = RH2
 443       ICOUNT = 20
 444       METH = 2
 445       MITER = JTYP
 446       PDLAST = 0.0D0
 447       NQ = NQM2
 448       L = NQ + 1
 449       GO TO 170
 450 C-----------------------------------------------------------------------
 451 C We are currently using a BDF method.  Consider switching to Adams.
 452 C Compute the step size we could have (ideally) used on this step,
 453 C with the current (BDF) method, and also that for the Adams.
 454 C If NQ .gt. MXORDN, we consider changing to order MXORDN on switching.
 455 C Compare the two step sizes to decide whether to switch.
 456 C The step size advantage must be at least 5/RATIO = 1 to switch.
 457 C If the step size for Adams would be so small as to cause
 458 C roundoff pollution, we stay with BDF.
 459 C-----------------------------------------------------------------------
 460  480  CONTINUE
 461       EXSM = 1.0D0/L
 462       IF (MXORDN .GE. NQ) GO TO 484
 463          NQM1 = MXORDN
 464          LM1 = MXORDN + 1
 465          EXM1 = 1.0D0/LM1
 466          LM1P1 = LM1 + 1
 467          DM1 = DMNORM (N, YH(1,LM1P1), EWT)/CM1(MXORDN)
 468          RH1 = 1.0D0/(1.2D0*DM1**EXM1 + 0.0000012D0)
 469          GO TO 486
 470  484  DM1 = DSM*(CM2(NQ)/CM1(NQ))
 471       RH1 = 1.0D0/(1.2D0*DM1**EXSM + 0.0000012D0)
 472       NQM1 = NQ
 473       EXM1 = EXSM
 474  486  RH1IT = 2.0D0*RH1
 475       PDH = PDNORM*ABS(H)
 476       IF (PDH*RH1 .GT. 0.00001D0) RH1IT = SM1(NQM1)/PDH
 477       RH1 = MIN(RH1,RH1IT)
 478       RH2 = 1.0D0/(1.2D0*DSM**EXSM + 0.0000012D0)
 479       IF (RH1*RATIO .LT. 5.0D0*RH2) GO TO 488
 480       ALPHA = MAX(0.001D0,RH1)
 481       DM1 = (ALPHA**EXM1)*DM1
 482       IF (DM1 .LE. 1000.0D0*UROUND*PNORM) GO TO 488
 483 C The switch test passed.  Reset relevant quantities for Adams. --------
 484       RH = RH1
 485       ICOUNT = 20
 486       METH = 1
 487       MITER = 0
 488       PDLAST = 0.0D0
 489       NQ = NQM1
 490       L = NQ + 1
 491       GO TO 170
 492 C
 493 C No method switch is being made.  Do the usual step/order selection. --
 494  488  CONTINUE
 495       IALTH = IALTH - 1
 496       IF (IALTH .EQ. 0) GO TO 520
 497       IF (IALTH .GT. 1) GO TO 700
 498       IF (L .EQ. LMAX) GO TO 700
 499       DO 490 I = 1,N
 500  490    YH(I,LMAX) = ACOR(I)
 501       GO TO 700
 502 C-----------------------------------------------------------------------
 503 C The error test failed.  KFLAG keeps track of multiple failures.
 504 C Restore TN and the YH array to their previous values, and prepare
 505 C to try the step again.  Compute the optimum step size for this or
 506 C one lower order.  After 2 or more failures, H is forced to decrease
 507 C by a factor of 0.2 or less.
 508 C-----------------------------------------------------------------------
 509  500  KFLAG = KFLAG - 1
 510       TN = TOLD
 511       I1 = NQNYH + 1
 512       DO 515 JB = 1,NQ
 513         I1 = I1 - NYH
 514 CDIR$ IVDEP
 515         DO 510 I = I1,NQNYH
 516  510      YH1(I) = YH1(I) - YH1(I+NYH)
 517  515    CONTINUE
 518       RMAX = 2.0D0
 519       IF (ABS(H) .LE. HMIN*1.00001D0) GO TO 660
 520       IF (KFLAG .LE. -3) GO TO 640
 521       IREDO = 2
 522       RHUP = 0.0D0
 523       GO TO 540
 524 C-----------------------------------------------------------------------
 525 C Regardless of the success or failure of the step, factors
 526 C RHDN, RHSM, and RHUP are computed, by which H could be multiplied
 527 C at order NQ - 1, order NQ, or order NQ + 1, respectively.
 528 C In the case of failure, RHUP = 0.0 to avoid an order increase.
 529 C The largest of these is determined and the new order chosen
 530 C accordingly.  If the order is to be increased, we compute one
 531 C additional scaled derivative.
 532 C-----------------------------------------------------------------------
 533  520  RHUP = 0.0D0
 534       IF (L .EQ. LMAX) GO TO 540
 535       DO 530 I = 1,N
 536  530    SAVF(I) = ACOR(I) - YH(I,LMAX)
 537       DUP = DMNORM (N, SAVF, EWT)/TESCO(3,NQ)
 538       EXUP = 1.0D0/(L+1)
 539       RHUP = 1.0D0/(1.4D0*DUP**EXUP + 0.0000014D0)
 540  540  EXSM = 1.0D0/L
 541       RHSM = 1.0D0/(1.2D0*DSM**EXSM + 0.0000012D0)
 542       RHDN = 0.0D0
 543       IF (NQ .EQ. 1) GO TO 550
 544       DDN = DMNORM (N, YH(1,L), EWT)/TESCO(1,NQ)
 545       EXDN = 1.0D0/NQ
 546       RHDN = 1.0D0/(1.3D0*DDN**EXDN + 0.0000013D0)
 547 C If METH = 1, limit RH according to the stability region also. --------
 548  550  IF (METH .EQ. 2) GO TO 560
 549       PDH = MAX(ABS(H)*PDLAST,0.000001D0)
 550       IF (L .LT. LMAX) RHUP = MIN(RHUP,SM1(L)/PDH)
 551       RHSM = MIN(RHSM,SM1(NQ)/PDH)
 552       IF (NQ .GT. 1) RHDN = MIN(RHDN,SM1(NQ-1)/PDH)
 553       PDEST = 0.0D0
 554  560  IF (RHSM .GE. RHUP) GO TO 570
 555       IF (RHUP .GT. RHDN) GO TO 590
 556       GO TO 580
 557  570  IF (RHSM .LT. RHDN) GO TO 580
 558       NEWQ = NQ
 559       RH = RHSM
 560       GO TO 620
 561  580  NEWQ = NQ - 1
 562       RH = RHDN
 563       IF (KFLAG .LT. 0 .AND. RH .GT. 1.0D0) RH = 1.0D0
 564       GO TO 620
 565  590  NEWQ = L
 566       RH = RHUP
 567       IF (RH .LT. 1.1D0) GO TO 610
 568       R = EL(L)/L
 569       DO 600 I = 1,N
 570  600    YH(I,NEWQ+1) = ACOR(I)*R
 571       GO TO 630
 572  610  IALTH = 3
 573       GO TO 700
 574 C If METH = 1 and H is restricted by stability, bypass 10 percent test.
 575  620  IF (METH .EQ. 2) GO TO 622
 576       IF (RH*PDH*1.00001D0 .GE. SM1(NEWQ)) GO TO 625
 577  622  IF (KFLAG .EQ. 0 .AND. RH .LT. 1.1D0) GO TO 610
 578  625  IF (KFLAG .LE. -2) RH = MIN(RH,0.2D0)
 579 C-----------------------------------------------------------------------
 580 C If there is a change of order, reset NQ, L, and the coefficients.
 581 C In any case H is reset according to RH and the YH array is rescaled.
 582 C Then exit from 690 if the step was OK, or redo the step otherwise.
 583 C-----------------------------------------------------------------------
 584       IF (NEWQ .EQ. NQ) GO TO 170
 585  630  NQ = NEWQ
 586       L = NQ + 1
 587       IRET = 2
 588       GO TO 150
 589 C-----------------------------------------------------------------------
 590 C Control reaches this section if 3 or more failures have occured.
 591 C If 10 failures have occurred, exit with KFLAG = -1.
 592 C It is assumed that the derivatives that have accumulated in the
 593 C YH array have errors of the wrong order.  Hence the first
 594 C derivative is recomputed, and the order is set to 1.  Then
 595 C H is reduced by a factor of 10, and the step is retried,
 596 C until it succeeds or H reaches HMIN.
 597 C-----------------------------------------------------------------------
 598  640  IF (KFLAG .EQ. -10) GO TO 660
 599       RH = 0.1D0
 600       RH = MAX(HMIN/ABS(H),RH)
 601       H = H*RH
 602       DO 645 I = 1,N
 603  645    Y(I) = YH(I,1)
 604       CALL F (NEQ, TN, Y, SAVF)
 605       NFE = NFE + 1
 606       DO 650 I = 1,N
 607  650    YH(I,2) = H*SAVF(I)
 608       IPUP = MITER
 609       IALTH = 5
 610       IF (NQ .EQ. 1) GO TO 200
 611       NQ = 1
 612       L = 2
 613       IRET = 3
 614       GO TO 150
 615 C-----------------------------------------------------------------------
 616 C All returns are made through this section.  H is saved in HOLD
 617 C to allow the caller to change H on the next step.
 618 C-----------------------------------------------------------------------
 619  660  KFLAG = -1
 620       GO TO 720
 621  670  KFLAG = -2
 622       GO TO 720
 623  680  KFLAG = -3
 624       GO TO 720
 625  690  RMAX = 10.0D0
 626  700  R = 1.0D0/TESCO(2,NQU)
 627       DO 710 I = 1,N
 628  710    ACOR(I) = ACOR(I)*R
 629  720  HOLD = H
 630       JSTART = 1
 631       RETURN
 632 C----------------------- End of Subroutine DSTODA ----------------------
 633       END