From oleg@nrcbsa.bio.nrc.ca Thu May 23 08:50:45 1991
Return-Path:
Received: from nrcnet0.nrc.ca by CS.UTK.EDU with SMTP (5.61++/2.5.1s-UTK)
id AA25956; Thu, 23 May 91 08:50:36 -0400
Message-Id: <9105231250.AA28656@ nrcbsa.bio.nrc.ca>
Date: Thu, 23 May 91 08:50:09 EDT
From: Dr. Oleg Keselyov
To: dongarra@cs.utk.edu
Subject: brent.shar
Status: RO
#--------------------------------CUT HERE-------------------------------------
#! /bin/sh
#
# This is a shell archive. Save this into a file, edit it
# and delete all lines above this comment. Then give this
# file to sh by executing the command 'sh file'. The files
# will be extracted into the current directory owned by
# you with default permissions.
# This file contains Brent's univariate minimizer and zero finder.
# C realization (with adaptation) of the algorithm
# G.Forsythe, M.Malcolm, C.Moler, Computer methods for
# mathematical computations.
# Contains the source code for the program fminbr.c and
# zeroin.c, test drivers for both and verifivation protocols.
# As to header files and service functions, see serv.shar
# Contents of the archive
#
# Source code for the programs being submitted
# fminbr.c Brent's univariate minimizer.
# zeroin.c Brent's univariate zero-finder.
#
# Validation routines
# vfminbr.c Verify fminbr.c function
# vzeroin.c Verify zeroin.c function
#
# Verification protocols
# vfminbr.dat On Silicon Graphics IRIS
# vzeroin.dat
# vfminbr.dat.vax On VAX station
# vzeroin.dat.vax
echo 'x - fminbr.c'
sed 's/^X//' << '________This_Is_The_END________' > fminbr.c
X/*
X ************************************************************************
X * C math library
X * function FMINBR - one-dimensional search for a function minimum
X * over the given range
X *
X * Input
X * double fminbr(a,b,f,tol)
X * double a; Minimum will be seeked for over
X * double b; a range [a,b], a being < b.
X * double (*f)(double x); Name of the function whose minimum
X * will be seeked for
X * double tol; Acceptable tolerance for the minimum
X * location. It have to be positive
X * (e.g. may be specified as EPSILON)
X *
X * Output
X * Fminbr returns an estimate for the minimum location with accuracy
X * 3*SQRT_EPSILON*abs(x) + tol.
X * The function always obtains a local minimum which coincides with
X * the global one only if a function under investigation being
X * unimodular.
X * If a function being examined possesses no local minimum within
X * the given range, Fminbr returns 'a' (if f(a) < f(b)), otherwise
X * it returns the right range boundary value b.
X *
X * Algorithm
X * G.Forsythe, M.Malcolm, C.Moler, Computer methods for mathematical
X * computations. M., Mir, 1980, p.202 of the Russian edition
X *
X * The function makes use of the "gold section" procedure combined with
X * the parabolic interpolation.
X * At every step program operates three abscissae - x,v, and w.
X * x - the last and the best approximation to the minimum location,
X * i.e. f(x) <= f(a) or/and f(x) <= f(b)
X * (if the function f has a local minimum in (a,b), then the both
X * conditions are fulfiled after one or two steps).
X * v,w are previous approximations to the minimum location. They may
X * coincide with a, b, or x (although the algorithm tries to make all
X * u, v, and w distinct). Points x, v, and w are used to construct
X * interpolating parabola whose minimum will be treated as a new
X * approximation to the minimum location if the former falls within
X * [a,b] and reduces the range enveloping minimum more efficient than
X * the gold section procedure.
X * When f(x) has a second derivative positive at the minimum location
X * (not coinciding with a or b) the procedure converges superlinearly
X * at a rate order about 1.324
X *
X ************************************************************************
X */
X
X#include "assert.h"
X#include "math.h"
X
Xdouble fminbr(a,b,f,tol) /* An estimate to the min location*/
Xdouble a; /* Left border | of the range */
Xdouble b; /* Right border| the min is seeked*/
Xdouble (*f)(double x); /* Function under investigation */
Xdouble tol; /* Acceptable tolerance */
X{
X double x,v,w; /* Abscissae, descr. see above */
X double fx; /* f(x) */
X double fv; /* f(v) */
X double fw; /* f(w) */
X const double r = (3.-sqrt(5.0))/2; /* Gold section ratio */
X
X assert( tol > 0 && b > a );
X
X v = a + r*(b-a); fv = (*f)(v); /* First step - always gold section*/
X x = v; w = v;
X fx=fv; fw=fv;
X
X for(;;) /* Main iteration loop */
X {
X double range = b-a; /* Range over which the minimum */
X /* is seeked for */
X double middle_range = (a+b)/2;
X double tol_act = /* Actual tolerance */
X SQRT_EPSILON*fabs(x) + tol/3;
X double new_step; /* Step at this iteration */
X
X
X
X if( fabs(x-middle_range) + range/2 <= 2*tol_act )
X return x; /* Acceptable approx. is found */
X
X /* Obtain the gold section step */
X new_step = r * ( x= tol_act ) /* If x and w are distinct */
X { /* interpolatiom may be tried */
X register double p; /* Interpolation step is calcula-*/
X register double q; /* ted as p/q; division operation*/
X /* is delayed until last moment */
X register double t;
X
X t = (x-w) * (fx-fv);
X q = (x-v) * (fx-fw);
X p = (x-v)*q - (x-w)*t;
X q = 2*(q-t);
X
X if( q>(double)0 ) /* q was calculated with the op-*/
X p = -p; /* posite sign; make q positive */
X else /* and assign possible minus to */
X q = -q; /* p */
X
X if( fabs(p) < fabs(new_step*q) && /* If x+p/q falls in [a,b]*/
X p > q*(a-x+2*tol_act) && /* not too close to a and */
X p < q*(b-x-2*tol_act) ) /* b, and isn't too large */
X new_step = p/q; /* it is accepted */
X /* If p/q is too large then the */
X /* gold section procedure can */
X /* reduce [a,b] range to more */
X /* extent */
X }
X
X if( fabs(new_step) < tol_act ) /* Adjust the step to be not less*/
X if( new_step > (double)0 ) /* than tolerance */
X new_step = tol_act;
X else
X new_step = -tol_act;
X
X /* Obtain the next approximation to min */
X { /* and reduce the enveloping range */
X register double t = x + new_step; /* Tentative point for the min */
X register double ft = (*f)(t);
X if( ft <= fx )
X { /* t is a better approximation */
X if( t < x ) /* Reduce the range so that */
X b = x; /* t would fall within it */
X else
X a = x;
X
X v = w; w = x; x = t; /* Assign the best approx to x */
X fv=fw; fw=fx; fx=ft;
X }
X else /* x remains the better approx */
X {
X if( t < x ) /* Reduce the range enclosing x */
X a = t;
X else
X b = t;
X
X if( ft <= fw || w==x )
X {
X v = w; w = t;
X fv=fw; fw=ft;
X }
X else if( ft<=fv || v==x || v==w )
X {
X v = t;
X fv=ft;
X }
X }
X
X } /* ----- end-of-block ----- */
X } /* ===== End of loop ===== */
X
X}
________This_Is_The_END________
if test `wc -l < fminbr.c` -ne 161; then
echo 'shar: fminbr.c was damaged during transit (should have had 161 lines)'
fi
echo 'x - zeroin.c'
sed 's/^X//' << '________This_Is_The_END________' > zeroin.c
X/*
X ************************************************************************
X * C math library
X * function ZEROIN - obtain a function zero within the given range
X *
X * Input
X * double zeroin(ax,bx,f,tol)
X * double ax; Root will be seeked for within
X * double bx; a range [ax,bx]
X * double (*f)(double x); Name of the function whose zero
X * will be seeked for
X * double tol; Acceptable tolerance for the root
X * value.
X * May be specified as 0.0 to cause
X * the program to find the root as
X * accurate as possible
X *
X * Output
X * Zeroin returns an estimate for the root with accuracy
X * 4*EPSILON*abs(x) + tol
X *
X * Algorithm
X * G.Forsythe, M.Malcolm, C.Moler, Computer methods for mathematical
X * computations. M., Mir, 1980, p.180 of the Russian edition
X *
X * The function makes use of the bissection procedure combined with
X * the linear or quadric inverse interpolation.
X * At every step program operates on three abscissae - a, b, and c.
X * b - the last and the best approximation to the root
X * a - the last but one approximation
X * c - the last but one or even earlier approximation than a that
X * 1) |f(b)| <= |f(c)|
X * 2) f(b) and f(c) have opposite signs, i.e. b and c confine
X * the root
X * At every step Zeroin selects one of the two new approximations, the
X * former being obtained by the bissection procedure and the latter
X * resulting in the interpolation (if a,b, and c are all different
X * the quadric interpolation is utilized, otherwise the linear one).
X * If the latter (i.e. obtained by the interpolation) point is
X * reasonable (i.e. lies within the current interval [b,c] not being
X * too close to the boundaries) it is accepted. The bissection result
X * is used in the other case. Therefore, the range of uncertainty is
X * ensured to be reduced at least by the factor 1.6
X *
X ************************************************************************
X */
X
X#include "math.h"
X
Xdouble zeroin(ax,bx,f,tol) /* An estimate to the root */
Xdouble ax; /* Left border | of the range */
Xdouble bx; /* Right border| the root is seeked*/
Xdouble (*f)(double x); /* Function under investigation */
Xdouble tol; /* Acceptable tolerance */
X{
X double a,b,c; /* Abscissae, descr. see above */
X double fa; /* f(a) */
X double fb; /* f(b) */
X double fc; /* f(c) */
X
X a = ax; b = bx; fa = (*f)(a); fb = (*f)(b);
X c = a; fc = fa;
X
X for(;;) /* Main iteration loop */
X {
X double prev_step = b-a; /* Distance from the last but one*/
X /* to the last approximation */
X double tol_act; /* Actual tolerance */
X double p; /* Interpolation step is calcu- */
X double q; /* lated in the form p/q; divi- */
X /* sion operations is delayed */
X /* until the last moment */
X double new_step; /* Step at this iteration */
X
X if( fabs(fc) < fabs(fb) )
X { /* Swap data for b to be the */
X a = b; b = c; c = a; /* best approximation */
X fa=fb; fb=fc; fc=fa;
X }
X tol_act = 2*EPSILON*fabs(b) + tol/2;
X new_step = (c-b)/2;
X
X if( fabs(new_step) <= tol_act || fb == (double)0 )
X return b; /* Acceptable approx. is found */
X
X /* Decide if the interpolation can be tried */
X if( fabs(prev_step) >= tol_act /* If prev_step was large enough*/
X && fabs(fa) > fabs(fb) ) /* and was in true direction, */
X { /* Interpolatiom may be tried */
X register double t1,cb,t2;
X cb = c-b;
X if( a==c ) /* If we have only two distinct */
X { /* points linear interpolation */
X t1 = fb/fa; /* can only be applied */
X p = cb*t1;
X q = 1.0 - t1;
X }
X else /* Quadric inverse interpolation*/
X {
X q = fa/fc; t1 = fb/fc; t2 = fb/fa;
X p = t2 * ( cb*q*(q-t1) - (b-a)*(t1-1.0) );
X q = (q-1.0) * (t1-1.0) * (t2-1.0);
X }
X if( p>(double)0 ) /* p was calculated with the op-*/
X q = -q; /* posite sign; make p positive */
X else /* and assign possible minus to */
X p = -p; /* q */
X
X if( p < (0.75*cb*q-fabs(tol_act*q)/2) /* If b+p/q falls in [b,c]*/
X && p < fabs(prev_step*q/2) ) /* and isn't too large */
X new_step = p/q; /* it is accepted */
X /* If p/q is too large then the */
X /* bissection procedure can */
X /* reduce [b,c] range to more */
X /* extent */
X }
X
X if( fabs(new_step) < tol_act ) /* Adjust the step to be not less*/
X if( new_step > (double)0 ) /* than tolerance */
X new_step = tol_act;
X else
X new_step = -tol_act;
X
X a = b; fa = fb; /* Save the previous approx. */
X b += new_step; fb = (*f)(b); /* Do step to a new approxim. */
X if( (fb > 0 && fc > 0) || (fb < 0 && fc < 0) )
X { /* Adjust c for it to have a sign*/
X c = a; fc = fa; /* opposite to that of b */
X }
X }
X
X}
________This_Is_The_END________
if test `wc -l < zeroin.c` -ne 132; then
echo 'shar: zeroin.c was damaged during transit (should have had 132 lines)'
fi
echo 'x - vfminbr.c'
sed 's/^X//' << '________This_Is_The_END________' > vfminbr.c
X/*
X *=======================================================================
X * Verify FMINBR routine
X */
X
X#include "math.h"
X
Xdouble fminbr(double a, double b, double (*f)(), double tol);
X
Xstatic int counter; /* Iteration counter */
X
Xtest(a,b,f,msg) /* Run a test */
Xdouble a,b; /* Range the root is seeked for */
Xdouble (*f)(double x); /* Functiom under examination */
Xchar * msg; /* Explanation message */
X{
X double minloc;
X counter = 0;
X printf("\nFor function %s\nin [%g,%g] min found is at\t%.9e\n",msg,a,b,
X (minloc=fminbr(a,b,f,EPSILON)) );
X printf("Min function value found\t%.4e\nNo. of iterations\t\t%d\n",
X (*f)(minloc), counter);
X}
X
Xdouble f1(x) /* Test of the Forsythe book */
Xdouble x;
X{
X counter++;
X return (powi(x,2)-2.0)*x - 5.0;
X}
X
Xdouble f2(x) /* Modified test 1 */
Xdouble x;
X{
X counter++;
X return powi( (powi(x,2)-2.0)*x - 5.0, 2 );
X}
X
Xdouble f3(x) /* My test */
Xdouble x;
X{
X counter++;
X return powi( cos(x) - x,2 ) - 2;
X}
X
Xdouble f4(x) /* My test */
Xdouble x;
X{
X counter++;
X return powi( sin(x) - x,2 ) + 1;
X}
X
X
Xmain()
X{
X test(0.0,1.0,f1,"x^3 - 2*x - 5");
X printf("Exact min is at\t\t0.81650\n");
X
X test(2.0,3.0,f2,"(x^3 - 2*x - 5)^2");
X printf("Exact root is \t\t2.0945514815\n");
X
X test(2.0,3.0,f3,"(cos(x)-x)^2 - 2");
X test(-1.0,3.0,f3,"(cos(x)-x)^2 - 2");
X test(-1.0,3.0,f4,"(sin(x)-x)^2 + 1");
X}
________This_Is_The_END________
if test `wc -l < vfminbr.c` -ne 65; then
echo 'shar: vfminbr.c was damaged during transit (should have had 65 lines)'
fi
echo 'x - vzeroin.c'
sed 's/^X//' << '________This_Is_The_END________' > vzeroin.c
X/*
X *=======================================================================
X * Verify ZEROIN routine
X */
X
X#include "math.h"
X
Xdouble zeroin(double ax, double bx, double (*f)(), double tol);
X
Xstatic int counter; /* Iteration counter */
X
Xtest(a,b,f,msg) /* Run a test */
Xdouble a,b; /* Range the root is seeked for */
Xdouble (*f)(double x); /* Functiom under examination */
Xchar * msg; /* Explanation message */
X{
X double root;
X counter = 0;
X printf("\nFor function %s\nin [%g,%g] root is\t%.9e\n",msg,a,b,
X (root=zeroin(a,b,f,0.0)) );
X printf("Function value at the root found\t%.4e\nNo. of iterations\t\t%d\n",
X (*f)(root), counter);
X}
X
Xdouble f1(x) /* Test of the Forsythe book */
Xdouble x;
X{
X counter++;
X return (pow(x,2)-2.0)*x - 5.0;
X}
X
Xdouble f2(x) /* My test */
Xdouble x;
X{
X counter++;
X return cos(x) - x;
X}
X
Xdouble f3(x) /* My test */
Xdouble x;
X{
X counter++;
X return sin(x) - x;
X}
X
X
Xmain()
X{
X test(2.0,3.0,f1,"x^3 - 2*x - 5");
X printf("Exact root is \t\t2.0945514815\n");
X
X test(2.0,3.0,f2,"cos(x)-x");
X test(-1.0,3.0,f2,"cos(x)-x");
X test(-1.0,3.0,f3,"sin(x)-x");
X}
________This_Is_The_END________
if test `wc -l < vzeroin.c` -ne 55; then
echo 'shar: vzeroin.c was damaged during transit (should have had 55 lines)'
fi
echo 'x - vfminbr.dat'
sed 's/^X//' << '________This_Is_The_END________' > vfminbr.dat
X
XFor function x^3 - 2*x - 5
Xin [0,1] min found is at 8.164965811e-01
XMin function value found -6.0887e+00
XNo. of iterations 12
XExact min is at 0.81650
X
XFor function (x^3 - 2*x - 5)^2
Xin [2,3] min found is at 2.094551483e+00
XMin function value found 2.7186e-16
XNo. of iterations 15
XExact root is 2.0945514815
X
XFor function (cos(x)-x)^2 - 2
Xin [2,3] min found is at 2.000000048e+00
XMin function value found 3.8378e+00
XNo. of iterations 36
X
XFor function (cos(x)-x)^2 - 2
Xin [-1,3] min found is at 7.390851269e-01
XMin function value found -2.0000e+00
XNo. of iterations 18
X
XFor function (sin(x)-x)^2 + 1
Xin [-1,3] min found is at -3.815499176e-03
XMin function value found 1.0000e+00
XNo. of iterations 60
________This_Is_The_END________
if test `wc -l < vfminbr.dat` -ne 27; then
echo 'shar: vfminbr.dat was damaged during transit (should have had 27 lines)'
fi
echo 'x - vzeroin.dat'
sed 's/^X//' << '________This_Is_The_END________' > vzeroin.dat
X
XFor function x^3 - 2*x - 5
Xin [2,3] root is 2.094551482e+00
XFunction value at the root found -1.7764e-15
XNo. of iterations 11
XExact root is 2.0945514815
X
XFor function cos(x)-x
Xin [2,3] root is 2.000000000e+00
XFunction value at the root found -2.4161e+00
XNo. of iterations 52
X
XFor function cos(x)-x
Xin [-1,3] root is 7.390851332e-01
XFunction value at the root found 0.0000e+00
XNo. of iterations 11
X
XFor function sin(x)-x
Xin [-1,3] root is -1.643737357e-08
XFunction value at the root found 0.0000e+00
XNo. of iterations 58
________This_Is_The_END________
if test `wc -l < vzeroin.dat` -ne 21; then
echo 'shar: vzeroin.dat was damaged during transit (should have had 21 lines)'
fi
echo 'x - vfminbr.dat.vax'
sed 's/^X//' << '________This_Is_The_END________' > vfminbr.dat.vax
X
XFor function x^3 - 2*x - 5
Xin [0,1] min found is at 8.164965778e-01
XMin function value found -6.0887e+00
XNo. of iterations 12
XExact min is at 0.81650
X
XFor function (x^3 - 2*x - 5)^2
Xin [2,3] min found is at 2.094551483e+00
XMin function value found 2.7186e-16
XNo. of iterations 15
XExact root is 2.0945514815
X
XFor function (cos(x)-x)^2 - 2
Xin [2,3] min found is at 2.000000011e+00
XMin function value found 3.8378e+00
XNo. of iterations 38
X
XFor function (cos(x)-x)^2 - 2
Xin [-1,3] min found is at 7.390851332e-01
XMin function value found -2.0000e+00
XNo. of iterations 18
X
XFor function (sin(x)-x)^2 + 1
Xin [-1,3] min found is at 2.859282036e-03
XMin function value found 1.0000e+00
XNo. of iterations 56
________This_Is_The_END________
if test `wc -l < vfminbr.dat.vax` -ne 27; then
echo 'shar: vfminbr.dat.vax was damaged during transit (should have had 27 lines)'
fi
echo 'x - vzeroin.dat.vax'
sed 's/^X//' << '________This_Is_The_END________' > vzeroin.dat.vax
X
XFor function x^3 - 2*x - 5
Xin [2,3] root is 2.094551482e+00
XFunction value at the root found 3.3307e-16
XNo. of iterations 10
XExact root is 2.0945514815
X
XFor function cos(x)-x
Xin [2,3] root is 2.000000000e+00
XFunction value at the root found -2.4161e+00
XNo. of iterations 55
X
XFor function cos(x)-x
Xin [-1,3] root is 7.390851332e-01
XFunction value at the root found 0.0000e+00
XNo. of iterations 10
X
XFor function sin(x)-x
Xin [-1,3] root is -5.916737541e-09
XFunction value at the root found 0.0000e+00
XNo. of iterations 59
________This_Is_The_END________
if test `wc -l < vzeroin.dat.vax` -ne 21; then
echo 'shar: vzeroin.dat.vax was damaged during transit (should have had 21 lines)'
fi