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

   1
   2 //g++-4.4 -DNOMTL  -Wl,-rpath /usr/local/lib/oski -L /usr/local/lib/oski/ -l oski -l oski_util -l oski_util_Tid  -DOSKI -I ~/Coding/LinearAlgebra/mtl4/  spmv.cpp  -I .. -O2 -DNDEBUG -lrt  -lm -l oski_mat_CSC_Tid  -loskilt && ./a.out r200000 c200000 n100 t1 p1
   3
   4 #define SCALAR double
   5
   6 #include <iostream>
   7 #include <algorithm>
   8 #include "BenchTimer.h"
   9 #include "BenchSparseUtil.h"
  10
  11 #define SPMV_BENCH(CODE) BENCH(t,tries,repeats,CODE);
  12
  13 // #ifdef MKL
  14 //
  15 // #include "mkl_types.h"
  16 // #include "mkl_spblas.h"
  17 //
  18 // template<typename Lhs,typename Rhs,typename Res>
  19 // void mkl_multiply(const Lhs& lhs, const Rhs& rhs, Res& res)
  20 // {
  21 //   char n = 'N';
  22 //   float alpha = 1;
  23 //   char matdescra[6];
  24 //   matdescra[0] = 'G';
  25 //   matdescra[1] = 0;
  26 //   matdescra[2] = 0;
  27 //   matdescra[3] = 'C';
  28 //   mkl_scscmm(&n, lhs.rows(), rhs.cols(), lhs.cols(), &alpha, matdescra,
  29 //              lhs._valuePtr(), lhs._innerIndexPtr(), lhs.outerIndexPtr(),
  30 //              pntre, b, &ldb, &beta, c, &ldc);
  31 // //   mkl_somatcopy('C', 'T', lhs.rows(), lhs.cols(), 1,
  32 // //                 lhs._valuePtr(), lhs.rows(), DST, dst_stride);
  33 // }
  34 //
  35 // #endif
  36
  37 int main(int argc, char *argv[])
  38 {
  39   int size = 10000;
  40   int rows = size;
  41   int cols = size;
  42   int nnzPerCol = 40;
  43   int tries = 2;
  44   int repeats = 2;
  45
  46   bool need_help = false;
  47   for(int i = 1; i < argc; i++)
  48   {
  49     if(argv[i][0] == 'r')
  50     {
  51       rows = atoi(argv[i]+1);
  52     }
  53     else if(argv[i][0] == 'c')
  54     {
  55       cols = atoi(argv[i]+1);
  56     }
  57     else if(argv[i][0] == 'n')
  58     {
  59       nnzPerCol = atoi(argv[i]+1);
  60     }
  61     else if(argv[i][0] == 't')
  62     {
  63       tries = atoi(argv[i]+1);
  64     }
  65     else if(argv[i][0] == 'p')
  66     {
  67       repeats = atoi(argv[i]+1);
  68     }
  69     else
  70     {
  71       need_help = true;
  72     }
  73   }
  74   if(need_help)
  75   {
  76     std::cout << argv[0] << " r<nb rows> c<nb columns> n<non zeros per column> t<nb tries> p<nb repeats>\n";
  77     return 1;
  78   }
  79
  80   std::cout << "SpMV " << rows << " x " << cols << " with " << nnzPerCol << " non zeros per column. (" << repeats << " repeats, and " << tries << " tries)\n\n";
  81
  82   EigenSparseMatrix sm(rows,cols);
  83   DenseVector dv(cols), res(rows);
  84   dv.setRandom();
  85
  86   BenchTimer t;
  87   while (nnzPerCol>=4)
  88   {
  89     std::cout << "nnz: " << nnzPerCol << "\n";
  90     sm.setZero();
  91     fillMatrix2(nnzPerCol, rows, cols, sm);
  92
  93     // dense matrices
  94     #ifdef DENSEMATRIX
  95     {
  96       DenseMatrix dm(rows,cols), (rows,cols);
  97       eiToDense(sm, dm);
  98
  99       SPMV_BENCH(res = dm * sm);
 100       std::cout << "Dense       " << t.value()/repeats << "\t";
 101
 102       SPMV_BENCH(res = dm.transpose() * sm);
 103       std::cout << t.value()/repeats << endl;
 104     }
 105     #endif
 106
 107     // eigen sparse matrices
 108     {
 109       SPMV_BENCH(res.noalias() += sm * dv; )
 110       std::cout << "Eigen       " << t.value()/repeats << "\t";
 111
 112       SPMV_BENCH(res.noalias() += sm.transpose() * dv; )
 113       std::cout << t.value()/repeats << endl;
 114     }
 115
 116     // CSparse
 117     #ifdef CSPARSE
 118     {
 119       std::cout << "CSparse \n";
 120       cs *csm;
 121       eiToCSparse(sm, csm);
 122
 123 //       BENCH();
 124 //       timer.stop();
 125 //       std::cout << "   a * b:\t" << timer.value() << endl;
 126
 127 //       BENCH( { m3 = cs_sorted_multiply2(m1, m2); cs_spfree(m3); } );
 128 //       std::cout << "   a * b:\t" << timer.value() << endl;
 129     }
 130     #endif
 131
 132     #ifdef OSKI
 133     {
 134       oski_matrix_t om;
 135       oski_vecview_t ov, ores;
 136       oski_Init();
 137       om = oski_CreateMatCSC(sm._outerIndexPtr(), sm._innerIndexPtr(), sm._valuePtr(), rows, cols,
 138                              SHARE_INPUTMAT, 1, INDEX_ZERO_BASED);
 139       ov = oski_CreateVecView(dv.data(), cols, STRIDE_UNIT);
 140       ores = oski_CreateVecView(res.data(), rows, STRIDE_UNIT);
 141
 142       SPMV_BENCH( oski_MatMult(om, OP_NORMAL, 1, ov, 0, ores) );
 143       std::cout << "OSKI        " << t.value()/repeats << "\t";
 144
 145       SPMV_BENCH( oski_MatMult(om, OP_TRANS, 1, ov, 0, ores) );
 146       std::cout << t.value()/repeats << "\n";
 147
 148       // tune
 149       t.reset();
 150       t.start();
 151       oski_SetHintMatMult(om, OP_NORMAL, 1.0, SYMBOLIC_VEC, 0.0, SYMBOLIC_VEC, ALWAYS_TUNE_AGGRESSIVELY);
 152       oski_TuneMat(om);
 153       t.stop();
 154       double tuning = t.value();
 155
 156       SPMV_BENCH( oski_MatMult(om, OP_NORMAL, 1, ov, 0, ores) );
 157       std::cout << "OSKI tuned  " << t.value()/repeats << "\t";
 158
 159       SPMV_BENCH( oski_MatMult(om, OP_TRANS, 1, ov, 0, ores) );
 160       std::cout << t.value()/repeats << "\t(" << tuning <<  ")\n";
 161
 162
 163       oski_DestroyMat(om);
 164       oski_DestroyVecView(ov);
 165       oski_DestroyVecView(ores);
 166       oski_Close();
 167     }
 168     #endif
 169
 170     #ifndef NOUBLAS
 171     {
 172       using namespace boost::numeric;
 173       UblasMatrix um(rows,cols);
 174       eiToUblas(sm, um);
 175
 176       boost::numeric::ublas::vector<Scalar> uv(cols), ures(rows);
 177       Map<Matrix<Scalar,Dynamic,1> >(&uv[0], cols) = dv;
 178       Map<Matrix<Scalar,Dynamic,1> >(&ures[0], rows) = res;
 179
 180       SPMV_BENCH(ublas::axpy_prod(um, uv, ures, true));
 181       std::cout << "ublas       " << t.value()/repeats << "\t";
 182
 183       SPMV_BENCH(ublas::axpy_prod(boost::numeric::ublas::trans(um), uv, ures, true));
 184       std::cout << t.value()/repeats << endl;
 185     }
 186     #endif
 187
 188     // GMM++
 189     #ifndef NOGMM
 190     {
 191       GmmSparse gm(rows,cols);
 192       eiToGmm(sm, gm);
 193
 194       std::vector<Scalar> gv(cols), gres(rows);
 195       Map<Matrix<Scalar,Dynamic,1> >(&gv[0], cols) = dv;
 196       Map<Matrix<Scalar,Dynamic,1> >(&gres[0], rows) = res;
 197
 198       SPMV_BENCH(gmm::mult(gm, gv, gres));
 199       std::cout << "GMM++       " << t.value()/repeats << "\t";
 200
 201       SPMV_BENCH(gmm::mult(gmm::transposed(gm), gv, gres));
 202       std::cout << t.value()/repeats << endl;
 203     }
 204     #endif
 205
 206     // MTL4
 207     #ifndef NOMTL
 208     {
 209       MtlSparse mm(rows,cols);
 210       eiToMtl(sm, mm);
 211       mtl::dense_vector<Scalar> mv(cols, 1.0);
 212       mtl::dense_vector<Scalar> mres(rows, 1.0);
 213
 214       SPMV_BENCH(mres = mm * mv);
 215       std::cout << "MTL4        " << t.value()/repeats << "\t";
 216
 217       SPMV_BENCH(mres = trans(mm) * mv);
 218       std::cout << t.value()/repeats << endl;
 219     }
 220     #endif
 221
 222     std::cout << "\n";
 223
 224     if(nnzPerCol==1)
 225       break;
 226     nnzPerCol -= nnzPerCol/2;
 227   }
 228
 229   return 0;
 230 }
 231
 232
 233