GNU Octave  4.2.1
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
dnsk.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 dnsk(X,Y,YPRIME,NEQ,RES,PSOL,WT,RPAR,IPAR,
6  * savr,delta,e,wm,iwm,cj,sqrtn,rsqrtn,eplin,epcon,
7  * s,confac,tolnew,muldel,maxit,ires,iersl,iernew)
8 C
9 C***BEGIN PROLOGUE DNSK
10 C***REFER TO DDASPK
11 C***DATE WRITTEN 891219 (YYMMDD)
12 C***REVISION DATE 900926 (YYMMDD)
13 C***REVISION DATE 950126 (YYMMDD)
14 C
15 C
16 C-----------------------------------------------------------------------
17 C***DESCRIPTION
18 C
19 C DNSK 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 PSOL -- External user-supplied routine to solve
35 C a linear system using preconditioning.
36 C See explanation inside DDASPK.
37 C WT -- Vector of weights for error criterion.
38 C RPAR,IPAR -- Real and integer arrays used for communication
39 C between the calling program and external user
40 C routines. They are not altered within DASPK.
41 C SAVR -- Work vector for DNSK of length NEQ.
42 C DELTA -- Work vector for DNSK of length NEQ.
43 C E -- Error accumulation vector for DNSK of length NEQ.
44 C WM,IWM -- Real and integer arrays storing
45 C matrix information such as the matrix
46 C of partial derivatives, permutation
47 C vector, and various other information.
48 C CJ -- Parameter always proportional to 1/H (step size).
49 C SQRTN -- Square root of NEQ.
50 C RSQRTN -- reciprical of square root of NEQ.
51 C EPLIN -- Tolerance for linear system solver.
52 C EPCON -- Tolerance to test for convergence of the Newton
53 C iteration.
54 C S -- Used for error convergence tests.
55 C In the Newton iteration: S = RATE/(1.D0-RATE),
56 C where RATE is the estimated rate of convergence
57 C of the Newton iteration.
58 C
59 C The closer RATE is to 0., the faster the Newton
60 C iteration is converging; the closer RATE is to 1.,
61 C the slower the Newton iteration is converging.
62 C
63 C The calling routine sends the initial value
64 C of S to the Newton iteration.
65 C CONFAC -- A residual scale factor to improve convergence.
66 C TOLNEW -- Tolerance on the norm of Newton correction in
67 C alternative Newton convergence test.
68 C MULDEL -- A flag indicating whether or not to multiply
69 C DELTA by CONFAC.
70 C 0 ==> do not scale DELTA by CONFAC.
71 C 1 ==> scale DELTA by CONFAC.
72 C MAXIT -- Maximum allowed number of Newton iterations.
73 C IRES -- Error flag returned from RES. See RES description
74 C in DDASPK prologue. If IRES = -1, then IERNEW
75 C will be set to 1.
76 C If IRES < -1, then IERNEW will be set to -1.
77 C IERSL -- Error flag for linear system solver.
78 C See IERSL description in subroutine DSLVK.
79 C If IERSL = 1, then IERNEW will be set to 1.
80 C If IERSL < 0, then IERNEW will be set to -1.
81 C IERNEW -- Error flag for Newton iteration.
82 C 0 ==> Newton iteration converged.
83 C 1 ==> recoverable error inside Newton iteration.
84 C -1 ==> unrecoverable error inside Newton iteration.
85 C-----------------------------------------------------------------------
86 C
87 C***ROUTINES CALLED
88 C RES, DSLVK, DDWNRM
89 C
90 C***END PROLOGUE DNSK
91 C
92 C
93  IMPLICIT DOUBLE PRECISION(a-h,o-z)
94  dimension y(*),yprime(*),wt(*),delta(*),e(*),savr(*)
95  dimension wm(*),iwm(*), rpar(*),ipar(*)
96  EXTERNAL res, psol
97 C
98  parameter(lnni=19, lnre=12)
99 C
100 C Initialize Newton counter M and accumulation vector E.
101 C
102  m = 0
103  DO 100 i=1,neq
104 100 e(i) = 0.0d0
105 C
106 C Corrector loop.
107 C
108 300 CONTINUE
109  iwm(lnni) = iwm(lnni) + 1
110 C
111 C If necessary, multiply residual by convergence factor.
112 C
113  IF (muldel .EQ. 1) THEN
114  DO 320 i = 1,neq
115 320 delta(i) = delta(i) * confac
116  ENDIF
117 C
118 C Save residual in SAVR.
119 C
120  DO 340 i = 1,neq
121 340 savr(i) = delta(i)
122 C
123 C Compute a new iterate. Store the correction in DELTA.
124 C
125  CALL dslvk(neq, y, x, yprime, savr, delta, wt, wm, iwm,
126  * res, ires, psol, iersl, cj, eplin, sqrtn, rsqrtn, rhok,
127  * rpar, ipar)
128  IF (ires .NE. 0 .OR. iersl .NE. 0) go to 380
129 C
130 C Update Y, E, and YPRIME.
131 C
132  DO 360 i=1,neq
133  y(i) = y(i) - delta(i)
134  e(i) = e(i) - delta(i)
135 360 yprime(i) = yprime(i) - cj*delta(i)
136 C
137 C Test for convergence of the iteration.
138 C
139  delnrm = ddwnrm(neq,delta,wt,rpar,ipar)
140  IF (delnrm .LE. tolnew) go to 370
141  IF (m .EQ. 0) THEN
142  oldnrm = delnrm
143  ELSE
144  rate = (delnrm/oldnrm)**(1.0d0/m)
145  IF (rate .GT. 0.9d0) go to 380
146  s = rate/(1.0d0 - rate)
147  ENDIF
148  IF (s*delnrm .LE. epcon) go to 370
149 C
150 C The corrector has not yet converged. Update M and test whether
151 C the maximum number of iterations have been tried.
152 C
153  m = m + 1
154  IF (m .GE. maxit) go to 380
155 C
156 C Evaluate the residual, and go back to do another iteration.
157 C
158  iwm(lnre) = iwm(lnre) + 1
159  CALL res(x,y,yprime,cj,delta,ires,rpar,ipar)
160  IF (ires .LT. 0) go to 380
161  go to 300
162 C
163 C The iteration has converged.
164 C
165 370 RETURN
166 C
167 C The iteration has not converged. Set IERNEW appropriately.
168 C
169 380 CONTINUE
170  IF (ires .LE. -2 .OR. iersl .LT. 0) THEN
171  iernew = -1
172  ELSE
173  iernew = 1
174  ENDIF
175  RETURN
176 C
177 C
178 C------END OF SUBROUTINE DNSK-------------------------------------------
179  END
dslvk
subroutine dslvk(NEQ, Y, TN, YPRIME, SAVR, X, EWT, WM, IWM, RES, IRES, PSOL, IERSL, CJ, EPLIN, SQRTN, RSQRTN, RHOK, RPAR, IPAR)
Definition: dslvk.f:5
dnsk
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
e
i e
Definition: data.cc:2724
h
double h
Definition: graphics.cc:11205
to
may be zero for pure relative error test tem the relative tolerance must be greater than or equal to
Definition: Quad-opts.cc:233
dimension
OCTAVE_EXPORT octave_value_list etc The functions then dimension(columns)
ddwnrm
double precision function ddwnrm(NEQ, V, RWT, RPAR, IPAR)
Definition: ddwnrm.f:5
Generated on Wed May 10 2017 15:42:57 for GNU Octave by   doxygen 1.8.8