projects / occt.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
0023360: Test cases for command mkoffset produce different results on different versi...
[occt.git] / src / math / math_FunctionSetRoot.cxx
CommitLineData
b311480e 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
7fd59977 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.
7fd59977 48// - Recadrage au niveau des bornes avec recalcul de la direction si une
49// inconnue a une valeur imposee.
7fd59977 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.
7fd59977 60// - Si on a progresse sur F
61// On fait les tests d'arret
62// Sinon
63// On change de direction
64//============================================================================
65
3e42bd70
J		 66#define FSR_DEBUG(arg)
67// Uncomment the following code to have debug output to cout
258ff83b 68/*
3e42bd70
J		 69static Standard_Boolean mydebug = Standard_False;
70#undef FSR_DEBUG
71#define FSR_DEBUG(arg) {if (mydebug) { cout << arg << endl; }}
258ff83b 72*/
fbadd2cc 73
7fd59977 74class 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
83public :
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
103MyDirFunction::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
116void MyDirFunction::Initialize(const math_Vector& p0,
117 const math_Vector& dir) const
118{
119 *P0 = p0;
120 *Dir = dir;
121}
122
123
124Standard_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
154Standard_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);
12f139fd 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 }
7fd59977 190 Gnr1 = GH.Norm2();
191 return Standard_True;
192 }
193 return Standard_False;
194}
195
196
197//--------------------------------------------------------------
198static 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
3e42bd70 232 FSR_DEBUG (" minimisation dans la direction")
7fd59977 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//----------------------------------------------------------------------
248static 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;
3e42bd70 273 FSR_DEBUG(" essai d interpolation")
fbadd2cc 274
7fd59977 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
fbadd2cc 294 Delta = Sqrt(Delta);
7fd59977 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;
3e42bd70 315 FSR_DEBUG("t= "<<tsol<<" F = " << fsol << " OldF = "<<PValue)
7fd59977 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 }
3e42bd70 328 FSR_DEBUG(" minimisation dans la direction")
7fd59977 329 math_BrentMinimum Sol(F, ax, bx, cx, tol1d, 100, tol1d);
330 if(Sol.IsDone()) {
331 if (Sol.Minimum() <= Result) {
3e42bd70
J		 332 tsol = Sol.Location();
333 good = Standard_True;
334 FSR_DEBUG("t= "<<tsol<<" F ="<< Sol.Minimum() << " OldF = "<<Result)
7fd59977 335 }
336 }
337 }
338 if (good) {
339 // mise a jour du Delta
340 Dir.Multiply(tsol);
341 }
342 return good;
343}
344
345//------------------------------------------------------
346static 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.
3e42bd70
J		 365 FSR_DEBUG(" Matrice singuliere : On prend SVD")
366 math_SVD SolvebySVD(DF);
7fd59977 367 if (SolvebySVD.IsDone()) SolvebySVD.Solve(-1*FF, Direction);
3e42bd70
J		 368 else ChangeDirection = Standard_True;
369 }
7fd59977 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//=====================================================================
411static 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//====================================================
3e42bd70
J		 486Standard_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)
7fd59977 494//
3e42bd70
J		 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
7fd59977 500//======================================================
501{
502 Standard_Boolean Out = Standard_False;
503 Standard_Integer i, Ninc = Sol.Length();
504 Standard_Real monratio = 1;
505
3e42bd70
J		 506 theIsNewSol = Standard_True;
507
7fd59977 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;
3e42bd70
J		 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) );
7fd59977 522 }
523 else if (Sol(i) > SupBound(i)) {
524 Constraints(i) = 1;
525 Out = Standard_True;
3e42bd70
J		 526 // Delta(i) is positive
527 if (Delta(i) > Tol(i))
528 monratio = Min(monratio, (SupBound(i) - SolSave(i))/Delta(i) );
7fd59977 529 }
530 }
531
532 if (Out){ // Troncature et derniers recadrage pour blinder (pb numeriques)
3e42bd70
J		 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 }
7fd59977 549 }
550 }
551 }
552 return Out;
553}
554
555
556
557
558
559math_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
591math_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
621math_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,
b659a6dc 626 const Standard_Integer NbIterations,
627 Standard_Boolean theStopOnDivergent) :
7fd59977 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;
b659a6dc 657 Perform(F, StartingPoint, infBound, supBound, theStopOnDivergent);
7fd59977 658}
659
660math_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
696void math_FunctionSetRoot::Delete()
697{}
698
699void 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
706void math_FunctionSetRoot::Perform(math_FunctionSetWithDerivatives& F,
707 const math_Vector& StartingPoint,
708 const math_Vector& InfBound,
b659a6dc 709 const math_Vector& SupBound,
710 Standard_Boolean theStopOnDivergent)
7fd59977 711{
7fd59977 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;
3e42bd70 720 Standard_Boolean ChangeDirection = Standard_False, Sort = Standard_False, isNewSol = Standard_False;
7fd59977 721 Standard_Boolean Good, Verif;
722 Standard_Boolean Stop;
3e42bd70 723 const Standard_Real EpsSqrt = 1.e-16, Eps = 1.e-32, Eps2 = 1.e-64, Progres = 0.005;
7fd59977 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
b659a6dc 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 }
7fd59977 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;
b659a6dc 761 if (!theStopOnDivergent || !myIsDivergent)
762 {
763 State = F.GetStateNumber();
764 }
7fd59977 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...
3e42bd70 771 Save(0) = Max (F2, EpsSqrt);
5368adff 772 Standard_Real aTol_Func = Epsilon(F2);
3e42bd70
J		 773 FSR_DEBUG("=== Mode Debug de Function Set Root" << endl)
774 FSR_DEBUG(" F2 Initial = " << F2)
7fd59977 775
3e42bd70 776 if ((F2 <= Eps) || (Gnr1 <= Eps2)) {
b659a6dc 777 Done = Standard_False;
778 if (!theStopOnDivergent || !myIsDivergent)
779 {
780 Done = Standard_True;
781 State = F.GetStateNumber();
782 }
7fd59977 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) {
b659a6dc 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 }
7fd59977 803 return;
804 }
805 if (ChangeDirection) {
3e42bd70 806 Ambda = Ambda2 / Sqrt(Abs(Dy));
7fd59977 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 }
b659a6dc 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 //
7fd59977 827 Sort = Bounds(InfBound, SupBound, Tol, Sol, SolSave,
3e42bd70 828 Constraints, Delta, isNewSol);
7fd59977 829
830
831 DHSave = GH;
3e42bd70
J		 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;
b659a6dc 836 if (!theStopOnDivergent || !myIsDivergent)
837 {
838 State = F.GetStateNumber();
839 }
3e42bd70
J		 840 return;
841 }
7fd59977 842 }
3e42bd70
J		 843
844 FSR_DEBUG("Kount = " << Kount)
845 FSR_DEBUG("Le premier F2 = " << F2)
846 FSR_DEBUG("Direction = " << ChangeDirection)
7fd59977 847
3e42bd70 848 if ((F2 <= Eps) || (Gnr1 <= Eps2)) {
b659a6dc 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 }
7fd59977 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.
3e42bd70 876 FSR_DEBUG(" iteration de descente = " << DescenteIter)
7fd59977 877 DescenteIter++;
878 SolSave = Sol;
879 OldF = F2;
880 for( i = Sol.Lower(); i <= Sol.Upper(); i++) {
3e42bd70
J		 881 Sol(i) = Sol(i) + Ambda * DH(i);
882 }
b659a6dc 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 //
7fd59977 893 Stop = Bounds(InfBound, SupBound, Tol, Sol, SolSave,
3e42bd70
J		 894 Constraints, Delta, isNewSol);
895 FSR_DEBUG(" Augmentation de lambda")
7fd59977 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 {
3e42bd70 920 Sol = SolSave+Delta;
b659a6dc 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 //
3e42bd70
J		 931 Sort = Bounds(InfBound, SupBound, Tol, Sol, SolSave,
932 Constraints, Delta, isNewSol);
7fd59977 933 }
7fd59977 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;
3e42bd70
J		 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;
b659a6dc 943 if (!theStopOnDivergent || !myIsDivergent)
944 {
945 State = F.GetStateNumber();
946 }
3e42bd70
J		 947 return;
948 }
949 }
7fd59977 950 Dy = GH*DH;
951 if (Abs(Dy) <= Eps) {
3e42bd70
J		 952 if (F2 > OldF)
953 Sol = SolSave;
b659a6dc 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 }
7fd59977 963 return;
964 }
965 if (DescenteIter >= 10) {
966 Stop = Standard_True;
967 }
968 }
3e42bd70
J		 969 FSR_DEBUG("--- Sortie du Traitement Standard")
970 FSR_DEBUG(" DescenteIter = "<<DescenteIter << " F2 = " << F2)
7fd59977 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);
3e42bd70 987 FSR_DEBUG(" Conditional Direction = " << ChangeDirection)
7fd59977 988 if (Dy<-Eps) { //Pour eviter des calculs inutiles et des /0...
989 if (ChangeDirection) {
990
3e42bd70
J		 991// Ambda = Ambda2 / Sqrt(Abs(Dy));
992 Ambda = Ambda2 / Sqrt(-Dy);
7fd59977 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 }
b659a6dc 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 //
7fd59977 1013 Sortbis = Bounds(InfBound, SupBound, Tol, Sol, SolSave,
3e42bd70 1014 Constraints, Delta, isNewSol);
7fd59977 1015
1016 DHSave = GH;
3e42bd70
J		 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;
b659a6dc 1021 if (!theStopOnDivergent || !myIsDivergent)
1022 {
1023 State = F.GetStateNumber();
1024 }
3e42bd70
J		 1025 return;
1026 }
1027 }
7fd59977 1028 Ambda2 = Gnr1;
3e42bd70
J		 1029 FSR_DEBUG("--- Iteration au bords : " << DescenteIter)
1030 FSR_DEBUG("--- F2 = " << F2)
7fd59977 1031 }
1032 else {
1033 Stop = Standard_True;
1034 }
1035
1036 while((F2/PreviousMinimum > Progres) && (F2<OldF) && (!Stop) ) {
1037 DescenteIter++;
3e42bd70 1038 FSR_DEBUG("--- Iteration de descente conditionnel = " << DescenteIter)
7fd59977 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 }
b659a6dc 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 //
7fd59977 1056 Sortbis = Bounds(InfBound, SupBound, Tol, Sol, SolSave,
3e42bd70 1057 Constraints, Delta, isNewSol);
7fd59977 1058 }
1059 DHSave = GH;
3e42bd70
J		 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;
b659a6dc 1064 if (!theStopOnDivergent || !myIsDivergent)
1065 {
1066 State = F.GetStateNumber();
1067 }
3e42bd70
J		 1068 return;
1069 }
1070 }
7fd59977 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;
3e42bd70 1085 Sort = Standard_False;
7fd59977 1086 }
1087 else {
1088 Sol = SolSave + Delta;
b659a6dc 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 //
3e42bd70
J		 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;
b659a6dc 1105 if (!theStopOnDivergent || !myIsDivergent)
1106 {
1107 State = F.GetStateNumber();
1108 }
3e42bd70
J		 1109 return;
1110 }
1111 }
7fd59977 1112 }
1113 Dy = GH*DH;
1114 }
3e42bd70
J		 1115 FSR_DEBUG("--- Sortie du Traitement des Bords")
1116 FSR_DEBUG("--- DescenteIter = "<<DescenteIter << " F2 = " << F2)
7fd59977 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 }
3e42bd70 1138
7fd59977 1139 if (IsSolutionReached(F)) {
1140 if (PreviousMinimum < F2) {
1141 Sol = SolSave;
1142 }
b659a6dc 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 }
7fd59977 1152 return;
1153 }
1154 }
1155 //fin du test solution
1156
1157 // Analyse de la progression...
5368adff 1158 //comparison of current minimum and previous minimum
1159 if ((F2 - PreviousMinimum) <= aTol_Func){
7fd59977 1160 if (Kount > 5) {
1161 // L'historique est il bon ?
1162 if (F2 >= 0.95*Save(Kount - 5)) {
1163 if (!ChangeDirection) ChangeDirection = Standard_True;
5368adff 1164 else
1165 {
b659a6dc 1166 Done = Standard_False;
1167 if (!theStopOnDivergent || !myIsDivergent)
1168 {
1169 Done = Standard_True;
1170 State = F.GetStateNumber();
1171 }
5368adff 1172 return; // si un gain inf a 5% on sort
1173 }
7fd59977 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)) {
b659a6dc 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 }
7fd59977 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;
b659a6dc 1214 if (!theStopOnDivergent || !myIsDivergent)
1215 {
1216 State = F.GetStateNumber();
1217 }
7fd59977 1218 return;
1219 }
1220 }
5368adff 1221 else
1222 {
1223
b659a6dc 1224 if (!theStopOnDivergent || !myIsDivergent)
1225 {
1226 State = F.GetStateNumber();
1227 }
5368adff 1228 return; // y a plus d'issues
1229 }
7fd59977 1230 }
1231 }
b659a6dc 1232 if (!theStopOnDivergent || !myIsDivergent)
1233 {
1234 State = F.GetStateNumber();
1235 }
7fd59977 1236}
1237
1238
1239
1240
1241Standard_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
1249void 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
1262void 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
1269void 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}
Open CASCADE Technology repository