00001 SUBROUTINE STODE (NEQ, Y, YH, NYH, YH1, EWT, SAVF, ACOR, 00002 1 WM, IWM, F, JAC, PJAC, SLVS, IERR) 00003 CLLL. OPTIMIZE 00004 EXTERNAL F, JAC, PJAC, SLVS 00005 INTEGER NEQ, NYH, IWM 00006 INTEGER IOWND, IALTH, IPUP, LMAX, MEO, NQNYH, NSLP, 00007 1 ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L, METH, MITER, 00008 2 MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU 00009 INTEGER I, I1, IREDO, IRET, J, JB, M, NCF, NEWQ 00010 DOUBLE PRECISION Y, YH, YH1, EWT, SAVF, ACOR, WM 00011 DOUBLE PRECISION CONIT, CRATE, EL, ELCO, HOLD, RMAX, TESCO, 00012 2 CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND 00013 DOUBLE PRECISION DCON, DDN, DEL, DELP, DSM, DUP, EXDN, EXSM, EXUP, 00014 1 R, RH, RHDN, RHSM, RHUP, TOLD, VNORM 00015 DIMENSION NEQ(*), Y(*), YH(NYH,*), YH1(*), EWT(*), SAVF(*), 00016 1 ACOR(*), WM(*), IWM(*) 00017 COMMON /LS0001/ CONIT, CRATE, EL(13), ELCO(13,12), 00018 1 HOLD, RMAX, TESCO(3,12), 00019 2 CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND, IOWND(14), 00020 3 IALTH, IPUP, LMAX, MEO, NQNYH, NSLP, 00021 4 ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L, METH, MITER, 00022 5 MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU 00023 C----------------------------------------------------------------------- 00024 C STODE PERFORMS ONE STEP OF THE INTEGRATION OF AN INITIAL VALUE 00025 C PROBLEM FOR A SYSTEM OF ORDINARY DIFFERENTIAL EQUATIONS. 00026 C NOTE.. STODE IS INDEPENDENT OF THE VALUE OF THE ITERATION METHOD 00027 C INDICATOR MITER, WHEN THIS IS .NE. 0, AND HENCE IS INDEPENDENT 00028 C OF THE TYPE OF CHORD METHOD USED, OR THE JACOBIAN STRUCTURE. 00029 C COMMUNICATION WITH STODE IS DONE WITH THE FOLLOWING VARIABLES.. 00030 C 00031 C NEQ = INTEGER ARRAY CONTAINING PROBLEM SIZE IN NEQ(1), AND 00032 C PASSED AS THE NEQ ARGUMENT IN ALL CALLS TO F AND JAC. 00033 C Y = AN ARRAY OF LENGTH .GE. N USED AS THE Y ARGUMENT IN 00034 C ALL CALLS TO F AND JAC. 00035 C YH = AN NYH BY LMAX ARRAY CONTAINING THE DEPENDENT VARIABLES 00036 C AND THEIR APPROXIMATE SCALED DERIVATIVES, WHERE 00037 C LMAX = MAXORD + 1. YH(I,J+1) CONTAINS THE APPROXIMATE 00038 C J-TH DERIVATIVE OF Y(I), SCALED BY H**J/FACTORIAL(J) 00039 C (J = 0,1,...,NQ). ON ENTRY FOR THE FIRST STEP, THE FIRST 00040 C TWO COLUMNS OF YH MUST BE SET FROM THE INITIAL VALUES. 00041 C NYH = A CONSTANT INTEGER .GE. N, THE FIRST DIMENSION OF YH. 00042 C YH1 = A ONE-DIMENSIONAL ARRAY OCCUPYING THE SAME SPACE AS YH. 00043 C EWT = AN ARRAY OF LENGTH N CONTAINING MULTIPLICATIVE WEIGHTS 00044 C FOR LOCAL ERROR MEASUREMENTS. LOCAL ERRORS IN Y(I) ARE 00045 C COMPARED TO 1.0/EWT(I) IN VARIOUS ERROR TESTS. 00046 C SAVF = AN ARRAY OF WORKING STORAGE, OF LENGTH N. 00047 C ALSO USED FOR INPUT OF YH(*,MAXORD+2) WHEN JSTART = -1 00048 C AND MAXORD .LT. THE CURRENT ORDER NQ. 00049 C ACOR = A WORK ARRAY OF LENGTH N, USED FOR THE ACCUMULATED 00050 C CORRECTIONS. ON A SUCCESSFUL RETURN, ACOR(I) CONTAINS 00051 C THE ESTIMATED ONE-STEP LOCAL ERROR IN Y(I). 00052 C WM,IWM = REAL AND INTEGER WORK ARRAYS ASSOCIATED WITH MATRIX 00053 C OPERATIONS IN CHORD ITERATION (MITER .NE. 0). 00054 C PJAC = NAME OF ROUTINE TO EVALUATE AND PREPROCESS JACOBIAN MATRIX 00055 C AND P = I - H*EL0*JAC, IF A CHORD METHOD IS BEING USED. 00056 C SLVS = NAME OF ROUTINE TO SOLVE LINEAR SYSTEM IN CHORD ITERATION. 00057 C CCMAX = MAXIMUM RELATIVE CHANGE IN H*EL0 BEFORE PJAC IS CALLED. 00058 C H = THE STEP SIZE TO BE ATTEMPTED ON THE NEXT STEP. 00059 C H IS ALTERED BY THE ERROR CONTROL ALGORITHM DURING THE 00060 C PROBLEM. H CAN BE EITHER POSITIVE OR NEGATIVE, BUT ITS 00061 C SIGN MUST REMAIN CONSTANT THROUGHOUT THE PROBLEM. 00062 C HMIN = THE MINIMUM ABSOLUTE VALUE OF THE STEP SIZE H TO BE USED. 00063 C HMXI = INVERSE OF THE MAXIMUM ABSOLUTE VALUE OF H TO BE USED. 00064 C HMXI = 0.0 IS ALLOWED AND CORRESPONDS TO AN INFINITE HMAX. 00065 C HMIN AND HMXI MAY BE CHANGED AT ANY TIME, BUT WILL NOT 00066 C TAKE EFFECT UNTIL THE NEXT CHANGE OF H IS CONSIDERED. 00067 C TN = THE INDEPENDENT VARIABLE. TN IS UPDATED ON EACH STEP TAKEN. 00068 C JSTART = AN INTEGER USED FOR INPUT ONLY, WITH THE FOLLOWING 00069 C VALUES AND MEANINGS.. 00070 C 0 PERFORM THE FIRST STEP. 00071 C .GT.0 TAKE A NEW STEP CONTINUING FROM THE LAST. 00072 C -1 TAKE THE NEXT STEP WITH A NEW VALUE OF H, MAXORD, 00073 C N, METH, MITER, AND/OR MATRIX PARAMETERS. 00074 C -2 TAKE THE NEXT STEP WITH A NEW VALUE OF H, 00075 C BUT WITH OTHER INPUTS UNCHANGED. 00076 C ON RETURN, JSTART IS SET TO 1 TO FACILITATE CONTINUATION. 00077 C KFLAG = A COMPLETION CODE WITH THE FOLLOWING MEANINGS.. 00078 C 0 THE STEP WAS SUCCESFUL. 00079 C -1 THE REQUESTED ERROR COULD NOT BE ACHIEVED. 00080 C -2 CORRECTOR CONVERGENCE COULD NOT BE ACHIEVED. 00081 C -3 FATAL ERROR IN PJAC OR SLVS. 00082 C A RETURN WITH KFLAG = -1 OR -2 MEANS EITHER 00083 C ABS(H) = HMIN OR 10 CONSECUTIVE FAILURES OCCURRED. 00084 C ON A RETURN WITH KFLAG NEGATIVE, THE VALUES OF TN AND 00085 C THE YH ARRAY ARE AS OF THE BEGINNING OF THE LAST 00086 C STEP, AND H IS THE LAST STEP SIZE ATTEMPTED. 00087 C MAXORD = THE MAXIMUM ORDER OF INTEGRATION METHOD TO BE ALLOWED. 00088 C MAXCOR = THE MAXIMUM NUMBER OF CORRECTOR ITERATIONS ALLOWED. 00089 C MSBP = MAXIMUM NUMBER OF STEPS BETWEEN PJAC CALLS (MITER .GT. 0). 00090 C MXNCF = MAXIMUM NUMBER OF CONVERGENCE FAILURES ALLOWED. 00091 C METH/MITER = THE METHOD FLAGS. SEE DESCRIPTION IN DRIVER. 00092 C N = THE NUMBER OF FIRST-ORDER DIFFERENTIAL EQUATIONS. 00093 C IERR = ERROR FLAG FROM USER-SUPPLIED FUNCTION 00094 C----------------------------------------------------------------------- 00095 KFLAG = 0 00096 TOLD = TN 00097 NCF = 0 00098 IERPJ = 0 00099 IERSL = 0 00100 JCUR = 0 00101 ICF = 0 00102 DELP = 0.0D0 00103 IF (JSTART .GT. 0) GO TO 200 00104 IF (JSTART .EQ. -1) GO TO 100 00105 IF (JSTART .EQ. -2) GO TO 160 00106 C----------------------------------------------------------------------- 00107 C ON THE FIRST CALL, THE ORDER IS SET TO 1, AND OTHER VARIABLES ARE 00108 C INITIALIZED. RMAX IS THE MAXIMUM RATIO BY WHICH H CAN BE INCREASED 00109 C IN A SINGLE STEP. IT IS INITIALLY 1.E4 TO COMPENSATE FOR THE SMALL 00110 C INITIAL H, BUT THEN IS NORMALLY EQUAL TO 10. IF A FAILURE 00111 C OCCURS (IN CORRECTOR CONVERGENCE OR ERROR TEST), RMAX IS SET AT 2 00112 C FOR THE NEXT INCREASE. 00113 C----------------------------------------------------------------------- 00114 LMAX = MAXORD + 1 00115 NQ = 1 00116 L = 2 00117 IALTH = 2 00118 RMAX = 10000.0D0 00119 RC = 0.0D0 00120 EL0 = 1.0D0 00121 CRATE = 0.7D0 00122 HOLD = H 00123 MEO = METH 00124 NSLP = 0 00125 IPUP = MITER 00126 IRET = 3 00127 GO TO 140 00128 C----------------------------------------------------------------------- 00129 C THE FOLLOWING BLOCK HANDLES PRELIMINARIES NEEDED WHEN JSTART = -1. 00130 C IPUP IS SET TO MITER TO FORCE A MATRIX UPDATE. 00131 C IF AN ORDER INCREASE IS ABOUT TO BE CONSIDERED (IALTH = 1), 00132 C IALTH IS RESET TO 2 TO POSTPONE CONSIDERATION ONE MORE STEP. 00133 C IF THE CALLER HAS CHANGED METH, CFODE IS CALLED TO RESET 00134 C THE COEFFICIENTS OF THE METHOD. 00135 C IF THE CALLER HAS CHANGED MAXORD TO A VALUE LESS THAN THE CURRENT 00136 C ORDER NQ, NQ IS REDUCED TO MAXORD, AND A NEW H CHOSEN ACCORDINGLY. 00137 C IF H IS TO BE CHANGED, YH MUST BE RESCALED. 00138 C IF H OR METH IS BEING CHANGED, IALTH IS RESET TO L = NQ + 1 00139 C TO PREVENT FURTHER CHANGES IN H FOR THAT MANY STEPS. 00140 C----------------------------------------------------------------------- 00141 100 IPUP = MITER 00142 LMAX = MAXORD + 1 00143 IF (IALTH .EQ. 1) IALTH = 2 00144 IF (METH .EQ. MEO) GO TO 110 00145 CALL CFODE (METH, ELCO, TESCO) 00146 MEO = METH 00147 IF (NQ .GT. MAXORD) GO TO 120 00148 IALTH = L 00149 IRET = 1 00150 GO TO 150 00151 110 IF (NQ .LE. MAXORD) GO TO 160 00152 120 NQ = MAXORD 00153 L = LMAX 00154 DO 125 I = 1,L 00155 125 EL(I) = ELCO(I,NQ) 00156 NQNYH = NQ*NYH 00157 RC = RC*EL(1)/EL0 00158 EL0 = EL(1) 00159 CONIT = 0.5D0/DBLE(NQ+2) 00160 DDN = VNORM (N, SAVF, EWT)/TESCO(1,L) 00161 EXDN = 1.0D0/DBLE(L) 00162 RHDN = 1.0D0/(1.3D0*DDN**EXDN + 0.0000013D0) 00163 RH = DMIN1(RHDN,1.0D0) 00164 IREDO = 3 00165 IF (H .EQ. HOLD) GO TO 170 00166 RH = DMIN1(RH,DABS(H/HOLD)) 00167 H = HOLD 00168 GO TO 175 00169 C----------------------------------------------------------------------- 00170 C CFODE IS CALLED TO GET ALL THE INTEGRATION COEFFICIENTS FOR THE 00171 C CURRENT METH. THEN THE EL VECTOR AND RELATED CONSTANTS ARE RESET 00172 C WHENEVER THE ORDER NQ IS CHANGED, OR AT THE START OF THE PROBLEM. 00173 C----------------------------------------------------------------------- 00174 140 CALL CFODE (METH, ELCO, TESCO) 00175 150 DO 155 I = 1,L 00176 155 EL(I) = ELCO(I,NQ) 00177 NQNYH = NQ*NYH 00178 RC = RC*EL(1)/EL0 00179 EL0 = EL(1) 00180 CONIT = 0.5D0/DBLE(NQ+2) 00181 GO TO (160, 170, 200), IRET 00182 C----------------------------------------------------------------------- 00183 C IF H IS BEING CHANGED, THE H RATIO RH IS CHECKED AGAINST 00184 C RMAX, HMIN, AND HMXI, AND THE YH ARRAY RESCALED. IALTH IS SET TO 00185 C L = NQ + 1 TO PREVENT A CHANGE OF H FOR THAT MANY STEPS, UNLESS 00186 C FORCED BY A CONVERGENCE OR ERROR TEST FAILURE. 00187 C----------------------------------------------------------------------- 00188 160 IF (H .EQ. HOLD) GO TO 200 00189 RH = H/HOLD 00190 H = HOLD 00191 IREDO = 3 00192 GO TO 175 00193 170 RH = DMAX1(RH,HMIN/DABS(H)) 00194 175 RH = DMIN1(RH,RMAX) 00195 RH = RH/DMAX1(1.0D0,DABS(H)*HMXI*RH) 00196 R = 1.0D0 00197 DO 180 J = 2,L 00198 R = R*RH 00199 DO 180 I = 1,N 00200 180 YH(I,J) = YH(I,J)*R 00201 H = H*RH 00202 RC = RC*RH 00203 IALTH = L 00204 IF (IREDO .EQ. 0) GO TO 690 00205 C----------------------------------------------------------------------- 00206 C THIS SECTION COMPUTES THE PREDICTED VALUES BY EFFECTIVELY 00207 C MULTIPLYING THE YH ARRAY BY THE PASCAL TRIANGLE MATRIX. 00208 C RC IS THE RATIO OF NEW TO OLD VALUES OF THE COEFFICIENT H*EL(1). 00209 C WHEN RC DIFFERS FROM 1 BY MORE THAN CCMAX, IPUP IS SET TO MITER 00210 C TO FORCE PJAC TO BE CALLED, IF A JACOBIAN IS INVOLVED. 00211 C IN ANY CASE, PJAC IS CALLED AT LEAST EVERY MSBP STEPS. 00212 C----------------------------------------------------------------------- 00213 200 IF (DABS(RC-1.0D0) .GT. CCMAX) IPUP = MITER 00214 IF (NST .GE. NSLP+MSBP) IPUP = MITER 00215 TN = TN + H 00216 I1 = NQNYH + 1 00217 DO 215 JB = 1,NQ 00218 I1 = I1 - NYH 00219 CDIR$ IVDEP 00220 DO 210 I = I1,NQNYH 00221 210 YH1(I) = YH1(I) + YH1(I+NYH) 00222 215 CONTINUE 00223 C----------------------------------------------------------------------- 00224 C UP TO MAXCOR CORRECTOR ITERATIONS ARE TAKEN. A CONVERGENCE TEST IS 00225 C MADE ON THE R.M.S. NORM OF EACH CORRECTION, WEIGHTED BY THE ERROR 00226 C WEIGHT VECTOR EWT. THE SUM OF THE CORRECTIONS IS ACCUMULATED IN THE 00227 C VECTOR ACOR(I). THE YH ARRAY IS NOT ALTERED IN THE CORRECTOR LOOP. 00228 C----------------------------------------------------------------------- 00229 220 M = 0 00230 DO 230 I = 1,N 00231 230 Y(I) = YH(I,1) 00232 IERR = 0 00233 CALL F (NEQ, TN, Y, SAVF, IERR) 00234 IF (IERR .LT. 0) RETURN 00235 NFE = NFE + 1 00236 IF (IPUP .LE. 0) GO TO 250 00237 C----------------------------------------------------------------------- 00238 C IF INDICATED, THE MATRIX P = I - H*EL(1)*J IS REEVALUATED AND 00239 C PREPROCESSED BEFORE STARTING THE CORRECTOR ITERATION. IPUP IS SET 00240 C TO 0 AS AN INDICATOR THAT THIS HAS BEEN DONE. 00241 C----------------------------------------------------------------------- 00242 IERR = 0 00243 CALL PJAC (NEQ, Y, YH, NYH, EWT, ACOR, SAVF, WM, IWM, F, JAC, 00244 1 IERR) 00245 IF (IERR .LT. 0) RETURN 00246 IPUP = 0 00247 RC = 1.0D0 00248 NSLP = NST 00249 CRATE = 0.7D0 00250 IF (IERPJ .NE. 0) GO TO 430 00251 250 DO 260 I = 1,N 00252 260 ACOR(I) = 0.0D0 00253 270 IF (MITER .NE. 0) GO TO 350 00254 C----------------------------------------------------------------------- 00255 C IN THE CASE OF FUNCTIONAL ITERATION, UPDATE Y DIRECTLY FROM 00256 C THE RESULT OF THE LAST FUNCTION EVALUATION. 00257 C----------------------------------------------------------------------- 00258 DO 290 I = 1,N 00259 SAVF(I) = H*SAVF(I) - YH(I,2) 00260 290 Y(I) = SAVF(I) - ACOR(I) 00261 DEL = VNORM (N, Y, EWT) 00262 DO 300 I = 1,N 00263 Y(I) = YH(I,1) + EL(1)*SAVF(I) 00264 300 ACOR(I) = SAVF(I) 00265 GO TO 400 00266 C----------------------------------------------------------------------- 00267 C IN THE CASE OF THE CHORD METHOD, COMPUTE THE CORRECTOR ERROR, 00268 C AND SOLVE THE LINEAR SYSTEM WITH THAT AS RIGHT-HAND SIDE AND 00269 C P AS COEFFICIENT MATRIX. 00270 C----------------------------------------------------------------------- 00271 350 DO 360 I = 1,N 00272 360 Y(I) = H*SAVF(I) - (YH(I,2) + ACOR(I)) 00273 CALL SLVS (WM, IWM, Y, SAVF) 00274 IF (IERSL .LT. 0) GO TO 430 00275 IF (IERSL .GT. 0) GO TO 410 00276 DEL = VNORM (N, Y, EWT) 00277 DO 380 I = 1,N 00278 ACOR(I) = ACOR(I) + Y(I) 00279 380 Y(I) = YH(I,1) + EL(1)*ACOR(I) 00280 C----------------------------------------------------------------------- 00281 C TEST FOR CONVERGENCE. IF M.GT.0, AN ESTIMATE OF THE CONVERGENCE 00282 C RATE CONSTANT IS STORED IN CRATE, AND THIS IS USED IN THE TEST. 00283 C----------------------------------------------------------------------- 00284 400 IF (M .NE. 0) CRATE = DMAX1(0.2D0*CRATE,DEL/DELP) 00285 DCON = DEL*DMIN1(1.0D0,1.5D0*CRATE)/(TESCO(2,NQ)*CONIT) 00286 IF (DCON .LE. 1.0D0) GO TO 450 00287 M = M + 1 00288 IF (M .EQ. MAXCOR) GO TO 410 00289 IF (M .GE. 2 .AND. DEL .GT. 2.0D0*DELP) GO TO 410 00290 DELP = DEL 00291 IERR = 0 00292 CALL F (NEQ, TN, Y, SAVF, IERR) 00293 IF (IERR .LT. 0) RETURN 00294 NFE = NFE + 1 00295 GO TO 270 00296 C----------------------------------------------------------------------- 00297 C THE CORRECTOR ITERATION FAILED TO CONVERGE. 00298 C IF MITER .NE. 0 AND THE JACOBIAN IS OUT OF DATE, PJAC IS CALLED FOR 00299 C THE NEXT TRY. OTHERWISE THE YH ARRAY IS RETRACTED TO ITS VALUES 00300 C BEFORE PREDICTION, AND H IS REDUCED, IF POSSIBLE. IF H CANNOT BE 00301 C REDUCED OR MXNCF FAILURES HAVE OCCURRED, EXIT WITH KFLAG = -2. 00302 C----------------------------------------------------------------------- 00303 410 IF (MITER .EQ. 0 .OR. JCUR .EQ. 1) GO TO 430 00304 ICF = 1 00305 IPUP = MITER 00306 GO TO 220 00307 430 ICF = 2 00308 NCF = NCF + 1 00309 RMAX = 2.0D0 00310 TN = TOLD 00311 I1 = NQNYH + 1 00312 DO 445 JB = 1,NQ 00313 I1 = I1 - NYH 00314 CDIR$ IVDEP 00315 DO 440 I = I1,NQNYH 00316 440 YH1(I) = YH1(I) - YH1(I+NYH) 00317 445 CONTINUE 00318 IF (IERPJ .LT. 0 .OR. IERSL .LT. 0) GO TO 680 00319 IF (DABS(H) .LE. HMIN*1.00001D0) GO TO 670 00320 IF (NCF .EQ. MXNCF) GO TO 670 00321 RH = 0.25D0 00322 IPUP = MITER 00323 IREDO = 1 00324 GO TO 170 00325 C----------------------------------------------------------------------- 00326 C THE CORRECTOR HAS CONVERGED. JCUR IS SET TO 0 00327 C TO SIGNAL THAT THE JACOBIAN INVOLVED MAY NEED UPDATING LATER. 00328 C THE LOCAL ERROR TEST IS MADE AND CONTROL PASSES TO STATEMENT 500 00329 C IF IT FAILS. 00330 C----------------------------------------------------------------------- 00331 450 JCUR = 0 00332 IF (M .EQ. 0) DSM = DEL/TESCO(2,NQ) 00333 IF (M .GT. 0) DSM = VNORM (N, ACOR, EWT)/TESCO(2,NQ) 00334 IF (DSM .GT. 1.0D0) GO TO 500 00335 C----------------------------------------------------------------------- 00336 C AFTER A SUCCESSFUL STEP, UPDATE THE YH ARRAY. 00337 C CONSIDER CHANGING H IF IALTH = 1. OTHERWISE DECREASE IALTH BY 1. 00338 C IF IALTH IS THEN 1 AND NQ .LT. MAXORD, THEN ACOR IS SAVED FOR 00339 C USE IN A POSSIBLE ORDER INCREASE ON THE NEXT STEP. 00340 C IF A CHANGE IN H IS CONSIDERED, AN INCREASE OR DECREASE IN ORDER 00341 C BY ONE IS CONSIDERED ALSO. A CHANGE IN H IS MADE ONLY IF IT IS BY A 00342 C FACTOR OF AT LEAST 1.1. IF NOT, IALTH IS SET TO 3 TO PREVENT 00343 C TESTING FOR THAT MANY STEPS. 00344 C----------------------------------------------------------------------- 00345 KFLAG = 0 00346 IREDO = 0 00347 NST = NST + 1 00348 HU = H 00349 NQU = NQ 00350 DO 470 J = 1,L 00351 DO 470 I = 1,N 00352 470 YH(I,J) = YH(I,J) + EL(J)*ACOR(I) 00353 IALTH = IALTH - 1 00354 IF (IALTH .EQ. 0) GO TO 520 00355 IF (IALTH .GT. 1) GO TO 700 00356 IF (L .EQ. LMAX) GO TO 700 00357 DO 490 I = 1,N 00358 490 YH(I,LMAX) = ACOR(I) 00359 GO TO 700 00360 C----------------------------------------------------------------------- 00361 C THE ERROR TEST FAILED. KFLAG KEEPS TRACK OF MULTIPLE FAILURES. 00362 C RESTORE TN AND THE YH ARRAY TO THEIR PREVIOUS VALUES, AND PREPARE 00363 C TO TRY THE STEP AGAIN. COMPUTE THE OPTIMUM STEP SIZE FOR THIS OR 00364 C ONE LOWER ORDER. AFTER 2 OR MORE FAILURES, H IS FORCED TO DECREASE 00365 C BY A FACTOR OF 0.2 OR LESS. 00366 C----------------------------------------------------------------------- 00367 500 KFLAG = KFLAG - 1 00368 TN = TOLD 00369 I1 = NQNYH + 1 00370 DO 515 JB = 1,NQ 00371 I1 = I1 - NYH 00372 CDIR$ IVDEP 00373 DO 510 I = I1,NQNYH 00374 510 YH1(I) = YH1(I) - YH1(I+NYH) 00375 515 CONTINUE 00376 RMAX = 2.0D0 00377 IF (DABS(H) .LE. HMIN*1.00001D0) GO TO 660 00378 IF (KFLAG .LE. -3) GO TO 640 00379 IREDO = 2 00380 RHUP = 0.0D0 00381 GO TO 540 00382 C----------------------------------------------------------------------- 00383 C REGARDLESS OF THE SUCCESS OR FAILURE OF THE STEP, FACTORS 00384 C RHDN, RHSM, AND RHUP ARE COMPUTED, BY WHICH H COULD BE MULTIPLIED 00385 C AT ORDER NQ - 1, ORDER NQ, OR ORDER NQ + 1, RESPECTIVELY. 00386 C IN THE CASE OF FAILURE, RHUP = 0.0 TO AVOID AN ORDER INCREASE. 00387 C THE LARGEST OF THESE IS DETERMINED AND THE NEW ORDER CHOSEN 00388 C ACCORDINGLY. IF THE ORDER IS TO BE INCREASED, WE COMPUTE ONE 00389 C ADDITIONAL SCALED DERIVATIVE. 00390 C----------------------------------------------------------------------- 00391 520 RHUP = 0.0D0 00392 IF (L .EQ. LMAX) GO TO 540 00393 DO 530 I = 1,N 00394 530 SAVF(I) = ACOR(I) - YH(I,LMAX) 00395 DUP = VNORM (N, SAVF, EWT)/TESCO(3,NQ) 00396 EXUP = 1.0D0/DBLE(L+1) 00397 RHUP = 1.0D0/(1.4D0*DUP**EXUP + 0.0000014D0) 00398 540 EXSM = 1.0D0/DBLE(L) 00399 RHSM = 1.0D0/(1.2D0*DSM**EXSM + 0.0000012D0) 00400 RHDN = 0.0D0 00401 IF (NQ .EQ. 1) GO TO 560 00402 DDN = VNORM (N, YH(1,L), EWT)/TESCO(1,NQ) 00403 EXDN = 1.0D0/DBLE(NQ) 00404 RHDN = 1.0D0/(1.3D0*DDN**EXDN + 0.0000013D0) 00405 560 IF (RHSM .GE. RHUP) GO TO 570 00406 IF (RHUP .GT. RHDN) GO TO 590 00407 GO TO 580 00408 570 IF (RHSM .LT. RHDN) GO TO 580 00409 NEWQ = NQ 00410 RH = RHSM 00411 GO TO 620 00412 580 NEWQ = NQ - 1 00413 RH = RHDN 00414 IF (KFLAG .LT. 0 .AND. RH .GT. 1.0D0) RH = 1.0D0 00415 GO TO 620 00416 590 NEWQ = L 00417 RH = RHUP 00418 IF (RH .LT. 1.1D0) GO TO 610 00419 R = EL(L)/DBLE(L) 00420 DO 600 I = 1,N 00421 600 YH(I,NEWQ+1) = ACOR(I)*R 00422 GO TO 630 00423 610 IALTH = 3 00424 GO TO 700 00425 620 IF ((KFLAG .EQ. 0) .AND. (RH .LT. 1.1D0)) GO TO 610 00426 IF (KFLAG .LE. -2) RH = DMIN1(RH,0.2D0) 00427 C----------------------------------------------------------------------- 00428 C IF THERE IS A CHANGE OF ORDER, RESET NQ, L, AND THE COEFFICIENTS. 00429 C IN ANY CASE H IS RESET ACCORDING TO RH AND THE YH ARRAY IS RESCALED. 00430 C THEN EXIT FROM 690 IF THE STEP WAS OK, OR REDO THE STEP OTHERWISE. 00431 C----------------------------------------------------------------------- 00432 IF (NEWQ .EQ. NQ) GO TO 170 00433 630 NQ = NEWQ 00434 L = NQ + 1 00435 IRET = 2 00436 GO TO 150 00437 C----------------------------------------------------------------------- 00438 C CONTROL REACHES THIS SECTION IF 3 OR MORE FAILURES HAVE OCCURED. 00439 C IF 10 FAILURES HAVE OCCURRED, EXIT WITH KFLAG = -1. 00440 C IT IS ASSUMED THAT THE DERIVATIVES THAT HAVE ACCUMULATED IN THE 00441 C YH ARRAY HAVE ERRORS OF THE WRONG ORDER. HENCE THE FIRST 00442 C DERIVATIVE IS RECOMPUTED, AND THE ORDER IS SET TO 1. THEN 00443 C H IS REDUCED BY A FACTOR OF 10, AND THE STEP IS RETRIED, 00444 C UNTIL IT SUCCEEDS OR H REACHES HMIN. 00445 C----------------------------------------------------------------------- 00446 640 IF (KFLAG .EQ. -10) GO TO 660 00447 RH = 0.1D0 00448 RH = DMAX1(HMIN/DABS(H),RH) 00449 H = H*RH 00450 DO 645 I = 1,N 00451 645 Y(I) = YH(I,1) 00452 IERR = 0 00453 CALL F (NEQ, TN, Y, SAVF, IERR) 00454 IF (IERR .LT. 0) RETURN 00455 NFE = NFE + 1 00456 DO 650 I = 1,N 00457 650 YH(I,2) = H*SAVF(I) 00458 IPUP = MITER 00459 IALTH = 5 00460 IF (NQ .EQ. 1) GO TO 200 00461 NQ = 1 00462 L = 2 00463 IRET = 3 00464 GO TO 150 00465 C----------------------------------------------------------------------- 00466 C ALL RETURNS ARE MADE THROUGH THIS SECTION. H IS SAVED IN HOLD 00467 C TO ALLOW THE CALLER TO CHANGE H ON THE NEXT STEP. 00468 C----------------------------------------------------------------------- 00469 660 KFLAG = -1 00470 GO TO 720 00471 670 KFLAG = -2 00472 GO TO 720 00473 680 KFLAG = -3 00474 GO TO 720 00475 690 RMAX = 10.0D0 00476 700 R = 1.0D0/TESCO(2,NQU) 00477 DO 710 I = 1,N 00478 710 ACOR(I) = ACOR(I)*R 00479 720 HOLD = H 00480 JSTART = 1 00481 RETURN 00482 C----------------------- END OF SUBROUTINE STODE ----------------------- 00483 END
1.7.1
