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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 89 - (hide annotations) (download) (as text)
Wed Dec 7 15:44:43 2005 UTC (18 years, 7 months ago) by johnpye
File MIME type: text/x-csrc
File size: 5092 byte(s)
Small changes to eliminate GCC warnings
1 aw0a 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 johnpye 89 * See page 360 of nr in C. We will later convert the netlib code
35 aw0a 1 * 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 johnpye 89
41     #ifdef STAND_ALONE
42     # include "codegen_support.h"
43 aw0a 1 #else
44 johnpye 89 # error "DON'T COMPILE THIS FILE UNLESS YOU ARE IN STANDALONE MODE"
45     # include "compiler/extfunc.h" /* ExtEvalFunc */
46     #endif
47    
48 aw0a 1 #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 johnpye 62 error_reporter(ASC_USER_ERROR,NULL,0,"Solver: zbrent: Root must be bracketed.");
104 aw0a 1 *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 johnpye 62 error_reporter(ASC_USER_ERROR,NULL,0,"Solver: zbrent: Maximum number of iterations exceeded.");
174 aw0a 1 *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