src/math/math_BrentMinimum.cxx

   1 //static const char* sccsid = "@(#)math_BrentMinimum.cxx        3.2 95/01/10"; // Do not delete this line. Used by sccs.
   2 #include <math_BrentMinimum.ixx>
   3 #include <math_Function.hxx>
   4
   5 #define CGOLD         0.3819660
   6 #ifdef MAX
   7 #undef MAX
   8 #endif
   9 #define MAX(a,b)      ((a) > (b) ? (a) : (b))
  10 #define SIGN(a,b)     ((b) > 0.0 ? fabs(a) : -fabs(a))
  11 #define SHFT(a,b,c,d) (a)=(b);(b)=(c);(c)=(d)
  12
  13 void math_BrentMinimum::Perform(math_Function& F,
  14                                 const Standard_Real    ax,
  15                                 const Standard_Real    bx,
  16                                 const Standard_Real    cx) {
  17
  18   Standard_Boolean OK;
  19   Standard_Real etemp, fu, p, q, r;
  20   Standard_Real tol1, tol2, u, v, w, xm;
  21   Standard_Real e = 0.0;
  22   Standard_Real d = RealLast();
  23
  24   a = ((ax < cx) ? ax : cx);
  25   b = ((ax > cx) ? ax : cx);
  26   x = w = v = bx;
  27   if (!myF) {
  28     OK = F.Value(x, fx);
  29     if(!OK) return;
  30   }
  31   fw = fv = fx;
  32   for(iter = 1; iter <= Itermax; iter++) {
  33     xm = 0.5 * (a + b);
  34     tol1 = XTol * fabs(x) + EPSZ;
  35     tol2 = 2.0 * tol1;
  36     if(IsSolutionReached(F)) {
  37       Done = Standard_True;
  38       return;
  39     }
  40     if(fabs(e) > tol1) {
  41       r = (x - w) * (fx - fv);
  42       q = (x - v) * (fx - fw);
  43       p = (x - v) * q - (x - w) * r;
  44       q = 2.0 * (q - r);
  45       if(q > 0.0) p = -p;
  46       q = fabs(q);
  47       etemp = e;
  48       e = d;
  49       if(fabs(p) >= fabs(0.5 * q * etemp)
  50          || p <= q * ( a - x) || p >= q * (b - x)) {
  51         e = (x >= xm ? a - x : b - x);
  52         d = CGOLD * e;
  53       }
  54       else {
  55         d = p / q;
  56         u = x + d;
  57         if(u - a < tol2 || b - u < tol2) d = SIGN(tol1, xm - x);
  58       }
  59     }
  60     else {
  61       e = (x >= xm ? a - x : b - x);
  62       d = CGOLD * e;
  63     }
  64     u = (fabs(d) >= tol1 ? x + d : x + SIGN(tol1, d));
  65     OK = F.Value(u, fu);
  66     if(!OK) return;
  67     if(fu <= fx) {
  68       if(u >= x) a = x; else b = x;
  69       SHFT(v, w, x, u);
  70       SHFT(fv, fw, fx, fu);
  71     }
  72     else {
  73       if(u < x) a = u; else b = u;
  74       if(fu <= fw || w == x) {
  75         v = w;
  76         w = u;
  77         fv = fw;
  78         fw = fu;
  79       }
  80       else if(fu <= fv || v == x || v == w) {
  81         v = u;
  82         fv = fu;
  83       }
  84     }
  85   }
  86   Done = Standard_False;
  87   return;
  88 }
  89
  90
  91 math_BrentMinimum::math_BrentMinimum(math_Function& F,
  92                                      const Standard_Real    Ax,
  93                                      const Standard_Real    Bx,
  94                                      const Standard_Real    Cx,
  95                                      const Standard_Real    TolX,
  96                                      const Standard_Integer NbIterations,
  97                                      const Standard_Real    ZEPS) {
  98
  99   XTol = TolX;
 100   EPSZ = ZEPS;
 101   Itermax = NbIterations;
 102   myF = Standard_False;
 103   Perform(F, Ax, Bx, Cx);
 104 }
 105
 106
 107 // Constructeur d'initialisation des champs.
 108
 109 math_BrentMinimum::math_BrentMinimum(const Standard_Real    TolX,
 110                                      const Standard_Integer NbIterations,
 111                                      const Standard_Real    ZEPS) {
 112   myF = Standard_False;
 113   XTol = TolX;
 114   EPSZ = ZEPS;
 115   Itermax = NbIterations;
 116 }
 117
 118 math_BrentMinimum::math_BrentMinimum(const Standard_Real    TolX,
 119                                      const Standard_Real    Fbx,
 120                                      const Standard_Integer NbIterations,
 121                                      const Standard_Real    ZEPS) {
 122
 123   fx = Fbx;
 124   myF = Standard_True;
 125   XTol = TolX;
 126   EPSZ = ZEPS;
 127   Itermax = NbIterations;
 128 }
 129
 130
 131 //    Standard_Boolean math_BrentMinimum::IsSolutionReached(math_Function& F) {
 132     Standard_Boolean math_BrentMinimum::IsSolutionReached(math_Function& ) {
 133
 134        Standard_Real xm = 0.5 * (a + b);
 135        // modified by NIZHNY-MKK  Mon Oct  3 17:45:57 2005.BEGIN
 136 //        Standard_Real tol = XTol * fabs(x) + EPSZ;
 137 //        return fabs(x - xm) <= 2.0 * tol - 0.5 * (b - a);
 138        Standard_Real TwoTol = 2.0 *(XTol * fabs(x) + EPSZ);
 139        return ((x <= (TwoTol + a)) && (x >= (b - TwoTol)));
 140        // modified by NIZHNY-MKK  Mon Oct  3 17:46:00 2005.END
 141   }
 142
 143
 144
 145     void math_BrentMinimum::Dump(Standard_OStream& o) const {
 146        o << "math_BrentMinimum ";
 147        if(Done) {
 148          o << " Status = Done \n";
 149          o << " Location value = " << x <<"\n";
 150          o << " Minimum value = " << fx << "\n";
 151          o << " Number of iterations = " << iter <<"\n";
 152        }
 153        else {
 154          o << " Status = not Done \n";
 155        }
 156     }