Updates to the Sparse unsupported solvers module.

* change Sparse* specialization's signatures from <..., int Backend> to <..., typename Backend>. Update SparseExtra accordingly to use structs instead of the SparseBackend enum.
* add SparseLDLT Cholmod specialization
* for Cholmod and UmfPack, SparseLU, SparseLLT and SparseLDLT now use ei_solve_retval and have the new solve() method (to be closer to the 3.0 API).

* fix doc

doc/examples/Tutorial_ReductionsVisitorsBroadcasting_broadcast_1nn.cpp~

@@ -1,24 +0,0 @@
-#include <iostream>
-#include <Eigen/Dense>
-
-using namespace std;
-using namespace Eigen;
-
-int main()
-{
-  Eigen::MatrixXf m(2,4);
-  Eigen::VectorXf v(2);
-  
-  m << 1, 2, 6, 9,
-       3, 1, 7, 2;
-       
-  v << 2,
-       3;
-
-  MatrixXf::Index index;
-  // find nearest neighbour
-  (m.colwise() - v).colwise().squaredNorm().minCoeff(&index);
-
-  cout << "Nearest neighbour is column " << index << ":" << endl;
-  cout << m.col(index) << endl;
-}

doc/examples/Tutorial_ReductionsVisitorsBroadcasting_broadcast_simple.cpp~

@@ -1,21 +0,0 @@
-#include <iostream>
-#include <Eigen/Dense>
-
-using namespace std;
-int main()
-{
-  Eigen::MatrixXf mat(2,4);
-  Eigen::MatrixXf v(2);
-  
-  mat << 1, 2, 6, 9,
-         3, 1, 7, 2;
-         
-  v << 0,
-       1;
-       
-  //add v to each column of m
-  mat.colwise() += v;
-  
-  std::cout << "Broadcasting result: " << std::endl;
-  std::cout << mat << std::endl;
-}

doc/examples/Tutorial_ReductionsVisitorsBroadcasting_reductions_norm.cpp~

@@ -1,28 +0,0 @@
-#include <Eigen/Dense>
-#include <iostream>
-
-using namespace std;
-using namespace Eigen;
-
-int main()
-{
-  VectorXf v(2);
-  MatrixXf m(2,2), n(2,2);
-  
-  v << 2,
-       5;
-  
-  m << 0,2,
-       3,4;
-
-  n << 1,2,
-       3,4;
-  
-  cout << "v.norm() = " << m.norm() << endl;
-  cout << "m.norm() = " << m.norm() << endl;
-  cout << "n.norm() = " << m.norm() << endl;
-  cout << endl;
-  cout << "v.squaredNorm() = " << v.squaredNorm() << endl;
-  cout << "m.squaredNorm() = " << m.squaredNorm() << endl;
-  cout << "n.squaredNorm() = " << n.squaredNorm() << endl;
-}

doc/examples/Tutorial_ReductionsVisitorsBroadcasting_visitors.cpp~

@@ -1,26 +0,0 @@
-#include <iostream>
-#include <Eigen/Dense>
-
-using namespace std;
-using namespace Eigen;
-
-int main()
-{
-  Eigen::MatrixXf m(2,2);
-  
-  m << 1, 2,
-       3, 4;
-
-  //get location of maximum
-  MatrixXf::Index maxRow, maxCol;
-  m.maxCoeff(&maxRow, &maxCol);
-
-  //get location of minimum
-  MatrixXf::Index minRow, minCol;
-  m.minCoeff(&minRow, &minCol);
-
-  cout << "Max at: " <<
-     maxRow << "," << maxCol << endl;
-  cout << "Min at: " <<
-     minRow << "," << minCol << endl;
-}

unsupported/Eigen/CholmodSupport

@@ -3,7 +3,7 @@
 
 #include "SparseExtra"
 
-#include "src/Core/util/DisableMSVCWarnings.h"
+#include "../../Eigen/src/Core/util/DisableMSVCWarnings.h"
 
 extern "C" {
   #include <cholmod.h>
@@ -20,11 +20,13 @@ namespace Eigen {
   * \endcode
   */
 
+struct Cholmod {};
+
 #include "src/SparseExtra/CholmodSupport.h"
 
 } // namespace Eigen
 
-#include "src/Core/util/EnableMSVCWarnings.h"
+#include "../../Eigen/src/Core/util/EnableMSVCWarnings.h"
 
 #endif // EIGEN_CHOLMODSUPPORT_MODULE_H
 

unsupported/Eigen/SparseExtra

@@ -27,6 +27,9 @@ namespace Eigen {
   * \endcode
   */
 
+struct DefaultBackend {};
+
+/*
 enum SparseBackend {
   DefaultBackend,
   Taucs,
@@ -34,6 +37,8 @@ enum SparseBackend {
   SuperLU,
   UmfPack
 };
+*/
+
 
 // solver flags
 enum {

unsupported/Eigen/SuperLUSupport

@@ -3,7 +3,7 @@
 
 #include "SparseExtra"
 
-#include "src/Core/util/DisableMSVCWarnings.h"
+#include "../../Eigen/src/Core/util/DisableMSVCWarnings.h"
 
 typedef int int_t;
 #include <slu_Cnames.h>
@@ -36,10 +36,12 @@ namespace Eigen {
   * \endcode
   */
 
+struct SuperLU {};
+
 #include "src/SparseExtra/SuperLUSupport.h"
 
 } // namespace Eigen
 
-#include "src/Core/util/EnableMSVCWarnings.h"
+#include "../../Eigen/src/Core/util/EnableMSVCWarnings.h"
 
 #endif // EIGEN_SUPERLUSUPPORT_MODULE_H

unsupported/Eigen/UmfPackSupport

@@ -3,7 +3,7 @@
 
 #include "SparseExtra"
 
-#include "src/Core/util/DisableMSVCWarnings.h"
+#include "../../Eigen/src/Core/util/DisableMSVCWarnings.h"
 
 #include <umfpack.h>
 
@@ -20,10 +20,12 @@ namespace Eigen {
   * \endcode
   */
 
+struct UmfPack {};
+
 #include "src/SparseExtra/UmfPackSupport.h"
 
 } // namespace Eigen
 
-#include "src/Core/util/EnableMSVCWarnings.h"
+#include "../../Eigen/src/Core/util/EnableMSVCWarnings.h"
 
 #endif // EIGEN_UMFPACKSUPPORT_MODULE_H

unsupported/Eigen/src/SparseExtra/CholmodSupport.h

@@ -25,6 +25,7 @@
 #ifndef EIGEN_CHOLMODSUPPORT_H
 #define EIGEN_CHOLMODSUPPORT_H
 
+
 template<typename Scalar, typename CholmodType>
 void ei_cholmod_configure_matrix(CholmodType& mat)
 {
@@ -54,10 +55,10 @@ void ei_cholmod_configure_matrix(CholmodType& mat)
   }
 }
 
-template<typename MatrixType>
+template<typename _MatrixType>
-cholmod_sparse ei_asCholmodMatrix(MatrixType& mat)
+cholmod_sparse ei_cholmod_map_eigen_to_sparse(_MatrixType& mat)
 {
-  typedef typename MatrixType::Scalar Scalar;
+  typedef typename _MatrixType::Scalar Scalar;
   cholmod_sparse res;
   res.nzmax   = mat.nonZeros();
   res.nrow    = mat.rows();;
@@ -75,11 +76,11 @@ cholmod_sparse ei_asCholmodMatrix(MatrixType& mat)
   ei_cholmod_configure_matrix<Scalar>(res);
 
 
-  if (MatrixType::Flags & SelfAdjoint)
+  if (_MatrixType::Flags & SelfAdjoint)
   {
-    if (MatrixType::Flags & Upper)
+    if (_MatrixType::Flags & Upper)
       res.stype = 1;
-    else if (MatrixType::Flags & Lower)
+    else if (_MatrixType::Flags & Lower)
       res.stype = -1;
     else
       res.stype = 0;
@@ -110,22 +111,23 @@ cholmod_dense ei_cholmod_map_eigen_to_dense(MatrixBase<Derived>& mat)
 }
 
 template<typename Scalar, int Flags, typename Index>
-MappedSparseMatrix<Scalar,Flags,Index> ei_map_cholmod(cholmod_sparse& cm)
+MappedSparseMatrix<Scalar,Flags,Index> ei_map_cholmod_sparse_to_eigen(cholmod_sparse& cm)
 {
   return MappedSparseMatrix<Scalar,Flags,Index>
          (cm.nrow, cm.ncol, reinterpret_cast<Index*>(cm.p)[cm.ncol],
           reinterpret_cast<Index*>(cm.p), reinterpret_cast<Index*>(cm.i),reinterpret_cast<Scalar*>(cm.x) );
 }
 
-template<typename MatrixType>
+
-class SparseLLT<MatrixType,Cholmod> : public SparseLLT<MatrixType>
+
+template<typename _MatrixType>
+class SparseLLT<_MatrixType, Cholmod> : public SparseLLT<_MatrixType>
 {
   protected:
-    typedef SparseLLT<MatrixType> Base;
+    typedef SparseLLT<_MatrixType> Base;
     typedef typename Base::Scalar Scalar;
     typedef typename Base::RealScalar RealScalar;
     typedef typename Base::CholMatrixType CholMatrixType;
-    typedef typename MatrixType::Index Index;
     using Base::MatrixLIsDirty;
     using Base::SupernodalFactorIsDirty;
     using Base::m_flags;
@@ -133,6 +135,8 @@ class SparseLLT<MatrixType,Cholmod> : public SparseLLT<MatrixType>
     using Base::m_status;
 
   public:
+    typedef _MatrixType MatrixType;
+    typedef typename MatrixType::Index Index;
 
     SparseLLT(int flags = 0)
       : Base(flags), m_cholmodFactor(0)
@@ -159,15 +163,71 @@ class SparseLLT<MatrixType,Cholmod> : public SparseLLT<MatrixType>
     template<typename Derived>
     bool solveInPlace(MatrixBase<Derived> &b) const;
 
+    template<typename Rhs>
+    inline const ei_solve_retval<SparseLLT<MatrixType, Cholmod>, Rhs>
+    solve(const MatrixBase<Rhs>& b) const
+    {
+      ei_assert(true && "SparseLLT is not initialized.");
+      return ei_solve_retval<SparseLLT<MatrixType, Cholmod>, Rhs>(*this, b.derived());
+    }
+
     void compute(const MatrixType& matrix);
 
+    inline Index cols() const { return m_matrix.cols(); }
+    inline Index rows() const { return m_matrix.rows(); }
+
+    inline const cholmod_factor* cholmodFactor() const
+    { return m_cholmodFactor; }
+
+    inline cholmod_common* cholmodCommon() const
+    { return &m_cholmod; }
+
+    bool succeeded() const;
+
   protected:
     mutable cholmod_common m_cholmod;
     cholmod_factor* m_cholmodFactor;
 };
 
-template<typename MatrixType>
+
-void SparseLLT<MatrixType,Cholmod>::compute(const MatrixType& a)
+
+template<typename _MatrixType, typename Rhs>
+  struct ei_solve_retval<SparseLLT<_MatrixType, Cholmod>, Rhs>
+  : ei_solve_retval_base<SparseLLT<_MatrixType, Cholmod>, Rhs>
+{
+  typedef SparseLLT<_MatrixType, Cholmod> SpLLTDecType;
+  EIGEN_MAKE_SOLVE_HELPERS(SpLLTDecType,Rhs)
+
+  template<typename Dest> void evalTo(Dest& dst) const
+  {
+    //Index size = dec().cholmodFactor()->n;
+    ei_assert((Index)dec().cholmodFactor()->n==rhs().rows());
+    
+    cholmod_factor* cholmodFactor = const_cast<cholmod_factor*>(dec().cholmodFactor());
+    cholmod_common* cholmodCommon = const_cast<cholmod_common*>(dec().cholmodCommon());
+    // 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(rhs().const_cast_derived());
+    cholmod_dense* x = cholmod_solve(CHOLMOD_A, cholmodFactor, &cdb, cholmodCommon);
+
+    dst = Matrix<typename Base::Scalar,Dynamic,1>::Map(reinterpret_cast<typename Base::Scalar*>(x->x), rhs().rows());  
+
+    cholmod_free_dense(&x, cholmodCommon);
+
+  }
+    
+};
+
+
+
+
+
+template<typename _MatrixType>
+void SparseLLT<_MatrixType,Cholmod>::compute(const _MatrixType& a)
 {
   if (m_cholmodFactor)
   {
@@ -175,7 +235,7 @@ void SparseLLT<MatrixType,Cholmod>::compute(const MatrixType& a)
     m_cholmodFactor = 0;
   }
 
-  cholmod_sparse A = ei_asCholmodMatrix(const_cast<MatrixType&>(a));
+  cholmod_sparse A = ei_cholmod_map_eigen_to_sparse(const_cast<_MatrixType&>(a));
 //   m_cholmod.supernodal = CHOLMOD_AUTO;
   // TODO
 //   if (m_flags&IncompleteFactorization)
@@ -197,16 +257,25 @@ void SparseLLT<MatrixType,Cholmod>::compute(const MatrixType& a)
   m_status = (m_status & ~SupernodalFactorIsDirty) | MatrixLIsDirty;
 }
 
-template<typename MatrixType>
+
-inline const typename SparseLLT<MatrixType,Cholmod>::CholMatrixType&
+// TODO
-SparseLLT<MatrixType,Cholmod>::matrixL() const
+template<typename _MatrixType>
+bool SparseLLT<_MatrixType,Cholmod>::succeeded() const
+{ return true; }
+
+
+
+template<typename _MatrixType>
+inline const typename SparseLLT<_MatrixType,Cholmod>::CholMatrixType&
+SparseLLT<_MatrixType,Cholmod>::matrixL() const
 {
   if (m_status & MatrixLIsDirty)
   {
     ei_assert(!(m_status & SupernodalFactorIsDirty));
 
     cholmod_sparse* cmRes = cholmod_factor_to_sparse(m_cholmodFactor, &m_cholmod);
-    const_cast<typename Base::CholMatrixType&>(m_matrix) = ei_map_cholmod<Scalar,ColMajor,Index>(*cmRes);
+    const_cast<typename Base::CholMatrixType&>(m_matrix) = 
+      ei_map_cholmod_sparse_to_eigen<Scalar,ColMajor,Index>(*cmRes);
     free(cmRes);
 
     m_status = (m_status & ~MatrixLIsDirty);
@@ -214,30 +283,234 @@ SparseLLT<MatrixType,Cholmod>::matrixL() const
   return m_matrix;
 }
 
-template<typename MatrixType>
+
+
+
+template<typename _MatrixType>
 template<typename Derived>
-bool SparseLLT<MatrixType,Cholmod>::solveInPlace(MatrixBase<Derived> &b) const
+bool SparseLLT<_MatrixType,Cholmod>::solveInPlace(MatrixBase<Derived> &b) const
 {
-  const Index size = m_cholmodFactor->n;
+  //Index size = m_cholmodFactor->n;
-  ei_assert(size==b.rows());
+  ei_assert((Index)m_cholmodFactor->n==b.rows());
 
   // this uses Eigen's triangular sparse solver
-//   if (m_status & MatrixLIsDirty)
+  //   if (m_status & MatrixLIsDirty)
-//     matrixL();
+  //     matrixL();
-//   Base::solveInPlace(b);
+  //   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);
+
   cholmod_dense* x = cholmod_solve(CHOLMOD_A, m_cholmodFactor, &cdb, &m_cholmod);
-  if(!x)
+  ei_assert(x && "Eigen: cholmod_solve failed.");
-  {
+
-    std::cerr << "Eigen: cholmod_solve failed\n";
-    return false;
-  }
   b = Matrix<typename Base::Scalar,Dynamic,1>::Map(reinterpret_cast<typename Base::Scalar*>(x->x),b.rows());
   cholmod_free_dense(&x, &m_cholmod);
   return true;
 }
 
+
+
+
+
+
+
+
+
+
+
+template<typename _MatrixType>
+class SparseLDLT<_MatrixType,Cholmod> : public SparseLDLT<_MatrixType>
+{
+  protected:
+    typedef SparseLDLT<_MatrixType> Base;
+    typedef typename Base::Scalar Scalar;
+    typedef typename Base::RealScalar RealScalar;
+    using Base::MatrixLIsDirty;
+    using Base::SupernodalFactorIsDirty;
+    using Base::m_flags;
+    using Base::m_matrix;
+    using Base::m_status;
+
+  public:
+    typedef _MatrixType MatrixType;
+    typedef typename MatrixType::Index Index;
+
+    SparseLDLT(int flags = 0)
+      : Base(flags), m_cholmodFactor(0)
+    {
+      cholmod_start(&m_cholmod);
+    }
+
+    SparseLDLT(const _MatrixType& matrix, int flags = 0)
+      : Base(flags), m_cholmodFactor(0)
+    {
+      cholmod_start(&m_cholmod);
+      compute(matrix);
+    }
+
+    ~SparseLDLT()
+    {
+      if (m_cholmodFactor)
+        cholmod_free_factor(&m_cholmodFactor, &m_cholmod);
+      cholmod_finish(&m_cholmod);
+    }
+
+    inline const typename Base::CholMatrixType& matrixL(void) const;
+
+    template<typename Derived>
+    void solveInPlace(MatrixBase<Derived> &b) const;
+
+    template<typename Rhs>
+    inline const ei_solve_retval<SparseLDLT<MatrixType, Cholmod>, Rhs>
+    solve(const MatrixBase<Rhs>& b) const
+    {
+      ei_assert(true && "SparseLDLT is not initialized.");
+      return ei_solve_retval<SparseLDLT<MatrixType, Cholmod>, Rhs>(*this, b.derived());
+    }
+
+    void compute(const _MatrixType& matrix);
+
+    inline Index cols() const { return m_matrix.cols(); }
+    inline Index rows() const { return m_matrix.rows(); }
+
+    inline const cholmod_factor* cholmodFactor() const
+    { return m_cholmodFactor; }
+
+    inline cholmod_common* cholmodCommon() const
+    { return &m_cholmod; }
+
+    bool succeeded() const;
+
+  protected:
+    mutable cholmod_common m_cholmod;
+    cholmod_factor* m_cholmodFactor;
+};
+
+
+
+
+
+template<typename _MatrixType, typename Rhs>
+  struct ei_solve_retval<SparseLDLT<_MatrixType, Cholmod>, Rhs>
+  : ei_solve_retval_base<SparseLDLT<_MatrixType, Cholmod>, Rhs>
+{
+  typedef SparseLDLT<_MatrixType, Cholmod> SpLDLTDecType;
+  EIGEN_MAKE_SOLVE_HELPERS(SpLDLTDecType,Rhs)
+
+  template<typename Dest> void evalTo(Dest& dst) const
+  {
+    //Index size = dec().cholmodFactor()->n;
+    ei_assert((Index)dec().cholmodFactor()->n==rhs().rows());
+    
+    cholmod_factor* cholmodFactor = const_cast<cholmod_factor*>(dec().cholmodFactor());
+    cholmod_common* cholmodCommon = const_cast<cholmod_common*>(dec().cholmodCommon());
+    // 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(rhs().const_cast_derived());
+    cholmod_dense* x = cholmod_solve(CHOLMOD_LDLt, cholmodFactor, &cdb, cholmodCommon);
+
+    dst = Matrix<typename Base::Scalar,Dynamic,1>::Map(reinterpret_cast<typename Base::Scalar*>(x->x), rhs().rows());  
+    cholmod_free_dense(&x, cholmodCommon);
+
+  }
+    
+};
+
+
+
+
+
+template<typename _MatrixType>
+void SparseLDLT<_MatrixType,Cholmod>::compute(const _MatrixType& a)
+{
+  if (m_cholmodFactor)
+  {
+    cholmod_free_factor(&m_cholmodFactor, &m_cholmod);
+    m_cholmodFactor = 0;
+  }
+
+  cholmod_sparse A = ei_cholmod_map_eigen_to_sparse(const_cast<_MatrixType&>(a));
+ 
+  //m_cholmod.supernodal = CHOLMOD_AUTO;
+  m_cholmod.supernodal = CHOLMOD_SIMPLICIAL;
+  //m_cholmod.supernodal = CHOLMOD_SUPERNODAL;
+  // TODO
+  if (m_flags & IncompleteFactorization)
+  {
+    m_cholmod.nmethods = 1;
+    //m_cholmod.method[0].ordering = CHOLMOD_NATURAL;
+    m_cholmod.method[0].ordering = CHOLMOD_COLAMD;
+    m_cholmod.postorder = 1;
+  }
+  else
+  {
+    m_cholmod.nmethods = 1;
+    m_cholmod.method[0].ordering = CHOLMOD_NATURAL;
+    m_cholmod.postorder = 0;
+  }
+  m_cholmod.final_ll = 0;
+  m_cholmodFactor = cholmod_analyze(&A, &m_cholmod);
+  cholmod_factorize(&A, m_cholmodFactor, &m_cholmod);
+  
+  m_status = (m_status & ~SupernodalFactorIsDirty) | MatrixLIsDirty;
+}
+
+
+// TODO
+template<typename _MatrixType>
+bool SparseLDLT<_MatrixType,Cholmod>::succeeded() const
+{ return true; }
+
+
+template<typename _MatrixType>
+inline const typename SparseLDLT<_MatrixType>::CholMatrixType&
+SparseLDLT<_MatrixType,Cholmod>::matrixL() const
+{
+  if (m_status & MatrixLIsDirty)
+  {
+    ei_assert(!(m_status & SupernodalFactorIsDirty));
+
+    cholmod_sparse* cmRes = cholmod_factor_to_sparse(m_cholmodFactor, &m_cholmod);
+    const_cast<typename Base::CholMatrixType&>(m_matrix) = MappedSparseMatrix<Scalar>(*cmRes);
+    free(cmRes);
+
+    m_status = (m_status & ~MatrixLIsDirty);
+  }
+  return m_matrix;
+}
+
+
+
+
+
+
+template<typename _MatrixType>
+template<typename Derived>
+void SparseLDLT<_MatrixType,Cholmod>::solveInPlace(MatrixBase<Derived> &b) const
+{
+  //Index size = m_cholmodFactor->n;
+  ei_assert((Index)m_cholmodFactor->n == 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_A, 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

unsupported/Eigen/src/SparseExtra/SparseLDLT.h

@@ -75,15 +75,14 @@ LDL License:
   *
   * \sa class LDLT, class LDLT
   */
-template<typename MatrixType, int Backend = DefaultBackend>
+template<typename _MatrixType, typename Backend = DefaultBackend>
 class SparseLDLT
 {
   protected:
-    typedef typename MatrixType::Scalar Scalar;
+    typedef typename _MatrixType::Scalar Scalar;
-    typedef typename MatrixType::Index Index;
+    typedef typename NumTraits<typename _MatrixType::Scalar>::Real RealScalar;
-    typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar;
+    
-    typedef SparseMatrix<Scalar> CholMatrixType;
+    typedef Matrix<Scalar,_MatrixType::ColsAtCompileTime,1> VectorType;
-    typedef Matrix<Scalar,MatrixType::ColsAtCompileTime,1> VectorType;
 
     enum {
       SupernodalFactorIsDirty      = 0x10000,
@@ -91,6 +90,10 @@ class SparseLDLT
     };
 
   public:
+    typedef SparseMatrix<Scalar> CholMatrixType;
+    typedef _MatrixType MatrixType;
+    typedef typename MatrixType::Index Index;
+
 
     /** Creates a dummy LDLT factorization object with flags \a flags. */
     SparseLDLT(int flags = 0)
@@ -162,6 +165,19 @@ class SparseLDLT
     template<typename Derived>
     bool solveInPlace(MatrixBase<Derived> &b) const;
 
+    template<typename Rhs>
+    inline const ei_solve_retval<SparseLDLT<MatrixType>, Rhs>
+    solve(const MatrixBase<Rhs>& b) const
+    {
+      ei_assert(true && "SparseLDLT is not initialized.");
+      return ei_solve_retval<SparseLDLT<MatrixType>, Rhs>(*this, b.derived());
+    }
+
+    inline Index cols() const { return m_matrix.cols(); }
+    inline Index rows() const { return m_matrix.rows(); }
+
+    inline const VectorType& diag() const { return m_diag; }
+
     /** \returns true if the factorization succeeded */
     inline bool succeeded(void) const { return m_succeeded; }
 
@@ -177,18 +193,52 @@ class SparseLDLT
     bool m_succeeded;
 };
 
+
+
+
+
+template<typename _MatrixType, typename Rhs>
+struct ei_solve_retval<SparseLDLT<_MatrixType>, Rhs>
+  : ei_solve_retval_base<SparseLDLT<_MatrixType>, Rhs>
+{
+  typedef SparseLDLT<_MatrixType> SpLDLTDecType;
+  EIGEN_MAKE_SOLVE_HELPERS(SpLDLTDecType,Rhs)
+
+  template<typename Dest> void evalTo(Dest& dst) const
+  {
+    //Index size = dec().matrixL().rows();
+    ei_assert(dec().matrixL().rows()==rhs().rows());
+
+    Rhs b(rhs().rows(), rhs().cols());
+    b = rhs();
+
+    if (dec().matrixL().nonZeros()>0) // otherwise L==I
+      dec().matrixL().template triangularView<UnitLower>().solveInPlace(b);
+
+    b = b.cwiseQuotient(dec().diag());
+    if (dec().matrixL().nonZeros()>0) // otherwise L==I
+      dec().matrixL().adjoint().template triangularView<UnitUpper>().solveInPlace(b);
+    
+    dst = b;
+
+  }
+    
+};
+
+
+
 /** Computes / recomputes the LDLT decomposition of matrix \a a
   * using the default algorithm.
   */
-template<typename MatrixType, int Backend>
+template<typename _MatrixType, typename Backend>
-void SparseLDLT<MatrixType,Backend>::compute(const MatrixType& a)
+void SparseLDLT<_MatrixType,Backend>::compute(const _MatrixType& a)
 {
   _symbolic(a);
   m_succeeded = _numeric(a);
 }
 
-template<typename MatrixType, int Backend>
+template<typename _MatrixType, typename Backend>
-void SparseLDLT<MatrixType,Backend>::_symbolic(const MatrixType& a)
+void SparseLDLT<_MatrixType,Backend>::_symbolic(const _MatrixType& a)
 {
   assert(a.rows()==a.cols());
   const Index size = a.rows();
@@ -244,8 +294,8 @@ void SparseLDLT<MatrixType,Backend>::_symbolic(const MatrixType& a)
   ei_aligned_stack_delete(Index, tags, size);
 }
 
-template<typename MatrixType, int Backend>
+template<typename _MatrixType, typename Backend>
-bool SparseLDLT<MatrixType,Backend>::_numeric(const MatrixType& a)
+bool SparseLDLT<_MatrixType,Backend>::_numeric(const _MatrixType& a)
 {
   assert(a.rows()==a.cols());
   const Index size = a.rows();
@@ -327,12 +377,12 @@ bool SparseLDLT<MatrixType,Backend>::_numeric(const MatrixType& a)
 }
 
 /** Computes b = L^-T D^-1 L^-1 b */
-template<typename MatrixType, int Backend>
+template<typename _MatrixType, typename Backend>
 template<typename Derived>
-bool SparseLDLT<MatrixType, Backend>::solveInPlace(MatrixBase<Derived> &b) const
+bool SparseLDLT<_MatrixType, Backend>::solveInPlace(MatrixBase<Derived> &b) const
 {
-  const Index size = m_matrix.rows();
+  //Index size = m_matrix.rows();
-  ei_assert(size==b.rows());
+  ei_assert(m_matrix.rows()==b.rows());
   if (!m_succeeded)
     return false;
 

unsupported/Eigen/src/SparseExtra/SparseLLT.h

@@ -35,14 +35,12 @@
   *
   * \sa class LLT, class LDLT
   */
-template<typename MatrixType, int Backend = DefaultBackend>
+template<typename _MatrixType, typename Backend = DefaultBackend>
 class SparseLLT
 {
   protected:
-    typedef typename MatrixType::Scalar Scalar;
+    typedef typename _MatrixType::Scalar Scalar;
-    typedef typename MatrixType::Index Index;
+    typedef typename NumTraits<typename _MatrixType::Scalar>::Real RealScalar;
-    typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar;
-    typedef SparseMatrix<Scalar> CholMatrixType;
 
     enum {
       SupernodalFactorIsDirty      = 0x10000,
@@ -50,6 +48,9 @@ class SparseLLT
     };
 
   public:
+    typedef SparseMatrix<Scalar> CholMatrixType;
+    typedef _MatrixType MatrixType;
+    typedef typename MatrixType::Index Index;
 
     /** Creates a dummy LLT factorization object with flags \a flags. */
     SparseLLT(int flags = 0)
@@ -110,6 +111,17 @@ class SparseLLT
     template<typename Derived>
     bool solveInPlace(MatrixBase<Derived> &b) const;
 
+    template<typename Rhs>
+    inline const ei_solve_retval<SparseLLT<MatrixType>, Rhs>
+    solve(const MatrixBase<Rhs>& b) const
+    {
+      ei_assert(true && "SparseLLT is not initialized.");
+      return ei_solve_retval<SparseLLT<MatrixType>, Rhs>(*this, b.derived());
+    }
+
+    inline Index cols() const { return m_matrix.cols(); }
+    inline Index rows() const { return m_matrix.rows(); }
+
     /** \returns true if the factorization succeeded */
     inline bool succeeded(void) const { return m_succeeded; }
 
@@ -121,11 +133,43 @@ class SparseLLT
     bool m_succeeded;
 };
 
+
+
+
+
+
+template<typename _MatrixType, typename Rhs>
+struct ei_solve_retval<SparseLLT<_MatrixType>, Rhs>
+  : ei_solve_retval_base<SparseLLT<_MatrixType>, Rhs>
+{
+  typedef SparseLLT<_MatrixType> SpLLTDecType;
+  EIGEN_MAKE_SOLVE_HELPERS(SpLLTDecType,Rhs)
+
+  template<typename Dest> void evalTo(Dest& dst) const
+  {
+    const Index size = dec().matrixL().rows();
+    ei_assert(size==rhs().rows());
+    
+    Rhs b(rhs().rows(), rhs().cols());
+    b = rhs();
+    
+    dec().matrixL().template triangularView<Lower>().solveInPlace(b);
+    dec().matrixL().adjoint().template triangularView<Upper>().solveInPlace(b);
+    
+    dst = b;
+    
+  }
+    
+};
+
+
+
+
 /** Computes / recomputes the LLT decomposition of matrix \a a
   * using the default algorithm.
   */
-template<typename MatrixType, int Backend>
+template<typename _MatrixType, typename Backend>
-void SparseLLT<MatrixType,Backend>::compute(const MatrixType& a)
+void SparseLLT<_MatrixType,Backend>::compute(const _MatrixType& a)
 {
   assert(a.rows()==a.cols());
   const Index size = a.rows();
@@ -148,7 +192,7 @@ void SparseLLT<MatrixType,Backend>::compute(const MatrixType& a)
     tempVector.setZero();
     // init with current matrix a
     {
-      typename MatrixType::InnerIterator it(a,j);
+      typename _MatrixType::InnerIterator it(a,j);
       ei_assert(it.index()==j &&
         "matrix must has non zero diagonal entries and only the lower triangular part must be stored");
       ++it; // skip diagonal element
@@ -187,9 +231,9 @@ void SparseLLT<MatrixType,Backend>::compute(const MatrixType& a)
 }
 
 /** Computes b = L^-T L^-1 b */
-template<typename MatrixType, int Backend>
+template<typename _MatrixType, typename Backend>
 template<typename Derived>
-bool SparseLLT<MatrixType, Backend>::solveInPlace(MatrixBase<Derived> &b) const
+bool SparseLLT<_MatrixType, Backend>::solveInPlace(MatrixBase<Derived> &b) const
 {
   const Index size = m_matrix.rows();
   ei_assert(size==b.rows());

unsupported/Eigen/src/SparseExtra/SparseLU.h

@@ -37,16 +37,16 @@ enum {
   *
   * \brief LU decomposition of a sparse matrix and associated features
   *
-  * \param MatrixType the type of the matrix of which we are computing the LU factorization
+  * \param _MatrixType the type of the matrix of which we are computing the LU factorization
   *
   * \sa class FullPivLU, class SparseLLT
   */
-template<typename MatrixType, int Backend = DefaultBackend>
+template<typename _MatrixType, typename Backend = DefaultBackend>
 class SparseLU
-{
+  {
   protected:
-    typedef typename MatrixType::Scalar Scalar;
+    typedef typename _MatrixType::Scalar Scalar;
-    typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar;
+    typedef typename NumTraits<typename _MatrixType::Scalar>::Real RealScalar;
     typedef SparseMatrix<Scalar> LUMatrixType;
 
     enum {
@@ -54,6 +54,7 @@ class SparseLU
     };
 
   public:
+    typedef _MatrixType MatrixType;
 
     /** Creates a dummy LU factorization object with flags \a flags. */
     SparseLU(int flags = 0)
@@ -64,7 +65,7 @@ class SparseLU
 
     /** Creates a LU object and compute the respective factorization of \a matrix using
       * flags \a flags. */
-    SparseLU(const MatrixType& matrix, int flags = 0)
+    SparseLU(const _MatrixType& matrix, int flags = 0)
       : /*m_matrix(matrix.rows(), matrix.cols()),*/ m_flags(flags), m_status(0)
     {
       m_precision = RealScalar(0.1) * Eigen::NumTraits<RealScalar>::dummy_precision();
@@ -112,13 +113,13 @@ class SparseLU
     }
 
     /** Computes/re-computes the LU factorization */
-    void compute(const MatrixType& matrix);
+    void compute(const _MatrixType& matrix);
 
     /** \returns the lower triangular matrix L */
-    //inline const MatrixType& matrixL() const { return m_matrixL; }
+    //inline const _MatrixType& matrixL() const { return m_matrixL; }
 
     /** \returns the upper triangular matrix U */
-    //inline const MatrixType& matrixU() const { return m_matrixU; }
+    //inline const _MatrixType& matrixU() const { return m_matrixU; }
 
     template<typename BDerived, typename XDerived>
     bool solve(const MatrixBase<BDerived> &b, MatrixBase<XDerived>* x,
@@ -137,8 +138,8 @@ class SparseLU
 /** Computes / recomputes the LU decomposition of matrix \a a
   * using the default algorithm.
   */
-template<typename MatrixType, int Backend>
+template<typename _MatrixType, typename Backend>
-void SparseLU<MatrixType,Backend>::compute(const MatrixType& )
+void SparseLU<_MatrixType,Backend>::compute(const _MatrixType& )
 {
   ei_assert(false && "not implemented yet");
 }
@@ -151,9 +152,9 @@ void SparseLU<MatrixType,Backend>::compute(const MatrixType& )
   * Not all backends implement the solution of the transposed or
   * adjoint system.
   */
-template<typename MatrixType, int Backend>
+template<typename _MatrixType, typename Backend>
 template<typename BDerived, typename XDerived>
-bool SparseLU<MatrixType,Backend>::solve(const MatrixBase<BDerived> &, MatrixBase<XDerived>* , const int ) const
+bool SparseLU<_MatrixType,Backend>::solve(const MatrixBase<BDerived> &, MatrixBase<XDerived>* , const int ) const
 {
   ei_assert(false && "not implemented yet");
   return false;

unsupported/Eigen/src/SparseExtra/UmfPackSupport.h

@@ -117,22 +117,24 @@ inline int umfpack_get_determinant(std::complex<double> *Mx, double *Ex, void *N
 }
 
 
-template<typename MatrixType>
+template<typename _MatrixType>
-class SparseLU<MatrixType,UmfPack> : public SparseLU<MatrixType>
+class SparseLU<_MatrixType,UmfPack> : public SparseLU<_MatrixType>
 {
   protected:
-    typedef SparseLU<MatrixType> Base;
+    typedef SparseLU<_MatrixType> Base;
     typedef typename Base::Scalar Scalar;
     typedef typename Base::RealScalar RealScalar;
     typedef Matrix<Scalar,Dynamic,1> Vector;
-    typedef Matrix<int, 1, MatrixType::ColsAtCompileTime> IntRowVectorType;
+    typedef Matrix<int, 1, _MatrixType::ColsAtCompileTime> IntRowVectorType;
-    typedef Matrix<int, MatrixType::RowsAtCompileTime, 1> IntColVectorType;
+    typedef Matrix<int, _MatrixType::RowsAtCompileTime, 1> IntColVectorType;
     typedef SparseMatrix<Scalar,Lower|UnitDiag> LMatrixType;
     typedef SparseMatrix<Scalar,Upper> UMatrixType;
     using Base::m_flags;
     using Base::m_status;
 
   public:
+    typedef _MatrixType MatrixType;
+    typedef typename MatrixType::Index Index;
 
     SparseLU(int flags = NaturalOrdering)
       : Base(flags), m_numeric(0)
@@ -180,8 +182,30 @@ class SparseLU<MatrixType,UmfPack> : public SparseLU<MatrixType>
     template<typename BDerived, typename XDerived>
     bool solve(const MatrixBase<BDerived> &b, MatrixBase<XDerived>* x) const;
 
+    template<typename Rhs>
+      inline const ei_solve_retval<SparseLU<MatrixType, UmfPack>, Rhs>
+    solve(const MatrixBase<Rhs>& b) const
+    {
+      ei_assert(true && "SparseLU is not initialized.");
+      return ei_solve_retval<SparseLU<MatrixType, UmfPack>, Rhs>(*this, b.derived());
+    }
+
     void compute(const MatrixType& matrix);
 
+    inline Index cols() const { return m_matrixRef->cols(); }
+    inline Index rows() const { return m_matrixRef->rows(); }
+
+    inline const MatrixType& matrixLU() const
+    {
+      //ei_assert(m_isInitialized && "LU is not initialized.");
+      return *m_matrixRef;
+    }
+
+    const void* numeric() const
+    {
+      return m_numeric;
+    }
+
   protected:
 
     void extractData() const;
@@ -197,6 +221,37 @@ class SparseLU<MatrixType,UmfPack> : public SparseLU<MatrixType>
     mutable bool m_extractedDataAreDirty;
 };
 
+
+template<typename _MatrixType, typename Rhs>
+  struct ei_solve_retval<SparseLU<_MatrixType, UmfPack>, Rhs>
+  : ei_solve_retval_base<SparseLU<_MatrixType, UmfPack>, Rhs>
+{
+  typedef SparseLU<_MatrixType, UmfPack> SpLUDecType;
+  EIGEN_MAKE_SOLVE_HELPERS(SpLUDecType,Rhs)
+
+  template<typename Dest> void evalTo(Dest& dst) const
+  {
+    const int rhsCols = rhs().cols();
+
+    ei_assert((Rhs::Flags&RowMajorBit)==0 && "UmfPack backend does not support non col-major rhs yet");
+    ei_assert((Dest::Flags&RowMajorBit)==0 && "UmfPack backend does not support non col-major result yet");
+
+    void* numeric = const_cast<void*>(dec().numeric());
+
+    int errorCode = 0;
+    for (int j=0; j<rhsCols; ++j)
+      {
+ errorCode = umfpack_solve(UMFPACK_A,
+         dec().matrixLU()._outerIndexPtr(), dec().matrixLU()._innerIndexPtr(), dec().matrixLU()._valuePtr(),
+         &dst.col(j).coeffRef(0), &rhs().const_cast_derived().col(j).coeffRef(0), numeric, 0, 0);
+ ei_assert(!errorCode && "UmfPack could not solve the system.");
+      }
+  }
+    
+};
+
+
+
 template<typename MatrixType>
 void SparseLU<MatrixType,UmfPack>::compute(const MatrixType& a)
 {

unsupported/test/sparse_ldlt.cpp

@@ -25,6 +25,10 @@
 #include "sparse.h"
 #include <Eigen/SparseExtra>
 
+#ifdef EIGEN_CHOLMOD_SUPPORT
+#include <Eigen/CholmodSupport>
+#endif
+
 #ifdef EIGEN_TAUCS_SUPPORT
 #include <Eigen/TaucsSupport>
 #endif
@@ -56,6 +60,29 @@ template<typename Scalar> void sparse_ldlt(int rows, int cols)
 
   VERIFY_IS_APPROX(refMat2.template selfadjointView<Upper>() * x, b);
   VERIFY(refX.isApprox(x,test_precision<Scalar>()) && "LDLT: default");
+
+#ifdef EIGEN_CHOLMOD_SUPPORT
+  x = b;
+  SparseLDLT<SparseSelfAdjointMatrix, Cholmod> ldlt2(m2);
+  if (ldlt2.succeeded())
+    ldlt2.solveInPlace(x);
+  else
+    std::cerr << "warning LDLT failed\n";
+
+  VERIFY_IS_APPROX(refMat2.template selfadjointView<Upper>() * x, b);
+  VERIFY(refX.isApprox(x,test_precision<Scalar>()) && "LDLT: cholmod solveInPlace");
+
+
+  SparseLDLT<SparseSelfAdjointMatrix, Cholmod> ldlt3(m2);
+  if (ldlt3.succeeded())
+    x = ldlt3.solve(b);
+  else
+    std::cerr << "warning LDLT failed\n";
+
+  VERIFY_IS_APPROX(refMat2.template selfadjointView<Upper>() * x, b);
+  VERIFY(refX.isApprox(x,test_precision<Scalar>()) && "LDLT: cholmod solve");
+
+#endif
 }
 
 void test_sparse_ldlt()

unsupported/test/sparse_llt.cpp

@@ -47,6 +47,7 @@ template<typename Scalar> void sparse_llt(int rows, int cols)
     DenseVector refX(cols), x(cols);
 
     initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag|MakeLowerTriangular, 0, 0);
+
     for(int i=0; i<rows; ++i)
       m2.coeffRef(i,i) = refMat2(i,i) = ei_abs(ei_real(refMat2(i,i)));
 
@@ -57,11 +58,24 @@ template<typename Scalar> void sparse_llt(int rows, int cols)
       SparseLLT<SparseMatrix<Scalar> > (m2).solveInPlace(x);
       VERIFY(refX.isApprox(x,test_precision<Scalar>()) && "LLT: default");
     }
-    #ifdef EIGEN_CHOLMOD_SUPPORT
+        
-    x = b;
+#ifdef EIGEN_CHOLMOD_SUPPORT
-    SparseLLT<SparseMatrix<Scalar> ,Cholmod>(m2).solveInPlace(x);
+    {
-    VERIFY(refX.isApprox(x,test_precision<Scalar>()) && "LLT: cholmod");
+      // Cholmod, as configured in CholmodSupport.h, only supports self-adjoint matrices
-    #endif
+      SparseMatrix<Scalar> m3 = m2.adjoint()*m2;
+      DenseMatrix refMat3 = refMat2.adjoint()*refMat2;
+      
+      refX = refMat3.template selfadjointView<Lower>().llt().solve(b);
+      
+      x = b;
+      SparseLLT<SparseMatrix<Scalar>, Cholmod>(m3).solveInPlace(x);
+      VERIFY((m3*x).isApprox(b,test_precision<Scalar>()) && "LLT: cholmod solveInPlace");
+      
+      x = SparseLLT<SparseMatrix<Scalar>, Cholmod>(m3).solve(b);
+      VERIFY(refX.isApprox(x,test_precision<Scalar>()) && "LLT: cholmod solve");
+    }
+#endif
+
 
     #ifdef EIGEN_TAUCS_SUPPORT
     // TODO fix TAUCS with complexes
@@ -72,10 +86,10 @@ template<typename Scalar> void sparse_llt(int rows, int cols)
 //       VERIFY(refX.isApprox(x,test_precision<Scalar>()) && "LLT: taucs (IncompleteFactorization)");
 
       x = b;
-      SparseLLT<SparseMatrix<Scalar> ,Taucs>(m2,SupernodalMultifrontal).solveInPlace(x);
+      SparseLLT<SparseMatrix<Scalar>, Taucs>(m2,SupernodalMultifrontal).solveInPlace(x);
       VERIFY(refX.isApprox(x,test_precision<Scalar>()) && "LLT: taucs (SupernodalMultifrontal)");
       x = b;
-      SparseLLT<SparseMatrix<Scalar> ,Taucs>(m2,SupernodalLeftLooking).solveInPlace(x);
+      SparseLLT<SparseMatrix<Scalar>, Taucs>(m2,SupernodalLeftLooking).solveInPlace(x);
       VERIFY(refX.isApprox(x,test_precision<Scalar>()) && "LLT: taucs (SupernodalLeftLooking)");
     }
     #endif