-
Notifications
You must be signed in to change notification settings - Fork 2
/
nleq1e.f
231 lines (231 loc) · 9.24 KB
/
nleq1e.f
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
SUBROUTINE NLEQ1E(N,X,RTOL,IERR)
C* Begin Prologue NLEQ1E
INTEGER N
DOUBLE PRECISION X(N), RTOL
INTEGER IERR
C ------------------------------------------------------------
C
C* Title
C
C Numerical solution of nonlinear (NL) equations (EQ)
C especially designed for numerically sensitive problems.
C (E)asy-to-use driver routine for NLEQ1.
C
C* Written by U. Nowak, L. Weimann
C* Purpose Solution of systems of highly nonlinear equations
C* Method Damped affine invariant Newton method
C (see references below)
C* Category F2a. - Systems of nonlinear equations
C* Keywords Nonlinear equations, Newton methods
C* Version 2.3
C* Revision November 1991
C* Latest Change November 1991
C* Library CodeLib
C* Code Fortran 77, Double Precision
C* Environment Standard Fortran 77 environment on PC's,
C workstations and hosts.
C* Copyright (c) Konrad-Zuse-Zentrum fuer
C Informationstechnik Berlin (ZIB)
C Takustrasse 7, D-14195 Berlin-Dahlem
C phone : + 49/30/84185-0
C fax : + 49/30/84185-125
C* Contact Lutz Weimann
C ZIB, Division Scientific Computing,
C Department Numerical Analysis and Modelling
C phone : + 49/30/84185-185
C fax : + 49/30/84185-107
C e-mail: [email protected]
C
C* References:
C
C /1/ P. Deuflhard:
C Newton Methods for Nonlinear Problems. -
C Affine Invariance and Adaptive Algorithms.
C Series Computational Mathematics 35, Springer (2004)
C
C /2/ U. Nowak, L. Weimann:
C A Family of Newton Codes for Systems of Highly Nonlinear
C Equations - Algorithm, Implementation, Application.
C ZIB, Technical Report TR 90-10 (December 1990)
C
C ---------------------------------------------------------------
C
C* Licence
C You may use or modify this code for your own non commercial
C purposes for an unlimited time.
C In any case you should not deliver this code without a special
C permission of ZIB.
C In case you intend to use the code commercially, we oblige you
C to sign an according licence agreement with ZIB.
C
C* Warranty
C This code has been tested up to a certain level. Defects and
C weaknesses, which may be included in the code, do not establish
C any warranties by ZIB. ZIB does not take over any liabilities
C which may follow from aquisition or application of this code.
C
C* Software status
C This code is under care of ZIB and belongs to ZIB software class 1.
C
C ------------------------------------------------------------
C
C* Summary:
C ========
C Damped Newton-algorithm for systems of highly nonlinear
C equations - damping strategy due to Ref. (1).
C
C (The iteration is done by subroutine N1INT actually. NLEQ1E
C calls the standard interface driver NLEQ1, which itself does
C some house keeping and builds up workspace.)
C
C Jacobian approximation by numerical differences.
C
C The numerical solution of the arising linear equations is
C done by means of the subroutines *GEFA and *GESL ( Gauss-
C algorithm with column-pivoting and row-interchange;
C replace '*' by 'S' or 'D' for single or double precision
C version respectively).
C
C ------------------------------------------------------------
C
C* External subroutine to be supplied by the user
C ==============================================
C
C FCN(N,X,F,IFAIL) Ext Function subroutine - the problem function
C (The routine must be named exactly FCN !)
C N Int Number of vector components (input)
C Must not be altered!
C X(N) Dble Vector of unknowns (input)
C Must not be altered!
C F(N) Dble Vector of function values (output)
C IFAIL Int FCN evaluation-failure indicator. (output)
C On input: Has always value 0 (zero).
C On output: Indicates failure of FCN eval-
C uation, if having a nonzero value.
C If <0: NLEQ1E will be terminated with
C IFAIL returned via IERR.
C If =1: A new trial Newton iterate will be
C computed, with the damping factor
C reduced to it's half.
C If =2: A new trial Newton iterate will
C computed, with the damping factor
C reduced by a reduction factor,
C which must be output through F(1)
C by FCN, and it's value must be
C >0 and < 1.
C
C* Parameters list description
C ===========================
C
C Name Type In/Out Meaning
C
C N Int In Number of unknowns ( N .LE. 50 )
C X(N) Dble In Initial estimate of the solution
C Out Solution values ( or final values,
C respectively )
C RTOL Dble In Required relative precision of
C solution components -
C RTOL.GE.EPMACH*TEN*N
C Out Finally achieved accuracy
C IERR Int Out Return value parameter
C < 0 Termination forced by user function FCN
C due to IFAIL<0 on output, IERR is set
C to IFAIL
C = 0 successfull completion of the iteration,
C solution has been computed
C > 0 see list of error/warning messages below
C
C* Error and warning messages:
C ===========================
C
C 1 Termination, since Jacobian matrix became singular
C 2 Termination after 100 iterations
C 3 Termination, since damping factor became to small
C 4 Warning: Superlinear or quadratic convergence slowed down
C near the solution.
C Iteration has been stopped therefore with an approximation
C of the solution not such accurate as requested by RTOL,
C because possibly the RTOL requirement may be too stringent
C (i.e. the nonlinear problem is ill-conditioned)
C 5 Warning: Iteration stopped with termination criterion
C (using RTOL as requested precision) satisfied, but no
C superlinear or quadratic convergence has been indicated yet.
C Therefore, possibly the error estimate for the solution may
C not match good enough the really achieved accuracy.
C 20 Bad input value to parameter N; 1 .LE. N .LE. 50 required
C 21 Nonpositive value for RTOL supplied
C 82 Termination, because user routine FCN returned with IFAIL>0
C
C Note : in case of failure:
C - use better initial guess
C - or refine model
C - or use non-standard options and/or analytical Jacobian
C via the standard interface NLEQ1
C
C* Machine dependent constants used:
C =================================
C
C DOUBLE PRECISION EPMACH in N1PCHK, N1INT
C DOUBLE PRECISION GREAT in N1PCHK
C DOUBLE PRECISION SMALL in N1PCHK, N1INT, N1SCAL
C
C* Subroutines called: NLEQ1
C
C ------------------------------------------------------------
C* End Prologue
INTEGER NMAX, LIOPT, LIWK, LRWK, LUPRT
PARAMETER ( NMAX=50 )
PARAMETER ( LIOPT=50, LIWK=NMAX+50, LRWK=(NMAX+13)*NMAX+60 )
PARAMETER ( LUPRT=6 )
EXTERNAL FCN, NLEQ1
DOUBLE PRECISION XSCAL(NMAX)
INTEGER IOPT(LIOPT), IWK(LIWK)
DOUBLE PRECISION RWK(LRWK)
INTEGER I, NIW, NRW
CHARACTER CHGDAT*20, PRODCT*8
C
C Version: 2.3 Latest change:
C -----------------------------------------
C
DATA CHGDAT /'November 15, 1991 '/
DATA PRODCT /'NLEQ1E '/
C* Begin
NIW = N+50
NRW = (N+13)*N+60
C Checking dimensional parameter N
IF ( N.LT.1 .OR. N.GT.NMAX ) THEN
WRITE(LUPRT,1001) NMAX, N
1001 FORMAT(/,' Error: Bad input to parameter N supplied',
$ /,8X,'choose 1 .LE. N .LE. ',I3,
$ ' , your input is: N = ',I5)
IERR = 20
RETURN
ENDIF
DO 10 I=1,LIOPT
IOPT(I) = 0
10 CONTINUE
DO 20 I=1,NIW
IWK(I) = 0
20 CONTINUE
DO 30 I=1,NRW
RWK(I) = 0.0D0
30 CONTINUE
DO 40 I=1,N
XSCAL(I) = 0.0D0
40 CONTINUE
C Print errors, warnings, monitor and time monitor
C to standard output
IOPT(11) = 3
IOPT(12) = LUPRT
IOPT(13) = 3
IOPT(14) = LUPRT
IOPT(19) = 1
IOPT(20) = LUPRT
C Maximum number of Newton iterations
IWK(31) = 100
C
CALL NLEQ1(N,FCN,DUMMY,X,XSCAL,RTOL,IOPT,IERR,LIWK,IWK,LRWK,RWK)
IF (IERR.EQ.82 .AND. IWK(23).LT.0) IERR=IWK(23)
C End of subroutine NLEQ1E
RETURN
END