/[ascend]/trunk/base/generic/solver/rootfind.c
ViewVC logotype

Contents of /trunk/base/generic/solver/rootfind.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 89 - (show annotations) (download) (as text)
Wed Dec 7 15:44:43 2005 UTC (14 years ago) by johnpye
File MIME type: text/x-csrc
File size: 5092 byte(s)
Small changes to eliminate GCC warnings
1 /*
2 * SLV: Ascend Nonlinear Solver
3 * by Kirk Andre' Abbott
4 * Created: 10/06/95
5 * Version: $Revision: 1.4 $
6 * Version control file: $RCSfile: rootfind.c,v $
7 * Date last modified: $Date: 1997/07/18 12:15:37 $
8 * Last modified by: $Author: mthomas $
9 *
10 * This file is part of the SLV solver.
11 *
12 * Copyright (C) 1990 Karl Michael Westerberg
13 * Copyright (C) 1993 Joseph Zaher
14 * Copyright (C) 1994 Joseph Zaher, Benjamin Andrew Allan
15 *
16 * The SLV solver is free software; you can redistribute
17 * it and/or modify it under the terms of the GNU General Public License as
18 * published by the Free Software Foundation; either version 2 of the
19 * License, or (at your option) any later version.
20 *
21 * The SLV solver is distributed in hope that it will be
22 * useful, but WITHOUT ANY WARRANTY; without even the implied warranty of
23 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
24 * General Public License for more details.
25 *
26 * You should have received a copy of the GNU General Public License along with
27 * the program; if not, write to the Free Software Foundation, Inc., 675
28 * Mass Ave, Cambridge, MA 02139 USA. Check the file named COPYING.
29 * COPYING is found in ../compiler.
30 */
31
32 /*
33 * This is the first pass implemenation of some rootfinding codes.
34 * See page 360 of nr in C. We will later convert the netlib code
35 * and use this to avoid copyright problems with the nr in C boys.
36 */
37 #include <math.h>
38 #include "utilities/ascConfig.h"
39 #include "compiler/compiler.h"
40
41 #ifdef STAND_ALONE
42 # include "codegen_support.h"
43 #else
44 # error "DON'T COMPILE THIS FILE UNLESS YOU ARE IN STANDALONE MODE"
45 # include "compiler/extfunc.h" /* ExtEvalFunc */
46 #endif
47
48 #include "compiler/rootfind.h"
49
50 #define ITMAX 100
51 #define EPS 1.0e-08
52 #define ZBIGNUM 1.0e08
53
54 /*
55 * We have maintained a consistent calling protocol between
56 * the (possibly) different versions of the rootfinding code.
57 * In this version, m is not used; n is the index into to
58 * the x-vector, of the variable that we are solving for.
59 * Our residuals are always written to f[0].
60 */
61 double zbrent(ExtEvalFunc *func, /* the evaluation function */
62 double *lowbound, /* low bound */
63 double *upbound, /* up bound */
64 int *mode, /* to pass to the eval func */
65 int *m, /* the relation index */
66 int *n, /* the variable index */
67 double *x, /* the x vector -- needed by eval func */
68 double *u, /* the u vector -- needed by eval func */
69 double *f, /* vector of residuals */
70 double *g, /* vector of gradients */
71 double *tolerance, /* accuracy of solution */
72 int *status) /* success or failure */
73
74 {
75 int iter, result;
76 double x1, x2;
77 double a, b, c;
78 double fa, fb;
79 double d, e, min1, min2;
80 double fc, p, q, r, s, tol1, xm;
81
82 x1 = *lowbound; /* initialization */
83 x2 = *upbound;
84 a = x1, b = x2, c = x2;
85
86 x[*n] = a;
87 result = (*func)(mode,m,n,x,u,f,g);
88 fa = f[*m];
89 if (result) {
90 *status = -1;
91 return ZBIGNUM;
92 }
93
94 x[*n] = b;
95 result = (*func)(mode,m,n,x,u,f,g);
96 fb = f[*m];
97 if (result) {
98 *status = -1;
99 return ZBIGNUM;
100 }
101
102 if ((fa > 0.0 && fb > 0.0) || (fa < 0.0 && fb < 0.0)) {
103 error_reporter(ASC_USER_ERROR,NULL,0,"Solver: zbrent: Root must be bracketed.");
104 *status = -1; /* cannot invert */
105 return ZBIGNUM;
106 }
107 fc = fb;
108 for (iter=1;iter<=ITMAX;iter++) {
109 if ((fb > 0.0 && fc > 0.0) || (fb < 0.0 && fc < 0.0)) {
110 c = a; /* rename a, b, c and adjust bounding interval d. */
111 fc = fa;
112 e = d = b - a;
113 }
114 if (fabs(fc) < fabs(fb)) {
115 a = b;
116 b = c;
117 c = a;
118 fa = fb;
119 fb = fc;
120 fc = fa;
121 }
122 tol1 = 2.0*EPS*fabs(b) + 0.5 * (*tolerance); /* Convergence check */
123 xm = 0.5*(c-b);
124 if (fabs(xm) <= tol1 || fb == 0.0) {
125 x[*n] = b;
126 *status = 0;
127 return b;
128 }
129 if (fabs(e) >= tol1 && fabs(fa) > fabs(fb)) {
130 s = fb/fa; /* Attempt inverse quadratic interpolation. */
131 if (a == c) {
132 p = 2.0*xm*s;
133 q = 1.0 - s;
134 }
135 else{
136 q = fa/fc;
137 r = fb/fc;
138 p = s*(2.0*xm*q*(q-r) - (b-a)*(r-1.0));
139 q = (q-1.0)*(r-1.0)*(s-1.0);
140 }
141 if (p > 0.0) q = -q; /* check whether in bounds */
142 p = fabs(p);
143 min1 = 3.0*xm*q - fabs(tol1*q);
144 min2 = fabs(e*q);
145 if (2.0*p < (min1 < min2 ? min1 : min2)) {
146 e = d; /* accept interpolation */
147 d = p/q;
148 }
149 else{
150 d = xm;
151 e = d; /* interpolation failed; use bisection */
152 }
153 }
154 else{ /* bounds decreasing too slowly; use bisection */
155 d = xm;
156 e = d;
157 }
158 a = b; /* move last best quess to a. */
159 fa = fb;
160 if (fabs(d) > tol1) /* evaluate new trial root */
161 b += d;
162 else
163 b += (xm > 0.0 ? fabs(tol1) : -fabs(tol1));
164
165 x[*n] = b;
166 result = (*func)(mode,m,n,x,u,f,g);
167 fb = f[*m] ;
168 if (result) {
169 *status = -1;
170 return ZBIGNUM;
171 }
172 }
173 error_reporter(ASC_USER_ERROR,NULL,0,"Solver: zbrent: Maximum number of iterations exceeded.");
174 *status = -1; /* cannot invert */
175 return ZBIGNUM; /* NOTREACHED */
176 }
177
178
179

john.pye@anu.edu.au
ViewVC Help
Powered by ViewVC 1.1.22