GNU Octave  3.8.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
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