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 //#ifndef OCCT_DEBUG
16 #define No_Standard_RangeError
17 #define No_Standard_OutOfRange
18 #define No_Standard_DimensionError
19
20 //#endif
21
22 #include <math_GaussLeastSquare.hxx>
23 #include <math_Matrix.hxx>
24 #include <math_Recipes.hxx>
25 #include <Standard_DimensionError.hxx>
26 #include <StdFail_NotDone.hxx>
27
math_GaussLeastSquare(const math_Matrix & A,const Standard_Real MinPivot)28 math_GaussLeastSquare::math_GaussLeastSquare (const math_Matrix& A,
29 const Standard_Real MinPivot) :
30 LU(1, A.ColNumber(),
31 1, A.ColNumber()),
32 A2(1, A.ColNumber(),
33 1, A.RowNumber()),
34 Index(1, A.ColNumber()) {
35 A2 = A.Transposed();
36 LU.Multiply(A2, A);
37
38 Standard_Integer Error = LU_Decompose(LU, Index, D, MinPivot);
39 Done = (!Error) ? Standard_True : Standard_False;
40
41 }
42
Solve(const math_Vector & B,math_Vector & X) const43 void math_GaussLeastSquare::Solve(const math_Vector& B, math_Vector& X) const{
44 StdFail_NotDone_Raise_if(!Done, " ");
45 Standard_DimensionError_Raise_if((B.Length() != A2.ColNumber()) ||
46 (X.Length() != A2.RowNumber()), " ");
47
48 X.Multiply(A2, B);
49
50 LU_Solve(LU, Index, X);
51
52 return;
53 }
54
55
Dump(Standard_OStream & o) const56 void math_GaussLeastSquare::Dump(Standard_OStream& o) const {
57
58 o <<"math_GaussLeastSquare ";
59 if (Done) {
60 o << " Status = Done \n";
61 }
62 else {
63 o << "Status = not Done \n";
64 }
65 }
66
67
68