share/odepack/fortran/dstodpk.f

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