sparse_setter.cpp
Go to the documentation of this file.
00001 
00002 //g++ -O3 -g0 -DNDEBUG  sparse_product.cpp -I.. -I/home/gael/Coding/LinearAlgebra/mtl4/ -DDENSITY=0.005 -DSIZE=10000 && ./a.out
00003 //g++ -O3 -g0 -DNDEBUG  sparse_product.cpp -I.. -I/home/gael/Coding/LinearAlgebra/mtl4/ -DDENSITY=0.05 -DSIZE=2000 && ./a.out
00004 // -DNOGMM -DNOMTL -DCSPARSE
00005 // -I /home/gael/Coding/LinearAlgebra/CSparse/Include/ /home/gael/Coding/LinearAlgebra/CSparse/Lib/libcsparse.a
00006 #ifndef SIZE
00007 #define SIZE 100000
00008 #endif
00009 
00010 #ifndef NBPERROW
00011 #define NBPERROW 24
00012 #endif
00013 
00014 #ifndef REPEAT
00015 #define REPEAT 2
00016 #endif
00017 
00018 #ifndef NBTRIES
00019 #define NBTRIES 2
00020 #endif
00021 
00022 #ifndef KK
00023 #define KK 10
00024 #endif
00025 
00026 #ifndef NOGOOGLE
00027 #define EIGEN_GOOGLEHASH_SUPPORT
00028 #include <google/sparse_hash_map>
00029 #endif
00030 
00031 #include "BenchSparseUtil.h"
00032 
00033 #define CHECK_MEM
00034 // #define CHECK_MEM  std/**/::cout << "check mem\n"; getchar();
00035 
00036 #define BENCH(X) \
00037   timer.reset(); \
00038   for (int _j=0; _j<NBTRIES; ++_j) { \
00039     timer.start(); \
00040     for (int _k=0; _k<REPEAT; ++_k) { \
00041         X  \
00042   } timer.stop(); }
00043 
00044 typedef std::vector<Vector2i> Coordinates;
00045 typedef std::vector<float> Values;
00046 
00047 EIGEN_DONT_INLINE Scalar* setinnerrand_eigen(const Coordinates& coords, const Values& vals);
00048 EIGEN_DONT_INLINE Scalar* setrand_eigen_dynamic(const Coordinates& coords, const Values& vals);
00049 EIGEN_DONT_INLINE Scalar* setrand_eigen_compact(const Coordinates& coords, const Values& vals);
00050 EIGEN_DONT_INLINE Scalar* setrand_eigen_sumeq(const Coordinates& coords, const Values& vals);
00051 EIGEN_DONT_INLINE Scalar* setrand_eigen_gnu_hash(const Coordinates& coords, const Values& vals);
00052 EIGEN_DONT_INLINE Scalar* setrand_eigen_google_dense(const Coordinates& coords, const Values& vals);
00053 EIGEN_DONT_INLINE Scalar* setrand_eigen_google_sparse(const Coordinates& coords, const Values& vals);
00054 EIGEN_DONT_INLINE Scalar* setrand_scipy(const Coordinates& coords, const Values& vals);
00055 EIGEN_DONT_INLINE Scalar* setrand_ublas_mapped(const Coordinates& coords, const Values& vals);
00056 EIGEN_DONT_INLINE Scalar* setrand_ublas_coord(const Coordinates& coords, const Values& vals);
00057 EIGEN_DONT_INLINE Scalar* setrand_ublas_compressed(const Coordinates& coords, const Values& vals);
00058 EIGEN_DONT_INLINE Scalar* setrand_ublas_genvec(const Coordinates& coords, const Values& vals);
00059 EIGEN_DONT_INLINE Scalar* setrand_mtl(const Coordinates& coords, const Values& vals);
00060 
00061 int main(int argc, char *argv[])
00062 {
00063   int rows = SIZE;
00064   int cols = SIZE;
00065   bool fullyrand = true;
00066 
00067   BenchTimer timer;
00068   Coordinates coords;
00069   Values values;
00070   if(fullyrand)
00071   {
00072     Coordinates pool;
00073     pool.reserve(cols*NBPERROW);
00074     std::cerr << "fill pool" << "\n";
00075     for (int i=0; i<cols*NBPERROW; )
00076     {
00077 //       DynamicSparseMatrix<int> stencil(SIZE,SIZE);
00078       Vector2i ij(internal::random<int>(0,rows-1),internal::random<int>(0,cols-1));
00079 //       if(stencil.coeffRef(ij.x(), ij.y())==0)
00080       {
00081 //         stencil.coeffRef(ij.x(), ij.y()) = 1;
00082         pool.push_back(ij);
00083 
00084       }
00085       ++i;
00086     }
00087     std::cerr << "pool ok" << "\n";
00088     int n = cols*NBPERROW*KK;
00089     coords.reserve(n);
00090     values.reserve(n);
00091     for (int i=0; i<n; ++i)
00092     {
00093       int i = internal::random<int>(0,pool.size());
00094       coords.push_back(pool[i]);
00095       values.push_back(internal::random<Scalar>());
00096     }
00097   }
00098   else
00099   {
00100     for (int j=0; j<cols; ++j)
00101     for (int i=0; i<NBPERROW; ++i)
00102     {
00103       coords.push_back(Vector2i(internal::random<int>(0,rows-1),j));
00104       values.push_back(internal::random<Scalar>());
00105     }
00106   }
00107   std::cout << "nnz = " << coords.size()  << "\n";
00108   CHECK_MEM
00109 
00110     // dense matrices
00111     #ifdef DENSEMATRIX
00112     {
00113       BENCH(setrand_eigen_dense(coords,values);)
00114       std::cout << "Eigen Dense\t" << timer.value() << "\n";
00115     }
00116     #endif
00117 
00118     // eigen sparse matrices
00119 //     if (!fullyrand)
00120 //     {
00121 //       BENCH(setinnerrand_eigen(coords,values);)
00122 //       std::cout << "Eigen fillrand\t" << timer.value() << "\n";
00123 //     }
00124     {
00125       BENCH(setrand_eigen_dynamic(coords,values);)
00126       std::cout << "Eigen dynamic\t" << timer.value() << "\n";
00127     }
00128 //     {
00129 //       BENCH(setrand_eigen_compact(coords,values);)
00130 //       std::cout << "Eigen compact\t" << timer.value() << "\n";
00131 //     }
00132     {
00133       BENCH(setrand_eigen_sumeq(coords,values);)
00134       std::cout << "Eigen sumeq\t" << timer.value() << "\n";
00135     }
00136     {
00137 //       BENCH(setrand_eigen_gnu_hash(coords,values);)
00138 //       std::cout << "Eigen std::map\t" << timer.value() << "\n";
00139     }
00140     {
00141       BENCH(setrand_scipy(coords,values);)
00142       std::cout << "scipy\t" << timer.value() << "\n";
00143     }
00144     #ifndef NOGOOGLE
00145     {
00146       BENCH(setrand_eigen_google_dense(coords,values);)
00147       std::cout << "Eigen google dense\t" << timer.value() << "\n";
00148     }
00149     {
00150       BENCH(setrand_eigen_google_sparse(coords,values);)
00151       std::cout << "Eigen google sparse\t" << timer.value() << "\n";
00152     }
00153     #endif
00154 
00155     #ifndef NOUBLAS
00156     {
00157 //       BENCH(setrand_ublas_mapped(coords,values);)
00158 //       std::cout << "ublas mapped\t" << timer.value() << "\n";
00159     }
00160     {
00161       BENCH(setrand_ublas_genvec(coords,values);)
00162       std::cout << "ublas vecofvec\t" << timer.value() << "\n";
00163     }
00164     /*{
00165       timer.reset();
00166       timer.start();
00167       for (int k=0; k<REPEAT; ++k)
00168         setrand_ublas_compressed(coords,values);
00169       timer.stop();
00170       std::cout << "ublas comp\t" << timer.value() << "\n";
00171     }
00172     {
00173       timer.reset();
00174       timer.start();
00175       for (int k=0; k<REPEAT; ++k)
00176         setrand_ublas_coord(coords,values);
00177       timer.stop();
00178       std::cout << "ublas coord\t" << timer.value() << "\n";
00179     }*/
00180     #endif
00181 
00182 
00183     // MTL4
00184     #ifndef NOMTL
00185     {
00186       BENCH(setrand_mtl(coords,values));
00187       std::cout << "MTL\t" << timer.value() << "\n";
00188     }
00189     #endif
00190 
00191   return 0;
00192 }
00193 
00194 EIGEN_DONT_INLINE Scalar* setinnerrand_eigen(const Coordinates& coords, const Values& vals)
00195 {
00196   using namespace Eigen;
00197   SparseMatrix<Scalar> mat(SIZE,SIZE);
00198   //mat.startFill(2000000/*coords.size()*/);
00199   for (int i=0; i<coords.size(); ++i)
00200   {
00201     mat.insert(coords[i].x(), coords[i].y()) = vals[i];
00202   }
00203   mat.finalize();
00204   CHECK_MEM;
00205   return 0;
00206 }
00207 
00208 EIGEN_DONT_INLINE Scalar* setrand_eigen_dynamic(const Coordinates& coords, const Values& vals)
00209 {
00210   using namespace Eigen;
00211   DynamicSparseMatrix<Scalar> mat(SIZE,SIZE);
00212   mat.reserve(coords.size()/10);
00213   for (int i=0; i<coords.size(); ++i)
00214   {
00215     mat.coeffRef(coords[i].x(), coords[i].y()) += vals[i];
00216   }
00217   mat.finalize();
00218   CHECK_MEM;
00219   return &mat.coeffRef(coords[0].x(), coords[0].y());
00220 }
00221 
00222 EIGEN_DONT_INLINE Scalar* setrand_eigen_sumeq(const Coordinates& coords, const Values& vals)
00223 {
00224   using namespace Eigen;
00225   int n = coords.size()/KK;
00226   DynamicSparseMatrix<Scalar> mat(SIZE,SIZE);
00227   for (int j=0; j<KK; ++j)
00228   {
00229     DynamicSparseMatrix<Scalar> aux(SIZE,SIZE);
00230     mat.reserve(n);
00231     for (int i=j*n; i<(j+1)*n; ++i)
00232     {
00233       aux.insert(coords[i].x(), coords[i].y()) += vals[i];
00234     }
00235     aux.finalize();
00236     mat += aux;
00237   }
00238   return &mat.coeffRef(coords[0].x(), coords[0].y());
00239 }
00240 
00241 EIGEN_DONT_INLINE Scalar* setrand_eigen_compact(const Coordinates& coords, const Values& vals)
00242 {
00243   using namespace Eigen;
00244   DynamicSparseMatrix<Scalar> setter(SIZE,SIZE);
00245   setter.reserve(coords.size()/10);
00246   for (int i=0; i<coords.size(); ++i)
00247   {
00248     setter.coeffRef(coords[i].x(), coords[i].y()) += vals[i];
00249   }
00250   SparseMatrix<Scalar> mat = setter;
00251   CHECK_MEM;
00252   return &mat.coeffRef(coords[0].x(), coords[0].y());
00253 }
00254 
00255 EIGEN_DONT_INLINE Scalar* setrand_eigen_gnu_hash(const Coordinates& coords, const Values& vals)
00256 {
00257   using namespace Eigen;
00258   SparseMatrix<Scalar> mat(SIZE,SIZE);
00259   {
00260     RandomSetter<SparseMatrix<Scalar>, StdMapTraits > setter(mat);
00261     for (int i=0; i<coords.size(); ++i)
00262     {
00263       setter(coords[i].x(), coords[i].y()) += vals[i];
00264     }
00265     CHECK_MEM;
00266   }
00267   return &mat.coeffRef(coords[0].x(), coords[0].y());
00268 }
00269 
00270 #ifndef NOGOOGLE
00271 EIGEN_DONT_INLINE Scalar* setrand_eigen_google_dense(const Coordinates& coords, const Values& vals)
00272 {
00273   using namespace Eigen;
00274   SparseMatrix<Scalar> mat(SIZE,SIZE);
00275   {
00276     RandomSetter<SparseMatrix<Scalar>, GoogleDenseHashMapTraits> setter(mat);
00277     for (int i=0; i<coords.size(); ++i)
00278       setter(coords[i].x(), coords[i].y()) += vals[i];
00279     CHECK_MEM;
00280   }
00281   return &mat.coeffRef(coords[0].x(), coords[0].y());
00282 }
00283 
00284 EIGEN_DONT_INLINE Scalar* setrand_eigen_google_sparse(const Coordinates& coords, const Values& vals)
00285 {
00286   using namespace Eigen;
00287   SparseMatrix<Scalar> mat(SIZE,SIZE);
00288   {
00289     RandomSetter<SparseMatrix<Scalar>, GoogleSparseHashMapTraits> setter(mat);
00290     for (int i=0; i<coords.size(); ++i)
00291       setter(coords[i].x(), coords[i].y()) += vals[i];
00292     CHECK_MEM;
00293   }
00294   return &mat.coeffRef(coords[0].x(), coords[0].y());
00295 }
00296 #endif
00297 
00298 
00299 template <class T>
00300 void coo_tocsr(const int n_row,
00301                const int n_col,
00302                const int nnz,
00303                const Coordinates Aij,
00304                const Values Ax,
00305                      int Bp[],
00306                      int Bj[],
00307                      T Bx[])
00308 {
00309     //compute number of non-zero entries per row of A coo_tocsr
00310     std::fill(Bp, Bp + n_row, 0);
00311 
00312     for (int n = 0; n < nnz; n++){
00313         Bp[Aij[n].x()]++;
00314     }
00315 
00316     //cumsum the nnz per row to get Bp[]
00317     for(int i = 0, cumsum = 0; i < n_row; i++){
00318         int temp = Bp[i];
00319         Bp[i] = cumsum;
00320         cumsum += temp;
00321     }
00322     Bp[n_row] = nnz;
00323 
00324     //write Aj,Ax into Bj,Bx
00325     for(int n = 0; n < nnz; n++){
00326         int row  = Aij[n].x();
00327         int dest = Bp[row];
00328 
00329         Bj[dest] = Aij[n].y();
00330         Bx[dest] = Ax[n];
00331 
00332         Bp[row]++;
00333     }
00334 
00335     for(int i = 0, last = 0; i <= n_row; i++){
00336         int temp = Bp[i];
00337         Bp[i]  = last;
00338         last   = temp;
00339     }
00340 
00341     //now Bp,Bj,Bx form a CSR representation (with possible duplicates)
00342 }
00343 
00344 template< class T1, class T2 >
00345 bool kv_pair_less(const std::pair<T1,T2>& x, const std::pair<T1,T2>& y){
00346     return x.first < y.first;
00347 }
00348 
00349 
00350 template<class I, class T>
00351 void csr_sort_indices(const I n_row,
00352                       const I Ap[],
00353                             I Aj[],
00354                             T Ax[])
00355 {
00356     std::vector< std::pair<I,T> > temp;
00357 
00358     for(I i = 0; i < n_row; i++){
00359         I row_start = Ap[i];
00360         I row_end   = Ap[i+1];
00361 
00362         temp.clear();
00363 
00364         for(I jj = row_start; jj < row_end; jj++){
00365             temp.push_back(std::make_pair(Aj[jj],Ax[jj]));
00366         }
00367 
00368         std::sort(temp.begin(),temp.end(),kv_pair_less<I,T>);
00369 
00370         for(I jj = row_start, n = 0; jj < row_end; jj++, n++){
00371             Aj[jj] = temp[n].first;
00372             Ax[jj] = temp[n].second;
00373         }
00374     }
00375 }
00376 
00377 template <class I, class T>
00378 void csr_sum_duplicates(const I n_row,
00379                         const I n_col,
00380                               I Ap[],
00381                               I Aj[],
00382                               T Ax[])
00383 {
00384     I nnz = 0;
00385     I row_end = 0;
00386     for(I i = 0; i < n_row; i++){
00387         I jj = row_end;
00388         row_end = Ap[i+1];
00389         while( jj < row_end ){
00390             I j = Aj[jj];
00391             T x = Ax[jj];
00392             jj++;
00393             while( jj < row_end && Aj[jj] == j ){
00394                 x += Ax[jj];
00395                 jj++;
00396             }
00397             Aj[nnz] = j;
00398             Ax[nnz] = x;
00399             nnz++;
00400         }
00401         Ap[i+1] = nnz;
00402     }
00403 }
00404 
00405 EIGEN_DONT_INLINE Scalar* setrand_scipy(const Coordinates& coords, const Values& vals)
00406 {
00407   using namespace Eigen;
00408   SparseMatrix<Scalar> mat(SIZE,SIZE);
00409   mat.resizeNonZeros(coords.size());
00410 //   std::cerr << "setrand_scipy...\n";
00411   coo_tocsr<Scalar>(SIZE,SIZE, coords.size(), coords, vals, mat._outerIndexPtr(), mat._innerIndexPtr(), mat._valuePtr());
00412 //   std::cerr << "coo_tocsr ok\n";
00413 
00414   csr_sort_indices(SIZE, mat._outerIndexPtr(), mat._innerIndexPtr(), mat._valuePtr());
00415 
00416   csr_sum_duplicates(SIZE, SIZE, mat._outerIndexPtr(), mat._innerIndexPtr(), mat._valuePtr());
00417 
00418   mat.resizeNonZeros(mat._outerIndexPtr()[SIZE]);
00419 
00420   return &mat.coeffRef(coords[0].x(), coords[0].y());
00421 }
00422 
00423 
00424 #ifndef NOUBLAS
00425 EIGEN_DONT_INLINE Scalar* setrand_ublas_mapped(const Coordinates& coords, const Values& vals)
00426 {
00427   using namespace boost;
00428   using namespace boost::numeric;
00429   using namespace boost::numeric::ublas;
00430   mapped_matrix<Scalar> aux(SIZE,SIZE);
00431   for (int i=0; i<coords.size(); ++i)
00432   {
00433     aux(coords[i].x(), coords[i].y()) += vals[i];
00434   }
00435   CHECK_MEM;
00436   compressed_matrix<Scalar> mat(aux);
00437   return 0;// &mat(coords[0].x(), coords[0].y());
00438 }
00439 /*EIGEN_DONT_INLINE Scalar* setrand_ublas_coord(const Coordinates& coords, const Values& vals)
00440 {
00441   using namespace boost;
00442   using namespace boost::numeric;
00443   using namespace boost::numeric::ublas;
00444   coordinate_matrix<Scalar> aux(SIZE,SIZE);
00445   for (int i=0; i<coords.size(); ++i)
00446   {
00447     aux(coords[i].x(), coords[i].y()) = vals[i];
00448   }
00449   compressed_matrix<Scalar> mat(aux);
00450   return 0;//&mat(coords[0].x(), coords[0].y());
00451 }
00452 EIGEN_DONT_INLINE Scalar* setrand_ublas_compressed(const Coordinates& coords, const Values& vals)
00453 {
00454   using namespace boost;
00455   using namespace boost::numeric;
00456   using namespace boost::numeric::ublas;
00457   compressed_matrix<Scalar> mat(SIZE,SIZE);
00458   for (int i=0; i<coords.size(); ++i)
00459   {
00460     mat(coords[i].x(), coords[i].y()) = vals[i];
00461   }
00462   return 0;//&mat(coords[0].x(), coords[0].y());
00463 }*/
00464 EIGEN_DONT_INLINE Scalar* setrand_ublas_genvec(const Coordinates& coords, const Values& vals)
00465 {
00466   using namespace boost;
00467   using namespace boost::numeric;
00468   using namespace boost::numeric::ublas;
00469 
00470 //   ublas::vector<coordinate_vector<Scalar> > foo;
00471   generalized_vector_of_vector<Scalar, row_major, ublas::vector<coordinate_vector<Scalar> > > aux(SIZE,SIZE);
00472   for (int i=0; i<coords.size(); ++i)
00473   {
00474     aux(coords[i].x(), coords[i].y()) += vals[i];
00475   }
00476   CHECK_MEM;
00477   compressed_matrix<Scalar,row_major> mat(aux);
00478   return 0;//&mat(coords[0].x(), coords[0].y());
00479 }
00480 #endif
00481 
00482 #ifndef NOMTL
00483 EIGEN_DONT_INLINE void setrand_mtl(const Coordinates& coords, const Values& vals);
00484 #endif
00485 


libicr
Author(s): Robert Krug
autogenerated on Mon Jan 6 2014 11:33:32