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