src/math/math_BracketMinimum.cxx

   1 //static const char* sccsid = "@(#)math_BracketMinimum.cxx      3.2 95/01/10"; // Do not delete this line. Used by sccs.
   2 #include <math_BracketMinimum.ixx>
   3
   4 #include <StdFail_NotDone.hxx>   // waiting for NotDone Exception
   5 #include <math_Function.hxx>
   6
   7 #define GOLD           1.618034
   8 #define CGOLD          0.3819660
   9 #define GLIMIT         100.0
  10 #define TINY           1.0e-20
  11 #ifdef MAX
  12 #undef MAX
  13 #endif
  14 #define MAX(a,b)       ((a) > (b) ? (a) : (b))
  15 #define SIGN(a,b)      ((b) > 0.0 ? fabs(a) : -fabs(a))
  16 #define SHFT(a,b,c,d)  (a)=(b);(b)=(c);(c)=(d)
  17
  18
  19     void math_BracketMinimum::Perform(math_Function& F,
  20                                       const Standard_Real A,
  21                                       const Standard_Real B) {
  22
  23      Standard_Boolean OK;
  24      Standard_Real ulim, u, r, q, f, fu, dum;
  25
  26      Done = Standard_False;
  27      Ax = A;
  28      Bx = B;
  29      Standard_Real Lambda = GOLD;
  30      if (!myFA) {
  31        OK = F.Value(Ax, FAx);
  32        if(!OK) return;
  33      }
  34      if (!myFB) {
  35        OK = F.Value(Bx, FBx);
  36        if(!OK) return;
  37      }
  38      if(FBx > FAx) {
  39        SHFT(dum, Ax, Bx, dum);
  40        SHFT(dum, FBx, FAx, dum);
  41      }
  42      Cx = Bx + Lambda * (Bx - Ax);
  43      OK = F.Value(Cx, FCx);
  44      if(!OK) return;
  45      while(FBx > FCx) {
  46        r = (Bx - Ax) * (FBx -FCx);
  47        q = (Bx - Cx) * (FBx -FAx);
  48        u = Bx - ((Bx - Cx) * q - (Bx - Ax) * r) /
  49            (2.0 * SIGN(MAX(fabs(q - r), TINY), q - r));
  50        ulim = Bx + GLIMIT * (Cx - Bx);
  51        if((Bx - u) * (u - Cx) > 0.0) {
  52          OK = F.Value(u, fu);
  53          if(!OK) return;
  54          if(fu < FCx) {
  55            Ax = Bx;
  56            Bx = u;
  57            FAx = FBx;
  58            FBx = fu;
  59            Done = Standard_True;
  60            return;
  61          }
  62          else if(fu > FBx) {
  63            Cx = u;
  64            FCx = fu;
  65            Done = Standard_True;
  66            return;
  67          }
  68          u = Cx + Lambda * (Cx - Bx);
  69          OK = F.Value(u, fu);
  70          if(!OK) return;
  71        }
  72        else if((Cx - u) * (u - ulim) > 0.0) {
  73          OK = F.Value(u, fu);
  74          if(!OK) return;
  75          if(fu < FCx) {
  76            SHFT(Bx, Cx, u, Cx + GOLD * (Cx - Bx));
  77            OK = F.Value(u, f);
  78            if(!OK) return;
  79            SHFT(FBx, FCx, fu, f);
  80          }
  81        }
  82        else if ((u - ulim) * (ulim - Cx) >= 0.0) {
  83          u = ulim;
  84          OK = F.Value(u, fu);
  85          if(!OK) return;
  86        }
  87        else {
  88          u = Cx + GOLD * (Cx - Bx);
  89          OK = F.Value(u, fu);
  90          if(!OK) return;
  91        }
  92        SHFT(Ax, Bx, Cx, u);
  93        SHFT(FAx, FBx, FCx, fu);
  94      }
  95      Done = Standard_True;
  96    }
  97
  98
  99
 100
 101     math_BracketMinimum::math_BracketMinimum(math_Function& F,
 102                                              const Standard_Real A,
 103                                              const Standard_Real B) {
 104
 105       myFA = Standard_False;
 106       myFB = Standard_False;
 107       Perform(F, A, B);
 108     }
 109
 110     math_BracketMinimum::math_BracketMinimum(math_Function& F,
 111                                              const Standard_Real A,
 112                                              const Standard_Real B,
 113                                              const Standard_Real FA) {
 114       FAx = FA;
 115       myFA = Standard_True;
 116       myFB = Standard_False;
 117       Perform(F, A, B);
 118     }
 119
 120     math_BracketMinimum::math_BracketMinimum(math_Function& F,
 121                                              const Standard_Real A,
 122                                              const Standard_Real B,
 123                                              const Standard_Real FA,
 124                                              const Standard_Real FB) {
 125       FAx = FA;
 126       FBx = FB;
 127       myFA = Standard_True;
 128       myFB = Standard_True;
 129       Perform(F, A, B);
 130     }
 131
 132
 133     void math_BracketMinimum::Values(Standard_Real& A, Standard_Real& B, Standard_Real& C) const{
 134
 135       StdFail_NotDone_Raise_if(!Done, " ");
 136       A = Ax;
 137       B = Bx;
 138       C = Cx;
 139     }
 140
 141     void math_BracketMinimum::FunctionValues(Standard_Real& FA, Standard_Real& FB, Standard_Real& FC) const{
 142
 143       StdFail_NotDone_Raise_if(!Done, " ");
 144       FA = FAx;
 145       FB = FBx;
 146       FC = FCx;
 147     }
 148
 149     void math_BracketMinimum::Dump(Standard_OStream& o) const {
 150
 151        o << "math_BracketMinimum ";
 152        if(Done) {
 153          o << " Status = Done \n";
 154          o << " The bracketed triplet is: " << endl;
 155          o << Ax << ", " << Bx << ", " << Cx << endl;
 156          o << " The corresponding function values are: "<< endl;
 157          o << FAx << ", " << FBx << ", " << FCx << endl;
 158        }
 159        else {
 160          o << " Status = not Done \n";
 161        }
 162 }
 163