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