bench/sparse_product.cpp
Go to the documentation of this file.
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 
7 #include <typeinfo>
8 
9 #ifndef SIZE
10 #define SIZE 1000000
11 #endif
12 
13 #ifndef NNZPERCOL
14 #define NNZPERCOL 6
15 #endif
16 
17 #ifndef REPEAT
18 #define REPEAT 1
19 #endif
20 
21 #include <algorithm>
22 #include "BenchTimer.h"
23 #include "BenchUtil.h"
24 #include "BenchSparseUtil.h"
25 
26 #ifndef NBTRIES
27 #define NBTRIES 1
28 #endif
29 
30 #define BENCH(X) \
31  timer.reset(); \
32  for (int _j=0; _j<NBTRIES; ++_j) { \
33  timer.start(); \
34  for (int _k=0; _k<REPEAT; ++_k) { \
35  X \
36  } timer.stop(); }
37 
38 // #ifdef MKL
39 //
40 // #include "mkl_types.h"
41 // #include "mkl_spblas.h"
42 //
43 // template<typename Lhs,typename Rhs,typename Res>
44 // void mkl_multiply(const Lhs& lhs, const Rhs& rhs, Res& res)
45 // {
46 // char n = 'N';
47 // float alpha = 1;
48 // char matdescra[6];
49 // matdescra[0] = 'G';
50 // matdescra[1] = 0;
51 // matdescra[2] = 0;
52 // matdescra[3] = 'C';
53 // mkl_scscmm(&n, lhs.rows(), rhs.cols(), lhs.cols(), &alpha, matdescra,
54 // lhs._valuePtr(), lhs._innerIndexPtr(), lhs.outerIndexPtr(),
55 // pntre, b, &ldb, &beta, c, &ldc);
56 // // mkl_somatcopy('C', 'T', lhs.rows(), lhs.cols(), 1,
57 // // lhs._valuePtr(), lhs.rows(), DST, dst_stride);
58 // }
59 //
60 // #endif
61 
62 
63 #ifdef CSPARSE
64 cs* cs_sorted_multiply(const cs* a, const cs* b)
65 {
66 // return cs_multiply(a,b);
67 
68  cs* A = cs_transpose(a, 1);
69  cs* B = cs_transpose(b, 1);
70  cs* D = cs_multiply(B,A); /* D = B'*A' */
71  cs_spfree (A) ;
72  cs_spfree (B) ;
73  cs_dropzeros (D) ; /* drop zeros from D */
74  cs* C = cs_transpose (D, 1) ; /* C = D', so that C is sorted */
75  cs_spfree (D) ;
76  return C;
77 
78 // cs* A = cs_transpose(a, 1);
79 // cs* C = cs_transpose(A, 1);
80 // return C;
81 }
82 
83 cs* cs_sorted_multiply2(const cs* a, const cs* b)
84 {
85  cs* D = cs_multiply(a,b);
86  cs* E = cs_transpose(D,1);
87  cs_spfree(D);
88  cs* C = cs_transpose(E,1);
89  cs_spfree(E);
90  return C;
91 }
92 #endif
93 
94 void bench_sort();
95 
96 int main(int argc, char *argv[])
97 {
98 // bench_sort();
99 
100  int rows = SIZE;
101  int cols = SIZE;
102  float density = DENSITY;
103 
104  EigenSparseMatrix sm1(rows,cols), sm2(rows,cols), sm3(rows,cols), sm4(rows,cols);
105 
106  BenchTimer timer;
107  for (int nnzPerCol = NNZPERCOL; nnzPerCol>1; nnzPerCol/=1.1)
108  {
109  sm1.setZero();
110  sm2.setZero();
111  fillMatrix2(nnzPerCol, rows, cols, sm1);
112  fillMatrix2(nnzPerCol, rows, cols, sm2);
113 // std::cerr << "filling OK\n";
114 
115  // dense matrices
116  #ifdef DENSEMATRIX
117  {
118  std::cout << "Eigen Dense\t" << nnzPerCol << "%\n";
119  DenseMatrix m1(rows,cols), m2(rows,cols), m3(rows,cols);
120  eiToDense(sm1, m1);
121  eiToDense(sm2, m2);
122 
123  timer.reset();
124  timer.start();
125  for (int k=0; k<REPEAT; ++k)
126  m3 = m1 * m2;
127  timer.stop();
128  std::cout << " a * b:\t" << timer.value() << endl;
129 
130  timer.reset();
131  timer.start();
132  for (int k=0; k<REPEAT; ++k)
133  m3 = m1.transpose() * m2;
134  timer.stop();
135  std::cout << " a' * b:\t" << timer.value() << endl;
136 
137  timer.reset();
138  timer.start();
139  for (int k=0; k<REPEAT; ++k)
140  m3 = m1.transpose() * m2.transpose();
141  timer.stop();
142  std::cout << " a' * b':\t" << timer.value() << endl;
143 
144  timer.reset();
145  timer.start();
146  for (int k=0; k<REPEAT; ++k)
147  m3 = m1 * m2.transpose();
148  timer.stop();
149  std::cout << " a * b':\t" << timer.value() << endl;
150  }
151  #endif
152 
153  // eigen sparse matrices
154  {
155  std::cout << "Eigen sparse\t" << sm1.nonZeros()/(float(sm1.rows())*float(sm1.cols()))*100 << "% * "
156  << sm2.nonZeros()/(float(sm2.rows())*float(sm2.cols()))*100 << "%\n";
157 
158  BENCH(sm3 = sm1 * sm2; )
159  std::cout << " a * b:\t" << timer.value() << endl;
160 
161 // BENCH(sm3 = sm1.transpose() * sm2; )
162 // std::cout << " a' * b:\t" << timer.value() << endl;
163 // //
164 // BENCH(sm3 = sm1.transpose() * sm2.transpose(); )
165 // std::cout << " a' * b':\t" << timer.value() << endl;
166 // //
167 // BENCH(sm3 = sm1 * sm2.transpose(); )
168 // std::cout << " a * b' :\t" << timer.value() << endl;
169 
170 
171 // std::cout << "\n";
172 //
173 // BENCH( sm3._experimentalNewProduct(sm1, sm2); )
174 // std::cout << " a * b:\t" << timer.value() << endl;
175 //
176 // BENCH(sm3._experimentalNewProduct(sm1.transpose(),sm2); )
177 // std::cout << " a' * b:\t" << timer.value() << endl;
178 // //
179 // BENCH(sm3._experimentalNewProduct(sm1.transpose(),sm2.transpose()); )
180 // std::cout << " a' * b':\t" << timer.value() << endl;
181 // //
182 // BENCH(sm3._experimentalNewProduct(sm1, sm2.transpose());)
183 // std::cout << " a * b' :\t" << timer.value() << endl;
184  }
185 
186  // eigen dyn-sparse matrices
187  /*{
188  DynamicSparseMatrix<Scalar> m1(sm1), m2(sm2), m3(sm3);
189  std::cout << "Eigen dyn-sparse\t" << m1.nonZeros()/(float(m1.rows())*float(m1.cols()))*100 << "% * "
190  << m2.nonZeros()/(float(m2.rows())*float(m2.cols()))*100 << "%\n";
191 
192 // timer.reset();
193 // timer.start();
194  BENCH(for (int k=0; k<REPEAT; ++k) m3 = m1 * m2;)
195 // timer.stop();
196  std::cout << " a * b:\t" << timer.value() << endl;
197 // std::cout << sm3 << "\n";
198 
199  timer.reset();
200  timer.start();
201 // std::cerr << "transpose...\n";
202 // EigenSparseMatrix sm4 = sm1.transpose();
203 // std::cout << sm4.nonZeros() << " == " << sm1.nonZeros() << "\n";
204 // exit(1);
205 // std::cerr << "transpose OK\n";
206 // std::cout << sm1 << "\n\n" << sm1.transpose() << "\n\n" << sm4.transpose() << "\n\n";
207  BENCH(for (int k=0; k<REPEAT; ++k) m3 = m1.transpose() * m2;)
208 // timer.stop();
209  std::cout << " a' * b:\t" << timer.value() << endl;
210 
211 // timer.reset();
212 // timer.start();
213  BENCH( for (int k=0; k<REPEAT; ++k) m3 = m1.transpose() * m2.transpose(); )
214 // timer.stop();
215  std::cout << " a' * b':\t" << timer.value() << endl;
216 
217 // timer.reset();
218 // timer.start();
219  BENCH( for (int k=0; k<REPEAT; ++k) m3 = m1 * m2.transpose(); )
220 // timer.stop();
221  std::cout << " a * b' :\t" << timer.value() << endl;
222  }*/
223 
224  // CSparse
225  #ifdef CSPARSE
226  {
227  std::cout << "CSparse \t" << nnzPerCol << "%\n";
228  cs *m1, *m2, *m3;
229  eiToCSparse(sm1, m1);
230  eiToCSparse(sm2, m2);
231 
232  BENCH(
233  {
234  m3 = cs_sorted_multiply(m1, m2);
235  if (!m3)
236  {
237  std::cerr << "cs_multiply failed\n";
238  }
239 // cs_print(m3, 0);
240  cs_spfree(m3);
241  }
242  );
243 // timer.stop();
244  std::cout << " a * b:\t" << timer.value() << endl;
245 
246 // BENCH( { m3 = cs_sorted_multiply2(m1, m2); cs_spfree(m3); } );
247 // std::cout << " a * b:\t" << timer.value() << endl;
248  }
249  #endif
250 
251  #ifndef NOUBLAS
252  {
253  std::cout << "ublas\t" << nnzPerCol << "%\n";
254  UBlasSparse m1(rows,cols), m2(rows,cols), m3(rows,cols);
255  eiToUblas(sm1, m1);
256  eiToUblas(sm2, m2);
257 
258  BENCH(boost::numeric::ublas::prod(m1, m2, m3););
259  std::cout << " a * b:\t" << timer.value() << endl;
260  }
261  #endif
262 
263  // GMM++
264  #ifndef NOGMM
265  {
266  std::cout << "GMM++ sparse\t" << nnzPerCol << "%\n";
267  GmmDynSparse gmmT3(rows,cols);
268  GmmSparse m1(rows,cols), m2(rows,cols), m3(rows,cols);
269  eiToGmm(sm1, m1);
270  eiToGmm(sm2, m2);
271 
272  BENCH(gmm::mult(m1, m2, gmmT3););
273  std::cout << " a * b:\t" << timer.value() << endl;
274 
275 // BENCH(gmm::mult(gmm::transposed(m1), m2, gmmT3););
276 // std::cout << " a' * b:\t" << timer.value() << endl;
277 //
278 // if (rows<500)
279 // {
280 // BENCH(gmm::mult(gmm::transposed(m1), gmm::transposed(m2), gmmT3););
281 // std::cout << " a' * b':\t" << timer.value() << endl;
282 //
283 // BENCH(gmm::mult(m1, gmm::transposed(m2), gmmT3););
284 // std::cout << " a * b':\t" << timer.value() << endl;
285 // }
286 // else
287 // {
288 // std::cout << " a' * b':\t" << "forever" << endl;
289 // std::cout << " a * b':\t" << "forever" << endl;
290 // }
291  }
292  #endif
293 
294  // MTL4
295  #ifndef NOMTL
296  {
297  std::cout << "MTL4\t" << nnzPerCol << "%\n";
298  MtlSparse m1(rows,cols), m2(rows,cols), m3(rows,cols);
299  eiToMtl(sm1, m1);
300  eiToMtl(sm2, m2);
301 
302  BENCH(m3 = m1 * m2;);
303  std::cout << " a * b:\t" << timer.value() << endl;
304 
305 // BENCH(m3 = trans(m1) * m2;);
306 // std::cout << " a' * b:\t" << timer.value() << endl;
307 //
308 // BENCH(m3 = trans(m1) * trans(m2););
309 // std::cout << " a' * b':\t" << timer.value() << endl;
310 //
311 // BENCH(m3 = m1 * trans(m2););
312 // std::cout << " a * b' :\t" << timer.value() << endl;
313  }
314  #endif
315 
316  std::cout << "\n\n";
317  }
318 
319  return 0;
320 }
321 
322 
323 
REPEAT
#define REPEAT
Definition: bench/sparse_product.cpp:18
Eigen::SparseMatrix< Scalar >
B
Matrix< SCALARB, Dynamic, Dynamic, opt_B > B
Definition: bench_gemm.cpp:49
MtlSparse
mtl::compressed2D< Scalar, mtl::matrix::parameters< mtl::tag::col_major > > MtlSparse
Definition: BenchSparseUtil.h:86
D
MatrixXcd D
Definition: EigenSolver_EigenSolver_MatrixType.cpp:14
eiToDense
void eiToDense(const EigenSparseMatrix &src, DenseMatrix &dst)
Definition: BenchSparseUtil.h:62
bench_sort
void bench_sort()
timer
static BenchTimer timer
Definition: benchmark-blocking-sizes.cpp:31
b
Scalar * b
Definition: benchVecAdd.cpp:17
m1
Matrix3d m1
Definition: IOFormat.cpp:2
main
int main(int argc, char *argv[])
Definition: bench/sparse_product.cpp:96
DENSITY
#define DENSITY
Definition: BenchSparseUtil.h:15
eiToMtl
void eiToMtl(const EigenSparseMatrix &src, MtlSparse &dst)
Definition: BenchSparseUtil.h:88
rows
int rows
Definition: Tutorial_commainit_02.cpp:1
SIZE
#define SIZE
Definition: bench/sparse_product.cpp:10
A
Matrix< SCALARA, Dynamic, Dynamic, opt_A > A
Definition: bench_gemm.cpp:48
Eigen::BenchTimer::value
double value(int TIMER=CPU_TIMER) const
Definition: BenchTimer.h:104
BenchUtil.h
m2
MatrixType m2(n_dims)
Eigen::BenchTimer::start
void start()
Definition: BenchTimer.h:81
GmmDynSparse
gmm::col_matrix< gmm::wsvector< Scalar > > GmmDynSparse
Definition: BenchSparseUtil.h:73
Eigen::BenchTimer
Definition: BenchTimer.h:59
density
Definition: testGaussianConditional.cpp:127
NNZPERCOL
#define NNZPERCOL
Definition: bench/sparse_product.cpp:14
BENCH
#define BENCH(X)
Definition: bench/sparse_product.cpp:30
eiToUblas
void eiToUblas(const EigenSparseMatrix &src, UBlasSparse &dst)
Definition: BenchSparseUtil.h:128
Eigen::BenchTimer::reset
void reset()
Definition: BenchTimer.h:75
E
DiscreteKey E(5, 2)
m3
static const DiscreteKey m3(M(3), 2)
a
ArrayXXi a
Definition: Array_initializer_list_23_cxx11.cpp:1
C
Matrix< Scalar, Dynamic, Dynamic > C
Definition: bench_gemm.cpp:50
BenchTimer.h
GmmSparse
gmm::csc_matrix< Scalar > GmmSparse
Definition: BenchSparseUtil.h:72
DenseMatrix
Matrix< Scalar, Dynamic, Dynamic > DenseMatrix
Definition: BenchSparseUtil.h:23
Eigen::BenchTimer::stop
void stop()
Definition: BenchTimer.h:86
gtsam.examples.DogLegOptimizerExample.float
float
Definition: DogLegOptimizerExample.py:113
eiToGmm
void eiToGmm(const EigenSparseMatrix &src, GmmSparse &dst)
Definition: BenchSparseUtil.h:74
UBlasSparse
boost::numeric::ublas::compressed_matrix< Scalar, boost::numeric::ublas::column_major > UBlasSparse
Definition: BenchSparseUtil.h:126
cols
int cols
Definition: Tutorial_commainit_02.cpp:1
prod
EIGEN_DONT_INLINE void prod(const Lhs &a, const Rhs &b, Res &c)
Definition: product_threshold.cpp:39
BenchSparseUtil.h
fillMatrix2
void fillMatrix2(int nnzPerCol, int rows, int cols, EigenSparseMatrix &dst)
Definition: BenchSparseUtil.h:42

gtsam
Author(s):
autogenerated on Fri Nov 1 2024 03:35:33