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