1
2//g++ -O3 -g0 -DNDEBUG  sparse_product.cpp -I.. -I/home/gael/Coding/LinearAlgebra/mtl4/ -DDENSITY=0.005 -DSIZE=10000 && ./a.out
3//g++ -O3 -g0 -DNDEBUG  sparse_product.cpp -I.. -I/home/gael/Coding/LinearAlgebra/mtl4/ -DDENSITY=0.05 -DSIZE=2000 && ./a.out
4// -DNOGMM -DNOMTL -DCSPARSE
5// -I /home/gael/Coding/LinearAlgebra/CSparse/Include/ /home/gael/Coding/LinearAlgebra/CSparse/Lib/libcsparse.a
6#ifndef SIZE
7#define SIZE 650000
8#endif
9
10#ifndef DENSITY
11#define DENSITY 0.01
12#endif
13
14#ifndef REPEAT
15#define REPEAT 1
16#endif
17
18#include "BenchSparseUtil.h"
19
20#ifndef MINDENSITY
21#define MINDENSITY 0.0004
22#endif
23
24#ifndef NBTRIES
25#define NBTRIES 10
26#endif
27
28#define BENCH(X) \
29  timer.reset(); \
30  for (int _j=0; _j<NBTRIES; ++_j) { \
31    timer.start(); \
32    for (int _k=0; _k<REPEAT; ++_k) { \
33        X  \
34  } timer.stop(); }
35
36
37#ifdef CSPARSE
38cs* cs_sorted_multiply(const cs* a, const cs* b)
39{
40  cs* A = cs_transpose (a, 1) ;
41  cs* B = cs_transpose (b, 1) ;
42  cs* D = cs_multiply (B,A) ;   /* D = B'*A' */
43  cs_spfree (A) ;
44  cs_spfree (B) ;
45  cs_dropzeros (D) ;      /* drop zeros from D */
46  cs* C = cs_transpose (D, 1) ;   /* C = D', so that C is sorted */
47  cs_spfree (D) ;
48  return C;
49}
50#endif
51
52int main(int argc, char *argv[])
53{
54  int rows = SIZE;
55  int cols = SIZE;
56  float density = DENSITY;
57
58  EigenSparseMatrix sm1(rows,cols);
59  DenseVector v1(cols), v2(cols);
60  v1.setRandom();
61
62  BenchTimer timer;
63  for (float density = DENSITY; density>=MINDENSITY; density*=0.5)
64  {
65    //fillMatrix(density, rows, cols, sm1);
66    fillMatrix2(7, rows, cols, sm1);
67
68    // dense matrices
69    #ifdef DENSEMATRIX
70    {
71      std::cout << "Eigen Dense\t" << density*100 << "%\n";
72      DenseMatrix m1(rows,cols);
73      eiToDense(sm1, m1);
74
75      timer.reset();
76      timer.start();
77      for (int k=0; k<REPEAT; ++k)
78        v2 = m1 * v1;
79      timer.stop();
80      std::cout << "   a * v:\t" << timer.best() << "  " << double(REPEAT)/timer.best() << " * / sec " << endl;
81
82      timer.reset();
83      timer.start();
84      for (int k=0; k<REPEAT; ++k)
85        v2 = m1.transpose() * v1;
86      timer.stop();
87      std::cout << "   a' * v:\t" << timer.best() << endl;
88    }
89    #endif
90
91    // eigen sparse matrices
92    {
93      std::cout << "Eigen sparse\t" << sm1.nonZeros()/float(sm1.rows()*sm1.cols())*100 << "%\n";
94
95      BENCH(asm("#myc"); v2 = sm1 * v1; asm("#myd");)
96      std::cout << "   a * v:\t" << timer.best()/REPEAT << "  " << double(REPEAT)/timer.best(REAL_TIMER) << " * / sec " << endl;
97
98
99      BENCH( { asm("#mya"); v2 = sm1.transpose() * v1; asm("#myb"); })
100
101      std::cout << "   a' * v:\t" << timer.best()/REPEAT << endl;
102    }
103
104//     {
105//       DynamicSparseMatrix<Scalar> m1(sm1);
106//       std::cout << "Eigen dyn-sparse\t" << m1.nonZeros()/float(m1.rows()*m1.cols())*100 << "%\n";
107//
108//       BENCH(for (int k=0; k<REPEAT; ++k) v2 = m1 * v1;)
109//       std::cout << "   a * v:\t" << timer.value() << endl;
110//
111//       BENCH(for (int k=0; k<REPEAT; ++k) v2 = m1.transpose() * v1;)
112//       std::cout << "   a' * v:\t" << timer.value() << endl;
113//     }
114
115    // GMM++
116    #ifndef NOGMM
117    {
118      std::cout << "GMM++ sparse\t" << density*100 << "%\n";
119      //GmmDynSparse  gmmT3(rows,cols);
120      GmmSparse m1(rows,cols);
121      eiToGmm(sm1, m1);
122
123      std::vector<Scalar> gmmV1(cols), gmmV2(cols);
124      Map<Matrix<Scalar,Dynamic,1> >(&gmmV1[0], cols) = v1;
125      Map<Matrix<Scalar,Dynamic,1> >(&gmmV2[0], cols) = v2;
126
127      BENCH( asm("#myx"); gmm::mult(m1, gmmV1, gmmV2); asm("#myy"); )
128      std::cout << "   a * v:\t" << timer.value() << endl;
129
130      BENCH( gmm::mult(gmm::transposed(m1), gmmV1, gmmV2); )
131      std::cout << "   a' * v:\t" << timer.value() << endl;
132    }
133    #endif
134
135    #ifndef NOUBLAS
136    {
137      std::cout << "ublas sparse\t" << density*100 << "%\n";
138      UBlasSparse m1(rows,cols);
139      eiToUblas(sm1, m1);
140
141      boost::numeric::ublas::vector<Scalar> uv1, uv2;
142      eiToUblasVec(v1,uv1);
143      eiToUblasVec(v2,uv2);
144
145//       std::vector<Scalar> gmmV1(cols), gmmV2(cols);
146//       Map<Matrix<Scalar,Dynamic,1> >(&gmmV1[0], cols) = v1;
147//       Map<Matrix<Scalar,Dynamic,1> >(&gmmV2[0], cols) = v2;
148
149      BENCH( uv2 = boost::numeric::ublas::prod(m1, uv1); )
150      std::cout << "   a * v:\t" << timer.value() << endl;
151
152//       BENCH( boost::ublas::prod(gmm::transposed(m1), gmmV1, gmmV2); )
153//       std::cout << "   a' * v:\t" << timer.value() << endl;
154    }
155    #endif
156
157    // MTL4
158    #ifndef NOMTL
159    {
160      std::cout << "MTL4\t" << density*100 << "%\n";
161      MtlSparse m1(rows,cols);
162      eiToMtl(sm1, m1);
163      mtl::dense_vector<Scalar> mtlV1(cols, 1.0);
164      mtl::dense_vector<Scalar> mtlV2(cols, 1.0);
165
166      timer.reset();
167      timer.start();
168      for (int k=0; k<REPEAT; ++k)
169        mtlV2 = m1 * mtlV1;
170      timer.stop();
171      std::cout << "   a * v:\t" << timer.value() << endl;
172
173      timer.reset();
174      timer.start();
175      for (int k=0; k<REPEAT; ++k)
176        mtlV2 = trans(m1) * mtlV1;
177      timer.stop();
178      std::cout << "   a' * v:\t" << timer.value() << endl;
179    }
180    #endif
181
182    std::cout << "\n\n";
183  }
184
185  return 0;
186}
187
188