GNU Octave  4.4.1
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
OCTAVE_EXPORT octave_value_list etc The functions then dimension(columns)
charNDArray min(char d, const charNDArray &m)
Definition: chNDArray.cc:204