src/math/math_SVD.cxx

   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_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>
  28
  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())
  33 {
  34    U.Init(0.0);
  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 ;
  39 }
  40
  41 void math_SVD::Solve (const math_Vector& B,
  42                       math_Vector& X,
  43                       const Standard_Real Eps)
  44 {
  45    StdFail_NotDone_Raise_if(!Done, " ");
  46    Standard_DimensionError_Raise_if((RowA != B.Length()) ||
  47                                     (X.Length() != Diag.Length()), " ");
  48
  49    math_Vector BB (1, U.RowNumber()) ;
  50    BB.Init(0.0);
  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++)
  54    {
  55       if (Diag(I) < wmin) Diag(I) = 0.0 ;
  56    }
  57    SVD_Solve (U, Diag, V, BB, X) ;
  58 }
  59
  60 void  math_SVD::PseudoInverse (math_Matrix& Result,
  61                                const Standard_Real Eps)
  62 {
  63    Standard_Integer i, j ;
  64
  65    StdFail_NotDone_Raise_if(!Done, " ");
  66
  67    Standard_Real wmin = Eps * Diag(Diag.Max());
  68    for (i=1; i<=Diag.Upper(); i++)
  69    {
  70       if (Diag(i) < wmin) Diag(i) = 0.0 ;
  71    }
  72
  73    Standard_Integer ColA = Diag.Length () ;
  74    math_Vector VNorme (1, U.RowNumber()) ;
  75    math_Vector Column (1, ColA) ;
  76
  77    for (j=1; j<=RowA; j++)
  78    {
  79       for (i=1; i<=VNorme.Upper(); i++) VNorme(i) = 0.0 ;
  80       VNorme(j) = 1.0 ;
  81       SVD_Solve (U, Diag, V, VNorme, Column) ;
  82       for (i=1; i<=ColA; i++) Result(i,j) = Column(i) ;
  83    }
  84 }
  85
  86
  87 void math_SVD::Dump(Standard_OStream& o) const {
  88
  89   o << "math_SVD";
  90   if(Done) {
  91     o << " Status = Done \n";
  92   }
  93   else {
  94     o << " Status = not Done \n";
  95   }
  96 }
  97