diff --git a/Eigen/src/Sparse/CholmodSupport.h b/Eigen/src/Sparse/CholmodSupport.h
index 6689a4669..f77edba42 100644
--- a/Eigen/src/Sparse/CholmodSupport.h
+++ b/Eigen/src/Sparse/CholmodSupport.h
@@ -25,6 +25,35 @@
 #ifndef EIGEN_CHOLMODSUPPORT_H
 #define EIGEN_CHOLMODSUPPORT_H
 
+template<typename Scalar, typename CholmodType>
+void ei_cholmod_configure_matrix(CholmodType& mat)
+{
+  if (ei_is_same_type<Scalar,float>::ret)
+  {
+    mat.xtype = CHOLMOD_REAL;
+    mat.dtype = 1;
+  }
+  else if (ei_is_same_type<Scalar,double>::ret)
+  {
+    mat.xtype = CHOLMOD_REAL;
+    mat.dtype = 0;
+  }
+  else if (ei_is_same_type<Scalar,std::complex<float> >::ret)
+  {
+    mat.xtype = CHOLMOD_COMPLEX;
+    mat.dtype = 1;
+  }
+  else if (ei_is_same_type<Scalar,std::complex<double> >::ret)
+  {
+    mat.xtype = CHOLMOD_COMPLEX;
+    mat.dtype = 0;
+  }
+  else
+  {
+    ei_assert(false && "Scalar type not supported by CHOLMOD");
+  }
+}
+
 template<typename Scalar, int Flags>
 cholmod_sparse SparseMatrix<Scalar,Flags>::asCholmodMatrix()
 {
@@ -42,30 +71,7 @@ cholmod_sparse SparseMatrix<Scalar,Flags>::asCholmodMatrix()
   res.dtype = 0;
   res.stype = -1;
 
-  if (ei_is_same_type<Scalar,float>::ret)
-  {
-    res.xtype = CHOLMOD_REAL;
-    res.dtype = 1;
-  }
-  else if (ei_is_same_type<Scalar,double>::ret)
-  {
-    res.xtype = CHOLMOD_REAL;
-    res.dtype = 0;
-  }
-  else if (ei_is_same_type<Scalar,std::complex<float> >::ret)
-  {
-    res.xtype = CHOLMOD_COMPLEX;
-    res.dtype = 1;
-  }
-  else if (ei_is_same_type<Scalar,std::complex<double> >::ret)
-  {
-    res.xtype = CHOLMOD_COMPLEX;
-    res.dtype = 0;
-  }
-  else
-  {
-    ei_assert(false && "Scalar type not supported by CHOLMOD");
-  }
+  ei_cholmod_configure_matrix<Scalar>(res);
 
   if (Flags & SelfAdjoint)
   {
@@ -82,6 +88,24 @@ cholmod_sparse SparseMatrix<Scalar,Flags>::asCholmodMatrix()
   return res;
 }
 
+template<typename Derived>
+cholmod_dense ei_cholmod_map_eigen_to_dense(MatrixBase<Derived>& mat)
+{
+  typedef typename Derived::Scalar Scalar;
+
+  cholmod_dense res;
+  res.nrow   = mat.rows();
+  res.ncol   = mat.cols();
+  res.nzmax  = res.nrow * res.ncol;
+  res.d      = mat.derived().stride();
+  res.x      = mat.derived().data();
+  res.z      = 0;
+
+  ei_cholmod_configure_matrix<Scalar>(res);
+
+  return res;
+}
+
 template<typename Scalar, int Flags>
 SparseMatrix<Scalar,Flags> SparseMatrix<Scalar,Flags>::Map(cholmod_sparse& cm)
 {
@@ -103,7 +127,7 @@ class SparseLLT<MatrixType,Cholmod> : public SparseLLT<MatrixType>
 {
   protected:
     typedef SparseLLT<MatrixType> Base;
-    using Base::Scalar;
+    using typename Base::Scalar;
     using Base::RealScalar;
     using Base::MatrixLIsDirty;
     using Base::SupernodalFactorIsDirty;
@@ -155,11 +179,12 @@ void SparseLLT<MatrixType,Cholmod>::compute(const MatrixType& a)
   }
 
   cholmod_sparse A = const_cast<MatrixType&>(a).asCholmodMatrix();
+  m_cholmod.supernodal = CHOLMOD_AUTO;
   // TODO
   if (m_flags&IncompleteFactorization)
   {
     m_cholmod.nmethods = 1;
-    m_cholmod.method [0].ordering = CHOLMOD_NATURAL;
+    m_cholmod.method[0].ordering = CHOLMOD_NATURAL;
     m_cholmod.postorder = 0;
   }
   else
@@ -196,13 +221,19 @@ template<typename MatrixType>
 template<typename Derived>
 void SparseLLT<MatrixType,Cholmod>::solveInPlace(MatrixBase<Derived> &b) const
 {
-  if (m_status & MatrixLIsDirty)
-    matrixL();
-
-  const int size = m_matrix.rows();
+  const int size = m_cholmodFactor->n;
   ei_assert(size==b.rows());
 
+  // this uses Eigen's triangular sparse solver
+  if (m_status & MatrixLIsDirty)
+    matrixL();
   Base::solveInPlace(b);
+  // as long as our own triangular sparse solver is not fully optimal,
+  // let's use CHOLMOD's one:
+//   cholmod_dense cdb = ei_cholmod_map_eigen_to_dense(b);
+//   cholmod_dense* x = cholmod_solve(CHOLMOD_LDLt, m_cholmodFactor, &cdb, &m_cholmod);
+//   b = Matrix<typename Base::Scalar,Dynamic,1>::Map(reinterpret_cast<typename Base::Scalar*>(x->x),b.rows());
+//   cholmod_free_dense(&x, &m_cholmod);
 }
 
 #endif // EIGEN_CHOLMODSUPPORT_H
diff --git a/Eigen/src/Sparse/SparseLLT.h b/Eigen/src/Sparse/SparseLLT.h
index 0b63f80ab..e7c314c2c 100644
--- a/Eigen/src/Sparse/SparseLLT.h
+++ b/Eigen/src/Sparse/SparseLLT.h
@@ -190,8 +190,14 @@ bool SparseLLT<MatrixType, Backend>::solveInPlace(MatrixBase<Derived> &b) const
   ei_assert(size==b.rows());
 
   m_matrix.solveTriangularInPlace(b);
-  // FIXME should be .adjoint() but it fails to compile...
-  m_matrix.transpose().solveTriangularInPlace(b);
+  // FIXME should be simply .adjoint() but it fails to compile...
+  if (NumTraits<Scalar>::IsComplex)
+  {
+    CholMatrixType aux = m_matrix.conjugate();
+    aux.transpose().solveTriangularInPlace(b);
+  }
+  else
+    m_matrix.transpose().solveTriangularInPlace(b);
 
   return true;
 }
diff --git a/Eigen/src/Sparse/SparseMatrix.h b/Eigen/src/Sparse/SparseMatrix.h
index 5a6b2e0d4..f41a9586d 100644
--- a/Eigen/src/Sparse/SparseMatrix.h
+++ b/Eigen/src/Sparse/SparseMatrix.h
@@ -148,7 +148,7 @@ class SparseMatrix
       */
     inline void startFill(int reserveSize = 1000)
     {
-      std::cerr << this << " startFill\n";
+//       std::cerr << this << " startFill\n";
       setZero();
       m_data.reserve(reserveSize);
     }
@@ -214,7 +214,7 @@ class SparseMatrix
       m_data.index(id+1) = inner;
       //return (m_data.value(id+1) = 0);
       m_data.value(id+1) = 0;
-      std::cerr << m_outerIndex[outer] << " " << m_outerIndex[outer+1] << "\n";
+//       std::cerr << m_outerIndex[outer] << " " << m_outerIndex[outer+1] << "\n";
       return m_data.value(id+1);
     }
 
@@ -222,7 +222,7 @@ class SparseMatrix
 
     inline void endFill()
     {
-      std::cerr << this << " endFill\n";
+//       std::cerr << this << " endFill\n";
       int size = m_data.size();
       int i = m_outerSize;
       // find the last filled column
@@ -238,7 +238,7 @@ class SparseMatrix
 
     void resize(int rows, int cols)
     {
-      std::cerr << this << " resize " << rows << "x" << cols << "\n";
+//       std::cerr << this << " resize " << rows << "x" << cols << "\n";
       const int outerSize = RowMajor ? rows : cols;
       m_innerSize = RowMajor ? cols : rows;
       m_data.clear();
@@ -273,6 +273,12 @@ class SparseMatrix
       *this = other.derived();
     }
 
+    inline SparseMatrix(const SparseMatrix& other)
+      : m_outerSize(0), m_innerSize(0), m_outerIndex(0)
+    {
+      *this = other.derived();
+    }
+
     inline void swap(SparseMatrix& other)
     {
       //EIGEN_DBG_SPARSE(std::cout << "SparseMatrix:: swap\n");
diff --git a/test/sparse_solvers.cpp b/test/sparse_solvers.cpp
index da9a245a2..8ce4e2264 100644
--- a/test/sparse_solvers.cpp
+++ b/test/sparse_solvers.cpp
@@ -24,6 +24,27 @@
 
 #include "sparse.h"
 
+template<typename Scalar> void
+initSPD(double density,
+        Matrix<Scalar,Dynamic,Dynamic>& refMat,
+        SparseMatrix<Scalar>& sparseMat)
+{
+  Matrix<Scalar,Dynamic,Dynamic> aux(refMat.rows(),refMat.cols());
+  initSparse(density,refMat,sparseMat);
+  refMat = refMat * refMat.adjoint();
+  for (int k=0; k<2; ++k)
+  {
+    initSparse(density,aux,sparseMat,ForceNonZeroDiag);
+    refMat += aux * aux.adjoint();
+  }
+  sparseMat.startFill();
+  for (int j=0 ; j<sparseMat.cols(); ++j)
+    for (int i=j ; i<sparseMat.rows(); ++i)
+      if (refMat(i,j)!=Scalar(0))
+        sparseMat.fill(i,j) = refMat(i,j);
+  sparseMat.endFill();
+}
+
 template<typename Scalar> void sparse_solvers(int rows, int cols)
 {
   double density = std::max(8./(rows*cols), 0.01);
@@ -56,7 +77,6 @@ template<typename Scalar> void sparse_solvers(int rows, int cols)
   }
 
   // test LLT
-  if (!NumTraits<Scalar>::IsComplex)
   {
     // TODO fix the issue with complex (see SparseLLT::solveInPlace)
     SparseMatrix<Scalar> m2(rows, cols);
@@ -65,49 +85,54 @@ template<typename Scalar> void sparse_solvers(int rows, int cols)
     DenseVector b = DenseVector::Random(cols);
     DenseVector refX(cols), x(cols);
 
-    initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag|MakeLowerTriangular, &zeroCoords, &nonzeroCoords);
-    refMat2 += refMat2.adjoint();
-    refMat2.diagonal() *= 0.5;
+    initSPD(density, refMat2, m2);
 
     refMat2.llt().solve(b, &refX);
     typedef SparseMatrix<Scalar,LowerTriangular|SelfAdjoint> SparseSelfAdjointMatrix;
-    x = b;
-    SparseLLT<SparseSelfAdjointMatrix> (m2).solveInPlace(x);
-    //VERIFY(refX.isApprox(x,test_precision<Scalar>()) && "LLT: default");
+    if (!NumTraits<Scalar>::IsComplex)
+    {
+      x = b;
+      SparseLLT<SparseSelfAdjointMatrix> (m2).solveInPlace(x);
+      VERIFY(refX.isApprox(x,test_precision<Scalar>()) && "LLT: default");
+    }
     #ifdef EIGEN_CHOLMOD_SUPPORT
     x = b;
     SparseLLT<SparseSelfAdjointMatrix,Cholmod>(m2).solveInPlace(x);
     VERIFY(refX.isApprox(x,test_precision<Scalar>()) && "LLT: cholmod");
     #endif
-    #ifdef EIGEN_TAUCS_SUPPORT
-    x = b;
-    SparseLLT<SparseSelfAdjointMatrix,Taucs>(m2,IncompleteFactorization).solveInPlace(x);
-    VERIFY(refX.isApprox(x,test_precision<Scalar>()) && "LLT: taucs (IncompleteFactorization)");
-    x = b;
-    SparseLLT<SparseSelfAdjointMatrix,Taucs>(m2,SupernodalMultifrontal).solveInPlace(x);
-    VERIFY(refX.isApprox(x,test_precision<Scalar>()) && "LLT: taucs (SupernodalMultifrontal)");
-    x = b;
-    SparseLLT<SparseSelfAdjointMatrix,Taucs>(m2,SupernodalLeftLooking).solveInPlace(x);
-    VERIFY(refX.isApprox(x,test_precision<Scalar>()) && "LLT: taucs (SupernodalLeftLooking)");
-    #endif
+    if (!NumTraits<Scalar>::IsComplex)
+    {
+      #ifdef EIGEN_TAUCS_SUPPORT
+      x = b;
+      SparseLLT<SparseSelfAdjointMatrix,Taucs>(m2,IncompleteFactorization).solveInPlace(x);
+      VERIFY(refX.isApprox(x,test_precision<Scalar>()) && "LLT: taucs (IncompleteFactorization)");
+      x = b;
+      SparseLLT<SparseSelfAdjointMatrix,Taucs>(m2,SupernodalMultifrontal).solveInPlace(x);
+      VERIFY(refX.isApprox(x,test_precision<Scalar>()) && "LLT: taucs (SupernodalMultifrontal)");
+      x = b;
+      SparseLLT<SparseSelfAdjointMatrix,Taucs>(m2,SupernodalLeftLooking).solveInPlace(x);
+      VERIFY(refX.isApprox(x,test_precision<Scalar>()) && "LLT: taucs (SupernodalLeftLooking)");
+      #endif
+    }
   }
 
   // test LDLT
   if (!NumTraits<Scalar>::IsComplex)
   {
-    // TODO fix the issue with complex (see SparseLDT::solveInPlace)
+    // TODO fix the issue with complex (see SparseLDLT::solveInPlace)
     SparseMatrix<Scalar> m2(rows, cols);
     DenseMatrix refMat2(rows, cols);
 
     DenseVector b = DenseVector::Random(cols);
     DenseVector refX(cols), x(cols);
 
-    initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag|MakeUpperTriangular, &zeroCoords, &nonzeroCoords);
+    //initSPD(density, refMat2, m2);
+    initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag|MakeUpperTriangular, 0, 0);
     refMat2 += refMat2.adjoint();
     refMat2.diagonal() *= 0.5;
 
     refMat2.ldlt().solve(b, &refX);
-    typedef SparseMatrix<Scalar,LowerTriangular|SelfAdjoint> SparseSelfAdjointMatrix;
+    typedef SparseMatrix<Scalar,UpperTriangular|SelfAdjoint> SparseSelfAdjointMatrix;
     x = b;
     SparseLDLT<SparseSelfAdjointMatrix> ldlt(m2);
     if (ldlt.succeeded())
@@ -177,6 +202,6 @@ void test_sparse_solvers()
   for(int i = 0; i < g_repeat; i++) {
     CALL_SUBTEST( sparse_solvers<double>(8, 8) );
     CALL_SUBTEST( sparse_solvers<std::complex<double> >(16, 16) );
-    CALL_SUBTEST( sparse_solvers<double>(33, 33) );
+    CALL_SUBTEST( sparse_solvers<double>(101, 101) );
   }
 }