From f7aa9873caab20d49afd622b5ef76ecff8bfef06 Mon Sep 17 00:00:00 2001
From: Gael Guennebaud <g.gael@free.fr>
Date: Wed, 24 Feb 2010 10:40:16 +0100
Subject: [PATCH] * fix LDLT's default ctor use * add a reconstructedMatrix()
 function to LDLT for debug purpose

---
 Eigen/src/Cholesky/LDLT.h | 55 ++++++++++++++++++++++++++++++++-------
 Eigen/src/Cholesky/LLT.h  |  2 +-
 2 files changed, 47 insertions(+), 10 deletions(-)

diff --git a/Eigen/src/Cholesky/LDLT.h b/Eigen/src/Cholesky/LDLT.h
index 7c8e1eb04..8cfc256bb 100644
--- a/Eigen/src/Cholesky/LDLT.h
+++ b/Eigen/src/Cholesky/LDLT.h
@@ -62,14 +62,21 @@ template<typename _MatrixType> class LDLT
     typedef Matrix<int, MatrixType::RowsAtCompileTime, 1> IntColVectorType;
     typedef Matrix<int, 1, MatrixType::RowsAtCompileTime> IntRowVectorType;
 
-    /**
-    * \brief Default Constructor.
-    *
-    * The default constructor is useful in cases in which the user intends to
-    * perform decompositions via LDLT::compute(const MatrixType&).
-    */
+    /** \brief Default Constructor.
+      *
+      * The default constructor is useful in cases in which the user intends to
+      * perform decompositions via LDLT::compute(const MatrixType&).
+      */
     LDLT() : m_matrix(), m_p(), m_transpositions(), m_isInitialized(false) {}
 
+    /** \brief Default Constructor with memory preallocation
+      *
+      * Like the default constructor but with preallocation of the internal data
+      * according to the specified problem \a size.
+      * \sa LDLT()
+      */
+    LDLT(int size) : m_matrix(size,size), m_p(size), m_transpositions(size), m_isInitialized(false) {}
+
     LDLT(const MatrixType& matrix)
       : m_matrix(matrix.rows(), matrix.cols()),
         m_p(matrix.rows()),
@@ -148,6 +155,8 @@ template<typename _MatrixType> class LDLT
       return m_matrix;
     }
 
+    const MatrixType reconstructedMatrix() const;
+
     inline int rows() const { return m_matrix.rows(); }
     inline int cols() const { return m_matrix.cols(); }
 
@@ -175,6 +184,10 @@ LDLT<MatrixType>& LDLT<MatrixType>::compute(const MatrixType& a)
 
   m_matrix = a;
 
+  m_p.resize(size);
+  m_transpositions.resize(size);
+  m_isInitialized = false;
+
   if (size <= 1) {
     m_p.setZero();
     m_transpositions.setZero();
@@ -228,8 +241,7 @@ LDLT<MatrixType>& LDLT<MatrixType>::compute(const MatrixType& a)
       continue;
     }
 
-    RealScalar Djj = ei_real(m_matrix.coeff(j,j) -  m_matrix.row(j).head(j)
-                                               .dot(m_matrix.col(j).head(j)));
+    RealScalar Djj = ei_real(m_matrix.coeff(j,j) -  m_matrix.row(j).head(j).dot(m_matrix.col(j).head(j)));
     m_matrix.coeffRef(j,j) = Djj;
 
     int endSize = size - j - 1;
@@ -238,7 +250,7 @@ LDLT<MatrixType>& LDLT<MatrixType>::compute(const MatrixType& a)
                                 * m_matrix.col(j).head(j).conjugate();
 
       m_matrix.row(j).tail(endSize) = m_matrix.row(j).tail(endSize).conjugate()
-                                   - _temporary.tail(endSize).transpose();
+                                    - _temporary.tail(endSize).transpose();
 
       if(ei_abs(Djj) > cutoff)
       {
@@ -308,6 +320,31 @@ bool LDLT<MatrixType>::solveInPlace(MatrixBase<Derived> &bAndX) const
   return true;
 }
 
+/** \returns the matrix represented by the decomposition,
+ * 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
+{
+  ei_assert(m_isInitialized && "LDLT is not initialized.");
+  const int size = m_matrix.rows();
+  MatrixType res(size,size);
+  res.setIdentity();
+
+  // PI
+  for(int i = 0; i < size; ++i) res.row(m_transpositions.coeff(i)).swap(res.row(i));
+  // L^* P
+  res = matrixL().adjoint() * res;
+  // D(L^*P)
+  res = vectorD().asDiagonal() * res;
+  // L(DL^*P)
+  res = matrixL() * res;
+  // P^T (LDL^*P)
+  for (int i = size-1; i >= 0; --i) res.row(m_transpositions.coeff(i)).swap(res.row(i));
+
+  return res;
+}
+
 /** \cholesky_module
   * \returns the Cholesky decomposition with full pivoting without square root of \c *this
   */
diff --git a/Eigen/src/Cholesky/LLT.h b/Eigen/src/Cholesky/LLT.h
index 474b82406..96e1e5f73 100644
--- a/Eigen/src/Cholesky/LLT.h
+++ b/Eigen/src/Cholesky/LLT.h
@@ -117,7 +117,7 @@ template<typename _MatrixType, int _UpLo> class LLT
                 && "LLT::solve(): invalid number of rows of the right hand side matrix b");
       return ei_solve_retval<LLT, Rhs>(*this, b.derived());
     }
-    
+
     template<typename Derived>
     bool solveInPlace(MatrixBase<Derived> &bAndX) const;