/[ascend]/trunk/ascend4/solver/rootfind.c
ViewVC logotype

Contents of /trunk/ascend4/solver/rootfind.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1 - (show annotations) (download) (as text)
Fri Oct 29 20:54:12 2004 UTC (17 years, 6 months ago) by aw0a
File MIME type: text/x-csrc
File size: 4985 byte(s)
Setting up web subdirectory in repository
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 or 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 #ifndef STAND_ALONE
41 #include "compiler/extfunc.h" /* ExtEvalFunc */
42 #else
43 #include "codegen_support.h"
44 #endif /* STAND_ALONE */
45 #include "compiler/rootfind.h"
46
47 #define ITMAX 100
48 #define EPS 1.0e-08
49 #define ZBIGNUM 1.0e08
50
51 /*
52 * We have maintained a consistent calling protocol between
53 * the (possibly) different versions of the rootfinding code.
54 * In this version, m is not used; n is the index into to
55 * the x-vector, of the variable that we are solving for.
56 * Our residuals are always written to f[0].
57 */
58 double zbrent(ExtEvalFunc *func, /* the evaluation function */
59 double *lowbound, /* low bound */
60 double *upbound, /* up bound */
61 int *mode, /* to pass to the eval func */
62 int *m, /* the relation index */
63 int *n, /* the variable index */
64 double *x, /* the x vector -- needed by eval func */
65 double *u, /* the u vector -- needed by eval func */
66 double *f, /* vector of residuals */
67 double *g, /* vector of gradients */
68 double *tolerance, /* accuracy of solution */
69 int *status) /* success or failure */
70
71 {
72 int iter, result;
73 double x1, x2;
74 double a, b, c;
75 double fa, fb;
76 double d, e, min1, min2;
77 double fc, p, q, r, s, tol1, xm;
78
79 x1 = *lowbound; /* initialization */
80 x2 = *upbound;
81 a = x1, b = x2, c = x2;
82
83 x[*n] = a;
84 result = (*func)(mode,m,n,x,u,f,g);
85 fa = f[*m];
86 if (result) {
87 *status = -1;
88 return ZBIGNUM;
89 }
90
91 x[*n] = b;
92 result = (*func)(mode,m,n,x,u,f,g);
93 fb = f[*m];
94 if (result) {
95 *status = -1;
96 return ZBIGNUM;
97 }
98
99 if ((fa > 0.0 && fb > 0.0) || (fa < 0.0 && fb < 0.0)) {
100 FPRINTF(stderr,"root must be bracketed in zbrent\n");
101 *status = -1; /* cannot invert */
102 return ZBIGNUM;
103 }
104 fc = fb;
105 for (iter=1;iter<=ITMAX;iter++) {
106 if ((fb > 0.0 && fc > 0.0) || (fb < 0.0 && fc < 0.0)) {
107 c = a; /* rename a, b, c and adjust bounding interval d. */
108 fc = fa;
109 e = d = b - a;
110 }
111 if (fabs(fc) < fabs(fb)) {
112 a = b;
113 b = c;
114 c = a;
115 fa = fb;
116 fb = fc;
117 fc = fa;
118 }
119 tol1 = 2.0*EPS*fabs(b) + 0.5 * (*tolerance); /* Convergence check */
120 xm = 0.5*(c-b);
121 if (fabs(xm) <= tol1 || fb == 0.0) {
122 x[*n] = b;
123 *status = 0;
124 return b;
125 }
126 if (fabs(e) >= tol1 && fabs(fa) > fabs(fb)) {
127 s = fb/fa; /* Attempt inverse quadratic interpolation. */
128 if (a == c) {
129 p = 2.0*xm*s;
130 q = 1.0 - s;
131 }
132 else{
133 q = fa/fc;
134 r = fb/fc;
135 p = s*(2.0*xm*q*(q-r) - (b-a)*(r-1.0));
136 q = (q-1.0)*(r-1.0)*(s-1.0);
137 }
138 if (p > 0.0) q = -q; /* check whether in bounds */
139 p = fabs(p);
140 min1 = 3.0*xm*q - fabs(tol1*q);
141 min2 = fabs(e*q);
142 if (2.0*p < (min1 < min2 ? min1 : min2)) {
143 e = d; /* accept interpolation */
144 d = p/q;
145 }
146 else{
147 d = xm;
148 e = d; /* interpolation failed; use bisection */
149 }
150 }
151 else{ /* bounds decreasing too slowly; use bisection */
152 d = xm;
153 e = d;
154 }
155 a = b; /* move last best quess to a. */
156 fa = fb;
157 if (fabs(d) > tol1) /* evaluate new trial root */
158 b += d;
159 else
160 b += (xm > 0.0 ? fabs(tol1) : -fabs(tol1));
161
162 x[*n] = b;
163 result = (*func)(mode,m,n,x,u,f,g);
164 fb = f[*m] ;
165 if (result) {
166 *status = -1;
167 return ZBIGNUM;
168 }
169 }
170 FPRINTF(stderr,"Maximum number of iterations exceeded in zbrent\n");
171 *status = -1; /* cannot invert */
172 return ZBIGNUM; /* NOTREACHED */
173 }
174
175
176

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