sparse_product.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 
00007 #include <typeinfo>
00008 
00009 #ifndef SIZE
00010 #define SIZE 1000000
00011 #endif
00012 
00013 #ifndef NNZPERCOL
00014 #define NNZPERCOL 6
00015 #endif
00016 
00017 #ifndef REPEAT
00018 #define REPEAT 1
00019 #endif
00020 
00021 #include <algorithm>
00022 #include "BenchTimer.h"
00023 #include "BenchSparseUtil.h"
00024 
00025 #ifndef NBTRIES
00026 #define NBTRIES 1
00027 #endif
00028 
00029 #define BENCH(X) \
00030   timer.reset(); \
00031   for (int _j=0; _j<NBTRIES; ++_j) { \
00032     timer.start(); \
00033     for (int _k=0; _k<REPEAT; ++_k) { \
00034         X  \
00035   } timer.stop(); }
00036 
00037 // #ifdef MKL
00038 //
00039 // #include "mkl_types.h"
00040 // #include "mkl_spblas.h"
00041 //
00042 // template<typename Lhs,typename Rhs,typename Res>
00043 // void mkl_multiply(const Lhs& lhs, const Rhs& rhs, Res& res)
00044 // {
00045 //   char n = 'N';
00046 //   float alpha = 1;
00047 //   char matdescra[6];
00048 //   matdescra[0] = 'G';
00049 //   matdescra[1] = 0;
00050 //   matdescra[2] = 0;
00051 //   matdescra[3] = 'C';
00052 //   mkl_scscmm(&n, lhs.rows(), rhs.cols(), lhs.cols(), &alpha, matdescra,
00053 //              lhs._valuePtr(), lhs._innerIndexPtr(), lhs.outerIndexPtr(),
00054 //              pntre, b, &ldb, &beta, c, &ldc);
00055 // //   mkl_somatcopy('C', 'T', lhs.rows(), lhs.cols(), 1,
00056 // //                 lhs._valuePtr(), lhs.rows(), DST, dst_stride);
00057 // }
00058 //
00059 // #endif
00060 
00061 
00062 #ifdef CSPARSE
00063 cs* cs_sorted_multiply(const cs* a, const cs* b)
00064 {
00065 //   return cs_multiply(a,b);
00066 
00067   cs* A = cs_transpose(a, 1);
00068   cs* B = cs_transpose(b, 1);
00069   cs* D = cs_multiply(B,A);   /* D = B'*A' */
00070   cs_spfree (A) ;
00071   cs_spfree (B) ;
00072   cs_dropzeros (D) ;      /* drop zeros from D */
00073   cs* C = cs_transpose (D, 1) ;   /* C = D', so that C is sorted */
00074   cs_spfree (D) ;
00075   return C;
00076 
00077 //   cs* A = cs_transpose(a, 1);
00078 //   cs* C = cs_transpose(A, 1);
00079 //   return C;
00080 }
00081 
00082 cs* cs_sorted_multiply2(const cs* a, const cs* b)
00083 {
00084   cs* D = cs_multiply(a,b);
00085   cs* E = cs_transpose(D,1);
00086   cs_spfree(D);
00087   cs* C = cs_transpose(E,1);
00088   cs_spfree(E);
00089   return C;
00090 }
00091 #endif
00092 
00093 void bench_sort();
00094 
00095 int main(int argc, char *argv[])
00096 {
00097 //   bench_sort();
00098 
00099   int rows = SIZE;
00100   int cols = SIZE;
00101   float density = DENSITY;
00102 
00103   EigenSparseMatrix sm1(rows,cols), sm2(rows,cols), sm3(rows,cols), sm4(rows,cols);
00104 
00105   BenchTimer timer;
00106   for (int nnzPerCol = NNZPERCOL; nnzPerCol>1; nnzPerCol/=1.1)
00107   {
00108     sm1.setZero();
00109     sm2.setZero();
00110     fillMatrix2(nnzPerCol, rows, cols, sm1);
00111     fillMatrix2(nnzPerCol, rows, cols, sm2);
00112 //     std::cerr << "filling OK\n";
00113 
00114     // dense matrices
00115     #ifdef DENSEMATRIX
00116     {
00117       std::cout << "Eigen Dense\t" << nnzPerCol << "%\n";
00118       DenseMatrix m1(rows,cols), m2(rows,cols), m3(rows,cols);
00119       eiToDense(sm1, m1);
00120       eiToDense(sm2, m2);
00121 
00122       timer.reset();
00123       timer.start();
00124       for (int k=0; k<REPEAT; ++k)
00125         m3 = m1 * m2;
00126       timer.stop();
00127       std::cout << "   a * b:\t" << timer.value() << endl;
00128 
00129       timer.reset();
00130       timer.start();
00131       for (int k=0; k<REPEAT; ++k)
00132         m3 = m1.transpose() * m2;
00133       timer.stop();
00134       std::cout << "   a' * b:\t" << timer.value() << endl;
00135 
00136       timer.reset();
00137       timer.start();
00138       for (int k=0; k<REPEAT; ++k)
00139         m3 = m1.transpose() * m2.transpose();
00140       timer.stop();
00141       std::cout << "   a' * b':\t" << timer.value() << endl;
00142 
00143       timer.reset();
00144       timer.start();
00145       for (int k=0; k<REPEAT; ++k)
00146         m3 = m1 * m2.transpose();
00147       timer.stop();
00148       std::cout << "   a * b':\t" << timer.value() << endl;
00149     }
00150     #endif
00151 
00152     // eigen sparse matrices
00153     {
00154       std::cout << "Eigen sparse\t" << sm1.nonZeros()/(float(sm1.rows())*float(sm1.cols()))*100 << "% * "
00155                 << sm2.nonZeros()/(float(sm2.rows())*float(sm2.cols()))*100 << "%\n";
00156 
00157       BENCH(sm3 = sm1 * sm2; )
00158       std::cout << "   a * b:\t" << timer.value() << endl;
00159 
00160 //       BENCH(sm3 = sm1.transpose() * sm2; )
00161 //       std::cout << "   a' * b:\t" << timer.value() << endl;
00162 // //
00163 //       BENCH(sm3 = sm1.transpose() * sm2.transpose(); )
00164 //       std::cout << "   a' * b':\t" << timer.value() << endl;
00165 // //
00166 //       BENCH(sm3 = sm1 * sm2.transpose(); )
00167 //       std::cout << "   a * b' :\t" << timer.value() << endl;
00168 
00169 
00170 //       std::cout << "\n";
00171 //
00172 //       BENCH( sm3._experimentalNewProduct(sm1, sm2); )
00173 //       std::cout << "   a * b:\t" << timer.value() << endl;
00174 //
00175 //       BENCH(sm3._experimentalNewProduct(sm1.transpose(),sm2); )
00176 //       std::cout << "   a' * b:\t" << timer.value() << endl;
00177 // //
00178 //       BENCH(sm3._experimentalNewProduct(sm1.transpose(),sm2.transpose()); )
00179 //       std::cout << "   a' * b':\t" << timer.value() << endl;
00180 // //
00181 //       BENCH(sm3._experimentalNewProduct(sm1, sm2.transpose());)
00182 //       std::cout << "   a * b' :\t" << timer.value() << endl;
00183     }
00184 
00185     // eigen dyn-sparse matrices
00186     /*{
00187       DynamicSparseMatrix<Scalar> m1(sm1), m2(sm2), m3(sm3);
00188       std::cout << "Eigen dyn-sparse\t" << m1.nonZeros()/(float(m1.rows())*float(m1.cols()))*100 << "% * "
00189                 << m2.nonZeros()/(float(m2.rows())*float(m2.cols()))*100 << "%\n";
00190 
00191 //       timer.reset();
00192 //       timer.start();
00193       BENCH(for (int k=0; k<REPEAT; ++k) m3 = m1 * m2;)
00194 //       timer.stop();
00195       std::cout << "   a * b:\t" << timer.value() << endl;
00196 //       std::cout << sm3 << "\n";
00197 
00198       timer.reset();
00199       timer.start();
00200 //       std::cerr << "transpose...\n";
00201 //       EigenSparseMatrix sm4 = sm1.transpose();
00202 //       std::cout << sm4.nonZeros() << " == " << sm1.nonZeros() << "\n";
00203 //       exit(1);
00204 //       std::cerr << "transpose OK\n";
00205 //       std::cout << sm1 << "\n\n" << sm1.transpose() << "\n\n" << sm4.transpose() << "\n\n";
00206       BENCH(for (int k=0; k<REPEAT; ++k) m3 = m1.transpose() * m2;)
00207 //       timer.stop();
00208       std::cout << "   a' * b:\t" << timer.value() << endl;
00209 
00210 //       timer.reset();
00211 //       timer.start();
00212       BENCH( for (int k=0; k<REPEAT; ++k) m3 = m1.transpose() * m2.transpose(); )
00213 //       timer.stop();
00214       std::cout << "   a' * b':\t" << timer.value() << endl;
00215 
00216 //       timer.reset();
00217 //       timer.start();
00218       BENCH( for (int k=0; k<REPEAT; ++k) m3 = m1 * m2.transpose(); )
00219 //       timer.stop();
00220       std::cout << "   a * b' :\t" << timer.value() << endl;
00221     }*/
00222 
00223     // CSparse
00224     #ifdef CSPARSE
00225     {
00226       std::cout << "CSparse \t" << nnzPerCol << "%\n";
00227       cs *m1, *m2, *m3;
00228       eiToCSparse(sm1, m1);
00229       eiToCSparse(sm2, m2);
00230 
00231 //       timer.reset();
00232 //       timer.start();
00233 //       for (int k=0; k<REPEAT; ++k)
00234       BENCH(
00235       {
00236         m3 = cs_sorted_multiply(m1, m2);
00237         if (!m3)
00238         {
00239           std::cerr << "cs_multiply failed\n";
00240 //           break;
00241         }
00242 //         cs_print(m3, 0);
00243         cs_spfree(m3);
00244       }
00245       );
00246 //       timer.stop();
00247       std::cout << "   a * b:\t" << timer.value() << endl;
00248 
00249 //       BENCH( { m3 = cs_sorted_multiply2(m1, m2); cs_spfree(m3); } );
00250 //       std::cout << "   a * b:\t" << timer.value() << endl;
00251     }
00252     #endif
00253 
00254     #ifndef NOUBLAS
00255     {
00256       std::cout << "ublas\t" << nnzPerCol << "%\n";
00257       UblasMatrix m1(rows,cols), m2(rows,cols), m3(rows,cols);
00258       eiToUblas(sm1, m1);
00259       eiToUblas(sm2, m2);
00260 
00261       BENCH(boost::numeric::ublas::prod(m1, m2, m3););
00262 //       timer.reset();
00263 //       timer.start();
00264 //       for (int k=0; k<REPEAT; ++k)
00265 //         gmm::mult(m1, m2, gmmT3);
00266 //       timer.stop();
00267       std::cout << "   a * b:\t" << timer.value() << endl;
00268     }
00269     #endif
00270 
00271     // GMM++
00272     #ifndef NOGMM
00273     {
00274       std::cout << "GMM++ sparse\t" << nnzPerCol << "%\n";
00275       GmmDynSparse  gmmT3(rows,cols);
00276       GmmSparse m1(rows,cols), m2(rows,cols), m3(rows,cols);
00277       eiToGmm(sm1, m1);
00278       eiToGmm(sm2, m2);
00279 
00280       timer.reset();
00281       timer.start();
00282       for (int k=0; k<REPEAT; ++k)
00283         gmm::mult(m1, m2, gmmT3);
00284       timer.stop();
00285       std::cout << "   a * b:\t" << timer.value() << endl;
00286 
00287 //       timer.reset();
00288 //       timer.start();
00289 //       for (int k=0; k<REPEAT; ++k)
00290 //         gmm::mult(gmm::transposed(m1), m2, gmmT3);
00291 //       timer.stop();
00292 //       std::cout << "   a' * b:\t" << timer.value() << endl;
00293 //
00294 //       if (rows<500)
00295 //       {
00296 //         timer.reset();
00297 //         timer.start();
00298 //         for (int k=0; k<REPEAT; ++k)
00299 //           gmm::mult(gmm::transposed(m1), gmm::transposed(m2), gmmT3);
00300 //         timer.stop();
00301 //         std::cout << "   a' * b':\t" << timer.value() << endl;
00302 //
00303 //         timer.reset();
00304 //         timer.start();
00305 //         for (int k=0; k<REPEAT; ++k)
00306 //           gmm::mult(m1, gmm::transposed(m2), gmmT3);
00307 //         timer.stop();
00308 //         std::cout << "   a * b':\t" << timer.value() << endl;
00309 //       }
00310 //       else
00311 //       {
00312 //         std::cout << "   a' * b':\t" << "forever" << endl;
00313 //         std::cout << "   a * b':\t" << "forever" << endl;
00314 //       }
00315     }
00316     #endif
00317 
00318     // MTL4
00319     #ifndef NOMTL
00320     {
00321       std::cout << "MTL4\t" << nnzPerCol << "%\n";
00322       MtlSparse m1(rows,cols), m2(rows,cols), m3(rows,cols);
00323       eiToMtl(sm1, m1);
00324       eiToMtl(sm2, m2);
00325 
00326       timer.reset();
00327       timer.start();
00328       for (int k=0; k<REPEAT; ++k)
00329         m3 = m1 * m2;
00330       timer.stop();
00331       std::cout << "   a * b:\t" << timer.value() << endl;
00332 
00333 //       timer.reset();
00334 //       timer.start();
00335 //       for (int k=0; k<REPEAT; ++k)
00336 //         m3 = trans(m1) * m2;
00337 //       timer.stop();
00338 //       std::cout << "   a' * b:\t" << timer.value() << endl;
00339 //
00340 //       timer.reset();
00341 //       timer.start();
00342 //       for (int k=0; k<REPEAT; ++k)
00343 //         m3 = trans(m1) * trans(m2);
00344 //       timer.stop();
00345 //       std::cout << "  a' * b':\t" << timer.value() << endl;
00346 //
00347 //       timer.reset();
00348 //       timer.start();
00349 //       for (int k=0; k<REPEAT; ++k)
00350 //         m3 = m1 * trans(m2);
00351 //       timer.stop();
00352 //       std::cout << "   a * b' :\t" << timer.value() << endl;
00353     }
00354     #endif
00355 
00356     std::cout << "\n\n";
00357   }
00358 
00359   return 0;
00360 }
00361 
00362 
00363 


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