1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
|
C Work perfored under the auspices of the U.S. Department of Energy
C by Lawrence Livermore National Laboratory under contract number
C W-7405-Eng-48.
C
SUBROUTINE DDASIK(X,Y,YPRIME,NEQ,ICOPT,ID,RES,JACK,PSOL,H,WT,
* JSKIP,RPAR,IPAR,SAVR,DELTA,R,YIC,YPIC,PWK,WM,IWM,CJ,UROUND,
* EPLI,SQRTN,RSQRTN,EPCON,RATEMX,STPTOL,JFLG,
* ICNFLG,ICNSTR,IERNLS)
C
C***BEGIN PROLOGUE DDASIK
C***REFER TO DDASPK
C***DATE WRITTEN 941026 (YYMMDD)
C***REVISION DATE 950808 (YYMMDD)
C***REVISION DATE 951110 Removed unreachable block 390.
C
C
C-----------------------------------------------------------------------
C***DESCRIPTION
C
C
C DDASIK solves a nonlinear system of algebraic equations of the
C form G(X,Y,YPRIME) = 0 for the unknown parts of Y and YPRIME in
C the initial conditions.
C
C An initial value for Y and initial guess for YPRIME are input.
C
C The method used is a Newton scheme with Krylov iteration and a
C linesearch algorithm.
C
C The parameters represent
C
C X -- Independent variable.
C Y -- Solution vector at x.
C YPRIME -- Derivative of solution vector.
C NEQ -- Number of equations to be integrated.
C ICOPT -- Initial condition option chosen (1 or 2).
C ID -- Array of dimension NEQ, which must be initialized
C if ICOPT = 1. See DDASIC.
C RES -- External user-supplied subroutine
C to evaluate the residual. See RES description
C in DDASPK prologue.
C JACK -- External user-supplied routine to update
C the preconditioner. (This is optional).
C See JAC description for the case
C INFO(12) = 1 in the DDASPK prologue.
C PSOL -- External user-supplied routine to solve
C a linear system using preconditioning.
C (This is optional). See explanation inside DDASPK.
C H -- Scaling factor for this initial condition calc.
C WT -- Vector of weights for error criterion.
C JSKIP -- input flag to signal if initial JAC call is to be
C skipped. 1 => skip the call, 0 => do not skip call.
C RPAR,IPAR -- Real and integer arrays used for communication
C between the calling program and external user
C routines. They are not altered within DASPK.
C SAVR -- Work vector for DDASIK of length NEQ.
C DELTA -- Work vector for DDASIK of length NEQ.
C R -- Work vector for DDASIK of length NEQ.
C YIC,YPIC -- Work vectors for DDASIK, each of length NEQ.
C PWK -- Work vector for DDASIK of length NEQ.
C WM,IWM -- Real and integer arrays storing
C matrix information for linear system
C solvers, and various other information.
C CJ -- Matrix parameter = 1/H (ICOPT = 1) or 0 (ICOPT = 2).
C UROUND -- Unit roundoff.
C EPLI -- convergence test constant.
C See DDASPK prologue for more details.
C SQRTN -- Square root of NEQ.
C RSQRTN -- reciprical of square root of NEQ.
C EPCON -- Tolerance to test for convergence of the Newton
C iteration.
C RATEMX -- Maximum convergence rate for which Newton iteration
C is considered converging.
C JFLG -- Flag showing whether a Jacobian routine is supplied.
C ICNFLG -- Integer scalar. If nonzero, then constraint
C violations in the proposed new approximate solution
C will be checked for, and the maximum step length
C will be adjusted accordingly.
C ICNSTR -- Integer array of length NEQ containing flags for
C checking constraints.
C IERNLS -- Error flag for nonlinear solver.
C 0 ==> nonlinear solver converged.
C 1,2 ==> recoverable error inside nonlinear solver.
C 1 => retry with current Y, YPRIME
C 2 => retry with original Y, YPRIME
C -1 ==> unrecoverable error in nonlinear solver.
C
C-----------------------------------------------------------------------
C
C***ROUTINES CALLED
C RES, JACK, DNSIK, DCOPY
C
C***END PROLOGUE DDASIK
C
C
IMPLICIT DOUBLE PRECISION(A-H,O-Z)
DIMENSION Y(*),YPRIME(*),ID(*),WT(*),ICNSTR(*)
DIMENSION SAVR(*),DELTA(*),R(*),YIC(*),YPIC(*),PWK(*)
DIMENSION WM(*),IWM(*), RPAR(*),IPAR(*)
EXTERNAL RES, JACK, PSOL
C
PARAMETER (LNRE=12, LNJE=13, LLOCWP=29, LLCIWP=30)
PARAMETER (LMXNIT=32, LMXNJ=33)
C
C
C Perform initializations.
C
LWP = IWM(LLOCWP)
LIWP = IWM(LLCIWP)
MXNIT = IWM(LMXNIT)
MXNJ = IWM(LMXNJ)
IERNLS = 0
NJ = 0
EPLIN = EPLI*EPCON
C
C Call RES to initialize DELTA.
C
IRES = 0
IWM(LNRE) = IWM(LNRE) + 1
CALL RES(X,Y,YPRIME,CJ,DELTA,IRES,RPAR,IPAR)
IF (IRES .LT. 0) GO TO 370
C
C Looping point for updating the preconditioner.
C
300 CONTINUE
C
C Initialize all error flags to zero.
C
IERPJ = 0
IRES = 0
IERNEW = 0
C
C If a Jacobian routine was supplied, call it.
C
IF (JFLG .EQ. 1 .AND. JSKIP .EQ. 0) THEN
NJ = NJ + 1
IWM(LNJE)=IWM(LNJE)+1
CALL JACK (RES, IRES, NEQ, X, Y, YPRIME, WT, DELTA, R, H, CJ,
* WM(LWP), IWM(LIWP), IERPJ, RPAR, IPAR)
IF (IRES .LT. 0 .OR. IERPJ .NE. 0) GO TO 370
ENDIF
JSKIP = 0
C
C Call the nonlinear Newton solver for up to MXNIT iterations.
C
CALL DNSIK(X,Y,YPRIME,NEQ,ICOPT,ID,RES,PSOL,WT,RPAR,IPAR,
* SAVR,DELTA,R,YIC,YPIC,PWK,WM,IWM,CJ,SQRTN,RSQRTN,
* EPLIN,EPCON,RATEMX,MXNIT,STPTOL,ICNFLG,ICNSTR,IERNEW)
C
IF (IERNEW .EQ. 1 .AND. NJ .LT. MXNJ .AND. JFLG .EQ. 1) THEN
C
C Up to MXNIT iterations were done, the convergence rate is < 1,
C a Jacobian routine is supplied, and the number of JACK calls
C is less than MXNJ.
C Copy the residual SAVR to DELTA, call JACK, and try again.
C
CALL DCOPY (NEQ, SAVR, 1, DELTA, 1)
GO TO 300
ENDIF
C
IF (IERNEW .NE. 0) GO TO 380
RETURN
C
C
C Unsuccessful exits from nonlinear solver.
C Set IERNLS accordingly.
C
370 IERNLS = 2
IF (IRES .LE. -2) IERNLS = -1
RETURN
C
380 IERNLS = MIN(IERNEW,2)
RETURN
C
C----------------------- END OF SUBROUTINE DDASIK-----------------------
END
|
