Main Page   Class Hierarchy   Compound List   File List   Compound Members   File Members  

brent.c

Go to the documentation of this file.
00001 
00002 /*
00003  *   THIS FILE IS UNDER RCS - DO NOT MODIFY UNLESS YOU HAVE
00004  *   CHECKED IT OUT USING THE COMMAND CHECKOUT.
00005  *
00006  *    $Id: brent_8c-source.html 2161 2006-05-19 16:55:03Z paulf $
00007  *
00008  *    Revision history:
00009  *     $Log$
00009  *     Revision 1.1  2006/05/19 16:55:01  paulf
00009  *     first inclusion
00009  *
00010  *     Revision 1.1  2000/02/14 18:51:48  lucky
00011  *     Initial revision
00012  *
00013  *
00014  */
00015 
00016 #include <math.h>
00017 
00018 #define ITMAX 100
00019 #define CGOLD 0.3819660
00020 #define ZEPS 1.0e-10
00021 #define SIGN(a,b) ((b) > 0.0 ? fabs(a) : -fabs(a))
00022 #define SHFT(a,b,c,d) (a)=(b);(b)=(c);(c)=(d);
00023 
00024 float brent(float ax, float bx, float cx, float (*f)(float), float tol, float *xmin)
00025 {
00026         int iter;
00027         float a,b,etemp,fu,fv,fw,fx,p,q,r,tol1,tol2,u,v,w,x,xm;
00028         float d=0.0;
00029         float e=0.0;
00030 
00031         a=((ax < cx) ? ax : cx);
00032         b=((ax > cx) ? ax : cx);
00033         x=w=v=bx;
00034         fw=fv=fx=(*f)(x);
00035 
00036         for (iter=1;iter<=ITMAX;iter++)
00037         {
00038            xm=(float)0.5*(a+b);
00039 
00040 
00041            tol2=(float)2.0*(tol1=(float)tol*(float)fabs(x)+(float)ZEPS);
00042 
00043            if (fabs(x-xm) <= (tol2-0.5*(b-a)))
00044            {
00045               *xmin=x;
00046               return fx;
00047            }
00048            if (fabs(e) > tol1)
00049            {
00050               r=(x-w)*(fx-fv);
00051               q=(x-v)*(fx-fw);
00052               p=(x-v)*q-(x-w)*r;
00053               q=(float)2.0*(q-r);
00054               if (q > 0.0) p = -p;
00055               q=(float)fabs(q);
00056               etemp=e;
00057               e=d;
00058               if (fabs(p) >= fabs(0.5*q*etemp) || p <= q*(a-x) || p >= q*(b-x))
00059                  d=(float)CGOLD*(e=(x >= xm ? a-x : b-x));
00060               else
00061               {
00062                  d=p/q;
00063                  u=x+d;
00064                  if (u-a < tol2 || b-u < tol2)
00065                     d=(float)SIGN(tol1,xm-x);
00066               }
00067            }
00068            else
00069               d=(float)CGOLD*(e=(x >= xm ? a-x : b-x));
00070 
00071            u=(float)(fabs(d) >= tol1 ? x+d : x+SIGN(tol1,d));
00072            fu=(*f)(u);
00073            if (fu <= fx)
00074            {
00075               if (u >= x) a=x; else b=x;
00076               SHFT(v,w,x,u)
00077               SHFT(fv,fw,fx,fu)
00078            }
00079            else
00080            {
00081               if (u < x) a=u; else b=u;
00082               if (fu <= fw || w == x)
00083               {
00084                  v=w;
00085                  w=u;
00086                  fv=fw;
00087                  fw=fu;
00088               }
00089               else if (fu <= fv || v == x || v == w)
00090               {
00091                  v=u;
00092                  fv=fu;
00093               }
00094            }
00095         }
00096         *xmin=x;
00097         return fx;
00098 }
00099 
00100 #undef ITMAX
00101 #undef CGOLD
00102 #undef ZEPS
00103 #undef SIGN

Generated on Tue May 6 09:16:00 2003 for Earthworm Libs by doxygen1.3-rc3