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_BENCHres = 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
