* added innerSize / outerSize functions to MatrixBase

* added complete implementation of sparse matrix product
  (with a little glue in Eigen/Core)
* added an exhaustive bench of sparse products including GMM++ and MTL4
  => Eigen outperforms in all transposed/density configurations !

3 files changed, 274 insertions, 211 deletions

diff --git a/bench/BenchSparseUtil.h b/bench/BenchSparseUtil.h
new file mode 100644
index 000000000..97838f4b1
--- /dev/null
+++ b/bench/BenchSparseUtil.h
@@ -0,0 +1,75 @@
+
+#include <Eigen/Array>
+#include <Eigen/Sparse>
+#include <bench/BenchTimer.h>
+
+
+
+using namespace std;
+using namespace Eigen;
+USING_PART_OF_NAMESPACE_EIGEN
+
+#ifndef SIZE
+#define SIZE 1024
+#endif
+
+#ifndef DENSITY
+#define DENSITY 0.01
+#endif
+
+#ifndef SCALAR
+#define SCALAR float
+#endif
+
+typedef SCALAR Scalar;
+typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix;
+typedef SparseMatrix<Scalar> EigenSparseMatrix;
+
+void fillMatrix(float density, int rows, int cols,  EigenSparseMatrix& dst)
+{
+  dst.startFill(rows*cols*density);
+  for(int j = 0; j < cols; j++)
+  {
+    for(int i = 0; i < rows; i++)
+    {
+      Scalar v = (ei_random<float>(0,1) < density) ? ei_random<Scalar>() : 0;
+      if (v!=0)
+        dst.fill(i,j) = v;
+    }
+  }
+  dst.endFill();
+}
+
+void eiToDense(const EigenSparseMatrix& src, DenseMatrix& dst)
+{
+  dst.setZero();
+  for (int j=0; j<src.cols(); ++j)
+    for (EigenSparseMatrix::InnerIterator it(src.derived(), j); it; ++it)
+      dst(it.index(),j) = it.value();
+}
+
+#ifndef NOGMM
+#include "gmm/gmm.h"
+typedef gmm::csc_matrix<Scalar> GmmSparse;
+typedef gmm::col_matrix< gmm::wsvector<Scalar> > GmmDynSparse;
+void eiToGmm(const EigenSparseMatrix& src, GmmSparse& dst)
+{
+  GmmDynSparse tmp(src.rows(), src.cols());
+  for (int j=0; j<src.cols(); ++j)
+    for (EigenSparseMatrix::InnerIterator it(src.derived(), j); it; ++it)
+      tmp(it.index(),j) = it.value();
+  gmm::copy(tmp, dst);
+}
+#endif
+
+#ifndef NOMTL
+#include <boost/numeric/mtl/mtl.hpp>
+typedef mtl::compressed2D<Scalar> MtlSparse;
+void eiToMtl(const EigenSparseMatrix& src, MtlSparse& dst)
+{
+  mtl::matrix::inserter<MtlSparse> ins(dst);
+  for (int j=0; j<src.cols(); ++j)
+    for (EigenSparseMatrix::InnerIterator it(src.derived(), j); it; ++it)
+      ins[it.index()][j] = it.value();
+}
+#endif
diff --git a/bench/sparse_01.cpp b/bench/sparse_01.cpp
deleted file mode 100644
index 69ab4dda3..000000000
--- a/bench/sparse_01.cpp
+++ /dev/null
@@ -1,211 +0,0 @@
-
-// 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 10
-#endif
-
-#define REPEATPRODUCT 1
-
-#define SIZE 10
-#define DENSITY 0.2
-
-// #define NODENSEMATRIX
-
-typedef MatrixXf DenseMatrix;
-// typedef Matrix<float,SIZE,SIZE> 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, DenseMatrix* 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 = SIZE;
-  int cols = SIZE;
-  float density = DENSITY;
-
-  // dense matrices
-  #ifndef NODENSEMATRIX
-  DenseMatrix m1(rows,cols), m2(rows,cols), m3(rows,cols), m4(rows,cols);
-  #endif
-
-  // sparse matrices
-  EigenSparseMatrix sm1(rows,cols), sm2(rows,cols), sm3(rows,cols), sm4(rows,cols);
-
-  HashMatrix<float> hm4(rows,cols);
-
-  // GMM++ matrices
-
-  GmmDynSparse  gmmT4(rows,cols);
-  GmmSparse gmmM1(rows,cols), gmmM2(rows,cols), gmmM3(rows,cols), gmmM4(rows,cols);
-
-  #ifndef NODENSEMATRIX
-  fillMatrix(density, rows, cols, &m1, &sm1, &gmmM1);
-  fillMatrix(density, rows, cols, &m2, &sm2, &gmmM2);
-  fillMatrix(density, rows, cols, &m3, &sm3, &gmmM3);
-  #else
-  fillMatrix(density, rows, cols, 0, &sm1, &gmmM1);
-  fillMatrix(density, rows, cols, 0, &sm2, &gmmM2);
-  fillMatrix(density, rows, cols, 0, &sm3, &gmmM3);
-  #endif
-
-  BenchTimer timer;
-
-  //--------------------------------------------------------------------------------
-  //  COEFF WISE OPERATORS
-  //--------------------------------------------------------------------------------
-#if 1
-  std::cout << "\n\n\"m4 = m1 + m2 + 2 * m3\":\n\n";
-
-  timer.reset();
-  timer.start();
-  asm("#begin");
-  for (int k=0; k<REPEAT; ++k)
-    m4 = m1 + m2 + 2 * m3;
-  asm("#end");
-  timer.stop();
-  std::cout << "Eigen dense = " << timer.value() << endl;
-
-  timer.reset();
-  timer.start();
-  for (int k=0; k<REPEAT; ++k)
-    sm4 = sm1 + sm2 + 2 * sm3;
-  timer.stop();
-  std::cout << "Eigen sparse = " << timer.value() << endl;
-
-  timer.reset();
-  timer.start();
-  for (int k=0; k<REPEAT; ++k)
-    hm4 = sm1 + sm2 + 2 * sm3;
-  timer.stop();
-  std::cout << "Eigen hash = " << timer.value() << endl;
-
-  LinkedVectorMatrix<float> lm4(rows, cols);
-  timer.reset();
-  timer.start();
-  for (int k=0; k<REPEAT; ++k)
-    lm4 = sm1 + sm2 + 2 * sm3;
-  timer.stop();
-  std::cout << "Eigen linked vector = " << timer.value() << endl;
-
-  timer.reset();
-  timer.start();
-  for (int k=0; k<REPEAT; ++k)
-  {
-    gmm::add(gmmM1, gmmM2, gmmT4);
-    gmm::add(gmm::scaled(gmmM3,2), gmmT4);
-  }
-  timer.stop();
-  std::cout << "GMM++ sparse = " << timer.value() << endl;
-#endif
-  //--------------------------------------------------------------------------------
-  //  PRODUCT
-  //--------------------------------------------------------------------------------
-#if 0
-  std::cout << "\n\nProduct:\n\n";
-
-  #ifndef NODENSEMATRIX
-  timer.reset();
-  timer.start();
-  asm("#begin");
-  for (int k=0; k<REPEATPRODUCT; ++k)
-    m1 = m1 * m2;
-  asm("#end");
-  timer.stop();
-  std::cout << "Eigen dense = " << timer.value() << endl;
-  #endif
-
-  timer.reset();
-  timer.start();
-  for (int k=0; k<REPEATPRODUCT; ++k)
-    sm4 = sm1 * sm2;
-  timer.stop();
-  std::cout << "Eigen sparse = " << timer.value() << endl;
-
-//   timer.reset();
-//   timer.start();
-//   for (int k=0; k<REPEATPRODUCT; ++k)
-//     hm4 = sm1 * sm2;
-//   timer.stop();
-//   std::cout << "Eigen hash = " << timer.value() << endl;
-
-  timer.reset();
-  timer.start();
-  for (int k=0; k<REPEATPRODUCT; ++k)
-  {
-    gmm::csr_matrix<float> R(rows,cols);
-    gmm::copy(gmmM1, R);
-    //gmm::mult(gmmM1, gmmM2, gmmT4);
-  }
-  timer.stop();
-  std::cout << "GMM++ sparse = " << timer.value() << endl;
-#endif
-  //--------------------------------------------------------------------------------
-  //  VARIOUS
-  //--------------------------------------------------------------------------------
-#if 1
-//   sm3 = sm1 + m2;
-  cout << m4.transpose() << "\n\n";
-//   sm4 = sm1+sm2;
-  cout << sm4 << "\n\n";
-  cout << lm4 << endl;
-
-  LinkedVectorMatrix<float,RowMajorBit> lm5(rows, cols);
-  lm5 = lm4;
-  lm5 = sm4;
-  cout << endl << lm5 << endl;
-
-  sm3 = sm4.transpose();
-  cout << endl << lm5 << endl;
-
-  cout << endl << "SM1 before random editing: " << endl << sm1 << endl;
-  {
-    SparseSetter<EigenSparseMatrix,RandomAccessPattern> w1(sm1);
-    w1->coeffRef(4,2) = ei_random<float>();
-    w1->coeffRef(2,6) = ei_random<float>();
-    w1->coeffRef(0,4) = ei_random<float>();
-    w1->coeffRef(9,3) = ei_random<float>();
-  }
-  cout << endl << "SM1 after random editing: " << endl << sm1 << endl;
-#endif
-
-  return 0;
-}
-
diff --git a/bench/sparse_product.cpp b/bench/sparse_product.cpp
new file mode 100644
index 000000000..846301fa5
--- /dev/null
+++ b/bench/sparse_product.cpp
@@ -0,0 +1,199 @@
+
+//g++ -O3 -g0 -DNDEBUG  sparse_product.cpp -I.. -I/home/gael/Coding/LinearAlgebra/mtl4/ -DDENSITY=0.005 -DSIZE=10000 && ./a.out
+//g++ -O3 -g0 -DNDEBUG  sparse_product.cpp -I.. -I/home/gael/Coding/LinearAlgebra/mtl4/ -DDENSITY=0.05 -DSIZE=2000 && ./a.out
+// -DNOGMM -DNOMTL
+
+#ifndef SIZE
+#define SIZE 10000
+#endif
+
+#ifndef DENSITY
+#define DENSITY 0.01
+#endif
+
+#ifndef REPEAT
+#define REPEAT 1
+#endif
+
+#include "BenchSparseUtil.h"
+
+#ifndef MINDENSITY
+#define MINDENSITY 0.0004
+#endif
+
+int main(int argc, char *argv[])
+{
+  int rows = SIZE;
+  int cols = SIZE;
+  float density = DENSITY;
+
+  EigenSparseMatrix sm1(rows,cols), sm2(rows,cols), sm3(rows,cols), sm4(rows,cols);
+
+  BenchTimer timer;
+  for (float density = DENSITY; density>=MINDENSITY; density*=0.5)
+  {
+    fillMatrix(density, rows, cols, sm1);
+    fillMatrix(density, rows, cols, sm2);
+
+    // dense matrices
+    #ifdef DENSEMATRIX
+    {
+      std::cout << "Eigen Dense\t" << density*100 << "%\n";
+      DenseMatrix m1(rows,cols), m2(rows,cols), m3(rows,cols);
+      eiToDense(sm1, m1);
+      eiToDense(sm2, m2);
+
+      timer.reset();
+      timer.start();
+      for (int k=0; k<REPEAT; ++k)
+        m3 = m1 * m2;
+      timer.stop();
+      std::cout << "   a * b:\t" << timer.value() << endl;
+
+      timer.reset();
+      timer.start();
+      for (int k=0; k<REPEAT; ++k)
+        m3 = m1.transpose() * m2;
+      timer.stop();
+      std::cout << "   a' * b:\t" << timer.value() << endl;
+
+      timer.reset();
+      timer.start();
+      for (int k=0; k<REPEAT; ++k)
+        m3 = m1.transpose() * m2.transpose();
+      timer.stop();
+      std::cout << "   a' * b':\t" << timer.value() << endl;
+
+      timer.reset();
+      timer.start();
+      for (int k=0; k<REPEAT; ++k)
+        m3 = m1 * m2.transpose();
+      timer.stop();
+      std::cout << "   a * b':\t" << timer.value() << endl;
+    }
+    #endif
+
+    // eigen sparse matrices
+    {
+      std::cout << "Eigen sparse\t" << density*100 << "%\n";
+
+      timer.reset();
+      timer.start();
+      for (int k=0; k<REPEAT; ++k)
+        sm3 = sm1 * sm2;
+      timer.stop();
+      std::cout << "   a * b:\t" << timer.value() << endl;
+
+      timer.reset();
+      timer.start();
+      for (int k=0; k<REPEAT; ++k)
+        sm3 = sm1.transpose() * sm2;
+      timer.stop();
+      std::cout << "   a' * b:\t" << timer.value() << endl;
+
+      timer.reset();
+      timer.start();
+      for (int k=0; k<REPEAT; ++k)
+        sm3 = sm1.transpose() * sm2.transpose();
+      timer.stop();
+      std::cout << "   a' * b':\t" << timer.value() << endl;
+
+      timer.reset();
+      timer.start();
+      for (int k=0; k<REPEAT; ++k)
+        sm3 = sm1 * sm2.transpose();
+      timer.stop();
+      std::cout << "   a * b' :\t" << timer.value() << endl;
+    }
+
+    // GMM++
+    #ifndef NOGMM
+    {
+      std::cout << "GMM++ sparse\t" << density*100 << "%\n";
+      GmmDynSparse  gmmT3(rows,cols);
+      GmmSparse m1(rows,cols), m2(rows,cols), m3(rows,cols);
+      eiToGmm(sm1, m1);
+      eiToGmm(sm2, m2);
+
+      timer.reset();
+      timer.start();
+      for (int k=0; k<REPEAT; ++k)
+        gmm::mult(m1, m2, gmmT3);
+      timer.stop();
+      std::cout << "   a * b:\t" << timer.value() << endl;
+
+      timer.reset();
+      timer.start();
+      for (int k=0; k<REPEAT; ++k)
+        gmm::mult(gmm::transposed(m1), m2, gmmT3);
+      timer.stop();
+      std::cout << "   a' * b:\t" << timer.value() << endl;
+
+      if (rows<500)
+      {
+        timer.reset();
+        timer.start();
+        for (int k=0; k<REPEAT; ++k)
+          gmm::mult(gmm::transposed(m1), gmm::transposed(m2), gmmT3);
+        timer.stop();
+        std::cout << "   a' * b':\t" << timer.value() << endl;
+
+        timer.reset();
+        timer.start();
+        for (int k=0; k<REPEAT; ++k)
+          gmm::mult(m1, gmm::transposed(m2), gmmT3);
+        timer.stop();
+        std::cout << "   a * b':\t" << timer.value() << endl;
+      }
+      else
+      {
+        std::cout << "   a' * b':\t" << "forever" << endl;
+        std::cout << "   a * b':\t" << "forever" << endl;
+      }
+    }
+    #endif
+
+    // MTL4
+    #ifndef NOMTL
+    {
+      std::cout << "MTL4\t" << density*100 << "%\n";
+      MtlSparse m1(rows,cols), m2(rows,cols), m3(rows,cols);
+      eiToMtl(sm1, m1);
+      eiToMtl(sm2, m2);
+
+      timer.reset();
+      timer.start();
+      for (int k=0; k<REPEAT; ++k)
+        m3 = m1 * m2;
+      timer.stop();
+      std::cout << "   a * b:\t" << timer.value() << endl;
+
+      timer.reset();
+      timer.start();
+      for (int k=0; k<REPEAT; ++k)
+        m3 = trans(m1) * m2;
+      timer.stop();
+      std::cout << "   a' * b:\t" << timer.value() << endl;
+
+      timer.reset();
+      timer.start();
+      for (int k=0; k<REPEAT; ++k)
+        m3 = trans(m1) * trans(m2);
+      timer.stop();
+      std::cout << "  a' * b':\t" << timer.value() << endl;
+
+      timer.reset();
+      timer.start();
+      for (int k=0; k<REPEAT; ++k)
+        m3 = m1 * trans(m2);
+      timer.stop();
+      std::cout << "   a * b' :\t" << timer.value() << endl;
+    }
+    #endif
+
+    std::cout << "\n\n";
+  }
+
+  return 0;
+}
+