Add a Transpositions class to ease the representation and

manipulation of permutations as a sequence of transpositions.
Make LDLT use it.

5 files changed, 351 insertions, 50 deletions

diff --git a/Eigen/Core b/Eigen/Core
index 595e1cb89..70c2c2207 100644
--- a/Eigen/Core
+++ b/Eigen/Core
@@ -260,6 +260,7 @@ using std::size_t;
 #include "src/Core/Diagonal.h"
 #include "src/Core/DiagonalProduct.h"
 #include "src/Core/PermutationMatrix.h"
+#include "src/Core/Transpositions.h"
 #include "src/Core/Redux.h"
 #include "src/Core/Visitor.h"
 #include "src/Core/Fuzzy.h"
diff --git a/Eigen/src/Cholesky/LDLT.h b/Eigen/src/Cholesky/LDLT.h
index 60cb98307..79657f086 100644
--- a/Eigen/src/Cholesky/LDLT.h
+++ b/Eigen/src/Cholesky/LDLT.h
@@ -1,7 +1,7 @@
 // This file is part of Eigen, a lightweight C++ template library
 // for linear algebra.
 //
-// Copyright (C) 2008 Gael Guennebaud <g.gael@free.fr>
+// Copyright (C) 2008-2010 Gael Guennebaud <g.gael@free.fr>
 // Copyright (C) 2009 Keir Mierle <mierle@gmail.com>
 // Copyright (C) 2009 Benoit Jacob <jacob.benoit.1@gmail.com>
 //
@@ -69,9 +69,12 @@ template<typename _MatrixType, int _UpLo> class LDLT
     typedef typename MatrixType::Scalar Scalar;
     typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar;
     typedef typename MatrixType::Index Index;
-    typedef typename ei_plain_col_type<MatrixType, Index>::type IntColVectorType;
+//     typedef typename ei_plain_col_type<MatrixType, Index>::type IntColVectorType;
     typedef Matrix<Scalar, RowsAtCompileTime, 1, Options, MaxRowsAtCompileTime, 1> TmpMatrixType;
 
+    typedef Transpositions<RowsAtCompileTime, MaxRowsAtCompileTime> TranspositionType;
+    typedef PermutationMatrix<RowsAtCompileTime, MaxRowsAtCompileTime> PermutationType;
+
     typedef LDLT_Traits<MatrixType,UpLo> Traits;
 
     /** \brief Default Constructor.
@@ -79,7 +82,7 @@ template<typename _MatrixType, int _UpLo> class LDLT
       * 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) {}
+    LDLT() : m_matrix(), m_transpositions(), m_isInitialized(false) {}
 
     /** \brief Default Constructor with memory preallocation
       *
@@ -89,7 +92,6 @@ template<typename _MatrixType, int _UpLo> class LDLT
       */
     LDLT(Index size)
       : m_matrix(size, size),
-        m_p(size),
         m_transpositions(size),
         m_temporary(size),
         m_isInitialized(false)
@@ -97,7 +99,6 @@ template<typename _MatrixType, int _UpLo> class LDLT
 
     LDLT(const MatrixType& matrix)
       : m_matrix(matrix.rows(), matrix.cols()),
-        m_p(matrix.rows()),
         m_transpositions(matrix.rows()),
         m_temporary(matrix.rows()),
         m_isInitialized(false)
@@ -119,14 +120,12 @@ template<typename _MatrixType, int _UpLo> class LDLT
       return Traits::getL(m_matrix);
     }
 
-    /** \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,
-      * see the examples given in the documentation of class FullPivLU.
+    /** \returns the permutation matrix P as a transposition sequence.
       */
-    inline const IntColVectorType& permutationP() const
+    inline const TranspositionType& transpositionsP() const
     {
       ei_assert(m_isInitialized && "LDLT is not initialized.");
-      return m_p;
+      return m_transpositions;
     }
 
     /** \returns the coefficients of the diagonal matrix D */
@@ -195,8 +194,7 @@ template<typename _MatrixType, int _UpLo> class LDLT
       * is not stored), and the diagonal entries correspond to D.
       */
     MatrixType m_matrix;
-    IntColVectorType m_p;
-    IntColVectorType m_transpositions; // FIXME do we really need to store permanently the transpositions?
+    TranspositionType m_transpositions;
     TmpMatrixType m_temporary;
     int m_sign;
     bool m_isInitialized;
@@ -206,18 +204,18 @@ template<int UpLo> struct ei_ldlt_inplace;
 
 template<> struct ei_ldlt_inplace<Lower>
 {
-  template<typename MatrixType, typename Transpositions, typename Workspace>
-  static bool unblocked(MatrixType& mat, Transpositions& transpositions, Workspace& temp, int* sign=0)
+  template<typename MatrixType, typename TranspositionType, typename Workspace>
+  static bool unblocked(MatrixType& mat, TranspositionType& transpositions, Workspace& temp, int* sign=0)
   {
     typedef typename MatrixType::Scalar Scalar;
     typedef typename MatrixType::RealScalar RealScalar;
     typedef typename MatrixType::Index Index;
     ei_assert(mat.rows()==mat.cols());
     const Index size = mat.rows();
-    
+
     if (size <= 1)
     {
-      transpositions.setZero();
+      transpositions.setIdentity();
       if(sign)
         *sign = ei_real(mat.coeff(0,0))>0 ? 1:-1;
       return true;
@@ -272,15 +270,15 @@ template<> struct ei_ldlt_inplace<Lower>
       if((rs>0) && (ei_abs(mat.coeffRef(k,k)) > cutoff))
         A21 /= mat.coeffRef(k,k);
     }
-    
+
     return true;
   }
 };
 
 template<> struct ei_ldlt_inplace<Upper>
 {
-  template<typename MatrixType, typename Transpositions, typename Workspace>
-  static EIGEN_STRONG_INLINE bool unblocked(MatrixType& mat, Transpositions& transpositions, Workspace& temp, int* sign=0)
+  template<typename MatrixType, typename TranspositionType, typename Workspace>
+  static EIGEN_STRONG_INLINE bool unblocked(MatrixType& mat, TranspositionType& transpositions, Workspace& temp, int* sign=0)
   {
     Transpose<MatrixType> matt(mat);
     return ei_ldlt_inplace<Lower>::unblocked(matt, transpositions, temp, sign);
@@ -313,22 +311,12 @@ LDLT<MatrixType,_UpLo>& LDLT<MatrixType,_UpLo>::compute(const MatrixType& a)
 
   m_matrix = a;
 
-  m_p.resize(size);
   m_transpositions.resize(size);
   m_isInitialized = false;
   m_temporary.resize(size);
 
   ei_ldlt_inplace<UpLo>::unblocked(m_matrix, m_transpositions, m_temporary, &m_sign);
 
-  if(size==0)
-    m_p.setZero();
-
-  // Reverse applied swaps to get P matrix.
-  for(Index k = 0; k < size; ++k) m_p.coeffRef(k) = k;
-  for(Index k = size-1; k >= 0; --k) {
-    std::swap(m_p.coeffRef(k), m_p.coeffRef(m_transpositions.coeff(k)));
-  }
-    
   m_isInitialized = true;
   return *this;
 }
@@ -342,8 +330,21 @@ struct ei_solve_retval<LDLT<_MatrixType,_UpLo>, Rhs>
 
   template<typename Dest> void evalTo(Dest& dst) const
   {
-    dst = rhs();
-    dec().solveInPlace(dst);
+    ei_assert(rhs().rows() == dec().matrixLDLT().rows());
+    // dst = P b
+    dst = dec().transpositionsP() * rhs();
+
+    // dst = L^-1 (P b)
+    dec().matrixL().solveInPlace(dst);
+
+    // dst = D^-1 (L^-1 P b)
+    dst = dec().vectorD().asDiagonal().inverse() * dst;
+
+    // dst = L^-T (D^-1 L^-1 P b)
+    dec().matrixU().solveInPlace(dst);
+
+    // dst = P^-1 (L^-T D^-1 L^-1 P b) = A^-1 b
+    dst = dec().transpositionsP().transpose() * dst;
   }
 };
 
@@ -366,21 +367,7 @@ bool LDLT<MatrixType,_UpLo>::solveInPlace(MatrixBase<Derived> &bAndX) const
   const Index size = m_matrix.rows();
   ei_assert(size == bAndX.rows());
 
-  // z = P b
-  for(Index i = 0; i < size; ++i) bAndX.row(m_transpositions.coeff(i)).swap(bAndX.row(i));
-
-  // y = L^-1 z
-  //matrixL().solveInPlace(bAndX);
-  matrixL().solveInPlace(bAndX);
-
-  // w = D^-1 y
-  bAndX = m_matrix.diagonal().asDiagonal().inverse() * bAndX;
-
-  // u = L^-T w
-  matrixU().solveInPlace(bAndX);
-
-  // x = P^T u
-  for (Index i = size-1; i >= 0; --i) bAndX.row(m_transpositions.coeff(i)).swap(bAndX.row(i));
+  bAndX = this->solve(bAndX);
 
   return true;
 }
@@ -394,10 +381,10 @@ MatrixType LDLT<MatrixType,_UpLo>::reconstructedMatrix() const
   ei_assert(m_isInitialized && "LDLT is not initialized.");
   const Index size = m_matrix.rows();
   MatrixType res(size,size);
-  res.setIdentity();
 
-  // PI
-  for(Index i = 0; i < size; ++i) res.row(m_transpositions.coeff(i)).swap(res.row(i));
+  // P
+  res.setIdentity();
+  res = transpositionsP() * res;
   // L^* P
   res = matrixU() * res;
   // D(L^*P)
@@ -405,7 +392,7 @@ MatrixType LDLT<MatrixType,_UpLo>::reconstructedMatrix() const
   // L(DL^*P)
   res = matrixL() * res;
   // P^T (LDL^*P)
-  for (Index i = size-1; i >= 0; --i) res.row(m_transpositions.coeff(i)).swap(res.row(i));
+  res = transpositionsP().transpose() * res;
 
   return res;
 }
diff --git a/Eigen/src/Core/PermutationMatrix.h b/Eigen/src/Core/PermutationMatrix.h
index 46884dc3f..d3e36c73a 100644
--- a/Eigen/src/Core/PermutationMatrix.h
+++ b/Eigen/src/Core/PermutationMatrix.h
@@ -2,6 +2,7 @@
 // for linear algebra.
 //
 // Copyright (C) 2009 Benoit Jacob <jacob.benoit.1@gmail.com>
+// Copyright (C) 2009 Gael Guennebaud <g.gael@free.fr>
 //
 // Eigen is free software; you can redistribute it and/or
 // modify it under the terms of the GNU Lesser General Public
diff --git a/Eigen/src/Core/Transpositions.h b/Eigen/src/Core/Transpositions.h
new file mode 100644
index 000000000..39cb24fd7
--- /dev/null
+++ b/Eigen/src/Core/Transpositions.h
@@ -0,0 +1,285 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra.
+//
+// Copyright (C) 2010 Gael Guennebaud <g.gael@free.fr>
+//
+// Eigen is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 3 of the License, or (at your option) any later version.
+//
+// Alternatively, you can redistribute it and/or
+// modify it under the terms of the GNU General Public License as
+// published by the Free Software Foundation; either version 2 of
+// the License, or (at your option) any later version.
+//
+// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
+// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License and a copy of the GNU General Public License along with
+// Eigen. If not, see <http://www.gnu.org/licenses/>.
+
+#ifndef EIGEN_TRANSPOSITIONS_H
+#define EIGEN_TRANSPOSITIONS_H
+
+/** \class Transpositions
+  *
+  * \brief Represents a sequence of transpositions (row/column interchange)
+  *
+  * \param SizeAtCompileTime the number of transpositions, or Dynamic
+  * \param MaxSizeAtCompileTime the maximum number of transpositions, or Dynamic. This optional parameter defaults to SizeAtCompileTime. Most of the time, you should not have to specify it.
+  *
+  * This class represents a permutation transformation as a sequence of \em n transpositions
+  * \f$[T_{n-1} \ldots T_{i} \ldots T_{0}]\f$. It is internally stored as a vector of integers \c indices.
+  * Each transposition \f$ T_{i} \f$ applied on the left of a matrix (\f$ T_{i} M\f$) interchanges
+  * the rows \c i and \c indices[i] of the matrix \c M.
+  * A transposition applied on the right (e.g., \f$ M T_{i}\f$) yields a column interchange.
+  *
+  * Compared to the class PermutationMatrix, such a sequence of transpositions is what is
+  * computed during a decomposition with pivoting, and it is faster when applying the permutation in-place.
+  * 
+  * To apply a sequence of transpositions to a matrix, simply use the operator * as in the following example:
+  * \code
+  * Transpositions tr;
+  * MatrixXf mat;
+  * mat = tr * mat;
+  * \endcode
+  * In this example, we detect that the matrix appears on both side, and so the transpositions
+  * are applied in-place without any temporary or extra copy.
+  *
+  * \sa class PermutationMatrix
+  */
+template<int SizeAtCompileTime, int MaxSizeAtCompileTime = SizeAtCompileTime> class Transpositions;
+template<typename TranspositionType, typename MatrixType, int Side, bool Transposed=false> struct ei_transposition_matrix_product_retval;
+
+template<int SizeAtCompileTime, int MaxSizeAtCompileTime>
+class Transpositions
+{
+  public:
+
+    typedef Matrix<DenseIndex, SizeAtCompileTime, 1, 0, MaxSizeAtCompileTime, 1> IndicesType;
+    typedef typename IndicesType::Index Index;
+
+    inline Transpositions() {}
+
+    /** Copy constructor. */
+    template<int OtherSize, int OtherMaxSize>
+    inline Transpositions(const Transpositions<OtherSize, OtherMaxSize>& other)
+      : m_indices(other.indices()) {}
+
+    #ifndef EIGEN_PARSED_BY_DOXYGEN
+    /** Standard copy constructor. Defined only to prevent a default copy constructor
+      * from hiding the other templated constructor */
+    inline Transpositions(const Transpositions& other) : m_indices(other.indices()) {}
+    #endif
+
+    /** Generic constructor from expression of the transposition indices. */
+    template<typename Other>
+    explicit inline Transpositions(const MatrixBase<Other>& indices) : m_indices(indices)
+    {}
+
+    /** Copies the \a other transpositions into \c *this */
+    template<int OtherSize, int OtherMaxSize>
+    Transpositions& operator=(const Transpositions<OtherSize, OtherMaxSize>& other)
+    {
+      m_indices = other.indices();
+      return *this;
+    }
+
+    #ifndef EIGEN_PARSED_BY_DOXYGEN
+    /** This is a special case of the templated operator=. Its purpose is to
+      * prevent a default operator= from hiding the templated operator=.
+      */
+    Transpositions& operator=(const Transpositions& other)
+    {
+      m_indices = other.m_indices;
+      return *this;
+    }
+    #endif
+
+    /** Constructs an uninitialized permutation matrix of given size.
+      */
+    inline Transpositions(Index size) : m_indices(size)
+    {}
+
+    /** \returns the number of transpositions */
+    inline Index size() const { return m_indices.size(); }
+
+    inline const Index& coeff(Index i) const { return m_indices.coeff(i); }
+    inline Index& coeffRef(Index i) { return m_indices.coeffRef(i); }
+    inline const Index& operator()(Index i) const { return m_indices(i); }
+    inline Index& operator()(Index i) { return m_indices(i); }
+
+    /** const version of indices(). */
+    const IndicesType& indices() const { return m_indices; }
+    /** \returns a reference to the stored array representing the transpositions. */
+    IndicesType& indices() { return m_indices; }
+
+    /** Resizes to given size. */
+    inline void resize(int size)
+    {
+      m_indices.resize(size);
+    }
+
+    /** Sets \c *this to represents an identity transformation */
+    void setIdentity()
+    {
+      for(int i = 0; i < m_indices.size(); ++i)
+        m_indices.coeffRef(i) = i;
+    }
+
+    // FIXME: do we want such methods ?
+    // might be usefull when the target matrix expression is complex, e.g.:
+    // object.matrix().block(..,..,..,..) = trans * object.matrix().block(..,..,..,..);
+    /*
+    template<typename MatrixType>
+    void applyForwardToRows(MatrixType& mat) const
+    {
+      for(Index k=0 ; k<size() ; ++k)
+        if(m_indices(k)!=k)
+          mat.row(k).swap(mat.row(m_indices(k)));
+    }
+
+    template<typename MatrixType>
+    void applyBackwardToRows(MatrixType& mat) const
+    {
+      for(Index k=size()-1 ; k>=0 ; --k)
+        if(m_indices(k)!=k)
+          mat.row(k).swap(mat.row(m_indices(k)));
+    }
+    */
+
+    /** \returns the inverse transformation */
+    inline Transpose<Transpositions> inverse() const
+    { return *this; }
+
+    /** \returns the tranpose transformation */
+    inline Transpose<Transpositions> transpose() const
+    { return *this; }
+
+#ifndef EIGEN_PARSED_BY_DOXYGEN
+    template<int OtherSize, int OtherMaxSize>
+    Transpositions(const Transpose<Transpositions<OtherSize,OtherMaxSize> >& other)
+      : m_indices(other.size())
+    {
+      Index n = size();
+      Index j = size-1;
+      for(Index i=0; i<n;++i,--j)
+        m_indices.coeffRef(j) = other.nestedTranspositions().indices().coeff(i);
+    }
+#endif
+
+  protected:
+
+    IndicesType m_indices;
+};
+
+/** \returns the \a matrix with the \a transpositions applied to the columns.
+  */
+template<typename Derived, int SizeAtCompileTime, int MaxSizeAtCompileTime>
+inline const ei_transposition_matrix_product_retval<Transpositions<SizeAtCompileTime, MaxSizeAtCompileTime>, Derived, OnTheRight>
+operator*(const MatrixBase<Derived>& matrix,
+          const Transpositions<SizeAtCompileTime, MaxSizeAtCompileTime> &transpositions)
+{
+  return ei_transposition_matrix_product_retval
+           <Transpositions<SizeAtCompileTime, MaxSizeAtCompileTime>, Derived, OnTheRight>
+           (transpositions, matrix.derived());
+}
+
+/** \returns the \a matrix with the \a transpositions applied to the rows.
+  */
+template<typename Derived, int SizeAtCompileTime, int MaxSizeAtCompileTime>
+inline const ei_transposition_matrix_product_retval
+               <Transpositions<SizeAtCompileTime, MaxSizeAtCompileTime>, Derived, OnTheLeft>
+operator*(const Transpositions<SizeAtCompileTime, MaxSizeAtCompileTime> &transpositions,
+          const MatrixBase<Derived>& matrix)
+{
+  return ei_transposition_matrix_product_retval
+           <Transpositions<SizeAtCompileTime, MaxSizeAtCompileTime>, Derived, OnTheLeft>
+           (transpositions, matrix.derived());
+}
+
+template<typename TranspositionType, typename MatrixType, int Side, bool Transposed>
+struct ei_traits<ei_transposition_matrix_product_retval<TranspositionType, MatrixType, Side, Transposed> >
+{
+  typedef typename MatrixType::PlainObject ReturnType;
+};
+
+template<typename TranspositionType, typename MatrixType, int Side, bool Transposed>
+struct ei_transposition_matrix_product_retval
+ : public ReturnByValue<ei_transposition_matrix_product_retval<TranspositionType, MatrixType, Side, Transposed> >
+{
+    typedef typename ei_cleantype<typename MatrixType::Nested>::type MatrixTypeNestedCleaned;
+    typedef typename TranspositionType::Index Index;
+
+    ei_transposition_matrix_product_retval(const TranspositionType& tr, const MatrixType& matrix)
+      : m_transpositions(tr), m_matrix(matrix)
+    {}
+
+    inline int rows() const { return m_matrix.rows(); }
+    inline int cols() const { return m_matrix.cols(); }
+
+    template<typename Dest> inline void evalTo(Dest& dst) const
+    {
+      const int size = m_transpositions.size();
+      Index j = 0;
+
+      if(!(ei_is_same_type<MatrixTypeNestedCleaned,Dest>::ret && ei_extract_data(dst) == ei_extract_data(m_matrix)))
+        dst = m_matrix;
+
+      for(int k=(Transposed?size-1:0) ; Transposed?k>=0:k<size ; Transposed?--k:++k)
+        if((j=m_transpositions.coeff(k))!=k)
+        {
+          if(Side==OnTheLeft)
+            dst.row(k).swap(dst.row(j));
+          else if(Side==OnTheRight)
+            dst.col(k).swap(dst.col(j));
+        }
+    }
+
+  protected:
+    const TranspositionType& m_transpositions;
+    const typename MatrixType::Nested m_matrix;
+};
+
+/* Template partial specialization for transposed/inverse transpositions */
+
+template<int SizeAtCompileTime, int MaxSizeAtCompileTime>
+class Transpose<Transpositions<SizeAtCompileTime, MaxSizeAtCompileTime> >
+{
+    typedef Transpositions<SizeAtCompileTime, MaxSizeAtCompileTime> TranspositionType;
+    typedef typename TranspositionType::IndicesType IndicesType;
+  public:
+
+    Transpose(const TranspositionType& t) : m_transpositions(t) {}
+
+    inline int size() const { return m_transpositions.size(); }
+
+    /** \returns the \a matrix with the inverse transpositions applied to the columns.
+      */
+    template<typename Derived> friend
+    inline const ei_transposition_matrix_product_retval<TranspositionType, Derived, OnTheRight, true>
+    operator*(const MatrixBase<Derived>& matrix, const Transpose& trt)
+    {
+      return ei_transposition_matrix_product_retval<TranspositionType, Derived, OnTheRight, true>(trt.m_transpositions, matrix.derived());
+    }
+
+    /** \returns the \a matrix with the inverse transpositions applied to the rows.
+      */
+    template<typename Derived>
+    inline const ei_transposition_matrix_product_retval<TranspositionType, Derived, OnTheLeft, true>
+    operator*(const MatrixBase<Derived>& matrix) const
+    {
+      return ei_transposition_matrix_product_retval<TranspositionType, Derived, OnTheLeft, true>(m_transpositions, matrix.derived());
+    }
+
+    const TranspositionType& nestedTranspositions() const { return m_transpositions; }
+
+  protected:
+    const TranspositionType& m_transpositions;
+};
+
+#endif // EIGEN_TRANSPOSITIONS_H
diff --git a/test/cholesky.cpp b/test/cholesky.cpp
index d403af7ba..feb7be289 100644
--- a/test/cholesky.cpp
+++ b/test/cholesky.cpp
@@ -26,10 +26,21 @@
 #define EIGEN_NO_ASSERTION_CHECKING
 #endif
 
+static int nb_temporaries;
+
+#define EIGEN_DEBUG_MATRIX_CTOR { if(size!=0) nb_temporaries++; }
+
 #include "main.h"
 #include <Eigen/Cholesky>
 #include <Eigen/QR>
 
+#define VERIFY_EVALUATION_COUNT(XPR,N) {\
+    nb_temporaries = 0; \
+    XPR; \
+    if(nb_temporaries!=N) std::cerr << "nb_temporaries == " << nb_temporaries << "\n"; \
+    VERIFY( (#XPR) && nb_temporaries==N ); \
+  }
+
 #ifdef HAS_GSL
 #include "gsl_helper.h"
 #endif
@@ -131,6 +142,21 @@ template<typename MatrixType> void cholesky(const MatrixType& m)
     VERIFY_IS_APPROX(symm * vecX, vecB);
     matX = ldltup.solve(matB);
     VERIFY_IS_APPROX(symm * matX, matB);
+
+    if(MatrixType::RowsAtCompileTime==Dynamic)
+    {
+      // note : each inplace permutation requires a small temporary vector (mask)
+
+      // check inplace solve
+      matX = matB;
+      VERIFY_EVALUATION_COUNT(matX = ldltlo.solve(matX), 0);
+      VERIFY_IS_APPROX(matX, ldltlo.solve(matB).eval());
+
+
+      matX = matB;
+      VERIFY_EVALUATION_COUNT(matX = ldltup.solve(matX), 0);
+      VERIFY_IS_APPROX(matX, ldltup.solve(matB).eval());
+    }
   }
 
 }
@@ -141,6 +167,7 @@ template<typename MatrixType> void cholesky_verify_assert()
 
   LLT<MatrixType> llt;
   VERIFY_RAISES_ASSERT(llt.matrixL())
+  VERIFY_RAISES_ASSERT(llt.matrixU())
   VERIFY_RAISES_ASSERT(llt.solve(tmp))
   VERIFY_RAISES_ASSERT(llt.solveInPlace(&tmp))