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