src/math/math_FunctionSetRoot.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 // pmn 15/05/97 pas de Gauss avec des pivot trop petit. SVD fait mieux
  16 // l'affaire + limitation de la longeur du pas + qq comentaire issus d'EUCLID3
  17 // pmn 10/06/97 refonte totale du traitement des bords + ajustement des init
  18 // et des tolerances pour brent...
  19
  20 //#ifndef OCCT_DEBUG
  21 #define No_Standard_RangeError
  22 #define No_Standard_OutOfRange
  23 #define No_Standard_DimensionError
  24
  25 //#endif
  26 //math_FunctionSetRoot.cxx
  27
  28 #include <math_BrentMinimum.hxx>
  29 #include <math_Function.hxx>
  30 #include <math_FunctionSetRoot.hxx>
  31 #include <math_FunctionSetWithDerivatives.hxx>
  32 #include <math_Gauss.hxx>
  33 #include <math_GaussLeastSquare.hxx>
  34 #include <math_IntegerVector.hxx>
  35 #include <math_Matrix.hxx>
  36 #include <math_SVD.hxx>
  37 #include <Precision.hxx>
  38 #include <Standard_DimensionError.hxx>
  39 #include <StdFail_NotDone.hxx>
  40
  41 //===========================================================================
  42 // - A partir d une solution de depart, recherche d une direction.( Newton la
  43 // plupart du temps, gradient si Newton echoue.
  44 // - Recadrage au niveau des bornes avec recalcul de la direction si une
  45 // inconnue a une valeur imposee.
  46 // -Si On ne sort pas des bornes
  47 //   Tant que (On ne progresse pas assez ou on ne change pas de direction)
  48 //    . Si (Progression encore possible)
  49 //        Si (On ne sort pas des bornes)
  50 //          On essaye de progresser dans cette meme direction.
  51 //        Sinon
  52 //          On diminue le pas d'avancement ou on change de direction.
  53 //      Sinon
  54 //        Si on depasse le minimum
  55 //          On fait une interpolation parabolique.
  56 // - Si on a progresse sur F
  57 //     On fait les tests d'arret
  58 //   Sinon
  59 //     On change de direction
  60 //============================================================================
  61 #define FSR_DEBUG(arg)
  62 // Uncomment the following code to have debug output to cout
  63 //==========================================================
  64 //static Standard_Boolean mydebug = Standard_True;
  65 //#undef FSR_DEBUG
  66 //#define FSR_DEBUG(arg) {if (mydebug) { std::cout << arg << std::endl; }}
  67 //===========================================================
  68
  69 class MyDirFunction : public math_Function
  70 {
  71
  72   math_Vector *P0;
  73   math_Vector *Dir;
  74   math_Vector *P;
  75   math_Vector *FV;
  76   math_FunctionSetWithDerivatives *F;
  77
  78 public :
  79
  80   MyDirFunction(math_Vector& V1,
  81     math_Vector& V2,
  82     math_Vector& V3,
  83     math_Vector& V4,
  84     math_FunctionSetWithDerivatives& f) ;
  85
  86   void Initialize(const math_Vector& p0, const math_Vector& dir) const;
  87   //For hp :
  88   Standard_Boolean Value(const math_Vector& Sol, math_Vector& FF,
  89     math_Matrix& DF, math_Vector& GH,
  90     Standard_Real& F2, Standard_Real& Gnr1);
  91   //     Standard_Boolean MyDirFunction::Value(const math_Vector& Sol, math_Vector& FF,
  92   //                                       math_Matrix& DF, math_Vector& GH,
  93   //                                       Standard_Real& F2, Standard_Real& Gnr1);
  94   Standard_Boolean Value(const Standard_Real x, Standard_Real& fval) ;
  95
  96 };
  97
  98 MyDirFunction::MyDirFunction(math_Vector& V1,
  99                              math_Vector& V2,
 100                              math_Vector& V3,
 101                              math_Vector& V4,
 102                              math_FunctionSetWithDerivatives& f) {
 103
 104                                P0  = &V1;
 105                                Dir = &V2;
 106                                P   = &V3;
 107                                FV =  &V4;
 108                                F   = &f;
 109 }
 110
 111 void MyDirFunction::Initialize(const math_Vector& p0,
 112                                const math_Vector& dir)  const
 113 {
 114   *P0 = p0;
 115   *Dir = dir;
 116 }
 117
 118
 119 Standard_Boolean MyDirFunction::Value(const Standard_Real x,
 120                                       Standard_Real& fval)
 121 {
 122   Standard_Real p;
 123   for(Standard_Integer i = P->Lower(); i <= P->Upper(); i++) {
 124     p = Dir->Value(i);
 125     P->Value(i) = p * x + P0->Value(i);
 126   }
 127   if( F->Value(*P, *FV) )
 128   {
 129
 130     Standard_Real aVal = 0.0;
 131
 132     for(Standard_Integer i = FV->Lower(); i <= FV->Upper(); i++)
 133     {
 134       aVal = FV->Value(i);
 135       if(aVal <= -1.e+100) // Precision::HalfInfinite() later
 136         return Standard_False;
 137       else if(aVal >= 1.e+100) // Precision::HalfInfinite() later
 138         return Standard_False;
 139     }
 140
 141     fval = 0.5 * (FV->Norm2());
 142     return Standard_True;
 143   }
 144   return Standard_False;
 145 }
 146
 147 Standard_Boolean  MyDirFunction::Value(const math_Vector& Sol,
 148                                        math_Vector& FF,
 149                                        math_Matrix& DF,
 150                                        math_Vector& GH,
 151                                        Standard_Real& F2,
 152                                        Standard_Real& Gnr1)
 153 {
 154   if( F->Values(Sol, FF, DF) ) {
 155
 156     Standard_Real aVal = 0.;
 157
 158     for(Standard_Integer i = FF.Lower(); i <= FF.Upper(); i++) {
 159       // modified by NIZHNY-MKK  Mon Oct  3 17:56:50 2005.BEGIN
 160       aVal = FF.Value(i);
 161       if(aVal < 0.) {
 162         if(aVal <= -1.e+100) // Precision::HalfInfinite() later
 163           //       if(Precision::IsInfinite(Abs(FF.Value(i)))) {
 164           //    F2 = Precision::Infinite();
 165           //    Gnr1 = Precision::Infinite();
 166           return Standard_False;
 167       }
 168       else if(aVal >= 1.e+100) // Precision::HalfInfinite() later
 169         return Standard_False;
 170       // modified by NIZHNY-MKK  Mon Oct  3 17:57:05 2005.END
 171     }
 172
 173
 174     F2 = 0.5 * (FF.Norm2());
 175     GH.TMultiply(DF, FF);
 176     for(Standard_Integer i = GH.Lower(); i <= GH.Upper(); i++)
 177     {
 178       if(Precision::IsInfinite((GH.Value(i))))
 179       {
 180         return Standard_False;
 181       }
 182     }
 183     Gnr1 = GH.Norm2();
 184     return Standard_True;
 185   }
 186   return Standard_False;
 187 }
 188
 189
 190 //--------------------------------------------------------------
 191 static Standard_Boolean MinimizeDirection(const math_Vector&   P0,
 192                                           const math_Vector&   P1,
 193                                           const math_Vector&   P2,
 194                                           const Standard_Real  F1,
 195                                           math_Vector&         Delta,
 196                                           const math_Vector&   Tol,
 197                                           MyDirFunction& F)
 198                                           // Purpose : minimisation a partir de 3 points
 199                                           //-------------------------------------------------------
 200 {
 201   // (1) Evaluation d'un tolerance parametrique 1D
 202   Standard_Real tol1d = 2.1 , invnorme, tsol;
 203   Standard_Real Eps = 1.e-16;
 204   Standard_Real ax, bx, cx;
 205
 206   for (Standard_Integer ii =1; ii<=Tol.Length(); ii++) {
 207     invnorme = Abs(Delta(ii));
 208     if  (invnorme>Eps) tol1d = Min(tol1d, Tol(ii) / invnorme);
 209   }
 210   if (tol1d > 1.9) return Standard_False; //Pas la peine de se fatiguer
 211   tol1d /= 3;
 212
 213   //JR/Hp :
 214   const math_Vector& PP0 = P0 ;
 215   const math_Vector& PP1 = P1 ;
 216   Delta = PP1 - PP0;
 217   //  Delta = P1 - P0;
 218   invnorme = Delta.Norm();
 219   if (invnorme <= Eps) return Standard_False;
 220   invnorme = ((Standard_Real) 1) / invnorme;
 221
 222   F.Initialize(P1, Delta);
 223
 224   // (2) On minimise
 225   FSR_DEBUG ("      minimisation dans la direction")
 226     ax = -1; bx = 0;
 227   cx = (P2-P1).Norm()*invnorme;
 228   if (cx < 1.e-2)
 229     return Standard_False;
 230
 231   math_BrentMinimum Sol(tol1d, 100, tol1d);
 232   Sol.Perform(F, ax, bx, cx);
 233
 234   if(Sol.IsDone()) {
 235     tsol = Sol.Location();
 236     if (Sol.Minimum() < F1) {
 237       Delta.Multiply(tsol);
 238       return Standard_True;
 239     }
 240   }
 241   return Standard_False;
 242 }
 243
 244 //----------------------------------------------------------------------
 245 static Standard_Boolean MinimizeDirection(const math_Vector&   P,
 246                                           math_Vector&   Dir,
 247                                           const Standard_Real& PValue,
 248                                           const Standard_Real& PDirValue,
 249                                           const math_Vector&   Gradient,
 250                                           const math_Vector&   DGradient,
 251                                           const math_Vector&   Tol,
 252                                           MyDirFunction& F)
 253                                           // Purpose: minimisation a partir de 2 points et une derives
 254                                           //----------------------------------------------------------------------
 255
 256 {
 257   if(Precision::IsInfinite(PValue) || Precision::IsInfinite(PDirValue))
 258   {
 259     return Standard_False;
 260   }
 261   // (0) Evaluation d'un tolerance parametrique 1D
 262   Standard_Boolean good = Standard_False;
 263   Standard_Real Eps = 1.e-20;
 264   Standard_Real tol1d = 1.1, Result = PValue, absdir;
 265
 266   for (Standard_Integer ii =1; ii<=Tol.Length(); ii++) {
 267     absdir = Abs(Dir(ii));
 268     if (absdir >Eps) tol1d = Min(tol1d, Tol(ii) / absdir);
 269   }
 270   if (tol1d > 0.9) return Standard_False;
 271
 272   // (1) On realise une premiere interpolation quadratique
 273   Standard_Real ax, bx, cx, df1, df2, Delta, tsol, fsol, tsolbis;
 274   FSR_DEBUG("     essai d interpolation");
 275
 276   df1 = Gradient*Dir;
 277   df2 = DGradient*Dir;
 278
 279   if (df1<-Eps && df2>Eps) { // cuvette
 280     tsol = - df1 / (df2 - df1);
 281   }
 282   else {
 283     cx = PValue;
 284     bx = df1;
 285     ax = PDirValue - (bx+cx);
 286
 287     if (Abs(ax) <= Eps) { // cas lineaire
 288       if ((Abs(bx) >= Eps)) tsol = - cx/bx;
 289       else                  tsol = 0;
 290     }
 291     else { // cas quadratique
 292       Delta = bx*bx - 4*ax*cx;
 293       if (Delta > 1.e-9) {
 294         // il y a des racines, on prend la plus proche de 0
 295         Delta = Sqrt(Delta);
 296         tsol = -(bx + Delta);
 297         tsolbis = (Delta - bx);
 298         if (Abs(tsolbis) < Abs(tsol)) tsol = tsolbis;
 299         tsol /= 2*ax;
 300       }
 301       else {
 302         // pas ou peu de racine : on "extremise"
 303         tsol = -(0.5*bx)/ax;
 304       }
 305     }
 306   }
 307
 308   if (Abs(tsol) >= 1) return Standard_False; //resultat sans interet
 309
 310   F.Initialize(P, Dir);
 311   F.Value(tsol, fsol);
 312
 313   if (fsol<PValue) {
 314     good = Standard_True;
 315     Result = fsol;
 316     FSR_DEBUG("t= "<<tsol<<" F = " << fsol << " OldF = "<<PValue);
 317   }
 318
 319   // (2) Si l'on a pas assez progresser on realise une recherche
 320   //     en bonne et due forme, a partir des inits precedents
 321   if ((fsol > 0.2*PValue) && (tol1d < 0.5))
 322   {
 323
 324     if (tsol <0) {
 325       ax = tsol; bx = 0.0; cx = 1.0;
 326     }
 327     else {
 328       ax = 0.0; bx = tsol; cx = 1.0;
 329     }
 330     FSR_DEBUG(" minimisation dans la direction");
 331
 332     math_BrentMinimum Sol(tol1d, 100, tol1d);
 333
 334     // Base invocation.
 335     Sol.Perform(F, ax, bx, cx);
 336     if(Sol.IsDone())
 337     {
 338       if (Sol.Minimum() <= Result)
 339       {
 340         tsol = Sol.Location();
 341         good = Standard_True;
 342         Result = Sol.Minimum();
 343
 344         // Objective function changes too fast ->
 345         // perform additional computations.
 346         if (Gradient.Norm2() > 1.0 / Precision::SquareConfusion() &&
 347             tsol > ax &&
 348             tsol < cx) // Solution inside of (ax, cx) interval.
 349         {
 350           // First and second part invocation.
 351           Sol.Perform(F, ax, (ax + tsol) / 2.0, tsol);
 352           if(Sol.IsDone())
 353           {
 354             if (Sol.Minimum() <= Result)
 355             {
 356               tsol = Sol.Location();
 357               good = Standard_True;
 358               Result = Sol.Minimum();
 359             }
 360           }
 361
 362           Sol.Perform(F, tsol, (cx + tsol) / 2.0, cx);
 363           if(Sol.IsDone())
 364           {
 365             if (Sol.Minimum() <= Result)
 366             {
 367               tsol = Sol.Location();
 368               good = Standard_True;
 369               Result = Sol.Minimum();
 370             }
 371           }
 372         }
 373       } // if (Sol.Minimum() <= Result)
 374     } // if(Sol.IsDone())
 375   }
 376
 377   if (good)
 378   {
 379     // mise a jour du Delta
 380     Dir.Multiply(tsol);
 381   }
 382   return good;
 383 }
 384
 385 //------------------------------------------------------
 386 static void SearchDirection(const math_Matrix& DF,
 387                             const math_Vector& GH,
 388                             const math_Vector& FF,
 389                             Standard_Boolean ChangeDirection,
 390                             const math_Vector& InvLengthMax,
 391                             math_Vector& Direction,
 392                             Standard_Real& Dy)
 393
 394 {
 395   Standard_Integer Ninc = DF.ColNumber(), Neq = DF.RowNumber();
 396   Standard_Real Eps = 1.e-32;
 397   if (!ChangeDirection) {
 398     if (Ninc == Neq) {
 399       for (Standard_Integer i = FF.Lower(); i <= FF.Upper(); i++) {
 400         Direction(i) = -FF(i);
 401       }
 402       math_Gauss Solut(DF, 1.e-9);
 403       if (Solut.IsDone()) Solut.Solve(Direction);
 404       else { // we have to "forget" singular directions.
 405         FSR_DEBUG(" Matrice singuliere : On prend SVD");
 406         math_SVD SolvebySVD(DF);
 407         if (SolvebySVD.IsDone()) SolvebySVD.Solve(-1*FF, Direction);
 408         else ChangeDirection = Standard_True;
 409       }
 410     }
 411     else if (Ninc > Neq) {
 412       math_SVD Solut(DF);
 413       if (Solut.IsDone()) Solut.Solve(-1*FF, Direction);
 414       else ChangeDirection = Standard_True;
 415     }
 416     else if (Ninc < Neq) {         // Calcul par GaussLeastSquare
 417       math_GaussLeastSquare Solut(DF);
 418       if (Solut.IsDone()) Solut.Solve(-1*FF, Direction);
 419       else ChangeDirection = Standard_True;
 420     }
 421   }
 422   // Il vaut mieux interdire des directions trops longue
 423   // Afin de blinder les cas trop mal conditionne
 424   // PMN 12/05/97 Traitement des singularite dans les conges
 425   // Sur des surfaces periodiques
 426
 427   Standard_Real ratio = Abs( Direction(Direction.Lower())
 428     *InvLengthMax(Direction.Lower()) );
 429   Standard_Integer i;
 430   for (i = Direction.Lower()+1; i<=Direction.Upper(); i++) {
 431     ratio = Max(ratio,  Abs( Direction(i)*InvLengthMax(i)) );
 432   }
 433   if (ratio > 1) {
 434     Direction /= ratio;
 435   }
 436
 437   Dy = Direction*GH;
 438   if (Dy >= -Eps) { // newton "ne descend pas" on prend le gradient
 439     ChangeDirection = Standard_True;
 440   }
 441   if (ChangeDirection) { // On va faire un gradient !
 442     for (i = Direction.Lower(); i <= Direction.Upper(); i++) {
 443       Direction(i) = - GH(i);
 444     }
 445     Dy = - (GH.Norm2());
 446   }
 447 }
 448
 449
 450 //=====================================================================
 451 static void SearchDirection(const math_Matrix& DF,
 452                             const math_Vector& GH,
 453                             const math_Vector& FF,
 454                             const math_IntegerVector& Constraints,
 455                             //                      const math_Vector& X, // Le point d'init
 456                             const math_Vector& , // Le point d'init
 457                             Standard_Boolean ChangeDirection,
 458                             const math_Vector& InvLengthMax,
 459                             math_Vector& Direction,
 460                             Standard_Real& Dy)
 461                             //Purpose : Recherche une direction (et un pas si Newton Fonctionne) le long
 462                             //          d'une frontiere
 463                             //=====================================================================
 464 {
 465   Standard_Integer Ninc = DF.ColNumber(), Neq = DF.RowNumber();
 466   Standard_Integer i, j, k, Cons = 0;
 467
 468   // verification sur les bornes imposees:
 469
 470   for (i = 1; i <= Ninc; i++) {
 471     if (Constraints(i) != 0) Cons++;
 472     // sinon le systeme a resoudre ne change pas.
 473   }
 474
 475   if (Cons == 0) {
 476     SearchDirection(DF, GH, FF, ChangeDirection, InvLengthMax,
 477       Direction, Dy);
 478   }
 479   else if (Cons == Ninc) { // il n'y a plus rien a faire...
 480     for(i = Direction.Lower(); i <= Direction.Upper(); i++) {
 481       Direction(i) = 0;
 482     }
 483     Dy = 0;
 484   }
 485   else { //(1) Cas general : On definit un sous probleme
 486     math_Matrix DF2(1, Neq, 1, Ninc-Cons);
 487     math_Vector MyGH(1, Ninc-Cons);
 488     math_Vector MyDirection(1, Ninc-Cons);
 489     math_Vector MyInvLengthMax(1, Ninc);
 490
 491     for (k=1, i = 1; i <= Ninc; i++) {
 492       if (Constraints(i) == 0) {
 493         MyGH(k) = GH(i);
 494         MyInvLengthMax(k) = InvLengthMax(i);
 495         MyDirection(k) = Direction(i);
 496         for (j = 1; j <= Neq; j++) {
 497           DF2(j, k) = DF(j, i);
 498         }
 499         k++; //on passe a l'inconnue suivante
 500       }
 501     }
 502     //(2) On le resoud
 503     SearchDirection(DF2, MyGH, FF, ChangeDirection, MyInvLengthMax,
 504       MyDirection, Dy);
 505
 506     // (3) On l'interprete...
 507     // Reconstruction de Direction:
 508     for (i = 1, k = 1; i <= Ninc; i++) {
 509       if (Constraints(i) == 0) {
 510         if (!ChangeDirection) {
 511           Direction(i) = MyDirection(k);
 512         }
 513         else Direction(i) = - GH(i);
 514         k++;
 515       }
 516       else {
 517         Direction(i) = 0.0;
 518       }
 519     }
 520   }
 521 }
 522
 523
 524
 525 //====================================================
 526 Standard_Boolean Bounds(const math_Vector&  InfBound,
 527                         const math_Vector&  SupBound,
 528                         const math_Vector&  Tol,
 529                         math_Vector&        Sol,
 530                         const math_Vector&  SolSave,
 531                         math_IntegerVector& Constraints,
 532                         math_Vector&        Delta,
 533                         Standard_Boolean&   theIsNewSol)
 534                         //
 535                         // Purpose: Trims an initial solution Sol to be within a domain defined by
 536                         //   InfBound and SupBound. Delta will contain a distance between final Sol and
 537                         //   SolSave.
 538                         //   IsNewSol returns False, if final Sol fully coincides with SolSave, i.e.
 539                         //   if SolSave already lied on a boundary and initial Sol was fully beyond it
 540                         //======================================================
 541 {
 542   Standard_Boolean Out = Standard_False;
 543   Standard_Integer i, Ninc = Sol.Length();
 544   Standard_Real    monratio = 1;
 545
 546   theIsNewSol = Standard_True;
 547
 548   // Calcul du ratio de recadrage
 549   for (i = 1; i <= Ninc; i++) {
 550     Constraints(i) = 0;
 551     Delta(i) = Sol(i) - SolSave(i);
 552     if (InfBound(i) == SupBound(i)) {
 553       Constraints(i) = 1;
 554       Out = Standard_True; // Ok mais, cela devrait etre eviter
 555     }
 556     else if(Sol(i) < InfBound(i)) {
 557       Constraints(i) = 1;
 558       Out = Standard_True;
 559       // Delta(i) is negative
 560       if (-Delta(i) > Tol(i)) // Afin d'eviter des ratio nulles pour rien
 561         monratio = Min(monratio, (InfBound(i) - SolSave(i))/Delta(i) );
 562     }
 563     else if (Sol(i) > SupBound(i)) {
 564       Constraints(i) = 1;
 565       Out = Standard_True;
 566       // Delta(i) is positive
 567       if (Delta(i) > Tol(i))
 568         monratio = Min(monratio, (SupBound(i) - SolSave(i))/Delta(i) );
 569     }
 570   }
 571
 572   if (Out){ // Troncature et derniers recadrage pour blinder (pb numeriques)
 573     if (monratio == 0.0) {
 574       theIsNewSol = Standard_False;
 575       Sol = SolSave;
 576       Delta.Init (0.0);
 577     } else {
 578       Delta *= monratio;
 579       Sol = SolSave+Delta;
 580       for (i = 1; i <= Ninc; i++) {
 581         if(Sol(i) < InfBound(i))  {
 582           Sol(i) = InfBound(i);
 583           Delta(i) = Sol(i) - SolSave(i);
 584         }
 585         else if (Sol(i) > SupBound(i)) {
 586           Sol(i) = SupBound(i);
 587           Delta(i) = Sol(i) - SolSave(i);
 588         }
 589       }
 590     }
 591   }
 592   return Out;
 593 }
 594
 595
 596
 597
 598 //=======================================================================
 599 //function : math_FunctionSetRoot
 600 //purpose  : Constructor
 601 //=======================================================================
 602 math_FunctionSetRoot::math_FunctionSetRoot(
 603   math_FunctionSetWithDerivatives& theFunction,
 604   const math_Vector&               theTolerance,
 605   const Standard_Integer           theNbIterations)
 606
 607 : Delta(1, theFunction.NbVariables()),
 608   Sol  (1, theFunction.NbVariables()),
 609   DF   (1, theFunction.NbEquations() , 1, theFunction.NbVariables()),
 610   Tol  (1, theFunction.NbVariables()),
 611   Done    (Standard_False),
 612   Kount   (0),
 613   State   (0),
 614   Itermax (theNbIterations),
 615   InfBound(1, theFunction.NbVariables(), RealFirst()),
 616   SupBound(1, theFunction.NbVariables(), RealLast ()),
 617   SolSave (1, theFunction.NbVariables()),
 618   GH      (1, theFunction.NbVariables()),
 619   DH      (1, theFunction.NbVariables()),
 620   DHSave  (1, theFunction.NbVariables()),
 621   FF      (1, theFunction.NbEquations()),
 622   PreviousSolution(1, theFunction.NbVariables()),
 623   Save    (0, theNbIterations),
 624   Constraints(1, theFunction.NbVariables()),
 625   Temp1   (1, theFunction.NbVariables()),
 626   Temp2   (1, theFunction.NbVariables()),
 627   Temp3   (1, theFunction.NbVariables()),
 628   Temp4   (1, theFunction.NbEquations()),
 629   myIsDivergent(Standard_False)
 630 {
 631   SetTolerance(theTolerance);
 632 }
 633
 634 //=======================================================================
 635 //function : math_FunctionSetRoot
 636 //purpose  : Constructor
 637 //=======================================================================
 638 math_FunctionSetRoot::math_FunctionSetRoot(math_FunctionSetWithDerivatives& theFunction,
 639                                            const Standard_Integer           theNbIterations)
 640
 641 : Delta(1, theFunction.NbVariables()),
 642   Sol  (1, theFunction.NbVariables()),
 643   DF   (1, theFunction.NbEquations() , 1, theFunction.NbVariables()),
 644   Tol  (1, theFunction.NbVariables()),
 645   Done    (Standard_False),
 646   Kount   (0),
 647   State   (0),
 648   Itermax (theNbIterations),
 649   InfBound(1, theFunction.NbVariables(), RealFirst()),
 650   SupBound(1, theFunction.NbVariables(), RealLast ()),
 651   SolSave (1, theFunction.NbVariables()),
 652   GH      (1, theFunction.NbVariables()),
 653   DH      (1, theFunction.NbVariables()),
 654   DHSave  (1, theFunction.NbVariables()),
 655   FF      (1, theFunction.NbEquations()),
 656   PreviousSolution(1, theFunction.NbVariables()),
 657   Save    (0, theNbIterations),
 658   Constraints(1, theFunction.NbVariables()),
 659   Temp1   (1, theFunction.NbVariables()),
 660   Temp2   (1, theFunction.NbVariables()),
 661   Temp3   (1, theFunction.NbVariables()),
 662   Temp4   (1, theFunction.NbEquations()),
 663   myIsDivergent(Standard_False)
 664 {
 665 }
 666
 667 //=======================================================================
 668 //function : ~math_FunctionSetRoot
 669 //purpose  : Destructor
 670 //=======================================================================
 671 math_FunctionSetRoot::~math_FunctionSetRoot()
 672 {
 673 }
 674
 675 //=======================================================================
 676 //function : SetTolerance
 677 //purpose  :
 678 //=======================================================================
 679 void math_FunctionSetRoot::SetTolerance(const math_Vector& theTolerance)
 680 {
 681   for (Standard_Integer i = 1; i <= Tol.Length(); ++i)
 682     Tol(i) = theTolerance(i);
 683 }
 684
 685 //=======================================================================
 686 //function : Perform
 687 //purpose  :
 688 //=======================================================================
 689 void math_FunctionSetRoot::Perform(math_FunctionSetWithDerivatives& theFunction,
 690                                    const math_Vector&               theStartingPoint,
 691                                    const Standard_Boolean           theStopOnDivergent)
 692 {
 693   Perform(theFunction, theStartingPoint, InfBound, SupBound, theStopOnDivergent);
 694 }
 695
 696 //=======================================================================
 697 //function : Perform
 698 //purpose  :
 699 //=======================================================================
 700 void math_FunctionSetRoot::Perform(math_FunctionSetWithDerivatives& F,
 701                                    const math_Vector& StartingPoint,
 702                                    const math_Vector& theInfBound,
 703                                    const math_Vector& theSupBound,
 704                                    Standard_Boolean theStopOnDivergent)
 705 {
 706   Standard_Integer Ninc = F.NbVariables(), Neq = F.NbEquations();
 707
 708   if ((Neq <= 0)                      ||
 709     (StartingPoint.Length()!= Ninc) ||
 710     (theInfBound.Length() != Ninc)     ||
 711     (theSupBound.Length() != Ninc))  { throw Standard_DimensionError(); }
 712
 713   Standard_Integer i;
 714   Standard_Boolean ChangeDirection = Standard_False, Sort = Standard_False, isNewSol = Standard_False;
 715   Standard_Boolean Good, Verif;
 716   Standard_Boolean Stop;
 717   const Standard_Real EpsSqrt = 1.e-16, Eps = 1.e-32, Eps2 = 1.e-64, Progres = 0.005;
 718   Standard_Real F2, PreviousMinimum, Dy, OldF;
 719   Standard_Real Ambda, Ambda2, Gnr1, Oldgr;
 720   math_Vector InvLengthMax(1, Ninc); // Pour bloquer les pas a 1/4 du domaine
 721   math_IntegerVector aConstraints(1, Ninc); // Pour savoir sur quels bord on se trouve
 722   for (i = 1; i <= Ninc ; i++) {
 723     const Standard_Real aSupBound  = Min (theSupBound(i),  Precision::Infinite());
 724     const Standard_Real anInfBound = Max (theInfBound(i), -Precision::Infinite());
 725     InvLengthMax(i) = 1. / Max((aSupBound - anInfBound)/4, 1.e-9);
 726   }
 727
 728   MyDirFunction F_Dir(Temp1, Temp2, Temp3, Temp4, F);
 729   Standard_Integer  DescenteIter;
 730
 731   Done = Standard_False;
 732   Sol = StartingPoint;
 733   Kount = 0;
 734
 735   //
 736   myIsDivergent = Standard_False;
 737   for (i = 1; i <= Ninc; i++)
 738   {
 739     myIsDivergent = myIsDivergent
 740                  || Sol(i) < theInfBound(i)
 741                  || Sol(i) > theSupBound(i);
 742   }
 743   if (theStopOnDivergent && myIsDivergent)
 744   {
 745     return;
 746   }
 747
 748   // Verification de la validite des inconnues par rapport aux bornes.
 749   // Recentrage sur les bornes si pas valide.
 750   for ( i = 1; i <= Ninc; i++) {
 751     if (Sol(i) <= theInfBound(i)) Sol(i) = theInfBound(i);
 752     else if (Sol(i) > theSupBound(i)) Sol(i) = theSupBound(i);
 753   }
 754
 755   // Calcul de la premiere valeur de F et de son gradient
 756   if(!F_Dir.Value(Sol, FF, DF, GH, F2, Gnr1)) {
 757     Done = Standard_False;
 758     if (!theStopOnDivergent || !myIsDivergent)
 759     {
 760       State = F.GetStateNumber();
 761     }
 762     return;
 763   }
 764   Ambda2 = Gnr1;
 765   // Le rang 0 de Save ne doit servir q'au test accelarteur en fin de boucle
 766   // s'il on est dejas sur la solution, il faut leurer ce test pour eviter
 767   // de faire une seconde iteration...
 768   Save(0) = Max (F2, EpsSqrt);
 769   Standard_Real aTol_Func = Epsilon(F2);
 770   FSR_DEBUG("=== Mode Debug de Function Set Root" << std::endl);
 771   FSR_DEBUG("    F2 Initial = " << F2);
 772
 773   if ((F2 <= Eps) || (Gnr1 <= Eps2)) {
 774     Done = Standard_False;
 775     if (!theStopOnDivergent || !myIsDivergent)
 776     {
 777       Done = Standard_True;
 778       State = F.GetStateNumber();
 779     }
 780     return;
 781   }
 782
 783   for (Kount = 1; Kount <= Itermax; Kount++) {
 784     PreviousMinimum = F2;
 785     Oldgr = Gnr1;
 786     PreviousSolution = Sol;
 787     SolSave = Sol;
 788
 789     SearchDirection(DF, GH, FF, ChangeDirection, InvLengthMax, DH, Dy);
 790     if (Abs(Dy) <= Eps) {
 791       Done = Standard_False;
 792       if (!theStopOnDivergent || !myIsDivergent)
 793       {
 794         Done = Standard_True;
 795         ////modified by jgv, 31.08.2011////
 796         F.Value(Sol, FF); //update F before GetStateNumber
 797         ///////////////////////////////////
 798         State = F.GetStateNumber();
 799       }
 800       return;
 801     }
 802     if (ChangeDirection) {
 803       Ambda = Ambda2 / Sqrt(Abs(Dy));
 804       if (Ambda > 1.0) Ambda = 1.0;
 805     }
 806     else {
 807       Ambda = 1.0;
 808       Ambda2 = 0.5*Ambda/DH.Norm();
 809     }
 810
 811     for( i = Sol.Lower(); i <= Sol.Upper(); i++) {
 812       Sol(i) = Sol(i) + Ambda * DH(i);
 813     }
 814     //
 815     for (i = 1; i <= Ninc; i++)
 816     {
 817       myIsDivergent = myIsDivergent
 818                    || Sol(i) < theInfBound(i)
 819                    || Sol(i) > theSupBound(i);
 820     }
 821     if (theStopOnDivergent && myIsDivergent)
 822     {
 823       return;
 824     }
 825     //
 826     Sort = Bounds(theInfBound, theSupBound, Tol, Sol, SolSave,
 827       aConstraints, Delta, isNewSol);
 828
 829
 830     DHSave = GH;
 831     if (isNewSol) {
 832       //      F_Dir.Value(Sol, FF, DF, GH, F2, Gnr1);
 833       if(!F_Dir.Value(Sol, FF, DF, GH, F2, Gnr1)) {
 834         Done = Standard_False;
 835         if (!theStopOnDivergent || !myIsDivergent)
 836         {
 837           State = F.GetStateNumber();
 838         }
 839         return;
 840       }
 841     }
 842
 843     FSR_DEBUG("Kount         = " << Kount);
 844     FSR_DEBUG("Le premier F2 = " << F2);
 845     FSR_DEBUG("Direction     = " << ChangeDirection);
 846
 847     if ((F2 <= Eps) || (Gnr1 <= Eps2)) {
 848       Done = Standard_False;
 849       if (!theStopOnDivergent || !myIsDivergent)
 850       {
 851         Done = Standard_True;
 852         ////modified by jgv, 31.08.2011////
 853         F.Value(Sol, FF); //update F before GetStateNumber
 854         ///////////////////////////////////
 855         State = F.GetStateNumber();
 856       }
 857       return;
 858     }
 859
 860     if (Sort || (F2/PreviousMinimum > Progres)) {
 861       Dy = GH*DH;
 862       OldF = PreviousMinimum;
 863       Stop = Standard_False;
 864       Good = Standard_False;
 865       DescenteIter = 0;
 866       Standard_Boolean Sortbis;
 867
 868       // -------------------------------------------------
 869       // Traitement standard sans traitement des bords
 870       // -------------------------------------------------
 871       if (!Sort) { // si l'on est pas sortie on essaye de progresser en avant
 872         while((F2/PreviousMinimum > Progres) && !Stop) {
 873           if (F2 < OldF && (Dy < 0.0)) {
 874             // On essaye de progresser dans cette direction.
 875             FSR_DEBUG(" iteration de descente = " << DescenteIter);
 876             DescenteIter++;
 877             SolSave = Sol;
 878             OldF = F2;
 879             for( i = Sol.Lower(); i <= Sol.Upper(); i++) {
 880               Sol(i) = Sol(i) + Ambda * DH(i);
 881             }
 882             //
 883             for (i = 1; i <= Ninc; i++)
 884             {
 885               myIsDivergent = myIsDivergent
 886                            || Sol(i) < theInfBound(i)
 887                            || Sol(i) > theSupBound(i);
 888             }
 889             if (theStopOnDivergent && myIsDivergent)
 890             {
 891               return;
 892             }
 893             //
 894             Stop = Bounds(theInfBound, theSupBound, Tol, Sol, SolSave,
 895               aConstraints, Delta, isNewSol);
 896             FSR_DEBUG(" Augmentation de lambda");
 897             Ambda *= 1.7;
 898           }
 899           else {
 900             if ((F2 >= OldF)||(F2 >= PreviousMinimum)) {
 901               Good = Standard_False;
 902               if (DescenteIter == 0) {
 903                 // C'est le premier pas qui flanche, on fait une interpolation.
 904                 // et on minimise si necessaire.
 905                 DescenteIter++;
 906                 Good = MinimizeDirection(SolSave, Delta, OldF, F2, DHSave, GH,
 907                   Tol, F_Dir);
 908               }
 909               else if (ChangeDirection || (DescenteIter>1)
 910                 || (OldF>PreviousMinimum) ) {
 911                   // La progression a ete utile, on minimise...
 912                   DescenteIter++;
 913                   Good = MinimizeDirection(PreviousSolution, SolSave, Sol,
 914                     OldF, Delta, Tol, F_Dir);
 915               }
 916               if (!Good) {
 917                 Sol = SolSave;
 918                 F2 = OldF;
 919               }
 920               else {
 921                 Sol = SolSave+Delta;
 922                 //
 923                 for (i = 1; i <= Ninc; i++)
 924                 {
 925                   myIsDivergent = myIsDivergent
 926                                || Sol(i) < theInfBound(i)
 927                                || Sol(i) > theSupBound(i);
 928                 }
 929                 if (theStopOnDivergent && myIsDivergent)
 930                 {
 931                   return;
 932                 }
 933                 //
 934                 Sort = Bounds(theInfBound, theSupBound, Tol, Sol, SolSave,
 935                   aConstraints, Delta, isNewSol);
 936               }
 937               Sort = Standard_False; // On a rejete le point sur la frontiere
 938             }
 939             Stop = Standard_True; // et on sort dans tous les cas...
 940           }
 941           DHSave = GH;
 942           if (isNewSol) {
 943             //            F_Dir.Value(Sol, FF, DF, GH, F2, Gnr1);
 944             if(!F_Dir.Value(Sol, FF, DF, GH, F2, Gnr1)) {
 945               Done = Standard_False;
 946               if (!theStopOnDivergent || !myIsDivergent)
 947               {
 948                 State = F.GetStateNumber();
 949               }
 950               return;
 951             }
 952           }
 953           Dy = GH*DH;
 954           if (Abs(Dy) <= Eps) {
 955             if (F2 > OldF)
 956               Sol = SolSave;
 957             Done = Standard_False;
 958             if (!theStopOnDivergent || !myIsDivergent)
 959             {
 960               Done = Standard_True;
 961               ////modified by jgv, 31.08.2011////
 962               F.Value(Sol, FF); //update F before GetStateNumber
 963               ///////////////////////////////////
 964               State = F.GetStateNumber();
 965             }
 966             return;
 967           }
 968           if (DescenteIter >= 100) {
 969             Stop = Standard_True;
 970           }
 971         }
 972         FSR_DEBUG("--- Sortie du Traitement Standard");
 973         FSR_DEBUG("    DescenteIter = "<<DescenteIter << " F2 = " << F2);
 974       }
 975       // ------------------------------------
 976       //  on passe au traitement des bords
 977       // ------------------------------------
 978       if (Sort) {
 979         Stop = (F2 > 1.001*OldF); // Pour ne pas progresser sur le bord
 980         Sortbis = Sort;
 981         DescenteIter = 0;
 982         while (Sortbis && ((F2<OldF)|| (DescenteIter == 0))
 983           && (!Stop)) {
 984             DescenteIter++;
 985             // On essaye de progresser sur le bord
 986             SolSave = Sol;
 987             OldF = F2;
 988             SearchDirection(DF, GH, FF,  aConstraints, Sol,
 989               ChangeDirection, InvLengthMax, DH, Dy);
 990             FSR_DEBUG(" Conditional Direction = " << ChangeDirection);
 991             if (Dy<-Eps) { //Pour eviter des calculs inutiles et des /0...
 992               if (ChangeDirection) {
 993
 994                 //            Ambda = Ambda2 / Sqrt(Abs(Dy));
 995                 Ambda = Ambda2 / Sqrt(-Dy);
 996                 if (Ambda > 1.0) Ambda = 1.0;
 997               }
 998               else {
 999                 Ambda = 1.0;
1000                 Ambda2 = 0.5*Ambda/DH.Norm();
1001               }
1002
1003               for( i = Sol.Lower(); i <= Sol.Upper(); i++) {
1004                 Sol(i) = Sol(i) + Ambda * DH(i);
1005               }
1006               //
1007               for (i = 1; i <= Ninc; i++)
1008               {
1009                 myIsDivergent = myIsDivergent
1010                              || Sol(i) < theInfBound(i)
1011                              || Sol(i) > theSupBound(i);
1012               }
1013               if (theStopOnDivergent && myIsDivergent)
1014               {
1015                 return;
1016               }
1017               //
1018               Sortbis = Bounds(theInfBound, theSupBound, Tol, Sol, SolSave,
1019                 aConstraints, Delta, isNewSol);
1020
1021               DHSave = GH;
1022               if (isNewSol) {
1023                 //              F_Dir.Value(Sol, FF, DF, GH, F2, Gnr1);
1024                 if(!F_Dir.Value(Sol, FF, DF, GH, F2, Gnr1)) {
1025                   Done = Standard_False;
1026                   if (!theStopOnDivergent || !myIsDivergent)
1027                   {
1028                     State = F.GetStateNumber();
1029                   }
1030                   return;
1031                 }
1032               }
1033               Ambda2 = Gnr1;
1034               FSR_DEBUG("---  Iteration au bords : " << DescenteIter);
1035               FSR_DEBUG("---  F2 = " << F2);
1036             }
1037             else {
1038               Stop = Standard_True;
1039             }
1040
1041             while((F2/PreviousMinimum > Progres) && (F2<OldF) && (!Stop) ) {
1042               DescenteIter++;
1043               FSR_DEBUG("--- Iteration de descente conditionnel = " << DescenteIter);
1044               if (F2 < OldF && Dy < 0.0) {
1045                 // On essaye de progresser dans cette direction.
1046                 SolSave = Sol;
1047                 OldF = F2;
1048                 for( i = Sol.Lower(); i <= Sol.Upper(); i++) {
1049                   Sol(i) = Sol(i) + Ambda * DH(i);
1050                 }
1051                 //
1052                 for (i = 1; i <= Ninc; i++)
1053                 {
1054                   myIsDivergent = myIsDivergent
1055                                || Sol(i) < theInfBound(i)
1056                                || Sol(i) > theSupBound(i);
1057                 }
1058                 if (theStopOnDivergent && myIsDivergent)
1059                 {
1060                   return;
1061                 }
1062                 //
1063                 Sortbis = Bounds(theInfBound, theSupBound, Tol, Sol, SolSave,
1064                   aConstraints, Delta, isNewSol);
1065               }
1066               DHSave = GH;
1067               if (isNewSol) {
1068                 //              F_Dir.Value(Sol, FF, DF, GH, F2, Gnr1);
1069                 if(!F_Dir.Value(Sol, FF, DF, GH, F2, Gnr1)) {
1070                   Done = Standard_False;
1071                   if (!theStopOnDivergent || !myIsDivergent)
1072                   {
1073                     State = F.GetStateNumber();
1074                   }
1075                   return;
1076                 }
1077               }
1078               Ambda2 = Gnr1;
1079               Dy = GH*DH;
1080               Stop = ((Dy >=0) || (DescenteIter >= 10) || Sortbis);
1081             }
1082             Stop = ((Dy >=0) || (DescenteIter >= 10));
1083         }
1084         if (((F2/PreviousMinimum > Progres) &&
1085           (F2>=OldF))||(F2>=PreviousMinimum)) {
1086             // On minimise par Brent
1087             DescenteIter++;
1088             Good = MinimizeDirection(SolSave, Delta, OldF, F2,
1089               DHSave, GH, Tol, F_Dir);
1090             if (!Good) {
1091               Sol = SolSave;
1092               Sort = Standard_False;
1093             }
1094             else {
1095               Sol = SolSave + Delta;
1096               //
1097               for (i = 1; i <= Ninc; i++)
1098               {
1099                 myIsDivergent = myIsDivergent
1100                              || Sol(i) < theInfBound(i)
1101                              || Sol(i) > theSupBound(i);
1102               }
1103               if (theStopOnDivergent && myIsDivergent)
1104               {
1105                 return;
1106               }
1107               //
1108               Sort = Bounds(theInfBound, theSupBound, Tol, Sol, SolSave,
1109                 aConstraints, Delta, isNewSol);
1110               if (isNewSol) {
1111                 //              F_Dir.Value(Sol, FF, DF, GH, F2, Gnr1);
1112                 if(!F_Dir.Value(Sol, FF, DF, GH, F2, Gnr1)) {
1113                   Done = Standard_False;
1114                   if (!theStopOnDivergent || !myIsDivergent)
1115                   {
1116                     State = F.GetStateNumber();
1117                   }
1118                   return;
1119                 }
1120               }
1121             }
1122             Dy = GH*DH;
1123         }
1124         FSR_DEBUG("--- Sortie du Traitement des Bords");
1125         FSR_DEBUG("--- DescenteIter = "<<DescenteIter << " F2 = " << F2);
1126       }
1127     }
1128
1129     // ---------------------------------------------
1130     //  on passe aux tests d'ARRET
1131     // ---------------------------------------------
1132     Save(Kount) = F2;
1133     // Est ce la solution ?
1134     if (ChangeDirection) Verif = Standard_True;
1135     // Gradient : Il faut eviter de boucler
1136     else {
1137       Verif = Standard_False;
1138       if (Kount > 1) { // Pour accelerer les cas quasi-quadratique
1139         if (Save(Kount-1)<1.e-4*Save(Kount-2)) Verif = Standard_True;
1140       }
1141       else Verif = (F2 < 1.e-6*Save(0)); //Pour les cas dejas solutions
1142     }
1143     if (Verif) {
1144       for(i = Delta.Lower(); i <= Delta.Upper(); i++) {
1145         Delta(i) = PreviousSolution(i) - Sol(i);
1146       }
1147
1148       if (IsSolutionReached(F)) {
1149         if (PreviousMinimum < F2) {
1150           Sol = SolSave;
1151         }
1152         Done = Standard_False;
1153         if (!theStopOnDivergent || !myIsDivergent)
1154         {
1155           Done = Standard_True;
1156           ////modified by jgv, 31.08.2011////
1157           F.Value(Sol, FF); //update F before GetStateNumber
1158           ///////////////////////////////////
1159           State = F.GetStateNumber();
1160         }
1161         return;
1162       }
1163     }
1164     //fin du test solution
1165
1166     // Analyse de la progression...
1167     //comparison of current minimum and previous minimum
1168     if ((F2 - PreviousMinimum) <= aTol_Func){
1169       if (Kount > 5) {
1170         // L'historique est il bon ?
1171         if (F2 >= 0.95*Save(Kount - 5)) {
1172           if (!ChangeDirection) ChangeDirection = Standard_True;
1173           else
1174           {
1175             Done = Standard_False;
1176             if (!theStopOnDivergent || !myIsDivergent)
1177             {
1178               Done = Standard_True;
1179               State = F.GetStateNumber();
1180             }
1181             return; //  si un gain inf a 5% on sort
1182           }
1183         }
1184         else ChangeDirection = Standard_False; //Si oui on recommence
1185       }
1186       else  ChangeDirection = Standard_False; //Pas d'historique on continue
1187       // Si le gradient ne diminue pas suffisemment par newton on essaie
1188       // le gradient sauf si f diminue (aussi bizarre que cela puisse
1189       // paraitre avec NEWTON le gradient de f peut augmenter alors que f
1190       // diminue: dans ce cas il faut garder NEWTON)
1191       if ((Gnr1 > 0.9*Oldgr) &&
1192         (F2 > 0.5*PreviousMinimum)) {
1193           ChangeDirection = Standard_True;
1194       }
1195
1196       // Si l'on ne decide pas de changer de strategie, on verifie,
1197       // si ce n'est dejas fait
1198       if ((!ChangeDirection) && (!Verif)) {
1199         for(i = Delta.Lower(); i <= Delta.Upper(); i++) {
1200           Delta(i) = PreviousSolution(i) - Sol(i);
1201         }
1202         if (IsSolutionReached(F)) {
1203           Done = Standard_False;
1204           if (!theStopOnDivergent || !myIsDivergent)
1205           {
1206             Done = Standard_True;
1207             ////modified by jgv, 31.08.2011////
1208             F.Value(Sol, FF); //update F before GetStateNumber
1209             ///////////////////////////////////
1210             State = F.GetStateNumber();
1211           }
1212           return;
1213         }
1214       }
1215     }
1216     else { // Cas de regression
1217       if (!ChangeDirection) { // On passe au gradient
1218         ChangeDirection = Standard_True;
1219         Sol = PreviousSolution;
1220         //      F_Dir.Value(Sol, FF, DF, GH, F2, Gnr1);
1221         if(!F_Dir.Value(Sol, FF, DF, GH, F2, Gnr1)) {
1222           Done = Standard_False;
1223           if (!theStopOnDivergent || !myIsDivergent)
1224           {
1225             State = F.GetStateNumber();
1226           }
1227           return;
1228         }
1229       }
1230       else
1231       {
1232
1233         if (!theStopOnDivergent || !myIsDivergent)
1234         {
1235           State = F.GetStateNumber();
1236         }
1237         return; // y a plus d'issues
1238       }
1239     }
1240   }
1241   if (!theStopOnDivergent || !myIsDivergent)
1242   {
1243     State = F.GetStateNumber();
1244   }
1245 }
1246
1247 //=======================================================================
1248 //function : Dump
1249 //purpose  :
1250 //=======================================================================
1251 void math_FunctionSetRoot::Dump(Standard_OStream& o) const
1252 {
1253   o << " math_FunctionSetRoot";
1254   if (Done) {
1255     o << " Status = Done\n";
1256     o << " Location value = " << Sol << "\n";
1257     o << " Number of iterations = " << Kount << "\n";
1258   }
1259   else {
1260     o << "Status = Not Done\n";
1261   }
1262 }
1263
1264 //=======================================================================
1265 //function : Root
1266 //purpose  :
1267 //=======================================================================
1268 void math_FunctionSetRoot::Root(math_Vector& Root) const
1269 {
1270   StdFail_NotDone_Raise_if(!Done, " ");
1271   Standard_DimensionError_Raise_if(Root.Length() != Sol.Length(), " ");
1272   Root = Sol;
1273 }
1274
1275 //=======================================================================
1276 //function : FunctionSetErrors
1277 //purpose  :
1278 //=======================================================================
1279 void math_FunctionSetRoot::FunctionSetErrors(math_Vector& Err) const
1280 {
1281   StdFail_NotDone_Raise_if(!Done, " ");
1282   Standard_DimensionError_Raise_if(Err.Length() != Sol.Length(), " ");
1283   Err = Delta;
1284 }