00001 C Work perfored 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 DDASIK(X,Y,YPRIME,NEQ,ICOPT,ID,RES,JACK,PSOL,H,WT, 00006 * JSKIP,RPAR,IPAR,SAVR,DELTA,R,YIC,YPIC,PWK,WM,IWM,CJ,UROUND, 00007 * EPLI,SQRTN,RSQRTN,EPCON,RATEMX,STPTOL,JFLG, 00008 * ICNFLG,ICNSTR,IERNLS) 00009 C 00010 C***BEGIN PROLOGUE DDASIK 00011 C***REFER TO DDASPK 00012 C***DATE WRITTEN 941026 (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 DDASIK 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 An initial value for Y and initial guess for YPRIME are input. 00026 C 00027 C The method used is a Newton scheme with Krylov iteration and a 00028 C linesearch algorithm. 00029 C 00030 C The parameters represent 00031 C 00032 C X -- Independent variable. 00033 C Y -- Solution vector at x. 00034 C YPRIME -- Derivative of solution vector. 00035 C NEQ -- Number of equations to be integrated. 00036 C ICOPT -- Initial condition option chosen (1 or 2). 00037 C ID -- Array of dimension NEQ, which must be initialized 00038 C if ICOPT = 1. See DDASIC. 00039 C RES -- External user-supplied subroutine 00040 C to evaluate the residual. See RES description 00041 C in DDASPK prologue. 00042 C JACK -- External user-supplied routine to update 00043 C the preconditioner. (This is optional). 00044 C See JAC description for the case 00045 C INFO(12) = 1 in the DDASPK prologue. 00046 C PSOL -- External user-supplied routine to solve 00047 C a linear system using preconditioning. 00048 C (This is optional). See explanation inside DDASPK. 00049 C H -- Scaling factor for this initial condition calc. 00050 C WT -- Vector of weights for error criterion. 00051 C JSKIP -- input flag to signal if initial JAC call is to be 00052 C skipped. 1 => skip the call, 0 => do not skip call. 00053 C RPAR,IPAR -- Real and integer arrays used for communication 00054 C between the calling program and external user 00055 C routines. They are not altered within DASPK. 00056 C SAVR -- Work vector for DDASIK of length NEQ. 00057 C DELTA -- Work vector for DDASIK of length NEQ. 00058 C R -- Work vector for DDASIK of length NEQ. 00059 C YIC,YPIC -- Work vectors for DDASIK, each of length NEQ. 00060 C PWK -- Work vector for DDASIK of length NEQ. 00061 C WM,IWM -- Real and integer arrays storing 00062 C matrix information for linear system 00063 C solvers, and various other information. 00064 C CJ -- Matrix parameter = 1/H (ICOPT = 1) or 0 (ICOPT = 2). 00065 C UROUND -- Unit roundoff. 00066 C EPLI -- convergence test constant. 00067 C See DDASPK prologue for more details. 00068 C SQRTN -- Square root of NEQ. 00069 C RSQRTN -- reciprical of square root of NEQ. 00070 C EPCON -- Tolerance to test for convergence of the Newton 00071 C iteration. 00072 C RATEMX -- Maximum convergence rate for which Newton iteration 00073 C is considered converging. 00074 C JFLG -- Flag showing whether a Jacobian routine is supplied. 00075 C ICNFLG -- Integer scalar. If nonzero, then constraint 00076 C violations in the proposed new approximate solution 00077 C will be checked for, and the maximum step length 00078 C will be adjusted accordingly. 00079 C ICNSTR -- Integer array of length NEQ containing flags for 00080 C checking constraints. 00081 C IERNLS -- Error flag for nonlinear solver. 00082 C 0 ==> nonlinear solver converged. 00083 C 1,2 ==> recoverable error inside nonlinear solver. 00084 C 1 => retry with current Y, YPRIME 00085 C 2 => retry with original Y, YPRIME 00086 C -1 ==> unrecoverable error in nonlinear solver. 00087 C 00088 C----------------------------------------------------------------------- 00089 C 00090 C***ROUTINES CALLED 00091 C RES, JACK, DNSIK, DCOPY 00092 C 00093 C***END PROLOGUE DDASIK 00094 C 00095 C 00096 IMPLICIT DOUBLE PRECISION(A-H,O-Z) 00097 DIMENSION Y(*),YPRIME(*),ID(*),WT(*),ICNSTR(*) 00098 DIMENSION SAVR(*),DELTA(*),R(*),YIC(*),YPIC(*),PWK(*) 00099 DIMENSION WM(*),IWM(*), RPAR(*),IPAR(*) 00100 EXTERNAL RES, JACK, PSOL 00101 C 00102 PARAMETER (LNRE=12, LNJE=13, LLOCWP=29, LLCIWP=30) 00103 PARAMETER (LMXNIT=32, LMXNJ=33) 00104 C 00105 C 00106 C Perform initializations. 00107 C 00108 LWP = IWM(LLOCWP) 00109 LIWP = IWM(LLCIWP) 00110 MXNIT = IWM(LMXNIT) 00111 MXNJ = IWM(LMXNJ) 00112 IERNLS = 0 00113 NJ = 0 00114 EPLIN = EPLI*EPCON 00115 C 00116 C Call RES to initialize DELTA. 00117 C 00118 IRES = 0 00119 IWM(LNRE) = IWM(LNRE) + 1 00120 CALL RES(X,Y,YPRIME,CJ,DELTA,IRES,RPAR,IPAR) 00121 IF (IRES .LT. 0) GO TO 370 00122 C 00123 C Looping point for updating the preconditioner. 00124 C 00125 300 CONTINUE 00126 C 00127 C Initialize all error flags to zero. 00128 C 00129 IERPJ = 0 00130 IRES = 0 00131 IERNEW = 0 00132 C 00133 C If a Jacobian routine was supplied, call it. 00134 C 00135 IF (JFLG .EQ. 1 .AND. JSKIP .EQ. 0) THEN 00136 NJ = NJ + 1 00137 IWM(LNJE)=IWM(LNJE)+1 00138 CALL JACK (RES, IRES, NEQ, X, Y, YPRIME, WT, DELTA, R, H, CJ, 00139 * WM(LWP), IWM(LIWP), IERPJ, RPAR, IPAR) 00140 IF (IRES .LT. 0 .OR. IERPJ .NE. 0) GO TO 370 00141 ENDIF 00142 JSKIP = 0 00143 C 00144 C Call the nonlinear Newton solver for up to MXNIT iterations. 00145 C 00146 CALL DNSIK(X,Y,YPRIME,NEQ,ICOPT,ID,RES,PSOL,WT,RPAR,IPAR, 00147 * SAVR,DELTA,R,YIC,YPIC,PWK,WM,IWM,CJ,SQRTN,RSQRTN, 00148 * EPLIN,EPCON,RATEMX,MXNIT,STPTOL,ICNFLG,ICNSTR,IERNEW) 00149 C 00150 IF (IERNEW .EQ. 1 .AND. NJ .LT. MXNJ .AND. JFLG .EQ. 1) THEN 00151 C 00152 C Up to MXNIT iterations were done, the convergence rate is < 1, 00153 C a Jacobian routine is supplied, and the number of JACK calls 00154 C is less than MXNJ. 00155 C Copy the residual SAVR to DELTA, call JACK, and try again. 00156 C 00157 CALL DCOPY (NEQ, SAVR, 1, DELTA, 1) 00158 GO TO 300 00159 ENDIF 00160 C 00161 IF (IERNEW .NE. 0) GO TO 380 00162 RETURN 00163 C 00164 C 00165 C Unsuccessful exits from nonlinear solver. 00166 C Set IERNLS accordingly. 00167 C 00168 370 IERNLS = 2 00169 IF (IRES .LE. -2) IERNLS = -1 00170 RETURN 00171 C 00172 380 IERNLS = MIN(IERNEW,2) 00173 RETURN 00174 C 00175 C----------------------- END OF SUBROUTINE DDASIK----------------------- 00176 END
1.7.1
