ddaini.f

Go to the documentation of this file.
00001       SUBROUTINE DDAINI (X, Y, YPRIME, NEQ, RES, JAC, H, WT, IDID, RPAR,
00002      +   IPAR, PHI, DELTA, E, WM, IWM, HMIN, UROUND, NONNEG, NTEMP)
00003 C***BEGIN PROLOGUE  DDAINI
00004 C***SUBSIDIARY
00005 C***PURPOSE  Initialization routine for DDASSL.
00006 C***LIBRARY   SLATEC (DASSL)
00007 C***TYPE      DOUBLE PRECISION (SDAINI-S, DDAINI-D)
00008 C***AUTHOR  PETZOLD, LINDA R., (LLNL)
00009 C***DESCRIPTION
00010 C-----------------------------------------------------------------
00011 C     DDAINI TAKES ONE STEP OF SIZE H OR SMALLER
00012 C     WITH THE BACKWARD EULER METHOD, TO
00013 C     FIND YPRIME.  X AND Y ARE UPDATED TO BE CONSISTENT WITH THE
00014 C     NEW STEP.  A MODIFIED DAMPED NEWTON ITERATION IS USED TO
00015 C     SOLVE THE CORRECTOR ITERATION.
00016 C
00017 C     THE INITIAL GUESS FOR YPRIME IS USED IN THE
00018 C     PREDICTION, AND IN FORMING THE ITERATION
00019 C     MATRIX, BUT IS NOT INVOLVED IN THE
00020 C     ERROR TEST. THIS MAY HAVE TROUBLE
00021 C     CONVERGING IF THE INITIAL GUESS IS NO
00022 C     GOOD, OR IF G(X,Y,YPRIME) DEPENDS
00023 C     NONLINEARLY ON YPRIME.
00024 C
00025 C     THE PARAMETERS REPRESENT:
00026 C     X --         INDEPENDENT VARIABLE
00027 C     Y --         SOLUTION VECTOR AT X
00028 C     YPRIME --    DERIVATIVE OF SOLUTION VECTOR
00029 C     NEQ --       NUMBER OF EQUATIONS
00030 C     H --         STEPSIZE. IMDER MAY USE A STEPSIZE
00031 C                  SMALLER THAN H.
00032 C     WT --        VECTOR OF WEIGHTS FOR ERROR
00033 C                  CRITERION
00034 C     IDID --      COMPLETION CODE WITH THE FOLLOWING MEANINGS
00035 C                  IDID= 1 -- YPRIME WAS FOUND SUCCESSFULLY
00036 C                  IDID=-12 -- DDAINI FAILED TO FIND YPRIME
00037 C     RPAR,IPAR -- REAL AND INTEGER PARAMETER ARRAYS
00038 C                  THAT ARE NOT ALTERED BY DDAINI
00039 C     PHI --       WORK SPACE FOR DDAINI
00040 C     DELTA,E --   WORK SPACE FOR DDAINI
00041 C     WM,IWM --    REAL AND INTEGER ARRAYS STORING
00042 C                  MATRIX INFORMATION
00043 C
00044 C-----------------------------------------------------------------
00045 C***ROUTINES CALLED  DDAJAC, DDANRM, DDASLV
00046 C***REVISION HISTORY  (YYMMDD)
00047 C   830315  DATE WRITTEN
00048 C   901009  Finished conversion to SLATEC 4.0 format (F.N.Fritsch)
00049 C   901019  Merged changes made by C. Ulrich with SLATEC 4.0 format.
00050 C   901026  Added explicit declarations for all variables and minor
00051 C           cosmetic changes to prologue.  (FNF)
00052 C   901030  Minor corrections to declarations.  (FNF)
00053 C***END PROLOGUE  DDAINI
00054 C
00055       INTEGER  NEQ, IDID, IPAR(*), IWM(*), NONNEG, NTEMP
00056       DOUBLE PRECISION
00057      *   X, Y(*), YPRIME(*), H, WT(*), RPAR(*), PHI(NEQ,*), DELTA(*),
00058      *   E(*), WM(*), HMIN, UROUND
00059       EXTERNAL  RES, JAC
00060 C
00061       EXTERNAL  DDAJAC, DDANRM, DDASLV
00062       DOUBLE PRECISION  DDANRM
00063 C
00064       INTEGER  I, IER, IRES, JCALC, LNJE, LNRE, M, MAXIT, MJAC, NCF,
00065      *   NEF, NSF
00066       DOUBLE PRECISION
00067      *   CJ, DAMP, DELNRM, ERR, OLDNRM, R, RATE, S, XOLD, YNORM
00068       LOGICAL  CONVGD
00069 C
00070       PARAMETER (LNRE=12)
00071       PARAMETER (LNJE=13)
00072 C
00073       DATA MAXIT/10/,MJAC/5/
00074       DATA DAMP/0.75D0/
00075 C
00076 C
00077 C---------------------------------------------------
00078 C     BLOCK 1.
00079 C     INITIALIZATIONS.
00080 C---------------------------------------------------
00081 C
00082 C***FIRST EXECUTABLE STATEMENT  DDAINI
00083       IDID=1
00084       NEF=0
00085       NCF=0
00086       NSF=0
00087       XOLD=X
00088       YNORM=DDANRM(NEQ,Y,WT,RPAR,IPAR)
00089 C
00090 C     SAVE Y AND YPRIME IN PHI
00091       DO 100 I=1,NEQ
00092          PHI(I,1)=Y(I)
00093 100      PHI(I,2)=YPRIME(I)
00094 C
00095 C
00096 C----------------------------------------------------
00097 C     BLOCK 2.
00098 C     DO ONE BACKWARD EULER STEP.
00099 C----------------------------------------------------
00100 C
00101 C     SET UP FOR START OF CORRECTOR ITERATION
00102 200   CJ=1.0D0/H
00103       X=X+H
00104 C
00105 C     PREDICT SOLUTION AND DERIVATIVE
00106       DO 250 I=1,NEQ
00107 250     Y(I)=Y(I)+H*YPRIME(I)
00108 C
00109       JCALC=-1
00110       M=0
00111       CONVGD=.TRUE.
00112 C
00113 C
00114 C     CORRECTOR LOOP.
00115 300   IWM(LNRE)=IWM(LNRE)+1
00116       IRES=0
00117 C
00118       CALL RES(X,Y,YPRIME,DELTA,IRES,RPAR,IPAR)
00119       IF (IRES.LT.0) GO TO 430
00120 C
00121 C
00122 C     EVALUATE THE ITERATION MATRIX
00123       IF (JCALC.NE.-1) GO TO 310
00124       IWM(LNJE)=IWM(LNJE)+1
00125       JCALC=0
00126       CALL DDAJAC(NEQ,X,Y,YPRIME,DELTA,CJ,H,
00127      *   IER,WT,E,WM,IWM,RES,IRES,
00128      *   UROUND,JAC,RPAR,IPAR,NTEMP)
00129 C
00130       S=1000000.D0
00131       IF (IRES.LT.0) GO TO 430
00132       IF (IER.NE.0) GO TO 430
00133       NSF=0
00134 C
00135 C
00136 C
00137 C     MULTIPLY RESIDUAL BY DAMPING FACTOR
00138 310   CONTINUE
00139       DO 320 I=1,NEQ
00140 320      DELTA(I)=DELTA(I)*DAMP
00141 C
00142 C     COMPUTE A NEW ITERATE (BACK SUBSTITUTION)
00143 C     STORE THE CORRECTION IN DELTA
00144 C
00145       CALL DDASLV(NEQ,DELTA,WM,IWM)
00146 C
00147 C     UPDATE Y AND YPRIME
00148       DO 330 I=1,NEQ
00149          Y(I)=Y(I)-DELTA(I)
00150 330      YPRIME(I)=YPRIME(I)-CJ*DELTA(I)
00151 C
00152 C     TEST FOR CONVERGENCE OF THE ITERATION.
00153 C
00154       DELNRM=DDANRM(NEQ,DELTA,WT,RPAR,IPAR)
00155       IF (DELNRM.LE.100.D0*UROUND*YNORM)
00156      *   GO TO 400
00157 C
00158       IF (M.GT.0) GO TO 340
00159          OLDNRM=DELNRM
00160          GO TO 350
00161 C
00162 340   RATE=(DELNRM/OLDNRM)**(1.0D0/M)
00163       IF (RATE.GT.0.90D0) GO TO 430
00164       S=RATE/(1.0D0-RATE)
00165 C
00166 350   IF (S*DELNRM .LE. 0.33D0) GO TO 400
00167 C
00168 C
00169 C     THE CORRECTOR HAS NOT YET CONVERGED. UPDATE
00170 C     M AND AND TEST WHETHER THE MAXIMUM
00171 C     NUMBER OF ITERATIONS HAVE BEEN TRIED.
00172 C     EVERY MJAC ITERATIONS, GET A NEW
00173 C     ITERATION MATRIX.
00174 C
00175       M=M+1
00176       IF (M.GE.MAXIT) GO TO 430
00177 C
00178       IF ((M/MJAC)*MJAC.EQ.M) JCALC=-1
00179       GO TO 300
00180 C
00181 C
00182 C     THE ITERATION HAS CONVERGED.
00183 C     CHECK NONNEGATIVITY CONSTRAINTS
00184 400   IF (NONNEG.EQ.0) GO TO 450
00185       DO 410 I=1,NEQ
00186 410      DELTA(I)=MIN(Y(I),0.0D0)
00187 C
00188       DELNRM=DDANRM(NEQ,DELTA,WT,RPAR,IPAR)
00189       IF (DELNRM.GT.0.33D0) GO TO 430
00190 C
00191       DO 420 I=1,NEQ
00192          Y(I)=Y(I)-DELTA(I)
00193 420      YPRIME(I)=YPRIME(I)-CJ*DELTA(I)
00194       GO TO 450
00195 C
00196 C
00197 C     EXITS FROM CORRECTOR LOOP.
00198 430   CONVGD=.FALSE.
00199 450   IF (.NOT.CONVGD) GO TO 600
00200 C
00201 C
00202 C
00203 C-----------------------------------------------------
00204 C     BLOCK 3.
00205 C     THE CORRECTOR ITERATION CONVERGED.
00206 C     DO ERROR TEST.
00207 C-----------------------------------------------------
00208 C
00209       DO 510 I=1,NEQ
00210 510      E(I)=Y(I)-PHI(I,1)
00211       ERR=DDANRM(NEQ,E,WT,RPAR,IPAR)
00212 C
00213       IF (ERR.LE.1.0D0) RETURN
00214 C
00215 C
00216 C
00217 C--------------------------------------------------------
00218 C     BLOCK 4.
00219 C     THE BACKWARD EULER STEP FAILED. RESTORE X, Y
00220 C     AND YPRIME TO THEIR ORIGINAL VALUES.
00221 C     REDUCE STEPSIZE AND TRY AGAIN, IF
00222 C     POSSIBLE.
00223 C---------------------------------------------------------
00224 C
00225 600   CONTINUE
00226       X = XOLD
00227       DO 610 I=1,NEQ
00228          Y(I)=PHI(I,1)
00229 610      YPRIME(I)=PHI(I,2)
00230 C
00231       IF (CONVGD) GO TO 640
00232       IF (IER.EQ.0) GO TO 620
00233          NSF=NSF+1
00234          H=H*0.25D0
00235          IF (NSF.LT.3.AND.ABS(H).GE.HMIN) GO TO 690
00236          IDID=-12
00237          RETURN
00238 620   IF (IRES.GT.-2) GO TO 630
00239          IDID=-12
00240          RETURN
00241 630   NCF=NCF+1
00242       H=H*0.25D0
00243       IF (NCF.LT.10.AND.ABS(H).GE.HMIN) GO TO 690
00244          IDID=-12
00245          RETURN
00246 C
00247 640   NEF=NEF+1
00248       R=0.90D0/(2.0D0*ERR+0.0001D0)
00249       R=MAX(0.1D0,MIN(0.5D0,R))
00250       H=H*R
00251       IF (ABS(H).GE.HMIN.AND.NEF.LT.10) GO TO 690
00252          IDID=-12
00253          RETURN
00254 690      GO TO 200
00255 C
00256 C-------------END OF SUBROUTINE DDAINI----------------------
00257       END
 All Classes Files Functions Variables Typedefs Enumerations Enumerator Friends Defines