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
main(int argc,char * argv[])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
234