1 SUBROUTINE STODE (NEQ, Y, YH, NYH, YH1, EWT, SAVF, ACOR, 2 1 WM, IWM, F, JAC, PJAC, SLVS, IERR) 3CLLL. OPTIMIZE 4 EXTERNAL F, JAC, PJAC, SLVS 5 INTEGER NEQ, NYH, IWM 6 INTEGER ILLIN, INIT, LYH, LEWT, LACOR, LSAVF, LWM, LIWM, 7 1 MXSTEP, MXHNIL, NHNIL, NTREP, NSLAST, CNYH, 8 2 IALTH, IPUP, LMAX, MEO, NQNYH, NSLP 9 INTEGER ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L, METH, MITER, 10 1 MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU 11 INTEGER I, I1, IREDO, IRET, J, JB, M, NCF, NEWQ 12 DOUBLE PRECISION Y, YH, YH1, EWT, SAVF, ACOR, WM 13 DOUBLE PRECISION CONIT, CRATE, EL, ELCO, HOLD, RMAX, TESCO, 14 2 CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND 15 DOUBLE PRECISION DCON, DDN, DEL, DELP, DSM, DUP, EXDN, EXSM, EXUP, 16 1 R, RH, RHDN, RHSM, RHUP, TOLD, VNORM 17 DIMENSION NEQ(*), Y(*), YH(NYH,*), YH1(*), EWT(*), SAVF(*), 18 1 ACOR(*), WM(*), IWM(*) 19 COMMON /LS0001/ CONIT, CRATE, EL(13), ELCO(13,12), 20 1 HOLD, RMAX, TESCO(3,12), 21 2 CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND, 22 2 ILLIN, INIT, LYH, LEWT, LACOR, LSAVF, LWM, LIWM, 23 3 MXSTEP, MXHNIL, NHNIL, NTREP, NSLAST, CNYH, 24 3 IALTH, IPUP, LMAX, MEO, NQNYH, NSLP, 25 4 ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L, METH, MITER, 26 5 MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU 27C----------------------------------------------------------------------- 28C STODE PERFORMS ONE STEP OF THE INTEGRATION OF AN INITIAL VALUE 29C PROBLEM FOR A SYSTEM OF ORDINARY DIFFERENTIAL EQUATIONS. 30C NOTE.. STODE IS INDEPENDENT OF THE VALUE OF THE ITERATION METHOD 31C INDICATOR MITER, WHEN THIS IS .NE. 0, AND HENCE IS INDEPENDENT 32C OF THE TYPE OF CHORD METHOD USED, OR THE JACOBIAN STRUCTURE. 33C COMMUNICATION WITH STODE IS DONE WITH THE FOLLOWING VARIABLES.. 34C 35C NEQ = INTEGER ARRAY CONTAINING PROBLEM SIZE IN NEQ(1), AND 36C PASSED AS THE NEQ ARGUMENT IN ALL CALLS TO F AND JAC. 37C Y = AN ARRAY OF LENGTH .GE. N USED AS THE Y ARGUMENT IN 38C ALL CALLS TO F AND JAC. 39C YH = AN NYH BY LMAX ARRAY CONTAINING THE DEPENDENT VARIABLES 40C AND THEIR APPROXIMATE SCALED DERIVATIVES, WHERE 41C LMAX = MAXORD + 1. YH(I,J+1) CONTAINS THE APPROXIMATE 42C J-TH DERIVATIVE OF Y(I), SCALED BY H**J/FACTORIAL(J) 43C (J = 0,1,...,NQ). ON ENTRY FOR THE FIRST STEP, THE FIRST 44C TWO COLUMNS OF YH MUST BE SET FROM THE INITIAL VALUES. 45C NYH = A CONSTANT INTEGER .GE. N, THE FIRST DIMENSION OF YH. 46C YH1 = A ONE-DIMENSIONAL ARRAY OCCUPYING THE SAME SPACE AS YH. 47C EWT = AN ARRAY OF LENGTH N CONTAINING MULTIPLICATIVE WEIGHTS 48C FOR LOCAL ERROR MEASUREMENTS. LOCAL ERRORS IN Y(I) ARE 49C COMPARED TO 1.0/EWT(I) IN VARIOUS ERROR TESTS. 50C SAVF = AN ARRAY OF WORKING STORAGE, OF LENGTH N. 51C ALSO USED FOR INPUT OF YH(*,MAXORD+2) WHEN JSTART = -1 52C AND MAXORD .LT. THE CURRENT ORDER NQ. 53C ACOR = A WORK ARRAY OF LENGTH N, USED FOR THE ACCUMULATED 54C CORRECTIONS. ON A SUCCESSFUL RETURN, ACOR(I) CONTAINS 55C THE ESTIMATED ONE-STEP LOCAL ERROR IN Y(I). 56C WM,IWM = REAL AND INTEGER WORK ARRAYS ASSOCIATED WITH MATRIX 57C OPERATIONS IN CHORD ITERATION (MITER .NE. 0). 58C PJAC = NAME OF ROUTINE TO EVALUATE AND PREPROCESS JACOBIAN MATRIX 59C AND P = I - H*EL0*JAC, IF A CHORD METHOD IS BEING USED. 60C SLVS = NAME OF ROUTINE TO SOLVE LINEAR SYSTEM IN CHORD ITERATION. 61C CCMAX = MAXIMUM RELATIVE CHANGE IN H*EL0 BEFORE PJAC IS CALLED. 62C H = THE STEP SIZE TO BE ATTEMPTED ON THE NEXT STEP. 63C H IS ALTERED BY THE ERROR CONTROL ALGORITHM DURING THE 64C PROBLEM. H CAN BE EITHER POSITIVE OR NEGATIVE, BUT ITS 65C SIGN MUST REMAIN CONSTANT THROUGHOUT THE PROBLEM. 66C HMIN = THE MINIMUM ABSOLUTE VALUE OF THE STEP SIZE H TO BE USED. 67C HMXI = INVERSE OF THE MAXIMUM ABSOLUTE VALUE OF H TO BE USED. 68C HMXI = 0.0 IS ALLOWED AND CORRESPONDS TO AN INFINITE HMAX. 69C HMIN AND HMXI MAY BE CHANGED AT ANY TIME, BUT WILL NOT 70C TAKE EFFECT UNTIL THE NEXT CHANGE OF H IS CONSIDERED. 71C TN = THE INDEPENDENT VARIABLE. TN IS UPDATED ON EACH STEP TAKEN. 72C JSTART = AN INTEGER USED FOR INPUT ONLY, WITH THE FOLLOWING 73C VALUES AND MEANINGS.. 74C 0 PERFORM THE FIRST STEP. 75C .GT.0 TAKE A NEW STEP CONTINUING FROM THE LAST. 76C -1 TAKE THE NEXT STEP WITH A NEW VALUE OF H, MAXORD, 77C N, METH, MITER, AND/OR MATRIX PARAMETERS. 78C -2 TAKE THE NEXT STEP WITH A NEW VALUE OF H, 79C BUT WITH OTHER INPUTS UNCHANGED. 80C ON RETURN, JSTART IS SET TO 1 TO FACILITATE CONTINUATION. 81C KFLAG = A COMPLETION CODE WITH THE FOLLOWING MEANINGS.. 82C 0 THE STEP WAS SUCCESFUL. 83C -1 THE REQUESTED ERROR COULD NOT BE ACHIEVED. 84C -2 CORRECTOR CONVERGENCE COULD NOT BE ACHIEVED. 85C -3 FATAL ERROR IN PJAC OR SLVS. 86C A RETURN WITH KFLAG = -1 OR -2 MEANS EITHER 87C ABS(H) = HMIN OR 10 CONSECUTIVE FAILURES OCCURRED. 88C ON A RETURN WITH KFLAG NEGATIVE, THE VALUES OF TN AND 89C THE YH ARRAY ARE AS OF THE BEGINNING OF THE LAST 90C STEP, AND H IS THE LAST STEP SIZE ATTEMPTED. 91C MAXORD = THE MAXIMUM ORDER OF INTEGRATION METHOD TO BE ALLOWED. 92C MAXCOR = THE MAXIMUM NUMBER OF CORRECTOR ITERATIONS ALLOWED. 93C MSBP = MAXIMUM NUMBER OF STEPS BETWEEN PJAC CALLS (MITER .GT. 0). 94C MXNCF = MAXIMUM NUMBER OF CONVERGENCE FAILURES ALLOWED. 95C METH/MITER = THE METHOD FLAGS. SEE DESCRIPTION IN DRIVER. 96C N = THE NUMBER OF FIRST-ORDER DIFFERENTIAL EQUATIONS. 97C IERR = ERROR FLAG FROM USER-SUPPLIED FUNCTION 98C----------------------------------------------------------------------- 99 KFLAG = 0 100 TOLD = TN 101 NCF = 0 102 IERPJ = 0 103 IERSL = 0 104 JCUR = 0 105 ICF = 0 106 DELP = 0.0D0 107 IF (JSTART .GT. 0) GO TO 200 108 IF (JSTART .EQ. -1) GO TO 100 109 IF (JSTART .EQ. -2) GO TO 160 110C----------------------------------------------------------------------- 111C ON THE FIRST CALL, THE ORDER IS SET TO 1, AND OTHER VARIABLES ARE 112C INITIALIZED. RMAX IS THE MAXIMUM RATIO BY WHICH H CAN BE INCREASED 113C IN A SINGLE STEP. IT IS INITIALLY 1.E4 TO COMPENSATE FOR THE SMALL 114C INITIAL H, BUT THEN IS NORMALLY EQUAL TO 10. IF A FAILURE 115C OCCURS (IN CORRECTOR CONVERGENCE OR ERROR TEST), RMAX IS SET AT 2 116C FOR THE NEXT INCREASE. 117C----------------------------------------------------------------------- 118 LMAX = MAXORD + 1 119 NQ = 1 120 L = 2 121 IALTH = 2 122 RMAX = 10000.0D0 123 RC = 0.0D0 124 EL0 = 1.0D0 125 CRATE = 0.7D0 126 HOLD = H 127 MEO = METH 128 NSLP = 0 129 IPUP = MITER 130 IRET = 3 131 GO TO 140 132C----------------------------------------------------------------------- 133C THE FOLLOWING BLOCK HANDLES PRELIMINARIES NEEDED WHEN JSTART = -1. 134C IPUP IS SET TO MITER TO FORCE A MATRIX UPDATE. 135C IF AN ORDER INCREASE IS ABOUT TO BE CONSIDERED (IALTH = 1), 136C IALTH IS RESET TO 2 TO POSTPONE CONSIDERATION ONE MORE STEP. 137C IF THE CALLER HAS CHANGED METH, CFODE IS CALLED TO RESET 138C THE COEFFICIENTS OF THE METHOD. 139C IF THE CALLER HAS CHANGED MAXORD TO A VALUE LESS THAN THE CURRENT 140C ORDER NQ, NQ IS REDUCED TO MAXORD, AND A NEW H CHOSEN ACCORDINGLY. 141C IF H IS TO BE CHANGED, YH MUST BE RESCALED. 142C IF H OR METH IS BEING CHANGED, IALTH IS RESET TO L = NQ + 1 143C TO PREVENT FURTHER CHANGES IN H FOR THAT MANY STEPS. 144C----------------------------------------------------------------------- 145 100 IPUP = MITER 146 LMAX = MAXORD + 1 147 IF (IALTH .EQ. 1) IALTH = 2 148 IF (METH .EQ. MEO) GO TO 110 149 CALL CFODE (METH, ELCO, TESCO) 150 MEO = METH 151 IF (NQ .GT. MAXORD) GO TO 120 152 IALTH = L 153 IRET = 1 154 GO TO 150 155 110 IF (NQ .LE. MAXORD) GO TO 160 156 120 NQ = MAXORD 157 L = LMAX 158 DO 125 I = 1,L 159 125 EL(I) = ELCO(I,NQ) 160 NQNYH = NQ*NYH 161 RC = RC*EL(1)/EL0 162 EL0 = EL(1) 163 CONIT = 0.5D0/DBLE(NQ+2) 164 DDN = VNORM (N, SAVF, EWT)/TESCO(1,L) 165 EXDN = 1.0D0/DBLE(L) 166 RHDN = 1.0D0/(1.3D0*DDN**EXDN + 0.0000013D0) 167 RH = DMIN1(RHDN,1.0D0) 168 IREDO = 3 169 IF (H .EQ. HOLD) GO TO 170 170 RH = DMIN1(RH,DABS(H/HOLD)) 171 H = HOLD 172 GO TO 175 173C----------------------------------------------------------------------- 174C CFODE IS CALLED TO GET ALL THE INTEGRATION COEFFICIENTS FOR THE 175C CURRENT METH. THEN THE EL VECTOR AND RELATED CONSTANTS ARE RESET 176C WHENEVER THE ORDER NQ IS CHANGED, OR AT THE START OF THE PROBLEM. 177C----------------------------------------------------------------------- 178 140 CALL CFODE (METH, ELCO, TESCO) 179 150 DO 155 I = 1,L 180 155 EL(I) = ELCO(I,NQ) 181 NQNYH = NQ*NYH 182 RC = RC*EL(1)/EL0 183 EL0 = EL(1) 184 CONIT = 0.5D0/DBLE(NQ+2) 185 GO TO (160, 170, 200), IRET 186C----------------------------------------------------------------------- 187C IF H IS BEING CHANGED, THE H RATIO RH IS CHECKED AGAINST 188C RMAX, HMIN, AND HMXI, AND THE YH ARRAY RESCALED. IALTH IS SET TO 189C L = NQ + 1 TO PREVENT A CHANGE OF H FOR THAT MANY STEPS, UNLESS 190C FORCED BY A CONVERGENCE OR ERROR TEST FAILURE. 191C----------------------------------------------------------------------- 192 160 IF (H .EQ. HOLD) GO TO 200 193 RH = H/HOLD 194 H = HOLD 195 IREDO = 3 196 GO TO 175 197 170 RH = DMAX1(RH,HMIN/DABS(H)) 198 175 RH = DMIN1(RH,RMAX) 199 RH = RH/DMAX1(1.0D0,DABS(H)*HMXI*RH) 200 R = 1.0D0 201 DO 180 J = 2,L 202 R = R*RH 203 DO 180 I = 1,N 204 180 YH(I,J) = YH(I,J)*R 205 H = H*RH 206 RC = RC*RH 207 IALTH = L 208 IF (IREDO .EQ. 0) GO TO 690 209C----------------------------------------------------------------------- 210C THIS SECTION COMPUTES THE PREDICTED VALUES BY EFFECTIVELY 211C MULTIPLYING THE YH ARRAY BY THE PASCAL TRIANGLE MATRIX. 212C RC IS THE RATIO OF NEW TO OLD VALUES OF THE COEFFICIENT H*EL(1). 213C WHEN RC DIFFERS FROM 1 BY MORE THAN CCMAX, IPUP IS SET TO MITER 214C TO FORCE PJAC TO BE CALLED, IF A JACOBIAN IS INVOLVED. 215C IN ANY CASE, PJAC IS CALLED AT LEAST EVERY MSBP STEPS. 216C----------------------------------------------------------------------- 217 200 IF (DABS(RC-1.0D0) .GT. CCMAX) IPUP = MITER 218 IF (NST .GE. NSLP+MSBP) IPUP = MITER 219 TN = TN + H 220 I1 = NQNYH + 1 221 DO 215 JB = 1,NQ 222 I1 = I1 - NYH 223CDIR$ IVDEP 224 DO 210 I = I1,NQNYH 225 210 YH1(I) = YH1(I) + YH1(I+NYH) 226 215 CONTINUE 227C----------------------------------------------------------------------- 228C UP TO MAXCOR CORRECTOR ITERATIONS ARE TAKEN. A CONVERGENCE TEST IS 229C MADE ON THE R.M.S. NORM OF EACH CORRECTION, WEIGHTED BY THE ERROR 230C WEIGHT VECTOR EWT. THE SUM OF THE CORRECTIONS IS ACCUMULATED IN THE 231C VECTOR ACOR(I). THE YH ARRAY IS NOT ALTERED IN THE CORRECTOR LOOP. 232C----------------------------------------------------------------------- 233 220 M = 0 234 DO 230 I = 1,N 235 230 Y(I) = YH(I,1) 236 IERR = 0 237 CALL F (NEQ, TN, Y, SAVF, IERR) 238 IF (IERR .LT. 0) RETURN 239 NFE = NFE + 1 240 IF (IPUP .LE. 0) GO TO 250 241C----------------------------------------------------------------------- 242C IF INDICATED, THE MATRIX P = I - H*EL(1)*J IS REEVALUATED AND 243C PREPROCESSED BEFORE STARTING THE CORRECTOR ITERATION. IPUP IS SET 244C TO 0 AS AN INDICATOR THAT THIS HAS BEEN DONE. 245C----------------------------------------------------------------------- 246 IERR = 0 247 CALL PJAC (NEQ, Y, YH, NYH, EWT, ACOR, SAVF, WM, IWM, F, JAC, 248 1 IERR) 249 IF (IERR .LT. 0) RETURN 250 IPUP = 0 251 RC = 1.0D0 252 NSLP = NST 253 CRATE = 0.7D0 254 IF (IERPJ .NE. 0) GO TO 430 255 250 DO 260 I = 1,N 256 260 ACOR(I) = 0.0D0 257 270 IF (MITER .NE. 0) GO TO 350 258C----------------------------------------------------------------------- 259C IN THE CASE OF FUNCTIONAL ITERATION, UPDATE Y DIRECTLY FROM 260C THE RESULT OF THE LAST FUNCTION EVALUATION. 261C----------------------------------------------------------------------- 262 DO 290 I = 1,N 263 SAVF(I) = H*SAVF(I) - YH(I,2) 264 290 Y(I) = SAVF(I) - ACOR(I) 265 DEL = VNORM (N, Y, EWT) 266 DO 300 I = 1,N 267 Y(I) = YH(I,1) + EL(1)*SAVF(I) 268 300 ACOR(I) = SAVF(I) 269 GO TO 400 270C----------------------------------------------------------------------- 271C IN THE CASE OF THE CHORD METHOD, COMPUTE THE CORRECTOR ERROR, 272C AND SOLVE THE LINEAR SYSTEM WITH THAT AS RIGHT-HAND SIDE AND 273C P AS COEFFICIENT MATRIX. 274C----------------------------------------------------------------------- 275 350 DO 360 I = 1,N 276 360 Y(I) = H*SAVF(I) - (YH(I,2) + ACOR(I)) 277 CALL SLVS (WM, IWM, Y, SAVF) 278 IF (IERSL .LT. 0) GO TO 430 279 IF (IERSL .GT. 0) GO TO 410 280 DEL = VNORM (N, Y, EWT) 281 DO 380 I = 1,N 282 ACOR(I) = ACOR(I) + Y(I) 283 380 Y(I) = YH(I,1) + EL(1)*ACOR(I) 284C----------------------------------------------------------------------- 285C TEST FOR CONVERGENCE. IF M.GT.0, AN ESTIMATE OF THE CONVERGENCE 286C RATE CONSTANT IS STORED IN CRATE, AND THIS IS USED IN THE TEST. 287C----------------------------------------------------------------------- 288 400 IF (M .NE. 0) CRATE = DMAX1(0.2D0*CRATE,DEL/DELP) 289 DCON = DEL*DMIN1(1.0D0,1.5D0*CRATE)/(TESCO(2,NQ)*CONIT) 290 IF (DCON .LE. 1.0D0) GO TO 450 291 M = M + 1 292 IF (M .EQ. MAXCOR) GO TO 410 293 IF (M .GE. 2 .AND. DEL .GT. 2.0D0*DELP) GO TO 410 294 DELP = DEL 295 IERR = 0 296 CALL F (NEQ, TN, Y, SAVF, IERR) 297 IF (IERR .LT. 0) RETURN 298 NFE = NFE + 1 299 GO TO 270 300C----------------------------------------------------------------------- 301C THE CORRECTOR ITERATION FAILED TO CONVERGE. 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 MXNCF FAILURES HAVE OCCURRED, EXIT WITH KFLAG = -2. 306C----------------------------------------------------------------------- 307 410 IF (MITER .EQ. 0 .OR. JCUR .EQ. 1) GO TO 430 308 ICF = 1 309 IPUP = MITER 310 GO TO 220 311 430 ICF = 2 312 NCF = NCF + 1 313 RMAX = 2.0D0 314 TN = TOLD 315 I1 = NQNYH + 1 316 DO 445 JB = 1,NQ 317 I1 = I1 - NYH 318CDIR$ IVDEP 319 DO 440 I = I1,NQNYH 320 440 YH1(I) = YH1(I) - YH1(I+NYH) 321 445 CONTINUE 322 IF (IERPJ .LT. 0 .OR. IERSL .LT. 0) GO TO 680 323 IF (DABS(H) .LE. HMIN*1.00001D0) GO TO 670 324 IF (NCF .EQ. MXNCF) GO TO 670 325 RH = 0.25D0 326 IPUP = MITER 327 IREDO = 1 328 GO TO 170 329C----------------------------------------------------------------------- 330C THE CORRECTOR HAS CONVERGED. JCUR IS SET TO 0 331C TO SIGNAL THAT THE JACOBIAN INVOLVED MAY NEED UPDATING LATER. 332C THE LOCAL ERROR TEST IS MADE AND CONTROL PASSES TO STATEMENT 500 333C IF IT FAILS. 334C----------------------------------------------------------------------- 335 450 JCUR = 0 336 IF (M .EQ. 0) DSM = DEL/TESCO(2,NQ) 337 IF (M .GT. 0) DSM = VNORM (N, ACOR, EWT)/TESCO(2,NQ) 338 IF (DSM .GT. 1.0D0) GO TO 500 339C----------------------------------------------------------------------- 340C AFTER A SUCCESSFUL STEP, UPDATE THE YH ARRAY. 341C CONSIDER CHANGING H IF IALTH = 1. OTHERWISE DECREASE IALTH BY 1. 342C IF IALTH IS THEN 1 AND NQ .LT. MAXORD, THEN ACOR IS SAVED FOR 343C USE IN A POSSIBLE ORDER INCREASE ON THE NEXT STEP. 344C IF A CHANGE IN H IS CONSIDERED, AN INCREASE OR DECREASE IN ORDER 345C BY ONE IS CONSIDERED ALSO. A CHANGE IN H IS MADE ONLY IF IT IS BY A 346C FACTOR OF AT LEAST 1.1. IF NOT, IALTH IS SET TO 3 TO PREVENT 347C TESTING FOR THAT MANY STEPS. 348C----------------------------------------------------------------------- 349 KFLAG = 0 350 IREDO = 0 351 NST = NST + 1 352 HU = H 353 NQU = NQ 354 DO 470 J = 1,L 355 DO 470 I = 1,N 356 470 YH(I,J) = YH(I,J) + EL(J)*ACOR(I) 357 IALTH = IALTH - 1 358 IF (IALTH .EQ. 0) GO TO 520 359 IF (IALTH .GT. 1) GO TO 700 360 IF (L .EQ. LMAX) GO TO 700 361 DO 490 I = 1,N 362 490 YH(I,LMAX) = ACOR(I) 363 GO TO 700 364C----------------------------------------------------------------------- 365C THE ERROR TEST FAILED. KFLAG KEEPS TRACK OF MULTIPLE FAILURES. 366C RESTORE TN AND THE YH ARRAY TO THEIR PREVIOUS VALUES, AND PREPARE 367C TO TRY THE STEP AGAIN. COMPUTE THE OPTIMUM STEP SIZE FOR THIS OR 368C ONE LOWER ORDER. AFTER 2 OR MORE FAILURES, H IS FORCED TO DECREASE 369C BY A FACTOR OF 0.2 OR LESS. 370C----------------------------------------------------------------------- 371 500 KFLAG = KFLAG - 1 372 TN = TOLD 373 I1 = NQNYH + 1 374 DO 515 JB = 1,NQ 375 I1 = I1 - NYH 376CDIR$ IVDEP 377 DO 510 I = I1,NQNYH 378 510 YH1(I) = YH1(I) - YH1(I+NYH) 379 515 CONTINUE 380 RMAX = 2.0D0 381 IF (DABS(H) .LE. HMIN*1.00001D0) GO TO 660 382 IF (KFLAG .LE. -3) GO TO 640 383 IREDO = 2 384 RHUP = 0.0D0 385 GO TO 540 386C----------------------------------------------------------------------- 387C REGARDLESS OF THE SUCCESS OR FAILURE OF THE STEP, FACTORS 388C RHDN, RHSM, AND RHUP ARE COMPUTED, BY WHICH H COULD BE MULTIPLIED 389C AT ORDER NQ - 1, ORDER NQ, OR ORDER NQ + 1, RESPECTIVELY. 390C IN THE CASE OF FAILURE, RHUP = 0.0 TO AVOID AN ORDER INCREASE. 391C THE LARGEST OF THESE IS DETERMINED AND THE NEW ORDER CHOSEN 392C ACCORDINGLY. IF THE ORDER IS TO BE INCREASED, WE COMPUTE ONE 393C ADDITIONAL SCALED DERIVATIVE. 394C----------------------------------------------------------------------- 395 520 RHUP = 0.0D0 396 IF (L .EQ. LMAX) GO TO 540 397 DO 530 I = 1,N 398 530 SAVF(I) = ACOR(I) - YH(I,LMAX) 399 DUP = VNORM (N, SAVF, EWT)/TESCO(3,NQ) 400 EXUP = 1.0D0/DBLE(L+1) 401 RHUP = 1.0D0/(1.4D0*DUP**EXUP + 0.0000014D0) 402 540 EXSM = 1.0D0/DBLE(L) 403 RHSM = 1.0D0/(1.2D0*DSM**EXSM + 0.0000012D0) 404 RHDN = 0.0D0 405 IF (NQ .EQ. 1) GO TO 560 406 DDN = VNORM (N, YH(1,L), EWT)/TESCO(1,NQ) 407 EXDN = 1.0D0/DBLE(NQ) 408 RHDN = 1.0D0/(1.3D0*DDN**EXDN + 0.0000013D0) 409 560 IF (RHSM .GE. RHUP) GO TO 570 410 IF (RHUP .GT. RHDN) GO TO 590 411 GO TO 580 412 570 IF (RHSM .LT. RHDN) GO TO 580 413 NEWQ = NQ 414 RH = RHSM 415 GO TO 620 416 580 NEWQ = NQ - 1 417 RH = RHDN 418 IF (KFLAG .LT. 0 .AND. RH .GT. 1.0D0) RH = 1.0D0 419 GO TO 620 420 590 NEWQ = L 421 RH = RHUP 422 IF (RH .LT. 1.1D0) GO TO 610 423 R = EL(L)/DBLE(L) 424 DO 600 I = 1,N 425 600 YH(I,NEWQ+1) = ACOR(I)*R 426 GO TO 630 427 610 IALTH = 3 428 GO TO 700 429 620 IF ((KFLAG .EQ. 0) .AND. (RH .LT. 1.1D0)) GO TO 610 430 IF (KFLAG .LE. -2) RH = DMIN1(RH,0.2D0) 431C----------------------------------------------------------------------- 432C IF THERE IS A CHANGE OF ORDER, RESET NQ, L, AND THE COEFFICIENTS. 433C IN ANY CASE H IS RESET ACCORDING TO RH AND THE YH ARRAY IS RESCALED. 434C THEN EXIT FROM 690 IF THE STEP WAS OK, OR REDO THE STEP OTHERWISE. 435C----------------------------------------------------------------------- 436 IF (NEWQ .EQ. NQ) GO TO 170 437 630 NQ = NEWQ 438 L = NQ + 1 439 IRET = 2 440 GO TO 150 441C----------------------------------------------------------------------- 442C CONTROL REACHES THIS SECTION IF 3 OR MORE FAILURES HAVE OCCURRED. 443C IF 10 FAILURES HAVE OCCURRED, EXIT WITH KFLAG = -1. 444C IT IS ASSUMED THAT THE DERIVATIVES THAT HAVE ACCUMULATED IN THE 445C YH ARRAY HAVE ERRORS OF THE WRONG ORDER. HENCE THE FIRST 446C DERIVATIVE IS RECOMPUTED, AND THE ORDER IS SET TO 1. THEN 447C H IS REDUCED BY A FACTOR OF 10, AND THE STEP IS RETRIED, 448C UNTIL IT SUCCEEDS OR H REACHES HMIN. 449C----------------------------------------------------------------------- 450 640 IF (KFLAG .EQ. -10) GO TO 660 451 RH = 0.1D0 452 RH = DMAX1(HMIN/DABS(H),RH) 453 H = H*RH 454 DO 645 I = 1,N 455 645 Y(I) = YH(I,1) 456 IERR = 0 457 CALL F (NEQ, TN, Y, SAVF, IERR) 458 IF (IERR .LT. 0) RETURN 459 NFE = NFE + 1 460 DO 650 I = 1,N 461 650 YH(I,2) = H*SAVF(I) 462 IPUP = MITER 463 IALTH = 5 464 IF (NQ .EQ. 1) GO TO 200 465 NQ = 1 466 L = 2 467 IRET = 3 468 GO TO 150 469C----------------------------------------------------------------------- 470C ALL RETURNS ARE MADE THROUGH THIS SECTION. H IS SAVED IN HOLD 471C TO ALLOW THE CALLER TO CHANGE H ON THE NEXT STEP. 472C----------------------------------------------------------------------- 473 660 KFLAG = -1 474 GO TO 720 475 670 KFLAG = -2 476 GO TO 720 477 680 KFLAG = -3 478 GO TO 720 479 690 RMAX = 10.0D0 480 700 R = 1.0D0/TESCO(2,NQU) 481 DO 710 I = 1,N 482 710 ACOR(I) = ACOR(I)*R 483 720 HOLD = H 484 JSTART = 1 485 RETURN 486C----------------------- END OF SUBROUTINE STODE ----------------------- 487 END 488