1C Work performed under the auspices of the U.S. Department of Energy 2C by Lawrence Livermore National Laboratory under contract number 3C W-7405-Eng-48. 4C 5 SUBROUTINE DDASID(X,Y,YPRIME,NEQ,ICOPT,ID,RES,JACD,PDUM,H,WT, 6 * JSDUM,RPAR,IPAR,DUMSVR,DELTA,R,YIC,YPIC,DUMPWK,WM,IWM,CJ,UROUND, 7 * DUME,DUMS,DUMR,EPCON,RATEMX,STPTOL,JFDUM, 8 * ICNFLG,ICNSTR,IERNLS) 9C 10C***BEGIN PROLOGUE DDASID 11C***REFER TO DDASPK 12C***DATE WRITTEN 940701 (YYMMDD) 13C***REVISION DATE 950808 (YYMMDD) 14C***REVISION DATE 951110 Removed unreachable block 390. 15C 16C 17C----------------------------------------------------------------------- 18C***DESCRIPTION 19C 20C 21C DDASID solves a nonlinear system of algebraic equations of the 22C form G(X,Y,YPRIME) = 0 for the unknown parts of Y and YPRIME in 23C the initial conditions. 24C 25C The method used is a modified Newton scheme. 26C 27C The parameters represent 28C 29C X -- Independent variable. 30C Y -- Solution vector. 31C YPRIME -- Derivative of solution vector. 32C NEQ -- Number of unknowns. 33C ICOPT -- Initial condition option chosen (1 or 2). 34C ID -- Array of dimension NEQ, which must be initialized 35C if ICOPT = 1. See DDASIC. 36C RES -- External user-supplied subroutine to evaluate the 37C residual. See RES description in DDASPK prologue. 38C JACD -- External user-supplied routine to evaluate the 39C Jacobian. See JAC description for the case 40C INFO(12) = 0 in the DDASPK prologue. 41C PDUM -- Dummy argument. 42C H -- Scaling factor for this initial condition calc. 43C WT -- Vector of weights for error criterion. 44C JSDUM -- Dummy argument. 45C RPAR,IPAR -- Real and integer arrays used for communication 46C between the calling program and external user 47C routines. They are not altered within DASPK. 48C DUMSVR -- Dummy argument. 49C DELTA -- Work vector for NLS of length NEQ. 50C R -- Work vector for NLS of length NEQ. 51C YIC,YPIC -- Work vectors for NLS, each of length NEQ. 52C DUMPWK -- Dummy argument. 53C WM,IWM -- Real and integer arrays storing matrix information 54C such as the matrix of partial derivatives, 55C permutation vector, and various other information. 56C CJ -- Matrix parameter = 1/H (ICOPT = 1) or 0 (ICOPT = 2). 57C UROUND -- Unit roundoff. 58C DUME -- Dummy argument. 59C DUMS -- Dummy argument. 60C DUMR -- Dummy argument. 61C EPCON -- Tolerance to test for convergence of the Newton 62C iteration. 63C RATEMX -- Maximum convergence rate for which Newton iteration 64C is considered converging. 65C JFDUM -- Dummy argument. 66C STPTOL -- Tolerance used in calculating the minimum lambda 67C value allowed. 68C ICNFLG -- Integer scalar. If nonzero, then constraint 69C violations in the proposed new approximate solution 70C will be checked for, and the maximum step length 71C will be adjusted accordingly. 72C ICNSTR -- Integer array of length NEQ containing flags for 73C checking constraints. 74C IERNLS -- Error flag for nonlinear solver. 75C 0 ==> nonlinear solver converged. 76C 1,2 ==> recoverable error inside nonlinear solver. 77C 1 => retry with current Y, YPRIME 78C 2 => retry with original Y, YPRIME 79C -1 ==> unrecoverable error in nonlinear solver. 80C 81C All variables with "DUM" in their names are dummy variables 82C which are not used in this routine. 83C 84C----------------------------------------------------------------------- 85C 86C***ROUTINES CALLED 87C RES, DMATD, DNSID 88C 89C***END PROLOGUE DDASID 90C 91C 92 IMPLICIT DOUBLE PRECISION(A-H,O-Z) 93 DIMENSION Y(*),YPRIME(*),ID(*),WT(*),ICNSTR(*) 94 DIMENSION DELTA(*),R(*),YIC(*),YPIC(*) 95 DIMENSION WM(*),IWM(*), RPAR(*),IPAR(*) 96 EXTERNAL RES, JACD 97C 98 PARAMETER (LNRE=12, LNJE=13, LMXNIT=32, LMXNJ=33) 99C 100C 101C Perform initializations. 102C 103 MXNIT = IWM(LMXNIT) 104 MXNJ = IWM(LMXNJ) 105 IERNLS = 0 106 NJ = 0 107C 108C Call RES to initialize DELTA. 109C 110 IRES = 0 111 IWM(LNRE) = IWM(LNRE) + 1 112 CALL RES(X,Y,YPRIME,CJ,DELTA,IRES,RPAR,IPAR) 113 IF (IRES .LT. 0) GO TO 370 114C 115C Looping point for updating the Jacobian. 116C 117300 CONTINUE 118C 119C Initialize all error flags to zero. 120C 121 IERJ = 0 122 IRES = 0 123 IERNEW = 0 124C 125C Reevaluate the iteration matrix, J = dG/dY + CJ*dG/dYPRIME, 126C where G(X,Y,YPRIME) = 0. 127C 128 NJ = NJ + 1 129 IWM(LNJE)=IWM(LNJE)+1 130 CALL DMATD(NEQ,X,Y,YPRIME,DELTA,CJ,H,IERJ,WT,R, 131 * WM,IWM,RES,IRES,UROUND,JACD,RPAR,IPAR) 132 IF (IRES .LT. 0 .OR. IERJ .NE. 0) GO TO 370 133C 134C Call the nonlinear Newton solver for up to MXNIT iterations. 135C 136 CALL DNSID(X,Y,YPRIME,NEQ,ICOPT,ID,RES,WT,RPAR,IPAR,DELTA,R, 137 * YIC,YPIC,WM,IWM,CJ,EPCON,RATEMX,MXNIT,STPTOL, 138 * ICNFLG,ICNSTR,IERNEW) 139C 140 IF (IERNEW .EQ. 1 .AND. NJ .LT. MXNJ) THEN 141C 142C MXNIT iterations were done, the convergence rate is < 1, 143C and the number of Jacobian evaluations is less than MXNJ. 144C Call RES, reevaluate the Jacobian, and try again. 145C 146 IWM(LNRE)=IWM(LNRE)+1 147 CALL RES(X,Y,YPRIME,CJ,DELTA,IRES,RPAR,IPAR) 148 IF (IRES .LT. 0) GO TO 370 149 GO TO 300 150 ENDIF 151C 152 IF (IERNEW .NE. 0) GO TO 380 153 154 RETURN 155C 156C 157C Unsuccessful exits from nonlinear solver. 158C Compute IERNLS accordingly. 159C 160370 IERNLS = 2 161 IF (IRES .LE. -2) IERNLS = -1 162 RETURN 163C 164380 IERNLS = MIN(IERNEW,2) 165 RETURN 166C 167C------END OF SUBROUTINE DDASID----------------------------------------- 168 END 169