GNU Octave  6.2.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
dnedd.f
Go to the documentation of this file.
1 C Work performed under the auspices of the U.S. Department of Energy
2 C by Lawrence Livermore National Laboratory under contract number
3 C W-7405-Eng-48.
4 C
5  SUBROUTINE dnedd(X,Y,YPRIME,NEQ,RES,JACD,PDUM,H,WT,
6  * JSTART,IDID,RPAR,IPAR,PHI,GAMMA,DUMSVR,DELTA,E,
7  * WM,IWM,CJ,CJOLD,CJLAST,S,UROUND,DUME,DUMS,DUMR,
8  * EPCON,JCALC,JFDUM,KP1,NONNEG,NTYPE,IERNLS)
9 C
10 C***BEGIN PROLOGUE DNEDD
11 C***REFER TO DDASPK
12 C***DATE WRITTEN 891219 (YYMMDD)
13 C***REVISION DATE 900926 (YYMMDD)
14 C
15 C
16 C-----------------------------------------------------------------------
17 C***DESCRIPTION
18 C
19 C DNEDD solves a nonlinear system of
20 C algebraic equations of the form
21 C G(X,Y,YPRIME) = 0 for the unknown Y.
22 C
23 C The method used is a modified Newton scheme.
24 C
25 C The parameters represent
26 C
27 C X -- Independent variable.
28 C Y -- Solution vector.
29 C YPRIME -- Derivative of solution vector.
30 C NEQ -- Number of unknowns.
31 C RES -- External user-supplied subroutine
32 C to evaluate the residual. See RES description
33 C in DDASPK prologue.
34 C JACD -- External user-supplied routine to evaluate the
35 C Jacobian. See JAC description for the case
36 C INFO(12) = 0 in the DDASPK prologue.
37 C PDUM -- Dummy argument.
38 C H -- Appropriate step size for next step.
39 C WT -- Vector of weights for error criterion.
40 C JSTART -- Indicates first call to this routine.
41 C If JSTART = 0, then this is the first call,
42 C otherwise it is not.
43 C IDID -- Completion flag, output by DNEDD.
44 C See IDID description in DDASPK prologue.
45 C RPAR,IPAR -- Real and integer arrays used for communication
46 C between the calling program and external user
47 C routines. They are not altered within DASPK.
48 C PHI -- Array of divided differences used by
49 C DNEDD. The length is NEQ*(K+1),where
50 C K is the maximum order.
51 C GAMMA -- Array used to predict Y and YPRIME. The length
52 C is MAXORD+1 where MAXORD is the maximum order.
53 C DUMSVR -- Dummy argument.
54 C DELTA -- Work vector for NLS of length NEQ.
55 C E -- Error accumulation vector for NLS of length NEQ.
56 C WM,IWM -- Real and integer arrays storing
57 C matrix information such as the matrix
58 C of partial derivatives, permutation
59 C vector, and various other information.
60 C CJ -- Parameter always proportional to 1/H.
61 C CJOLD -- Saves the value of CJ as of the last call to DMATD.
62 C Accounts for changes in CJ needed to
63 C decide whether to call DMATD.
64 C CJLAST -- Previous value of CJ.
65 C S -- A scalar determined by the approximate rate
66 C of convergence of the Newton iteration and used
67 C in the convergence test for the Newton iteration.
68 C
69 C If RATE is defined to be an estimate of the
70 C rate of convergence of the Newton iteration,
71 C then S = RATE/(1.D0-RATE).
72 C
73 C The closer RATE is to 0., the faster the Newton
74 C iteration is converging; the closer RATE is to 1.,
75 C the slower the Newton iteration is converging.
76 C
77 C On the first Newton iteration with an up-dated
78 C preconditioner S = 100.D0, Thus the initial
79 C RATE of convergence is approximately 1.
80 C
81 C S is preserved from call to call so that the rate
82 C estimate from a previous step can be applied to
83 C the current step.
84 C UROUND -- Unit roundoff.
85 C DUME -- Dummy argument.
86 C DUMS -- Dummy argument.
87 C DUMR -- Dummy argument.
88 C EPCON -- Tolerance to test for convergence of the Newton
89 C iteration.
90 C JCALC -- Flag used to determine when to update
91 C the Jacobian matrix. In general:
92 C
93 C JCALC = -1 ==> Call the DMATD routine to update
94 C the Jacobian matrix.
95 C JCALC = 0 ==> Jacobian matrix is up-to-date.
96 C JCALC = 1 ==> Jacobian matrix is out-dated,
97 C but DMATD will not be called unless
98 C JCALC is set to -1.
99 C JFDUM -- Dummy argument.
100 C KP1 -- The current order(K) + 1; updated across calls.
101 C NONNEG -- Flag to determine nonnegativity constraints.
102 C NTYPE -- Identification code for the NLS routine.
103 C 0 ==> modified Newton; direct solver.
104 C IERNLS -- Error flag for nonlinear solver.
105 C 0 ==> nonlinear solver converged.
106 C 1 ==> recoverable error inside nonlinear solver.
107 C -1 ==> unrecoverable error inside nonlinear solver.
108 C
109 C All variables with "DUM" in their names are dummy variables
110 C which are not used in this routine.
111 C
112 C Following is a list and description of local variables which
113 C may not have an obvious usage. They are listed in roughly the
114 C order they occur in this subroutine.
115 C
116 C The following group of variables are passed as arguments to
117 C the Newton iteration solver. They are explained in greater detail
118 C in DNSD:
119 C TOLNEW, MULDEL, MAXIT, IERNEW
120 C
121 C IERTYP -- Flag which tells whether this subroutine is correct.
122 C 0 ==> correct subroutine.
123 C 1 ==> incorrect subroutine.
124 C
125 C-----------------------------------------------------------------------
126 C***ROUTINES CALLED
127 C DDWNRM, RES, DMATD, DNSD
128 C
129 C***END PROLOGUE DNEDD
130 C
131 C
132  IMPLICIT DOUBLE PRECISION(a-h,o-z)
133  dimension y(*),yprime(*),wt(*)
134  dimension delta(*),e(*)
135  dimension wm(*),iwm(*), rpar(*),ipar(*)
136  dimension phi(neq,*),gamma(*)
137  EXTERNAL res, jacd
138 C
139  parameter(lnre=12, lnje=13)
140 C
141  SAVE muldel, maxit, xrate
142  DATA muldel/1/, maxit/4/, xrate/0.25d0/
143 C
144 C Verify that this is the correct subroutine.
145 C
146  iertyp = 0
147  IF (ntype .NE. 0) THEN
148  iertyp = 1
149  GO TO 380
150  ENDIF
151 C
152 C If this is the first step, perform initializations.
153 C
154  IF (jstart .EQ. 0) THEN
155  cjold = cj
156  jcalc = -1
157  ENDIF
158 C
159 C Perform all other initializations.
160 C
161  iernls = 0
162 C
163 C Decide whether new Jacobian is needed.
164 C
165  temp1 = (1.0d0 - xrate)/(1.0d0 + xrate)
166  temp2 = 1.0d0/temp1
167  IF (cj/cjold .LT. temp1 .OR. cj/cjold .GT. temp2) jcalc = -1
168  IF (cj .NE. cjlast) s = 100.d0
169 C
170 C-----------------------------------------------------------------------
171 C Entry point for updating the Jacobian with current
172 C stepsize.
173 C-----------------------------------------------------------------------
174 300 CONTINUE
175 C
176 C Initialize all error flags to zero.
177 C
178  ierj = 0
179  ires = 0
180  iernew = 0
181 C
182 C Predict the solution and derivative and compute the tolerance
183 C for the Newton iteration.
184 C
185  DO 310 i=1,neq
186  y(i)=phi(i,1)
187 310 yprime(i)=0.0d0
188  DO 330 j=2,kp1
189  DO 320 i=1,neq
190  y(i)=y(i)+phi(i,j)
191 320 yprime(i)=yprime(i)+gamma(j)*phi(i,j)
192 330 CONTINUE
193  pnorm = ddwnrm(neq,y,wt,rpar,ipar)
194  tolnew = 100.d0*uround*pnorm
195 C
196 C Call RES to initialize DELTA.
197 C
198  iwm(lnre)=iwm(lnre)+1
199  CALL res(x,y,yprime,cj,delta,ires,rpar,ipar)
200  IF (ires .LT. 0) GO TO 380
201 C
202 C If indicated, reevaluate the iteration matrix
203 C J = dG/dY + CJ*dG/dYPRIME (where G(X,Y,YPRIME)=0).
204 C Set JCALC to 0 as an indicator that this has been done.
205 C
206  IF(jcalc .EQ. -1) THEN
207  iwm(lnje)=iwm(lnje)+1
208  jcalc=0
209  CALL dmatd(neq,x,y,yprime,delta,cj,h,ierj,wt,e,wm,iwm,
210  * res,ires,uround,jacd,rpar,ipar)
211  cjold=cj
212  s = 100.d0
213  IF (ires .LT. 0) GO TO 380
214  IF(ierj .NE. 0)GO TO 380
215  ENDIF
216 C
217 C Call the nonlinear Newton solver.
218 C
219  temp1 = 2.0d0/(1.0d0 + cj/cjold)
220  CALL dnsd(x,y,yprime,neq,res,pdum,wt,rpar,ipar,dumsvr,
221  * delta,e,wm,iwm,cj,dums,dumr,dume,epcon,s,temp1,
222  * tolnew,muldel,maxit,ires,idum,iernew)
223 C
224  IF (iernew .GT. 0 .AND. jcalc .NE. 0) THEN
225 C
226 C The Newton iteration had a recoverable failure with an old
227 C iteration matrix. Retry the step with a new iteration matrix.
228 C
229  jcalc = -1
230  GO TO 300
231  ENDIF
232 C
233  IF (iernew .NE. 0) GO TO 380
234 C
235 C The Newton iteration has converged. If nonnegativity of
236 C solution is required, set the solution nonnegative, if the
237 C perturbation to do it is small enough. If the change is too
238 C large, then consider the corrector iteration to have failed.
239 C
240 375 IF(nonneg .EQ. 0) GO TO 390
241  DO 377 i = 1,neq
242 377 delta(i) = min(y(i),0.0d0)
243  delnrm = ddwnrm(neq,delta,wt,rpar,ipar)
244  IF(delnrm .GT. epcon) GO TO 380
245  DO 378 i = 1,neq
246 378 e(i) = e(i) - delta(i)
247  GO TO 390
248 C
249 C
250 C Exits from nonlinear solver.
251 C No convergence with current iteration
252 C matrix, or singular iteration matrix.
253 C Compute IERNLS and IDID accordingly.
254 C
255 380 CONTINUE
256  IF (ires .LE. -2 .OR. iertyp .NE. 0) THEN
257  iernls = -1
258  IF (ires .LE. -2) idid = -11
259  IF (iertyp .NE. 0) idid = -15
260  ELSE
261  iernls = 1
262  IF (ires .LT. 0) idid = -10
263  IF (ierj .NE. 0) idid = -8
264  ENDIF
265 C
266 390 jcalc = 1
267  RETURN
268 C
269 C------END OF SUBROUTINE DNEDD------------------------------------------
270  END
charNDArray min(char d, const charNDArray &m)
Definition: chNDArray.cc:207
double precision function ddwnrm(NEQ, V, RWT, RPAR, IPAR)
Definition: ddwnrm.f:6
subroutine dmatd(NEQ, X, Y, YPRIME, DELTA, CJ, H, IER, EWT, E, WM, IWM, RES, IRES, UROUND, JACD, RPAR, IPAR)
Definition: dmatd.f:7
subroutine dnedd(X, Y, YPRIME, NEQ, RES, JACD, PDUM, H, WT, JSTART, IDID, RPAR, IPAR, PHI, GAMMA, DUMSVR, DELTA, E, WM, IWM, CJ, CJOLD, CJLAST, S, UROUND, DUME, DUMS, DUMR, EPCON, JCALC, JFDUM, KP1, NONNEG, NTYPE, IERNLS)
Definition: dnedd.f:9
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)
Definition: dnsd.f:8
function gamma(X)
Definition: gamma.f:3