spmv.cpp
Go to the documentation of this file.
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 
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 
Eigen::SparseMatrix< Scalar >
MtlSparse
mtl::compressed2D< Scalar, mtl::matrix::parameters< mtl::tag::col_major > > MtlSparse
Definition: BenchSparseUtil.h:86
eiToDense
void eiToDense(const EigenSparseMatrix &src, DenseMatrix &dst)
Definition: BenchSparseUtil.h:62
trans
static char trans
Definition: blas_interface.hh:58
SPMV_BENCH
#define SPMV_BENCH(CODE)
Definition: spmv.cpp:11
eiToMtl
void eiToMtl(const EigenSparseMatrix &src, MtlSparse &dst)
Definition: BenchSparseUtil.h:88
res
cout<< "Here is the matrix m:"<< endl<< m<< endl;Matrix< ptrdiff_t, 3, 1 > res
Definition: PartialRedux_count.cpp:3
rows
int rows
Definition: Tutorial_commainit_02.cpp:1
size
Scalar Scalar int size
Definition: benchVecAdd.cpp:17
Eigen::BenchTimer
Definition: BenchTimer.h:59
eiToUblas
void eiToUblas(const EigenSparseMatrix &src, UBlasSparse &dst)
Definition: BenchSparseUtil.h:128
Eigen::Map
A matrix or vector expression mapping an existing array of data.
Definition: Map.h:94
DenseVector
Matrix< Scalar, Dynamic, 1 > DenseVector
Definition: BenchSparseUtil.h:24
main
int main(int argc, char *argv[])
Definition: spmv.cpp:37
BenchTimer.h
GmmSparse
gmm::csc_matrix< Scalar > GmmSparse
Definition: BenchSparseUtil.h:72
DenseMatrix
Matrix< Scalar, Dynamic, Dynamic > DenseMatrix
Definition: BenchSparseUtil.h:23
Eigen::SparseMatrix::setZero
void setZero()
Definition: SparseMatrix.h:253
eiToGmm
void eiToGmm(const EigenSparseMatrix &src, GmmSparse &dst)
Definition: BenchSparseUtil.h:74
cols
int cols
Definition: Tutorial_commainit_02.cpp:1
align_3::t
Point2 t(10, 10)
i
int i
Definition: BiCGSTAB_step_by_step.cpp:9
BenchSparseUtil.h
fillMatrix2
void fillMatrix2(int nnzPerCol, int rows, int cols, EigenSparseMatrix &dst)
Definition: BenchSparseUtil.h:42


gtsam
Author(s):
autogenerated on Fri Nov 1 2024 03:36:24