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
37int 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
234