ascend4/solver/rootfind.c

/*
 *  SLV: Ascend Nonlinear Solver
 *  by Kirk Andre' Abbott
 *  Created: 10/06/95
 *  Version: $Revision: 1.4 $
 *  Version control file: $RCSfile: rootfind.c,v $
 *  Date last modified: $Date: 1997/07/18 12:15:37 $
 *  Last modified by: $Author: mthomas $
 *
 *  This file is part of the SLV solver.
 *
 *  Copyright (C) 1990 Karl Michael Westerberg
 *  Copyright (C) 1993 Joseph Zaher
 *  Copyright (C) 1994 Joseph Zaher, Benjamin Andrew Allan
 *
 *  The SLV solver is free software; you can redistribute
 *  it and/or modify it under the terms of the GNU General Public License as
 *  published by the Free Software Foundation; either version 2 of the
 *  License, or (at your option) any later version.
 *
 *  The SLV solver is distributed in hope that it will be
 *  useful, but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 *  General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public License along with
 *  the program; if not, write to the Free Software Foundation, Inc., 675
 *  Mass Ave, Cambridge, MA 02139 USA.  Check the file named COPYING.
 *  COPYING is found in ../compiler.
 */

/*
 * This is the first pass implemenation of some rootfinding codes.
 * See page 360 or nr in C. We will later convert the netlib code
 * and use this to avoid copyright problems with the nr in C boys.
 */
#include <math.h>
#include "utilities/ascConfig.h"
#include "compiler/compiler.h"
#ifndef STAND_ALONE
#include "compiler/extfunc.h"   /* ExtEvalFunc */
#else
#include "codegen_support.h"
#endif /* STAND_ALONE */
#include "compiler/rootfind.h"

#define ITMAX 100
#define EPS 1.0e-08
#define ZBIGNUM 1.0e08

/*
 * We have maintained a consistent calling protocol between
 * the (possibly) different versions of the rootfinding code.
 * In this version, m is not used; n is the index into to
 * the x-vector, of the variable that we are solving for.
 * Our residuals are always written to f[0].
 */
double zbrent(ExtEvalFunc *func,    /* the evaluation function */
          double *lowbound,     /* low bound */
          double *upbound,      /* up bound */
          int *mode,        /* to pass to the eval func */
          int *m,           /* the relation index */
          int *n,           /* the variable index */
          double *x,    /* the x vector -- needed by eval func */
          double *u,    /* the u vector -- needed by eval func */
          double *f,        /* vector of residuals */
          double *g,        /* vector of gradients */
          double *tolerance,    /* accuracy of solution */
          int *status)      /* success or failure */

{
  int iter, result;
  double x1, x2;
  double a, b, c;
  double fa, fb;
  double d, e, min1, min2;
  double fc, p, q, r, s, tol1, xm;

  x1 = *lowbound;       /* initialization */
  x2 = *upbound;
  a = x1, b = x2, c = x2;

  x[*n] = a;
  result = (*func)(mode,m,n,x,u,f,g);
  fa = f[*m];
  if (result) {
    *status = -1;
    return ZBIGNUM;
  }

  x[*n] = b;
  result = (*func)(mode,m,n,x,u,f,g);
  fb = f[*m];
  if (result) {
    *status = -1;
    return ZBIGNUM;
  }

  if ((fa > 0.0 && fb > 0.0) || (fa < 0.0 && fb < 0.0)) {
    FPRINTF(stderr,"root must be bracketed in zbrent\n");
    *status = -1;       /* cannot invert */
    return ZBIGNUM;
  }
  fc = fb;
  for (iter=1;iter<=ITMAX;iter++) {
    if ((fb > 0.0 && fc > 0.0) || (fb < 0.0 && fc < 0.0)) {
      c = a;        /* rename a, b, c and adjust bounding interval d. */
      fc = fa;
      e = d = b - a;
    }
    if (fabs(fc) < fabs(fb)) {
      a = b;
      b = c;
      c = a;
      fa = fb;
      fb = fc;
      fc = fa;
    }
    tol1 = 2.0*EPS*fabs(b) + 0.5 * (*tolerance); /* Convergence check */
    xm = 0.5*(c-b);
    if (fabs(xm) <= tol1 || fb == 0.0) {
      x[*n] = b;
      *status = 0;
      return b;
    }
    if (fabs(e) >= tol1 && fabs(fa) > fabs(fb)) {
      s = fb/fa;    /* Attempt inverse quadratic interpolation. */
      if (a == c) {
    p = 2.0*xm*s;
    q = 1.0 - s;
      }
      else{
    q = fa/fc;
    r = fb/fc;
    p = s*(2.0*xm*q*(q-r) - (b-a)*(r-1.0));
    q = (q-1.0)*(r-1.0)*(s-1.0);
      }
      if (p > 0.0) q = -q;  /* check whether in bounds */
      p = fabs(p);
      min1 = 3.0*xm*q - fabs(tol1*q);
      min2 = fabs(e*q);
      if (2.0*p < (min1 < min2 ? min1 : min2)) {
    e = d;          /* accept interpolation */
    d = p/q;
      }
      else{
    d = xm;
    e = d;      /* interpolation failed; use bisection */
      }
    }
    else{       /* bounds decreasing too slowly; use bisection */
      d = xm;
      e = d;
    }
    a = b;          /* move last best quess to a. */
    fa = fb;
    if (fabs(d) > tol1)     /* evaluate new trial root */
      b += d;
    else
      b += (xm > 0.0 ? fabs(tol1) : -fabs(tol1));

    x[*n] = b;
    result = (*func)(mode,m,n,x,u,f,g);
    fb = f[*m]  ;
    if (result) {
      *status = -1;
      return ZBIGNUM;
    }
  }
  FPRINTF(stderr,"Maximum number of iterations exceeded in zbrent\n");
  *status = -1;     /* cannot invert */
  return ZBIGNUM;   /* NOTREACHED */
}


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.0EPSfabs(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.0xms;
130	q = 1.0 - s;
131	}
132	else{
133	q = fa/fc;
134	r = fb/fc;
135	p = s(2.0xmq(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.0xmq - 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