1 // Copyright (c) 1997-1999 Matra Datavision
2 // Copyright (c) 1999-2014 OPEN CASCADE SAS
4 // This file is part of Open CASCADE Technology software library.
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.
12 // Alternatively, this file may be used under the terms of Open CASCADE
13 // commercial license or contractual agreement.
16 #define No_Standard_RangeError
17 #define No_Standard_OutOfRange
18 #define No_Standard_DimensionError
22 #include <math_Matrix.hxx>
23 #include <math_Recipes.hxx>
24 #include <math_SVD.hxx>
25 #include <Standard_DimensionError.hxx>
26 #include <Standard_NotImplemented.hxx>
27 #include <StdFail_NotDone.hxx>
29 math_SVD::math_SVD (const math_Matrix& A) :
30 U (1, Max(A.RowNumber(),A.ColNumber()), 1, A.ColNumber()) ,
31 V (1, A.ColNumber(), 1, A.ColNumber()),
32 Diag (1, A.ColNumber())
35 RowA = A.RowNumber() ;
36 U.Set(1, A.RowNumber(), 1, A.ColNumber(), A) ;
37 Standard_Integer Error = SVD_Decompose(U, Diag, V) ;
38 Done = (!Error) ? Standard_True : Standard_False ;
41 void math_SVD::Solve (const math_Vector& B,
43 const Standard_Real Eps)
45 StdFail_NotDone_Raise_if(!Done, " ");
46 Standard_DimensionError_Raise_if((RowA != B.Length()) ||
47 (X.Length() != Diag.Length()), " ");
49 math_Vector BB (1, U.RowNumber()) ;
51 BB.Set (1, B.Length(), B) ;
52 Standard_Real wmin = Eps * Diag(Diag.Max());
53 for (Standard_Integer I = 1; I <= Diag.Upper(); I++)
55 if (Diag(I) < wmin) Diag(I) = 0.0 ;
57 SVD_Solve (U, Diag, V, BB, X) ;
60 void math_SVD::PseudoInverse (math_Matrix& Result,
61 const Standard_Real Eps)
63 Standard_Integer i, j ;
65 StdFail_NotDone_Raise_if(!Done, " ");
67 Standard_Real wmin = Eps * Diag(Diag.Max());
68 for (i=1; i<=Diag.Upper(); i++)
70 if (Diag(i) < wmin) Diag(i) = 0.0 ;
73 Standard_Integer ColA = Diag.Length () ;
74 math_Vector VNorme (1, U.RowNumber()) ;
75 math_Vector Column (1, ColA) ;
77 for (j=1; j<=RowA; j++)
79 for (i=1; i<=VNorme.Upper(); i++) VNorme(i) = 0.0 ;
81 SVD_Solve (U, Diag, V, VNorme, Column) ;
82 for (i=1; i<=ColA; i++) Result(i,j) = Column(i) ;
87 void math_SVD::Dump(Standard_OStream& o) const {
91 o << " Status = Done \n";
94 o << " Status = not Done \n";