add experimental code for sparse matrix:

- uses the common "Compressed Column Storage" scheme
 - supports every unary and binary operators with xpr template
   assuming binaryOp(0,0) == 0 and unaryOp(0) = 0 (otherwise a sparse
   matrix doesnot make sense)
 - this is the first commit, so of course, there are still several shorcommings !

Eigen/Sparse

@@ -0,0 +1,16 @@
+#ifndef EIGEN_SPARSE_MODULE_H
+#define EIGEN_SPARSE_MODULE_H
+
+#include "Core"
+#include <vector>
+#include <map>
+
+namespace Eigen {
+
+#include "src/Sparse/SparseArray.h"
+#include "src/Sparse/SparseMatrix.h"
+#include "src/Sparse/CoreIterators.h"
+
+} // namespace Eigen
+
+#endif // EIGEN_SPARSE_MODULE_H

Eigen/src/Core/CwiseBinaryOp.h

@@ -84,6 +84,8 @@ class CwiseBinaryOp : ei_no_assignment_operator,
     typedef typename ei_traits<CwiseBinaryOp>::LhsNested LhsNested;
     typedef typename ei_traits<CwiseBinaryOp>::RhsNested RhsNested;
 
+    class InnerIterator;
+
     inline CwiseBinaryOp(const Lhs& lhs, const Rhs& rhs, const BinaryOp& func = BinaryOp())
       : m_lhs(lhs), m_rhs(rhs), m_functor(func)
     {

Eigen/src/Core/CwiseUnaryOp.h

@@ -69,6 +69,8 @@ class CwiseUnaryOp : ei_no_assignment_operator,
 
     EIGEN_GENERIC_PUBLIC_INTERFACE(CwiseUnaryOp)
 
+    class InnerIterator;
+
     inline CwiseUnaryOp(const MatrixType& mat, const UnaryOp& func = UnaryOp())
       : m_matrix(mat), m_functor(func) {}
 

Eigen/src/Core/MatrixBase.h

@@ -55,6 +55,8 @@ template<typename Derived> class MatrixBase
 
   public:
 
+    class InnerIterator;
+
     typedef typename ei_traits<Derived>::Scalar Scalar;
     typedef typename ei_packet_traits<Scalar>::type PacketScalar;
 

Eigen/src/Sparse/CoreIterators.h

@@ -0,0 +1,152 @@
+// 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_COREITERATORS_H
+#define EIGEN_COREITERATORS_H
+
+template<typename Derived>
+class MatrixBase<Derived>::InnerIterator
+{
+    typedef typename Derived::Scalar Scalar;
+  public:
+    InnerIterator(const Derived& mat, int col)
+      : m_matrix(mat), m_row(0), m_col(col), m_end(mat.rows())
+    {}
+
+    Scalar value() { return m_matrix.coeff(m_row, m_col); }
+
+    InnerIterator& operator++() { m_row++; return *this; }
+
+    int index() const { return m_row; }
+
+    operator bool() const { return m_row < m_end && m_row>=0; }
+
+  protected:
+    const Derived& m_matrix;
+    int m_row;
+    const int m_col;
+    const int m_end;
+};
+
+template<typename UnaryOp, typename MatrixType>
+class CwiseUnaryOp<UnaryOp,MatrixType>::InnerIterator
+{
+    typedef typename CwiseUnaryOp::Scalar Scalar;
+    typedef typename ei_traits<CwiseUnaryOp>::_MatrixTypeNested _MatrixTypeNested;
+    typedef typename _MatrixTypeNested::InnerIterator MatrixTypeIterator;
+  public:
+
+    InnerIterator(const CwiseUnaryOp& unaryOp, int col)
+      : m_iter(unaryOp.m_matrix,col), m_functor(unaryOp.m_functor), m_id(-1)
+    {
+      this->operator++();
+    }
+
+    InnerIterator& operator++()
+    {
+      if (m_iter)
+      {
+        m_id = m_iter.index();
+        m_value = m_functor(m_iter.value());
+        ++m_iter;
+      }
+      else
+      {
+        m_id = -1;
+      }
+      return *this;
+    }
+
+    Scalar value() const { return m_value; }
+
+    int index() const { return m_id; }
+
+    operator bool() const { return m_id>=0; }
+
+  protected:
+    MatrixTypeIterator m_iter;
+    const UnaryOp& m_functor;
+    Scalar m_value;
+    int m_id;
+};
+
+template<typename BinaryOp, typename Lhs, typename Rhs>
+class CwiseBinaryOp<BinaryOp,Lhs,Rhs>::InnerIterator
+{
+    typedef typename CwiseBinaryOp::Scalar Scalar;
+    typedef typename ei_traits<CwiseBinaryOp>::_LhsNested _LhsNested;
+    typedef typename _LhsNested::InnerIterator LhsIterator;
+    typedef typename ei_traits<CwiseBinaryOp>::_RhsNested _RhsNested;
+    typedef typename _RhsNested::InnerIterator RhsIterator;
+  public:
+
+    InnerIterator(const CwiseBinaryOp& binOp, int col)
+      : m_lhsIter(binOp.m_lhs,col), m_rhsIter(binOp.m_rhs,col), m_functor(binOp.m_functor), m_id(-1)
+    {
+      this->operator++();
+    }
+
+    InnerIterator& operator++()
+    {
+      if (m_lhsIter && m_rhsIter && (m_lhsIter.index() == m_rhsIter.index()))
+      {
+        m_id = m_lhsIter.index();
+        m_value = m_functor(m_lhsIter.value(), m_rhsIter.value());
+        ++m_lhsIter;
+        ++m_rhsIter;
+      }
+      else if (m_lhsIter && ((!m_rhsIter) || m_lhsIter.index() < m_rhsIter.index()))
+      {
+        m_id = m_lhsIter.index();
+        m_value = m_functor(m_lhsIter.value(), Scalar(0));
+        ++m_lhsIter;
+      }
+      else if (m_rhsIter && ((!m_lhsIter) || m_lhsIter.index() > m_rhsIter.index()))
+      {
+        m_id = m_rhsIter.index();
+        m_value = m_functor(Scalar(0), m_rhsIter.value());
+        ++m_rhsIter;
+      }
+      else
+      {
+        m_id = -1;
+      }
+      return *this;
+    }
+
+    Scalar value() const { return m_value; }
+
+    int index() const { return m_id; }
+
+    operator bool() const { return m_id>=0; }
+
+  protected:
+    LhsIterator m_lhsIter;
+    RhsIterator m_rhsIter;
+    const BinaryOp& m_functor;
+    Scalar m_value;
+    int m_id;
+};
+
+#endif // EIGEN_COREITERATORS_H

Eigen/src/Sparse/SparseArray.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_SPARSE_ARRAY_H
+#define EIGEN_SPARSE_ARRAY_H
+
+/** Stores a sparse set of values as a list of values and a list of indices.
+  *
+  */
+template<typename Scalar> class SparseArray
+{
+  public:
+    SparseArray()
+      : m_values(0), m_indices(0), m_size(0), m_allocatedSize(0)
+    {}
+
+    SparseArray(int size)
+      : m_values(0), m_indices(0), m_size(0), m_allocatedSize(0)
+    {
+      resize(size);
+    }
+
+    SparseArray(const SparseArray& other)
+    {
+      *this = other;
+    }
+
+    SparseArray& operator=(const SparseArray& other)
+    {
+      resize(other.size());
+      memcpy(m_values, other.m_values, m_size * sizeof(Scalar));
+      memcpy(m_indices, other.m_indices, m_size * sizeof(int));
+    }
+
+    void reserve(int size)
+    {
+      int newAllocatedSize = m_size + size;
+      if (newAllocatedSize > m_allocatedSize)
+      {
+        Scalar* newValues  = new Scalar[newAllocatedSize];
+        int* newIndices = new int[newAllocatedSize];
+        // copy
+        memcpy(newValues,  m_values,  m_size * sizeof(Scalar));
+        memcpy(newIndices, m_indices, m_size * sizeof(int));
+        // delete old stuff
+        delete[] m_values;
+        delete[] m_indices;
+        m_values = newValues;
+        m_indices = newIndices;
+        m_allocatedSize = newAllocatedSize;
+      }
+    }
+
+    void resize(int size, int reserveSizeFactor = 0)
+    {
+      if (m_allocatedSize<size)
+      {
+        int newAllocatedSize = size + reserveSizeFactor*size;
+        Scalar* newValues  = new Scalar[newAllocatedSize];
+        int* newIndices = new int[newAllocatedSize];
+        // copy
+        memcpy(newValues,  m_values,  m_size * sizeof(Scalar));
+        memcpy(newIndices, m_indices, m_size * sizeof(int));
+        // delete old stuff
+        delete[] m_values;
+        delete[] m_indices;
+        m_values = newValues;
+        m_indices = newIndices;
+        m_allocatedSize = newAllocatedSize;
+      }
+      m_size = size;
+    }
+
+    void append(const Scalar& v, int i)
+    {
+      int id = m_size;
+      resize(m_size+1, 1);
+      m_values[id] = v;
+      m_indices[id] = i;
+    }
+
+    int size() const { return m_size; }
+
+    void clear()
+    {
+      m_size = 0;
+    }
+
+    Scalar& value(int i) { return m_values[i]; }
+    Scalar value(int i) const { return m_values[i]; }
+
+    int& index(int i) { return m_indices[i]; }
+    int index(int i) const { return m_indices[i]; }
+
+  protected:
+    Scalar* m_values;
+    int* m_indices;
+    int m_size;
+    int m_allocatedSize;
+};
+
+#endif // EIGEN_SPARSE_ARRAY_H

Eigen/src/Sparse/SparseMatrix.h

@@ -0,0 +1,331 @@
+// 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_SPARSEMATRIX_H
+#define EIGEN_SPARSEMATRIX_H
+
+template<typename _Scalar> class SparseMatrix;
+
+/** \class SparseMatrix
+  *
+  * \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>
+struct ei_traits<SparseMatrix<_Scalar> >
+{
+  typedef _Scalar Scalar;
+  enum {
+    RowsAtCompileTime = Dynamic,
+    ColsAtCompileTime = Dynamic,
+    MaxRowsAtCompileTime = Dynamic,
+    MaxColsAtCompileTime = Dynamic,
+    Flags = 0,
+    CoeffReadCost = NumTraits<Scalar>::ReadCost
+  };
+};
+
+template<typename _Scalar>
+class SparseMatrix : public MatrixBase<SparseMatrix<_Scalar> >
+{
+  public:
+    EIGEN_GENERIC_PUBLIC_INTERFACE(SparseMatrix)
+
+  protected:
+
+    int* m_colPtrs;
+    SparseArray<Scalar> m_data;
+    int m_rows;
+    int m_cols;
+
+    inline int _rows() const { return m_rows; }
+    inline int _cols() const { return m_cols; }
+
+    inline const Scalar& _coeff(int row, int col) const
+    {
+      int id = m_colPtrs[col];
+      int end = m_colPtrs[col+1];
+      while (id<end && m_data.index(id)!=row)
+      {
+        ++id;
+      }
+      if (id==end)
+        return 0;
+      return m_data.value(id);
+    }
+
+    inline Scalar& _coeffRef(int row, int col)
+    {
+      int id = m_colPtrs[cols];
+      int end = m_colPtrs[cols+1];
+      while (id<end && m_data.index(id)!=row)
+      {
+        ++id;
+      }
+      ei_assert(id!=end);
+      return m_data.value(id);
+    }
+
+  public:
+
+    class InnerIterator;
+
+    inline int rows() const { return _rows(); }
+    inline int cols() const { return _cols(); }
+    /** \returns the number of non zero coefficients */
+    inline int nonZeros() const  { return m_data.size(); }
+
+    inline const Scalar& operator() (int row, int col) const
+    {
+      return _coeff(row, col);
+    }
+
+    inline Scalar& operator() (int row, int col)
+    {
+      return _coeffRef(row, col);
+    }
+
+    inline void startFill(int reserveSize = 1000)
+    {
+      m_data.clear();
+      m_data.reserve(reserveSize);
+      for (int i=0; i<=m_cols; ++i)
+        m_colPtrs[i] = 0;
+    }
+
+    inline Scalar& fill(int row, int col)
+    {
+      if (m_colPtrs[col+1]==0)
+      {
+        int i=col;
+        while (i>=0 && m_colPtrs[i]==0)
+        {
+          m_colPtrs[i] = m_data.size();
+          --i;
+        }
+        m_colPtrs[col+1] = m_colPtrs[col];
+      }
+      assert(m_colPtrs[col+1] == m_data.size());
+      int id = m_colPtrs[col+1];
+      m_colPtrs[col+1]++;
+
+      m_data.append(0, row);
+      return m_data.value(id);
+    }
+
+    inline void endFill()
+    {
+      int size = m_data.size();
+      int i = m_cols;
+      // find the last filled column
+      while (i>=0 && m_colPtrs[i]==0)
+        --i;
+      i++;
+      while (i<=m_cols)
+      {
+        m_colPtrs[i] = size;
+        ++i;
+      }
+    }
+
+    void resize(int rows, int cols)
+    {
+      if (m_cols != cols)
+      {
+        delete[] m_colPtrs;
+        m_colPtrs = new int [cols+1];
+        m_rows = rows;
+        m_cols = cols;
+      }
+    }
+
+    inline SparseMatrix(int rows, int cols)
+      : m_rows(0), m_cols(0), m_colPtrs(0)
+    {
+      resize(rows, cols);
+    }
+
+    inline SparseMatrix& operator=(const SparseMatrix& other)
+    {
+      resize(other.rows(), other.cols());
+      m_colPtrs = other.m_colPtrs;
+      for (int col=0; col<=cols(); ++col)
+        m_colPtrs[col] = other.m_colPtrs[col];
+      m_data = other.m_data;
+      return *this;
+    }
+
+    template<typename OtherDerived>
+    inline SparseMatrix& operator=(const MatrixBase<OtherDerived>& other)
+    {
+      resize(other.rows(), other.cols());
+      startFill(std::max(m_rows,m_cols)*2);
+      for (int col=0; col<cols(); ++col)
+      {
+        for (typename OtherDerived::InnerIterator it(other.derived(), col); it; ++it)
+        {
+          Scalar v = it.value();
+          if (v!=Scalar(0))
+            fill(it.index(),col) = v;
+        }
+      }
+      endFill();
+      return *this;
+    }
+
+
+    // old explicit operator+
+//     template<typename Other>
+//     SparseMatrix operator+(const Other& other)
+//     {
+//       SparseMatrix res(rows(), cols());
+//       res.startFill(nonZeros()*3);
+//       for (int col=0; col<cols(); ++col)
+//       {
+//         InnerIterator row0(*this,col);
+//         typename Other::InnerIterator row1(other,col);
+//         while (row0 && row1)
+//         {
+//           if (row0.index()==row1.index())
+//           {
+//             std::cout << "both " << col << " " << row0.index() << "\n";
+//             Scalar v = row0.value() + row1.value();
+//             if (v!=Scalar(0))
+//               res.fill(row0.index(),col) = v;
+//             ++row0;
+//             ++row1;
+//           }
+//           else if (row0.index()<row1.index())
+//           {
+//             std::cout << "row0 " << col << " " << row0.index() << "\n";
+//             Scalar v = row0.value();
+//             if (v!=Scalar(0))
+//               res.fill(row0.index(),col) = v;
+//             ++row0;
+//           }
+//           else if (row1)
+//           {
+//             std::cout << "row1 " << col << " " << row0.index() << "\n";
+//             Scalar v = row1.value();
+//             if (v!=Scalar(0))
+//               res.fill(row1.index(),col) = v;
+//             ++row1;
+//           }
+//         }
+//         while (row0)
+//         {
+//           std::cout << "row0 " << col << " " << row0.index() << "\n";
+//           Scalar v = row0.value();
+//           if (v!=Scalar(0))
+//             res.fill(row0.index(),col) = v;
+//           ++row0;
+//         }
+//         while (row1)
+//         {
+//           std::cout << "row1 " << col << " " << row1.index() << "\n";
+//           Scalar v = row1.value();
+//           if (v!=Scalar(0))
+//             res.fill(row1.index(),col) = v;
+//           ++row1;
+//         }
+//       }
+//       res.endFill();
+//       return res;
+// //       return binaryOp(other, ei_scalar_sum_op<Scalar>());
+//     }
+
+
+    // WARNING for efficiency reason it currently outputs the transposed matrix
+    friend std::ostream & operator << (std::ostream & s, const SparseMatrix& m)
+    {
+      s << "Nonzero entries:\n";
+      for (uint i=0; i<m.nonZeros(); ++i)
+      {
+        s << "(" << m.m_data.value(i) << "," << m.m_data.index(i) << ") ";
+      }
+      s << std::endl;
+      s << std::endl;
+      s << "Column pointers:\n";
+      for (uint i=0; i<m.cols(); ++i)
+      {
+        s << m.m_colPtrs[i] << " ";
+      }
+      s << std::endl;
+      s << std::endl;
+      s << "Matrix (transposed):\n";
+      for (int j=0; j<m.cols(); j++ )
+      {
+        int end = m.m_colPtrs[j+1];
+        int i=0;
+        for (int id=m.m_colPtrs[j]; id<end; id++)
+        {
+          int row = m.m_data.index(id);
+          // fill with zeros
+          for (int k=i; k<row; ++k)
+            s << "0 ";
+          i = row+1;
+          s << m.m_data.value(id) << " ";
+        }
+        for (int k=i; k<m.rows(); ++k)
+          s << "0 ";
+        s << std::endl;
+      }
+      return s;
+    }
+
+    /** Destructor */
+    inline ~SparseMatrix()
+    {
+      delete[] m_colPtrs;
+    }
+};
+
+template<typename Scalar>
+class SparseMatrix<Scalar>::InnerIterator
+{
+  public:
+    InnerIterator(const SparseMatrix& mat, int col)
+      : m_matrix(mat), m_id(mat.m_colPtrs[col]), m_start(m_id), m_end(mat.m_colPtrs[col+1])
+    {}
+
+    InnerIterator& operator++() { m_id++; return *this; }
+
+    Scalar value() { return m_matrix.m_data.value(m_id); }
+
+    int index() const { return m_matrix.m_data.index(m_id); }
+
+    operator bool() const { return (m_id < m_end) && (m_id>=m_start); }
+
+  protected:
+    const SparseMatrix& m_matrix;
+    int m_id;
+    const int m_start;
+    const int m_end;
+};
+
+#endif // EIGEN_SPARSEMATRIX_H

bench/sparse_01.cpp

@@ -0,0 +1,100 @@
+
+// g++ -O3 -DNDEBUG sparse_01.cpp -I .. -o sparse_01 && ./sparse_01
+
+#include <Eigen/Array>
+#include <Eigen/Sparse>
+#include <bench/BenchTimer.h>
+
+#include "gmm/gmm.h"
+
+using namespace std;
+using namespace Eigen;
+USING_PART_OF_NAMESPACE_EIGEN
+
+#ifndef REPEAT
+#define REPEAT 40000000
+#endif
+
+typedef MatrixXf DenseMatrix;
+typedef SparseMatrix<float> EigenSparseMatrix;
+typedef gmm::csc_matrix<float> GmmSparse;
+typedef gmm::col_matrix< gmm::wsvector<float> > GmmDynSparse;
+
+void fillMatrix(float density, int rows, int cols, MatrixXf* pDenseMatrix, EigenSparseMatrix* pSparseMatrix, GmmSparse* pGmmMatrix=0)
+{
+  GmmDynSparse gmmT(rows, cols);
+  if (pSparseMatrix)
+    pSparseMatrix->startFill(rows*cols*density);
+  for(int j = 0; j < cols; j++)
+  {
+    for(int i = 0; i < rows; i++)
+    {
+      float v = (ei_random<float>(0,1) < density) ? ei_random<float>() : 0;
+      if (pDenseMatrix)
+        (*pDenseMatrix)(i,j) = v;
+      if (v!=0)
+      {
+        if (pSparseMatrix)
+          pSparseMatrix->fill(i,j) = v;
+        if (pGmmMatrix)
+          gmmT(i,j) = v;
+      }
+    }
+  }
+  if (pSparseMatrix)
+    pSparseMatrix->endFill();
+  if (pGmmMatrix)
+    gmm::copy(gmmT, *pGmmMatrix);
+}
+
+int main(int argc, char *argv[])
+{
+  int rows = 4000;
+  int cols = 4000;
+  float density = 0.1;
+
+  // dense matrices
+  DenseMatrix m1(rows,cols), m2(rows,cols), m3(rows,cols), m4(rows,cols);
+
+  // sparse matrices
+  EigenSparseMatrix sm1(rows,cols), sm2(rows,cols), sm3(rows,cols), sm4(rows,cols);
+
+  // GMM++ matrices
+
+  GmmDynSparse  gmmT4(rows,cols);
+  GmmSparse gmmM1(rows,cols), gmmM2(rows,cols), gmmM3(rows,cols), gmmM4(rows,cols);
+
+  fillMatrix(density, rows, cols, &m1, &sm1, &gmmM1);
+  fillMatrix(density, rows, cols, &m2, &sm2, &gmmM2);
+  fillMatrix(density, rows, cols, &m3, &sm3, &gmmM3);
+
+  BenchTimer timer;
+
+  timer.start();
+  for (int k=0; k<10; ++k)
+    m4 = m1 + m2 + 2 * m3;
+  timer.stop();
+  std::cout << "Eigen dense = " << timer.value() << endl;
+
+  timer.reset();
+  timer.start();
+  for (int k=0; k<10; ++k)
+    sm4 = sm1 + sm2 + 2 * sm3;
+  timer.stop();
+  std::cout << "Eigen sparse = " << timer.value() << endl;
+
+  timer.reset();
+  timer.start();
+  for (int k=0; k<10; ++k)
+  {
+    gmm::add(gmmM1, gmmM2, gmmT4);
+    gmm::add(gmm::scaled(gmmM3,2), gmmT4);
+  }
+  timer.stop();
+  std::cout << "GMM++ sparse = " << timer.value() << endl;
+
+//   sm3 = sm1 + m2;
+//   cout << m4.transpose() << "\n\n" << sm4 << endl;
+  return 0;
+}
+