0024510: Remove unused local variables
[occt.git] / src / math / math_BracketedRoot.cxx
1 // Copyright (c) 1997-1999 Matra Datavision
2 // Copyright (c) 1999-2014 OPEN CASCADE SAS
3 //
4 // This file is part of Open CASCADE Technology software library.
5 //
6 // This library is free software; you can redistribute it and / or modify it
7 // under the terms of the GNU Lesser General Public version 2.1 as published
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.
11 //
12 // Alternatively, this file may be used under the terms of Open CASCADE
13 // commercial license or contractual agreement.
14
15 #include <math_BracketedRoot.ixx>
16 #include <math_Function.hxx>
17
18 // reference algorithme:  
19 //                   Brent method 
20 //                   numerical recipes in C  p 269
21
22 math_BracketedRoot::math_BracketedRoot (math_Function& F, 
23                                            const Standard_Real Bound1, 
24                                            const Standard_Real Bound2, 
25                                            const Standard_Real Tolerance, 
26                                            const Standard_Integer NbIterations,
27                                            const Standard_Real ZEPS ) {
28
29     Standard_Real Fa,Fc,a,c=0,d=0,e=0;
30     Standard_Real min1,min2,p,q,r,s,tol1,xm;
31   
32     a = Bound1;
33     TheRoot = Bound2;
34     F.Value(a,Fa);
35     F.Value(TheRoot,TheError);
36     if (Fa*TheError > 0.) { Done = Standard_False;}
37     else {
38       Fc = TheError ;
39       for (NbIter = 1; NbIter <= NbIterations; NbIter++) {
40           if (TheError*Fc > 0.) {
41              c = a;      // rename a TheRoot c and adjust bounding interval d
42              Fc = Fa;
43              d = TheRoot - a;
44              e = d;
45           } 
46           if ( Abs(Fc) < Abs(Fa) ) {
47              a = TheRoot;
48              TheRoot = c;
49              c = a;
50              Fa = TheError;
51              TheError = Fc;
52              Fc = Fa;
53           }
54           tol1 = 2.*ZEPS * Abs(TheRoot) + 0.5 * Tolerance; // convergence check
55           xm = 0.5 * ( c - TheRoot );
56           if (Abs(xm) <= tol1 || TheError == 0. ) {
57                Done = Standard_True;
58                return;
59           }
60           if (Abs(e) >= tol1 && Abs(Fa) > Abs(TheError) ) {
61              s = TheError / Fa; // attempt inverse quadratic interpolation
62              if (a == c) {
63                 p = 2.*xm*s;
64                 q = 1. - s;
65              }
66              else {
67                 q = Fa / Fc;
68                 r = TheError / Fc;
69                 p = s * (2.*xm *q * (q - r) - (TheRoot - a)*(r - 1.)); 
70                 q = (q -1.) * (r - 1.) * (s - 1.);
71              }
72             if ( p > 0. ) { q = -q;} // check whether in bounds
73             p = Abs(p);
74             min1 = 3.* xm* q - Abs(tol1 *q);
75             min2 = Abs(e * q);
76             if (2.* p < (min1 < min2 ? min1 : min2) ) {
77                e = d ;  // accept interpolation
78                d = p / q;
79             }
80             else {
81                d = xm;  // interpolation failed,use bissection
82                e = d;
83             }
84           }
85           else {   // bounds decreasing too slowly ,use bissection
86               d = xm;
87               e =d;            
88           }
89           a = TheRoot ;   // move last best guess to a
90           Fa = TheError;
91           if (Abs(d) > tol1) {  // evaluate new trial root
92              TheRoot += d;
93           }
94           else {
95              TheRoot += (xm > 0. ? Abs(tol1) : -Abs(tol1));
96           }
97           F.Value(TheRoot,TheError);
98       }  
99      Done = Standard_False;
100     }  
101   }
102
103
104     void math_BracketedRoot::Dump(Standard_OStream& o) const {
105
106        o << "math_BracketedRoot ";
107        if(Done) {
108          o << " Status = Done \n";
109          o << " Number of iterations = " << NbIter << endl;
110          o << " The Root is: " << TheRoot << endl;
111          o << " The value at the root is: " << TheError << endl;
112        }
113        else {
114          o << " Status = not Done \n";
115        }
116 }