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