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 DDASID(X,Y,YPRIME,NEQ,ICOPT,ID,RES,JACD,PDUM,H,WT, 00006 * JSDUM,RPAR,IPAR,DUMSVR,DELTA,R,YIC,YPIC,DUMPWK,WM,IWM,CJ,UROUND, 00007 * DUME,DUMS,DUMR,EPCON,RATEMX,STPTOL,JFDUM, 00008 * ICNFLG,ICNSTR,IERNLS) 00009 C 00010 C***BEGIN PROLOGUE DDASID 00011 C***REFER TO DDASPK 00012 C***DATE WRITTEN 940701 (YYMMDD) 00013 C***REVISION DATE 950808 (YYMMDD) 00014 C***REVISION DATE 951110 Removed unreachable block 390. 00015 C 00016 C 00017 C----------------------------------------------------------------------- 00018 C***DESCRIPTION 00019 C 00020 C 00021 C DDASID solves a nonlinear system of algebraic equations of the 00022 C form G(X,Y,YPRIME) = 0 for the unknown parts of Y and YPRIME in 00023 C the initial conditions. 00024 C 00025 C The method used is a modified Newton scheme. 00026 C 00027 C The parameters represent 00028 C 00029 C X -- Independent variable. 00030 C Y -- Solution vector. 00031 C YPRIME -- Derivative of solution vector. 00032 C NEQ -- Number of unknowns. 00033 C ICOPT -- Initial condition option chosen (1 or 2). 00034 C ID -- Array of dimension NEQ, which must be initialized 00035 C if ICOPT = 1. See DDASIC. 00036 C RES -- External user-supplied subroutine to evaluate the 00037 C residual. See RES description in DDASPK prologue. 00038 C JACD -- External user-supplied routine to evaluate the 00039 C Jacobian. See JAC description for the case 00040 C INFO(12) = 0 in the DDASPK prologue. 00041 C PDUM -- Dummy argument. 00042 C H -- Scaling factor for this initial condition calc. 00043 C WT -- Vector of weights for error criterion. 00044 C JSDUM -- Dummy argument. 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 DUMSVR -- Dummy argument. 00049 C DELTA -- Work vector for NLS of length NEQ. 00050 C R -- Work vector for NLS of length NEQ. 00051 C YIC,YPIC -- Work vectors for NLS, each of length NEQ. 00052 C DUMPWK -- Dummy argument. 00053 C WM,IWM -- Real and integer arrays storing matrix information 00054 C such as the matrix of partial derivatives, 00055 C permutation vector, and various other information. 00056 C CJ -- Matrix parameter = 1/H (ICOPT = 1) or 0 (ICOPT = 2). 00057 C UROUND -- Unit roundoff. 00058 C DUME -- Dummy argument. 00059 C DUMS -- Dummy argument. 00060 C DUMR -- Dummy argument. 00061 C EPCON -- Tolerance to test for convergence of the Newton 00062 C iteration. 00063 C RATEMX -- Maximum convergence rate for which Newton iteration 00064 C is considered converging. 00065 C JFDUM -- Dummy argument. 00066 C STPTOL -- Tolerance used in calculating the minimum lambda 00067 C value allowed. 00068 C ICNFLG -- Integer scalar. If nonzero, then constraint 00069 C violations in the proposed new approximate solution 00070 C will be checked for, and the maximum step length 00071 C will be adjusted accordingly. 00072 C ICNSTR -- Integer array of length NEQ containing flags for 00073 C checking constraints. 00074 C IERNLS -- Error flag for nonlinear solver. 00075 C 0 ==> nonlinear solver converged. 00076 C 1,2 ==> recoverable error inside nonlinear solver. 00077 C 1 => retry with current Y, YPRIME 00078 C 2 => retry with original Y, YPRIME 00079 C -1 ==> unrecoverable error in nonlinear solver. 00080 C 00081 C All variables with "DUM" in their names are dummy variables 00082 C which are not used in this routine. 00083 C 00084 C----------------------------------------------------------------------- 00085 C 00086 C***ROUTINES CALLED 00087 C RES, DMATD, DNSID 00088 C 00089 C***END PROLOGUE DDASID 00090 C 00091 C 00092 IMPLICIT DOUBLE PRECISION(A-H,O-Z) 00093 DIMENSION Y(*),YPRIME(*),ID(*),WT(*),ICNSTR(*) 00094 DIMENSION DELTA(*),R(*),YIC(*),YPIC(*) 00095 DIMENSION WM(*),IWM(*), RPAR(*),IPAR(*) 00096 EXTERNAL RES, JACD 00097 C 00098 PARAMETER (LNRE=12, LNJE=13, LMXNIT=32, LMXNJ=33) 00099 C 00100 C 00101 C Perform initializations. 00102 C 00103 MXNIT = IWM(LMXNIT) 00104 MXNJ = IWM(LMXNJ) 00105 IERNLS = 0 00106 NJ = 0 00107 C 00108 C Call RES to initialize DELTA. 00109 C 00110 IRES = 0 00111 IWM(LNRE) = IWM(LNRE) + 1 00112 CALL RES(X,Y,YPRIME,CJ,DELTA,IRES,RPAR,IPAR) 00113 IF (IRES .LT. 0) GO TO 370 00114 C 00115 C Looping point for updating the Jacobian. 00116 C 00117 300 CONTINUE 00118 C 00119 C Initialize all error flags to zero. 00120 C 00121 IERJ = 0 00122 IRES = 0 00123 IERNEW = 0 00124 C 00125 C Reevaluate the iteration matrix, J = dG/dY + CJ*dG/dYPRIME, 00126 C where G(X,Y,YPRIME) = 0. 00127 C 00128 NJ = NJ + 1 00129 IWM(LNJE)=IWM(LNJE)+1 00130 CALL DMATD(NEQ,X,Y,YPRIME,DELTA,CJ,H,IERJ,WT,R, 00131 * WM,IWM,RES,IRES,UROUND,JACD,RPAR,IPAR) 00132 IF (IRES .LT. 0 .OR. IERJ .NE. 0) GO TO 370 00133 C 00134 C Call the nonlinear Newton solver for up to MXNIT iterations. 00135 C 00136 CALL DNSID(X,Y,YPRIME,NEQ,ICOPT,ID,RES,WT,RPAR,IPAR,DELTA,R, 00137 * YIC,YPIC,WM,IWM,CJ,EPCON,RATEMX,MXNIT,STPTOL, 00138 * ICNFLG,ICNSTR,IERNEW) 00139 C 00140 IF (IERNEW .EQ. 1 .AND. NJ .LT. MXNJ) THEN 00141 C 00142 C MXNIT iterations were done, the convergence rate is < 1, 00143 C and the number of Jacobian evaluations is less than MXNJ. 00144 C Call RES, reevaluate the Jacobian, and try again. 00145 C 00146 IWM(LNRE)=IWM(LNRE)+1 00147 CALL RES(X,Y,YPRIME,CJ,DELTA,IRES,RPAR,IPAR) 00148 IF (IRES .LT. 0) GO TO 370 00149 GO TO 300 00150 ENDIF 00151 C 00152 IF (IERNEW .NE. 0) GO TO 380 00153 00154 RETURN 00155 C 00156 C 00157 C Unsuccessful exits from nonlinear solver. 00158 C Compute IERNLS accordingly. 00159 C 00160 370 IERNLS = 2 00161 IF (IRES .LE. -2) IERNLS = -1 00162 RETURN 00163 C 00164 380 IERNLS = MIN(IERNEW,2) 00165 RETURN 00166 C 00167 C------END OF SUBROUTINE DDASID----------------------------------------- 00168 END
1.7.1
