source: pacpussensors/trunk/Vislab/lib3dv/eigen/bench/sparse_product.cpp@ 136

Last change on this file since 136 was 136, checked in by ldecherf, 7 years ago

Doc
File size: 8.8 KB
Line 
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
64cs* 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
83cs* 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
94void bench_sort();
95
96int 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
Note: See TracBrowser for help on using the repository browser.