sparse_setter.cpp
Go to the documentation of this file.
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);
54 EIGEN_DONT_INLINE Scalar* setrand_scipy(const Coordinates& coords, const Values& vals);
59 EIGEN_DONT_INLINE Scalar* setrand_mtl(const Coordinates& coords, const Values& vals);
60 
61 int main(int argc, char *argv[])
62 {
63  int rows = SIZE;
64  int cols = SIZE;
65  bool fullyrand = true;
66 
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 
195 {
196  using namespace Eigen;
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 
209 {
210  using namespace Eigen;
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 
223 {
224  using namespace Eigen;
225  int n = coords.size()/KK;
227  for (int j=0; j<KK; ++j)
228  {
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 
242 {
243  using namespace Eigen;
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 
256 {
257  using namespace Eigen;
259  {
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
272 {
273  using namespace Eigen;
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 
285 {
286  using namespace Eigen;
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>
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 >
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>
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>
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 
406 {
407  using namespace Eigen;
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
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 }*/
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 
SCALAR Scalar
Definition: bench_gemm.cpp:46
#define I
Definition: main.h:112
EIGEN_DONT_INLINE Scalar * setrand_ublas_compressed(const Coordinates &coords, const Values &vals)
EIGEN_DONT_INLINE Scalar * setrand_eigen_compact(const Coordinates &coords, const Values &vals)
Scalar * y
EIGEN_DONT_INLINE Scalar * setinnerrand_eigen(const Coordinates &coords, const Values &vals)
EIGEN_DONT_INLINE Scalar * setrand_ublas_genvec(const Coordinates &coords, const Values &vals)
#define SIZE
void 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[])
bool kv_pair_less(const std::pair< T1, T2 > &x, const std::pair< T1, T2 > &y)
int n
leaf::MyValues values
EIGEN_DONT_INLINE Scalar * setrand_ublas_mapped(const Coordinates &coords, const Values &vals)
double value(int TIMER=CPU_TIMER) const
Definition: BenchTimer.h:104
Namespace containing all symbols from the Eigen library.
Definition: jet.h:637
EIGEN_DONT_INLINE Scalar * setrand_eigen_google_dense(const Coordinates &coords, const Values &vals)
EIGEN_DONT_INLINE Scalar * setrand_eigen_google_sparse(const Coordinates &coords, const Values &vals)
void resizeNonZeros(Index size)
Definition: SparseMatrix.h:649
static const symbolic::SymbolExpr< internal::symbolic_last_tag > last
EIGEN_DONT_INLINE Scalar * setrand_eigen_sumeq(const Coordinates &coords, const Values &vals)
#define EIGEN_DONT_INLINE
Definition: Macros.h:940
EIGEN_DONT_INLINE Scalar * setrand_ublas_coord(const Coordinates &coords, const Values &vals)
m row(1)
EIGEN_DONT_INLINE Scalar * setrand_eigen_gnu_hash(const Coordinates &coords, const Values &vals)
Scalar & insert(Index row, Index col)
EIGEN_DONT_INLINE Scalar * setrand_mtl(const Coordinates &coords, const Values &vals)
A sparse matrix class designed for matrix assembly purpose.
Definition: SparseUtil.h:53
void reserve(Index reserveSize=1000)
int main(int argc, char *argv[])
Scalar & coeffRef(Index row, Index col)
Definition: SparseMatrix.h:208
std::vector< float > Values
void csr_sort_indices(const I n_row, const I Ap[], I Aj[], T Ax[])
std::vector< Vector2i > Coordinates
#define CHECK_MEM
The RandomSetter is a wrapper object allowing to set/update a sparse matrix with random access...
Definition: RandomSetter.h:176
EIGEN_DONT_INLINE Scalar * setrand_scipy(const Coordinates &coords, const Values &vals)
Scalar & coeffRef(Index row, Index col)
Scalar & insert(Index row, Index col)
set noclip points set clip one set noclip two set bar set border lt lw set xdata set ydata set zdata set x2data set y2data set boxwidth set dummy x
#define BENCH(X)
static BenchTimer timer
std::ptrdiff_t j
void csr_sum_duplicates(const I n_row, const I n_col, I Ap[], I Aj[], T Ax[])
EIGEN_DONT_INLINE Scalar * setrand_eigen_dynamic(const Coordinates &coords, const Values &vals)
#define KK
#define NBPERROW


gtsam
Author(s):
autogenerated on Tue Jul 4 2023 02:35:54