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
|
C Work performed 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 DNSD(X,Y,YPRIME,NEQ,RES,PDUM,WT,RPAR,IPAR,
* DUMSVR,DELTA,E,WM,IWM,CJ,DUMS,DUMR,DUME,EPCON,
* S,CONFAC,TOLNEW,MULDEL,MAXIT,IRES,IDUM,IERNEW)
C
C***BEGIN PROLOGUE DNSD
C***REFER TO DDASPK
C***DATE WRITTEN 891219 (YYMMDD)
C***REVISION DATE 900926 (YYMMDD)
C***REVISION DATE 950126 (YYMMDD)
C
C
C-----------------------------------------------------------------------
C***DESCRIPTION
C
C DNSD solves a nonlinear system of
C algebraic equations of the form
C G(X,Y,YPRIME) = 0 for the unknown Y.
C
C The method used is a modified Newton scheme.
C
C The parameters represent
C
C X -- Independent variable.
C Y -- Solution vector.
C YPRIME -- Derivative of solution vector.
C NEQ -- Number of unknowns.
C RES -- External user-supplied subroutine
C to evaluate the residual. See RES description
C in DDASPK prologue.
C PDUM -- Dummy argument.
C WT -- Vector of weights for error criterion.
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 DUMSVR -- Dummy argument.
C DELTA -- Work vector for DNSD of length NEQ.
C E -- Error accumulation vector for DNSD of length NEQ.
C WM,IWM -- Real and integer arrays storing
C matrix information such as the matrix
C of partial derivatives, permutation
C vector, and various other information.
C CJ -- Parameter always proportional to 1/H (step size).
C DUMS -- Dummy argument.
C DUMR -- Dummy argument.
C DUME -- Dummy argument.
C EPCON -- Tolerance to test for convergence of the Newton
C iteration.
C S -- Used for error convergence tests.
C In the Newton iteration: S = RATE/(1 - RATE),
C where RATE is the estimated rate of convergence
C of the Newton iteration.
C The calling routine passes the initial value
C of S to the Newton iteration.
C CONFAC -- A residual scale factor to improve convergence.
C TOLNEW -- Tolerance on the norm of Newton correction in
C alternative Newton convergence test.
C MULDEL -- A flag indicating whether or not to multiply
C DELTA by CONFAC.
C 0 ==> do not scale DELTA by CONFAC.
C 1 ==> scale DELTA by CONFAC.
C MAXIT -- Maximum allowed number of Newton iterations.
C IRES -- Error flag returned from RES. See RES description
C in DDASPK prologue. If IRES = -1, then IERNEW
C will be set to 1.
C If IRES < -1, then IERNEW will be set to -1.
C IDUM -- Dummy argument.
C IERNEW -- Error flag for Newton iteration.
C 0 ==> Newton iteration converged.
C 1 ==> recoverable error inside Newton iteration.
C -1 ==> unrecoverable error inside Newton iteration.
C
C All arguments with "DUM" in their names are dummy arguments
C which are not used in this routine.
C-----------------------------------------------------------------------
C
C***ROUTINES CALLED
C DSLVD, DDWNRM, RES
C
C***END PROLOGUE DNSD
C
C
IMPLICIT DOUBLE PRECISION(A-H,O-Z)
DIMENSION Y(*),YPRIME(*),WT(*),DELTA(*),E(*)
DIMENSION WM(*),IWM(*), RPAR(*),IPAR(*)
EXTERNAL RES
C
PARAMETER (LNRE=12, LNNI=19)
C
C Initialize Newton counter M and accumulation vector E.
C
M = 0
DO 100 I=1,NEQ
100 E(I)=0.0D0
C
C Corrector loop.
C
300 CONTINUE
IWM(LNNI) = IWM(LNNI) + 1
C
C If necessary, multiply residual by convergence factor.
C
IF (MULDEL .EQ. 1) THEN
DO 320 I = 1,NEQ
320 DELTA(I) = DELTA(I) * CONFAC
ENDIF
C
C Compute a new iterate (back-substitution).
C Store the correction in DELTA.
C
CALL DSLVD(NEQ,DELTA,WM,IWM)
C
C Update Y, E, and YPRIME.
C
DO 340 I=1,NEQ
Y(I)=Y(I)-DELTA(I)
E(I)=E(I)-DELTA(I)
340 YPRIME(I)=YPRIME(I)-CJ*DELTA(I)
C
C Test for convergence of the iteration.
C
DELNRM=DDWNRM(NEQ,DELTA,WT,RPAR,IPAR)
IF (DELNRM .LE. TOLNEW) GO TO 370
IF (M .EQ. 0) THEN
OLDNRM = DELNRM
ELSE
RATE = (DELNRM/OLDNRM)**(1.0D0/M)
IF (RATE .GT. 0.9D0) GO TO 380
S = RATE/(1.0D0 - RATE)
ENDIF
IF (S*DELNRM .LE. EPCON) GO TO 370
C
C The corrector has not yet converged.
C Update M and test whether the
C maximum number of iterations have
C been tried.
C
M=M+1
IF(M.GE.MAXIT) GO TO 380
C
C Evaluate the residual,
C and go back to do another iteration.
C
IWM(LNRE)=IWM(LNRE)+1
CALL RES(X,Y,YPRIME,CJ,DELTA,IRES,RPAR,IPAR)
IF (IRES .LT. 0) GO TO 380
GO TO 300
C
C The iteration has converged.
C
370 RETURN
C
C The iteration has not converged. Set IERNEW appropriately.
C
380 CONTINUE
IF (IRES .LE. -2 ) THEN
IERNEW = -1
ELSE
IERNEW = 1
ENDIF
RETURN
C
C
C------END OF SUBROUTINE DNSD-------------------------------------------
END
|
