add reconstructedMatrix() to LLT, and LUs

=> they show that some improvements have still to be done
   for permutations, tr*tr, trapezoidal matrices

4 files changed, 63 insertions, 2 deletions

diff --git a/Eigen/src/Cholesky/LDLT.h b/Eigen/src/Cholesky/LDLT.h
index 8cfc256bb..8699fe7e0 100644
--- a/Eigen/src/Cholesky/LDLT.h
+++ b/Eigen/src/Cholesky/LDLT.h
@@ -155,7 +155,7 @@ template<typename _MatrixType> class LDLT
       return m_matrix;
     }
 
-    const MatrixType reconstructedMatrix() const;
+    MatrixType reconstructedMatrix() const;
 
     inline int rows() const { return m_matrix.rows(); }
     inline int cols() const { return m_matrix.cols(); }
@@ -324,7 +324,7 @@ bool LDLT<MatrixType>::solveInPlace(MatrixBase<Derived> &bAndX) const
  * i.e., it returns the product: P^T L D L^* P.
  * This function is provided for debug purpose. */
 template<typename MatrixType>
-const MatrixType LDLT<MatrixType>::reconstructedMatrix() const
+MatrixType LDLT<MatrixType>::reconstructedMatrix() const
 {
   ei_assert(m_isInitialized && "LDLT is not initialized.");
   const int size = m_matrix.rows();
diff --git a/Eigen/src/Cholesky/LLT.h b/Eigen/src/Cholesky/LLT.h
index 96e1e5f73..2e8df7661 100644
--- a/Eigen/src/Cholesky/LLT.h
+++ b/Eigen/src/Cholesky/LLT.h
@@ -133,6 +133,8 @@ template<typename _MatrixType, int _UpLo> class LLT
       return m_matrix;
     }
 
+    MatrixType reconstructedMatrix() const;
+
     inline int rows() const { return m_matrix.rows(); }
     inline int cols() const { return m_matrix.cols(); }
 
@@ -295,6 +297,16 @@ bool LLT<MatrixType,_UpLo>::solveInPlace(MatrixBase<Derived> &bAndX) const
   return true;
 }
 
+/** \returns the matrix represented by the decomposition,
+ * i.e., it returns the product: L L^*.
+ * This function is provided for debug purpose. */
+template<typename MatrixType, int _UpLo>
+MatrixType LLT<MatrixType,_UpLo>::reconstructedMatrix() const
+{
+  ei_assert(m_isInitialized && "LLT is not initialized.");
+  return matrixL() * matrixL().adjoint().toDenseMatrix();
+}
+
 /** \cholesky_module
   * \returns the LLT decomposition of \c *this
   */
diff --git a/Eigen/src/LU/FullPivLU.h b/Eigen/src/LU/FullPivLU.h
index 0a305d52b..cd63b9ec7 100644
--- a/Eigen/src/LU/FullPivLU.h
+++ b/Eigen/src/LU/FullPivLU.h
@@ -361,6 +361,8 @@ template<typename _MatrixType> class FullPivLU
                (*this, MatrixType::Identity(m_lu.rows(), m_lu.cols()));
     }
 
+    MatrixType reconstructedMatrix() const;
+
     inline int rows() const { return m_lu.rows(); }
     inline int cols() const { return m_lu.cols(); }
 
@@ -487,6 +489,33 @@ typename ei_traits<MatrixType>::Scalar FullPivLU<MatrixType>::determinant() cons
   return Scalar(m_det_pq) * Scalar(m_lu.diagonal().prod());
 }
 
+/** \returns the matrix represented by the decomposition,
+ * i.e., it returns the product: P^{-1} L U Q^{-1}.
+ * This function is provided for debug purpose. */
+template<typename MatrixType>
+MatrixType FullPivLU<MatrixType>::reconstructedMatrix() const
+{
+  ei_assert(m_isInitialized && "LU is not initialized.");
+  const int smalldim = std::min(m_lu.rows(), m_lu.cols());
+  // LU
+  MatrixType res(m_lu.rows(),m_lu.cols());
+  // FIXME the .toDenseMatrix() should not be needed...
+  res = m_lu.corner(TopLeft,m_lu.rows(),smalldim)
+            .template triangularView<UnitLower>().toDenseMatrix()
+      * m_lu.corner(TopLeft,smalldim,m_lu.cols())
+            .template triangularView<Upper>().toDenseMatrix();
+
+  // P^{-1}(LU)
+  // FIXME implement inplace permutation
+  res = (m_p.inverse() * res).eval();
+
+  // (P^{-1}LU)Q^{-1}
+  // FIXME implement inplace permutation
+  res = (res * m_q.inverse()).eval();
+
+  return res;
+}
+
 /********* Implementation of kernel() **************************************************/
 
 template<typename _MatrixType>
diff --git a/Eigen/src/LU/PartialPivLU.h b/Eigen/src/LU/PartialPivLU.h
index ed2354d78..fcffc2458 100644
--- a/Eigen/src/LU/PartialPivLU.h
+++ b/Eigen/src/LU/PartialPivLU.h
@@ -165,6 +165,8 @@ template<typename _MatrixType> class PartialPivLU
       */
     typename ei_traits<MatrixType>::Scalar determinant() const;
 
+    MatrixType reconstructedMatrix() const;
+
     inline int rows() const { return m_lu.rows(); }
     inline int cols() const { return m_lu.cols(); }
 
@@ -400,6 +402,24 @@ typename ei_traits<MatrixType>::Scalar PartialPivLU<MatrixType>::determinant() c
   return Scalar(m_det_p) * m_lu.diagonal().prod();
 }
 
+/** \returns the matrix represented by the decomposition,
+ * i.e., it returns the product: P^{-1} L U.
+ * This function is provided for debug purpose. */
+template<typename MatrixType>
+MatrixType PartialPivLU<MatrixType>::reconstructedMatrix() const
+{
+  ei_assert(m_isInitialized && "LU is not initialized.");
+  // LU
+  MatrixType res = m_lu.template triangularView<UnitLower>().toDenseMatrix()
+                 * m_lu.template triangularView<Upper>();
+  
+  // P^{-1}(LU)
+  // FIXME implement inplace permutation
+  res = (m_p.inverse() * res).eval();
+
+  return res;
+}
+
 /***** Implementation of solve() *****************************************************/
 
 template<typename _MatrixType, typename Rhs>