Sparse module: add experimental support for TAUCS and CHOLMOD with:

 * bidirectionnal mapping
 * full cholesky factorization

6 files changed, 312 insertions, 18 deletions

diff --git a/Eigen/Sparse b/Eigen/Sparse
index e5126a2d1..d1418943b 100644
--- a/Eigen/Sparse
+++ b/Eigen/Sparse
@@ -8,6 +8,27 @@
 #include <cstring>
 #include <algorithm>
 
+#ifdef EIGEN_CHOLMOD_SUPPORT
+  extern "C" {
+    #include "cholmod.h"
+  }
+#endif
+
+#ifdef EIGEN_TAUCS_SUPPORT
+
+  extern "C" {
+    #include "taucs.h"
+  }
+
+  #ifdef min
+    #undef min
+  #endif
+  #ifdef max
+    #undef max
+  #endif
+
+#endif
+
 namespace Eigen {
 
 #include "src/Sparse/SparseUtil.h"
@@ -22,7 +43,15 @@ namespace Eigen {
 #include "src/Sparse/SparseSetter.h"
 #include "src/Sparse/SparseProduct.h"
 #include "src/Sparse/TriangularSolver.h"
-#include "src/Sparse/BasicSparseCholesky.h"
+#include "src/Sparse/SparseCholesky.h"
+
+#ifdef EIGEN_CHOLMOD_SUPPORT
+# include "src/Sparse/CholmodSupport.h"
+#endif
+
+#ifdef EIGEN_TAUCS_SUPPORT
+# include "src/Sparse/TaucsSupport.h"
+#endif
 
 } // namespace Eigen
 
diff --git a/Eigen/src/Sparse/CholmodSupport.h b/Eigen/src/Sparse/CholmodSupport.h
new file mode 100644
index 000000000..f451dbcc0
--- /dev/null
+++ b/Eigen/src/Sparse/CholmodSupport.h
@@ -0,0 +1,123 @@
+// 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_CHOLMODSUPPORT_H
+#define EIGEN_CHOLMODSUPPORT_H
+
+template<typename Scalar, int Flags>
+cholmod_sparse SparseMatrix<Scalar,Flags>::asCholmodMatrix()
+{
+  cholmod_sparse res;
+  res.nzmax   = nonZeros();
+  res.nrow    = rows();;
+  res.ncol    = cols();
+  res.p       = _outerIndexPtr();
+  res.i       = _innerIndexPtr();
+  res.x       = _valuePtr();
+  res.xtype = CHOLMOD_REAL;
+  res.itype = CHOLMOD_INT;
+  res.sorted = 1;
+  res.packed = 1;
+  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");
+  }
+
+  if (Flags & SelfAdjoint)
+  {
+    if (Flags & Upper)
+      res.stype = 1;
+    else if (Flags & Lower)
+      res.stype = -1;
+    else
+      res.stype = 0;
+  }
+  else
+    res.stype = 0;
+
+  return res;
+}
+
+template<typename Scalar, int Flags>
+SparseMatrix<Scalar,Flags> SparseMatrix<Scalar,Flags>::Map(cholmod_sparse& cm)
+{
+  SparseMatrix res;
+  res.m_innerSize = cm.nrow;
+  res.m_outerSize = cm.ncol;
+  res.m_outerIndex = reinterpret_cast<int*>(cm.p);
+  SparseArray<Scalar> data = SparseArray<Scalar>::Map(
+                                reinterpret_cast<int*>(cm.i),
+                                reinterpret_cast<Scalar*>(cm.x),
+                                res.m_outerIndex[cm.ncol]);
+  res.m_data.swap(data);
+  // res.markAsRValue();
+  return res;
+}
+
+template<typename MatrixType>
+void SparseCholesky<MatrixType>::computeUsingCholmod(const MatrixType& a)
+{
+  cholmod_common c;
+  cholmod_start(&c);
+  cholmod_sparse A = const_cast<MatrixType&>(a).asCholmodMatrix();
+  std::vector<int> perm(a.cols());
+  for (int i=0; i<a.cols(); ++i)
+    perm[i] = i;
+  c.nmethods = 1;
+  c.method [0].ordering = CHOLMOD_NATURAL;
+  c.postorder = 0;
+  c.final_ll = 1;
+  cholmod_factor *L = cholmod_analyze_p(&A, &perm[0], &perm[0], a.cols(), &c);
+  cholmod_factorize(&A, L, &c);
+  cholmod_sparse* cmRes = cholmod_factor_to_sparse(L, &c);
+  m_matrix = CholMatrixType::Map(*cmRes);
+  free(cmRes);
+  cholmod_free_factor(&L, &c);
+  cholmod_finish(&c);
+}
+
+#endif // EIGEN_CHOLMODSUPPORT_H
diff --git a/Eigen/src/Sparse/SparseArray.h b/Eigen/src/Sparse/SparseArray.h
index 0315c93d9..f47ad39a9 100644
--- a/Eigen/src/Sparse/SparseArray.h
+++ b/Eigen/src/Sparse/SparseArray.h
@@ -28,7 +28,8 @@
 /** Stores a sparse set of values as a list of values and a list of indices.
   *
   */
-template<typename Scalar> class SparseArray
+template<typename Scalar>
+class SparseArray
 {
   public:
     SparseArray()
@@ -105,6 +106,15 @@ template<typename Scalar> class SparseArray
     int& index(int i) { return m_indices[i]; }
     const int& index(int i) const { return m_indices[i]; }
 
+    static SparseArray Map(int* indices, Scalar* values, int size)
+    {
+      SparseArray res;
+      res.m_indices = indices;
+      res.m_values = values;
+      res.m_allocatedSize = res.m_size = size;
+      return res;
+    }
+
   protected:
 
     void reallocate(int size)
diff --git a/Eigen/src/Sparse/BasicSparseCholesky.h b/Eigen/src/Sparse/SparseCholesky.h
index 92193a3bc..5af84ebed 100644
--- a/Eigen/src/Sparse/BasicSparseCholesky.h
+++ b/Eigen/src/Sparse/SparseCholesky.h
@@ -22,12 +22,20 @@
 // License and a copy of the GNU General Public License along with
 // Eigen. If not, see <http://www.gnu.org/licenses/>.
 
-#ifndef EIGEN_BASICSPARSECHOLESKY_H
-#define EIGEN_BASICSPARSECHOLESKY_H
+#ifndef EIGEN_SPARSECHOLESKY_H
+#define EIGEN_SPARSECHOLESKY_H
+
+enum {
+  CholFull        = 0x0,  // full is the default
+  CholPartial     = 0x1,
+  CholUseEigen    = 0x0,  // Eigen's impl is the default
+  CholUseTaucs    = 0x2,
+  CholUseCholmod  = 0x4,
+};
 
 /** \ingroup Sparse_Module
   *
-  * \class BasicSparseCholesky
+  * \class SparseCholesky
   *
   * \brief Standard Cholesky decomposition of a matrix and associated features
   *
@@ -35,12 +43,13 @@
   *
   * \sa class Cholesky, class CholeskyWithoutSquareRoot
   */
-template<typename MatrixType> class BasicSparseCholesky
+template<typename MatrixType> class SparseCholesky
 {
   private:
     typedef typename MatrixType::Scalar Scalar;
     typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar;
     typedef Matrix<Scalar, MatrixType::ColsAtCompileTime, 1> VectorType;
+    typedef SparseMatrix<Scalar,Lower> CholMatrixType;
 
     enum {
       PacketSize = ei_packet_traits<Scalar>::size,
@@ -49,13 +58,13 @@ template<typename MatrixType> class BasicSparseCholesky
 
   public:
 
-    BasicSparseCholesky(const MatrixType& matrix)
-      : m_matrix(matrix.rows(), matrix.cols())
+    SparseCholesky(const MatrixType& matrix, int flags = 0)
+      : m_matrix(matrix.rows(), matrix.cols()), m_flags(flags)
     {
       compute(matrix);
     }
 
-    inline const MatrixType& matrixL(void) const { return m_matrix; }
+    inline const CholMatrixType& matrixL(void) const { return m_matrix; }
 
     /** \returns true if the matrix is positive definite */
     inline bool isPositiveDefinite(void) const { return m_isPositiveDefinite; }
@@ -67,25 +76,35 @@ template<typename MatrixType> class BasicSparseCholesky
     void compute(const MatrixType& matrix);
 
   protected:
+    void computeUsingEigen(const MatrixType& matrix);
+    void computeUsingTaucs(const MatrixType& matrix);
+    void computeUsingCholmod(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;
+    CholMatrixType m_matrix;
+    int m_flags;
     bool m_isPositiveDefinite;
-
-    struct ListEl
-    {
-      int next;
-      int index;
-      Scalar value;
-    };
 };
 
 /** Computes / recomputes the Cholesky decomposition A = LL^* = U^*U of \a matrix
   */
 template<typename MatrixType>
-void BasicSparseCholesky<MatrixType>::compute(const MatrixType& a)
+void SparseCholesky<MatrixType>::compute(const MatrixType& a)
+{
+  if (m_flags&CholUseTaucs)
+    computeUsingTaucs(a);
+  else if (m_flags&CholUseCholmod)
+    computeUsingCholmod(a);
+  else
+    computeUsingEigen(a);
+}
+
+template<typename MatrixType>
+void SparseCholesky<MatrixType>::computeUsingEigen(const MatrixType& a)
 {
   assert(a.rows()==a.cols());
   const int size = a.rows();
diff --git a/Eigen/src/Sparse/SparseMatrix.h b/Eigen/src/Sparse/SparseMatrix.h
index 480c92dcb..e3b224d7f 100644
--- a/Eigen/src/Sparse/SparseMatrix.h
+++ b/Eigen/src/Sparse/SparseMatrix.h
@@ -80,6 +80,15 @@ class SparseMatrix
     inline int outerSize() const { return m_outerSize; }
     inline int innerNonZeros(int j) const { return m_outerIndex[j+1]-m_outerIndex[j]; }
 
+    inline const Scalar* _valuePtr() const { return &m_data.value(0); }
+    inline Scalar* _valuePtr() { return &m_data.value(0); }
+
+    inline const int* _innerIndexPtr() const { return &m_data.index(0); }
+    inline int* _innerIndexPtr() { return &m_data.index(0); }
+
+    inline const int* _outerIndexPtr() const { return m_outerIndex; }
+    inline int* _outerIndexPtr() { return m_outerIndex; }
+
     inline Scalar coeff(int row, int col) const
     {
       const int outer = RowMajor ? row : col;
@@ -180,6 +189,10 @@ class SparseMatrix
       }
     }
 
+    inline SparseMatrix()
+      : m_outerSize(0), m_innerSize(0), m_outerIndex(0)
+    {}
+
     inline SparseMatrix(int rows, int cols)
       : m_outerSize(0), m_innerSize(0), m_outerIndex(0)
     {
@@ -248,6 +261,16 @@ class SparseMatrix
       return s;
     }
 
+    #ifdef EIGEN_TAUCS_SUPPORT
+    static SparseMatrix Map(taucs_ccs_matrix& taucsMatrix);
+    taucs_ccs_matrix asTaucsMatrix();
+    #endif
+
+    #ifdef EIGEN_CHOLMOD_SUPPORT
+    static SparseMatrix Map(cholmod_sparse& cholmodMatrix);
+    cholmod_sparse asCholmodMatrix();
+    #endif
+
     /** Destructor */
     inline ~SparseMatrix()
     {
diff --git a/Eigen/src/Sparse/TaucsSupport.h b/Eigen/src/Sparse/TaucsSupport.h
new file mode 100644
index 000000000..26bf19c5b
--- /dev/null
+++ b/Eigen/src/Sparse/TaucsSupport.h
@@ -0,0 +1,90 @@
+// 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_TAUCSSUPPORT_H
+#define EIGEN_TAUCSSUPPORT_H
+
+template<typename Scalar, int Flags>
+taucs_ccs_matrix SparseMatrix<Scalar,Flags>::asTaucsMatrix()
+{
+  taucs_ccs_matrix res;
+  res.n         = cols();
+  res.m         = rows();
+  res.flags     = 0;
+  res.colptr    = _outerIndexPtr();
+  res.rowind    = _innerIndexPtr();
+  res.values.v  = _valuePtr();
+  if (ei_is_same_type<Scalar,int>::ret)
+    res.flags |= TAUCS_INT;
+  else if (ei_is_same_type<Scalar,float>::ret)
+    res.flags |= TAUCS_SINGLE;
+  else if (ei_is_same_type<Scalar,double>::ret)
+    res.flags |= TAUCS_DOUBLE;
+  else if (ei_is_same_type<Scalar,std::complex<float> >::ret)
+    res.flags |= TAUCS_SCOMPLEX;
+  else if (ei_is_same_type<Scalar,std::complex<double> >::ret)
+    res.flags |= TAUCS_DCOMPLEX;
+  else
+  {
+    ei_assert(false && "Scalar type not supported by TAUCS");
+  }
+
+  if (Flags & Upper)
+    res.flags |= TAUCS_UPPER;
+  if (Flags & Lower)
+    res.flags |= TAUCS_LOWER;
+  if (Flags & SelfAdjoint)
+    res.flags |= (NumTraits<Scalar>::IsComplex ? TAUCS_HERMITIAN : TAUCS_SYMMETRIC);
+  else if ((Flags & Upper) || (Flags & Lower))
+    res.flags |= TAUCS_TRIANGULAR;
+
+  return res;
+}
+
+template<typename Scalar, int Flags>
+SparseMatrix<Scalar,Flags> SparseMatrix<Scalar,Flags>::Map(taucs_ccs_matrix& taucsMat)
+{
+  SparseMatrix res;
+  res.m_innerSize = taucsMat.m;
+  res.m_outerSize = taucsMat.n;
+  res.m_outerIndex = taucsMat.colptr;
+  SparseArray<Scalar> data = SparseArray<Scalar>::Map(
+                                taucsMat.rowind,
+                                reinterpret_cast<Scalar*>(taucsMat.values.v),
+                                taucsMat.colptr[taucsMat.n]);
+  res.m_data.swap(data);
+  // res.markAsRValue();
+  return res;
+}
+
+template<typename MatrixType>
+void SparseCholesky<MatrixType>::computeUsingTaucs(const MatrixType& a)
+{
+  taucs_ccs_matrix taucsMatA = const_cast<MatrixType&>(a).asTaucsMatrix();
+  taucs_ccs_matrix* taucsRes = taucs_ccs_factor_llt(&taucsMatA, 0, 0);
+  m_matrix = CholMatrixType::Map(*taucsRes);
+  free(taucsRes);
+}
+
+#endif // EIGEN_TAUCSSUPPORT_H