1 2 //g++ -O3 -g0 -DNDEBUG sparse_product.cpp -I.. -I/home/gael/Coding/LinearAlgebra/mtl4/ -DDENSITY=0.005 -DSIZE=10000 && ./a.out 3 //g++ -O3 -g0 -DNDEBUG sparse_product.cpp -I.. -I/home/gael/Coding/LinearAlgebra/mtl4/ -DDENSITY=0.05 -DSIZE=2000 && ./a.out 4 // -DNOGMM -DNOMTL -DCSPARSE 5 // -I /home/gael/Coding/LinearAlgebra/CSparse/Include/ /home/gael/Coding/LinearAlgebra/CSparse/Lib/libcsparse.a 6 7 #include <typeinfo> 8 9 #ifndef SIZE 10 #define SIZE 1000000 11 #endif 12 13 #ifndef NNZPERCOL 14 #define NNZPERCOL 6 15 #endif 16 17 #ifndef REPEAT 18 #define REPEAT 1 19 #endif 20 21 #include <algorithm> 22 #include "BenchTimer.h" 23 #include "BenchUtil.h" 24 #include "BenchSparseUtil.h" 25 26 #ifndef NBTRIES 27 #define NBTRIES 1 28 #endif 29 30 #define BENCH(X) \ 31 timer.reset(); \ 32 for (int _j=0; _j<NBTRIES; ++_j) { \ 33 timer.start(); \ 34 for (int _k=0; _k<REPEAT; ++_k) { \ 35 X \ 36 } timer.stop(); } 37 38 // #ifdef MKL 39 // 40 // #include "mkl_types.h" 41 // #include "mkl_spblas.h" 42 // 43 // template<typename Lhs,typename Rhs,typename Res> 44 // void mkl_multiply(const Lhs& lhs, const Rhs& rhs, Res& res) 45 // { 46 // char n = 'N'; 47 // float alpha = 1; 48 // char matdescra[6]; 49 // matdescra[0] = 'G'; 50 // matdescra[1] = 0; 51 // matdescra[2] = 0; 52 // matdescra[3] = 'C'; 53 // mkl_scscmm(&n, lhs.rows(), rhs.cols(), lhs.cols(), &alpha, matdescra, 54 // lhs._valuePtr(), lhs._innerIndexPtr(), lhs.outerIndexPtr(), 55 // pntre, b, &ldb, &beta, c, &ldc); 56 // // mkl_somatcopy('C', 'T', lhs.rows(), lhs.cols(), 1, 57 // // lhs._valuePtr(), lhs.rows(), DST, dst_stride); 58 // } 59 // 60 // #endif 61 62 63 #ifdef CSPARSE 64 cs* cs_sorted_multiply(const cs* a, const cs* b) 65 { 66 // return cs_multiply(a,b); 67 68 cs* A = cs_transpose(a, 1); 69 cs* B = cs_transpose(b, 1); 70 cs* D = cs_multiply(B,A); /* D = B'*A' */ 71 cs_spfree (A) ; 72 cs_spfree (B) ; 73 cs_dropzeros (D) ; /* drop zeros from D */ 74 cs* C = cs_transpose (D, 1) ; /* C = D', so that C is sorted */ 75 cs_spfree (D) ; 76 return C; 77 78 // cs* A = cs_transpose(a, 1); 79 // cs* C = cs_transpose(A, 1); 80 // return C; 81 } 82 83 cs* cs_sorted_multiply2(const cs* a, const cs* b) 84 { 85 cs* D = cs_multiply(a,b); 86 cs* E = cs_transpose(D,1); 87 cs_spfree(D); 88 cs* C = cs_transpose(E,1); 89 cs_spfree(E); 90 return C; 91 } 92 #endif 93 94 void bench_sort(); 95 96 int main(int argc, char *argv[]) 97 { 98 // bench_sort(); 99 100 int rows = SIZE; 101 int cols = SIZE; 102 float density = DENSITY; 103 104 EigenSparseMatrix sm1(rows,cols), sm2(rows,cols), sm3(rows,cols), sm4(rows,cols); 105 106 BenchTimer timer; 107 for (int nnzPerCol = NNZPERCOL; nnzPerCol>1; nnzPerCol/=1.1) 108 { 109 sm1.setZero(); 110 sm2.setZero(); 111 fillMatrix2(nnzPerCol, rows, cols, sm1); 112 fillMatrix2(nnzPerCol, rows, cols, sm2); 113 // std::cerr << "filling OK\n"; 114 115 // dense matrices 116 #ifdef DENSEMATRIX 117 { 118 std::cout << "Eigen Dense\t" << nnzPerCol << "%\n"; 119 DenseMatrix m1(rows,cols), m2(rows,cols), m3(rows,cols); 120 eiToDense(sm1, m1); 121 eiToDense(sm2, m2); 122 123 timer.reset(); 124 timer.start(); 125 for (int k=0; k<REPEAT; ++k) 126 m3 = m1 * m2; 127 timer.stop(); 128 std::cout << " a * b:\t" << timer.value() << endl; 129 130 timer.reset(); 131 timer.start(); 132 for (int k=0; k<REPEAT; ++k) 133 m3 = m1.transpose() * m2; 134 timer.stop(); 135 std::cout << " a' * b:\t" << timer.value() << endl; 136 137 timer.reset(); 138 timer.start(); 139 for (int k=0; k<REPEAT; ++k) 140 m3 = m1.transpose() * m2.transpose(); 141 timer.stop(); 142 std::cout << " a' * b':\t" << timer.value() << endl; 143 144 timer.reset(); 145 timer.start(); 146 for (int k=0; k<REPEAT; ++k) 147 m3 = m1 * m2.transpose(); 148 timer.stop(); 149 std::cout << " a * b':\t" << timer.value() << endl; 150 } 151 #endif 152 153 // eigen sparse matrices 154 { 155 std::cout << "Eigen sparse\t" << sm1.nonZeros()/(float(sm1.rows())*float(sm1.cols()))*100 << "% * " 156 << sm2.nonZeros()/(float(sm2.rows())*float(sm2.cols()))*100 << "%\n"; 157 158 BENCH(sm3 = sm1 * sm2; ) 159 std::cout << " a * b:\t" << timer.value() << endl; 160 161 // BENCH(sm3 = sm1.transpose() * sm2; ) 162 // std::cout << " a' * b:\t" << timer.value() << endl; 163 // // 164 // BENCH(sm3 = sm1.transpose() * sm2.transpose(); ) 165 // std::cout << " a' * b':\t" << timer.value() << endl; 166 // // 167 // BENCH(sm3 = sm1 * sm2.transpose(); ) 168 // std::cout << " a * b' :\t" << timer.value() << endl; 169 170 171 // std::cout << "\n"; 172 // 173 // BENCH( sm3._experimentalNewProduct(sm1, sm2); ) 174 // std::cout << " a * b:\t" << timer.value() << endl; 175 // 176 // BENCH(sm3._experimentalNewProduct(sm1.transpose(),sm2); ) 177 // std::cout << " a' * b:\t" << timer.value() << endl; 178 // // 179 // BENCH(sm3._experimentalNewProduct(sm1.transpose(),sm2.transpose()); ) 180 // std::cout << " a' * b':\t" << timer.value() << endl; 181 // // 182 // BENCH(sm3._experimentalNewProduct(sm1, sm2.transpose());) 183 // std::cout << " a * b' :\t" << timer.value() << endl; 184 } 185 186 // eigen dyn-sparse matrices 187 /*{ 188 DynamicSparseMatrix<Scalar> m1(sm1), m2(sm2), m3(sm3); 189 std::cout << "Eigen dyn-sparse\t" << m1.nonZeros()/(float(m1.rows())*float(m1.cols()))*100 << "% * " 190 << m2.nonZeros()/(float(m2.rows())*float(m2.cols()))*100 << "%\n"; 191 192 // timer.reset(); 193 // timer.start(); 194 BENCH(for (int k=0; k<REPEAT; ++k) m3 = m1 * m2;) 195 // timer.stop(); 196 std::cout << " a * b:\t" << timer.value() << endl; 197 // std::cout << sm3 << "\n"; 198 199 timer.reset(); 200 timer.start(); 201 // std::cerr << "transpose...\n"; 202 // EigenSparseMatrix sm4 = sm1.transpose(); 203 // std::cout << sm4.nonZeros() << " == " << sm1.nonZeros() << "\n"; 204 // exit(1); 205 // std::cerr << "transpose OK\n"; 206 // std::cout << sm1 << "\n\n" << sm1.transpose() << "\n\n" << sm4.transpose() << "\n\n"; 207 BENCH(for (int k=0; k<REPEAT; ++k) m3 = m1.transpose() * m2;) 208 // timer.stop(); 209 std::cout << " a' * b:\t" << timer.value() << endl; 210 211 // timer.reset(); 212 // timer.start(); 213 BENCH( for (int k=0; k<REPEAT; ++k) m3 = m1.transpose() * m2.transpose(); ) 214 // timer.stop(); 215 std::cout << " a' * b':\t" << timer.value() << endl; 216 217 // timer.reset(); 218 // timer.start(); 219 BENCH( for (int k=0; k<REPEAT; ++k) m3 = m1 * m2.transpose(); ) 220 // timer.stop(); 221 std::cout << " a * b' :\t" << timer.value() << endl; 222 }*/ 223 224 // CSparse 225 #ifdef CSPARSE 226 { 227 std::cout << "CSparse \t" << nnzPerCol << "%\n"; 228 cs *m1, *m2, *m3; 229 eiToCSparse(sm1, m1); 230 eiToCSparse(sm2, m2); 231 232 BENCH( 233 { 234 m3 = cs_sorted_multiply(m1, m2); 235 if (!m3) 236 { 237 std::cerr << "cs_multiply failed\n"; 238 } 239 // cs_print(m3, 0); 240 cs_spfree(m3); 241 } 242 ); 243 // timer.stop(); 244 std::cout << " a * b:\t" << timer.value() << endl; 245 246 // BENCH( { m3 = cs_sorted_multiply2(m1, m2); cs_spfree(m3); } ); 247 // std::cout << " a * b:\t" << timer.value() << endl; 248 } 249 #endif 250 251 #ifndef NOUBLAS 252 { 253 std::cout << "ublas\t" << nnzPerCol << "%\n"; 254 UBlasSparse m1(rows,cols), m2(rows,cols), m3(rows,cols); 255 eiToUblas(sm1, m1); 256 eiToUblas(sm2, m2); 257 258 BENCH(boost::numeric::ublas::prod(m1, m2, m3);); 259 std::cout << " a * b:\t" << timer.value() << endl; 260 } 261 #endif 262 263 // GMM++ 264 #ifndef NOGMM 265 { 266 std::cout << "GMM++ sparse\t" << nnzPerCol << "%\n"; 267 GmmDynSparse gmmT3(rows,cols); 268 GmmSparse m1(rows,cols), m2(rows,cols), m3(rows,cols); 269 eiToGmm(sm1, m1); 270 eiToGmm(sm2, m2); 271 272 BENCH(gmm::mult(m1, m2, gmmT3);); 273 std::cout << " a * b:\t" << timer.value() << endl; 274 275 // BENCH(gmm::mult(gmm::transposed(m1), m2, gmmT3);); 276 // std::cout << " a' * b:\t" << timer.value() << endl; 277 // 278 // if (rows<500) 279 // { 280 // BENCH(gmm::mult(gmm::transposed(m1), gmm::transposed(m2), gmmT3);); 281 // std::cout << " a' * b':\t" << timer.value() << endl; 282 // 283 // BENCH(gmm::mult(m1, gmm::transposed(m2), gmmT3);); 284 // std::cout << " a * b':\t" << timer.value() << endl; 285 // } 286 // else 287 // { 288 // std::cout << " a' * b':\t" << "forever" << endl; 289 // std::cout << " a * b':\t" << "forever" << endl; 290 // } 291 } 292 #endif 293 294 // MTL4 295 #ifndef NOMTL 296 { 297 std::cout << "MTL4\t" << nnzPerCol << "%\n"; 298 MtlSparse m1(rows,cols), m2(rows,cols), m3(rows,cols); 299 eiToMtl(sm1, m1); 300 eiToMtl(sm2, m2); 301 302 BENCH(m3 = m1 * m2;); 303 std::cout << " a * b:\t" << timer.value() << endl; 304 305 // BENCH(m3 = trans(m1) * m2;); 306 // std::cout << " a' * b:\t" << timer.value() << endl; 307 // 308 // BENCH(m3 = trans(m1) * trans(m2);); 309 // std::cout << " a' * b':\t" << timer.value() << endl; 310 // 311 // BENCH(m3 = m1 * trans(m2);); 312 // std::cout << " a * b' :\t" << timer.value() << endl; 313 } 314 #endif 315 316 std::cout << "\n\n"; 317 } 318 319 return 0; 320 } 321 322 323 324
