QR and SVD decomposition interface unification.

Added default ctor and public compute method as
well as safe-guards against uninitialized usage.
Added unit tests for the safe-guards.

Eigen/src/QR/QR.h

@@ -49,11 +49,20 @@ template<typename MatrixType> class QR
     typedef Matrix<Scalar, MatrixType::ColsAtCompileTime, MatrixType::ColsAtCompileTime> MatrixTypeR;
     typedef Matrix<Scalar, MatrixType::ColsAtCompileTime, 1> VectorType;
 
+    /** 
+    * \brief Default Constructor.
+    *
+    * The default constructor is useful in cases in which the user intends to
+    * perform decompositions via QR::compute(const MatrixType&).
+    */
+    QR() : m_qr(), m_hCoeffs(), m_isInitialized(false) {}
+
     QR(const MatrixType& matrix)
       : m_qr(matrix.rows(), matrix.cols()),
-        m_hCoeffs(matrix.cols())
+        m_hCoeffs(matrix.cols()),
+        m_isInitialized(false)
     {
-      _compute(matrix);
+      compute(matrix);
     }
     
     /** \deprecated use isInjective()
@@ -62,7 +71,11 @@ template<typename MatrixType> class QR
       * \note Since the rank is computed only once, i.e. the first time it is needed, this
       *       method almost does not perform any further computation.
       */
-    EIGEN_DEPRECATED bool isFullRank() const { return rank() == m_qr.cols(); }
+    EIGEN_DEPRECATED bool isFullRank() const 
+    { 
+      ei_assert(m_isInitialized && "QR is not initialized.");
+      return rank() == m_qr.cols(); 
+    }
     
     /** \returns the rank of the matrix of which *this is the QR decomposition.
       *
@@ -78,6 +91,7 @@ template<typename MatrixType> class QR
       */
     inline int dimensionOfKernel() const
     {
+      ei_assert(m_isInitialized && "QR is not initialized.");
       return m_qr.cols() - rank();
     }
     
@@ -89,6 +103,7 @@ template<typename MatrixType> class QR
       */
     inline bool isInjective() const
     {
+      ei_assert(m_isInitialized && "QR is not initialized.");
       return rank() == m_qr.cols();
     }
     
@@ -100,6 +115,7 @@ template<typename MatrixType> class QR
       */
     inline bool isSurjective() const
     {
+      ei_assert(m_isInitialized && "QR is not initialized.");
       return rank() == m_qr.rows();
     }
 
@@ -110,6 +126,7 @@ template<typename MatrixType> class QR
       */
     inline bool isInvertible() const
     {
+      ei_assert(m_isInitialized && "QR is not initialized.");
       return isInjective() && isSurjective();
     }
     
@@ -117,6 +134,7 @@ template<typename MatrixType> class QR
     const Part<NestByValue<MatrixRBlockType>, UpperTriangular>
     matrixR(void) const
     {
+      ei_assert(m_isInitialized && "QR is not initialized.");
       int cols = m_qr.cols();
       return MatrixRBlockType(m_qr, 0, 0, cols, cols).nestByValue().template part<UpperTriangular>();
     }
@@ -149,21 +167,21 @@ template<typename MatrixType> class QR
 
     MatrixType matrixQ(void) const;
 
-  private:
+    void compute(const MatrixType& matrix);
-
-    void _compute(const MatrixType& matrix);
 
   protected:
     MatrixType m_qr;
     VectorType m_hCoeffs;
     mutable int m_rank;
     mutable bool m_rankIsUptodate;
+    bool m_isInitialized;
 };
 
 /** \returns the rank of the matrix of which *this is the QR decomposition. */
 template<typename MatrixType>
 int QR<MatrixType>::rank() const
 {
+  ei_assert(m_isInitialized && "QR is not initialized.");
   if (!m_rankIsUptodate)
   {
     RealScalar maxCoeff = m_qr.diagonal().cwise().abs().maxCoeff();
@@ -179,10 +197,12 @@ int QR<MatrixType>::rank() const
 #ifndef EIGEN_HIDE_HEAVY_CODE
 
 template<typename MatrixType>
-void QR<MatrixType>::_compute(const MatrixType& matrix)
+void QR<MatrixType>::compute(const MatrixType& matrix)
 { 
   m_rankIsUptodate = false;
   m_qr = matrix;
+  m_hCoeffs.resize(matrix.cols());
+
   int rows = matrix.rows();
   int cols = matrix.cols();
   RealScalar eps2 = precision<RealScalar>()*precision<RealScalar>();
@@ -237,6 +257,7 @@ void QR<MatrixType>::_compute(const MatrixType& matrix)
       m_hCoeffs.coeffRef(k) = 0;
     }
   }
+  m_isInitialized = true;
 }
 
 template<typename MatrixType>
@@ -246,6 +267,7 @@ bool QR<MatrixType>::solve(
   ResultType *result
 ) const
 {
+  ei_assert(m_isInitialized && "QR is not initialized.");
   const int rows = m_qr.rows();
   ei_assert(b.rows() == rows);
   result->resize(rows, b.cols());
@@ -274,6 +296,7 @@ bool QR<MatrixType>::solve(
 template<typename MatrixType>
 MatrixType QR<MatrixType>::matrixQ() const
 {
+  ei_assert(m_isInitialized && "QR is not initialized.");
   // compute the product Q_0 Q_1 ... Q_n-1,
   // where Q_k is the k-th Householder transformation I - h_k v_k v_k'
   // and v_k is the k-th Householder vector [1,m_qr(k+1,k), m_qr(k+2,k), ...]

Eigen/src/SVD/SVD.h

@@ -61,10 +61,19 @@ template<typename MatrixType> class SVD
 
   public:
 
+    /** 
+    * \brief Default Constructor.
+    *
+    * The default constructor is useful in cases in which the user intends to
+    * perform decompositions via QR::compute(const MatrixType&).
+    */
+    SVD() : m_matU(), m_matV(), m_sigma(), m_isInitialized(false) {}
+
     SVD(const MatrixType& matrix)
       : m_matU(matrix.rows(), std::min(matrix.rows(), matrix.cols())),
         m_matV(matrix.cols(),matrix.cols()),
-        m_sigma(std::min(matrix.rows(),matrix.cols()))
+        m_sigma(std::min(matrix.rows(),matrix.cols())),
+        m_isInitialized(false)
     {
       compute(matrix);
     }
@@ -72,9 +81,23 @@ template<typename MatrixType> class SVD
     template<typename OtherDerived, typename ResultType>
     bool solve(const MatrixBase<OtherDerived> &b, ResultType* result) const;
 
-    const MatrixUType& matrixU() const { return m_matU; }
+    const MatrixUType& matrixU() const 
-    const SingularValuesType& singularValues() const { return m_sigma; }
+    { 
-    const MatrixVType& matrixV() const { return m_matV; }
+      ei_assert(m_isInitialized && "SVD is not initialized.");
+      return m_matU; 
+    }
+
+    const SingularValuesType& singularValues() const 
+    {
+      ei_assert(m_isInitialized && "SVD is not initialized.");
+      return m_sigma; 
+    }
+
+    const MatrixVType& matrixV() const 
+    {
+      ei_assert(m_isInitialized && "SVD is not initialized.");
+      return m_matV; 
+    }
 
     void compute(const MatrixType& matrix);
     SVD& sort();
@@ -95,6 +118,7 @@ template<typename MatrixType> class SVD
     MatrixVType m_matV;
     /** \internal */
     SingularValuesType m_sigma;
+    bool m_isInitialized;
 };
 
 /** Computes / recomputes the SVD decomposition A = U S V^* of \a matrix
@@ -473,11 +497,15 @@ void SVD<MatrixType>::compute(const MatrixType& matrix)
       break;
     } // end big switch
   } // end iterations
+
+  m_isInitialized = true;
 }
 
 template<typename MatrixType>
 SVD<MatrixType>& SVD<MatrixType>::sort()
 {
+  ei_assert(m_isInitialized && "SVD is not initialized.");
+
   int mu = m_matU.rows();
   int mv = m_matV.rows();
   int n  = m_matU.cols();
@@ -521,6 +549,8 @@ template<typename MatrixType>
 template<typename OtherDerived, typename ResultType>
 bool SVD<MatrixType>::solve(const MatrixBase<OtherDerived> &b, ResultType* result) const
 {
+  ei_assert(m_isInitialized && "SVD is not initialized.");
+
   const int rows = m_matU.rows();
   ei_assert(b.rows() == rows);
 
@@ -556,6 +586,7 @@ template<typename UnitaryType, typename PositiveType>
 void SVD<MatrixType>::computeUnitaryPositive(UnitaryType *unitary,
                                              PositiveType *positive) const
 {
+  ei_assert(m_isInitialized && "SVD is not initialized.");
   ei_assert(m_matU.cols() == m_matV.cols() && "Polar decomposition is only for square matrices");
   if(unitary) *unitary = m_matU * m_matV.adjoint();
   if(positive) *positive = m_matV * m_sigma.asDiagonal() * m_matV.adjoint();
@@ -574,6 +605,7 @@ template<typename UnitaryType, typename PositiveType>
 void SVD<MatrixType>::computePositiveUnitary(UnitaryType *positive,
                                              PositiveType *unitary) const
 {
+  ei_assert(m_isInitialized && "SVD is not initialized.");
   ei_assert(m_matU.rows() == m_matV.rows() && "Polar decomposition is only for square matrices");
   if(unitary) *unitary = m_matU * m_matV.adjoint();
   if(positive) *positive = m_matU * m_sigma.asDiagonal() * m_matU.adjoint();
@@ -592,6 +624,7 @@ template<typename MatrixType>
 template<typename RotationType, typename ScalingType>
 void SVD<MatrixType>::computeRotationScaling(RotationType *rotation, ScalingType *scaling) const
 {
+  ei_assert(m_isInitialized && "SVD is not initialized.");
   ei_assert(m_matU.rows() == m_matV.rows() && "Polar decomposition is only for square matrices");
   Scalar x = (m_matU * m_matV.adjoint()).determinant(); // so x has absolute value 1
   Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> sv(m_sigma);
@@ -618,6 +651,7 @@ template<typename MatrixType>
 template<typename ScalingType, typename RotationType>
 void SVD<MatrixType>::computeScalingRotation(ScalingType *scaling, RotationType *rotation) const
 {
+  ei_assert(m_isInitialized && "SVD is not initialized.");
   ei_assert(m_matU.rows() == m_matV.rows() && "Polar decomposition is only for square matrices");
   Scalar x = (m_matU * m_matV.adjoint()).determinant(); // so x has absolute value 1
   Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> sv(m_sigma);

test/qr.cpp

@@ -121,6 +121,22 @@ template<typename MatrixType> void qr_invertible()
   VERIFY(lu.solve(m3, &m2));
 }
 
+template<typename MatrixType> void qr_verify_assert()
+{
+  MatrixType tmp;
+
+  QR<MatrixType> qr;
+  VERIFY_RAISES_ASSERT(qr.isFullRank())
+  VERIFY_RAISES_ASSERT(qr.rank())
+  VERIFY_RAISES_ASSERT(qr.dimensionOfKernel())
+  VERIFY_RAISES_ASSERT(qr.isInjective())
+  VERIFY_RAISES_ASSERT(qr.isSurjective())
+  VERIFY_RAISES_ASSERT(qr.isInvertible())
+  VERIFY_RAISES_ASSERT(qr.matrixR())
+  VERIFY_RAISES_ASSERT(qr.solve(tmp,&tmp))
+  VERIFY_RAISES_ASSERT(qr.matrixQ())
+}
+
 void test_qr()
 {
   for(int i = 0; i < 1; i++) {
@@ -144,4 +160,11 @@ void test_qr()
 //     CALL_SUBTEST( qr_invertible<MatrixXcf>() );
 //     CALL_SUBTEST( qr_invertible<MatrixXcd>() );
   }
+
+  CALL_SUBTEST(qr_verify_assert<Matrix3f>());
+  CALL_SUBTEST(qr_verify_assert<Matrix3d>());
+  CALL_SUBTEST(qr_verify_assert<MatrixXf>());
+  CALL_SUBTEST(qr_verify_assert<MatrixXd>());
+  CALL_SUBTEST(qr_verify_assert<MatrixXcf>());
+  CALL_SUBTEST(qr_verify_assert<MatrixXcd>());
 }

test/svd.cpp

@@ -24,6 +24,7 @@
 
 #include "main.h"
 #include <Eigen/SVD>
+#include <Eigen/LU>
 
 template<typename MatrixType> void svd(const MatrixType& m)
 {
@@ -85,6 +86,22 @@ template<typename MatrixType> void svd(const MatrixType& m)
   }
 }
 
+template<typename MatrixType> void svd_verify_assert()
+{
+  MatrixType tmp;
+
+  SVD<MatrixType> svd;
+  VERIFY_RAISES_ASSERT(svd.solve(tmp, &tmp))
+  VERIFY_RAISES_ASSERT(svd.matrixU())
+  VERIFY_RAISES_ASSERT(svd.singularValues())
+  VERIFY_RAISES_ASSERT(svd.matrixV())
+  VERIFY_RAISES_ASSERT(svd.sort())
+  VERIFY_RAISES_ASSERT(svd.computeUnitaryPositive(&tmp,&tmp))
+  VERIFY_RAISES_ASSERT(svd.computePositiveUnitary(&tmp,&tmp))
+  VERIFY_RAISES_ASSERT(svd.computeRotationScaling(&tmp,&tmp))
+  VERIFY_RAISES_ASSERT(svd.computeScalingRotation(&tmp,&tmp))
+}
+
 void test_svd()
 {
   for(int i = 0; i < g_repeat; i++) {
@@ -96,4 +113,9 @@ void test_svd()
 //     CALL_SUBTEST( svd(MatrixXcd(6,6)) );
 //     CALL_SUBTEST( svd(MatrixXcf(3,3)) );
   }
+
+  CALL_SUBTEST( svd_verify_assert<Matrix3f>() );
+  CALL_SUBTEST( svd_verify_assert<Matrix3d>() );
+  CALL_SUBTEST( svd_verify_assert<MatrixXf>() );
+  CALL_SUBTEST( svd_verify_assert<MatrixXd>() );
 }