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
cs_sorted_multiply(const cs * a,const cs * b)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
cs_sorted_multiply2(const cs * a,const cs * b)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
main(int argc,char * argv[])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
324