GNU Octave  6.2.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
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
ddwnrm
double precision function ddwnrm(NEQ, V, RWT, RPAR, IPAR)
Definition: ddwnrm.f:6
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:8
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:8