GNU Octave  9.1.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
dnsd.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 dnsd(X,Y,YPRIME,NEQ,RES,PDUM,WT,RPAR,IPAR,
6  * DUMSVR,DELTA,E,WM,IWM,CJ,DUMS,DUMR,DUME,EPCON,
7  * S,CONFAC,TOLNEW,MULDEL,MAXIT,IRES,IDUM,IERNEW)
8 C
9 C***BEGIN PROLOGUE DNSD
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 DNSD 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 PDUM -- Dummy argument.
35 C WT -- Vector of weights for error criterion.
36 C RPAR,IPAR -- Real and integer arrays used for communication
37 C between the calling program and external user
38 C routines. They are not altered within DASPK.
39 C DUMSVR -- Dummy argument.
40 C DELTA -- Work vector for DNSD of length NEQ.
41 C E -- Error accumulation vector for DNSD of length NEQ.
42 C WM,IWM -- Real and integer arrays storing
43 C matrix information such as the matrix
44 C of partial derivatives, permutation
45 C vector, and various other information.
46 C CJ -- Parameter always proportional to 1/H (step size).
47 C DUMS -- Dummy argument.
48 C DUMR -- Dummy argument.
49 C DUME -- Dummy argument.
50 C EPCON -- Tolerance to test for convergence of the Newton
51 C iteration.
52 C S -- Used for error convergence tests.
53 C In the Newton iteration: S = RATE/(1 - RATE),
54 C where RATE is the estimated rate of convergence
55 C of the Newton iteration.
56 C The calling routine passes the initial value
57 C of S to the Newton iteration.
58 C CONFAC -- A residual scale factor to improve convergence.
59 C TOLNEW -- Tolerance on the norm of Newton correction in
60 C alternative Newton convergence test.
61 C MULDEL -- A flag indicating whether or not to multiply
62 C DELTA by CONFAC.
63 C 0 ==> do not scale DELTA by CONFAC.
64 C 1 ==> scale DELTA by CONFAC.
65 C MAXIT -- Maximum allowed number of Newton iterations.
66 C IRES -- Error flag returned from RES. See RES description
67 C in DDASPK prologue. If IRES = -1, then IERNEW
68 C will be set to 1.
69 C If IRES < -1, then IERNEW will be set to -1.
70 C IDUM -- Dummy argument.
71 C IERNEW -- Error flag for Newton iteration.
72 C 0 ==> Newton iteration converged.
73 C 1 ==> recoverable error inside Newton iteration.
74 C -1 ==> unrecoverable error inside Newton iteration.
75 C
76 C All arguments with "DUM" in their names are dummy arguments
77 C which are not used in this routine.
78 C-----------------------------------------------------------------------
79 C
80 C***ROUTINES CALLED
81 C DSLVD, DDWNRM, RES
82 C
83 C***END PROLOGUE DNSD
84 C
85 C
86  IMPLICIT DOUBLE PRECISION(a-h,o-z)
87  dimension y(*),yprime(*),wt(*),delta(*),e(*)
88  dimension wm(*),iwm(*), rpar(*),ipar(*)
89  EXTERNAL res
90 C
91  parameter(lnre=12, lnni=19)
92 C
93 C Initialize Newton counter M and accumulation vector E.
94 C
95  m = 0
96  DO 100 i=1,neq
97 100 e(i)=0.0d0
98 C
99 C Corrector loop.
100 C
101 300 CONTINUE
102  iwm(lnni) = iwm(lnni) + 1
103 C
104 C If necessary, multiply residual by convergence factor.
105 C
106  IF (muldel .EQ. 1) THEN
107  DO 320 i = 1,neq
108 320 delta(i) = delta(i) * confac
109  ENDIF
110 C
111 C Compute a new iterate (back-substitution).
112 C Store the correction in DELTA.
113 C
114  CALL dslvd(neq,delta,wm,iwm)
115 C
116 C Update Y, E, and YPRIME.
117 C
118  DO 340 i=1,neq
119  y(i)=y(i)-delta(i)
120  e(i)=e(i)-delta(i)
121 340 yprime(i)=yprime(i)-cj*delta(i)
122 C
123 C Test for convergence of the iteration.
124 C
125  delnrm=ddwnrm(neq,delta,wt,rpar,ipar)
126  IF (delnrm .LE. tolnew) GO TO 370
127  IF (m .EQ. 0) THEN
128  oldnrm = delnrm
129  ELSE
130  rate = (delnrm/oldnrm)**(1.0d0/m)
131  IF (rate .GT. 0.9d0) GO TO 380
132  s = rate/(1.0d0 - rate)
133  ENDIF
134  IF (s*delnrm .LE. epcon) GO TO 370
135 C
136 C The corrector has not yet converged.
137 C Update M and test whether the
138 C maximum number of iterations have
139 C been tried.
140 C
141  m=m+1
142  IF(m.GE.maxit) GO TO 380
143 C
144 C Evaluate the residual,
145 C and go back to do another iteration.
146 C
147  iwm(lnre)=iwm(lnre)+1
148  CALL res(x,y,yprime,cj,delta,ires,rpar,ipar)
149  IF (ires .LT. 0) GO TO 380
150  GO TO 300
151 C
152 C The iteration has converged.
153 C
154 370 RETURN
155 C
156 C The iteration has not converged. Set IERNEW appropriately.
157 C
158 380 CONTINUE
159  IF (ires .LE. -2 ) THEN
160  iernew = -1
161  ELSE
162  iernew = 1
163  ENDIF
164  RETURN
165 C
166 C
167 C------END OF SUBROUTINE DNSD-------------------------------------------
168  END
ddwnrm
double precision function ddwnrm(NEQ, V, RWT, RPAR, IPAR)
Definition: ddwnrm.f:6
dnsd
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
dslvd
subroutine dslvd(NEQ, DELTA, WM, IWM)
Definition: dslvd.f:6