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