ddastp.f

Go to the documentation of this file.
00001       SUBROUTINE DDASTP (X, Y, YPRIME, NEQ, RES, JAC, H, WT, JSTART,
00002      +   IDID, RPAR, IPAR, PHI, DELTA, E, WM, IWM, ALPHA, BETA, GAMMA,
00003      +   PSI, SIGMA, CJ, CJOLD, HOLD, S, HMIN, UROUND, IPHASE, JCALC,
00004      +   K, KOLD, NS, NONNEG, NTEMP)
00005 C***BEGIN PROLOGUE  DDASTP
00006 C***SUBSIDIARY
00007 C***PURPOSE  Perform one step of the DDASSL integration.
00008 C***LIBRARY   SLATEC (DASSL)
00009 C***TYPE      DOUBLE PRECISION (SDASTP-S, DDASTP-D)
00010 C***AUTHOR  PETZOLD, LINDA R., (LLNL)
00011 C***DESCRIPTION
00012 C-----------------------------------------------------------------------
00013 C     DDASTP SOLVES A SYSTEM OF DIFFERENTIAL/
00014 C     ALGEBRAIC EQUATIONS OF THE FORM
00015 C     G(X,Y,YPRIME) = 0,  FOR ONE STEP (NORMALLY
00016 C     FROM X TO X+H).
00017 C
00018 C     THE METHODS USED ARE MODIFIED DIVIDED
00019 C     DIFFERENCE,FIXED LEADING COEFFICIENT
00020 C     FORMS OF BACKWARD DIFFERENTIATION
00021 C     FORMULAS. THE CODE ADJUSTS THE STEPSIZE
00022 C     AND ORDER TO CONTROL THE LOCAL ERROR PER
00023 C     STEP.
00024 C
00025 C
00026 C     THE PARAMETERS REPRESENT
00027 C     X  --        INDEPENDENT VARIABLE
00028 C     Y  --        SOLUTION VECTOR AT X
00029 C     YPRIME --    DERIVATIVE OF SOLUTION VECTOR
00030 C                  AFTER SUCCESSFUL STEP
00031 C     NEQ --       NUMBER OF EQUATIONS TO BE INTEGRATED
00032 C     RES --       EXTERNAL USER-SUPPLIED SUBROUTINE
00033 C                  TO EVALUATE THE RESIDUAL.  THE CALL IS
00034 C                  CALL RES(X,Y,YPRIME,DELTA,IRES,RPAR,IPAR)
00035 C                  X,Y,YPRIME ARE INPUT.  DELTA IS OUTPUT.
00036 C                  ON INPUT, IRES=0.  RES SHOULD ALTER IRES ONLY
00037 C                  IF IT ENCOUNTERS AN ILLEGAL VALUE OF Y OR A
00038 C                  STOP CONDITION.  SET IRES=-1 IF AN INPUT VALUE
00039 C                  OF Y IS ILLEGAL, AND DDASTP WILL TRY TO SOLVE
00040 C                  THE PROBLEM WITHOUT GETTING IRES = -1.  IF
00041 C                  IRES=-2, DDASTP RETURNS CONTROL TO THE CALLING
00042 C                  PROGRAM WITH IDID = -11.
00043 C     JAC --       EXTERNAL USER-SUPPLIED ROUTINE TO EVALUATE
00044 C                  THE ITERATION MATRIX (THIS IS OPTIONAL)
00045 C                  THE CALL IS OF THE FORM
00046 C                  CALL JAC(X,Y,YPRIME,PD,CJ,RPAR,IPAR)
00047 C                  PD IS THE MATRIX OF PARTIAL DERIVATIVES,
00048 C                  PD=DG/DY+CJ*DG/DYPRIME
00049 C     H --         APPROPRIATE STEP SIZE FOR NEXT STEP.
00050 C                  NORMALLY DETERMINED BY THE CODE
00051 C     WT --        VECTOR OF WEIGHTS FOR ERROR CRITERION.
00052 C     JSTART --    INTEGER VARIABLE SET 0 FOR
00053 C                  FIRST STEP, 1 OTHERWISE.
00054 C     IDID --      COMPLETION CODE WITH THE FOLLOWING MEANINGS:
00055 C                  IDID= 1 -- THE STEP WAS COMPLETED SUCCESSFULLY
00056 C                  IDID=-6 -- THE ERROR TEST FAILED REPEATEDLY
00057 C                  IDID=-7 -- THE CORRECTOR COULD NOT CONVERGE
00058 C                  IDID=-8 -- THE ITERATION MATRIX IS SINGULAR
00059 C                  IDID=-9 -- THE CORRECTOR COULD NOT CONVERGE.
00060 C                             THERE WERE REPEATED ERROR TEST
00061 C                             FAILURES ON THIS STEP.
00062 C                  IDID=-10-- THE CORRECTOR COULD NOT CONVERGE
00063 C                             BECAUSE IRES WAS EQUAL TO MINUS ONE
00064 C                  IDID=-11-- IRES EQUAL TO -2 WAS ENCOUNTERED,
00065 C                             AND CONTROL IS BEING RETURNED TO
00066 C                             THE CALLING PROGRAM
00067 C     RPAR,IPAR -- REAL AND INTEGER PARAMETER ARRAYS THAT
00068 C                  ARE USED FOR COMMUNICATION BETWEEN THE
00069 C                  CALLING PROGRAM AND EXTERNAL USER ROUTINES
00070 C                  THEY ARE NOT ALTERED BY DDASTP
00071 C     PHI --       ARRAY OF DIVIDED DIFFERENCES USED BY
00072 C                  DDASTP. THE LENGTH IS NEQ*(K+1),WHERE
00073 C                  K IS THE MAXIMUM ORDER
00074 C     DELTA,E --   WORK VECTORS FOR DDASTP OF LENGTH NEQ
00075 C     WM,IWM --    REAL AND INTEGER ARRAYS STORING
00076 C                  MATRIX INFORMATION SUCH AS THE MATRIX
00077 C                  OF PARTIAL DERIVATIVES,PERMUTATION
00078 C                  VECTOR,AND VARIOUS OTHER INFORMATION.
00079 C
00080 C     THE OTHER PARAMETERS ARE INFORMATION
00081 C     WHICH IS NEEDED INTERNALLY BY DDASTP TO
00082 C     CONTINUE FROM STEP TO STEP.
00083 C
00084 C-----------------------------------------------------------------------
00085 C***ROUTINES CALLED  DDAJAC, DDANRM, DDASLV, DDATRP
00086 C***REVISION HISTORY  (YYMMDD)
00087 C   830315  DATE WRITTEN
00088 C   901009  Finished conversion to SLATEC 4.0 format (F.N.Fritsch)
00089 C   901019  Merged changes made by C. Ulrich with SLATEC 4.0 format.
00090 C   901026  Added explicit declarations for all variables and minor
00091 C           cosmetic changes to prologue.  (FNF)
00092 C***END PROLOGUE  DDASTP
00093 C
00094       INTEGER  NEQ, JSTART, IDID, IPAR(*), IWM(*), IPHASE, JCALC, K,
00095      *   KOLD, NS, NONNEG, NTEMP
00096       DOUBLE PRECISION
00097      *   X, Y(*), YPRIME(*), H, WT(*), RPAR(*), PHI(NEQ,*), DELTA(*),
00098      *   E(*), WM(*), ALPHA(*), BETA(*), GAMMA(*), PSI(*), SIGMA(*), CJ,
00099      *   CJOLD, HOLD, S, HMIN, UROUND
00100       EXTERNAL  RES, JAC
00101 C
00102       EXTERNAL  DDAJAC, DDANRM, DDASLV, DDATRP
00103       DOUBLE PRECISION  DDANRM
00104 C
00105       INTEGER  I, IER, IRES, J, J1, KDIFF, KM1, KNEW, KP1, KP2, LCTF,
00106      *   LETF, LMXORD, LNJE, LNRE, LNST, M, MAXIT, NCF, NEF, NSF, NSP1
00107       DOUBLE PRECISION
00108      *   ALPHA0, ALPHAS, CJLAST, CK, DELNRM, ENORM, ERK, ERKM1,
00109      *   ERKM2, ERKP1, ERR, EST, HNEW, OLDNRM, PNORM, R, RATE, TEMP1,
00110      *   TEMP2, TERK, TERKM1, TERKM2, TERKP1, XOLD, XRATE
00111       LOGICAL  CONVGD
00112 C
00113       PARAMETER (LMXORD=3)
00114       PARAMETER (LNST=11)
00115       PARAMETER (LNRE=12)
00116       PARAMETER (LNJE=13)
00117       PARAMETER (LETF=14)
00118       PARAMETER (LCTF=15)
00119 C
00120       DATA MAXIT/4/
00121       DATA XRATE/0.25D0/
00122 C
00123 C
00124 C
00125 C
00126 C
00127 C-----------------------------------------------------------------------
00128 C     BLOCK 1.
00129 C     INITIALIZE. ON THE FIRST CALL,SET
00130 C     THE ORDER TO 1 AND INITIALIZE
00131 C     OTHER VARIABLES.
00132 C-----------------------------------------------------------------------
00133 C
00134 C     INITIALIZATIONS FOR ALL CALLS
00135 C***FIRST EXECUTABLE STATEMENT  DDASTP
00136       IDID=1
00137       XOLD=X
00138       NCF=0
00139       NSF=0
00140       NEF=0
00141       IF(JSTART .NE. 0) GO TO 120
00142 C
00143 C     IF THIS IS THE FIRST STEP,PERFORM
00144 C     OTHER INITIALIZATIONS
00145       IWM(LETF) = 0
00146       IWM(LCTF) = 0
00147       K=1
00148       KOLD=0
00149       HOLD=0.0D0
00150       JSTART=1
00151       PSI(1)=H
00152       CJOLD = 1.0D0/H
00153       CJ = CJOLD
00154       S = 100.D0
00155       JCALC = -1
00156       DELNRM=1.0D0
00157       IPHASE = 0
00158       NS=0
00159 120   CONTINUE
00160 C
00161 C
00162 C
00163 C
00164 C
00165 C-----------------------------------------------------------------------
00166 C     BLOCK 2
00167 C     COMPUTE COEFFICIENTS OF FORMULAS FOR
00168 C     THIS STEP.
00169 C-----------------------------------------------------------------------
00170 200   CONTINUE
00171       KP1=K+1
00172       KP2=K+2
00173       KM1=K-1
00174       XOLD=X
00175       IF(H.NE.HOLD.OR.K .NE. KOLD) NS = 0
00176       NS=MIN(NS+1,KOLD+2)
00177       NSP1=NS+1
00178       IF(KP1 .LT. NS)GO TO 230
00179 C
00180       BETA(1)=1.0D0
00181       ALPHA(1)=1.0D0
00182       TEMP1=H
00183       GAMMA(1)=0.0D0
00184       SIGMA(1)=1.0D0
00185       DO 210 I=2,KP1
00186          TEMP2=PSI(I-1)
00187          PSI(I-1)=TEMP1
00188          BETA(I)=BETA(I-1)*PSI(I-1)/TEMP2
00189          TEMP1=TEMP2+H
00190          ALPHA(I)=H/TEMP1
00191          SIGMA(I)=(I-1)*SIGMA(I-1)*ALPHA(I)
00192          GAMMA(I)=GAMMA(I-1)+ALPHA(I-1)/H
00193 210      CONTINUE
00194       PSI(KP1)=TEMP1
00195 230   CONTINUE
00196 C
00197 C     COMPUTE ALPHAS, ALPHA0
00198       ALPHAS = 0.0D0
00199       ALPHA0 = 0.0D0
00200       DO 240 I = 1,K
00201         ALPHAS = ALPHAS - 1.0D0/I
00202         ALPHA0 = ALPHA0 - ALPHA(I)
00203 240     CONTINUE
00204 C
00205 C     COMPUTE LEADING COEFFICIENT CJ
00206       CJLAST = CJ
00207       CJ = -ALPHAS/H
00208 C
00209 C     COMPUTE VARIABLE STEPSIZE ERROR COEFFICIENT CK
00210       CK = ABS(ALPHA(KP1) + ALPHAS - ALPHA0)
00211       CK = MAX(CK,ALPHA(KP1))
00212 C
00213 C     DECIDE WHETHER NEW JACOBIAN IS NEEDED
00214       TEMP1 = (1.0D0 - XRATE)/(1.0D0 + XRATE)
00215       TEMP2 = 1.0D0/TEMP1
00216       IF (CJ/CJOLD .LT. TEMP1 .OR. CJ/CJOLD .GT. TEMP2) JCALC = -1
00217       IF (CJ .NE. CJLAST) S = 100.D0
00218 C
00219 C     CHANGE PHI TO PHI STAR
00220       IF(KP1 .LT. NSP1) GO TO 280
00221       DO 270 J=NSP1,KP1
00222          DO 260 I=1,NEQ
00223 260         PHI(I,J)=BETA(J)*PHI(I,J)
00224 270      CONTINUE
00225 280   CONTINUE
00226 C
00227 C     UPDATE TIME
00228       X=X+H
00229 C
00230 C
00231 C
00232 C
00233 C
00234 C-----------------------------------------------------------------------
00235 C     BLOCK 3
00236 C     PREDICT THE SOLUTION AND DERIVATIVE,
00237 C     AND SOLVE THE CORRECTOR EQUATION
00238 C-----------------------------------------------------------------------
00239 C
00240 C     FIRST,PREDICT THE SOLUTION AND DERIVATIVE
00241 300   CONTINUE
00242       DO 310 I=1,NEQ
00243          Y(I)=PHI(I,1)
00244 310      YPRIME(I)=0.0D0
00245       DO 330 J=2,KP1
00246          DO 320 I=1,NEQ
00247             Y(I)=Y(I)+PHI(I,J)
00248 320         YPRIME(I)=YPRIME(I)+GAMMA(J)*PHI(I,J)
00249 330   CONTINUE
00250       PNORM = DDANRM (NEQ,Y,WT,RPAR,IPAR)
00251 C
00252 C
00253 C
00254 C     SOLVE THE CORRECTOR EQUATION USING A
00255 C     MODIFIED NEWTON SCHEME.
00256       CONVGD= .TRUE.
00257       M=0
00258       IWM(LNRE)=IWM(LNRE)+1
00259       IRES = 0
00260       CALL RES(X,Y,YPRIME,DELTA,IRES,RPAR,IPAR)
00261       IF (IRES .LT. 0) GO TO 380
00262 C
00263 C
00264 C     IF INDICATED,REEVALUATE THE
00265 C     ITERATION MATRIX PD = DG/DY + CJ*DG/DYPRIME
00266 C     (WHERE G(X,Y,YPRIME)=0). SET
00267 C     JCALC TO 0 AS AN INDICATOR THAT
00268 C     THIS HAS BEEN DONE.
00269       IF(JCALC .NE. -1)GO TO 340
00270       IWM(LNJE)=IWM(LNJE)+1
00271       JCALC=0
00272       CALL DDAJAC(NEQ,X,Y,YPRIME,DELTA,CJ,H,
00273      * IER,WT,E,WM,IWM,RES,IRES,UROUND,JAC,RPAR,
00274      * IPAR,NTEMP)
00275       CJOLD=CJ
00276       S = 100.D0
00277       IF (IRES .LT. 0) GO TO 380
00278       IF(IER .NE. 0)GO TO 380
00279       NSF=0
00280 C
00281 C
00282 C     INITIALIZE THE ERROR ACCUMULATION VECTOR E.
00283 340   CONTINUE
00284       DO 345 I=1,NEQ
00285 345      E(I)=0.0D0
00286 C
00287 C
00288 C     CORRECTOR LOOP.
00289 350   CONTINUE
00290 C
00291 C     MULTIPLY RESIDUAL BY TEMP1 TO ACCELERATE CONVERGENCE
00292       TEMP1 = 2.0D0/(1.0D0 + CJ/CJOLD)
00293       DO 355 I = 1,NEQ
00294 355     DELTA(I) = DELTA(I) * TEMP1
00295 C
00296 C     COMPUTE A NEW ITERATE (BACK-SUBSTITUTION).
00297 C     STORE THE CORRECTION IN DELTA.
00298       CALL DDASLV(NEQ,DELTA,WM,IWM)
00299 C
00300 C     UPDATE Y,E,AND YPRIME
00301       DO 360 I=1,NEQ
00302          Y(I)=Y(I)-DELTA(I)
00303          E(I)=E(I)-DELTA(I)
00304 360      YPRIME(I)=YPRIME(I)-CJ*DELTA(I)
00305 C
00306 C     TEST FOR CONVERGENCE OF THE ITERATION
00307       DELNRM=DDANRM(NEQ,DELTA,WT,RPAR,IPAR)
00308       IF (DELNRM .LE. 100.D0*UROUND*PNORM) GO TO 375
00309       IF (M .GT. 0) GO TO 365
00310          OLDNRM = DELNRM
00311          GO TO 367
00312 365   RATE = (DELNRM/OLDNRM)**(1.0D0/M)
00313       IF (RATE .GT. 0.90D0) GO TO 370
00314       S = RATE/(1.0D0 - RATE)
00315 367   IF (S*DELNRM .LE. 0.33D0) GO TO 375
00316 C
00317 C     THE CORRECTOR HAS NOT YET CONVERGED.
00318 C     UPDATE M AND TEST WHETHER THE
00319 C     MAXIMUM NUMBER OF ITERATIONS HAVE
00320 C     BEEN TRIED.
00321       M=M+1
00322       IF(M.GE.MAXIT)GO TO 370
00323 C
00324 C     EVALUATE THE RESIDUAL
00325 C     AND GO BACK TO DO ANOTHER ITERATION
00326       IWM(LNRE)=IWM(LNRE)+1
00327       IRES = 0
00328       CALL RES(X,Y,YPRIME,DELTA,IRES,
00329      *  RPAR,IPAR)
00330       IF (IRES .LT. 0) GO TO 380
00331       GO TO 350
00332 C
00333 C
00334 C     THE CORRECTOR FAILED TO CONVERGE IN MAXIT
00335 C     ITERATIONS. IF THE ITERATION MATRIX
00336 C     IS NOT CURRENT,RE-DO THE STEP WITH
00337 C     A NEW ITERATION MATRIX.
00338 370   CONTINUE
00339       IF(JCALC.EQ.0)GO TO 380
00340       JCALC=-1
00341       GO TO 300
00342 C
00343 C
00344 C     THE ITERATION HAS CONVERGED.  IF NONNEGATIVITY OF SOLUTION IS
00345 C     REQUIRED, SET THE SOLUTION NONNEGATIVE, IF THE PERTURBATION
00346 C     TO DO IT IS SMALL ENOUGH.  IF THE CHANGE IS TOO LARGE, THEN
00347 C     CONSIDER THE CORRECTOR ITERATION TO HAVE FAILED.
00348 375   IF(NONNEG .EQ. 0) GO TO 390
00349       DO 377 I = 1,NEQ
00350 377      DELTA(I) = MIN(Y(I),0.0D0)
00351       DELNRM = DDANRM(NEQ,DELTA,WT,RPAR,IPAR)
00352       IF(DELNRM .GT. 0.33D0) GO TO 380
00353       DO 378 I = 1,NEQ
00354 378      E(I) = E(I) - DELTA(I)
00355       GO TO 390
00356 C
00357 C
00358 C     EXITS FROM BLOCK 3
00359 C     NO CONVERGENCE WITH CURRENT ITERATION
00360 C     MATRIX,OR SINGULAR ITERATION MATRIX
00361 380   CONVGD= .FALSE.
00362 390   JCALC = 1
00363       IF(.NOT.CONVGD)GO TO 600
00364 C
00365 C
00366 C
00367 C
00368 C
00369 C-----------------------------------------------------------------------
00370 C     BLOCK 4
00371 C     ESTIMATE THE ERRORS AT ORDERS K,K-1,K-2
00372 C     AS IF CONSTANT STEPSIZE WAS USED. ESTIMATE
00373 C     THE LOCAL ERROR AT ORDER K AND TEST
00374 C     WHETHER THE CURRENT STEP IS SUCCESSFUL.
00375 C-----------------------------------------------------------------------
00376 C
00377 C     ESTIMATE ERRORS AT ORDERS K,K-1,K-2
00378       ENORM = DDANRM(NEQ,E,WT,RPAR,IPAR)
00379       ERK = SIGMA(K+1)*ENORM
00380       TERK = (K+1)*ERK
00381       EST = ERK
00382       KNEW=K
00383       IF(K .EQ. 1)GO TO 430
00384       DO 405 I = 1,NEQ
00385 405     DELTA(I) = PHI(I,KP1) + E(I)
00386       ERKM1=SIGMA(K)*DDANRM(NEQ,DELTA,WT,RPAR,IPAR)
00387       TERKM1 = K*ERKM1
00388       IF(K .GT. 2)GO TO 410
00389       IF(TERKM1 .LE. 0.5D0*TERK)GO TO 420
00390       GO TO 430
00391 410   CONTINUE
00392       DO 415 I = 1,NEQ
00393 415     DELTA(I) = PHI(I,K) + DELTA(I)
00394       ERKM2=SIGMA(K-1)*DDANRM(NEQ,DELTA,WT,RPAR,IPAR)
00395       TERKM2 = (K-1)*ERKM2
00396       IF(MAX(TERKM1,TERKM2).GT.TERK)GO TO 430
00397 C     LOWER THE ORDER
00398 420   CONTINUE
00399       KNEW=K-1
00400       EST = ERKM1
00401 C
00402 C
00403 C     CALCULATE THE LOCAL ERROR FOR THE CURRENT STEP
00404 C     TO SEE IF THE STEP WAS SUCCESSFUL
00405 430   CONTINUE
00406       ERR = CK * ENORM
00407       IF(ERR .GT. 1.0D0)GO TO 600
00408 C
00409 C
00410 C
00411 C
00412 C
00413 C-----------------------------------------------------------------------
00414 C     BLOCK 5
00415 C     THE STEP IS SUCCESSFUL. DETERMINE
00416 C     THE BEST ORDER AND STEPSIZE FOR
00417 C     THE NEXT STEP. UPDATE THE DIFFERENCES
00418 C     FOR THE NEXT STEP.
00419 C-----------------------------------------------------------------------
00420       IDID=1
00421       IWM(LNST)=IWM(LNST)+1
00422       KDIFF=K-KOLD
00423       KOLD=K
00424       HOLD=H
00425 C
00426 C
00427 C     ESTIMATE THE ERROR AT ORDER K+1 UNLESS:
00428 C        ALREADY DECIDED TO LOWER ORDER, OR
00429 C        ALREADY USING MAXIMUM ORDER, OR
00430 C        STEPSIZE NOT CONSTANT, OR
00431 C        ORDER RAISED IN PREVIOUS STEP
00432       IF(KNEW.EQ.KM1.OR.K.EQ.IWM(LMXORD))IPHASE=1
00433       IF(IPHASE .EQ. 0)GO TO 545
00434       IF(KNEW.EQ.KM1)GO TO 540
00435       IF(K.EQ.IWM(LMXORD)) GO TO 550
00436       IF(KP1.GE.NS.OR.KDIFF.EQ.1)GO TO 550
00437       DO 510 I=1,NEQ
00438 510      DELTA(I)=E(I)-PHI(I,KP2)
00439       ERKP1 = (1.0D0/(K+2))*DDANRM(NEQ,DELTA,WT,RPAR,IPAR)
00440       TERKP1 = (K+2)*ERKP1
00441       IF(K.GT.1)GO TO 520
00442       IF(TERKP1.GE.0.5D0*TERK)GO TO 550
00443       GO TO 530
00444 520   IF(TERKM1.LE.MIN(TERK,TERKP1))GO TO 540
00445       IF(TERKP1.GE.TERK.OR.K.EQ.IWM(LMXORD))GO TO 550
00446 C
00447 C     RAISE ORDER
00448 530   K=KP1
00449       EST = ERKP1
00450       GO TO 550
00451 C
00452 C     LOWER ORDER
00453 540   K=KM1
00454       EST = ERKM1
00455       GO TO 550
00456 C
00457 C     IF IPHASE = 0, INCREASE ORDER BY ONE AND MULTIPLY STEPSIZE BY
00458 C     FACTOR TWO
00459 545   K = KP1
00460       HNEW = H*2.0D0
00461       H = HNEW
00462       GO TO 575
00463 C
00464 C
00465 C     DETERMINE THE APPROPRIATE STEPSIZE FOR
00466 C     THE NEXT STEP.
00467 550   HNEW=H
00468       TEMP2=K+1
00469       R=(2.0D0*EST+0.0001D0)**(-1.0D0/TEMP2)
00470       IF(R .LT. 2.0D0) GO TO 555
00471       HNEW = 2.0D0*H
00472       GO TO 560
00473 555   IF(R .GT. 1.0D0) GO TO 560
00474       R = MAX(0.5D0,MIN(0.9D0,R))
00475       HNEW = H*R
00476 560   H=HNEW
00477 C
00478 C
00479 C     UPDATE DIFFERENCES FOR NEXT STEP
00480 575   CONTINUE
00481       IF(KOLD.EQ.IWM(LMXORD))GO TO 585
00482       DO 580 I=1,NEQ
00483 580      PHI(I,KP2)=E(I)
00484 585   CONTINUE
00485       DO 590 I=1,NEQ
00486 590      PHI(I,KP1)=PHI(I,KP1)+E(I)
00487       DO 595 J1=2,KP1
00488          J=KP1-J1+1
00489          DO 595 I=1,NEQ
00490 595      PHI(I,J)=PHI(I,J)+PHI(I,J+1)
00491       RETURN
00492 C
00493 C
00494 C
00495 C
00496 C
00497 C-----------------------------------------------------------------------
00498 C     BLOCK 6
00499 C     THE STEP IS UNSUCCESSFUL. RESTORE X,PSI,PHI
00500 C     DETERMINE APPROPRIATE STEPSIZE FOR
00501 C     CONTINUING THE INTEGRATION, OR EXIT WITH
00502 C     AN ERROR FLAG IF THERE HAVE BEEN MANY
00503 C     FAILURES.
00504 C-----------------------------------------------------------------------
00505 600   IPHASE = 1
00506 C
00507 C     RESTORE X,PHI,PSI
00508       X=XOLD
00509       IF(KP1.LT.NSP1)GO TO 630
00510       DO 620 J=NSP1,KP1
00511          TEMP1=1.0D0/BETA(J)
00512          DO 610 I=1,NEQ
00513 610         PHI(I,J)=TEMP1*PHI(I,J)
00514 620      CONTINUE
00515 630   CONTINUE
00516       DO 640 I=2,KP1
00517 640      PSI(I-1)=PSI(I)-H
00518 C
00519 C
00520 C     TEST WHETHER FAILURE IS DUE TO CORRECTOR ITERATION
00521 C     OR ERROR TEST
00522       IF(CONVGD)GO TO 660
00523       IWM(LCTF)=IWM(LCTF)+1
00524 C
00525 C
00526 C     THE NEWTON ITERATION FAILED TO CONVERGE WITH
00527 C     A CURRENT ITERATION MATRIX.  DETERMINE THE CAUSE
00528 C     OF THE FAILURE AND TAKE APPROPRIATE ACTION.
00529       IF(IER.EQ.0)GO TO 650
00530 C
00531 C     THE ITERATION MATRIX IS SINGULAR. REDUCE
00532 C     THE STEPSIZE BY A FACTOR OF 4. IF
00533 C     THIS HAPPENS THREE TIMES IN A ROW ON
00534 C     THE SAME STEP, RETURN WITH AN ERROR FLAG
00535       NSF=NSF+1
00536       R = 0.25D0
00537       H=H*R
00538       IF (NSF .LT. 3 .AND. ABS(H) .GE. HMIN) GO TO 690
00539       IDID=-8
00540       GO TO 675
00541 C
00542 C
00543 C     THE NEWTON ITERATION FAILED TO CONVERGE FOR A REASON
00544 C     OTHER THAN A SINGULAR ITERATION MATRIX.  IF IRES = -2, THEN
00545 C     RETURN.  OTHERWISE, REDUCE THE STEPSIZE AND TRY AGAIN, UNLESS
00546 C     TOO MANY FAILURES HAVE OCCURED.
00547 650   CONTINUE
00548       IF (IRES .GT. -2) GO TO 655
00549       IDID = -11
00550       GO TO 675
00551 655   NCF = NCF + 1
00552       R = 0.25D0
00553       H = H*R
00554       IF (NCF .LT. 10 .AND. ABS(H) .GE. HMIN) GO TO 690
00555       IDID = -7
00556       IF (IRES .LT. 0) IDID = -10
00557       IF (NEF .GE. 3) IDID = -9
00558       GO TO 675
00559 C
00560 C
00561 C     THE NEWTON SCHEME CONVERGED,AND THE CAUSE
00562 C     OF THE FAILURE WAS THE ERROR ESTIMATE
00563 C     EXCEEDING THE TOLERANCE.
00564 660   NEF=NEF+1
00565       IWM(LETF)=IWM(LETF)+1
00566       IF (NEF .GT. 1) GO TO 665
00567 C
00568 C     ON FIRST ERROR TEST FAILURE, KEEP CURRENT ORDER OR LOWER
00569 C     ORDER BY ONE.  COMPUTE NEW STEPSIZE BASED ON DIFFERENCES
00570 C     OF THE SOLUTION.
00571       K = KNEW
00572       TEMP2 = K + 1
00573       R = 0.90D0*(2.0D0*EST+0.0001D0)**(-1.0D0/TEMP2)
00574       R = MAX(0.25D0,MIN(0.9D0,R))
00575       H = H*R
00576       IF (ABS(H) .GE. HMIN) GO TO 690
00577       IDID = -6
00578       GO TO 675
00579 C
00580 C     ON SECOND ERROR TEST FAILURE, USE THE CURRENT ORDER OR
00581 C     DECREASE ORDER BY ONE.  REDUCE THE STEPSIZE BY A FACTOR OF
00582 C     FOUR.
00583 665   IF (NEF .GT. 2) GO TO 670
00584       K = KNEW
00585       H = 0.25D0*H
00586       IF (ABS(H) .GE. HMIN) GO TO 690
00587       IDID = -6
00588       GO TO 675
00589 C
00590 C     ON THIRD AND SUBSEQUENT ERROR TEST FAILURES, SET THE ORDER TO
00591 C     ONE AND REDUCE THE STEPSIZE BY A FACTOR OF FOUR.
00592 670   K = 1
00593       H = 0.25D0*H
00594       IF (ABS(H) .GE. HMIN) GO TO 690
00595       IDID = -6
00596       GO TO 675
00597 C
00598 C
00599 C
00600 C
00601 C     FOR ALL CRASHES, RESTORE Y TO ITS LAST VALUE,
00602 C     INTERPOLATE TO FIND YPRIME AT LAST X, AND RETURN
00603 675   CONTINUE
00604       CALL DDATRP(X,X,Y,YPRIME,NEQ,K,PHI,PSI)
00605       RETURN
00606 C
00607 C
00608 C     GO BACK AND TRY THIS STEP AGAIN
00609 690   GO TO 200
00610 C
00611 C------END OF SUBROUTINE DDASTP------
00612       END
 All Classes Files Functions Variables Typedefs Enumerations Enumerator Friends Defines