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
486