LU and PartialLU 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/LU/LU.h

@@ -92,12 +92,22 @@ template<typename MatrixType> class LU
                    MatrixType::MaxColsAtCompileTime  // so it has the same number of rows and at most as many columns.
     > ImageResultType;
 
+    /** 
+    * \brief Default Constructor.
+    *
+    * The default constructor is useful in cases in which the user intends to
+    * perform decompositions via LU::compute(const MatrixType&).
+    */
+    LU();
+
     /** Constructor.
       *
       * \param matrix the matrix of which to compute the LU decomposition.
       */
     LU(const MatrixType& matrix);
 
+    void compute(const MatrixType& matrix);
+
     /** \returns the LU decomposition matrix: the upper-triangular part is U, the
       * unit-lower-triangular part is L (at least for square matrices; in the non-square
       * case, special care is needed, see the documentation of class LU).
@@ -106,6 +116,7 @@ template<typename MatrixType> class LU
       */
     inline const MatrixType& matrixLU() const
     {
+      ei_assert(m_originalMatrix != 0 && "LU is not initialized.");
       return m_lu;
     }
 
@@ -117,6 +128,7 @@ template<typename MatrixType> class LU
       */
     inline const IntColVectorType& permutationP() const
     {
+      ei_assert(m_originalMatrix != 0 && "LU is not initialized.");
       return m_p;
     }
 
@@ -128,6 +140,7 @@ template<typename MatrixType> class LU
       */
     inline const IntRowVectorType& permutationQ() const
     {
+      ei_assert(m_originalMatrix != 0 && "LU is not initialized.");
       return m_q;
     }
 
@@ -243,6 +256,7 @@ template<typename MatrixType> class LU
       */
     inline int rank() const
     {
+      ei_assert(m_originalMatrix != 0 && "LU is not initialized.");
       return m_rank;
     }
 
@@ -253,6 +267,7 @@ template<typename MatrixType> class LU
       */
     inline int dimensionOfKernel() const
     {
+      ei_assert(m_originalMatrix != 0 && "LU is not initialized.");
       return m_lu.cols() - m_rank;
     }
 
@@ -264,6 +279,7 @@ template<typename MatrixType> class LU
       */
     inline bool isInjective() const
     {
+      ei_assert(m_originalMatrix != 0 && "LU is not initialized.");
       return m_rank == m_lu.cols();
     }
 
@@ -275,6 +291,7 @@ template<typename MatrixType> class LU
       */
     inline bool isSurjective() const
     {
+      ei_assert(m_originalMatrix != 0 && "LU is not initialized.");
       return m_rank == m_lu.rows();
     }
 
@@ -285,6 +302,7 @@ template<typename MatrixType> class LU
       */
     inline bool isInvertible() const
     {
+      ei_assert(m_originalMatrix != 0 && "LU is not initialized.");
       return isInjective() && isSurjective();
     }
 
@@ -317,7 +335,7 @@ template<typename MatrixType> class LU
     }
 
   protected:
-    const MatrixType& m_originalMatrix;
+    const MatrixType* m_originalMatrix;
     MatrixType m_lu;
     IntColVectorType m_p;
     IntRowVectorType m_q;
@@ -327,12 +345,38 @@ template<typename MatrixType> class LU
 };
 
 template<typename MatrixType>
-LU<MatrixType>::LU(const MatrixType& matrix)
+LU<MatrixType>::LU()
-  : m_originalMatrix(matrix),
+  : m_originalMatrix(0),
-    m_lu(matrix),
+    m_lu(),
-    m_p(matrix.rows()),
+    m_p(),
-    m_q(matrix.cols())
+    m_q(),
+    m_det_pq(0),
+    m_rank(-1),
+    m_precision(precision<RealScalar>())
 {
+}
+
+template<typename MatrixType>
+LU<MatrixType>::LU(const MatrixType& matrix)
+  : m_originalMatrix(0),
+    m_lu(),
+    m_p(),
+    m_q(),
+    m_det_pq(0),
+    m_rank(-1),
+    m_precision(precision<RealScalar>())
+{
+  compute(matrix);
+}
+
+template<typename MatrixType>
+void LU<MatrixType>::compute(const MatrixType& matrix)
+{
+  m_originalMatrix = &matrix;
+  m_lu = matrix;
+  m_p.resize(matrix.rows());
+  m_q.resize(matrix.cols());
+
   const int size = matrix.diagonalSize();
   const int rows = matrix.rows();
   const int cols = matrix.cols();
@@ -409,6 +453,7 @@ LU<MatrixType>::LU(const MatrixType& matrix)
 template<typename MatrixType>
 typename ei_traits<MatrixType>::Scalar LU<MatrixType>::determinant() const
 {
+  ei_assert(m_originalMatrix != 0 && "LU is not initialized.");
   return Scalar(m_det_pq) * m_lu.diagonal().prod();
 }
 
@@ -462,17 +507,20 @@ template<typename MatrixType>
 template<typename ImageMatrixType>
 void LU<MatrixType>::computeImage(ImageMatrixType *result) const
 {
-  ei_assert(m_rank > 0);
+  // can be caused by a rank deficient matrix or by calling computeImage on an
-  result->resize(m_originalMatrix.rows(), m_rank);
+  // unitialized LU object
+  ei_assert(m_rank > 0 && "LU is not initialized or Matrix has rank zero.");
+  result->resize(m_originalMatrix->rows(), m_rank);
   for(int i = 0; i < m_rank; ++i)
-    result->col(i) = m_originalMatrix.col(m_q.coeff(i));
+    result->col(i) = m_originalMatrix->col(m_q.coeff(i));
 }
 
 template<typename MatrixType>
 const typename LU<MatrixType>::ImageResultType
 LU<MatrixType>::image() const
 {
-  ImageResultType result(m_originalMatrix.rows(), m_rank);
+  ei_assert(m_originalMatrix != 0 && "LU is not initialized.");
+  ImageResultType result(m_originalMatrix->rows(), m_rank);
   computeImage(&result);
   return result;
 }
@@ -484,6 +532,8 @@ bool LU<MatrixType>::solve(
   ResultType *result
 ) const
 {
+  ei_assert(m_originalMatrix != 0 && "LU is not initialized.");
+
   /* The decomposition PAQ = LU can be rewritten as A = P^{-1} L U Q^{-1}.
    * So we proceed as follows:
    * Step 1: compute c = Pb.

Eigen/src/LU/PartialLU.h

@@ -70,6 +70,14 @@ template<typename MatrixType> class PartialLU
              MatrixType::MaxRowsAtCompileTime)
     };
 
+    /** 
+    * \brief Default Constructor.
+    *
+    * The default constructor is useful in cases in which the user intends to
+    * perform decompositions via PartialLU::compute(const MatrixType&).
+    */
+    PartialLU();
+
     /** Constructor.
       *
       * \param matrix the matrix of which to compute the LU decomposition.
@@ -79,6 +87,8 @@ template<typename MatrixType> class PartialLU
       */
     PartialLU(const MatrixType& matrix);
 
+    void compute(const MatrixType& matrix);
+
     /** \returns the LU decomposition matrix: the upper-triangular part is U, the
       * unit-lower-triangular part is L (at least for square matrices; in the non-square
       * case, special care is needed, see the documentation of class LU).
@@ -87,8 +97,9 @@ template<typename MatrixType> class PartialLU
       */
     inline const MatrixType& matrixLU() const
     {
+      ei_assert(m_isInitialized && "PartialLU is not initialized.");
       return m_lu;
-    }
+    }    
 
     /** \returns a vector of integers, whose size is the number of rows of the matrix being decomposed,
       * representing the P permutation i.e. the permutation of the rows. For its precise meaning,
@@ -96,6 +107,7 @@ template<typename MatrixType> class PartialLU
       */
     inline const IntColVectorType& permutationP() const
     {
+      ei_assert(m_isInitialized && "PartialLU is not initialized.");
       return m_p;
     }
 
@@ -164,18 +176,37 @@ template<typename MatrixType> class PartialLU
     }
 
   protected:
-    const MatrixType& m_originalMatrix;
     MatrixType m_lu;
     IntColVectorType m_p;
     int m_det_p;
+    bool m_isInitialized;
 };
 
 template<typename MatrixType>
-PartialLU<MatrixType>::PartialLU(const MatrixType& matrix)
+PartialLU<MatrixType>::PartialLU()
-  : m_originalMatrix(matrix),
+  : m_lu(),
-    m_lu(matrix),
+    m_p(),
-    m_p(matrix.rows())
+    m_det_p(0),
+    m_isInitialized(false)
 {
+}
+
+template<typename MatrixType>
+PartialLU<MatrixType>::PartialLU(const MatrixType& matrix)
+  : m_lu(),
+    m_p(),
+    m_det_p(0),
+    m_isInitialized(false)
+{
+  compute(matrix);
+}
+
+template<typename MatrixType>
+void PartialLU<MatrixType>::compute(const MatrixType& matrix)
+{
+  m_lu = matrix;
+  m_p.resize(matrix.rows());
+
   ei_assert(matrix.rows() == matrix.cols() && "PartialLU is only for square (and moreover invertible) matrices");
   const int size = matrix.rows();
 
@@ -213,11 +244,14 @@ PartialLU<MatrixType>::PartialLU(const MatrixType& matrix)
     std::swap(m_p.coeffRef(k), m_p.coeffRef(rows_transpositions.coeff(k)));
 
   m_det_p = (number_of_transpositions%2) ? -1 : 1;
+
+  m_isInitialized = true;
 }
 
 template<typename MatrixType>
 typename ei_traits<MatrixType>::Scalar PartialLU<MatrixType>::determinant() const
 {
+  ei_assert(m_isInitialized && "PartialLU is not initialized.");
   return Scalar(m_det_p) * m_lu.diagonal().prod();
 }
 
@@ -228,6 +262,8 @@ void PartialLU<MatrixType>::solve(
   ResultType *result
 ) const
 {
+  ei_assert(m_isInitialized && "PartialLU is not initialized.");
+
   /* The decomposition PA = LU can be rewritten as A = P^{-1} L U.
    * So we proceed as follows:
    * Step 1: compute c = Pb.

test/lu.cpp

@@ -92,6 +92,37 @@ template<typename MatrixType> void lu_invertible()
   VERIFY(lu.solve(m3, &m2));
 }
 
+template<typename MatrixType> void lu_verify_assert()
+{
+  MatrixType tmp;
+
+  LU<MatrixType> lu;
+  VERIFY_RAISES_ASSERT(lu.matrixLU())
+  VERIFY_RAISES_ASSERT(lu.permutationP())
+  VERIFY_RAISES_ASSERT(lu.permutationQ())
+  VERIFY_RAISES_ASSERT(lu.computeKernel(&tmp))
+  VERIFY_RAISES_ASSERT(lu.computeImage(&tmp))
+  VERIFY_RAISES_ASSERT(lu.kernel())
+  VERIFY_RAISES_ASSERT(lu.image())
+  VERIFY_RAISES_ASSERT(lu.solve(tmp,&tmp))
+  VERIFY_RAISES_ASSERT(lu.determinant())
+  VERIFY_RAISES_ASSERT(lu.rank())
+  VERIFY_RAISES_ASSERT(lu.dimensionOfKernel())
+  VERIFY_RAISES_ASSERT(lu.isInjective())
+  VERIFY_RAISES_ASSERT(lu.isSurjective())
+  VERIFY_RAISES_ASSERT(lu.isInvertible())
+  VERIFY_RAISES_ASSERT(lu.computeInverse(&tmp))
+  VERIFY_RAISES_ASSERT(lu.inverse())
+
+  PartialLU<MatrixType> plu;
+  VERIFY_RAISES_ASSERT(plu.matrixLU())
+  VERIFY_RAISES_ASSERT(plu.permutationP())
+  VERIFY_RAISES_ASSERT(plu.solve(tmp,&tmp))
+  VERIFY_RAISES_ASSERT(plu.determinant())
+  VERIFY_RAISES_ASSERT(plu.computeInverse(&tmp))
+  VERIFY_RAISES_ASSERT(plu.inverse())
+}
+
 void test_lu()
 {
   for(int i = 0; i < g_repeat; i++) {
@@ -104,4 +135,11 @@ void test_lu()
     CALL_SUBTEST( lu_invertible<MatrixXcf>() );
     CALL_SUBTEST( lu_invertible<MatrixXcd>() );
   }
+
+  CALL_SUBTEST( lu_verify_assert<Matrix3f>() );
+  CALL_SUBTEST( lu_verify_assert<Matrix3d>() );
+  CALL_SUBTEST( lu_verify_assert<MatrixXf>() );
+  CALL_SUBTEST( lu_verify_assert<MatrixXd>() );
+  CALL_SUBTEST( lu_verify_assert<MatrixXcf>() );
+  CALL_SUBTEST( lu_verify_assert<MatrixXcd>() );
 }