-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathnleq2e.c
287 lines (242 loc) · 11.1 KB
/
nleq2e.c
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
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
/* nleq2e.f -- translated by f2c (version 20100827).
You must link the resulting object file with libf2c:
on Microsoft Windows system, link with libf2c.lib;
on Linux or Unix systems, link with .../path/to/libf2c.a -lm
or, if you install libf2c.a in a standard place, with -lf2c -lm
-- in that order, at the end of the command line, as in
cc *.o -lf2c -lm
Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
http://www.netlib.org/f2c/libf2c.zip
*/
#include "f2c_nleq.h"
/* Table of constant values */
static integer c__50 = 50;
static integer c__1 = 1;
static integer c__100 = 100;
static integer c__3210 = 3210;
/* Subroutine */ int nleq2e_(integer *n, doublereal *x, doublereal *rtol,
integer *ierr)
{
/* Format strings */
static char fmt_1001[] = "(/,\002 Error: Bad input to parameter N suppli"
"ed\002,/,8x,\002choose 1 .LE. N .LE. \002,i3,\002 , your input i"
"s: N = \002,i5)";
/* System generated locals */
integer i__1;
/* Builtin functions */
integer s_wsfe(cilist *), do_fio(integer *, char *, ftnlen), e_wsfe(void);
/* Local variables */
static integer i__;
extern /* Subroutine */ int fcn_();
static integer iwk[100], niw;
static doublereal rwk[3210];
static integer nrw, iopt[50];
extern /* Subroutine */ int nleq2_(integer *, U_fp, f2c_real *, doublereal *,
doublereal *, doublereal *, integer *, integer *, integer *,
integer *, integer *, doublereal *);
static doublereal xscal[50];
static integer luerr;
static f2c_real dummy;
static integer mprerr;
/* Fortran I/O blocks */
static cilist io___7 = { 0, 0, 0, fmt_1001, 0 };
/* * Begin Prologue NLEQ2E */
/* ------------------------------------------------------------ */
/* * Title */
/* Numerical solution of nonlinear (NL) equations (EQ) */
/* especially designed for numerically sensitive problems. */
/* (E)asy-to-use driver routine for NLEQ2. */
/* * Written by U. Nowak, L. Weimann */
/* * Purpose Solution of systems of highly nonlinear equations */
/* * Method Damped affine invariant Newton method */
/* (see references below) */
/* * Category F2a. - Systems of nonlinear equations */
/* * Keywords Nonlinear equations, Newton methods */
/* * Version 2.3 */
/* * Revision November 1991 */
/* * Latest Change November 1991 */
/* * Library CodeLib */
/* * Code Fortran 77, Double Precision */
/* * Environment Standard Fortran 77 environment on PC's, */
/* workstations and hosts. */
/* * Copyright (c) Konrad-Zuse-Zentrum fuer */
/* Informationstechnik Berlin (ZIB) */
/* Takustrasse 7, D-14195 Berlin-Dahlem */
/* phone : + 49/30/84185-0 */
/* fax : + 49/30/84185-125 */
/* * Contact Bodo Erdmann */
/* ZIB, Division Scientific Computing, */
/* Department Numerical Analysis and Modelling */
/* phone : + 49/30/84185-185 */
/* fax : + 49/30/84185-107 */
/* e-mail: [email protected] */
/* * References: */
/* /1/ P. Deuflhard: */
/* Newton Methods for Nonlinear Problems. - */
/* Affine Invariance and Adaptive Algorithms. */
/* Series Computational Mathematics 35, Springer (2004) */
/* /2/ U. Nowak, L. Weimann: */
/* A Family of Newton Codes for Systems of Highly Nonlinear */
/* Equations - Algorithm, Implementation, Application. */
/* ZIB, Technical Report TR 90-10 (December 1990) */
/* --------------------------------------------------------------- */
/* * Licence */
/* You may use or modify this code for your own non commercial */
/* purposes for an unlimited time. */
/* In any case you should not deliver this code without a special */
/* permission of ZIB. */
/* In case you intend to use the code commercially, we oblige you */
/* to sign an according licence agreement with ZIB. */
/* * Warranty */
/* This code has been tested up to a certain level. Defects and */
/* weaknesses, which may be included in the code, do not establish */
/* any warranties by ZIB. ZIB does not take over any liabilities */
/* which may follow from aquisition or application of this code. */
/* * Software status */
/* This code is under care of ZIB and belongs to ZIB software class 1. */
/* ------------------------------------------------------------ */
/* * Summary: */
/* ======== */
/* Damped Newton-algorithm for systems of highly nonlinear */
/* equations - damping strategy due to Ref. (1). */
/* (The iteration is done by subroutine N2INT actually. NLEQ2E */
/* calls the standard interface driver NLEQ2, which itself does */
/* some house keeping and builds up workspace.) */
/* Jacobian approximation by numerical differences. */
/* The numerical solution of the arising linear equations is */
/* done by means of the subroutines *GEFA and *GESL ( Gauss- */
/* algorithm with column-pivoting and row-interchange; */
/* replace '*' by 'S' or 'D' for single or double precision */
/* version respectively). */
/* ------------------------------------------------------------ */
/* * External subroutine to be supplied by the user */
/* ============================================== */
/* FCN(N,X,F,IFAIL) Ext Function subroutine - the problem function */
/* (The routine must be named exactly FCN !) */
/* N Int Number of vector components (input) */
/* Must not be altered! */
/* X(N) Dble Vector of unknowns (input) */
/* Must not be altered! */
/* F(N) Dble Vector of function values (output) */
/* IFAIL Int FCN evaluation-failure indicator. (output) */
/* On input: Has always value 0 (zero). */
/* On output: Indicates failure of FCN eval- */
/* uation, if having a nonzero value. */
/* If <0: NLEQ2E will be terminated with */
/* IFAIL returned via IERR. */
/* If =1: A new trial Newton iterate will be */
/* computed, with the damping factor */
/* reduced to it's half. */
/* If =2: A new trial Newton iterate will */
/* computed, with the damping factor */
/* reduced by a reduction factor, */
/* which must be output through F(1) */
/* by FCN, and it's value must be */
/* >0 and < 1. */
/* * Parameters list description */
/* =========================== */
/* Name Type In/Out Meaning */
/* N Int In Number of unknowns ( N .LE. 50 ) */
/* X(N) Dble In Initial estimate of the solution */
/* Out Solution values ( or final values, */
/* respectively ) */
/* RTOL Dble In Required relative precision of */
/* solution components - */
/* RTOL.GE.EPMACH*TEN*N */
/* Out Finally achieved accuracy */
/* IERR Int Out Return value parameter */
/* < 0 Termination forced by user function FCN */
/* due to IFAIL<0 on output, IERR is set */
/* to IFAIL */
/* = 0 successfull completion of the iteration, */
/* solution has been computed */
/* > 0 see list of error/warning messages below */
/* * Error and warning messages: */
/* =========================== */
/* 1 Termination, since Jacobian matrix became singular */
/* 2 Termination after 100 iterations */
/* 3 Termination, since damping factor became to small */
/* 4 Warning: Superlinear or quadratic convergence slowed down */
/* near the solution. */
/* Iteration has been stopped therefore with an approximation */
/* of the solution not such accurate as requested by RTOL, */
/* because possibly the RTOL requirement may be too stringent */
/* (i.e. the nonlinear problem is ill-conditioned) */
/* 5 Warning: Iteration stopped with termination criterion */
/* (using RTOL as requested precision) satisfied, but no */
/* superlinear or quadratic convergence has been indicated yet. */
/* Therefore, possibly the error estimate for the solution may */
/* not match good enough the really achieved accuracy. */
/* 20 Bad input value to parameter N; 1 .LE. N .LE. 50 required */
/* 21 Nonpositive value for RTOL supplied */
/* 82 Termination, because user routine FCN returned with IFAIL>0 */
/* Note : in case of failure: */
/* - use better initial guess */
/* - or refine model */
/* - or use non-standard options and/or analytical Jacobian */
/* via the standard interface NLEQ2 */
/* * Machine dependent constants used: */
/* ================================= */
/* DOUBLE PRECISION EPMACH in N2PCHK, N2INT */
/* DOUBLE PRECISION GREAT in N2PCHK */
/* DOUBLE PRECISION SMALL in N2PCHK, N2INT, N2SCAL */
/* * Subroutines called: NLEQ2 */
/* ------------------------------------------------------------ */
/* * End Prologue */
/* Version: 2.3 Latest change: */
/* ----------------------------------------- */
/* Parameter adjustments */
--x;
/* Function Body */
/* * Begin */
niw = *n + 50;
nrw = (*n + 13) * *n + 60;
/* Checking dimensional parameter N */
if (*n < 1 || *n > 50) {
if (mprerr >= 1) {
io___7.ciunit = luerr;
s_wsfe(&io___7);
do_fio(&c__1, (char *)&c__50, (ftnlen)sizeof(integer));
do_fio(&c__1, (char *)&(*n), (ftnlen)sizeof(integer));
e_wsfe();
}
*ierr = 20;
return 0;
}
for (i__ = 1; i__ <= 50; ++i__) {
iopt[i__ - 1] = 0;
/* L10: */
}
i__1 = niw;
for (i__ = 1; i__ <= i__1; ++i__) {
iwk[i__ - 1] = 0;
/* L20: */
}
i__1 = nrw;
for (i__ = 1; i__ <= i__1; ++i__) {
rwk[i__ - 1] = 0.;
/* L30: */
}
i__1 = *n;
for (i__ = 1; i__ <= i__1; ++i__) {
xscal[i__ - 1] = 0.;
/* L40: */
}
/* Print errors, warnings, monitor and time monitor */
/* to standard output */
iopt[10] = 3;
iopt[11] = 6;
iopt[12] = 3;
iopt[13] = 6;
iopt[18] = 1;
iopt[19] = 6;
/* Maximum number of Newton iterations */
iwk[30] = 100;
nleq2_(n, (U_fp)fcn_, &dummy, &x[1], xscal, rtol, iopt, ierr, &c__100,
iwk, &c__3210, rwk);
if (*ierr == 82 && iwk[22] < 0) {
*ierr = iwk[22];
}
/* End of subroutine NLEQ2E */
return 0;
} /* nleq2e_ */