Sparse module:

* add a MappedSparseMatrix class (like Eigen::Map but for sparse
  matrices)
* rename SparseArray to CompressedStorage

Eigen/Sparse

@@ -72,11 +72,12 @@ namespace Eigen {
 
 #include "src/Sparse/SparseUtil.h"
 #include "src/Sparse/SparseMatrixBase.h"
-#include "src/Sparse/SparseArray.h"
+#include "src/Sparse/CompressedStorage.h"
 #include "src/Sparse/AmbiVector.h"
 #include "src/Sparse/RandomSetter.h"
 #include "src/Sparse/SparseBlock.h"
 #include "src/Sparse/SparseMatrix.h"
+#include "src/Sparse/MappedSparseMatrix.h"
 #include "src/Sparse/SparseVector.h"
 #include "src/Sparse/CoreIterators.h"
 #include "src/Sparse/SparseTranspose.h"

Eigen/src/Sparse/CholmodSupport.h

@@ -55,7 +55,7 @@ void ei_cholmod_configure_matrix(CholmodType& mat)
 }
 
 template<typename Scalar, int Flags>
-cholmod_sparse SparseMatrix<Scalar,Flags>::asCholmodMatrix()
+cholmod_sparse SparseMatrixBase<Scalar,Flags>::asCholmodMatrix()
 {
   cholmod_sparse res;
   res.nzmax   = nonZeros();
@@ -108,19 +108,14 @@ cholmod_dense ei_cholmod_map_eigen_to_dense(MatrixBase<Derived>& mat)
 }
 
 template<typename Scalar, int Flags>
-SparseMatrix<Scalar,Flags> SparseMatrix<Scalar,Flags>::Map(cholmod_sparse& cm)
+MappedSparseMatrix<Scalar,Flags>::MappedSparseMatrix(taucs_ccs_matrix& taucsMat)
 {
-  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;
+  m_innerSize = cm.nrow;
+  m_outerSize = cm.ncol;
+  m_outerIndex = reinterpret_cast<int*>(cm.p);
+  m_innerIndices = reinterpret_cast<int*>(cm.i);
+  m_values = reinterpret_cast<Scalar*>(cm.x);
+  m_nnz = res.m_outerIndex[cm.ncol]);
 }
 
 template<typename MatrixType>

Eigen/src/Sparse/CompressedStorage.h

@@ -22,32 +22,32 @@
 // License and a copy of the GNU General Public License along with
 // Eigen. If not, see <http://www.gnu.org/licenses/>.
 
-#ifndef EIGEN_SPARSE_ARRAY_H
-#define EIGEN_SPARSE_ARRAY_H
+#ifndef EIGEN_COMPRESSED_STORAGE_H
+#define EIGEN_COMPRESSED_STORAGE_H
 
 /** Stores a sparse set of values as a list of values and a list of indices.
   *
   */
 template<typename Scalar>
-class SparseArray
+class CompressedStorage
 {
   public:
-    SparseArray()
+    CompressedStorage()
       : m_values(0), m_indices(0), m_size(0), m_allocatedSize(0)
     {}
 
-    SparseArray(int size)
+    CompressedStorage(int size)
       : m_values(0), m_indices(0), m_size(0), m_allocatedSize(0)
     {
       resize(size);
     }
 
-    SparseArray(const SparseArray& other)
+    CompressedStorage(const CompressedStorage& other)
     {
       *this = other;
     }
 
-    SparseArray& operator=(const SparseArray& other)
+    CompressedStorage& operator=(const CompressedStorage& other)
     {
       resize(other.size());
       memcpy(m_values, other.m_values, m_size * sizeof(Scalar));
@@ -55,7 +55,7 @@ class SparseArray
       return *this;
     }
 
-    void swap(SparseArray& other)
+    void swap(CompressedStorage& other)
     {
       std::swap(m_values, other.m_values);
       std::swap(m_indices, other.m_indices);
@@ -63,7 +63,7 @@ class SparseArray
       std::swap(m_allocatedSize, other.m_allocatedSize);
     }
 
-    ~SparseArray()
+    ~CompressedStorage()
     {
       delete[] m_values;
       delete[] m_indices;
@@ -106,9 +106,9 @@ 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)
+    static CompressedStorage Map(int* indices, Scalar* values, int size)
     {
-      SparseArray res;
+      CompressedStorage res;
       res.m_indices = indices;
       res.m_values = values;
       res.m_allocatedSize = res.m_size = size;
@@ -141,4 +141,4 @@ class SparseArray
 
 };
 
-#endif // EIGEN_SPARSE_ARRAY_H
+#endif // EIGEN_COMPRESSED_STORAGE_H

Eigen/src/Sparse/MappedSparseMatrix.h

@@ -0,0 +1,168 @@
+// 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_MAPPED_SPARSEMATRIX_H
+#define EIGEN_MAPPED_SPARSEMATRIX_H
+
+/** \class MappedSparseMatrix
+  *
+  * \brief Sparse matrix
+  *
+  * \param _Scalar the scalar type, i.e. the type of the coefficients
+  *
+  * See http://www.netlib.org/linalg/html_templates/node91.html for details on the storage scheme.
+  *
+  */
+template<typename _Scalar, int _Flags>
+struct ei_traits<MappedSparseMatrix<_Scalar, _Flags> > : ei_traits<SparseMatrix<_Scalar, _Flags> >
+{};
+
+template<typename _Scalar, int _Flags>
+class MappedSparseMatrix
+  : public SparseMatrixBase<MappedSparseMatrix<_Scalar, _Flags> >
+{
+  public:
+    EIGEN_SPARSE_GENERIC_PUBLIC_INTERFACE(MappedSparseMatrix)
+
+  protected:
+    enum { IsRowMajor = Base::IsRowMajor };
+
+    int m_outerSize;
+    int m_innerSize;
+    int m_nnz;
+    int* m_outerIndex;
+    int* m_innerIndices;
+    Scalar* m_values;
+
+  public:
+
+    inline int rows() const { return IsRowMajor ? m_outerSize : m_innerSize; }
+    inline int cols() const { return IsRowMajor ? m_innerSize : m_outerSize; }
+    inline int innerSize() const { return m_innerSize; }
+    inline int outerSize() const { return m_outerSize; }
+    inline int innerNonZeros(int j) const { return m_outerIndex[j+1]-m_outerIndex[j]; }
+
+    //----------------------------------------
+    // direct access interface
+    inline const Scalar* _valuePtr() const { return &m_values; }
+    inline Scalar* _valuePtr() { return &m_values; }
+
+    inline const int* _innerIndexPtr() const { return &m_innerIndices; }
+    inline int* _innerIndexPtr() { return m_innerIndices; }
+
+    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;
+      const int inner = RowMajor ? col : row;
+
+      int start = m_outerIndex[outer];
+      int end = m_outerIndex[outer+1];
+      if (start==end)
+        return Scalar(0);
+      else if (end>0 && inner==m_innerIndices[end-1])
+        return m_values[end-1];
+      // ^^  optimization: let's first check if it is the last coefficient
+      // (very common in high level algorithms)
+
+      const int* r = std::lower_bound(&m_innerIndices[start],&m_innerIndices[end-1],inner);
+      const int id = r-&m_innerIndices[0];
+      return ((*r==inner) && (id<end)) ? m_values[id] : Scalar(0);
+    }
+
+    inline Scalar& coeffRef(int row, int col)
+    {
+      const int outer = RowMajor ? row : col;
+      const int inner = RowMajor ? col : row;
+
+      int start = m_outerIndex[outer];
+      int end = m_outerIndex[outer+1];
+      ei_assert(end>=start && "you probably called coeffRef on a non finalized matrix");
+      ei_assert(end>start && "coeffRef cannot be called on a zero coefficient");
+      int* r = std::lower_bound(&m_innerIndices[start],&m_innerIndices[end],inner);
+      const int id = r-&m_innerIndices[0];
+      ei_assert((*r==inner) && (id<end) && "coeffRef cannot be called on a zero coefficient");
+      return m_values[id];
+    }
+    
+    class InnerIterator;
+
+    /** \returns the number of non zero coefficients */
+    inline int nonZeros() const  { return m_nnz; }
+
+    inline MappedSparseMatrix(int rows, int cols, int nnz, int* outerIndexPtr, int* innerIndexPtr, Scalar* valuePtr)
+      : m_outerSize(IsRowMajor?rows:cols), m_innerSize(IsRowMajor?cols:rows), m_nnz(nnz), m_outerIndex(outerIndexPtr),
+        m_innerIndices(innerIndexPtr), m_values(valuePtr)
+    {}
+
+    #ifdef EIGEN_TAUCS_SUPPORT
+    explicit MappedSparseMatrix(taucs_ccs_matrix& taucsMatrix);
+    #endif
+
+    #ifdef EIGEN_CHOLMOD_SUPPORT
+    explicit MappedSparseMatrix(cholmod_sparse& cholmodMatrix);
+    #endif
+
+    #ifdef EIGEN_SUPERLU_SUPPORT
+    explicit MappedSparseMatrix(SluMatrix& sluMatrix);
+    #endif
+
+    /** Empty destructor */
+    inline ~MappedSparseMatrix() {}
+};
+
+template<typename Scalar, int _Flags>
+class MappedSparseMatrix<Scalar,_Flags>::InnerIterator
+{
+  public:
+    InnerIterator(const MappedSparseMatrix& mat, int outer)
+      : m_matrix(mat), m_id(mat._outerIndexPtr[outer]), m_start(m_id), m_end(mat._outerIndexPtr[outer+1])
+    {}
+
+    template<unsigned int Added, unsigned int Removed>
+    InnerIterator(const Flagged<MappedSparseMatrix,Added,Removed>& mat, int outer)
+      : m_matrix(mat._expression()), m_id(m_matrix._outerIndexPtr[outer]),
+        m_start(m_id), m_end(m_matrix._outerIndexPtr[outer+1])
+    {}
+
+    inline InnerIterator& operator++() { m_id++; return *this; }
+
+    inline Scalar value() const { return m_matrix.m_valuePtr[m_id]; }
+    inline Scalar& valueRef() { return const_cast<Scalar&>(m_matrix._valuePtr[m_id]); }
+
+    inline int index() const { return m_matrix._innerIndexPtr(m_id); }
+
+    inline operator bool() const { return (m_id < m_end) && (m_id>=m_start); }
+
+  protected:
+    const MappedSparseMatrix& m_matrix;
+    int m_id;
+    const int m_start;
+    const int m_end;
+};
+
+#endif // EIGEN_MAPPED_SPARSEMATRIX_H

Eigen/src/Sparse/SparseMatrix.h

@@ -57,20 +57,17 @@ class SparseMatrix
 {
   public:
     EIGEN_SPARSE_GENERIC_PUBLIC_INTERFACE(SparseMatrix)
+    typedef MappedSparseMatrix<Scalar,Flags> Map;
 
   protected:
-  public:
 
-    typedef SparseMatrixBase<SparseMatrix> SparseBase;
-    enum {
-      RowMajor = SparseBase::IsRowMajor
-    };
+    enum { RowMajor = Base::IsRowMajor };
     typedef SparseMatrix<Scalar,(Flags&~RowMajorBit)|(RowMajor?RowMajorBit:0)> TransposedSparseMatrix;
 
     int m_outerSize;
     int m_innerSize;
     int* m_outerIndex;
-    SparseArray<Scalar> m_data;
+    CompressedStorage<Scalar> m_data;
 
 
   public:
@@ -380,21 +377,6 @@ 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
-
-    #ifdef EIGEN_SUPERLU_SUPPORT
-    static SparseMatrix Map(SluMatrix& sluMatrix);
-    SluMatrix asSluMatrix();
-    #endif
-
     /** Destructor */
     inline ~SparseMatrix()
     {

Eigen/src/Sparse/SparseMatrixBase.h

@@ -119,7 +119,7 @@ template<typename Derived> class SparseMatrixBase
     inline int size() const { return rows() * cols(); }
     /** \returns the number of nonzero coefficients which is in practice the number
       * of stored coefficients. */
-    inline int nonZeros() const { return derived.nonZeros(); }
+    inline int nonZeros() const { return derived().nonZeros(); }
     /** \returns true if either the number of rows or the number of columns is equal to 1.
       * In other words, this function returns
       * \code rows()==1 || cols()==1 \endcode
@@ -596,6 +596,18 @@ template<typename Derived> class SparseMatrixBase
 //       return res;
 //     }
 
+    #ifdef EIGEN_TAUCS_SUPPORT
+    taucs_ccs_matrix asTaucsMatrix();
+    #endif
+
+    #ifdef EIGEN_CHOLMOD_SUPPORT
+    cholmod_sparse asCholmodMatrix();
+    #endif
+
+    #ifdef EIGEN_SUPERLU_SUPPORT
+    SluMatrix asSluMatrix();
+    #endif
+
   protected:
 
     bool m_isRValue;

Eigen/src/Sparse/SparseUtil.h

@@ -103,6 +103,7 @@ enum {
 template<typename Derived> class SparseMatrixBase;
 template<typename _Scalar, int _Flags = 0> class SparseMatrix;
 template<typename _Scalar, int _Flags = 0> class SparseVector;
+template<typename _Scalar, int _Flags = 0> class MappedSparseMatrix;
 
 template<typename MatrixType> class SparseTranspose;
 template<typename MatrixType> class SparseInnerVector;

Eigen/src/Sparse/SparseVector.h

@@ -68,7 +68,7 @@ class SparseVector
       IsColVector = ei_traits<SparseVector>::IsColVector
     };
 
-    SparseArray<Scalar> m_data;
+    CompressedStorage<Scalar> m_data;
     int m_size;
 
 

Eigen/src/Sparse/SuperLUSupport.h

@@ -104,12 +104,58 @@ struct SluMatrix : SuperMatrix
       ei_assert(false && "Scalar type not supported by SuperLU");
     }
   }
+  
+  template<typename Scalar, int Rows, int Cols, int Options, int MRows, int MCols>
+  static SluMatrix Map(Matrix<Scalar,Rows,Cols,Options,MRows,MCols>& mat)
+  {
+    typedef Matrix<Scalar,Rows,Cols,Options,MRows,MCols> MatrixType;
+    ei_assert( ((Options&RowMajor)!=RowMajor) && "row-major dense matrices is not supported by SuperLU");
+    SluMatrix res;
+    res.setStorageType(SLU_DN);
+    res.setScalarType<Scalar>();
+    res.Mtype     = SLU_GE;
+
+    res.nrow      = mat.rows();
+    res.ncol      = mat.cols();
+
+    res.storage.lda       = mat.stride();
+    res.storage.values    = mat.data();
+    return res;
+  }
 
   template<typename MatrixType>
-  static SluMatrix Map(MatrixType& mat)
+  static SluMatrix Map(SparseMatrixBase<MatrixType>& mat)
   {
     SluMatrix res;
-    SluMatrixMapHelper<MatrixType>::run(mat, res);
+    if ((MatrixType::Flags&RowMajorBit)==RowMajorBit)
+    {
+      res.setStorageType(SLU_NR);
+      res.nrow      = mat.cols();
+      res.ncol      = mat.rows();
+    }
+    else
+    {
+      res.setStorageType(SLU_NC);
+      res.nrow      = mat.rows();
+      res.ncol      = mat.cols();
+    }
+
+    res.Mtype     = SLU_GE;
+
+    res.storage.nnz       = mat.nonZeros();
+    res.storage.values    = mat.derived()._valuePtr();
+    res.storage.innerInd  = mat.derived()._innerIndexPtr();
+    res.storage.outerInd  = mat.derived()._outerIndexPtr();
+
+    res.setScalarType<typename MatrixType::Scalar>();
+
+    // FIXME the following is not very accurate
+    if (MatrixType::Flags & UpperTriangular)
+      res.Mtype = SLU_TRU;
+    if (MatrixType::Flags & LowerTriangular)
+      res.Mtype = SLU_TRL;
+    if (MatrixType::Flags & SelfAdjoint)
+      ei_assert(false && "SelfAdjoint matrix shape not supported by SuperLU");
     return res;
   }
 };
@@ -133,13 +179,13 @@ struct SluMatrixMapHelper<Matrix<Scalar,Rows,Cols,Options,MRows,MCols> >
   }
 };
 
-template<typename Scalar, int Flags>
-struct SluMatrixMapHelper<SparseMatrix<Scalar,Flags> >
+template<typename Derived>
+struct SluMatrixMapHelper<SparseMatrixBase<Derived> >
 {
-  typedef SparseMatrix<Scalar,Flags> MatrixType;
+  typedef Derived MatrixType;
   static void run(MatrixType& mat, SluMatrix& res)
   {
-    if ((Flags&RowMajorBit)==RowMajorBit)
+    if ((MatrixType::Flags&RowMajorBit)==RowMajorBit)
     {
       res.setStorageType(SLU_NR);
       res.nrow      = mat.cols();
@@ -159,48 +205,43 @@ struct SluMatrixMapHelper<SparseMatrix<Scalar,Flags> >
     res.storage.innerInd  = mat._innerIndexPtr();
     res.storage.outerInd  = mat._outerIndexPtr();
 
-    res.setScalarType<Scalar>();
+    res.setScalarType<typename MatrixType::Scalar>();
 
     // FIXME the following is not very accurate
-    if (Flags & UpperTriangular)
+    if (MatrixType::Flags & UpperTriangular)
       res.Mtype = SLU_TRU;
-    if (Flags & LowerTriangular)
+    if (MatrixType::Flags & LowerTriangular)
       res.Mtype = SLU_TRL;
-    if (Flags & SelfAdjoint)
+    if (MatrixType::Flags & SelfAdjoint)
       ei_assert(false && "SelfAdjoint matrix shape not supported by SuperLU");
   }
 };
 
-template<typename Scalar, int Flags>
-SluMatrix SparseMatrix<Scalar,Flags>::asSluMatrix()
+template<typename Derived>
+SluMatrix SparseMatrixBase<Derived>::asSluMatrix()
 {
-  return SluMatrix::Map(*this);
+  return SluMatrix::Map(derived());
 }
 
 template<typename Scalar, int Flags>
-SparseMatrix<Scalar,Flags> SparseMatrix<Scalar,Flags>::Map(SluMatrix& sluMat)
+MappedSparseMatrix<Scalar,Flags>::MappedSparseMatrix(SluMatrix& sluMat)
 {
-  SparseMatrix res;
   if ((Flags&RowMajorBit)==RowMajorBit)
   {
     assert(sluMat.Stype == SLU_NR);
-    res.m_innerSize   = sluMat.ncol;
-    res.m_outerSize   = sluMat.nrow;
+    m_innerSize   = sluMat.ncol;
+    m_outerSize   = sluMat.nrow;
   }
   else
   {
     assert(sluMat.Stype == SLU_NC);
-    res.m_innerSize   = sluMat.nrow;
-    res.m_outerSize   = sluMat.ncol;
+    m_innerSize   = sluMat.nrow;
+    m_outerSize   = sluMat.ncol;
   }
-  res.m_outerIndex  = sluMat.storage.outerInd;
-  SparseArray<Scalar> data = SparseArray<Scalar>::Map(
-                                sluMat.storage.innerInd,
-                                reinterpret_cast<Scalar*>(sluMat.storage.values),
-                                sluMat.storage.outerInd[res.m_outerSize]);
-  res.m_data.swap(data);
-  res.markAsRValue();
-  return res;
+  m_outerIndex = sluMat.storage.outerInd;
+  m_innerIndices = sluMat.storage.innerInd;
+  m_values = reinterpret_cast<Scalar*>(sluMat.storage.values);
+  m_nnz = sluMat.storage.outerInd[m_outerSize];
 }
 
 template<typename MatrixType>

Eigen/src/Sparse/TaucsSupport.h

@@ -25,8 +25,8 @@
 #ifndef EIGEN_TAUCSSUPPORT_H
 #define EIGEN_TAUCSSUPPORT_H
 
-template<typename Scalar, int Flags>
-taucs_ccs_matrix SparseMatrix<Scalar,Flags>::asTaucsMatrix()
+template<typename Derived>
+taucs_ccs_matrix SparseMatrixBase<Derived>::asTaucsMatrix()
 {
   taucs_ccs_matrix res;
   res.n         = cols();
@@ -63,19 +63,14 @@ taucs_ccs_matrix SparseMatrix<Scalar,Flags>::asTaucsMatrix()
 }
 
 template<typename Scalar, int Flags>
-SparseMatrix<Scalar,Flags> SparseMatrix<Scalar,Flags>::Map(taucs_ccs_matrix& taucsMat)
+MappedSparseMatrix<Scalar,Flags>::MappedSparseMatrix(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;
+  m_innerSize = taucsMat.m;
+  m_outerSize = taucsMat.n;
+  m_outerIndex = taucsMat.colptr;
+  m_innerIndices = taucsMat.rowind;
+  m_values = reinterpret_cast<Scalar*>(taucsMat.values.v);
+  m_nnz = taucsMat.colptr[taucsMat.n];
 }
 
 template<typename MatrixType>