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 100000
8 #endif
9
10 #ifndef NBPERROW
11 #define NBPERROW 24
12 #endif
13
14 #ifndef REPEAT
15 #define REPEAT 2
16 #endif
17
18 #ifndef NBTRIES
19 #define NBTRIES 2
20 #endif
21
22 #ifndef KK
23 #define KK 10
24 #endif
25
26 #ifndef NOGOOGLE
27 #define EIGEN_GOOGLEHASH_SUPPORT
28 #include <google/sparse_hash_map>
29 #endif
30
31 #include "BenchSparseUtil.h"
32
33 #define CHECK_MEM
34 // #define CHECK_MEM std/**/::cout << "check mem\n"; getchar();
35
36 #define BENCH(X) \
37 timer.reset(); \
38 for (int _j=0; _j<NBTRIES; ++_j) { \
39 timer.start(); \
40 for (int _k=0; _k<REPEAT; ++_k) { \
41 X \
42 } timer.stop(); }
43
44 typedef std::vector<Vector2i> Coordinates;
45 typedef std::vector<float> Values;
46
47 EIGEN_DONT_INLINE Scalar* setinnerrand_eigen(const Coordinates& coords, const Values& vals);
48 EIGEN_DONT_INLINE Scalar* setrand_eigen_dynamic(const Coordinates& coords, const Values& vals);
49 EIGEN_DONT_INLINE Scalar* setrand_eigen_compact(const Coordinates& coords, const Values& vals);
50 EIGEN_DONT_INLINE Scalar* setrand_eigen_sumeq(const Coordinates& coords, const Values& vals);
51 EIGEN_DONT_INLINE Scalar* setrand_eigen_gnu_hash(const Coordinates& coords, const Values& vals);
52 EIGEN_DONT_INLINE Scalar* setrand_eigen_google_dense(const Coordinates& coords, const Values& vals);
53 EIGEN_DONT_INLINE Scalar* setrand_eigen_google_sparse(const Coordinates& coords, const Values& vals);
54 EIGEN_DONT_INLINE Scalar* setrand_scipy(const Coordinates& coords, const Values& vals);
55 EIGEN_DONT_INLINE Scalar* setrand_ublas_mapped(const Coordinates& coords, const Values& vals);
56 EIGEN_DONT_INLINE Scalar* setrand_ublas_coord(const Coordinates& coords, const Values& vals);
57 EIGEN_DONT_INLINE Scalar* setrand_ublas_compressed(const Coordinates& coords, const Values& vals);
58 EIGEN_DONT_INLINE Scalar* setrand_ublas_genvec(const Coordinates& coords, const Values& vals);
59 EIGEN_DONT_INLINE Scalar* setrand_mtl(const Coordinates& coords, const Values& vals);
60
main(int argc,char * argv[])61 int main(int argc, char *argv[])
62 {
63 int rows = SIZE;
64 int cols = SIZE;
65 bool fullyrand = true;
66
67 BenchTimer timer;
68 Coordinates coords;
69 Values values;
70 if(fullyrand)
71 {
72 Coordinates pool;
73 pool.reserve(cols*NBPERROW);
74 std::cerr << "fill pool" << "\n";
75 for (int i=0; i<cols*NBPERROW; )
76 {
77 // DynamicSparseMatrix<int> stencil(SIZE,SIZE);
78 Vector2i ij(internal::random<int>(0,rows-1),internal::random<int>(0,cols-1));
79 // if(stencil.coeffRef(ij.x(), ij.y())==0)
80 {
81 // stencil.coeffRef(ij.x(), ij.y()) = 1;
82 pool.push_back(ij);
83
84 }
85 ++i;
86 }
87 std::cerr << "pool ok" << "\n";
88 int n = cols*NBPERROW*KK;
89 coords.reserve(n);
90 values.reserve(n);
91 for (int i=0; i<n; ++i)
92 {
93 int i = internal::random<int>(0,pool.size());
94 coords.push_back(pool[i]);
95 values.push_back(internal::random<Scalar>());
96 }
97 }
98 else
99 {
100 for (int j=0; j<cols; ++j)
101 for (int i=0; i<NBPERROW; ++i)
102 {
103 coords.push_back(Vector2i(internal::random<int>(0,rows-1),j));
104 values.push_back(internal::random<Scalar>());
105 }
106 }
107 std::cout << "nnz = " << coords.size() << "\n";
108 CHECK_MEM
109
110 // dense matrices
111 #ifdef DENSEMATRIX
112 {
113 BENCH(setrand_eigen_dense(coords,values);)
114 std::cout << "Eigen Dense\t" << timer.value() << "\n";
115 }
116 #endif
117
118 // eigen sparse matrices
119 // if (!fullyrand)
120 // {
121 // BENCH(setinnerrand_eigen(coords,values);)
122 // std::cout << "Eigen fillrand\t" << timer.value() << "\n";
123 // }
124 {
125 BENCH(setrand_eigen_dynamic(coords,values);)
126 std::cout << "Eigen dynamic\t" << timer.value() << "\n";
127 }
128 // {
129 // BENCH(setrand_eigen_compact(coords,values);)
130 // std::cout << "Eigen compact\t" << timer.value() << "\n";
131 // }
132 {
133 BENCH(setrand_eigen_sumeq(coords,values);)
134 std::cout << "Eigen sumeq\t" << timer.value() << "\n";
135 }
136 {
137 // BENCH(setrand_eigen_gnu_hash(coords,values);)
138 // std::cout << "Eigen std::map\t" << timer.value() << "\n";
139 }
140 {
141 BENCH(setrand_scipy(coords,values);)
142 std::cout << "scipy\t" << timer.value() << "\n";
143 }
144 #ifndef NOGOOGLE
145 {
146 BENCH(setrand_eigen_google_dense(coords,values);)
147 std::cout << "Eigen google dense\t" << timer.value() << "\n";
148 }
149 {
150 BENCH(setrand_eigen_google_sparse(coords,values);)
151 std::cout << "Eigen google sparse\t" << timer.value() << "\n";
152 }
153 #endif
154
155 #ifndef NOUBLAS
156 {
157 // BENCH(setrand_ublas_mapped(coords,values);)
158 // std::cout << "ublas mapped\t" << timer.value() << "\n";
159 }
160 {
161 BENCH(setrand_ublas_genvec(coords,values);)
162 std::cout << "ublas vecofvec\t" << timer.value() << "\n";
163 }
164 /*{
165 timer.reset();
166 timer.start();
167 for (int k=0; k<REPEAT; ++k)
168 setrand_ublas_compressed(coords,values);
169 timer.stop();
170 std::cout << "ublas comp\t" << timer.value() << "\n";
171 }
172 {
173 timer.reset();
174 timer.start();
175 for (int k=0; k<REPEAT; ++k)
176 setrand_ublas_coord(coords,values);
177 timer.stop();
178 std::cout << "ublas coord\t" << timer.value() << "\n";
179 }*/
180 #endif
181
182
183 // MTL4
184 #ifndef NOMTL
185 {
186 BENCH(setrand_mtl(coords,values));
187 std::cout << "MTL\t" << timer.value() << "\n";
188 }
189 #endif
190
191 return 0;
192 }
193
setinnerrand_eigen(const Coordinates & coords,const Values & vals)194 EIGEN_DONT_INLINE Scalar* setinnerrand_eigen(const Coordinates& coords, const Values& vals)
195 {
196 using namespace Eigen;
197 SparseMatrix<Scalar> mat(SIZE,SIZE);
198 //mat.startFill(2000000/*coords.size()*/);
199 for (int i=0; i<coords.size(); ++i)
200 {
201 mat.insert(coords[i].x(), coords[i].y()) = vals[i];
202 }
203 mat.finalize();
204 CHECK_MEM;
205 return 0;
206 }
207
setrand_eigen_dynamic(const Coordinates & coords,const Values & vals)208 EIGEN_DONT_INLINE Scalar* setrand_eigen_dynamic(const Coordinates& coords, const Values& vals)
209 {
210 using namespace Eigen;
211 DynamicSparseMatrix<Scalar> mat(SIZE,SIZE);
212 mat.reserve(coords.size()/10);
213 for (int i=0; i<coords.size(); ++i)
214 {
215 mat.coeffRef(coords[i].x(), coords[i].y()) += vals[i];
216 }
217 mat.finalize();
218 CHECK_MEM;
219 return &mat.coeffRef(coords[0].x(), coords[0].y());
220 }
221
setrand_eigen_sumeq(const Coordinates & coords,const Values & vals)222 EIGEN_DONT_INLINE Scalar* setrand_eigen_sumeq(const Coordinates& coords, const Values& vals)
223 {
224 using namespace Eigen;
225 int n = coords.size()/KK;
226 DynamicSparseMatrix<Scalar> mat(SIZE,SIZE);
227 for (int j=0; j<KK; ++j)
228 {
229 DynamicSparseMatrix<Scalar> aux(SIZE,SIZE);
230 mat.reserve(n);
231 for (int i=j*n; i<(j+1)*n; ++i)
232 {
233 aux.insert(coords[i].x(), coords[i].y()) += vals[i];
234 }
235 aux.finalize();
236 mat += aux;
237 }
238 return &mat.coeffRef(coords[0].x(), coords[0].y());
239 }
240
setrand_eigen_compact(const Coordinates & coords,const Values & vals)241 EIGEN_DONT_INLINE Scalar* setrand_eigen_compact(const Coordinates& coords, const Values& vals)
242 {
243 using namespace Eigen;
244 DynamicSparseMatrix<Scalar> setter(SIZE,SIZE);
245 setter.reserve(coords.size()/10);
246 for (int i=0; i<coords.size(); ++i)
247 {
248 setter.coeffRef(coords[i].x(), coords[i].y()) += vals[i];
249 }
250 SparseMatrix<Scalar> mat = setter;
251 CHECK_MEM;
252 return &mat.coeffRef(coords[0].x(), coords[0].y());
253 }
254
setrand_eigen_gnu_hash(const Coordinates & coords,const Values & vals)255 EIGEN_DONT_INLINE Scalar* setrand_eigen_gnu_hash(const Coordinates& coords, const Values& vals)
256 {
257 using namespace Eigen;
258 SparseMatrix<Scalar> mat(SIZE,SIZE);
259 {
260 RandomSetter<SparseMatrix<Scalar>, StdMapTraits > setter(mat);
261 for (int i=0; i<coords.size(); ++i)
262 {
263 setter(coords[i].x(), coords[i].y()) += vals[i];
264 }
265 CHECK_MEM;
266 }
267 return &mat.coeffRef(coords[0].x(), coords[0].y());
268 }
269
270 #ifndef NOGOOGLE
setrand_eigen_google_dense(const Coordinates & coords,const Values & vals)271 EIGEN_DONT_INLINE Scalar* setrand_eigen_google_dense(const Coordinates& coords, const Values& vals)
272 {
273 using namespace Eigen;
274 SparseMatrix<Scalar> mat(SIZE,SIZE);
275 {
276 RandomSetter<SparseMatrix<Scalar>, GoogleDenseHashMapTraits> setter(mat);
277 for (int i=0; i<coords.size(); ++i)
278 setter(coords[i].x(), coords[i].y()) += vals[i];
279 CHECK_MEM;
280 }
281 return &mat.coeffRef(coords[0].x(), coords[0].y());
282 }
283
setrand_eigen_google_sparse(const Coordinates & coords,const Values & vals)284 EIGEN_DONT_INLINE Scalar* setrand_eigen_google_sparse(const Coordinates& coords, const Values& vals)
285 {
286 using namespace Eigen;
287 SparseMatrix<Scalar> mat(SIZE,SIZE);
288 {
289 RandomSetter<SparseMatrix<Scalar>, GoogleSparseHashMapTraits> setter(mat);
290 for (int i=0; i<coords.size(); ++i)
291 setter(coords[i].x(), coords[i].y()) += vals[i];
292 CHECK_MEM;
293 }
294 return &mat.coeffRef(coords[0].x(), coords[0].y());
295 }
296 #endif
297
298
299 template <class T>
coo_tocsr(const int n_row,const int n_col,const int nnz,const Coordinates Aij,const Values Ax,int Bp[],int Bj[],T Bx[])300 void coo_tocsr(const int n_row,
301 const int n_col,
302 const int nnz,
303 const Coordinates Aij,
304 const Values Ax,
305 int Bp[],
306 int Bj[],
307 T Bx[])
308 {
309 //compute number of non-zero entries per row of A coo_tocsr
310 std::fill(Bp, Bp + n_row, 0);
311
312 for (int n = 0; n < nnz; n++){
313 Bp[Aij[n].x()]++;
314 }
315
316 //cumsum the nnz per row to get Bp[]
317 for(int i = 0, cumsum = 0; i < n_row; i++){
318 int temp = Bp[i];
319 Bp[i] = cumsum;
320 cumsum += temp;
321 }
322 Bp[n_row] = nnz;
323
324 //write Aj,Ax into Bj,Bx
325 for(int n = 0; n < nnz; n++){
326 int row = Aij[n].x();
327 int dest = Bp[row];
328
329 Bj[dest] = Aij[n].y();
330 Bx[dest] = Ax[n];
331
332 Bp[row]++;
333 }
334
335 for(int i = 0, last = 0; i <= n_row; i++){
336 int temp = Bp[i];
337 Bp[i] = last;
338 last = temp;
339 }
340
341 //now Bp,Bj,Bx form a CSR representation (with possible duplicates)
342 }
343
344 template< class T1, class T2 >
kv_pair_less(const std::pair<T1,T2> & x,const std::pair<T1,T2> & y)345 bool kv_pair_less(const std::pair<T1,T2>& x, const std::pair<T1,T2>& y){
346 return x.first < y.first;
347 }
348
349
350 template<class I, class T>
csr_sort_indices(const I n_row,const I Ap[],I Aj[],T Ax[])351 void csr_sort_indices(const I n_row,
352 const I Ap[],
353 I Aj[],
354 T Ax[])
355 {
356 std::vector< std::pair<I,T> > temp;
357
358 for(I i = 0; i < n_row; i++){
359 I row_start = Ap[i];
360 I row_end = Ap[i+1];
361
362 temp.clear();
363
364 for(I jj = row_start; jj < row_end; jj++){
365 temp.push_back(std::make_pair(Aj[jj],Ax[jj]));
366 }
367
368 std::sort(temp.begin(),temp.end(),kv_pair_less<I,T>);
369
370 for(I jj = row_start, n = 0; jj < row_end; jj++, n++){
371 Aj[jj] = temp[n].first;
372 Ax[jj] = temp[n].second;
373 }
374 }
375 }
376
377 template <class I, class T>
csr_sum_duplicates(const I n_row,const I n_col,I Ap[],I Aj[],T Ax[])378 void csr_sum_duplicates(const I n_row,
379 const I n_col,
380 I Ap[],
381 I Aj[],
382 T Ax[])
383 {
384 I nnz = 0;
385 I row_end = 0;
386 for(I i = 0; i < n_row; i++){
387 I jj = row_end;
388 row_end = Ap[i+1];
389 while( jj < row_end ){
390 I j = Aj[jj];
391 T x = Ax[jj];
392 jj++;
393 while( jj < row_end && Aj[jj] == j ){
394 x += Ax[jj];
395 jj++;
396 }
397 Aj[nnz] = j;
398 Ax[nnz] = x;
399 nnz++;
400 }
401 Ap[i+1] = nnz;
402 }
403 }
404
setrand_scipy(const Coordinates & coords,const Values & vals)405 EIGEN_DONT_INLINE Scalar* setrand_scipy(const Coordinates& coords, const Values& vals)
406 {
407 using namespace Eigen;
408 SparseMatrix<Scalar> mat(SIZE,SIZE);
409 mat.resizeNonZeros(coords.size());
410 // std::cerr << "setrand_scipy...\n";
411 coo_tocsr<Scalar>(SIZE,SIZE, coords.size(), coords, vals, mat._outerIndexPtr(), mat._innerIndexPtr(), mat._valuePtr());
412 // std::cerr << "coo_tocsr ok\n";
413
414 csr_sort_indices(SIZE, mat._outerIndexPtr(), mat._innerIndexPtr(), mat._valuePtr());
415
416 csr_sum_duplicates(SIZE, SIZE, mat._outerIndexPtr(), mat._innerIndexPtr(), mat._valuePtr());
417
418 mat.resizeNonZeros(mat._outerIndexPtr()[SIZE]);
419
420 return &mat.coeffRef(coords[0].x(), coords[0].y());
421 }
422
423
424 #ifndef NOUBLAS
setrand_ublas_mapped(const Coordinates & coords,const Values & vals)425 EIGEN_DONT_INLINE Scalar* setrand_ublas_mapped(const Coordinates& coords, const Values& vals)
426 {
427 using namespace boost;
428 using namespace boost::numeric;
429 using namespace boost::numeric::ublas;
430 mapped_matrix<Scalar> aux(SIZE,SIZE);
431 for (int i=0; i<coords.size(); ++i)
432 {
433 aux(coords[i].x(), coords[i].y()) += vals[i];
434 }
435 CHECK_MEM;
436 compressed_matrix<Scalar> mat(aux);
437 return 0;// &mat(coords[0].x(), coords[0].y());
438 }
439 /*EIGEN_DONT_INLINE Scalar* setrand_ublas_coord(const Coordinates& coords, const Values& vals)
440 {
441 using namespace boost;
442 using namespace boost::numeric;
443 using namespace boost::numeric::ublas;
444 coordinate_matrix<Scalar> aux(SIZE,SIZE);
445 for (int i=0; i<coords.size(); ++i)
446 {
447 aux(coords[i].x(), coords[i].y()) = vals[i];
448 }
449 compressed_matrix<Scalar> mat(aux);
450 return 0;//&mat(coords[0].x(), coords[0].y());
451 }
452 EIGEN_DONT_INLINE Scalar* setrand_ublas_compressed(const Coordinates& coords, const Values& vals)
453 {
454 using namespace boost;
455 using namespace boost::numeric;
456 using namespace boost::numeric::ublas;
457 compressed_matrix<Scalar> mat(SIZE,SIZE);
458 for (int i=0; i<coords.size(); ++i)
459 {
460 mat(coords[i].x(), coords[i].y()) = vals[i];
461 }
462 return 0;//&mat(coords[0].x(), coords[0].y());
463 }*/
setrand_ublas_genvec(const Coordinates & coords,const Values & vals)464 EIGEN_DONT_INLINE Scalar* setrand_ublas_genvec(const Coordinates& coords, const Values& vals)
465 {
466 using namespace boost;
467 using namespace boost::numeric;
468 using namespace boost::numeric::ublas;
469
470 // ublas::vector<coordinate_vector<Scalar> > foo;
471 generalized_vector_of_vector<Scalar, row_major, ublas::vector<coordinate_vector<Scalar> > > aux(SIZE,SIZE);
472 for (int i=0; i<coords.size(); ++i)
473 {
474 aux(coords[i].x(), coords[i].y()) += vals[i];
475 }
476 CHECK_MEM;
477 compressed_matrix<Scalar,row_major> mat(aux);
478 return 0;//&mat(coords[0].x(), coords[0].y());
479 }
480 #endif
481
482 #ifndef NOMTL
483 EIGEN_DONT_INLINE void setrand_mtl(const Coordinates& coords, const Values& vals);
484 #endif
485
486