[136] | 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 |
|
---|