share/odepack/fortran/dstodi.f

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