Changeset 4286

Timestamp:

Jan 30, 2023, 3:08:38 AM (8 weeks ago)

Author:

Peter

Message:

provide the SVD inverse function as a member function so the decomposition can be recycled.

Location:

trunk/yat/utility

Files:

• Matrix.cc (modified) (1 diff)
• SVD.cc (modified) (2 diffs)
• SVD.h (modified) (1 diff)

trunk/yat/utility/Matrix.cc

@@ -260 +260 @@
     assert(A.rows() == A.columns());
     SVD svd(A);
-    // A = U S V'
     svd.decompose();
-
-    // result = A^-1 = V * S^-1 * U'
-
-    DiagonalMatrix D(A.rows(), A.columns());
-    // If eigenvalue is zero, A is not invertable, perhaps that should
-    // trigger an exception
-    for (size_t i=0; i<D.rows(); ++i)
-      D(i) = 1.0 / svd.s()(i);
-
-    result = svd.V() * D * transpose(svd.U());
-
+    svd.inverse(result);
     assert(result.columns() == result.rows());
     assert(result.rows() == A.rows());

trunk/yat/utility/SVD.cc

@@ -30 +30 @@
 #include "SVD.h"
 
+#include "DiagonalMatrix.h"
 #include "Exception.h"
 #include "Vector.h"
@@ -109 +110 @@
 
 
+  void SVD::inverse(Matrix& result) const
+  {
+    assert(U_.rows()==U_.columns() && "input matrix must be square");
+    // A = U S V'
+
+    // result = A^-1 = V * S^-1 * U'
+
+    DiagonalMatrix D(s_.size(), s_.size());
+    // If eigenvalue is zero, A is not invertable, perhaps that should
+    // trigger an exception
+    for (size_t i=0; i<D.rows(); ++i)
+      D(i) = 1.0 / s_(i);
+
+    result = V_ * D * transpose(U_);
+  }
+
+
   const Vector& SVD::s(void) const
   {

trunk/yat/utility/SVD.h

@@ -101 +101 @@
 
     /**
+       Calculate the inverse of the Matrix provided in the
+       constructor. Behaviour is undefined if the input matrix is not
+       square.
+
+       \since New in yat 0.21
+
+       \see inverse_svd(const MatrixBase&, Matrix&)
+     */
+    void inverse(Matrix& result) const;
+
+    /**
        \brief Access to the s vector.