src/math/math_NewtonFunctionSetRoot.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 under
   7 // the terms of the GNU Lesser General Public License 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 //#ifndef OCCT_DEBUG
  16 #define No_Standard_RangeError
  17 #define No_Standard_OutOfRange
  18 #define No_Standard_DimensionError
  19
  20 //#endif
  21
  22 #include <math_FunctionSetWithDerivatives.hxx>
  23 #include <math_Matrix.hxx>
  24 #include <math_NewtonFunctionSetRoot.hxx>
  25 #include <math_Recipes.hxx>
  26 #include <Standard_DimensionError.hxx>
  27 #include <StdFail_NotDone.hxx>
  28
  29 //=======================================================================
  30 //function : math_NewtonFunctionSetRoot
  31 //purpose  : Constructor
  32 //=======================================================================
  33 math_NewtonFunctionSetRoot::math_NewtonFunctionSetRoot(
  34   math_FunctionSetWithDerivatives& theFunction,
  35   const math_Vector&               theXTolerance,
  36   const Standard_Real              theFTolerance,
  37   const Standard_Integer           theNbIterations)
  38
  39 : TolX    (1, theFunction.NbVariables()),
  40   TolF    (theFTolerance),
  41   Indx    (1, theFunction.NbVariables()),
  42   Scratch (1, theFunction.NbVariables()),
  43   Sol     (1, theFunction.NbVariables()),
  44   DeltaX  (1, theFunction.NbVariables()),
  45   FValues (1, theFunction.NbVariables()),
  46   Jacobian(1, theFunction.NbVariables(), 1, theFunction.NbVariables()),
  47   Done    (Standard_False),
  48   State   (0),
  49   Iter    (0),
  50   Itermax (theNbIterations)
  51 {
  52   SetTolerance(theXTolerance);
  53 }
  54
  55 //=======================================================================
  56 //function : math_NewtonFunctionSetRoot
  57 //purpose  : Constructor
  58 //=======================================================================
  59 math_NewtonFunctionSetRoot::math_NewtonFunctionSetRoot(
  60   math_FunctionSetWithDerivatives& theFunction,
  61   const Standard_Real              theFTolerance,
  62   const Standard_Integer           theNbIterations)
  63
  64 : TolX    (1, theFunction.NbVariables()),
  65   TolF    (theFTolerance),
  66   Indx    (1, theFunction.NbVariables()),
  67   Scratch (1, theFunction.NbVariables()),
  68   Sol     (1, theFunction.NbVariables()),
  69   DeltaX  (1, theFunction.NbVariables()),
  70   FValues (1, theFunction.NbVariables()),
  71   Jacobian(1, theFunction.NbVariables(), 1, theFunction.NbVariables()),
  72   Done    (Standard_False),
  73   State   (0),
  74   Iter    (0),
  75   Itermax (theNbIterations)
  76 {
  77 }
  78
  79 //=======================================================================
  80 //function : ~math_NewtonFunctionSetRoot
  81 //purpose  : Destructor
  82 //=======================================================================
  83 math_NewtonFunctionSetRoot::~math_NewtonFunctionSetRoot()
  84 {
  85 }
  86
  87 //=======================================================================
  88 //function : SetTolerance
  89 //purpose  :
  90 //=======================================================================
  91 void math_NewtonFunctionSetRoot::SetTolerance(const math_Vector& theXTolerance)
  92 {
  93   for (Standard_Integer i = 1; i <= TolX.Length(); ++i)
  94     TolX(i) = theXTolerance(i);
  95 }
  96
  97 //=======================================================================
  98 //function : Perform
  99 //purpose  :
 100 //=======================================================================
 101 void math_NewtonFunctionSetRoot::Perform(
 102   math_FunctionSetWithDerivatives& theFunction,
 103   const math_Vector&               theStartingPoint)
 104 {
 105   const math_Vector anInf(1, theFunction.NbVariables(), RealFirst());
 106   const math_Vector aSup (1, theFunction.NbVariables(), RealLast ());
 107
 108   Perform(theFunction, theStartingPoint, anInf, aSup);
 109 }
 110
 111 //=======================================================================
 112 //function : Perform
 113 //purpose  :
 114 //=======================================================================
 115 void math_NewtonFunctionSetRoot::Perform(
 116                            math_FunctionSetWithDerivatives& F,
 117                            const math_Vector& StartingPoint,
 118                            const math_Vector& InfBound,
 119                            const math_Vector& SupBound)
 120 {
 121
 122   Standard_Real d;
 123   Standard_Boolean OK;
 124   Standard_Integer Error;
 125
 126   Done = Standard_False;
 127   Sol = StartingPoint;
 128   OK = F.Values(Sol, FValues, Jacobian);
 129   if(!OK) return;
 130   for(Iter = 1; Iter <= Itermax; Iter++) {
 131     for(Standard_Integer k = 1; k <= DeltaX.Length(); k++) {
 132       DeltaX(k) = -FValues(k);
 133     }
 134     Error = LU_Decompose(Jacobian, Indx, d, Scratch, 1.0e-30);
 135     if(Error) return;
 136     LU_Solve(Jacobian, Indx, DeltaX);
 137     for(Standard_Integer i = 1; i <= Sol.Length(); i++) {
 138       Sol(i) += DeltaX(i);
 139
 140       // Limitation de Sol dans les bornes [InfBound, SupBound] :
 141       if (Sol(i) <= InfBound(i)) Sol(i) = InfBound(i);
 142       if (Sol(i) >= SupBound(i)) Sol(i) = SupBound(i);
 143
 144     }
 145     OK = F.Values(Sol, FValues, Jacobian);
 146     if(!OK) return;
 147     if(IsSolutionReached(F)) {
 148       State = F.GetStateNumber();
 149       Done = Standard_True;
 150       return;
 151     }
 152   }
 153 }
 154
 155 //=======================================================================
 156 //function : Dump
 157 //purpose  :
 158 //=======================================================================
 159 void math_NewtonFunctionSetRoot::Dump(Standard_OStream& o) const
 160 {
 161   o <<"math_NewtonFunctionSetRoot ";
 162   if (Done) {
 163     o << " Status = Done \n";
 164     o << " Vector solution = " << Sol <<"\n";
 165     o << " Value of the function at this solution = \n";
 166     o << FValues <<"\n";
 167     o << " Number of iterations = " << Iter <<"\n";
 168   }
 169   else {
 170     o << "Status = not Done \n";
 171   }
 172 }