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 DNEDD(X,Y,YPRIME,NEQ,RES,JACD,PDUM,H,WT, 00006 * JSTART,IDID,RPAR,IPAR,PHI,GAMMA,DUMSVR,DELTA,E, 00007 * WM,IWM,CJ,CJOLD,CJLAST,S,UROUND,DUME,DUMS,DUMR, 00008 * EPCON,JCALC,JFDUM,KP1,NONNEG,NTYPE,IERNLS) 00009 C 00010 C***BEGIN PROLOGUE DNEDD 00011 C***REFER TO DDASPK 00012 C***DATE WRITTEN 891219 (YYMMDD) 00013 C***REVISION DATE 900926 (YYMMDD) 00014 C 00015 C 00016 C----------------------------------------------------------------------- 00017 C***DESCRIPTION 00018 C 00019 C DNEDD solves a nonlinear system of 00020 C algebraic equations of the form 00021 C G(X,Y,YPRIME) = 0 for the unknown Y. 00022 C 00023 C The method used is a modified Newton scheme. 00024 C 00025 C The parameters represent 00026 C 00027 C X -- Independent variable. 00028 C Y -- Solution vector. 00029 C YPRIME -- Derivative of solution vector. 00030 C NEQ -- Number of unknowns. 00031 C RES -- External user-supplied subroutine 00032 C to evaluate the residual. See RES description 00033 C in DDASPK prologue. 00034 C JACD -- External user-supplied routine to evaluate the 00035 C Jacobian. See JAC description for the case 00036 C INFO(12) = 0 in the DDASPK prologue. 00037 C PDUM -- Dummy argument. 00038 C H -- Appropriate step size for next step. 00039 C WT -- Vector of weights for error criterion. 00040 C JSTART -- Indicates first call to this routine. 00041 C If JSTART = 0, then this is the first call, 00042 C otherwise it is not. 00043 C IDID -- Completion flag, output by DNEDD. 00044 C See IDID description in DDASPK prologue. 00045 C RPAR,IPAR -- Real and integer arrays used for communication 00046 C between the calling program and external user 00047 C routines. They are not altered within DASPK. 00048 C PHI -- Array of divided differences used by 00049 C DNEDD. The length is NEQ*(K+1),where 00050 C K is the maximum order. 00051 C GAMMA -- Array used to predict Y and YPRIME. The length 00052 C is MAXORD+1 where MAXORD is the maximum order. 00053 C DUMSVR -- Dummy argument. 00054 C DELTA -- Work vector for NLS of length NEQ. 00055 C E -- Error accumulation vector for NLS of length NEQ. 00056 C WM,IWM -- Real and integer arrays storing 00057 C matrix information such as the matrix 00058 C of partial derivatives, permutation 00059 C vector, and various other information. 00060 C CJ -- Parameter always proportional to 1/H. 00061 C CJOLD -- Saves the value of CJ as of the last call to DMATD. 00062 C Accounts for changes in CJ needed to 00063 C decide whether to call DMATD. 00064 C CJLAST -- Previous value of CJ. 00065 C S -- A scalar determined by the approximate rate 00066 C of convergence of the Newton iteration and used 00067 C in the convergence test for the Newton iteration. 00068 C 00069 C If RATE is defined to be an estimate of the 00070 C rate of convergence of the Newton iteration, 00071 C then S = RATE/(1.D0-RATE). 00072 C 00073 C The closer RATE is to 0., the faster the Newton 00074 C iteration is converging; the closer RATE is to 1., 00075 C the slower the Newton iteration is converging. 00076 C 00077 C On the first Newton iteration with an up-dated 00078 C preconditioner S = 100.D0, Thus the initial 00079 C RATE of convergence is approximately 1. 00080 C 00081 C S is preserved from call to call so that the rate 00082 C estimate from a previous step can be applied to 00083 C the current step. 00084 C UROUND -- Unit roundoff. 00085 C DUME -- Dummy argument. 00086 C DUMS -- Dummy argument. 00087 C DUMR -- Dummy argument. 00088 C EPCON -- Tolerance to test for convergence of the Newton 00089 C iteration. 00090 C JCALC -- Flag used to determine when to update 00091 C the Jacobian matrix. In general: 00092 C 00093 C JCALC = -1 ==> Call the DMATD routine to update 00094 C the Jacobian matrix. 00095 C JCALC = 0 ==> Jacobian matrix is up-to-date. 00096 C JCALC = 1 ==> Jacobian matrix is out-dated, 00097 C but DMATD will not be called unless 00098 C JCALC is set to -1. 00099 C JFDUM -- Dummy argument. 00100 C KP1 -- The current order(K) + 1; updated across calls. 00101 C NONNEG -- Flag to determine nonnegativity constraints. 00102 C NTYPE -- Identification code for the NLS routine. 00103 C 0 ==> modified Newton; direct solver. 00104 C IERNLS -- Error flag for nonlinear solver. 00105 C 0 ==> nonlinear solver converged. 00106 C 1 ==> recoverable error inside nonlinear solver. 00107 C -1 ==> unrecoverable error inside nonlinear solver. 00108 C 00109 C All variables with "DUM" in their names are dummy variables 00110 C which are not used in this routine. 00111 C 00112 C Following is a list and description of local variables which 00113 C may not have an obvious usage. They are listed in roughly the 00114 C order they occur in this subroutine. 00115 C 00116 C The following group of variables are passed as arguments to 00117 C the Newton iteration solver. They are explained in greater detail 00118 C in DNSD: 00119 C TOLNEW, MULDEL, MAXIT, IERNEW 00120 C 00121 C IERTYP -- Flag which tells whether this subroutine is correct. 00122 C 0 ==> correct subroutine. 00123 C 1 ==> incorrect subroutine. 00124 C 00125 C----------------------------------------------------------------------- 00126 C***ROUTINES CALLED 00127 C DDWNRM, RES, DMATD, DNSD 00128 C 00129 C***END PROLOGUE DNEDD 00130 C 00131 C 00132 IMPLICIT DOUBLE PRECISION(A-H,O-Z) 00133 DIMENSION Y(*),YPRIME(*),WT(*) 00134 DIMENSION DELTA(*),E(*) 00135 DIMENSION WM(*),IWM(*), RPAR(*),IPAR(*) 00136 DIMENSION PHI(NEQ,*),GAMMA(*) 00137 EXTERNAL RES, JACD 00138 C 00139 PARAMETER (LNRE=12, LNJE=13) 00140 C 00141 SAVE MULDEL, MAXIT, XRATE 00142 DATA MULDEL/1/, MAXIT/4/, XRATE/0.25D0/ 00143 C 00144 C Verify that this is the correct subroutine. 00145 C 00146 IERTYP = 0 00147 IF (NTYPE .NE. 0) THEN 00148 IERTYP = 1 00149 GO TO 380 00150 ENDIF 00151 C 00152 C If this is the first step, perform initializations. 00153 C 00154 IF (JSTART .EQ. 0) THEN 00155 CJOLD = CJ 00156 JCALC = -1 00157 ENDIF 00158 C 00159 C Perform all other initializations. 00160 C 00161 IERNLS = 0 00162 C 00163 C Decide whether new Jacobian is needed. 00164 C 00165 TEMP1 = (1.0D0 - XRATE)/(1.0D0 + XRATE) 00166 TEMP2 = 1.0D0/TEMP1 00167 IF (CJ/CJOLD .LT. TEMP1 .OR. CJ/CJOLD .GT. TEMP2) JCALC = -1 00168 IF (CJ .NE. CJLAST) S = 100.D0 00169 C 00170 C----------------------------------------------------------------------- 00171 C Entry point for updating the Jacobian with current 00172 C stepsize. 00173 C----------------------------------------------------------------------- 00174 300 CONTINUE 00175 C 00176 C Initialize all error flags to zero. 00177 C 00178 IERJ = 0 00179 IRES = 0 00180 IERNEW = 0 00181 C 00182 C Predict the solution and derivative and compute the tolerance 00183 C for the Newton iteration. 00184 C 00185 DO 310 I=1,NEQ 00186 Y(I)=PHI(I,1) 00187 310 YPRIME(I)=0.0D0 00188 DO 330 J=2,KP1 00189 DO 320 I=1,NEQ 00190 Y(I)=Y(I)+PHI(I,J) 00191 320 YPRIME(I)=YPRIME(I)+GAMMA(J)*PHI(I,J) 00192 330 CONTINUE 00193 PNORM = DDWNRM (NEQ,Y,WT,RPAR,IPAR) 00194 TOLNEW = 100.D0*UROUND*PNORM 00195 C 00196 C Call RES to initialize DELTA. 00197 C 00198 IWM(LNRE)=IWM(LNRE)+1 00199 CALL RES(X,Y,YPRIME,CJ,DELTA,IRES,RPAR,IPAR) 00200 IF (IRES .LT. 0) GO TO 380 00201 C 00202 C If indicated, reevaluate the iteration matrix 00203 C J = dG/dY + CJ*dG/dYPRIME (where G(X,Y,YPRIME)=0). 00204 C Set JCALC to 0 as an indicator that this has been done. 00205 C 00206 IF(JCALC .EQ. -1) THEN 00207 IWM(LNJE)=IWM(LNJE)+1 00208 JCALC=0 00209 CALL DMATD(NEQ,X,Y,YPRIME,DELTA,CJ,H,IERJ,WT,E,WM,IWM, 00210 * RES,IRES,UROUND,JACD,RPAR,IPAR) 00211 CJOLD=CJ 00212 S = 100.D0 00213 IF (IRES .LT. 0) GO TO 380 00214 IF(IERJ .NE. 0)GO TO 380 00215 ENDIF 00216 C 00217 C Call the nonlinear Newton solver. 00218 C 00219 TEMP1 = 2.0D0/(1.0D0 + CJ/CJOLD) 00220 CALL DNSD(X,Y,YPRIME,NEQ,RES,PDUM,WT,RPAR,IPAR,DUMSVR, 00221 * DELTA,E,WM,IWM,CJ,DUMS,DUMR,DUME,EPCON,S,TEMP1, 00222 * TOLNEW,MULDEL,MAXIT,IRES,IDUM,IERNEW) 00223 C 00224 IF (IERNEW .GT. 0 .AND. JCALC .NE. 0) THEN 00225 C 00226 C The Newton iteration had a recoverable failure with an old 00227 C iteration matrix. Retry the step with a new iteration matrix. 00228 C 00229 JCALC = -1 00230 GO TO 300 00231 ENDIF 00232 C 00233 IF (IERNEW .NE. 0) GO TO 380 00234 C 00235 C The Newton iteration has converged. If nonnegativity of 00236 C solution is required, set the solution nonnegative, if the 00237 C perturbation to do it is small enough. If the change is too 00238 C large, then consider the corrector iteration to have failed. 00239 C 00240 375 IF(NONNEG .EQ. 0) GO TO 390 00241 DO 377 I = 1,NEQ 00242 377 DELTA(I) = MIN(Y(I),0.0D0) 00243 DELNRM = DDWNRM(NEQ,DELTA,WT,RPAR,IPAR) 00244 IF(DELNRM .GT. EPCON) GO TO 380 00245 DO 378 I = 1,NEQ 00246 378 E(I) = E(I) - DELTA(I) 00247 GO TO 390 00248 C 00249 C 00250 C Exits from nonlinear solver. 00251 C No convergence with current iteration 00252 C matrix, or singular iteration matrix. 00253 C Compute IERNLS and IDID accordingly. 00254 C 00255 380 CONTINUE 00256 IF (IRES .LE. -2 .OR. IERTYP .NE. 0) THEN 00257 IERNLS = -1 00258 IF (IRES .LE. -2) IDID = -11 00259 IF (IERTYP .NE. 0) IDID = -15 00260 ELSE 00261 IERNLS = 1 00262 IF (IRES .LT. 0) IDID = -10 00263 IF (IERJ .NE. 0) IDID = -8 00264 ENDIF 00265 C 00266 390 JCALC = 1 00267 RETURN 00268 C 00269 C------END OF SUBROUTINE DNEDD------------------------------------------ 00270 END
1.7.1
