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