ground/gcs/src/libs/eigen/bench/sparse_setter.cpp

   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
  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
 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
 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
 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
 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
 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
 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
 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>
 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
 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
 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 }*/
 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