source: pacpussensors/trunk/Vislab/lib3dv/eigen/bench/sparse_setter.cpp@ 141

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

Doc
File size: 13.4 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#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
44typedef std::vector<Vector2i> Coordinates;
45typedef std::vector<float> Values;
46
47EIGEN_DONT_INLINE Scalar* setinnerrand_eigen(const Coordinates& coords, const Values& vals);
48EIGEN_DONT_INLINE Scalar* setrand_eigen_dynamic(const Coordinates& coords, const Values& vals);
49EIGEN_DONT_INLINE Scalar* setrand_eigen_compact(const Coordinates& coords, const Values& vals);
50EIGEN_DONT_INLINE Scalar* setrand_eigen_sumeq(const Coordinates& coords, const Values& vals);
51EIGEN_DONT_INLINE Scalar* setrand_eigen_gnu_hash(const Coordinates& coords, const Values& vals);
52EIGEN_DONT_INLINE Scalar* setrand_eigen_google_dense(const Coordinates& coords, const Values& vals);
53EIGEN_DONT_INLINE Scalar* setrand_eigen_google_sparse(const Coordinates& coords, const Values& vals);
54EIGEN_DONT_INLINE Scalar* setrand_scipy(const Coordinates& coords, const Values& vals);
55EIGEN_DONT_INLINE Scalar* setrand_ublas_mapped(const Coordinates& coords, const Values& vals);
56EIGEN_DONT_INLINE Scalar* setrand_ublas_coord(const Coordinates& coords, const Values& vals);
57EIGEN_DONT_INLINE Scalar* setrand_ublas_compressed(const Coordinates& coords, const Values& vals);
58EIGEN_DONT_INLINE Scalar* setrand_ublas_genvec(const Coordinates& coords, const Values& vals);
59EIGEN_DONT_INLINE Scalar* setrand_mtl(const Coordinates& coords, const Values& vals);
60
61int 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
194EIGEN_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
208EIGEN_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
222EIGEN_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
241EIGEN_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
255EIGEN_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
271EIGEN_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
284EIGEN_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
299template <class T>
300void 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
344template< class T1, class T2 >
345bool kv_pair_less(const std::pair<T1,T2>& x, const std::pair<T1,T2>& y){
346 return x.first < y.first;
347}
348
349
350template<class I, class T>
351void 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
377template <class I, class T>
378void 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
405EIGEN_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
425EIGEN_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}
452EIGEN_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}*/
464EIGEN_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
483EIGEN_DONT_INLINE void setrand_mtl(const Coordinates& coords, const Values& vals);
484#endif
485
Note: See TracBrowser for help on using the repository browser.