add support for ublas

4 files changed, 70 insertions, 65 deletions

diff --git a/bench/BenchSparseUtil.h b/bench/BenchSparseUtil.h
index ff836bffe..13981f6b7 100644
--- a/bench/BenchSparseUtil.h
+++ b/bench/BenchSparseUtil.h
@@ -26,7 +26,7 @@ typedef SparseMatrix<Scalar> EigenSparseMatrix;
 
 void fillMatrix(float density, int rows, int cols,  EigenSparseMatrix& dst)
 {
-  dst.reserve(rows*cols*density);
+  dst.reserve(double(rows)*cols*density);
   for(int j = 0; j < cols; j++)
   {
     for(int i = 0; i < rows; i++)
@@ -122,22 +122,24 @@ void eiToCSparse(const EigenSparseMatrix& src, cs* &dst)
 #include <boost/numeric/ublas/matrix_sparse.hpp>
 #include <boost/numeric/ublas/vector_of_vector.hpp>
 #include <boost/numeric/ublas/operation.hpp>
-// #include <boost/numeric/ublas/sparse_prod.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> UBlasSparse;
 
-typedef boost::numeric::ublas::compressed_matrix<Scalar,boost::numeric::ublas::column_major> UblasMatrix;
-
-void eiToUblas(const EigenSparseMatrix& src, UblasMatrix& dst)
+void eiToUblas(const EigenSparseMatrix& src, UBlasSparse& dst)
 {
+  dst.resize(src.rows(), src.cols(), false);
   for (int j=0; j<src.cols(); ++j)
     for (EigenSparseMatrix::InnerIterator it(src.derived(), j); it; ++it)
       dst(it.index(),j) = it.value();
 }
 
-
+template <typename EigenType, typename UblasType>
+void eiToUblasVec(const EigenType& src, UblasType& dst)
+{
+  dst.resize(src.size());
+  for (int j=0; j<src.size(); ++j)
+      dst[j] = src.coeff(j);
+}
 #endif
 
 #ifdef OSKI
diff --git a/bench/BenchUtil.h b/bench/BenchUtil.h
index 9798fa385..8883a1380 100644
--- a/bench/BenchUtil.h
+++ b/bench/BenchUtil.h
@@ -69,4 +69,24 @@ void eiToGsl(const EigenMatrixType& src, gsl_matrix** dst)
 }
 #endif
 
+#ifdef BENCH_UBLAS
+#include <boost/numeric/ublas/matrix.hpp>
+#include <boost/numeric/ublas/vector.hpp>
+template <typename EigenMatrixType, typename UblasMatrixType>
+void eiToUblas(const EigenMatrixType& src, UblasMatrixType& dst)
+{
+  dst.resize(src.rows(),src.cols());
+  for (int j=0; j<src.cols(); ++j)
+    for (int i=0; i<src.rows(); ++i)
+      dst(i,j) = src.coeff(i,j);
+}
+template <typename EigenType, typename UblasType>
+void eiToUblasVec(const EigenType& src, UblasType& dst)
+{
+  dst.resize(src.size());
+  for (int j=0; j<src.size(); ++j)
+      dst[j] = src.coeff(j);
+}
+#endif
+
 #endif // EIGEN_BENCH_UTIL_H
diff --git a/bench/sparse_dense_product.cpp b/bench/sparse_dense_product.cpp
index bfe46122d..f3f519406 100644
--- a/bench/sparse_dense_product.cpp
+++ b/bench/sparse_dense_product.cpp
@@ -4,7 +4,7 @@
 // -DNOGMM -DNOMTL -DCSPARSE
 // -I /home/gael/Coding/LinearAlgebra/CSparse/Include/ /home/gael/Coding/LinearAlgebra/CSparse/Lib/libcsparse.a
 #ifndef SIZE
-#define SIZE 10000
+#define SIZE 650000
 #endif
 
 #ifndef DENSITY
@@ -62,7 +62,8 @@ int main(int argc, char *argv[])
   BenchTimer timer;
   for (float density = DENSITY; density>=MINDENSITY; density*=0.5)
   {
-    fillMatrix(density, rows, cols, sm1);
+    //fillMatrix(density, rows, cols, sm1);
+    fillMatrix2(7, rows, cols, sm1);
 
     // dense matrices
     #ifdef DENSEMATRIX
@@ -76,14 +77,14 @@ int main(int argc, char *argv[])
       for (int k=0; k<REPEAT; ++k)
         v2 = m1 * v1;
       timer.stop();
-      std::cout << "   a * v:\t" << timer.value() << endl;
+      std::cout << "   a * v:\t" << timer.best() << "  " << double(REPEAT)/timer.best() << " * / sec " << endl;
 
       timer.reset();
       timer.start();
       for (int k=0; k<REPEAT; ++k)
         v2 = m1.transpose() * v1;
       timer.stop();
-      std::cout << "   a' * v:\t" << timer.value() << endl;
+      std::cout << "   a' * v:\t" << timer.best() << endl;
     }
     #endif
 
@@ -92,12 +93,12 @@ int main(int argc, char *argv[])
       std::cout << "Eigen sparse\t" << sm1.nonZeros()/float(sm1.rows()*sm1.cols())*100 << "%\n";
 
       BENCH(asm("#myc"); v2 = sm1 * v1; asm("#myd");)
-      std::cout << "   a * v:\t" << timer.value() << endl;
+      std::cout << "   a * v:\t" << timer.best()/REPEAT << "  " << double(REPEAT)/timer.best(REAL_TIMER) << " * / sec " << endl;
 
 
       BENCH( { asm("#mya"); v2 = sm1.transpose() * v1; asm("#myb"); })
 
-      std::cout << "   a' * v:\t" << timer.value() << endl;
+      std::cout << "   a' * v:\t" << timer.best()/REPEAT << endl;
     }
 
 //     {
@@ -130,6 +131,28 @@ int main(int argc, char *argv[])
       std::cout << "   a' * v:\t" << timer.value() << endl;
     }
     #endif
+    
+    #ifndef NOUBLAS
+    {
+      std::cout << "ublas sparse\t" << density*100 << "%\n";
+      UBlasSparse m1(rows,cols);
+      eiToUblas(sm1, m1);
+      
+      boost::numeric::ublas::vector<Scalar> uv1, uv2;
+      eiToUblasVec(v1,uv1);
+      eiToUblasVec(v2,uv2);
+
+//       std::vector<Scalar> gmmV1(cols), gmmV2(cols);
+//       Map<Matrix<Scalar,Dynamic,1> >(&gmmV1[0], cols) = v1;
+//       Map<Matrix<Scalar,Dynamic,1> >(&gmmV2[0], cols) = v2;
+
+      BENCH( uv2 = boost::numeric::ublas::prod(m1, uv1); )
+      std::cout << "   a * v:\t" << timer.value() << endl;
+
+//       BENCH( boost::ublas::prod(gmm::transposed(m1), gmmV1, gmmV2); )
+//       std::cout << "   a' * v:\t" << timer.value() << endl;
+    }
+    #endif
 
     // MTL4
     #ifndef NOMTL
diff --git a/bench/sparse_product.cpp b/bench/sparse_product.cpp
index 0b5558b89..d2fc44f0d 100644
--- a/bench/sparse_product.cpp
+++ b/bench/sparse_product.cpp
@@ -20,6 +20,7 @@
 
 #include <algorithm>
 #include "BenchTimer.h"
+#include "BenchUtil.h"
 #include "BenchSparseUtil.h"
 
 #ifndef NBTRIES
@@ -228,16 +229,12 @@ int main(int argc, char *argv[])
       eiToCSparse(sm1, m1);
       eiToCSparse(sm2, m2);
 
-//       timer.reset();
-//       timer.start();
-//       for (int k=0; k<REPEAT; ++k)
       BENCH(
       {
         m3 = cs_sorted_multiply(m1, m2);
         if (!m3)
         {
           std::cerr << "cs_multiply failed\n";
-//           break;
         }
 //         cs_print(m3, 0);
         cs_spfree(m3);
@@ -254,16 +251,11 @@ int main(int argc, char *argv[])
     #ifndef NOUBLAS
     {
       std::cout << "ublas\t" << nnzPerCol << "%\n";
-      UblasMatrix m1(rows,cols), m2(rows,cols), m3(rows,cols);
+      UBlasSparse 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
@@ -277,34 +269,18 @@ int main(int argc, char *argv[])
       eiToGmm(sm1, m1);
       eiToGmm(sm2, m2);
 
-      timer.reset();
-      timer.start();
-      for (int k=0; k<REPEAT; ++k)
-        gmm::mult(m1, m2, gmmT3);
-      timer.stop();
+      BENCH(gmm::mult(m1, m2, gmmT3););
       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();
+//       BENCH(gmm::mult(gmm::transposed(m1), m2, gmmT3););
 //       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();
+//         BENCH(gmm::mult(gmm::transposed(m1), gmm::transposed(m2), gmmT3););
 //         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();
+//         BENCH(gmm::mult(m1, gmm::transposed(m2), gmmT3););
 //         std::cout << "   a * b':\t" << timer.value() << endl;
 //       }
 //       else
@@ -323,32 +299,16 @@ int main(int argc, char *argv[])
       eiToMtl(sm1, m1);
       eiToMtl(sm2, m2);
 
-      timer.reset();
-      timer.start();
-      for (int k=0; k<REPEAT; ++k)
-        m3 = m1 * m2;
-      timer.stop();
+      BENCH(m3 = m1 * m2;);
       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();
+//       BENCH(m3 = trans(m1) * m2;);
 //       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();
+//       BENCH(m3 = trans(m1) * trans(m2););
 //       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();
+//       BENCH(m3 = m1 * trans(m2););
 //       std::cout << "   a * b' :\t" << timer.value() << endl;
     }
     #endif