0027252: Implicit-implicit intersection (Cylinder-Plane) loses intersection curve
[occt.git] / src / math / math_BrentMinimum.cxx
CommitLineData
b311480e 1// Copyright (c) 1997-1999 Matra Datavision
973c2be1 2// Copyright (c) 1999-2014 OPEN CASCADE SAS
b311480e 3//
973c2be1 4// This file is part of Open CASCADE Technology software library.
b311480e 5//
d5f74e42 6// This library is free software; you can redistribute it and/or modify it under
7// the terms of the GNU Lesser General Public License version 2.1 as published
973c2be1 8// by the Free Software Foundation, with special exception defined in the file
9// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
10// distribution for complete text of the license and disclaimer of any warranty.
b311480e 11//
973c2be1 12// Alternatively, this file may be used under the terms of Open CASCADE
13// commercial license or contractual agreement.
b311480e 14
42cf5bc1 15
16#include <math_BrentMinimum.hxx>
7fd59977 17#include <math_Function.hxx>
42cf5bc1 18#include <StdFail_NotDone.hxx>
7fd59977 19
f542b7bb 20static const Standard_Real CGOLD = 0.3819660; //0.5*(3 - sqrt(5));
21
22
23//=======================================================================
24//function : SHFT
25//purpose : Shifts arguments
26//=======================================================================
27inline void SHFT(Standard_Real &theA, Standard_Real &theB,
28 Standard_Real &theC, Standard_Real &theD)
29{
30 theA = theB;
31 theB = theC;
32 theC = theD;
33}
7fd59977 34
859a47c3 35//=======================================================================
36//function : math_BrentMinimum
37//purpose : Constructor
38//=======================================================================
39math_BrentMinimum::math_BrentMinimum(const Standard_Real theTolX,
40 const Standard_Integer theNbIterations,
41 const Standard_Real theZEPS)
42: a (0.0),
43 b (0.0),
44 x (0.0),
45 fx (0.0),
46 fv (0.0),
47 fw (0.0),
48 XTol(theTolX),
49 EPSZ(theZEPS),
50 Done (Standard_False),
51 iter (0),
52 Itermax(theNbIterations),
53 myF (Standard_False)
54{
55}
56
57//=======================================================================
58//function : math_BrentMinimum
59//purpose : Constructor
60//=======================================================================
61math_BrentMinimum::math_BrentMinimum(const Standard_Real theTolX,
62 const Standard_Real theFbx,
63 const Standard_Integer theNbIterations,
64 const Standard_Real theZEPS)
65: a (0.0),
66 b (0.0),
67 x (0.0),
68 fx (theFbx),
69 fv (0.0),
70 fw (0.0),
71 XTol(theTolX),
72 EPSZ(theZEPS),
73 Done (Standard_False),
74 iter (0),
75 Itermax(theNbIterations),
76 myF (Standard_True)
77{
78}
79
80//=======================================================================
81//function : ~math_BrentMinimum
82//purpose : Destructor
83//=======================================================================
6da30ff1 84math_BrentMinimum::~math_BrentMinimum()
85{
86}
87
859a47c3 88//=======================================================================
89//function : Perform
90//purpose :
91//=======================================================================
7fd59977 92void math_BrentMinimum::Perform(math_Function& F,
859a47c3 93 const Standard_Real ax,
94 const Standard_Real bx,
95 const Standard_Real cx)
96{
7fd59977 97 Standard_Boolean OK;
98 Standard_Real etemp, fu, p, q, r;
99 Standard_Real tol1, tol2, u, v, w, xm;
100 Standard_Real e = 0.0;
101 Standard_Real d = RealLast();
f542b7bb 102
7fd59977 103 a = ((ax < cx) ? ax : cx);
104 b = ((ax > cx) ? ax : cx);
105 x = w = v = bx;
106 if (!myF) {
107 OK = F.Value(x, fx);
f542b7bb 108 if (!OK) return;
7fd59977 109 }
110 fw = fv = fx;
f542b7bb 111 for (iter = 1; iter <= Itermax; iter++) {
7fd59977 112 xm = 0.5 * (a + b);
113 tol1 = XTol * fabs(x) + EPSZ;
114 tol2 = 2.0 * tol1;
f542b7bb 115 if (IsSolutionReached(F)) {
7fd59977 116 Done = Standard_True;
117 return;
118 }
f542b7bb 119 if (fabs(e) > tol1) {
7fd59977 120 r = (x - w) * (fx - fv);
121 q = (x - v) * (fx - fw);
122 p = (x - v) * q - (x - w) * r;
123 q = 2.0 * (q - r);
f542b7bb 124 if (q > 0.0) p = -p;
7fd59977 125 q = fabs(q);
126 etemp = e;
127 e = d;
f542b7bb 128 if (fabs(p) >= fabs(0.5 * q * etemp)
129 || p <= q * (a - x) || p >= q * (b - x)) {
130 e = (x >= xm ? a - x : b - x);
131 d = CGOLD * e;
7fd59977 132 }
133 else {
f542b7bb 134 d = p / q;
135 u = x + d;
136 if (u - a < tol2 || b - u < tol2) d = Sign(tol1, xm - x);
7fd59977 137 }
138 }
139 else {
140 e = (x >= xm ? a - x : b - x);
141 d = CGOLD * e;
142 }
f542b7bb 143 u = (fabs(d) >= tol1 ? x + d : x + Sign(tol1, d));
7fd59977 144 OK = F.Value(u, fu);
f542b7bb 145 if (!OK) return;
146 if (fu <= fx) {
147 if (u >= x) a = x; else b = x;
7fd59977 148 SHFT(v, w, x, u);
149 SHFT(fv, fw, fx, fu);
150 }
151 else {
f542b7bb 152 if (u < x) a = u; else b = u;
153 if (fu <= fw || w == x) {
154 v = w;
155 w = u;
156 fv = fw;
157 fw = fu;
7fd59977 158 }
f542b7bb 159 else if (fu <= fv || v == x || v == w) {
160 v = u;
161 fv = fu;
7fd59977 162 }
163 }
164 }
165 Done = Standard_False;
166 return;
167}
168
859a47c3 169//=======================================================================
170//function : Dump
171//purpose :
172//=======================================================================
173void math_BrentMinimum::Dump(Standard_OStream& o) const
174{
175 o << "math_BrentMinimum ";
176 if(Done) {
177 o << " Status = Done \n";
178 o << " Location value = " << x <<"\n";
179 o << " Minimum value = " << fx << "\n";
180 o << " Number of iterations = " << iter <<"\n";
7fd59977 181 }
859a47c3 182 else {
183 o << " Status = not Done \n";
184 }
185}