From 905050b2392681ffb9b0f98b9fd6e738a562ae76 Mon Sep 17 00:00:00 2001
From: Gael Guennebaud <g.gael@free.fr>
Date: Tue, 9 Feb 2010 15:55:36 +0100
Subject: extend sparse product benchmark with ublas

---
 bench/BenchSparseUtil.h  |  14 +++++
 bench/sparse_product.cpp | 147 +++++++++++++++++++++++++++--------------------
 2 files changed, 98 insertions(+), 63 deletions(-)

(limited to 'bench')

diff --git a/bench/BenchSparseUtil.h b/bench/BenchSparseUtil.h
index 39db69345..a5ab10711 100644
--- a/bench/BenchSparseUtil.h
+++ b/bench/BenchSparseUtil.h
@@ -123,4 +123,18 @@ void eiToCSparse(const EigenSparseMatrix& src, cs* &dst)
 #include <boost/numeric/ublas/matrix_sparse.hpp>
 #include <boost/numeric/ublas/vector_of_vector.hpp>
 
+// using namespace boost;
+// using namespace boost::numeric;
+// using namespace boost::numeric::ublas;
+
+typedef boost::numeric::ublas::compressed_matrix<Scalar,boost::numeric::ublas::column_major> UblasMatrix;
+
+void eiToUblas(const EigenSparseMatrix& src, UblasMatrix& dst)
+{
+  for (int j=0; j<src.cols(); ++j)
+    for (EigenSparseMatrix::InnerIterator it(src.derived(), j); it; ++it)
+      dst(it.index(),j) = it.value();
+}
+
+
 #endif
diff --git a/bench/sparse_product.cpp b/bench/sparse_product.cpp
index 0c72993cd..0b5558b89 100644
--- a/bench/sparse_product.cpp
+++ b/bench/sparse_product.cpp
@@ -3,12 +3,15 @@
 //g++ -O3 -g0 -DNDEBUG  sparse_product.cpp -I.. -I/home/gael/Coding/LinearAlgebra/mtl4/ -DDENSITY=0.05 -DSIZE=2000 && ./a.out
 // -DNOGMM -DNOMTL -DCSPARSE
 // -I /home/gael/Coding/LinearAlgebra/CSparse/Include/ /home/gael/Coding/LinearAlgebra/CSparse/Lib/libcsparse.a
+
+#include <typeinfo>
+
 #ifndef SIZE
-#define SIZE 10000
+#define SIZE 1000000
 #endif
 
 #ifndef NNZPERCOL
-#define NNZPERCOL 32
+#define NNZPERCOL 6
 #endif
 
 #ifndef REPEAT
@@ -20,7 +23,7 @@
 #include "BenchSparseUtil.h"
 
 #ifndef NBTRIES
-#define NBTRIES 4
+#define NBTRIES 1
 #endif
 
 #define BENCH(X) \
@@ -60,6 +63,7 @@
 cs* cs_sorted_multiply(const cs* a, const cs* b)
 {
 //   return cs_multiply(a,b);
+
   cs* A = cs_transpose(a, 1);
   cs* B = cs_transpose(b, 1);
   cs* D = cs_multiply(B,A);   /* D = B'*A' */
@@ -155,26 +159,26 @@ int main(int argc, char *argv[])
 
 //       BENCH(sm3 = sm1.transpose() * sm2; )
 //       std::cout << "   a' * b:\t" << timer.value() << endl;
-//
+// //
 //       BENCH(sm3 = sm1.transpose() * sm2.transpose(); )
 //       std::cout << "   a' * b':\t" << timer.value() << endl;
-//
+// //
 //       BENCH(sm3 = sm1 * sm2.transpose(); )
 //       std::cout << "   a * b' :\t" << timer.value() << endl;
 
 
-//       std::cout << "\n\n";
-
-      BENCH( sm3.setprod(sm1, sm2); )
-      std::cout << "   a * b:\t" << timer.value() << endl;
-
-//       BENCH(sm3.setprod(sm1.transpose(),sm2); )
-//       std::cout << "   a' * b:\t" << timer.value() << endl;
+//       std::cout << "\n";
 //
-//       BENCH(sm3.setprod(sm1.transpose(),sm2.transpose()); )
-//       std::cout << "   a' * b':\t" << timer.value() << endl;
+//       BENCH( sm3._experimentalNewProduct(sm1, sm2); )
+//       std::cout << "   a * b:\t" << timer.value() << endl;
 //
-//       BENCH(sm3.setprod(sm1, sm2.transpose());)
+//       BENCH(sm3._experimentalNewProduct(sm1.transpose(),sm2); )
+//       std::cout << "   a' * b:\t" << timer.value() << endl;
+// //
+//       BENCH(sm3._experimentalNewProduct(sm1.transpose(),sm2.transpose()); )
+//       std::cout << "   a' * b':\t" << timer.value() << endl;
+// //
+//       BENCH(sm3._experimentalNewProduct(sm1, sm2.transpose());)
 //       std::cout << "   a * b' :\t" << timer.value() << endl;
     }
 
@@ -247,6 +251,23 @@ int main(int argc, char *argv[])
     }
     #endif
 
+    #ifndef NOUBLAS
+    {
+      std::cout << "ublas\t" << nnzPerCol << "%\n";
+      UblasMatrix m1(rows,cols), m2(rows,cols), m3(rows,cols);
+      eiToUblas(sm1, m1);
+      eiToUblas(sm2, m2);
+
+      BENCH(boost::numeric::ublas::prod(m1, m2, m3););
+//       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;
+    }
+    #endif
+
     // GMM++
     #ifndef NOGMM
     {
@@ -263,34 +284,34 @@ int main(int argc, char *argv[])
       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;
-      }
+//       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
 
@@ -309,26 +330,26 @@ int main(int argc, char *argv[])
       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;
+//       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
 
-- 
cgit v1.2.3