src/math/math_FunctionSetRoot.cxx

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