Big API change in Cholesky module:

 * rename Cholesky to LLT
 * rename CholeskyWithoutSquareRoot to LDLT
 * rename MatrixBase::cholesky() to llt()
 * rename MatrixBase::choleskyNoSqrt() to ldlt()
 * make {LLT,LDLT}::solve() API consistent with other modules

Note that we are going to keep a source compatibility untill the next beta release.
E.g., the "old" Cholesky* classes, etc are still available for some time.
To be clear, Eigen beta2 should be (hopefully) source compatible with beta1,
and so beta2 will contain all the deprecated API of beta1. Those features marked
as deprecated will be removed in beta3 (or in the final 2.0 if there is no beta 3 !).

Also includes various updated in sparse Cholesky.

23 files changed, 609 insertions, 184 deletions

diff --git a/Eigen/Cholesky b/Eigen/Cholesky
index f5bcec52d..ddde61a31 100644
--- a/Eigen/Cholesky
+++ b/Eigen/Cholesky
@@ -17,8 +17,8 @@ namespace Eigen {
 /** \defgroup Cholesky_Module Cholesky module
   * This module provides two variants of the Cholesky decomposition for selfadjoint (hermitian) matrices.
   * Those decompositions are accessible via the following MatrixBase methods:
-  *  - MatrixBase::cholesky(),
-  *  - MatrixBase::choleskyNoSqrt()
+  *  - MatrixBase::llt(),
+  *  - MatrixBase::ldlt()
   *
   * \code
   * #include <Eigen/Cholesky>
@@ -27,6 +27,8 @@ namespace Eigen {
 
 #include "src/Array/CwiseOperators.h"
 #include "src/Array/Functors.h"
+#include "src/Cholesky/LLT.h"
+#include "src/Cholesky/LDLT.h"
 #include "src/Cholesky/Cholesky.h"
 #include "src/Cholesky/CholeskyWithoutSquareRoot.h"
 
diff --git a/Eigen/src/Cholesky/Cholesky.h b/Eigen/src/Cholesky/Cholesky.h
index b94fea8dc..ada413b33 100644
--- a/Eigen/src/Cholesky/Cholesky.h
+++ b/Eigen/src/Cholesky/Cholesky.h
@@ -29,22 +29,7 @@
   *
   * \class Cholesky
   *
-  * \brief Standard Cholesky decomposition of a matrix and associated features
-  *
-  * \param MatrixType the type of the matrix of which we are computing the Cholesky decomposition
-  *
-  * This class performs a standard Cholesky decomposition of a symmetric, positive definite
-  * matrix A such that A = LL^* = U^*U, where L is lower triangular.
-  *
-  * While the Cholesky decomposition is particularly useful to solve selfadjoint problems like  D^*D x = b,
-  * for that purpose, we recommend the Cholesky decomposition without square root which is more stable
-  * and even faster. Nevertheless, this standard Cholesky decomposition remains useful in many other
-  * situations like generalised eigen problems with hermitian matrices.
-  *
-  * Note that during the decomposition, only the upper triangular part of A is considered. Therefore,
-  * the strict lower part does not have to store correct values.
-  *
-  * \sa MatrixBase::cholesky(), class CholeskyWithoutSquareRoot
+  * \deprecated this class has been renamed LLT
   */
 template<typename MatrixType> class Cholesky
 {
@@ -72,7 +57,10 @@ template<typename MatrixType> class Cholesky
     inline bool isPositiveDefinite(void) const { return m_isPositiveDefinite; }
 
     template<typename Derived>
-    typename Derived::Eval solve(const MatrixBase<Derived> &b) const;
+    typename Derived::Eval solve(const MatrixBase<Derived> &b) const EIGEN_DEPRECATED;
+
+    template<typename RhsDerived, typename ResDerived>
+    bool solve(const MatrixBase<RhsDerived> &b, MatrixBase<ResDerived> *result) const;
 
     template<typename Derived>
     bool solveInPlace(MatrixBase<Derived> &bAndX) const;
@@ -128,16 +116,7 @@ void Cholesky<MatrixType>::compute(const MatrixType& a)
   }
 }
 
-/** \returns the solution of \f$ A x = b \f$ using the current decomposition of A.
-  * In other words, it returns \f$ A^{-1} b \f$ computing
-  * \f$ {L^{*}}^{-1} L^{-1} b \f$ from right to left.
-  * \param b the column vector \f$ b \f$, which can also be a matrix.
-  *
-  * Example: \include Cholesky_solve.cpp
-  * Output: \verbinclude Cholesky_solve.out
-  *
-  * \sa MatrixBase::cholesky(), CholeskyWithoutSquareRoot::solve()
-  */
+/** \deprecated */
 template<typename MatrixType>
 template<typename Derived>
 typename Derived::Eval Cholesky<MatrixType>::solve(const MatrixBase<Derived> &b) const
@@ -147,7 +126,6 @@ typename Derived::Eval Cholesky<MatrixType>::solve(const MatrixBase<Derived> &b)
   typename ei_eval_to_column_major<Derived>::type x(b);
   solveInPlace(x);
   return x;
-  //return m_matrix.adjoint().template part<Upper>().solveTriangular(matrixL().solveTriangular(b));
 }
 
 /** Computes the solution x of \f$ A x = b \f$ using the current decomposition of A.
@@ -162,7 +140,25 @@ typename Derived::Eval Cholesky<MatrixType>::solve(const MatrixBase<Derived> &b)
   * Example: \include Cholesky_solve.cpp
   * Output: \verbinclude Cholesky_solve.out
   *
-  * \sa MatrixBase::cholesky(), Cholesky::solve()
+  * \sa MatrixBase::cholesky(), Cholesky::solveInPlace()
+  */
+template<typename MatrixType>
+template<typename RhsDerived, typename ResDerived>
+bool Cholesky<MatrixType>::solve(const MatrixBase<RhsDerived> &b, MatrixBase<ResDerived> *result) const
+{
+  const int size = m_matrix.rows();
+  ei_assert(size==b.rows() && "Cholesky::solve(): invalid number of rows of the right hand side matrix b");
+  return solveInPlace((*result) = b);
+}
+
+/** This is the \em in-place version of solve().
+  *
+  * \param bAndX represents both the right-hand side matrix b and result x.
+  *
+  * This version avoids a copy when the right hand side matrix b is not
+  * needed anymore.
+  *
+  * \sa Cholesky::solve(), MatrixBase::cholesky()
   */
 template<typename MatrixType>
 template<typename Derived>
@@ -178,7 +174,7 @@ bool Cholesky<MatrixType>::solveInPlace(MatrixBase<Derived> &bAndX) const
 }
 
 /** \cholesky_module
-  * \returns the Cholesky decomposition of \c *this
+  * \deprecated has been renamed llt()
   */
 template<typename Derived>
 inline const Cholesky<typename MatrixBase<Derived>::EvalType>
diff --git a/Eigen/src/Cholesky/CholeskyWithoutSquareRoot.h b/Eigen/src/Cholesky/CholeskyWithoutSquareRoot.h
index fd111fb1f..af44634a0 100644
--- a/Eigen/src/Cholesky/CholeskyWithoutSquareRoot.h
+++ b/Eigen/src/Cholesky/CholeskyWithoutSquareRoot.h
@@ -25,25 +25,11 @@
 #ifndef EIGEN_CHOLESKY_WITHOUT_SQUARE_ROOT_H
 #define EIGEN_CHOLESKY_WITHOUT_SQUARE_ROOT_H
 
-/** \ingroup Cholesky_Module
+/** \deprecated \ingroup Cholesky_Module
   *
   * \class CholeskyWithoutSquareRoot
   *
-  * \brief Robust Cholesky decomposition of a matrix and associated features
-  *
-  * \param MatrixType the type of the matrix of which we are computing the Cholesky decomposition
-  *
-  * This class performs a Cholesky decomposition without square root of a symmetric, positive definite
-  * matrix A such that A = L D L^* = U^* D U, where L is lower triangular with a unit diagonal
-  * and D is a diagonal matrix.
-  *
-  * Compared to a standard Cholesky decomposition, avoiding the square roots allows for faster and more
-  * stable computation.
-  *
-  * Note that during the decomposition, only the upper triangular part of A is considered. Therefore,
-  * the strict lower part does not have to store correct values.
-  *
-  * \sa MatrixBase::choleskyNoSqrt(), class Cholesky
+  * \deprecated this class has been renamed LDLT
   */
 template<typename MatrixType> class CholeskyWithoutSquareRoot
 {
@@ -69,7 +55,10 @@ template<typename MatrixType> class CholeskyWithoutSquareRoot
     inline bool isPositiveDefinite(void) const { return m_isPositiveDefinite; }
 
     template<typename Derived>
-    typename Derived::Eval solve(const MatrixBase<Derived> &b) const;
+    typename Derived::Eval solve(const MatrixBase<Derived> &b) const EIGEN_DEPRECATED;
+
+    template<typename RhsDerived, typename ResDerived>
+    bool solve(const MatrixBase<RhsDerived> &b, MatrixBase<ResDerived> *result) const;
 
 		template<typename Derived>
     bool solveInPlace(MatrixBase<Derived> &bAndX) const;
@@ -141,15 +130,7 @@ void CholeskyWithoutSquareRoot<MatrixType>::compute(const MatrixType& a)
   }
 }
 
-/** \returns the solution of \f$ A x = b \f$ using the current decomposition of A.
-  * In other words, it returns \f$ A^{-1} b \f$ computing
-  * \f$ {L^{*}}^{-1} D^{-1} L^{-1} b \f$ from right to left.
-  * \param b the column vector \f$ b \f$, which can also be a matrix.
-  *
-  * See Cholesky::solve() for a example.
-  *
-  * \sa CholeskyWithoutSquareRoot::solveInPlace(), MatrixBase::choleskyNoSqrt()
-  */
+/** \deprecated */
 template<typename MatrixType>
 template<typename Derived>
 typename Derived::Eval CholeskyWithoutSquareRoot<MatrixType>::solve(const MatrixBase<Derived> &b) const
@@ -173,10 +154,30 @@ typename Derived::Eval CholeskyWithoutSquareRoot<MatrixType>::solve(const Matrix
   * \f$ {L^{*}}^{-1} D^{-1} L^{-1} b \f$ from right to left.
   * \param bAndX stores both the matrix \f$ b \f$ and the result \f$ x \f$
   *
-  * Example: \include Cholesky_solve.cpp
-  * Output: \verbinclude Cholesky_solve.out
+  * Example: \include CholeskyCholeskyWithoutSquareRoot_solve.cpp
+  * Output: \verbinclude CholeskyCholeskyWithoutSquareRoot_solve.out
+  *
+  * \sa CholeskyWithoutSquareRoot::solveInPlace(), MatrixBase::choleskyNoSqrt()
+  */
+template<typename MatrixType>
+template<typename RhsDerived, typename ResDerived>
+bool CholeskyWithoutSquareRoot<MatrixType>
+::solve(const MatrixBase<RhsDerived> &b, MatrixBase<ResDerived> *result) const
+{
+  const int size = m_matrix.rows();
+  ei_assert(size==b.rows() && "Cholesky::solve(): invalid number of rows of the right hand side matrix b");
+  *result = b;
+  return solveInPlace(*result);
+}
+
+/** This is the \em in-place version of solve().
+  *
+  * \param bAndX represents both the right-hand side matrix b and result x.
+  *
+  * This version avoids a copy when the right hand side matrix b is not
+  * needed anymore.
   *
-  * \sa MatrixBase::cholesky(), CholeskyWithoutSquareRoot::solve()
+  * \sa CholeskyWithoutSquareRoot::solve(), MatrixBase::choleskyNoSqrt()
   */
 template<typename MatrixType>
 template<typename Derived>
@@ -187,13 +188,13 @@ bool CholeskyWithoutSquareRoot<MatrixType>::solveInPlace(MatrixBase<Derived> &bA
   if (!m_isPositiveDefinite)
     return false;
   matrixL().solveTriangularInPlace(bAndX);
-	bAndX *= m_matrix.cwise().inverse().template part<Diagonal>();
+  bAndX = (m_matrix.cwise().inverse().template part<Diagonal>() * bAndX).lazy();
   m_matrix.adjoint().template part<UnitUpper>().solveTriangularInPlace(bAndX);
   return true;
 }
 
-/** \cholesky_module
-  * \returns the Cholesky decomposition without square root of \c *this
+/** \deprecated \cholesky_module
+  * \deprecated has been renamed ldlt()
   */
 template<typename Derived>
 inline const CholeskyWithoutSquareRoot<typename MatrixBase<Derived>::EvalType>
diff --git a/Eigen/src/Cholesky/LDLT.h b/Eigen/src/Cholesky/LDLT.h
new file mode 100644
index 000000000..e70a324f6
--- /dev/null
+++ b/Eigen/src/Cholesky/LDLT.h
@@ -0,0 +1,202 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra. Eigen itself is part of the KDE project.
+//
+// Copyright (C) 2008 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_LDLT_H
+#define EIGEN_LDLT_H
+
+/** \ingroup cholesky_Module
+  *
+  * \class LDLT
+  *
+  * \brief Robust Cholesky decomposition of a matrix and associated features
+  *
+  * \param MatrixType the type of the matrix of which we are computing the LDL^T Cholesky decomposition
+  *
+  * This class performs a Cholesky decomposition without square root of a symmetric, positive definite
+  * matrix A such that A = L D L^* = U^* D U, where L is lower triangular with a unit diagonal
+  * and D is a diagonal matrix.
+  *
+  * Compared to a standard Cholesky decomposition, avoiding the square roots allows for faster and more
+  * stable computation.
+  *
+  * Note that during the decomposition, only the upper triangular part of A is considered. Therefore,
+  * the strict lower part does not have to store correct values.
+  *
+  * \sa MatrixBase::ldlt(), class LLT
+  */
+template<typename MatrixType> class LDLT
+{
+  public:
+
+    typedef typename MatrixType::Scalar Scalar;
+    typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar;
+    typedef Matrix<Scalar, MatrixType::ColsAtCompileTime, 1> VectorType;
+
+    LDLT(const MatrixType& matrix)
+      : m_matrix(matrix.rows(), matrix.cols())
+    {
+      compute(matrix);
+    }
+
+    /** \returns the lower triangular matrix L */
+    inline Part<MatrixType, UnitLower> matrixL(void) const { return m_matrix; }
+
+    /** \returns the coefficients of the diagonal matrix D */
+    inline DiagonalCoeffs<MatrixType> vectorD(void) const { return m_matrix.diagonal(); }
+
+    /** \returns true if the matrix is positive definite */
+    inline bool isPositiveDefinite(void) const { return m_isPositiveDefinite; }
+
+    template<typename RhsDerived, typename ResDerived>
+    bool solve(const MatrixBase<RhsDerived> &b, MatrixBase<ResDerived> *result) const;
+
+    template<typename Derived>
+    bool solveInPlace(MatrixBase<Derived> &bAndX) const;
+
+    void compute(const MatrixType& matrix);
+
+  protected:
+    /** \internal
+      * Used to compute and store the cholesky decomposition A = L D L^* = U^* D U.
+      * The strict upper part is used during the decomposition, the strict lower
+      * part correspond to the coefficients of L (its diagonal is equal to 1 and
+      * is not stored), and the diagonal entries correspond to D.
+      */
+    MatrixType m_matrix;
+
+    bool m_isPositiveDefinite;
+};
+
+/** Compute / recompute the LLT decomposition A = L D L^* = U^* D U of \a matrix
+  */
+template<typename MatrixType>
+void LDLT<MatrixType>::compute(const MatrixType& a)
+{
+  assert(a.rows()==a.cols());
+  const int size = a.rows();
+  m_matrix.resize(size, size);
+  m_isPositiveDefinite = true;
+  const RealScalar eps = ei_sqrt(precision<Scalar>());
+
+  if (size<=1)
+  {
+    m_matrix = a;
+    return;
+  }
+
+  // Let's preallocate a temporay vector to evaluate the matrix-vector product into it.
+  // Unlike the standard LLT decomposition, here we cannot evaluate it to the destination
+  // matrix because it a sub-row which is not compatible suitable for efficient packet evaluation.
+  // (at least if we assume the matrix is col-major)
+  Matrix<Scalar,MatrixType::RowsAtCompileTime,1> _temporary(size);
+
+  // Note that, in this algorithm the rows of the strict upper part of m_matrix is used to store
+  // column vector, thus the strange .conjugate() and .transpose()...
+
+  m_matrix.row(0) = a.row(0).conjugate();
+  m_matrix.col(0).end(size-1) = m_matrix.row(0).end(size-1) / m_matrix.coeff(0,0);
+  for (int j = 1; j < size; ++j)
+  {
+    RealScalar tmp = ei_real(a.coeff(j,j) - (m_matrix.row(j).start(j) * m_matrix.col(j).start(j).conjugate()).coeff(0,0));
+    m_matrix.coeffRef(j,j) = tmp;
+
+    if (tmp < eps)
+    {
+      m_isPositiveDefinite = false;
+      return;
+    }
+
+    int endSize = size-j-1;
+    if (endSize>0)
+    {
+      _temporary.end(endSize) = ( m_matrix.block(j+1,0, endSize, j)
+                                  * m_matrix.col(j).start(j).conjugate() ).lazy();
+
+      m_matrix.row(j).end(endSize) = a.row(j).end(endSize).conjugate()
+                                   - _temporary.end(endSize).transpose();
+
+      m_matrix.col(j).end(endSize) = m_matrix.row(j).end(endSize) / tmp;
+    }
+  }
+}
+
+/** Computes the solution x of \f$ A x = b \f$ using the current decomposition of A.
+  * The result is stored in \a bAndx
+  *
+  * \returns true in case of success, false otherwise.
+  *
+  * In other words, it computes \f$ b = A^{-1} b \f$ with
+  * \f$ {L^{*}}^{-1} D^{-1} L^{-1} b \f$ from right to left.
+  * \param bAndX stores both the matrix \f$ b \f$ and the result \f$ x \f$
+  *
+  * Example: \include LLTLDLT_solve.cpp
+  * Output: \verbinclude LLTLDLT_solve.out
+  *
+  * \sa LDLT::solveInPlace(), MatrixBase::ldlt()
+  */
+template<typename MatrixType>
+template<typename RhsDerived, typename ResDerived>
+bool LDLT<MatrixType>
+::solve(const MatrixBase<RhsDerived> &b, MatrixBase<ResDerived> *result) const
+{
+  const int size = m_matrix.rows();
+  ei_assert(size==b.rows() && "LLT::solve(): invalid number of rows of the right hand side matrix b");
+  *result = b;
+  return solveInPlace(*result);
+}
+
+/** This is the \em in-place version of solve().
+  *
+  * \param bAndX represents both the right-hand side matrix b and result x.
+  *
+  * This version avoids a copy when the right hand side matrix b is not
+  * needed anymore.
+  *
+  * \sa LDLT::solve(), MatrixBase::ldlt()
+  */
+template<typename MatrixType>
+template<typename Derived>
+bool LDLT<MatrixType>::solveInPlace(MatrixBase<Derived> &bAndX) const
+{
+  const int size = m_matrix.rows();
+  ei_assert(size==bAndX.rows());
+  if (!m_isPositiveDefinite)
+    return false;
+  matrixL().solveTriangularInPlace(bAndX);
+  bAndX = (m_matrix.cwise().inverse().template part<Diagonal>() * bAndX).lazy();
+  m_matrix.adjoint().template part<UnitUpper>().solveTriangularInPlace(bAndX);
+  return true;
+}
+
+/** \cholesky_module
+  * \returns the Cholesky decomposition without square root of \c *this
+  */
+template<typename Derived>
+inline const LDLT<typename MatrixBase<Derived>::EvalType>
+MatrixBase<Derived>::ldlt() const
+{
+  return derived();
+}
+
+#endif // EIGEN_LDLT_H
diff --git a/Eigen/src/Cholesky/LLT.h b/Eigen/src/Cholesky/LLT.h
new file mode 100644
index 000000000..8d4a1a752
--- /dev/null
+++ b/Eigen/src/Cholesky/LLT.h
@@ -0,0 +1,186 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra. Eigen itself is part of the KDE project.
+//
+// Copyright (C) 2008 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_LLT_H
+#define EIGEN_LLT_H
+
+/** \ingroup cholesky_Module
+  *
+  * \class LLT
+  *
+  * \brief Standard Cholesky decomposition (LL^T) of a matrix and associated features
+  *
+  * \param MatrixType the type of the matrix of which we are computing the LL^T Cholesky decomposition
+  *
+  * This class performs a LL^T Cholesky decomposition of a symmetric, positive definite
+  * matrix A such that A = LL^* = U^*U, where L is lower triangular.
+  *
+  * While the Cholesky decomposition is particularly useful to solve selfadjoint problems like  D^*D x = b,
+  * for that purpose, we recommend the Cholesky decomposition without square root which is more stable
+  * and even faster. Nevertheless, this standard Cholesky decomposition remains useful in many other
+  * situations like generalised eigen problems with hermitian matrices.
+  *
+  * Note that during the decomposition, only the upper triangular part of A is considered. Therefore,
+  * the strict lower part does not have to store correct values.
+  *
+  * \sa MatrixBase::llt(), class LDLT
+  */
+template<typename MatrixType> class LLT
+{
+  private:
+    typedef typename MatrixType::Scalar Scalar;
+    typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar;
+    typedef Matrix<Scalar, MatrixType::ColsAtCompileTime, 1> VectorType;
+
+    enum {
+      PacketSize = ei_packet_traits<Scalar>::size,
+      AlignmentMask = int(PacketSize)-1
+    };
+
+  public:
+
+    LLT(const MatrixType& matrix)
+      : m_matrix(matrix.rows(), matrix.cols())
+    {
+      compute(matrix);
+    }
+
+    inline Part<MatrixType, Lower> matrixL(void) const { return m_matrix; }
+
+    /** \returns true if the matrix is positive definite */
+    inline bool isPositiveDefinite(void) const { return m_isPositiveDefinite; }
+
+    template<typename RhsDerived, typename ResDerived>
+    bool solve(const MatrixBase<RhsDerived> &b, MatrixBase<ResDerived> *result) const;
+
+    template<typename Derived>
+    bool solveInPlace(MatrixBase<Derived> &bAndX) const;
+
+    void compute(const MatrixType& matrix);
+
+  protected:
+    /** \internal
+      * Used to compute and store L
+      * The strict upper part is not used and even not initialized.
+      */
+    MatrixType m_matrix;
+    bool m_isPositiveDefinite;
+};
+
+/** Computes / recomputes the Cholesky decomposition A = LL^* = U^*U of \a matrix
+  */
+template<typename MatrixType>
+void LLT<MatrixType>::compute(const MatrixType& a)
+{
+  assert(a.rows()==a.cols());
+  const int size = a.rows();
+  m_matrix.resize(size, size);
+  const RealScalar eps = ei_sqrt(precision<Scalar>());
+
+  RealScalar x;
+  x = ei_real(a.coeff(0,0));
+  m_isPositiveDefinite = x > eps && ei_isMuchSmallerThan(ei_imag(a.coeff(0,0)), RealScalar(1));
+  m_matrix.coeffRef(0,0) = ei_sqrt(x);
+  m_matrix.col(0).end(size-1) = a.row(0).end(size-1).adjoint() / ei_real(m_matrix.coeff(0,0));
+  for (int j = 1; j < size; ++j)
+  {
+    Scalar tmp = ei_real(a.coeff(j,j)) - m_matrix.row(j).start(j).norm2();
+    x = ei_real(tmp);
+    if (x < eps || (!ei_isMuchSmallerThan(ei_imag(tmp), RealScalar(1))))
+    {
+      m_isPositiveDefinite = false;
+      return;
+    }
+    m_matrix.coeffRef(j,j) = x = ei_sqrt(x);
+
+    int endSize = size-j-1;
+    if (endSize>0) {
+      // Note that when all matrix columns have good alignment, then the following
+      // product is guaranteed to be optimal with respect to alignment.
+      m_matrix.col(j).end(endSize) =
+        (m_matrix.block(j+1, 0, endSize, j) * m_matrix.row(j).start(j).adjoint()).lazy();
+
+      // FIXME could use a.col instead of a.row
+      m_matrix.col(j).end(endSize) = (a.row(j).end(endSize).adjoint()
+        - m_matrix.col(j).end(endSize) ) / x;
+    }
+  }
+}
+
+/** Computes the solution x of \f$ A x = b \f$ using the current decomposition of A.
+  * The result is stored in \a bAndx
+  *
+  * \returns true in case of success, false otherwise.
+  *
+  * In other words, it computes \f$ b = A^{-1} b \f$ with
+  * \f$ {L^{*}}^{-1} L^{-1} b \f$ from right to left.
+  * \param bAndX stores both the matrix \f$ b \f$ and the result \f$ x \f$
+  *
+  * Example: \include LLT_solve.cpp
+  * Output: \verbinclude LLT_solve.out
+  *
+  * \sa LLT::solveInPlace(), MatrixBase::llt()
+  */
+template<typename MatrixType>
+template<typename RhsDerived, typename ResDerived>
+bool LLT<MatrixType>::solve(const MatrixBase<RhsDerived> &b, MatrixBase<ResDerived> *result) const
+{
+  const int size = m_matrix.rows();
+  ei_assert(size==b.rows() && "LLT::solve(): invalid number of rows of the right hand side matrix b");
+  return solveInPlace((*result) = b);
+}
+
+/** This is the \em in-place version of solve().
+  *
+  * \param bAndX represents both the right-hand side matrix b and result x.
+  *
+  * This version avoids a copy when the right hand side matrix b is not
+  * needed anymore.
+  *
+  * \sa LLT::solve(), MatrixBase::llt()
+  */
+template<typename MatrixType>
+template<typename Derived>
+bool LLT<MatrixType>::solveInPlace(MatrixBase<Derived> &bAndX) const
+{
+  const int size = m_matrix.rows();
+  ei_assert(size==bAndX.rows());
+  if (!m_isPositiveDefinite)
+    return false;
+  matrixL().solveTriangularInPlace(bAndX);
+  m_matrix.adjoint().template part<Upper>().solveTriangularInPlace(bAndX);
+  return true;
+}
+
+/** \cholesky_module
+  * \returns the LLT decomposition of \c *this
+  */
+template<typename Derived>
+inline const LLT<typename MatrixBase<Derived>::EvalType>
+MatrixBase<Derived>::llt() const
+{
+  return LLT<typename ei_eval<Derived>::type>(derived());
+}
+
+#endif // EIGEN_LLT_H
diff --git a/Eigen/src/Core/MatrixBase.h b/Eigen/src/Core/MatrixBase.h
index 944d353d8..8b0129060 100644
--- a/Eigen/src/Core/MatrixBase.h
+++ b/Eigen/src/Core/MatrixBase.h
@@ -563,6 +563,9 @@ template<typename Derived> class MatrixBase
 
 /////////// Cholesky module ///////////
 
+    const LLT<EvalType>  llt() const;
+    const LDLT<EvalType> ldlt() const;
+    // deprecated:
     const Cholesky<EvalType> cholesky() const;
     const CholeskyWithoutSquareRoot<EvalType> choleskyNoSqrt() const;
 
diff --git a/Eigen/src/Core/util/ForwardDeclarations.h b/Eigen/src/Core/util/ForwardDeclarations.h
index 55016e05d..98d15b415 100644
--- a/Eigen/src/Core/util/ForwardDeclarations.h
+++ b/Eigen/src/Core/util/ForwardDeclarations.h
@@ -102,6 +102,9 @@ template<typename ExpressionType, int Direction> class PartialRedux;
 template<typename MatrixType> class LU;
 template<typename MatrixType> class QR;
 template<typename MatrixType> class SVD;
+template<typename MatrixType> class LLT;
+template<typename MatrixType> class LDLT;
+// deprecated:
 template<typename MatrixType> class Cholesky;
 template<typename MatrixType> class CholeskyWithoutSquareRoot;
 
diff --git a/Eigen/src/Core/util/Macros.h b/Eigen/src/Core/util/Macros.h
index e3429e146..348f313d2 100644
--- a/Eigen/src/Core/util/Macros.h
+++ b/Eigen/src/Core/util/Macros.h
@@ -98,6 +98,12 @@ using Eigen::ei_cos;
 #endif
 
 #if (defined __GNUC__)
+#define EIGEN_DEPRECATED __attribute__((deprecated))
+#else
+#define EIGEN_DEPRECATED
+#endif
+
+#if (defined __GNUC__)
 #define EIGEN_ALIGN_128 __attribute__ ((aligned(16)))
 #else
 #define EIGEN_ALIGN_128
diff --git a/Eigen/src/LU/LU.h b/Eigen/src/LU/LU.h
index f7ed1b412..554d8bd3c 100644
--- a/Eigen/src/LU/LU.h
+++ b/Eigen/src/LU/LU.h
@@ -181,7 +181,8 @@ template<typename MatrixType> class LU
       * \returns true if any solution exists, false if no solution exists.
       *
       * \note If there exist more than one solution, this method will arbitrarily choose one.
-      *       If you need a complete analysis of the space of solutions, take the one solution obtained  *       by this method and add to it elements of the kernel, as determined by kernel().
+      *       If you need a complete analysis of the space of solutions, take the one solution obtained
+      *       by this method and add to it elements of the kernel, as determined by kernel().
       *
       * Example: \include LU_solve.cpp
       * Output: \verbinclude LU_solve.out
diff --git a/Eigen/src/QR/QrInstantiations.cpp b/Eigen/src/QR/QrInstantiations.cpp
index 3b2594e34..dacb05d3d 100644
--- a/Eigen/src/QR/QrInstantiations.cpp
+++ b/Eigen/src/QR/QrInstantiations.cpp
@@ -26,8 +26,6 @@
 #define EIGEN_EXTERN_INSTANTIATIONS
 #endif
 #include "../../Core"
-// commented because of -pedantic
-// #include "../../Cholesky"
 #undef EIGEN_EXTERN_INSTANTIATIONS
 
 #include "../../QR"
diff --git a/Eigen/src/QR/SelfAdjointEigenSolver.h b/Eigen/src/QR/SelfAdjointEigenSolver.h
index 9e929620c..fdb2764c4 100644
--- a/Eigen/src/QR/SelfAdjointEigenSolver.h
+++ b/Eigen/src/QR/SelfAdjointEigenSolver.h
@@ -227,7 +227,7 @@ compute(const MatrixType& matA, const MatrixType& matB, bool computeEigenvectors
   ei_assert(matA.cols()==matA.rows() && matB.rows()==matA.rows() && matB.cols()==matB.rows());
   
   // Compute the cholesky decomposition of matB = L L'
-  Cholesky<MatrixType> cholB(matB);
+  LLT<MatrixType> cholB(matB);
 
   // compute C = inv(L) A inv(L')
   MatrixType matC = matA;
diff --git a/Eigen/src/SVD/SVD.h b/Eigen/src/SVD/SVD.h
index c3f3bb235..100ca9310 100644
--- a/Eigen/src/SVD/SVD.h
+++ b/Eigen/src/SVD/SVD.h
@@ -505,7 +505,7 @@ SVD<MatrixType>& SVD<MatrixType>::sort()
 /** \returns the solution of \f$ A x = b \f$ using the current SVD decomposition of A.
   * The parts of the solution corresponding to zero singular values are ignored.
   *
-  * \sa MatrixBase::svd(), LU::solve(), Cholesky::solve()
+  * \sa MatrixBase::svd(), LU::solve(), LLT::solve()
   */
 template<typename MatrixType>
 template<typename OtherDerived, typename ResultType>
diff --git a/Eigen/src/Sparse/AmbiVector.h b/Eigen/src/Sparse/AmbiVector.h
index 0c4bd1af1..44697bceb 100644
--- a/Eigen/src/Sparse/AmbiVector.h
+++ b/Eigen/src/Sparse/AmbiVector.h
@@ -28,7 +28,7 @@
 /** \internal
   * Hybrid sparse/dense vector class designed for intensive read-write operations.
   *
-  * See BasicSparseCholesky and SparseProduct for usage examples.
+  * See BasicSparseLLT and SparseProduct for usage examples.
   */
 template<typename _Scalar> class AmbiVector
 {
@@ -269,12 +269,14 @@ class AmbiVector<_Scalar>::Iterator
 
     /** Default constructor
       * \param vec the vector on which we iterate
-      * \param nonZeroReferenceValue reference value used to prune zero coefficients.
-      * In practice, the coefficient are compared to \a nonZeroReferenceValue * precision<Scalar>().
+      * \param epsilon the minimal value used to prune zero coefficients.
+      * In practice, all coefficients having a magnitude smaller than \a epsilon
+      * are skipped.
       */
-    Iterator(const AmbiVector& vec, RealScalar nonZeroReferenceValue = RealScalar(0.1)) : m_vector(vec)
+    Iterator(const AmbiVector& vec, RealScalar epsilon = RealScalar(0.1)*precision<RealScalar>())
+      : m_vector(vec)
     {
-      m_epsilon = nonZeroReferenceValue * precision<Scalar>();
+      m_epsilon = epsilon;
       m_isDense = m_vector.m_mode==IsDense;
       if (m_isDense)
       {
diff --git a/Eigen/src/Sparse/CholmodSupport.h b/Eigen/src/Sparse/CholmodSupport.h
index ed031b264..ba2161320 100644
--- a/Eigen/src/Sparse/CholmodSupport.h
+++ b/Eigen/src/Sparse/CholmodSupport.h
@@ -99,10 +99,10 @@ SparseMatrix<Scalar,Flags> SparseMatrix<Scalar,Flags>::Map(cholmod_sparse& cm)
 }
 
 template<typename MatrixType>
-class SparseCholesky<MatrixType,Cholmod> : public SparseCholesky<MatrixType>
+class SparseLLT<MatrixType,Cholmod> : public SparseLLT<MatrixType>
 {
   protected:
-    typedef SparseCholesky<MatrixType> Base;
+    typedef SparseLLT<MatrixType> Base;
     using Base::Scalar;
     using Base::RealScalar;
     using Base::MatrixLIsDirty;
@@ -113,14 +113,20 @@ class SparseCholesky<MatrixType,Cholmod> : public SparseCholesky<MatrixType>
 
   public:
 
-    SparseCholesky(const MatrixType& matrix, int flags = 0)
+    SparseLLT(int flags = 0)
+      : Base(flags), m_cholmodFactor(0)
+    {
+      cholmod_start(&m_cholmod);
+    }
+
+    SparseLLT(const MatrixType& matrix, int flags = 0)
       : Base(matrix, flags), m_cholmodFactor(0)
     {
       cholmod_start(&m_cholmod);
       compute(matrix);
     }
 
-    ~SparseCholesky()
+    ~SparseLLT()
     {
       if (m_cholmodFactor)
         cholmod_free_factor(&m_cholmodFactor, &m_cholmod);
@@ -140,7 +146,7 @@ class SparseCholesky<MatrixType,Cholmod> : public SparseCholesky<MatrixType>
 };
 
 template<typename MatrixType>
-void SparseCholesky<MatrixType,Cholmod>::compute(const MatrixType& a)
+void SparseLLT<MatrixType,Cholmod>::compute(const MatrixType& a)
 {
   if (m_cholmodFactor)
   {
@@ -169,8 +175,8 @@ void SparseCholesky<MatrixType,Cholmod>::compute(const MatrixType& a)
 }
 
 template<typename MatrixType>
-inline const typename SparseCholesky<MatrixType>::CholMatrixType&
-SparseCholesky<MatrixType,Cholmod>::matrixL() const
+inline const typename SparseLLT<MatrixType>::CholMatrixType&
+SparseLLT<MatrixType,Cholmod>::matrixL() const
 {
   if (m_status & MatrixLIsDirty)
   {
@@ -187,7 +193,7 @@ SparseCholesky<MatrixType,Cholmod>::matrixL() const
 
 template<typename MatrixType>
 template<typename Derived>
-void SparseCholesky<MatrixType,Cholmod>::solveInPlace(MatrixBase<Derived> &b) const
+void SparseLLT<MatrixType,Cholmod>::solveInPlace(MatrixBase<Derived> &b) const
 {
   const int size = m_matrix.rows();
   ei_assert(size==b.rows());
diff --git a/Eigen/src/Sparse/SparseCholesky.h b/Eigen/src/Sparse/SparseCholesky.h
index 839a84bf3..efe005709 100644
--- a/Eigen/src/Sparse/SparseCholesky.h
+++ b/Eigen/src/Sparse/SparseCholesky.h
@@ -38,25 +38,19 @@ enum {
   MemoryEfficient              = 0x2,
   SupernodalMultifrontal       = 0x4,
   SupernodalLeftLooking        = 0x8
-
-
-/*
-  CholUseEigen    = 0x0,  // Eigen's impl is the default
-  CholUseTaucs    = 0x2,
-  CholUseCholmod  = 0x4*/
 };
 
 /** \ingroup Sparse_Module
   *
-  * \class SparseCholesky
+  * \class SparseLLT
   *
-  * \brief Standard Cholesky decomposition of a matrix and associated features
+  * \brief Standard LLT decomposition of a matrix and associated features
   *
-  * \param MatrixType the type of the matrix of which we are computing the Cholesky decomposition
+  * \param MatrixType the type of the matrix of which we are computing the LLT decomposition
   *
-  * \sa class Cholesky, class CholeskyWithoutSquareRoot
+  * \sa class LLT, class LDLT
   */
-template<typename MatrixType, int Backend = DefaultBackend> class SparseCholesky
+template<typename MatrixType, int Backend = DefaultBackend> class SparseLLT
 {
   protected:
     typedef typename MatrixType::Scalar Scalar;
@@ -70,66 +64,66 @@ template<typename MatrixType, int Backend = DefaultBackend> class SparseCholesky
 
   public:
 
-    SparseCholesky(const MatrixType& matrix, int flags = 0)
+    SparseLLT(int flags = 0)
+      : m_flags(flags), m_status(0)
+    {
+      m_precision = RealScalar(0.1) * Eigen::precision<RealScalar>();
+    }
+
+    SparseLLT(const MatrixType& matrix, int flags = 0)
       : m_matrix(matrix.rows(), matrix.cols()), m_flags(flags), m_status(0)
     {
+      m_precision = RealScalar(0.1) * Eigen::precision<RealScalar>();
       compute(matrix);
     }
 
-    inline const CholMatrixType& matrixL(void) const { return m_matrix; }
+    void setPrecision(RealScalar v) { m_precision = v; }
+    RealScalar precision() const { return m_precision; }
 
-    /** \returns true if the matrix is positive definite */
-    inline bool isPositiveDefinite(void) const { return m_isPositiveDefinite; }
+    void setFlags(int f) { m_flags = f; }
+    int flags() const { return m_flags; }
+
+    void compute(const MatrixType& matrix);
+
+    inline const CholMatrixType& matrixL(void) const { return m_matrix; }
 
     template<typename Derived>
     void solveInPlace(MatrixBase<Derived> &b) const;
 
-    void compute(const MatrixType& matrix);
+    /** \returns true if the factorization succeeded */
+    inline bool succeeded(void) const { return m_succeeded; }
 
   protected:
-    /** \internal
-      * Used to compute and store L
-      * The strict upper part is not used and even not initialized.
-      */
     CholMatrixType m_matrix;
+    RealScalar m_precision;
     int m_flags;
     mutable int m_status;
-    bool m_isPositiveDefinite;
+    bool m_succeeded;
 };
 
-/** Computes / recomputes the Cholesky decomposition A = LL^* = U^*U of \a matrix
+/** Computes / recomputes the LLT decomposition A = LL^* = U^*U of \a matrix
   */
 template<typename MatrixType, int Backend>
-void SparseCholesky<MatrixType,Backend>::compute(const MatrixType& a)
+void SparseLLT<MatrixType,Backend>::compute(const MatrixType& a)
 {
   assert(a.rows()==a.cols());
   const int size = a.rows();
   m_matrix.resize(size, size);
-  const RealScalar eps = ei_sqrt(precision<Scalar>());
+//   const RealScalar eps = ei_sqrt(precision<Scalar>());
 
   // allocate a temporary vector for accumulations
   AmbiVector<Scalar> tempVector(size);
+  RealScalar density = a.nonZeros()/RealScalar(size*size);
 
   // TODO estimate the number of nnz
   m_matrix.startFill(a.nonZeros()*2);
   for (int j = 0; j < size; ++j)
   {
-//     std::cout << j << "\n";
     Scalar x = ei_real(a.coeff(j,j));
     int endSize = size-j-1;
 
-    // TODO estimate the number of non zero entries
-//       float ratioLhs = float(lhs.nonZeros())/float(lhs.rows()*lhs.cols());
-//       float avgNnzPerRhsColumn = float(rhs.nonZeros())/float(cols);
-//       float ratioRes = std::min(ratioLhs * avgNnzPerRhsColumn, 1.f);
-
-        // let's do a more accurate determination of the nnz ratio for the current column j of res
-        //float ratioColRes = std::min(ratioLhs * rhs.innerNonZeros(j), 1.f);
-        // FIXME find a nice way to get the number of nonzeros of a sub matrix (here an inner vector)
-//         float ratioColRes = ratioRes;
-//         if (ratioColRes>0.1)
-//     tempVector.init(IsSparse);
-    tempVector.init(IsDense);
+    // TODO better estimate the density !
+    tempVector.init(density>0.001? IsDense : IsSparse);
     tempVector.setBounds(j+1,size);
     tempVector.setZero();
     // init with current matrix a
@@ -161,7 +155,7 @@ void SparseCholesky<MatrixType,Backend>::compute(const MatrixType& a)
     RealScalar rx = ei_sqrt(ei_real(x));
     m_matrix.fill(j,j) = rx;
     Scalar y = Scalar(1)/rx;
-    for (typename AmbiVector<Scalar>::Iterator it(tempVector); it; ++it)
+    for (typename AmbiVector<Scalar>::Iterator it(tempVector, m_precision*rx); it; ++it)
     {
       m_matrix.fill(it.index(), j) = it.value() * y;
     }
@@ -171,7 +165,7 @@ void SparseCholesky<MatrixType,Backend>::compute(const MatrixType& a)
 
 template<typename MatrixType, int Backend>
 template<typename Derived>
-void SparseCholesky<MatrixType, Backend>::solveInPlace(MatrixBase<Derived> &b) const
+void SparseLLT<MatrixType, Backend>::solveInPlace(MatrixBase<Derived> &b) const
 {
   const int size = m_matrix.rows();
   ei_assert(size==b.rows());
diff --git a/Eigen/src/Sparse/TaucsSupport.h b/Eigen/src/Sparse/TaucsSupport.h
index a0a5b4b71..a92a63b57 100644
--- a/Eigen/src/Sparse/TaucsSupport.h
+++ b/Eigen/src/Sparse/TaucsSupport.h
@@ -79,10 +79,10 @@ SparseMatrix<Scalar,Flags> SparseMatrix<Scalar,Flags>::Map(taucs_ccs_matrix& tau
 }
 
 template<typename MatrixType>
-class SparseCholesky<MatrixType,Taucs> : public SparseCholesky<MatrixType>
+class SparseLLT<MatrixType,Taucs> : public SparseLLT<MatrixType>
 {
   protected:
-    typedef SparseCholesky<MatrixType> Base;
+    typedef SparseLLT<MatrixType> Base;
     using Base::Scalar;
     using Base::RealScalar;
     using Base::MatrixLIsDirty;
@@ -93,13 +93,18 @@ class SparseCholesky<MatrixType,Taucs> : public SparseCholesky<MatrixType>
 
   public:
 
-    SparseCholesky(const MatrixType& matrix, int flags = 0)
+    SparseLLT(int flags = 0)
+      : Base(flags), m_taucsSupernodalFactor(0)
+    {
+    }
+
+    SparseLLT(const MatrixType& matrix, int flags = 0)
       : Base(matrix, flags), m_taucsSupernodalFactor(0)
     {
       compute(matrix);
     }
 
-    ~SparseCholesky()
+    ~SparseLLT()
     {
       if (m_taucsSupernodalFactor)
         taucs_supernodal_factor_free(m_taucsSupernodalFactor);
@@ -117,7 +122,7 @@ class SparseCholesky<MatrixType,Taucs> : public SparseCholesky<MatrixType>
 };
 
 template<typename MatrixType>
-void SparseCholesky<MatrixType,Taucs>::compute(const MatrixType& a)
+void SparseLLT<MatrixType,Taucs>::compute(const MatrixType& a)
 {
   if (m_taucsSupernodalFactor)
   {
@@ -128,7 +133,7 @@ void SparseCholesky<MatrixType,Taucs>::compute(const MatrixType& a)
   if (m_flags & IncompleteFactorization)
   {
     taucs_ccs_matrix taucsMatA = const_cast<MatrixType&>(a).asTaucsMatrix();
-    taucs_ccs_matrix* taucsRes = taucs_ccs_factor_llt(&taucsMatA, 0, 0);
+    taucs_ccs_matrix* taucsRes = taucs_ccs_factor_llt(&taucsMatA, Base::m_precision, 0);
     m_matrix = Base::CholMatrixType::Map(*taucsRes);
     free(taucsRes);
     m_status = (m_status & ~(CompleteFactorization|MatrixLIsDirty))
@@ -153,8 +158,8 @@ void SparseCholesky<MatrixType,Taucs>::compute(const MatrixType& a)
 }
 
 template<typename MatrixType>
-inline const typename SparseCholesky<MatrixType>::CholMatrixType&
-SparseCholesky<MatrixType,Taucs>::matrixL() const
+inline const typename SparseLLT<MatrixType>::CholMatrixType&
+SparseLLT<MatrixType,Taucs>::matrixL() const
 {
   if (m_status & MatrixLIsDirty)
   {
@@ -170,7 +175,7 @@ SparseCholesky<MatrixType,Taucs>::matrixL() const
 
 template<typename MatrixType>
 template<typename Derived>
-void SparseCholesky<MatrixType,Taucs>::solveInPlace(MatrixBase<Derived> &b) const
+void SparseLLT<MatrixType,Taucs>::solveInPlace(MatrixBase<Derived> &b) const
 {
   const int size = m_matrix.rows();
   ei_assert(size==b.rows());
@@ -179,9 +184,9 @@ void SparseCholesky<MatrixType,Taucs>::solveInPlace(MatrixBase<Derived> &b) cons
   {
 //     ei_assert(!(m_status & SupernodalFactorIsDirty));
 //     taucs_supernodal_solve_llt(m_taucsSupernodalFactor,double* b);
-    matrixL();
+    //matrixL();
   }
-//   else
+  else
   {
     Base::solveInPlace(b);
   }
diff --git a/Eigen/src/Sparse/TriangularSolver.h b/Eigen/src/Sparse/TriangularSolver.h
index 3a1a33b3c..8e71c936d 100644
--- a/Eigen/src/Sparse/TriangularSolver.h
+++ b/Eigen/src/Sparse/TriangularSolver.h
@@ -102,7 +102,7 @@ struct ei_solve_triangular_selector<Lhs,Rhs,Lower,ColMajor|IsSparse>
         typename Lhs::InnerIterator it(lhs, i);
         if(!(Lhs::Flags & UnitDiagBit))
         {
-          std::cerr << it.value() << " ; " << it.index() << " == " << i << "\n";
+          // std::cerr << it.value() << " ; " << it.index() << " == " << i << "\n";
           ei_assert(it.index()==i);
           other.coeffRef(i,col) /= it.value();
         }
diff --git a/bench/BenchSparseUtil.h b/bench/BenchSparseUtil.h
index 2c24c29e6..35c9a5263 100644
--- a/bench/BenchSparseUtil.h
+++ b/bench/BenchSparseUtil.h
@@ -16,7 +16,7 @@ USING_PART_OF_NAMESPACE_EIGEN
 #endif
 
 #ifndef SCALAR
-#define SCALAR float
+#define SCALAR double
 #endif
 
 typedef SCALAR Scalar;
diff --git a/bench/benchCholesky.cpp b/bench/benchCholesky.cpp
index f64b61b71..e998d8536 100644
--- a/bench/benchCholesky.cpp
+++ b/bench/benchCholesky.cpp
@@ -1,5 +1,5 @@
 
-// g++ -DNDEBUG -O3 -I.. benchCholesky.cpp  -o benchCholesky && ./benchCholesky
+// g++ -DNDEBUG -O3 -I.. benchLLT.cpp  -o benchLLT && ./benchLLT
 // options:
 //  -DBENCH_GSL -lgsl /usr/lib/libcblas.so.3
 //  -DEIGEN_DONT_VECTORIZE
@@ -9,7 +9,7 @@
 //  -DSCALAR=double
 
 #include <Eigen/Array>
-#include <Eigen/Cholesky>
+#include <Eigen/LLT>
 #include <bench/BenchUtil.h>
 using namespace Eigen;
 
@@ -24,7 +24,7 @@ using namespace Eigen;
 typedef float Scalar;
 
 template <typename MatrixType>
-__attribute__ ((noinline)) void benchCholesky(const MatrixType& m)
+__attribute__ ((noinline)) void benchLLT(const MatrixType& m)
 {
   int rows = m.rows();
   int cols = m.cols();
@@ -54,7 +54,7 @@ __attribute__ ((noinline)) void benchCholesky(const MatrixType& m)
     timerNoSqrt.start();
     for (int k=0; k<repeats; ++k)
     {
-      CholeskyWithoutSquareRoot<SquareMatrixType> cholnosqrt(covMat);
+      LDLT<SquareMatrixType> cholnosqrt(covMat);
       acc += cholnosqrt.matrixL().coeff(r,c);
     }
     timerNoSqrt.stop();
@@ -65,7 +65,7 @@ __attribute__ ((noinline)) void benchCholesky(const MatrixType& m)
     timerSqrt.start();
     for (int k=0; k<repeats; ++k)
     {
-      Cholesky<SquareMatrixType> chol(covMat);
+      LLT<SquareMatrixType> chol(covMat);
       acc += chol.matrixL().coeff(r,c);
     }
     timerSqrt.stop();
@@ -124,17 +124,17 @@ int main(int argc, char* argv[])
   std::cout << "\n";
 
   for (uint i=0; dynsizes[i]>0; ++i)
-    benchCholesky(Matrix<Scalar,Dynamic,Dynamic>(dynsizes[i],dynsizes[i]));
-
-//   benchCholesky(Matrix<Scalar,2,2>());
-//   benchCholesky(Matrix<Scalar,3,3>());
-//   benchCholesky(Matrix<Scalar,4,4>());
-//   benchCholesky(Matrix<Scalar,5,5>());
-//   benchCholesky(Matrix<Scalar,6,6>());
-//   benchCholesky(Matrix<Scalar,7,7>());
-//   benchCholesky(Matrix<Scalar,8,8>());
-//   benchCholesky(Matrix<Scalar,12,12>());
-//   benchCholesky(Matrix<Scalar,16,16>());
+    benchLLT(Matrix<Scalar,Dynamic,Dynamic>(dynsizes[i],dynsizes[i]));
+
+//   benchLLT(Matrix<Scalar,2,2>());
+//   benchLLT(Matrix<Scalar,3,3>());
+//   benchLLT(Matrix<Scalar,4,4>());
+//   benchLLT(Matrix<Scalar,5,5>());
+//   benchLLT(Matrix<Scalar,6,6>());
+//   benchLLT(Matrix<Scalar,7,7>());
+//   benchLLT(Matrix<Scalar,8,8>());
+//   benchLLT(Matrix<Scalar,12,12>());
+//   benchLLT(Matrix<Scalar,16,16>());
   return 0;
 }
 
diff --git a/bench/sparse_product.cpp b/bench/sparse_product.cpp
index 846301fa5..edeb08c5d 100644
--- a/bench/sparse_product.cpp
+++ b/bench/sparse_product.cpp
@@ -21,6 +21,18 @@
 #define MINDENSITY 0.0004
 #endif
 
+#ifndef NBTRIES
+#define NBTRIES 10
+#endif
+
+#define BENCH(X) \
+  timer.reset(); \
+  for (int _j=0; _j<NBTRIES; ++_j) { \
+    timer.start(); \
+    for (int _k=0; _k<REPEAT; ++_k) { \
+        X  \
+  } timer.stop(); }
+
 int main(int argc, char *argv[])
 {
   int rows = SIZE;
@@ -77,32 +89,34 @@ int main(int argc, char *argv[])
     {
       std::cout << "Eigen sparse\t" << density*100 << "%\n";
 
-      timer.reset();
-      timer.start();
-      for (int k=0; k<REPEAT; ++k)
-        sm3 = sm1 * sm2;
-      timer.stop();
+//       timer.reset();
+//       timer.start();
+      BENCH(for (int k=0; k<REPEAT; ++k) sm3 = sm1 * sm2;)
+//       timer.stop();
       std::cout << "   a * b:\t" << timer.value() << endl;
 
       timer.reset();
       timer.start();
-      for (int k=0; k<REPEAT; ++k)
-        sm3 = sm1.transpose() * sm2;
-      timer.stop();
+//       std::cerr << "transpose...\n";
+//       EigenSparseMatrix sm4 = sm1.transpose();
+//       std::cout << sm4.nonZeros() << " == " << sm1.nonZeros() << "\n";
+//       exit(1);
+//       std::cerr << "transpose OK\n";
+//       std::cout << sm1 << "\n\n" << sm1.transpose() << "\n\n" << sm4.transpose() << "\n\n";
+      BENCH(for (int k=0; k<REPEAT; ++k) sm3 = sm1.transpose() * sm2;)
+//       timer.stop();
       std::cout << "   a' * b:\t" << timer.value() << endl;
 
-      timer.reset();
-      timer.start();
-      for (int k=0; k<REPEAT; ++k)
-        sm3 = sm1.transpose() * sm2.transpose();
-      timer.stop();
+//       timer.reset();
+//       timer.start();
+      BENCH( for (int k=0; k<REPEAT; ++k) sm3 = sm1.transpose() * sm2.transpose(); )
+//       timer.stop();
       std::cout << "   a' * b':\t" << timer.value() << endl;
 
-      timer.reset();
-      timer.start();
-      for (int k=0; k<REPEAT; ++k)
-        sm3 = sm1 * sm2.transpose();
-      timer.stop();
+//       timer.reset();
+//       timer.start();
+      BENCH( for (int k=0; k<REPEAT; ++k) sm3 = sm1 * sm2.transpose(); )
+//       timer.stop();
       std::cout << "   a * b' :\t" << timer.value() << endl;
     }
 
diff --git a/doc/snippets/Cholesky_solve.cpp b/doc/snippets/LLT_solve.cpp
index ac743cb55..76ab09ec5 100644
--- a/doc/snippets/Cholesky_solve.cpp
+++ b/doc/snippets/LLT_solve.cpp
@@ -2,5 +2,7 @@ typedef Matrix<float,Dynamic,2> DataMatrix;
 // let's generate some samples on the 3D plane of equation z = 2x+3y (with some noise)
 DataMatrix samples = DataMatrix::Random(12,2);
 VectorXf elevations = 2*samples.col(0) + 3*samples.col(1) + VectorXf::Random(12)*0.1;
-// and let's solve samples * x = elevations in least square sense:
-cout << (samples.adjoint() * samples).cholesky().solve((samples.adjoint()*elevations).eval()) << endl;
+// and let's solve samples * [x y]^T = elevations in least square sense:
+Matrix<float,2,1> xy;
+(samples.adjoint() * samples).llt().solve((samples.adjoint()*elevations), &xy);
+cout << xy << endl;
diff --git a/test/cholesky.cpp b/test/cholesky.cpp
index 80614346c..f7eb7800e 100644
--- a/test/cholesky.cpp
+++ b/test/cholesky.cpp
@@ -33,7 +33,7 @@
 template<typename MatrixType> void cholesky(const MatrixType& m)
 {
   /* this test covers the following files:
-     Cholesky.h CholeskyWithoutSquareRoot.h
+     LLT.h LDLT.h
   */
   int rows = m.rows();
   int cols = m.cols();
@@ -44,8 +44,8 @@ template<typename MatrixType> void cholesky(const MatrixType& m)
   typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> VectorType;
 
   MatrixType a0 = MatrixType::Random(rows,cols);
-  VectorType vecB = VectorType::Random(rows);
-  MatrixType matB = MatrixType::Random(rows,cols);
+  VectorType vecB = VectorType::Random(rows), vecX(rows);
+  MatrixType matB = MatrixType::Random(rows,cols), matX(rows,cols);
   SquareMatrixType symm =  a0 * a0.adjoint();
   // let's make sure the matrix is not singular or near singular
   MatrixType a1 = MatrixType::Random(rows,cols);
@@ -80,28 +80,32 @@ template<typename MatrixType> void cholesky(const MatrixType& m)
   #endif
 
   {
-    CholeskyWithoutSquareRoot<SquareMatrixType> cholnosqrt(symm);
-    VERIFY(cholnosqrt.isPositiveDefinite());
-    VERIFY_IS_APPROX(symm, cholnosqrt.matrixL() * cholnosqrt.vectorD().asDiagonal() * cholnosqrt.matrixL().adjoint());
-    VERIFY_IS_APPROX(symm * cholnosqrt.solve(vecB), vecB);
-    VERIFY_IS_APPROX(symm * cholnosqrt.solve(matB), matB);
+    LDLT<SquareMatrixType> ldlt(symm);
+    VERIFY(ldlt.isPositiveDefinite());
+    VERIFY_IS_APPROX(symm, ldlt.matrixL() * ldlt.vectorD().asDiagonal() * ldlt.matrixL().adjoint());
+    ldlt.solve(vecB, &vecX);
+    VERIFY_IS_APPROX(symm * vecX, vecB);
+    ldlt.solve(matB, &matX);
+    VERIFY_IS_APPROX(symm * matX, matB);
   }
 
   {
-    Cholesky<SquareMatrixType> chol(symm);
+    LLT<SquareMatrixType> chol(symm);
     VERIFY(chol.isPositiveDefinite());
     VERIFY_IS_APPROX(symm, chol.matrixL() * chol.matrixL().adjoint());
-    VERIFY_IS_APPROX(symm * chol.solve(vecB), vecB);
-    VERIFY_IS_APPROX(symm * chol.solve(matB), matB);
+    chol.solve(vecB, &vecX);
+    VERIFY_IS_APPROX(symm * vecX, vecB);
+    chol.solve(matB, &matX);
+    VERIFY_IS_APPROX(symm * matX, matB);
   }
 
   // test isPositiveDefinite on non definite matrix
   if (rows>4)
   {
     SquareMatrixType symm =  a0.block(0,0,rows,cols-4) * a0.block(0,0,rows,cols-4).adjoint();
-    Cholesky<SquareMatrixType> chol(symm);
+    LLT<SquareMatrixType> chol(symm);
     VERIFY(!chol.isPositiveDefinite());
-    CholeskyWithoutSquareRoot<SquareMatrixType> cholnosqrt(symm);
+    LDLT<SquareMatrixType> cholnosqrt(symm);
     VERIFY(!cholnosqrt.isPositiveDefinite());
   }
 }
diff --git a/test/sparse.cpp b/test/sparse.cpp
index 39ea05b8b..f54835972 100644
--- a/test/sparse.cpp
+++ b/test/sparse.cpp
@@ -217,7 +217,7 @@ template<typename Scalar> void sparse(int rows, int cols)
     // TODO test row major
   }
 
-  // test Cholesky
+  // test LLT
   {
   }