00001 SUBROUTINE SSTODE (NEQ, Y, YH, NYH, YH1, EWT, SAVF, ACOR, 00002 1 WM, IWM, F, JAC, PJAC, SLVS) 00003 C***BEGIN PROLOGUE SSTODE 00004 C***SUBSIDIARY 00005 C***PURPOSE Performs one step of an ODEPACK integration. 00006 C***TYPE SINGLE PRECISION (SSTODE-S, DSTODE-D) 00007 C***AUTHOR Hindmarsh, Alan C., (LLNL) 00008 C***DESCRIPTION 00009 C 00010 C SSTODE performs one step of the integration of an initial value 00011 C problem for a system of ordinary differential equations. 00012 C Note: SSTODE is independent of the value of the iteration method 00013 C indicator MITER, when this is .ne. 0, and hence is independent 00014 C of the type of chord method used, or the Jacobian structure. 00015 C Communication with SSTODE is done with the following variables: 00016 C 00017 C NEQ = integer array containing problem size in NEQ(1), and 00018 C passed as the NEQ argument in all calls to F and JAC. 00019 C Y = an array of length .ge. N used as the Y argument in 00020 C all calls to F and JAC. 00021 C YH = an NYH by LMAX array containing the dependent variables 00022 C and their approximate scaled derivatives, where 00023 C LMAX = MAXORD + 1. YH(i,j+1) contains the approximate 00024 C j-th derivative of y(i), scaled by h**j/factorial(j) 00025 C (j = 0,1,...,NQ). on entry for the first step, the first 00026 C two columns of YH must be set from the initial values. 00027 C NYH = a constant integer .ge. N, the first dimension of YH. 00028 C YH1 = a one-dimensional array occupying the same space as YH. 00029 C EWT = an array of length N containing multiplicative weights 00030 C for local error measurements. Local errors in Y(i) are 00031 C compared to 1.0/EWT(i) in various error tests. 00032 C SAVF = an array of working storage, of length N. 00033 C Also used for input of YH(*,MAXORD+2) when JSTART = -1 00034 C and MAXORD .lt. the current order NQ. 00035 C ACOR = a work array of length N, used for the accumulated 00036 C corrections. On a successful return, ACOR(i) contains 00037 C the estimated one-step local error in Y(i). 00038 C WM,IWM = real and integer work arrays associated with matrix 00039 C operations in chord iteration (MITER .ne. 0). 00040 C PJAC = name of routine to evaluate and preprocess Jacobian matrix 00041 C and P = I - h*el0*JAC, if a chord method is being used. 00042 C SLVS = name of routine to solve linear system in chord iteration. 00043 C CCMAX = maximum relative change in h*el0 before PJAC is called. 00044 C H = the step size to be attempted on the next step. 00045 C H is altered by the error control algorithm during the 00046 C problem. H can be either positive or negative, but its 00047 C sign must remain constant throughout the problem. 00048 C HMIN = the minimum absolute value of the step size h to be used. 00049 C HMXI = inverse of the maximum absolute value of h to be used. 00050 C HMXI = 0.0 is allowed and corresponds to an infinite hmax. 00051 C HMIN and HMXI may be changed at any time, but will not 00052 C take effect until the next change of h is considered. 00053 C TN = the independent variable. TN is updated on each step taken. 00054 C JSTART = an integer used for input only, with the following 00055 C values and meanings: 00056 C 0 perform the first step. 00057 C .gt.0 take a new step continuing from the last. 00058 C -1 take the next step with a new value of H, MAXORD, 00059 C N, METH, MITER, and/or matrix parameters. 00060 C -2 take the next step with a new value of H, 00061 C but with other inputs unchanged. 00062 C On return, JSTART is set to 1 to facilitate continuation. 00063 C KFLAG = a completion code with the following meanings: 00064 C 0 the step was succesful. 00065 C -1 the requested error could not be achieved. 00066 C -2 corrector convergence could not be achieved. 00067 C -3 fatal error in PJAC or SLVS. 00068 C A return with KFLAG = -1 or -2 means either 00069 C abs(H) = HMIN or 10 consecutive failures occurred. 00070 C On a return with KFLAG negative, the values of TN and 00071 C the YH array are as of the beginning of the last 00072 C step, and H is the last step size attempted. 00073 C MAXORD = the maximum order of integration method to be allowed. 00074 C MAXCOR = the maximum number of corrector iterations allowed. 00075 C MSBP = maximum number of steps between PJAC calls (MITER .gt. 0). 00076 C MXNCF = maximum number of convergence failures allowed. 00077 C METH/MITER = the method flags. See description in driver. 00078 C N = the number of first-order differential equations. 00079 C The values of CCMAX, H, HMIN, HMXI, TN, JSTART, KFLAG, MAXORD, 00080 C MAXCOR, MSBP, MXNCF, METH, MITER, and N are communicated via COMMON. 00081 C 00082 C***SEE ALSO SLSODE 00083 C***ROUTINES CALLED SCFODE, SVNORM 00084 C***COMMON BLOCKS SLS001 00085 C***REVISION HISTORY (YYMMDD) 00086 C 791129 DATE WRITTEN 00087 C 890501 Modified prologue to SLATEC/LDOC format. (FNF) 00088 C 890503 Minor cosmetic changes. (FNF) 00089 C 930809 Renamed to allow single/double precision versions. (ACH) 00090 C 010413 Reduced size of Common block /SLS001/. (ACH) 00091 C 031105 Restored 'own' variables to Common block /SLS001/, to 00092 C enable interrupt/restart feature. (ACH) 00093 C***END PROLOGUE SSTODE 00094 C**End 00095 EXTERNAL F, JAC, PJAC, SLVS 00096 INTEGER NEQ, NYH, IWM 00097 REAL Y, YH, YH1, EWT, SAVF, ACOR, WM 00098 DIMENSION NEQ(*), Y(*), YH(NYH,*), YH1(*), EWT(*), SAVF(*), 00099 1 ACOR(*), WM(*), IWM(*) 00100 INTEGER IOWND, IALTH, IPUP, LMAX, MEO, NQNYH, NSLP, 00101 1 ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L, 00102 2 LYH, LEWT, LACOR, LSAVF, LWM, LIWM, METH, MITER, 00103 3 MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU 00104 INTEGER I, I1, IREDO, IRET, J, JB, M, NCF, NEWQ 00105 REAL CONIT, CRATE, EL, ELCO, HOLD, RMAX, TESCO, 00106 2 CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND 00107 REAL DCON, DDN, DEL, DELP, DSM, DUP, EXDN, EXSM, EXUP, 00108 1 R, RH, RHDN, RHSM, RHUP, TOLD, SVNORM 00109 COMMON /SLS001/ CONIT, CRATE, EL(13), ELCO(13,12), 00110 1 HOLD, RMAX, TESCO(3,12), 00111 2 CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND, 00112 3 IOWND(6), IALTH, IPUP, LMAX, MEO, NQNYH, NSLP, 00113 3 ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L, 00114 4 LYH, LEWT, LACOR, LSAVF, LWM, LIWM, METH, MITER, 00115 5 MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU 00116 C 00117 C***FIRST EXECUTABLE STATEMENT SSTODE 00118 KFLAG = 0 00119 TOLD = TN 00120 NCF = 0 00121 IERPJ = 0 00122 IERSL = 0 00123 JCUR = 0 00124 ICF = 0 00125 DELP = 0.0E0 00126 IF (JSTART .GT. 0) GO TO 200 00127 IF (JSTART .EQ. -1) GO TO 100 00128 IF (JSTART .EQ. -2) GO TO 160 00129 C----------------------------------------------------------------------- 00130 C On the first call, the order is set to 1, and other variables are 00131 C initialized. RMAX is the maximum ratio by which H can be increased 00132 C in a single step. It is initially 1.E4 to compensate for the small 00133 C initial H, but then is normally equal to 10. If a failure 00134 C occurs (in corrector convergence or error test), RMAX is set to 2 00135 C for the next increase. 00136 C----------------------------------------------------------------------- 00137 LMAX = MAXORD + 1 00138 NQ = 1 00139 L = 2 00140 IALTH = 2 00141 RMAX = 10000.0E0 00142 RC = 0.0E0 00143 EL0 = 1.0E0 00144 CRATE = 0.7E0 00145 HOLD = H 00146 MEO = METH 00147 NSLP = 0 00148 IPUP = MITER 00149 IRET = 3 00150 GO TO 140 00151 C----------------------------------------------------------------------- 00152 C The following block handles preliminaries needed when JSTART = -1. 00153 C IPUP is set to MITER to force a matrix update. 00154 C If an order increase is about to be considered (IALTH = 1), 00155 C IALTH is reset to 2 to postpone consideration one more step. 00156 C If the caller has changed METH, SCFODE is called to reset 00157 C the coefficients of the method. 00158 C If the caller has changed MAXORD to a value less than the current 00159 C order NQ, NQ is reduced to MAXORD, and a new H chosen accordingly. 00160 C If H is to be changed, YH must be rescaled. 00161 C If H or METH is being changed, IALTH is reset to L = NQ + 1 00162 C to prevent further changes in H for that many steps. 00163 C----------------------------------------------------------------------- 00164 100 IPUP = MITER 00165 LMAX = MAXORD + 1 00166 IF (IALTH .EQ. 1) IALTH = 2 00167 IF (METH .EQ. MEO) GO TO 110 00168 CALL SCFODE (METH, ELCO, TESCO) 00169 MEO = METH 00170 IF (NQ .GT. MAXORD) GO TO 120 00171 IALTH = L 00172 IRET = 1 00173 GO TO 150 00174 110 IF (NQ .LE. MAXORD) GO TO 160 00175 120 NQ = MAXORD 00176 L = LMAX 00177 DO 125 I = 1,L 00178 125 EL(I) = ELCO(I,NQ) 00179 NQNYH = NQ*NYH 00180 RC = RC*EL(1)/EL0 00181 EL0 = EL(1) 00182 CONIT = 0.5E0/(NQ+2) 00183 DDN = SVNORM (N, SAVF, EWT)/TESCO(1,L) 00184 EXDN = 1.0E0/L 00185 RHDN = 1.0E0/(1.3E0*DDN**EXDN + 0.0000013E0) 00186 RH = MIN(RHDN,1.0E0) 00187 IREDO = 3 00188 IF (H .EQ. HOLD) GO TO 170 00189 RH = MIN(RH,ABS(H/HOLD)) 00190 H = HOLD 00191 GO TO 175 00192 C----------------------------------------------------------------------- 00193 C SCFODE is called to get all the integration coefficients for the 00194 C current METH. Then the EL vector and related constants are reset 00195 C whenever the order NQ is changed, or at the start of the problem. 00196 C----------------------------------------------------------------------- 00197 140 CALL SCFODE (METH, ELCO, TESCO) 00198 150 DO 155 I = 1,L 00199 155 EL(I) = ELCO(I,NQ) 00200 NQNYH = NQ*NYH 00201 RC = RC*EL(1)/EL0 00202 EL0 = EL(1) 00203 CONIT = 0.5E0/(NQ+2) 00204 GO TO (160, 170, 200), IRET 00205 C----------------------------------------------------------------------- 00206 C If H is being changed, the H ratio RH is checked against 00207 C RMAX, HMIN, and HMXI, and the YH array rescaled. IALTH is set to 00208 C L = NQ + 1 to prevent a change of H for that many steps, unless 00209 C forced by a convergence or error test failure. 00210 C----------------------------------------------------------------------- 00211 160 IF (H .EQ. HOLD) GO TO 200 00212 RH = H/HOLD 00213 H = HOLD 00214 IREDO = 3 00215 GO TO 175 00216 170 RH = MAX(RH,HMIN/ABS(H)) 00217 175 RH = MIN(RH,RMAX) 00218 RH = RH/MAX(1.0E0,ABS(H)*HMXI*RH) 00219 R = 1.0E0 00220 DO 180 J = 2,L 00221 R = R*RH 00222 DO 180 I = 1,N 00223 180 YH(I,J) = YH(I,J)*R 00224 H = H*RH 00225 RC = RC*RH 00226 IALTH = L 00227 IF (IREDO .EQ. 0) GO TO 690 00228 C----------------------------------------------------------------------- 00229 C This section computes the predicted values by effectively 00230 C multiplying the YH array by the Pascal Triangle matrix. 00231 C RC is the ratio of new to old values of the coefficient H*EL(1). 00232 C When RC differs from 1 by more than CCMAX, IPUP is set to MITER 00233 C to force PJAC to be called, if a Jacobian is involved. 00234 C In any case, PJAC is called at least every MSBP steps. 00235 C----------------------------------------------------------------------- 00236 200 IF (ABS(RC-1.0E0) .GT. CCMAX) IPUP = MITER 00237 IF (NST .GE. NSLP+MSBP) IPUP = MITER 00238 TN = TN + H 00239 I1 = NQNYH + 1 00240 DO 215 JB = 1,NQ 00241 I1 = I1 - NYH 00242 Cdir$ ivdep 00243 DO 210 I = I1,NQNYH 00244 210 YH1(I) = YH1(I) + YH1(I+NYH) 00245 215 CONTINUE 00246 C----------------------------------------------------------------------- 00247 C Up to MAXCOR corrector iterations are taken. A convergence test is 00248 C made on the R.M.S. norm of each correction, weighted by the error 00249 C weight vector EWT. The sum of the corrections is accumulated in the 00250 C vector ACOR(i). The YH array is not altered in the corrector loop. 00251 C----------------------------------------------------------------------- 00252 220 M = 0 00253 DO 230 I = 1,N 00254 230 Y(I) = YH(I,1) 00255 CALL F (NEQ, TN, Y, SAVF) 00256 NFE = NFE + 1 00257 IF (IPUP .LE. 0) GO TO 250 00258 C----------------------------------------------------------------------- 00259 C If indicated, the matrix P = I - h*el(1)*J is reevaluated and 00260 C preprocessed before starting the corrector iteration. IPUP is set 00261 C to 0 as an indicator that this has been done. 00262 C----------------------------------------------------------------------- 00263 CALL PJAC (NEQ, Y, YH, NYH, EWT, ACOR, SAVF, WM, IWM, F, JAC) 00264 IPUP = 0 00265 RC = 1.0E0 00266 NSLP = NST 00267 CRATE = 0.7E0 00268 IF (IERPJ .NE. 0) GO TO 430 00269 250 DO 260 I = 1,N 00270 260 ACOR(I) = 0.0E0 00271 270 IF (MITER .NE. 0) GO TO 350 00272 C----------------------------------------------------------------------- 00273 C In the case of functional iteration, update Y directly from 00274 C the result of the last function evaluation. 00275 C----------------------------------------------------------------------- 00276 DO 290 I = 1,N 00277 SAVF(I) = H*SAVF(I) - YH(I,2) 00278 290 Y(I) = SAVF(I) - ACOR(I) 00279 DEL = SVNORM (N, Y, EWT) 00280 DO 300 I = 1,N 00281 Y(I) = YH(I,1) + EL(1)*SAVF(I) 00282 300 ACOR(I) = SAVF(I) 00283 GO TO 400 00284 C----------------------------------------------------------------------- 00285 C In the case of the chord method, compute the corrector error, 00286 C and solve the linear system with that as right-hand side and 00287 C P as coefficient matrix. 00288 C----------------------------------------------------------------------- 00289 350 DO 360 I = 1,N 00290 360 Y(I) = H*SAVF(I) - (YH(I,2) + ACOR(I)) 00291 CALL SLVS (WM, IWM, Y, SAVF) 00292 IF (IERSL .LT. 0) GO TO 430 00293 IF (IERSL .GT. 0) GO TO 410 00294 DEL = SVNORM (N, Y, EWT) 00295 DO 380 I = 1,N 00296 ACOR(I) = ACOR(I) + Y(I) 00297 380 Y(I) = YH(I,1) + EL(1)*ACOR(I) 00298 C----------------------------------------------------------------------- 00299 C Test for convergence. If M.gt.0, an estimate of the convergence 00300 C rate constant is stored in CRATE, and this is used in the test. 00301 C----------------------------------------------------------------------- 00302 400 IF (M .NE. 0) CRATE = MAX(0.2E0*CRATE,DEL/DELP) 00303 DCON = DEL*MIN(1.0E0,1.5E0*CRATE)/(TESCO(2,NQ)*CONIT) 00304 IF (DCON .LE. 1.0E0) GO TO 450 00305 M = M + 1 00306 IF (M .EQ. MAXCOR) GO TO 410 00307 IF (M .GE. 2 .AND. DEL .GT. 2.0E0*DELP) GO TO 410 00308 DELP = DEL 00309 CALL F (NEQ, TN, Y, SAVF) 00310 NFE = NFE + 1 00311 GO TO 270 00312 C----------------------------------------------------------------------- 00313 C The corrector iteration failed to converge. 00314 C If MITER .ne. 0 and the Jacobian is out of date, PJAC is called for 00315 C the next try. Otherwise the YH array is retracted to its values 00316 C before prediction, and H is reduced, if possible. If H cannot be 00317 C reduced or MXNCF failures have occurred, exit with KFLAG = -2. 00318 C----------------------------------------------------------------------- 00319 410 IF (MITER .EQ. 0 .OR. JCUR .EQ. 1) GO TO 430 00320 ICF = 1 00321 IPUP = MITER 00322 GO TO 220 00323 430 ICF = 2 00324 NCF = NCF + 1 00325 RMAX = 2.0E0 00326 TN = TOLD 00327 I1 = NQNYH + 1 00328 DO 445 JB = 1,NQ 00329 I1 = I1 - NYH 00330 Cdir$ ivdep 00331 DO 440 I = I1,NQNYH 00332 440 YH1(I) = YH1(I) - YH1(I+NYH) 00333 445 CONTINUE 00334 IF (IERPJ .LT. 0 .OR. IERSL .LT. 0) GO TO 680 00335 IF (ABS(H) .LE. HMIN*1.00001E0) GO TO 670 00336 IF (NCF .EQ. MXNCF) GO TO 670 00337 RH = 0.25E0 00338 IPUP = MITER 00339 IREDO = 1 00340 GO TO 170 00341 C----------------------------------------------------------------------- 00342 C The corrector has converged. JCUR is set to 0 00343 C to signal that the Jacobian involved may need updating later. 00344 C The local error test is made and control passes to statement 500 00345 C if it fails. 00346 C----------------------------------------------------------------------- 00347 450 JCUR = 0 00348 IF (M .EQ. 0) DSM = DEL/TESCO(2,NQ) 00349 IF (M .GT. 0) DSM = SVNORM (N, ACOR, EWT)/TESCO(2,NQ) 00350 IF (DSM .GT. 1.0E0) GO TO 500 00351 C----------------------------------------------------------------------- 00352 C After a successful step, update the YH array. 00353 C Consider changing H if IALTH = 1. Otherwise decrease IALTH by 1. 00354 C If IALTH is then 1 and NQ .lt. MAXORD, then ACOR is saved for 00355 C use in a possible order increase on the next step. 00356 C If a change in H is considered, an increase or decrease in order 00357 C by one is considered also. A change in H is made only if it is by a 00358 C factor of at least 1.1. If not, IALTH is set to 3 to prevent 00359 C testing for that many steps. 00360 C----------------------------------------------------------------------- 00361 KFLAG = 0 00362 IREDO = 0 00363 NST = NST + 1 00364 HU = H 00365 NQU = NQ 00366 DO 470 J = 1,L 00367 DO 470 I = 1,N 00368 470 YH(I,J) = YH(I,J) + EL(J)*ACOR(I) 00369 IALTH = IALTH - 1 00370 IF (IALTH .EQ. 0) GO TO 520 00371 IF (IALTH .GT. 1) GO TO 700 00372 IF (L .EQ. LMAX) GO TO 700 00373 DO 490 I = 1,N 00374 490 YH(I,LMAX) = ACOR(I) 00375 GO TO 700 00376 C----------------------------------------------------------------------- 00377 C The error test failed. KFLAG keeps track of multiple failures. 00378 C Restore TN and the YH array to their previous values, and prepare 00379 C to try the step again. Compute the optimum step size for this or 00380 C one lower order. After 2 or more failures, H is forced to decrease 00381 C by a factor of 0.2 or less. 00382 C----------------------------------------------------------------------- 00383 500 KFLAG = KFLAG - 1 00384 TN = TOLD 00385 I1 = NQNYH + 1 00386 DO 515 JB = 1,NQ 00387 I1 = I1 - NYH 00388 Cdir$ ivdep 00389 DO 510 I = I1,NQNYH 00390 510 YH1(I) = YH1(I) - YH1(I+NYH) 00391 515 CONTINUE 00392 RMAX = 2.0E0 00393 IF (ABS(H) .LE. HMIN*1.00001E0) GO TO 660 00394 IF (KFLAG .LE. -3) GO TO 640 00395 IREDO = 2 00396 RHUP = 0.0E0 00397 GO TO 540 00398 C----------------------------------------------------------------------- 00399 C Regardless of the success or failure of the step, factors 00400 C RHDN, RHSM, and RHUP are computed, by which H could be multiplied 00401 C at order NQ - 1, order NQ, or order NQ + 1, respectively. 00402 C In the case of failure, RHUP = 0.0 to avoid an order increase. 00403 C The largest of these is determined and the new order chosen 00404 C accordingly. If the order is to be increased, we compute one 00405 C additional scaled derivative. 00406 C----------------------------------------------------------------------- 00407 520 RHUP = 0.0E0 00408 IF (L .EQ. LMAX) GO TO 540 00409 DO 530 I = 1,N 00410 530 SAVF(I) = ACOR(I) - YH(I,LMAX) 00411 DUP = SVNORM (N, SAVF, EWT)/TESCO(3,NQ) 00412 EXUP = 1.0E0/(L+1) 00413 RHUP = 1.0E0/(1.4E0*DUP**EXUP + 0.0000014E0) 00414 540 EXSM = 1.0E0/L 00415 RHSM = 1.0E0/(1.2E0*DSM**EXSM + 0.0000012E0) 00416 RHDN = 0.0E0 00417 IF (NQ .EQ. 1) GO TO 560 00418 DDN = SVNORM (N, YH(1,L), EWT)/TESCO(1,NQ) 00419 EXDN = 1.0E0/NQ 00420 RHDN = 1.0E0/(1.3E0*DDN**EXDN + 0.0000013E0) 00421 560 IF (RHSM .GE. RHUP) GO TO 570 00422 IF (RHUP .GT. RHDN) GO TO 590 00423 GO TO 580 00424 570 IF (RHSM .LT. RHDN) GO TO 580 00425 NEWQ = NQ 00426 RH = RHSM 00427 GO TO 620 00428 580 NEWQ = NQ - 1 00429 RH = RHDN 00430 IF (KFLAG .LT. 0 .AND. RH .GT. 1.0E0) RH = 1.0E0 00431 GO TO 620 00432 590 NEWQ = L 00433 RH = RHUP 00434 IF (RH .LT. 1.1E0) GO TO 610 00435 R = EL(L)/L 00436 DO 600 I = 1,N 00437 600 YH(I,NEWQ+1) = ACOR(I)*R 00438 GO TO 630 00439 610 IALTH = 3 00440 GO TO 700 00441 620 IF ((KFLAG .EQ. 0) .AND. (RH .LT. 1.1E0)) GO TO 610 00442 IF (KFLAG .LE. -2) RH = MIN(RH,0.2E0) 00443 C----------------------------------------------------------------------- 00444 C If there is a change of order, reset NQ, l, and the coefficients. 00445 C In any case H is reset according to RH and the YH array is rescaled. 00446 C Then exit from 690 if the step was OK, or redo the step otherwise. 00447 C----------------------------------------------------------------------- 00448 IF (NEWQ .EQ. NQ) GO TO 170 00449 630 NQ = NEWQ 00450 L = NQ + 1 00451 IRET = 2 00452 GO TO 150 00453 C----------------------------------------------------------------------- 00454 C Control reaches this section if 3 or more failures have occured. 00455 C If 10 failures have occurred, exit with KFLAG = -1. 00456 C It is assumed that the derivatives that have accumulated in the 00457 C YH array have errors of the wrong order. Hence the first 00458 C derivative is recomputed, and the order is set to 1. Then 00459 C H is reduced by a factor of 10, and the step is retried, 00460 C until it succeeds or H reaches HMIN. 00461 C----------------------------------------------------------------------- 00462 640 IF (KFLAG .EQ. -10) GO TO 660 00463 RH = 0.1E0 00464 RH = MAX(HMIN/ABS(H),RH) 00465 H = H*RH 00466 DO 645 I = 1,N 00467 645 Y(I) = YH(I,1) 00468 CALL F (NEQ, TN, Y, SAVF) 00469 NFE = NFE + 1 00470 DO 650 I = 1,N 00471 650 YH(I,2) = H*SAVF(I) 00472 IPUP = MITER 00473 IALTH = 5 00474 IF (NQ .EQ. 1) GO TO 200 00475 NQ = 1 00476 L = 2 00477 IRET = 3 00478 GO TO 150 00479 C----------------------------------------------------------------------- 00480 C All returns are made through this section. H is saved in HOLD 00481 C to allow the caller to change H on the next step. 00482 C----------------------------------------------------------------------- 00483 660 KFLAG = -1 00484 GO TO 720 00485 670 KFLAG = -2 00486 GO TO 720 00487 680 KFLAG = -3 00488 GO TO 720 00489 690 RMAX = 10.0E0 00490 700 R = 1.0E0/TESCO(2,NQU) 00491 DO 710 I = 1,N 00492 710 ACOR(I) = ACOR(I)*R 00493 720 HOLD = H 00494 JSTART = 1 00495 RETURN 00496 C----------------------- END OF SUBROUTINE SSTODE ---------------------- 00497 END
1.7.1
