GNU Octave: libcruft/daspk/ddstp.f Source File

Go to the documentation of this file.
00001 C Work performed under the auspices of the U.S. Department of Energy
00002 C by Lawrence Livermore National Laboratory under contract number 
00003 C W-7405-Eng-48.
00004 C
00005       SUBROUTINE DDSTP(X,Y,YPRIME,NEQ,RES,JAC,PSOL,H,WT,VT,
00006      *  JSTART,IDID,RPAR,IPAR,PHI,SAVR,DELTA,E,WM,IWM,
00007      *  ALPHA,BETA,GAMMA,PSI,SIGMA,CJ,CJOLD,HOLD,S,HMIN,UROUND,
00008      *  EPLI,SQRTN,RSQRTN,EPCON,IPHASE,JCALC,JFLG,K,KOLD,NS,NONNEG,
00009      *  NTYPE,NLS)
00010 C
00011 C***BEGIN PROLOGUE  DDSTP
00012 C***REFER TO  DDASPK
00013 C***DATE WRITTEN   890101   (YYMMDD)
00014 C***REVISION DATE  900926   (YYMMDD)
00015 C***REVISION DATE  940909   (YYMMDD) (Reset PSI(1), PHI(*,2) at 690)
00016 C
00017 C
00018 C-----------------------------------------------------------------------
00019 C***DESCRIPTION
00020 C
00021 C     DDSTP solves a system of differential/algebraic equations of 
00022 C     the form G(X,Y,YPRIME) = 0, for one step (normally from X to X+H).
00023 C
00024 C     The methods used are modified divided difference, fixed leading 
00025 C     coefficient forms of backward differentiation formulas.  
00026 C     The code adjusts the stepsize and order to control the local error
00027 C     per step.
00028 C
00029 C
00030 C     The parameters represent
00031 C     X  --        Independent variable.
00032 C     Y  --        Solution vector at X.
00033 C     YPRIME --    Derivative of solution vector
00034 C                  after successful step.
00035 C     NEQ --       Number of equations to be integrated.
00036 C     RES --       External user-supplied subroutine
00037 C                  to evaluate the residual.  See RES description
00038 C                  in DDASPK prologue.
00039 C     JAC --       External user-supplied routine to update
00040 C                  Jacobian or preconditioner information in the
00041 C                  nonlinear solver.  See JAC description in DDASPK
00042 C                  prologue.
00043 C     PSOL --      External user-supplied routine to solve
00044 C                  a linear system using preconditioning. 
00045 C                  (This is optional).  See PSOL in DDASPK prologue.
00046 C     H --         Appropriate step size for next step.
00047 C                  Normally determined by the code.
00048 C     WT --        Vector of weights for error criterion used in Newton test.
00049 C     VT --        Masked vector of weights used in error test.
00050 C     JSTART --    Integer variable set 0 for
00051 C                  first step, 1 otherwise.
00052 C     IDID --      Completion code returned from the nonlinear solver.
00053 C                  See IDID description in DDASPK prologue.
00054 C     RPAR,IPAR -- Real and integer parameter arrays that
00055 C                  are used for communication between the
00056 C                  calling program and external user routines.
00057 C                  They are not altered by DNSK
00058 C     PHI --       Array of divided differences used by
00059 C                  DDSTP. The length is NEQ*(K+1), where
00060 C                  K is the maximum order.
00061 C     SAVR --      Work vector for DDSTP of length NEQ.
00062 C     DELTA,E --   Work vectors for DDSTP of length NEQ.
00063 C     WM,IWM --    Real and integer arrays storing
00064 C                  information required by the linear solver.
00065 C
00066 C     The other parameters are information
00067 C     which is needed internally by DDSTP to
00068 C     continue from step to step.
00069 C
00070 C-----------------------------------------------------------------------
00071 C***ROUTINES CALLED
00072 C   NLS, DDWNRM, DDATRP
00073 C
00074 C***END PROLOGUE  DDSTP
00075 C
00076 C
00077       IMPLICIT DOUBLE PRECISION(A-H,O-Z)
00078       DIMENSION Y(*),YPRIME(*),WT(*),VT(*)
00079       DIMENSION PHI(NEQ,*),SAVR(*),DELTA(*),E(*)
00080       DIMENSION WM(*),IWM(*)
00081       DIMENSION PSI(*),ALPHA(*),BETA(*),GAMMA(*),SIGMA(*)
00082       DIMENSION RPAR(*),IPAR(*)
00083       EXTERNAL  RES, JAC, PSOL, NLS
00084 C
00085       PARAMETER (LMXORD=3)
00086       PARAMETER (LNST=11, LETF=14, LCFN=15)
00087 C
00088 C
00089 C-----------------------------------------------------------------------
00090 C     BLOCK 1.
00091 C     Initialize.  On the first call, set
00092 C     the order to 1 and initialize
00093 C     other variables.
00094 C-----------------------------------------------------------------------
00095 C
00096 C     Initializations for all calls
00097 C
00098       XOLD=X
00099       NCF=0
00100       NEF=0
00101       IF(JSTART .NE. 0) GO TO 120
00102 C
00103 C     If this is the first step, perform
00104 C     other initializations
00105 C
00106       K=1
00107       KOLD=0
00108       HOLD=0.0D0
00109       PSI(1)=H
00110       CJ = 1.D0/H
00111       IPHASE = 0
00112       NS=0
00113 120   CONTINUE
00114 C
00115 C
00116 C
00117 C
00118 C
00119 C-----------------------------------------------------------------------
00120 C     BLOCK 2
00121 C     Compute coefficients of formulas for
00122 C     this step.
00123 C-----------------------------------------------------------------------
00124 200   CONTINUE
00125       KP1=K+1
00126       KP2=K+2
00127       KM1=K-1
00128       IF(H.NE.HOLD.OR.K .NE. KOLD) NS = 0
00129       NS=MIN0(NS+1,KOLD+2)
00130       NSP1=NS+1
00131       IF(KP1 .LT. NS)GO TO 230
00132 C
00133       BETA(1)=1.0D0
00134       ALPHA(1)=1.0D0
00135       TEMP1=H
00136       GAMMA(1)=0.0D0
00137       SIGMA(1)=1.0D0
00138       DO 210 I=2,KP1
00139          TEMP2=PSI(I-1)
00140          PSI(I-1)=TEMP1
00141          BETA(I)=BETA(I-1)*PSI(I-1)/TEMP2
00142          TEMP1=TEMP2+H
00143          ALPHA(I)=H/TEMP1
00144          SIGMA(I)=(I-1)*SIGMA(I-1)*ALPHA(I)
00145          GAMMA(I)=GAMMA(I-1)+ALPHA(I-1)/H
00146 210      CONTINUE
00147       PSI(KP1)=TEMP1
00148 230   CONTINUE
00149 C
00150 C     Compute ALPHAS, ALPHA0
00151 C
00152       ALPHAS = 0.0D0
00153       ALPHA0 = 0.0D0
00154       DO 240 I = 1,K
00155         ALPHAS = ALPHAS - 1.0D0/I
00156         ALPHA0 = ALPHA0 - ALPHA(I)
00157 240     CONTINUE
00158 C
00159 C     Compute leading coefficient CJ
00160 C
00161       CJLAST = CJ
00162       CJ = -ALPHAS/H
00163 C
00164 C     Compute variable stepsize error coefficient CK
00165 C
00166       CK = ABS(ALPHA(KP1) + ALPHAS - ALPHA0)
00167       CK = MAX(CK,ALPHA(KP1))
00168 C
00169 C     Change PHI to PHI STAR
00170 C
00171       IF(KP1 .LT. NSP1) GO TO 280
00172       DO 270 J=NSP1,KP1
00173          DO 260 I=1,NEQ
00174 260         PHI(I,J)=BETA(J)*PHI(I,J)
00175 270      CONTINUE
00176 280   CONTINUE
00177 C
00178 C     Update time
00179 C
00180       X=X+H
00181 C
00182 C     Initialize IDID to 1
00183 C
00184       IDID = 1
00185 C
00186 C
00187 C
00188 C
00189 C
00190 C-----------------------------------------------------------------------
00191 C     BLOCK 3
00192 C     Call the nonlinear system solver to obtain the solution and
00193 C     derivative.
00194 C-----------------------------------------------------------------------
00195 C
00196       CALL NLS(X,Y,YPRIME,NEQ,
00197      *   RES,JAC,PSOL,H,WT,JSTART,IDID,RPAR,IPAR,PHI,GAMMA,
00198      *   SAVR,DELTA,E,WM,IWM,CJ,CJOLD,CJLAST,S,
00199      *   UROUND,EPLI,SQRTN,RSQRTN,EPCON,JCALC,JFLG,KP1,
00200      *   NONNEG,NTYPE,IERNLS)
00201 C
00202       IF(IERNLS .NE. 0)GO TO 600
00203 C
00204 C
00205 C
00206 C
00207 C
00208 C-----------------------------------------------------------------------
00209 C     BLOCK 4
00210 C     Estimate the errors at orders K,K-1,K-2
00211 C     as if constant stepsize was used. Estimate
00212 C     the local error at order K and test
00213 C     whether the current step is successful.
00214 C-----------------------------------------------------------------------
00215 C
00216 C     Estimate errors at orders K,K-1,K-2
00217 C
00218       ENORM = DDWNRM(NEQ,E,VT,RPAR,IPAR)
00219       ERK = SIGMA(K+1)*ENORM
00220       TERK = (K+1)*ERK
00221       EST = ERK
00222       KNEW=K
00223       IF(K .EQ. 1)GO TO 430
00224       DO 405 I = 1,NEQ
00225 405     DELTA(I) = PHI(I,KP1) + E(I)
00226       ERKM1=SIGMA(K)*DDWNRM(NEQ,DELTA,VT,RPAR,IPAR)
00227       TERKM1 = K*ERKM1
00228       IF(K .GT. 2)GO TO 410
00229       IF(TERKM1 .LE. 0.5*TERK)GO TO 420
00230       GO TO 430
00231 410   CONTINUE
00232       DO 415 I = 1,NEQ
00233 415     DELTA(I) = PHI(I,K) + DELTA(I)
00234       ERKM2=SIGMA(K-1)*DDWNRM(NEQ,DELTA,VT,RPAR,IPAR)
00235       TERKM2 = (K-1)*ERKM2
00236       IF(MAX(TERKM1,TERKM2).GT.TERK)GO TO 430
00237 C
00238 C     Lower the order
00239 C
00240 420   CONTINUE
00241       KNEW=K-1
00242       EST = ERKM1
00243 C
00244 C
00245 C     Calculate the local error for the current step
00246 C     to see if the step was successful
00247 C
00248 430   CONTINUE
00249       ERR = CK * ENORM
00250       IF(ERR .GT. 1.0D0)GO TO 600
00251 C
00252 C
00253 C
00254 C
00255 C
00256 C-----------------------------------------------------------------------
00257 C     BLOCK 5
00258 C     The step is successful. Determine
00259 C     the best order and stepsize for
00260 C     the next step. Update the differences
00261 C     for the next step.
00262 C-----------------------------------------------------------------------
00263       IDID=1
00264       IWM(LNST)=IWM(LNST)+1
00265       KDIFF=K-KOLD
00266       KOLD=K
00267       HOLD=H
00268 C
00269 C
00270 C     Estimate the error at order K+1 unless
00271 C        already decided to lower order, or
00272 C        already using maximum order, or
00273 C        stepsize not constant, or
00274 C        order raised in previous step
00275 C
00276       IF(KNEW.EQ.KM1.OR.K.EQ.IWM(LMXORD))IPHASE=1
00277       IF(IPHASE .EQ. 0)GO TO 545
00278       IF(KNEW.EQ.KM1)GO TO 540
00279       IF(K.EQ.IWM(LMXORD)) GO TO 550
00280       IF(KP1.GE.NS.OR.KDIFF.EQ.1)GO TO 550
00281       DO 510 I=1,NEQ
00282 510      DELTA(I)=E(I)-PHI(I,KP2)
00283       ERKP1 = (1.0D0/(K+2))*DDWNRM(NEQ,DELTA,VT,RPAR,IPAR)
00284       TERKP1 = (K+2)*ERKP1
00285       IF(K.GT.1)GO TO 520
00286       IF(TERKP1.GE.0.5D0*TERK)GO TO 550
00287       GO TO 530
00288 520   IF(TERKM1.LE.MIN(TERK,TERKP1))GO TO 540
00289       IF(TERKP1.GE.TERK.OR.K.EQ.IWM(LMXORD))GO TO 550
00290 C
00291 C     Raise order
00292 C
00293 530   K=KP1
00294       EST = ERKP1
00295       GO TO 550
00296 C
00297 C     Lower order
00298 C
00299 540   K=KM1
00300       EST = ERKM1
00301       GO TO 550
00302 C
00303 C     If IPHASE = 0, increase order by one and multiply stepsize by
00304 C     factor two
00305 C
00306 545   K = KP1
00307       HNEW = H*2.0D0
00308       H = HNEW
00309       GO TO 575
00310 C
00311 C
00312 C     Determine the appropriate stepsize for
00313 C     the next step.
00314 C
00315 550   HNEW=H
00316       TEMP2=K+1
00317       R=(2.0D0*EST+0.0001D0)**(-1.0D0/TEMP2)
00318       IF(R .LT. 2.0D0) GO TO 555
00319       HNEW = 2.0D0*H
00320       GO TO 560
00321 555   IF(R .GT. 1.0D0) GO TO 560
00322       R = MAX(0.5D0,MIN(0.9D0,R))
00323       HNEW = H*R
00324 560   H=HNEW
00325 C
00326 C
00327 C     Update differences for next step
00328 C
00329 575   CONTINUE
00330       IF(KOLD.EQ.IWM(LMXORD))GO TO 585
00331       DO 580 I=1,NEQ
00332 580      PHI(I,KP2)=E(I)
00333 585   CONTINUE
00334       DO 590 I=1,NEQ
00335 590      PHI(I,KP1)=PHI(I,KP1)+E(I)
00336       DO 595 J1=2,KP1
00337          J=KP1-J1+1
00338          DO 595 I=1,NEQ
00339 595      PHI(I,J)=PHI(I,J)+PHI(I,J+1)
00340       JSTART = 1
00341       RETURN
00342 C
00343 C
00344 C
00345 C
00346 C
00347 C-----------------------------------------------------------------------
00348 C     BLOCK 6
00349 C     The step is unsuccessful. Restore X,PSI,PHI
00350 C     Determine appropriate stepsize for
00351 C     continuing the integration, or exit with
00352 C     an error flag if there have been many
00353 C     failures.
00354 C-----------------------------------------------------------------------
00355 600   IPHASE = 1
00356 C
00357 C     Restore X,PHI,PSI
00358 C
00359       X=XOLD
00360       IF(KP1.LT.NSP1)GO TO 630
00361       DO 620 J=NSP1,KP1
00362          TEMP1=1.0D0/BETA(J)
00363          DO 610 I=1,NEQ
00364 610         PHI(I,J)=TEMP1*PHI(I,J)
00365 620      CONTINUE
00366 630   CONTINUE
00367       DO 640 I=2,KP1
00368 640      PSI(I-1)=PSI(I)-H
00369 C
00370 C
00371 C     Test whether failure is due to nonlinear solver
00372 C     or error test
00373 C
00374       IF(IERNLS .EQ. 0)GO TO 660
00375       IWM(LCFN)=IWM(LCFN)+1
00376 C
00377 C
00378 C     The nonlinear solver failed to converge.
00379 C     Determine the cause of the failure and take appropriate action.
00380 C     If IERNLS .LT. 0, then return.  Otherwise, reduce the stepsize
00381 C     and try again, unless too many failures have occurred.
00382 C
00383       IF (IERNLS .LT. 0) GO TO 675
00384       NCF = NCF + 1
00385       R = 0.25D0
00386       H = H*R
00387       IF (NCF .LT. 10 .AND. ABS(H) .GE. HMIN) GO TO 690
00388       IF (IDID .EQ. 1) IDID = -7
00389       IF (NEF .GE. 3) IDID = -9
00390       GO TO 675
00391 C
00392 C
00393 C     The nonlinear solver converged, and the cause
00394 C     of the failure was the error estimate
00395 C     exceeding the tolerance.
00396 C
00397 660   NEF=NEF+1
00398       IWM(LETF)=IWM(LETF)+1
00399       IF (NEF .GT. 1) GO TO 665
00400 C
00401 C     On first error test failure, keep current order or lower
00402 C     order by one.  Compute new stepsize based on differences
00403 C     of the solution.
00404 C
00405       K = KNEW
00406       TEMP2 = K + 1
00407       R = 0.90D0*(2.0D0*EST+0.0001D0)**(-1.0D0/TEMP2)
00408       R = MAX(0.25D0,MIN(0.9D0,R))
00409       H = H*R
00410       IF (ABS(H) .GE. HMIN) GO TO 690
00411       IDID = -6
00412       GO TO 675
00413 C
00414 C     On second error test failure, use the current order or
00415 C     decrease order by one.  Reduce the stepsize by a factor of
00416 C     one quarter.
00417 C
00418 665   IF (NEF .GT. 2) GO TO 670
00419       K = KNEW
00420       R = 0.25D0
00421       H = R*H
00422       IF (ABS(H) .GE. HMIN) GO TO 690
00423       IDID = -6
00424       GO TO 675
00425 C
00426 C     On third and subsequent error test failures, set the order to
00427 C     one, and reduce the stepsize by a factor of one quarter.
00428 C
00429 670   K = 1
00430       R = 0.25D0
00431       H = R*H
00432       IF (ABS(H) .GE. HMIN) GO TO 690
00433       IDID = -6
00434       GO TO 675
00435 C
00436 C
00437 C
00438 C
00439 C     For all crashes, restore Y to its last value,
00440 C     interpolate to find YPRIME at last X, and return.
00441 C
00442 C     Before returning, verify that the user has not set
00443 C     IDID to a nonnegative value.  If the user has set IDID
00444 C     to a nonnegative value, then reset IDID to be -7, indicating
00445 C     a failure in the nonlinear system solver.
00446 C
00447 675   CONTINUE
00448       CALL DDATRP(X,X,Y,YPRIME,NEQ,K,PHI,PSI)
00449       JSTART = 1
00450       IF (IDID .GE. 0) IDID = -7
00451       RETURN
00452 C
00453 C
00454 C     Go back and try this step again.  
00455 C     If this is the first step, reset PSI(1) and rescale PHI(*,2).
00456 C
00457 690   IF (KOLD .EQ. 0) THEN
00458         PSI(1) = H
00459         DO 695 I = 1,NEQ
00460 695       PHI(I,2) = R*PHI(I,2)
00461         ENDIF
00462       GO TO 200
00463 C
00464 C------END OF SUBROUTINE DDSTP------------------------------------------
00465       END