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