1*DECK STOD 2 SUBROUTINE STOD (NEQ, Y, YH, NYH, YH1, EWT, SAVF, ACOR, WM, IWM, 3 + F, JAC, RPAR, IPAR) 4C***BEGIN PROLOGUE STOD 5C***SUBSIDIARY 6C***PURPOSE Subsidiary to DEBDF 7C***LIBRARY SLATEC 8C***TYPE SINGLE PRECISION (STOD-S, DSTOD-D) 9C***AUTHOR Watts, H. A., (SNLA) 10C***DESCRIPTION 11C 12C STOD integrates a system of first order odes over one step in the 13C integrator package DEBDF. 14C ---------------------------------------------------------------------- 15C STOD performs one step of the integration of an initial value 16C problem for a system of ordinary differential equations. 17C Note.. STOD is independent of the value of the iteration method 18C indicator MITER, when this is .NE. 0, and hence is independent 19C of the type of chord method used, or the Jacobian structure. 20C Communication with STOD is done with the following variables.. 21C 22C Y = An array of length .GE. n used as the Y argument in 23C all calls to F and JAC. 24C NEQ = Integer array containing problem size in NEQ(1), and 25C passed as the NEQ argument in all calls to F and JAC. 26C YH = An NYH by LMAX array containing the dependent variables 27C and their approximate scaled derivatives, where 28C LMAX = MAXORD + 1. YH(I,J+1) contains the approximate 29C J-th derivative of Y(I), scaled by H**J/Factorial(j) 30C (J = 0,1,...,NQ). On entry for the first step, the first 31C two columns of YH must be set from the initial values. 32C NYH = A constant integer .GE. N, the first dimension of YH. 33C YH1 = A one-dimensional array occupying the same space as YH. 34C EWT = An array of N elements with which the estimated local 35C errors in YH are compared. 36C SAVF = An array of working storage, of length N. 37C ACOR = A work array of length N, used for the accumulated 38C corrections. On a successful return, ACOR(I) contains 39C the estimated one-step local error in Y(I). 40C WM,IWM = Real and integer work arrays associated with matrix 41C operations in chord iteration (MITER .NE. 0). 42C PJAC = Name of routine to evaluate and preprocess Jacobian matrix 43C if a chord method is being used. 44C SLVS = Name of routine to solve linear system in chord iteration. 45C H = The step size to be attempted on the next step. 46C H is altered by the error control algorithm during the 47C problem. H can be either positive or negative, but its 48C sign must remain constant throughout the problem. 49C HMIN = The minimum absolute value of the step size H to be used. 50C HMXI = Inverse of the maximum absolute value of H to be used. 51C HMXI = 0.0 is allowed and corresponds to an infinite HMAX. 52C HMIN and HMXI may be changed at any time, but will not 53C take effect until the next change of H is considered. 54C TN = The independent variable. TN is updated on each step taken. 55C JSTART = An integer used for input only, with the following 56C values and meanings.. 57C 0 Perform the first step. 58C .GT.0 Take a new step continuing from the last. 59C -1 Take the next step with a new value of H, MAXORD, 60C N, METH, MITER, and/or matrix parameters. 61C -2 Take the next step with a new value of H, 62C but with other inputs unchanged. 63C On return, JSTART is set to 1 to facilitate continuation. 64C KFLAG = a completion code with the following meanings.. 65C 0 The step was successful. 66C -1 The requested error could not be achieved. 67C -2 Corrector convergence could not be achieved. 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 METH/MITER = The method flags. See description in driver. 75C N = The number of first-order differential equations. 76C ---------------------------------------------------------------------- 77C 78C***SEE ALSO DEBDF 79C***ROUTINES CALLED CFOD, PJAC, SLVS, VNWRMS 80C***COMMON BLOCKS DEBDF1 81C***REVISION HISTORY (YYMMDD) 82C 800901 DATE WRITTEN 83C 890531 Changed all specific intrinsics to generic. (WRB) 84C 891214 Prologue converted to Version 4.0 format. (BAB) 85C 900328 Added TYPE section. (WRB) 86C 910722 Updated AUTHOR section. (ALS) 87C 920422 Changed DIMENSION statement. (WRB) 88C***END PROLOGUE STOD 89 EXTERNAL F, JAC 90C 91CLLL. OPTIMIZE 92 INTEGER NEQ, NYH, IWM, I, I1, IALTH, IER, IOWND, IREDO, IRET, 93 1 IPUP, J, JB, JSTART, KFLAG, L, LMAX, M, MAXORD, MEO, METH, 94 2 MITER, N, NCF, NEWQ, NFE, NJE, NQ, NQNYH, NQU, NST, NSTEPJ 95 REAL Y, YH, YH1, EWT, SAVF, ACOR, WM, 96 1 ROWND, CONIT, CRATE, EL, ELCO, HOLD, RC, RMAX, TESCO, 97 2 EL0, H, HMIN, HMXI, HU, TN, UROUND, 98 3 DCON, DDN, DEL, DELP, DSM, DUP, EXDN, EXSM, EXUP, 99 4 R, RH, RHDN, RHSM, RHUP, TOLD, VNWRMS 100 DIMENSION Y(*), YH(NYH,*), YH1(*), EWT(*), SAVF(*), 101 1 ACOR(*), WM(*), IWM(*), RPAR(*), IPAR(*) 102 COMMON /DEBDF1/ ROWND, CONIT, CRATE, EL(13), ELCO(13,12), 103 1 HOLD, RC, RMAX, TESCO(3,12), 104 2 EL0, H, HMIN, HMXI, HU, TN, UROUND, IOWND(7), KSTEPS, IOD(6), 105 3 IALTH, IPUP, LMAX, MEO, NQNYH, NSTEPJ, 106 4 IER, JSTART, KFLAG, L, METH, MITER, MAXORD, N, NQ, NST, NFE, 107 5 NJE, NQU 108C 109C 110C***FIRST EXECUTABLE STATEMENT STOD 111 KFLAG = 0 112 TOLD = TN 113 NCF = 0 114 IF (JSTART .GT. 0) GO TO 200 115 IF (JSTART .EQ. -1) GO TO 100 116 IF (JSTART .EQ. -2) GO TO 160 117C----------------------------------------------------------------------- 118C ON THE FIRST CALL, THE ORDER IS SET TO 1, AND OTHER VARIABLES ARE 119C INITIALIZED. RMAX IS THE MAXIMUM RATIO BY WHICH H CAN BE INCREASED 120C IN A SINGLE STEP. IT IS INITIALLY 1.E4 TO COMPENSATE FOR THE SMALL 121C INITIAL H, BUT THEN IS NORMALLY EQUAL TO 10. IF A FAILURE 122C OCCURS (IN CORRECTOR CONVERGENCE OR ERROR TEST), RMAX IS SET AT 2 123C FOR THE NEXT INCREASE. 124C----------------------------------------------------------------------- 125 LMAX = MAXORD + 1 126 NQ = 1 127 L = 2 128 IALTH = 2 129 RMAX = 10000.0E0 130 RC = 0.0E0 131 EL0 = 1.0E0 132 CRATE = 0.7E0 133 DELP = 0.0E0 134 HOLD = H 135 MEO = METH 136 NSTEPJ = 0 137 IRET = 3 138 GO TO 140 139C----------------------------------------------------------------------- 140C THE FOLLOWING BLOCK HANDLES PRELIMINARIES NEEDED WHEN JSTART = -1. 141C IPUP IS SET TO MITER TO FORCE A MATRIX UPDATE. 142C IF AN ORDER INCREASE IS ABOUT TO BE CONSIDERED (IALTH = 1), 143C IALTH IS RESET TO 2 TO POSTPONE CONSIDERATION ONE MORE STEP. 144C IF THE CALLER HAS CHANGED METH, CFOD IS CALLED TO RESET 145C THE COEFFICIENTS OF THE METHOD. 146C IF THE CALLER HAS CHANGED MAXORD TO A VALUE LESS THAN THE CURRENT 147C ORDER NQ, NQ IS REDUCED TO MAXORD, AND A NEW H CHOSEN ACCORDINGLY. 148C IF H IS TO BE CHANGED, YH MUST BE RESCALED. 149C IF H OR METH IS BEING CHANGED, IALTH IS RESET TO L = NQ + 1 150C TO PREVENT FURTHER CHANGES IN H FOR THAT MANY STEPS. 151C----------------------------------------------------------------------- 152 100 IPUP = MITER 153 LMAX = MAXORD + 1 154 IF (IALTH .EQ. 1) IALTH = 2 155 IF (METH .EQ. MEO) GO TO 110 156 CALL CFOD (METH, ELCO, TESCO) 157 MEO = METH 158 IF (NQ .GT. MAXORD) GO TO 120 159 IALTH = L 160 IRET = 1 161 GO TO 150 162 110 IF (NQ .LE. MAXORD) GO TO 160 163 120 NQ = MAXORD 164 L = LMAX 165 DO 125 I = 1,L 166 125 EL(I) = ELCO(I,NQ) 167 NQNYH = NQ*NYH 168 RC = RC*EL(1)/EL0 169 EL0 = EL(1) 170 CONIT = 0.5E0/(NQ+2) 171 DDN = VNWRMS (N, SAVF, EWT)/TESCO(1,L) 172 EXDN = 1.0E0/L 173 RHDN = 1.0E0/(1.3E0*DDN**EXDN + 0.0000013E0) 174 RH = MIN(RHDN,1.0E0) 175 IREDO = 3 176 IF (H .EQ. HOLD) GO TO 170 177 RH = MIN(RH,ABS(H/HOLD)) 178 H = HOLD 179 GO TO 175 180C----------------------------------------------------------------------- 181C CFOD IS CALLED TO GET ALL THE INTEGRATION COEFFICIENTS FOR THE 182C CURRENT METH. THEN THE EL VECTOR AND RELATED CONSTANTS ARE RESET 183C WHENEVER THE ORDER NQ IS CHANGED, OR AT THE START OF THE PROBLEM. 184C----------------------------------------------------------------------- 185 140 CALL CFOD (METH, ELCO, TESCO) 186 150 DO 155 I = 1,L 187 155 EL(I) = ELCO(I,NQ) 188 NQNYH = NQ*NYH 189 RC = RC*EL(1)/EL0 190 EL0 = EL(1) 191 CONIT = 0.5E0/(NQ+2) 192 GO TO (160, 170, 200), IRET 193C----------------------------------------------------------------------- 194C IF H IS BEING CHANGED, THE H RATIO RH IS CHECKED AGAINST 195C RMAX, HMIN, AND HMXI, AND THE YH ARRAY RESCALED. IALTH IS SET TO 196C L = NQ + 1 TO PREVENT A CHANGE OF H FOR THAT MANY STEPS, UNLESS 197C FORCED BY A CONVERGENCE OR ERROR TEST FAILURE. 198C----------------------------------------------------------------------- 199 160 IF (H .EQ. HOLD) GO TO 200 200 RH = H/HOLD 201 H = HOLD 202 IREDO = 3 203 GO TO 175 204 170 RH = MAX(RH,HMIN/ABS(H)) 205 175 RH = MIN(RH,RMAX) 206 RH = RH/MAX(1.0E0,ABS(H)*HMXI*RH) 207 R = 1.0E0 208 DO 180 J = 2,L 209 R = R*RH 210 DO 180 I = 1,N 211 180 YH(I,J) = YH(I,J)*R 212 H = H*RH 213 RC = RC*RH 214 IALTH = L 215 IF (IREDO .EQ. 0) GO TO 680 216C----------------------------------------------------------------------- 217C THIS SECTION COMPUTES THE PREDICTED VALUES BY EFFECTIVELY 218C MULTIPLYING THE YH ARRAY BY THE PASCAL TRIANGLE MATRIX. 219C RC IS THE RATIO OF NEW TO OLD VALUES OF THE COEFFICIENT H*EL(1). 220C WHEN RC DIFFERS FROM 1 BY MORE THAN 30 PERCENT, IPUP IS SET TO MITER 221C TO FORCE PJAC TO BE CALLED, IF A JACOBIAN IS INVOLVED. 222C IN ANY CASE, PJAC IS CALLED AT LEAST EVERY 20-TH STEP. 223C----------------------------------------------------------------------- 224 200 IF (ABS(RC-1.0E0) .GT. 0.3E0) IPUP = MITER 225 IF (NST .GE. NSTEPJ+20) IPUP = MITER 226 TN = TN + H 227 I1 = NQNYH + 1 228 DO 215 JB = 1,NQ 229 I1 = I1 - NYH 230 DO 210 I = I1,NQNYH 231 210 YH1(I) = YH1(I) + YH1(I+NYH) 232 215 CONTINUE 233 KSTEPS = KSTEPS + 1 234C----------------------------------------------------------------------- 235C UP TO 3 CORRECTOR ITERATIONS ARE TAKEN. A CONVERGENCE TEST IS 236C MADE ON THE R.M.S. NORM OF EACH CORRECTION, WEIGHTED BY THE ERROR 237C WEIGHT VECTOR EWT. THE SUM OF THE CORRECTIONS IS ACCUMULATED IN THE 238C VECTOR ACOR(I). THE YH ARRAY IS NOT ALTERED IN THE CORRECTOR LOOP. 239C----------------------------------------------------------------------- 240 220 M = 0 241 DO 230 I = 1,N 242 230 Y(I) = YH(I,1) 243 CALL F (TN, Y, SAVF, RPAR, IPAR) 244 NFE = NFE + 1 245 IF (IPUP .LE. 0) GO TO 250 246C----------------------------------------------------------------------- 247C IF INDICATED, THE MATRIX P = I - H*EL(1)*J IS REEVALUATED AND 248C PREPROCESSED BEFORE STARTING THE CORRECTOR ITERATION. IPUP IS SET 249C TO 0 AS AN INDICATOR THAT THIS HAS BEEN DONE. 250C----------------------------------------------------------------------- 251 IPUP = 0 252 RC = 1.0E0 253 NSTEPJ = NST 254 CRATE = 0.7E0 255 CALL PJAC (NEQ, Y, YH, NYH, EWT, ACOR, SAVF, WM, IWM, F, JAC, 256 1 RPAR, IPAR) 257 IF (IER .NE. 0) GO TO 430 258 250 DO 260 I = 1,N 259 260 ACOR(I) = 0.0E0 260 270 IF (MITER .NE. 0) GO TO 350 261C----------------------------------------------------------------------- 262C IN THE CASE OF FUNCTIONAL ITERATION, UPDATE Y DIRECTLY FROM 263C THE RESULT OF THE LAST FUNCTION EVALUATION. 264C----------------------------------------------------------------------- 265 DO 290 I = 1,N 266 SAVF(I) = H*SAVF(I) - YH(I,2) 267 290 Y(I) = SAVF(I) - ACOR(I) 268 DEL = VNWRMS (N, Y, EWT) 269 DO 300 I = 1,N 270 Y(I) = YH(I,1) + EL(1)*SAVF(I) 271 300 ACOR(I) = SAVF(I) 272 GO TO 400 273C----------------------------------------------------------------------- 274C IN THE CASE OF THE CHORD METHOD, COMPUTE THE CORRECTOR ERROR, 275C AND SOLVE THE LINEAR SYSTEM WITH THAT AS RIGHT-HAND SIDE AND 276C P AS COEFFICIENT MATRIX. 277C----------------------------------------------------------------------- 278 350 DO 360 I = 1,N 279 360 Y(I) = H*SAVF(I) - (YH(I,2) + ACOR(I)) 280 CALL SLVS (WM, IWM, Y, SAVF) 281 IF (IER .NE. 0) GO TO 410 282 DEL = VNWRMS (N, Y, EWT) 283 DO 380 I = 1,N 284 ACOR(I) = ACOR(I) + Y(I) 285 380 Y(I) = YH(I,1) + EL(1)*ACOR(I) 286C----------------------------------------------------------------------- 287C TEST FOR CONVERGENCE. IF M.GT.0, AN ESTIMATE OF THE CONVERGENCE 288C RATE CONSTANT IS STORED IN CRATE, AND THIS IS USED IN THE TEST. 289C----------------------------------------------------------------------- 290 400 IF (M .NE. 0) CRATE = MAX(0.2E0*CRATE,DEL/DELP) 291 DCON = DEL*MIN(1.0E0,1.5E0*CRATE)/(TESCO(2,NQ)*CONIT) 292 IF (DCON .LE. 1.0E0) GO TO 450 293 M = M + 1 294 IF (M .EQ. 3) GO TO 410 295 IF (M .GE. 2 .AND. DEL .GT. 2.0E0*DELP) GO TO 410 296 DELP = DEL 297 CALL F (TN, Y, SAVF, RPAR, IPAR) 298 NFE = NFE + 1 299 GO TO 270 300C----------------------------------------------------------------------- 301C THE CORRECTOR ITERATION FAILED TO CONVERGE IN 3 TRIES. 302C IF MITER .NE. 0 AND THE JACOBIAN IS OUT OF DATE, PJAC IS CALLED FOR 303C THE NEXT TRY. OTHERWISE THE YH ARRAY IS RETRACTED TO ITS VALUES 304C BEFORE PREDICTION, AND H IS REDUCED, IF POSSIBLE. IF H CANNOT BE 305C REDUCED OR 10 FAILURES HAVE OCCURRED, EXIT WITH KFLAG = -2. 306C----------------------------------------------------------------------- 307 410 IF (IPUP .EQ. 0) GO TO 430 308 IPUP = MITER 309 GO TO 220 310 430 TN = TOLD 311 NCF = NCF + 1 312 RMAX = 2.0E0 313 I1 = NQNYH + 1 314 DO 445 JB = 1,NQ 315 I1 = I1 - NYH 316 DO 440 I = I1,NQNYH 317 440 YH1(I) = YH1(I) - YH1(I+NYH) 318 445 CONTINUE 319 IF (ABS(H) .LE. HMIN*1.00001E0) GO TO 670 320 IF (NCF .EQ. 10) GO TO 670 321 RH = 0.25E0 322 IPUP = MITER 323 IREDO = 1 324 GO TO 170 325C----------------------------------------------------------------------- 326C THE CORRECTOR HAS CONVERGED. IPUP IS SET TO -1 IF MITER .NE. 0, 327C TO SIGNAL THAT THE JACOBIAN INVOLVED MAY NEED UPDATING LATER. 328C THE LOCAL ERROR TEST IS MADE AND CONTROL PASSES TO STATEMENT 500 329C IF IT FAILS. 330C----------------------------------------------------------------------- 331 450 IF (MITER .NE. 0) IPUP = -1 332 IF (M .EQ. 0) DSM = DEL/TESCO(2,NQ) 333 IF (M .GT. 0) DSM = VNWRMS (N, ACOR, EWT)/TESCO(2,NQ) 334 IF (DSM .GT. 1.0E0) GO TO 500 335C----------------------------------------------------------------------- 336C AFTER A SUCCESSFUL STEP, UPDATE THE YH ARRAY. 337C CONSIDER CHANGING H IF IALTH = 1. OTHERWISE DECREASE IALTH BY 1. 338C IF IALTH IS THEN 1 AND NQ .LT. MAXORD, THEN ACOR IS SAVED FOR 339C USE IN A POSSIBLE ORDER INCREASE ON THE NEXT STEP. 340C IF A CHANGE IN H IS CONSIDERED, AN INCREASE OR DECREASE IN ORDER 341C BY ONE IS CONSIDERED ALSO. A CHANGE IN H IS MADE ONLY IF IT IS BY A 342C FACTOR OF AT LEAST 1.1. IF NOT, IALTH IS SET TO 3 TO PREVENT 343C TESTING FOR THAT MANY STEPS. 344C----------------------------------------------------------------------- 345 KFLAG = 0 346 IREDO = 0 347 NST = NST + 1 348 HU = H 349 NQU = NQ 350 DO 470 J = 1,L 351 DO 470 I = 1,N 352 470 YH(I,J) = YH(I,J) + EL(J)*ACOR(I) 353 IALTH = IALTH - 1 354 IF (IALTH .EQ. 0) GO TO 520 355 IF (IALTH .GT. 1) GO TO 690 356 IF (L .EQ. LMAX) GO TO 690 357 DO 490 I = 1,N 358 490 YH(I,LMAX) = ACOR(I) 359 GO TO 690 360C----------------------------------------------------------------------- 361C THE ERROR TEST FAILED. KFLAG KEEPS TRACK OF MULTIPLE FAILURES. 362C RESTORE TN AND THE YH ARRAY TO THEIR PREVIOUS VALUES, AND PREPARE 363C TO TRY THE STEP AGAIN. COMPUTE THE OPTIMUM STEP SIZE FOR THIS OR 364C ONE LOWER ORDER. AFTER 2 OR MORE FAILURES, H IS FORCED TO DECREASE 365C BY A FACTOR OF 0.2 OR LESS. 366C----------------------------------------------------------------------- 367 500 KFLAG = KFLAG - 1 368 TN = TOLD 369 I1 = NQNYH + 1 370 DO 515 JB = 1,NQ 371 I1 = I1 - NYH 372 DO 510 I = I1,NQNYH 373 510 YH1(I) = YH1(I) - YH1(I+NYH) 374 515 CONTINUE 375 RMAX = 2.0E0 376 IF (ABS(H) .LE. HMIN*1.00001E0) GO TO 660 377 IF (KFLAG .LE. -3) GO TO 640 378 IREDO = 2 379 RHUP = 0.0E0 380 GO TO 540 381C----------------------------------------------------------------------- 382C REGARDLESS OF THE SUCCESS OR FAILURE OF THE STEP, FACTORS 383C RHDN, RHSM, AND RHUP ARE COMPUTED, BY WHICH H COULD BE MULTIPLIED 384C AT ORDER NQ - 1, ORDER NQ, OR ORDER NQ + 1, RESPECTIVELY. 385C IN THE CASE OF FAILURE, RHUP = 0.0 TO AVOID AN ORDER INCREASE. 386C THE LARGEST OF THESE IS DETERMINED AND THE NEW ORDER CHOSEN 387C ACCORDINGLY. IF THE ORDER IS TO BE INCREASED, WE COMPUTE ONE 388C ADDITIONAL SCALED DERIVATIVE. 389C----------------------------------------------------------------------- 390 520 RHUP = 0.0E0 391 IF (L .EQ. LMAX) GO TO 540 392 DO 530 I = 1,N 393 530 SAVF(I) = ACOR(I) - YH(I,LMAX) 394 DUP = VNWRMS (N, SAVF, EWT)/TESCO(3,NQ) 395 EXUP = 1.0E0/(L+1) 396 RHUP = 1.0E0/(1.4E0*DUP**EXUP + 0.0000014E0) 397 540 EXSM = 1.0E0/L 398 RHSM = 1.0E0/(1.2E0*DSM**EXSM + 0.0000012E0) 399 RHDN = 0.0E0 400 IF (NQ .EQ. 1) GO TO 560 401 DDN = VNWRMS (N, YH(1,L), EWT)/TESCO(1,NQ) 402 EXDN = 1.0E0/NQ 403 RHDN = 1.0E0/(1.3E0*DDN**EXDN + 0.0000013E0) 404 560 IF (RHSM .GE. RHUP) GO TO 570 405 IF (RHUP .GT. RHDN) GO TO 590 406 GO TO 580 407 570 IF (RHSM .LT. RHDN) GO TO 580 408 NEWQ = NQ 409 RH = RHSM 410 GO TO 620 411 580 NEWQ = NQ - 1 412 RH = RHDN 413 IF (KFLAG .LT. 0 .AND. RH .GT. 1.0E0) RH = 1.0E0 414 GO TO 620 415 590 NEWQ = L 416 RH = RHUP 417 IF (RH .LT. 1.1E0) GO TO 610 418 R = EL(L)/L 419 DO 600 I = 1,N 420 600 YH(I,NEWQ+1) = ACOR(I)*R 421 GO TO 630 422 610 IALTH = 3 423 GO TO 690 424 620 IF ((KFLAG .EQ. 0) .AND. (RH .LT. 1.1E0)) GO TO 610 425 IF (KFLAG .LE. -2) RH = MIN(RH,0.2E0) 426C----------------------------------------------------------------------- 427C IF THERE IS A CHANGE OF ORDER, RESET NQ, L, AND THE COEFFICIENTS. 428C IN ANY CASE H IS RESET ACCORDING TO RH AND THE YH ARRAY IS RESCALED. 429C THEN EXIT FROM 680 IF THE STEP WAS OK, OR REDO THE STEP OTHERWISE. 430C----------------------------------------------------------------------- 431 IF (NEWQ .EQ. NQ) GO TO 170 432 630 NQ = NEWQ 433 L = NQ + 1 434 IRET = 2 435 GO TO 150 436C----------------------------------------------------------------------- 437C CONTROL REACHES THIS SECTION IF 3 OR MORE FAILURES HAVE OCCURRED. 438C IF 10 FAILURES HAVE OCCURRED, EXIT WITH KFLAG = -1. 439C IT IS ASSUMED THAT THE DERIVATIVES THAT HAVE ACCUMULATED IN THE 440C YH ARRAY HAVE ERRORS OF THE WRONG ORDER. HENCE THE FIRST 441C DERIVATIVE IS RECOMPUTED, AND THE ORDER IS SET TO 1. THEN 442C H IS REDUCED BY A FACTOR OF 10, AND THE STEP IS RETRIED, 443C UNTIL IT SUCCEEDS OR H REACHES HMIN. 444C----------------------------------------------------------------------- 445 640 IF (KFLAG .EQ. -10) GO TO 660 446 RH = 0.1E0 447 RH = MAX(HMIN/ABS(H),RH) 448 H = H*RH 449 DO 645 I = 1,N 450 645 Y(I) = YH(I,1) 451 CALL F (TN, Y, SAVF, RPAR, IPAR) 452 NFE = NFE + 1 453 DO 650 I = 1,N 454 650 YH(I,2) = H*SAVF(I) 455 IPUP = MITER 456 IALTH = 5 457 IF (NQ .EQ. 1) GO TO 200 458 NQ = 1 459 L = 2 460 IRET = 3 461 GO TO 150 462C----------------------------------------------------------------------- 463C ALL RETURNS ARE MADE THROUGH THIS SECTION. H IS SAVED IN HOLD 464C TO ALLOW THE CALLER TO CHANGE H ON THE NEXT STEP. 465C----------------------------------------------------------------------- 466 660 KFLAG = -1 467 GO TO 700 468 670 KFLAG = -2 469 GO TO 700 470 680 RMAX = 10.0E0 471 690 R = 1.0E0/TESCO(2,NQU) 472 DO 695 I = 1,N 473 695 ACOR(I) = ACOR(I)*R 474 700 HOLD = H 475 JSTART = 1 476 RETURN 477C----------------------- END OF SUBROUTINE STOD ----------------------- 478 END 479