1*DECK DSTOD 2 SUBROUTINE DSTOD (NEQ, Y, YH, NYH, YH1, EWT, SAVF, ACOR, WM, IWM, 3 + DF, DJAC, RPAR, IPAR) 4C***BEGIN PROLOGUE DSTOD 5C***SUBSIDIARY 6C***PURPOSE Subsidiary to DDEBDF 7C***LIBRARY SLATEC 8C***TYPE DOUBLE PRECISION (STOD-S, DSTOD-D) 9C***AUTHOR Watts, H. A., (SNLA) 10C***DESCRIPTION 11C 12C DSTOD integrates a system of first order odes over one step in the 13C integrator package DDEBDF. 14C ---------------------------------------------------------------------- 15C DSTOD performs one step of the integration of an initial value 16C problem for a system of ordinary differential equations. 17C Note.. DSTOD 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 DSTOD 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 DF and DJAC. 24C NEQ = Integer array containing problem size in NEQ(1), and 25C passed as the NEQ argument in all calls to DF and DJAC. 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 = DOUBLE PRECISION and INTEGER work arrays associated with 41C matrix operations in chord iteration (MITER .NE. 0). 42C DPJAC = Name of routine to evaluate and preprocess Jacobian matrix 43C if a chord method is being used. 44C DSLVS = 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 DDEBDF 79C***ROUTINES CALLED DCFOD, DPJAC, DSLVS, DVNRMS 80C***COMMON BLOCKS DDEBD1 81C***REVISION HISTORY (YYMMDD) 82C 820301 DATE WRITTEN 83C 890531 Changed all specific intrinsics to generic. (WRB) 84C 890911 Removed unnecessary intrinsics. (WRB) 85C 891214 Prologue converted to Version 4.0 format. (BAB) 86C 900328 Added TYPE section. (WRB) 87C 910722 Updated AUTHOR section. (ALS) 88C 920422 Changed DIMENSION statement. (WRB) 89C***END PROLOGUE DSTOD 90C 91 INTEGER I, I1, IALTH, IER, IOD, IOWND, IPAR, IPUP, IREDO, IRET, 92 1 IWM, J, JB, JSTART, KFLAG, KSTEPS, L, LMAX, M, MAXORD, 93 2 MEO, METH, MITER, N, NCF, NEQ, NEWQ, NFE, NJE, NQ, NQNYH, 94 3 NQU, NST, NSTEPJ, NYH 95 DOUBLE PRECISION ACOR, CONIT, CRATE, DCON, DDN, 96 1 DEL, DELP, DSM, DUP, DVNRMS, EL, EL0, ELCO, 97 2 EWT, EXDN, EXSM, EXUP, H, HMIN, HMXI, HOLD, HU, R, RC, 98 3 RH, RHDN, RHSM, RHUP, RMAX, ROWND, RPAR, SAVF, TESCO, 99 4 TN, TOLD, UROUND, WM, Y, YH, YH1 100 EXTERNAL DF, DJAC 101C 102 DIMENSION Y(*),YH(NYH,*),YH1(*),EWT(*),SAVF(*),ACOR(*),WM(*), 103 1 IWM(*),RPAR(*),IPAR(*) 104 COMMON /DDEBD1/ ROWND,CONIT,CRATE,EL(13),ELCO(13,12),HOLD,RC,RMAX, 105 1 TESCO(3,12),EL0,H,HMIN,HMXI,HU,TN,UROUND,IOWND(7), 106 2 KSTEPS,IOD(6),IALTH,IPUP,LMAX,MEO,NQNYH,NSTEPJ, 107 3 IER,JSTART,KFLAG,L,METH,MITER,MAXORD,N,NQ,NST,NFE, 108 4 NJE,NQU 109C 110C 111C BEGIN BLOCK PERMITTING ...EXITS TO 690 112C BEGIN BLOCK PERMITTING ...EXITS TO 60 113C***FIRST EXECUTABLE STATEMENT DSTOD 114 KFLAG = 0 115 TOLD = TN 116 NCF = 0 117 IF (JSTART .GT. 0) GO TO 160 118 IF (JSTART .EQ. -1) GO TO 10 119 IF (JSTART .EQ. -2) GO TO 90 120C --------------------------------------------------------- 121C ON THE FIRST CALL, THE ORDER IS SET TO 1, AND OTHER 122C VARIABLES ARE INITIALIZED. RMAX IS THE MAXIMUM RATIO BY 123C WHICH H CAN BE INCREASED IN A SINGLE STEP. IT IS 124C INITIALLY 1.E4 TO COMPENSATE FOR THE SMALL INITIAL H, 125C BUT THEN IS NORMALLY EQUAL TO 10. IF A FAILURE OCCURS 126C (IN CORRECTOR CONVERGENCE OR ERROR TEST), RMAX IS SET AT 127C 2 FOR THE NEXT INCREASE. 128C --------------------------------------------------------- 129 LMAX = MAXORD + 1 130 NQ = 1 131 L = 2 132 IALTH = 2 133 RMAX = 10000.0D0 134 RC = 0.0D0 135 EL0 = 1.0D0 136 CRATE = 0.7D0 137 DELP = 0.0D0 138 HOLD = H 139 MEO = METH 140 NSTEPJ = 0 141 IRET = 3 142 GO TO 50 143 10 CONTINUE 144C BEGIN BLOCK PERMITTING ...EXITS TO 30 145C ------------------------------------------------------ 146C THE FOLLOWING BLOCK HANDLES PRELIMINARIES NEEDED WHEN 147C JSTART = -1. IPUP IS SET TO MITER TO FORCE A MATRIX 148C UPDATE. IF AN ORDER INCREASE IS ABOUT TO BE 149C CONSIDERED (IALTH = 1), IALTH IS RESET TO 2 TO 150C POSTPONE CONSIDERATION ONE MORE STEP. IF THE CALLER 151C HAS CHANGED METH, DCFOD IS CALLED TO RESET THE 152C COEFFICIENTS OF THE METHOD. IF THE CALLER HAS 153C CHANGED MAXORD TO A VALUE LESS THAN THE CURRENT 154C ORDER NQ, NQ IS REDUCED TO MAXORD, AND A NEW H CHOSEN 155C ACCORDINGLY. IF H IS TO BE CHANGED, YH MUST BE 156C RESCALED. IF H OR METH IS BEING CHANGED, IALTH IS 157C RESET TO L = NQ + 1 TO PREVENT FURTHER CHANGES IN H 158C FOR THAT MANY STEPS. 159C ------------------------------------------------------ 160 IPUP = MITER 161 LMAX = MAXORD + 1 162 IF (IALTH .EQ. 1) IALTH = 2 163 IF (METH .EQ. MEO) GO TO 20 164 CALL DCFOD(METH,ELCO,TESCO) 165 MEO = METH 166C ......EXIT 167 IF (NQ .GT. MAXORD) GO TO 30 168 IALTH = L 169 IRET = 1 170C ............EXIT 171 GO TO 60 172 20 CONTINUE 173 IF (NQ .LE. MAXORD) GO TO 90 174 30 CONTINUE 175 NQ = MAXORD 176 L = LMAX 177 DO 40 I = 1, L 178 EL(I) = ELCO(I,NQ) 179 40 CONTINUE 180 NQNYH = NQ*NYH 181 RC = RC*EL(1)/EL0 182 EL0 = EL(1) 183 CONIT = 0.5D0/(NQ+2) 184 DDN = DVNRMS(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 660 190 RH = MIN(RH,ABS(H/HOLD)) 191 H = HOLD 192 GO TO 100 193 50 CONTINUE 194C ------------------------------------------------------------ 195C DCFOD IS CALLED TO GET ALL THE INTEGRATION COEFFICIENTS 196C FOR THE CURRENT METH. THEN THE EL VECTOR AND RELATED 197C CONSTANTS ARE RESET WHENEVER THE ORDER NQ IS CHANGED, OR AT 198C THE START OF THE PROBLEM. 199C ------------------------------------------------------------ 200 CALL DCFOD(METH,ELCO,TESCO) 201 60 CONTINUE 202 70 CONTINUE 203C BEGIN BLOCK PERMITTING ...EXITS TO 680 204 DO 80 I = 1, L 205 EL(I) = ELCO(I,NQ) 206 80 CONTINUE 207 NQNYH = NQ*NYH 208 RC = RC*EL(1)/EL0 209 EL0 = EL(1) 210 CONIT = 0.5D0/(NQ+2) 211 GO TO (90,660,160), IRET 212C --------------------------------------------------------- 213C IF H IS BEING CHANGED, THE H RATIO RH IS CHECKED AGAINST 214C RMAX, HMIN, AND HMXI, AND THE YH ARRAY RESCALED. IALTH 215C IS SET TO L = NQ + 1 TO PREVENT A CHANGE OF H FOR THAT 216C MANY STEPS, UNLESS FORCED BY A CONVERGENCE OR ERROR TEST 217C FAILURE. 218C --------------------------------------------------------- 219 90 CONTINUE 220 IF (H .EQ. HOLD) GO TO 160 221 RH = H/HOLD 222 H = HOLD 223 IREDO = 3 224 100 CONTINUE 225 110 CONTINUE 226 RH = MIN(RH,RMAX) 227 RH = RH/MAX(1.0D0,ABS(H)*HMXI*RH) 228 R = 1.0D0 229 DO 130 J = 2, L 230 R = R*RH 231 DO 120 I = 1, N 232 YH(I,J) = YH(I,J)*R 233 120 CONTINUE 234 130 CONTINUE 235 H = H*RH 236 RC = RC*RH 237 IALTH = L 238 IF (IREDO .NE. 0) GO TO 150 239 RMAX = 10.0D0 240 R = 1.0D0/TESCO(2,NQU) 241 DO 140 I = 1, N 242 ACOR(I) = ACOR(I)*R 243 140 CONTINUE 244C ...............EXIT 245 GO TO 690 246 150 CONTINUE 247C ------------------------------------------------------ 248C THIS SECTION COMPUTES THE PREDICTED VALUES BY 249C EFFECTIVELY MULTIPLYING THE YH ARRAY BY THE PASCAL 250C TRIANGLE MATRIX. RC IS THE RATIO OF NEW TO OLD 251C VALUES OF THE COEFFICIENT H*EL(1). WHEN RC DIFFERS 252C FROM 1 BY MORE THAN 30 PERCENT, IPUP IS SET TO MITER 253C TO FORCE DPJAC TO BE CALLED, IF A JACOBIAN IS 254C INVOLVED. IN ANY CASE, DPJAC IS CALLED AT LEAST 255C EVERY 20-TH STEP. 256C ------------------------------------------------------ 257 160 CONTINUE 258 170 CONTINUE 259C BEGIN BLOCK PERMITTING ...EXITS TO 610 260C BEGIN BLOCK PERMITTING ...EXITS TO 490 261 IF (ABS(RC-1.0D0) .GT. 0.3D0) IPUP = MITER 262 IF (NST .GE. NSTEPJ + 20) IPUP = MITER 263 TN = TN + H 264 I1 = NQNYH + 1 265 DO 190 JB = 1, NQ 266 I1 = I1 - NYH 267 DO 180 I = I1, NQNYH 268 YH1(I) = YH1(I) + YH1(I+NYH) 269 180 CONTINUE 270 190 CONTINUE 271 KSTEPS = KSTEPS + 1 272C --------------------------------------------- 273C UP TO 3 CORRECTOR ITERATIONS ARE TAKEN. A 274C CONVERGENCE TEST IS MADE ON THE R.M.S. NORM 275C OF EACH CORRECTION, WEIGHTED BY THE ERROR 276C WEIGHT VECTOR EWT. THE SUM OF THE 277C CORRECTIONS IS ACCUMULATED IN THE VECTOR 278C ACOR(I). THE YH ARRAY IS NOT ALTERED IN THE 279C CORRECTOR LOOP. 280C --------------------------------------------- 281 200 CONTINUE 282 M = 0 283 DO 210 I = 1, N 284 Y(I) = YH(I,1) 285 210 CONTINUE 286 CALL DF(TN,Y,SAVF,RPAR,IPAR) 287 NFE = NFE + 1 288 IF (IPUP .LE. 0) GO TO 220 289C --------------------------------------- 290C IF INDICATED, THE MATRIX P = I - 291C H*EL(1)*J IS REEVALUATED AND 292C PREPROCESSED BEFORE STARTING THE 293C CORRECTOR ITERATION. IPUP IS SET TO 0 294C AS AN INDICATOR THAT THIS HAS BEEN 295C DONE. 296C --------------------------------------- 297 IPUP = 0 298 RC = 1.0D0 299 NSTEPJ = NST 300 CRATE = 0.7D0 301 CALL DPJAC(NEQ,Y,YH,NYH,EWT,ACOR,SAVF, 302 1 WM,IWM,DF,DJAC,RPAR,IPAR) 303C ......EXIT 304 IF (IER .NE. 0) GO TO 440 305 220 CONTINUE 306 DO 230 I = 1, N 307 ACOR(I) = 0.0D0 308 230 CONTINUE 309 240 CONTINUE 310 IF (MITER .NE. 0) GO TO 270 311C ------------------------------------ 312C IN THE CASE OF FUNCTIONAL 313C ITERATION, UPDATE Y DIRECTLY FROM 314C THE RESULT OF THE LAST FUNCTION 315C EVALUATION. 316C ------------------------------------ 317 DO 250 I = 1, N 318 SAVF(I) = H*SAVF(I) - YH(I,2) 319 Y(I) = SAVF(I) - ACOR(I) 320 250 CONTINUE 321 DEL = DVNRMS(N,Y,EWT) 322 DO 260 I = 1, N 323 Y(I) = YH(I,1) + EL(1)*SAVF(I) 324 ACOR(I) = SAVF(I) 325 260 CONTINUE 326 GO TO 300 327 270 CONTINUE 328C ------------------------------------ 329C IN THE CASE OF THE CHORD METHOD, 330C COMPUTE THE CORRECTOR ERROR, AND 331C SOLVE THE LINEAR SYSTEM WITH THAT 332C AS RIGHT-HAND SIDE AND P AS 333C COEFFICIENT MATRIX. 334C ------------------------------------ 335 DO 280 I = 1, N 336 Y(I) = H*SAVF(I) 337 1 - (YH(I,2) + ACOR(I)) 338 280 CONTINUE 339 CALL DSLVS(WM,IWM,Y,SAVF) 340C ......EXIT 341 IF (IER .NE. 0) GO TO 430 342 DEL = DVNRMS(N,Y,EWT) 343 DO 290 I = 1, N 344 ACOR(I) = ACOR(I) + Y(I) 345 Y(I) = YH(I,1) + EL(1)*ACOR(I) 346 290 CONTINUE 347 300 CONTINUE 348C --------------------------------------- 349C TEST FOR CONVERGENCE. IF M.GT.0, AN 350C ESTIMATE OF THE CONVERGENCE RATE 351C CONSTANT IS STORED IN CRATE, AND THIS 352C IS USED IN THE TEST. 353C --------------------------------------- 354 IF (M .NE. 0) 355 1 CRATE = MAX(0.2D0*CRATE,DEL/DELP) 356 DCON = DEL*MIN(1.0D0,1.5D0*CRATE) 357 1 /(TESCO(2,NQ)*CONIT) 358 IF (DCON .GT. 1.0D0) GO TO 420 359C ------------------------------------ 360C THE CORRECTOR HAS CONVERGED. IPUP 361C IS SET TO -1 IF MITER .NE. 0, TO 362C SIGNAL THAT THE JACOBIAN INVOLVED 363C MAY NEED UPDATING LATER. THE LOCAL 364C ERROR TEST IS MADE AND CONTROL 365C PASSES TO STATEMENT 500 IF IT 366C FAILS. 367C ------------------------------------ 368 IF (MITER .NE. 0) IPUP = -1 369 IF (M .EQ. 0) DSM = DEL/TESCO(2,NQ) 370 IF (M .GT. 0) 371 1 DSM = DVNRMS(N,ACOR,EWT) 372 2 /TESCO(2,NQ) 373 IF (DSM .GT. 1.0D0) GO TO 380 374C BEGIN BLOCK 375C PERMITTING ...EXITS TO 360 376C ------------------------------ 377C AFTER A SUCCESSFUL STEP, 378C UPDATE THE YH ARRAY. 379C CONSIDER CHANGING H IF IALTH 380C = 1. OTHERWISE DECREASE 381C IALTH BY 1. IF IALTH IS THEN 382C 1 AND NQ .LT. MAXORD, THEN 383C ACOR IS SAVED FOR USE IN A 384C POSSIBLE ORDER INCREASE ON 385C THE NEXT STEP. IF A CHANGE 386C IN H IS CONSIDERED, AN 387C INCREASE OR DECREASE IN ORDER 388C BY ONE IS CONSIDERED ALSO. A 389C CHANGE IN H IS MADE ONLY IF 390C IT IS BY A FACTOR OF AT LEAST 391C 1.1. IF NOT, IALTH IS SET TO 392C 3 TO PREVENT TESTING FOR THAT 393C MANY STEPS. 394C ------------------------------ 395 KFLAG = 0 396 IREDO = 0 397 NST = NST + 1 398 HU = H 399 NQU = NQ 400 DO 320 J = 1, L 401 DO 310 I = 1, N 402 YH(I,J) = YH(I,J) 403 1 + EL(J) 404 2 *ACOR(I) 405 310 CONTINUE 406 320 CONTINUE 407 IALTH = IALTH - 1 408 IF (IALTH .NE. 0) GO TO 340 409C --------------------------- 410C REGARDLESS OF THE SUCCESS 411C OR FAILURE OF THE STEP, 412C FACTORS RHDN, RHSM, AND 413C RHUP ARE COMPUTED, BY 414C WHICH H COULD BE 415C MULTIPLIED AT ORDER NQ - 416C 1, ORDER NQ, OR ORDER NQ + 417C 1, RESPECTIVELY. IN THE 418C CASE OF FAILURE, RHUP = 419C 0.0 TO AVOID AN ORDER 420C INCREASE. THE LARGEST OF 421C THESE IS DETERMINED AND 422C THE NEW ORDER CHOSEN 423C ACCORDINGLY. IF THE ORDER 424C IS TO BE INCREASED, WE 425C COMPUTE ONE ADDITIONAL 426C SCALED DERIVATIVE. 427C --------------------------- 428 RHUP = 0.0D0 429C .....................EXIT 430 IF (L .EQ. LMAX) GO TO 490 431 DO 330 I = 1, N 432 SAVF(I) = ACOR(I) 433 1 - YH(I,LMAX) 434 330 CONTINUE 435 DUP = DVNRMS(N,SAVF,EWT) 436 1 /TESCO(3,NQ) 437 EXUP = 1.0D0/(L+1) 438 RHUP = 1.0D0 439 1 /(1.4D0*DUP**EXUP 440 2 + 0.0000014D0) 441C .....................EXIT 442 GO TO 490 443 340 CONTINUE 444C ...EXIT 445 IF (IALTH .GT. 1) GO TO 360 446C ...EXIT 447 IF (L .EQ. LMAX) GO TO 360 448 DO 350 I = 1, N 449 YH(I,LMAX) = ACOR(I) 450 350 CONTINUE 451 360 CONTINUE 452 R = 1.0D0/TESCO(2,NQU) 453 DO 370 I = 1, N 454 ACOR(I) = ACOR(I)*R 455 370 CONTINUE 456C .................................EXIT 457 GO TO 690 458 380 CONTINUE 459C ------------------------------------ 460C THE ERROR TEST FAILED. KFLAG KEEPS 461C TRACK OF MULTIPLE FAILURES. 462C RESTORE TN AND THE YH ARRAY TO 463C THEIR PREVIOUS VALUES, AND PREPARE 464C TO TRY THE STEP AGAIN. COMPUTE THE 465C OPTIMUM STEP SIZE FOR THIS OR ONE 466C LOWER ORDER. AFTER 2 OR MORE 467C FAILURES, H IS FORCED TO DECREASE 468C BY A FACTOR OF 0.2 OR LESS. 469C ------------------------------------ 470 KFLAG = KFLAG - 1 471 TN = TOLD 472 I1 = NQNYH + 1 473 DO 400 JB = 1, NQ 474 I1 = I1 - NYH 475 DO 390 I = I1, NQNYH 476 YH1(I) = YH1(I) - YH1(I+NYH) 477 390 CONTINUE 478 400 CONTINUE 479 RMAX = 2.0D0 480 IF (ABS(H) .GT. HMIN*1.00001D0) 481 1 GO TO 410 482C --------------------------------- 483C ALL RETURNS ARE MADE THROUGH 484C THIS SECTION. H IS SAVED IN 485C HOLD TO ALLOW THE CALLER TO 486C CHANGE H ON THE NEXT STEP. 487C --------------------------------- 488 KFLAG = -1 489C .................................EXIT 490 GO TO 690 491 410 CONTINUE 492C ...............EXIT 493 IF (KFLAG .LE. -3) GO TO 610 494 IREDO = 2 495 RHUP = 0.0D0 496C ............EXIT 497 GO TO 490 498 420 CONTINUE 499 M = M + 1 500C ...EXIT 501 IF (M .EQ. 3) GO TO 430 502C ...EXIT 503 IF (M .GE. 2 .AND. DEL .GT. 2.0D0*DELP) 504 1 GO TO 430 505 DELP = DEL 506 CALL DF(TN,Y,SAVF,RPAR,IPAR) 507 NFE = NFE + 1 508 GO TO 240 509 430 CONTINUE 510C ------------------------------------------ 511C THE CORRECTOR ITERATION FAILED TO 512C CONVERGE IN 3 TRIES. IF MITER .NE. 0 AND 513C THE JACOBIAN IS OUT OF DATE, DPJAC IS 514C CALLED FOR THE NEXT TRY. OTHERWISE THE 515C YH ARRAY IS RETRACTED TO ITS VALUES 516C BEFORE PREDICTION, AND H IS REDUCED, IF 517C POSSIBLE. IF H CANNOT BE REDUCED OR 10 518C FAILURES HAVE OCCURRED, EXIT WITH KFLAG = 519C -2. 520C ------------------------------------------ 521C ...EXIT 522 IF (IPUP .EQ. 0) GO TO 440 523 IPUP = MITER 524 GO TO 200 525 440 CONTINUE 526 TN = TOLD 527 NCF = NCF + 1 528 RMAX = 2.0D0 529 I1 = NQNYH + 1 530 DO 460 JB = 1, NQ 531 I1 = I1 - NYH 532 DO 450 I = I1, NQNYH 533 YH1(I) = YH1(I) - YH1(I+NYH) 534 450 CONTINUE 535 460 CONTINUE 536 IF (ABS(H) .GT. HMIN*1.00001D0) GO TO 470 537 KFLAG = -2 538C ........................EXIT 539 GO TO 690 540 470 CONTINUE 541 IF (NCF .NE. 10) GO TO 480 542 KFLAG = -2 543C ........................EXIT 544 GO TO 690 545 480 CONTINUE 546 RH = 0.25D0 547 IPUP = MITER 548 IREDO = 1 549C .........EXIT 550 GO TO 650 551 490 CONTINUE 552 EXSM = 1.0D0/L 553 RHSM = 1.0D0/(1.2D0*DSM**EXSM + 0.0000012D0) 554 RHDN = 0.0D0 555 IF (NQ .EQ. 1) GO TO 500 556 DDN = DVNRMS(N,YH(1,L),EWT)/TESCO(1,NQ) 557 EXDN = 1.0D0/NQ 558 RHDN = 1.0D0/(1.3D0*DDN**EXDN + 0.0000013D0) 559 500 CONTINUE 560 IF (RHSM .GE. RHUP) GO TO 550 561 IF (RHUP .LE. RHDN) GO TO 540 562 NEWQ = L 563 RH = RHUP 564 IF (RH .GE. 1.1D0) GO TO 520 565 IALTH = 3 566 R = 1.0D0/TESCO(2,NQU) 567 DO 510 I = 1, N 568 ACOR(I) = ACOR(I)*R 569 510 CONTINUE 570C ...........................EXIT 571 GO TO 690 572 520 CONTINUE 573 R = EL(L)/L 574 DO 530 I = 1, N 575 YH(I,NEWQ+1) = ACOR(I)*R 576 530 CONTINUE 577 NQ = NEWQ 578 L = NQ + 1 579 IRET = 2 580C ..................EXIT 581 GO TO 680 582 540 CONTINUE 583 GO TO 580 584 550 CONTINUE 585 IF (RHSM .LT. RHDN) GO TO 580 586 NEWQ = NQ 587 RH = RHSM 588 IF (KFLAG .EQ. 0 .AND. RH .LT. 1.1D0) 589 1 GO TO 560 590 IF (KFLAG .LE. -2) RH = MIN(RH,0.2D0) 591C ------------------------------------------ 592C IF THERE IS A CHANGE OF ORDER, RESET NQ, 593C L, AND THE COEFFICIENTS. IN ANY CASE H 594C IS RESET ACCORDING TO RH AND THE YH ARRAY 595C IS RESCALED. THEN EXIT FROM 680 IF THE 596C STEP WAS OK, OR REDO THE STEP OTHERWISE. 597C ------------------------------------------ 598C ............EXIT 599 IF (NEWQ .EQ. NQ) GO TO 650 600 NQ = NEWQ 601 L = NQ + 1 602 IRET = 2 603C ..................EXIT 604 GO TO 680 605 560 CONTINUE 606 IALTH = 3 607 R = 1.0D0/TESCO(2,NQU) 608 DO 570 I = 1, N 609 ACOR(I) = ACOR(I)*R 610 570 CONTINUE 611C .....................EXIT 612 GO TO 690 613 580 CONTINUE 614 NEWQ = NQ - 1 615 RH = RHDN 616 IF (KFLAG .LT. 0 .AND. RH .GT. 1.0D0) RH = 1.0D0 617 IF (KFLAG .EQ. 0 .AND. RH .LT. 1.1D0) GO TO 590 618 IF (KFLAG .LE. -2) RH = MIN(RH,0.2D0) 619C --------------------------------------------- 620C IF THERE IS A CHANGE OF ORDER, RESET NQ, L, 621C AND THE COEFFICIENTS. IN ANY CASE H IS 622C RESET ACCORDING TO RH AND THE YH ARRAY IS 623C RESCALED. THEN EXIT FROM 680 IF THE STEP 624C WAS OK, OR REDO THE STEP OTHERWISE. 625C --------------------------------------------- 626C .........EXIT 627 IF (NEWQ .EQ. NQ) GO TO 650 628 NQ = NEWQ 629 L = NQ + 1 630 IRET = 2 631C ...............EXIT 632 GO TO 680 633 590 CONTINUE 634 IALTH = 3 635 R = 1.0D0/TESCO(2,NQU) 636 DO 600 I = 1, N 637 ACOR(I) = ACOR(I)*R 638 600 CONTINUE 639C ..................EXIT 640 GO TO 690 641 610 CONTINUE 642C --------------------------------------------------- 643C CONTROL REACHES THIS SECTION IF 3 OR MORE FAILURES 644C HAVE OCCURRED. IF 10 FAILURES HAVE OCCURRED, EXIT 645C WITH KFLAG = -1. IT IS ASSUMED THAT THE 646C DERIVATIVES THAT HAVE ACCUMULATED IN THE YH ARRAY 647C HAVE ERRORS OF THE WRONG ORDER. HENCE THE FIRST 648C DERIVATIVE IS RECOMPUTED, AND THE ORDER IS SET TO 649C 1. THEN H IS REDUCED BY A FACTOR OF 10, AND THE 650C STEP IS RETRIED, UNTIL IT SUCCEEDS OR H REACHES 651C HMIN. 652C --------------------------------------------------- 653 IF (KFLAG .NE. -10) GO TO 620 654C ------------------------------------------------ 655C ALL RETURNS ARE MADE THROUGH THIS SECTION. H 656C IS SAVED IN HOLD TO ALLOW THE CALLER TO CHANGE 657C H ON THE NEXT STEP. 658C ------------------------------------------------ 659 KFLAG = -1 660C ..................EXIT 661 GO TO 690 662 620 CONTINUE 663 RH = 0.1D0 664 RH = MAX(HMIN/ABS(H),RH) 665 H = H*RH 666 DO 630 I = 1, N 667 Y(I) = YH(I,1) 668 630 CONTINUE 669 CALL DF(TN,Y,SAVF,RPAR,IPAR) 670 NFE = NFE + 1 671 DO 640 I = 1, N 672 YH(I,2) = H*SAVF(I) 673 640 CONTINUE 674 IPUP = MITER 675 IALTH = 5 676C ......EXIT 677 IF (NQ .NE. 1) GO TO 670 678 GO TO 170 679 650 CONTINUE 680 660 CONTINUE 681 RH = MAX(RH,HMIN/ABS(H)) 682 GO TO 110 683 670 CONTINUE 684 NQ = 1 685 L = 2 686 IRET = 3 687 680 CONTINUE 688 GO TO 70 689 690 CONTINUE 690 HOLD = H 691 JSTART = 1 692 RETURN 693C ----------------------- END OF SUBROUTINE DSTOD 694C ----------------------- 695 END 696