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 #ifndef SIZE 7 #define SIZE 650000 8 #endif 9 10 #ifndef DENSITY 11 #define DENSITY 0.01 12 #endif 13 14 #ifndef REPEAT 15 #define REPEAT 1 16 #endif 17 18 #include "BenchSparseUtil.h" 19 20 #ifndef MINDENSITY 21 #define MINDENSITY 0.0004 22 #endif 23 24 #ifndef NBTRIES 25 #define NBTRIES 10 26 #endif 27 28 #define BENCH(X) \ 29 timer.reset(); \ 30 for (int _j=0; _j<NBTRIES; ++_j) { \ 31 timer.start(); \ 32 for (int _k=0; _k<REPEAT; ++_k) { \ 33 X \ 34 } timer.stop(); } 35 36 37 #ifdef CSPARSE 38 cs* cs_sorted_multiply(const cs* a, const cs* b) 39 { 40 cs* A = cs_transpose (a, 1) ; 41 cs* B = cs_transpose (b, 1) ; 42 cs* D = cs_multiply (B,A) ; /* D = B'*A' */ 43 cs_spfree (A) ; 44 cs_spfree (B) ; 45 cs_dropzeros (D) ; /* drop zeros from D */ 46 cs* C = cs_transpose (D, 1) ; /* C = D', so that C is sorted */ 47 cs_spfree (D) ; 48 return C; 49 } 50 #endif 51 52 int main(int argc, char *argv[]) 53 { 54 int rows = SIZE; 55 int cols = SIZE; 56 float density = DENSITY; 57 58 EigenSparseMatrix sm1(rows,cols); 59 DenseVector v1(cols), v2(cols); 60 v1.setRandom(); 61 62 BenchTimer timer; 63 for (float density = DENSITY; density>=MINDENSITY; density*=0.5) 64 { 65 //fillMatrix(density, rows, cols, sm1); 66 fillMatrix2(7, rows, cols, sm1); 67 68 // dense matrices 69 #ifdef DENSEMATRIX 70 { 71 std::cout << "Eigen Dense\t" << density*100 << "%\n"; 72 DenseMatrix m1(rows,cols); 73 eiToDense(sm1, m1); 74 75 timer.reset(); 76 timer.start(); 77 for (int k=0; k<REPEAT; ++k) 78 v2 = m1 * v1; 79 timer.stop(); 80 std::cout << " a * v:\t" << timer.best() << " " << double(REPEAT)/timer.best() << " * / sec " << endl; 81 82 timer.reset(); 83 timer.start(); 84 for (int k=0; k<REPEAT; ++k) 85 v2 = m1.transpose() * v1; 86 timer.stop(); 87 std::cout << " a' * v:\t" << timer.best() << endl; 88 } 89 #endif 90 91 // eigen sparse matrices 92 { 93 std::cout << "Eigen sparse\t" << sm1.nonZeros()/float(sm1.rows()*sm1.cols())*100 << "%\n"; 94 95 BENCH(asm("#myc"); v2 = sm1 * v1; asm("#myd");) 96 std::cout << " a * v:\t" << timer.best()/REPEAT << " " << double(REPEAT)/timer.best(REAL_TIMER) << " * / sec " << endl; 97 98 99 BENCH( { asm("#mya"); v2 = sm1.transpose() * v1; asm("#myb"); }) 100 101 std::cout << " a' * v:\t" << timer.best()/REPEAT << endl; 102 } 103 104 // { 105 // DynamicSparseMatrix<Scalar> m1(sm1); 106 // std::cout << "Eigen dyn-sparse\t" << m1.nonZeros()/float(m1.rows()*m1.cols())*100 << "%\n"; 107 // 108 // BENCH(for (int k=0; k<REPEAT; ++k) v2 = m1 * v1;) 109 // std::cout << " a * v:\t" << timer.value() << endl; 110 // 111 // BENCH(for (int k=0; k<REPEAT; ++k) v2 = m1.transpose() * v1;) 112 // std::cout << " a' * v:\t" << timer.value() << endl; 113 // } 114 115 // GMM++ 116 #ifndef NOGMM 117 { 118 std::cout << "GMM++ sparse\t" << density*100 << "%\n"; 119 //GmmDynSparse gmmT3(rows,cols); 120 GmmSparse m1(rows,cols); 121 eiToGmm(sm1, m1); 122 123 std::vector<Scalar> gmmV1(cols), gmmV2(cols); 124 Map<Matrix<Scalar,Dynamic,1> >(&gmmV1[0], cols) = v1; 125 Map<Matrix<Scalar,Dynamic,1> >(&gmmV2[0], cols) = v2; 126 127 BENCH( asm("#myx"); gmm::mult(m1, gmmV1, gmmV2); asm("#myy"); ) 128 std::cout << " a * v:\t" << timer.value() << endl; 129 130 BENCH( gmm::mult(gmm::transposed(m1), gmmV1, gmmV2); ) 131 std::cout << " a' * v:\t" << timer.value() << endl; 132 } 133 #endif 134 135 #ifndef NOUBLAS 136 { 137 std::cout << "ublas sparse\t" << density*100 << "%\n"; 138 UBlasSparse m1(rows,cols); 139 eiToUblas(sm1, m1); 140 141 boost::numeric::ublas::vector<Scalar> uv1, uv2; 142 eiToUblasVec(v1,uv1); 143 eiToUblasVec(v2,uv2); 144 145 // std::vector<Scalar> gmmV1(cols), gmmV2(cols); 146 // Map<Matrix<Scalar,Dynamic,1> >(&gmmV1[0], cols) = v1; 147 // Map<Matrix<Scalar,Dynamic,1> >(&gmmV2[0], cols) = v2; 148 149 BENCH( uv2 = boost::numeric::ublas::prod(m1, uv1); ) 150 std::cout << " a * v:\t" << timer.value() << endl; 151 152 // BENCH( boost::ublas::prod(gmm::transposed(m1), gmmV1, gmmV2); ) 153 // std::cout << " a' * v:\t" << timer.value() << endl; 154 } 155 #endif 156 157 // MTL4 158 #ifndef NOMTL 159 { 160 std::cout << "MTL4\t" << density*100 << "%\n"; 161 MtlSparse m1(rows,cols); 162 eiToMtl(sm1, m1); 163 mtl::dense_vector<Scalar> mtlV1(cols, 1.0); 164 mtl::dense_vector<Scalar> mtlV2(cols, 1.0); 165 166 timer.reset(); 167 timer.start(); 168 for (int k=0; k<REPEAT; ++k) 169 mtlV2 = m1 * mtlV1; 170 timer.stop(); 171 std::cout << " a * v:\t" << timer.value() << endl; 172 173 timer.reset(); 174 timer.start(); 175 for (int k=0; k<REPEAT; ++k) 176 mtlV2 = trans(m1) * mtlV1; 177 timer.stop(); 178 std::cout << " a' * v:\t" << timer.value() << endl; 179 } 180 #endif 181 182 std::cout << "\n\n"; 183 } 184 185 return 0; 186 } 187 188
