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 5 // -I /home/gael/Coding/LinearAlgebra/CSparse/Include/ /home/gael/Coding/LinearAlgebra/CSparse/Lib/libcsparse.a 6 7 #ifndef SIZE 8 #define SIZE 10000 9 #endif 10 11 #ifndef DENSITY 12 #define DENSITY 0.01 13 #endif 14 15 #ifndef REPEAT 16 #define REPEAT 1 17 #endif 18 19 #include "BenchSparseUtil.h" 20 21 #ifndef MINDENSITY 22 #define MINDENSITY 0.0004 23 #endif 24 25 #ifndef NBTRIES 26 #define NBTRIES 10 27 #endif 28 29 #define BENCH(X) \ 30 timer.reset(); \ 31 for (int _j=0; _j<NBTRIES; ++_j) { \ 32 timer.start(); \ 33 for (int _k=0; _k<REPEAT; ++_k) { \ 34 X \ 35 } timer.stop(); } 36 37 typedef SparseMatrix<Scalar,UpperTriangular> EigenSparseTriMatrix; 38 typedef SparseMatrix<Scalar,RowMajorBit|UpperTriangular> EigenSparseTriMatrixRow; 39 40 void fillMatrix(float density, int rows, int cols, EigenSparseTriMatrix& dst) 41 { 42 dst.startFill(rows*cols*density); 43 for(int j = 0; j < cols; j++) 44 { 45 for(int i = 0; i < j; i++) 46 { 47 Scalar v = (internal::random<float>(0,1) < density) ? internal::random<Scalar>() : 0; 48 if (v!=0) 49 dst.fill(i,j) = v; 50 } 51 dst.fill(j,j) = internal::random<Scalar>(); 52 } 53 dst.endFill(); 54 } 55 56 int main(int argc, char *argv[]) 57 { 58 int rows = SIZE; 59 int cols = SIZE; 60 float density = DENSITY; 61 BenchTimer timer; 62 #if 1 63 EigenSparseTriMatrix sm1(rows,cols); 64 typedef Matrix<Scalar,Dynamic,1> DenseVector; 65 DenseVector b = DenseVector::Random(cols); 66 DenseVector x = DenseVector::Random(cols); 67 68 bool densedone = false; 69 70 for (float density = DENSITY; density>=MINDENSITY; density*=0.5) 71 { 72 EigenSparseTriMatrix sm1(rows, cols); 73 fillMatrix(density, rows, cols, sm1); 74 75 // dense matrices 76 #ifdef DENSEMATRIX 77 if (!densedone) 78 { 79 densedone = true; 80 std::cout << "Eigen Dense\t" << density*100 << "%\n"; 81 DenseMatrix m1(rows,cols); 82 Matrix<Scalar,Dynamic,Dynamic,Dynamic,Dynamic,RowMajorBit> m2(rows,cols); 83 eiToDense(sm1, m1); 84 m2 = m1; 85 86 BENCH(x = m1.marked<UpperTriangular>().solveTriangular(b);) 87 std::cout << " colmajor^-1 * b:\t" << timer.value() << endl; 88 // std::cerr << x.transpose() << "\n"; 89 90 BENCH(x = m2.marked<UpperTriangular>().solveTriangular(b);) 91 std::cout << " rowmajor^-1 * b:\t" << timer.value() << endl; 92 // std::cerr << x.transpose() << "\n"; 93 } 94 #endif 95 96 // eigen sparse matrices 97 { 98 std::cout << "Eigen sparse\t" << density*100 << "%\n"; 99 EigenSparseTriMatrixRow sm2 = sm1; 100 101 BENCH(x = sm1.solveTriangular(b);) 102 std::cout << " colmajor^-1 * b:\t" << timer.value() << endl; 103 // std::cerr << x.transpose() << "\n"; 104 105 BENCH(x = sm2.solveTriangular(b);) 106 std::cout << " rowmajor^-1 * b:\t" << timer.value() << endl; 107 // std::cerr << x.transpose() << "\n"; 108 109 // x = b; 110 // BENCH(sm1.inverseProductInPlace(x);) 111 // std::cout << " colmajor^-1 * b:\t" << timer.value() << " (inplace)" << endl; 112 // std::cerr << x.transpose() << "\n"; 113 // 114 // x = b; 115 // BENCH(sm2.inverseProductInPlace(x);) 116 // std::cout << " rowmajor^-1 * b:\t" << timer.value() << " (inplace)" << endl; 117 // std::cerr << x.transpose() << "\n"; 118 } 119 120 121 122 // CSparse 123 #ifdef CSPARSE 124 { 125 std::cout << "CSparse \t" << density*100 << "%\n"; 126 cs *m1; 127 eiToCSparse(sm1, m1); 128 129 BENCH(x = b; if (!cs_lsolve (m1, x.data())){std::cerr << "cs_lsolve failed\n"; break;}; ) 130 std::cout << " colmajor^-1 * b:\t" << timer.value() << endl; 131 } 132 #endif 133 134 // GMM++ 135 #ifndef NOGMM 136 { 137 std::cout << "GMM++ sparse\t" << density*100 << "%\n"; 138 GmmSparse m1(rows,cols); 139 gmm::csr_matrix<Scalar> m2; 140 eiToGmm(sm1, m1); 141 gmm::copy(m1,m2); 142 std::vector<Scalar> gmmX(cols), gmmB(cols); 143 Map<Matrix<Scalar,Dynamic,1> >(&gmmX[0], cols) = x; 144 Map<Matrix<Scalar,Dynamic,1> >(&gmmB[0], cols) = b; 145 146 gmmX = gmmB; 147 BENCH(gmm::upper_tri_solve(m1, gmmX, false);) 148 std::cout << " colmajor^-1 * b:\t" << timer.value() << endl; 149 // std::cerr << Map<Matrix<Scalar,Dynamic,1> >(&gmmX[0], cols).transpose() << "\n"; 150 151 gmmX = gmmB; 152 BENCH(gmm::upper_tri_solve(m2, gmmX, false);) 153 timer.stop(); 154 std::cout << " rowmajor^-1 * b:\t" << timer.value() << endl; 155 // std::cerr << Map<Matrix<Scalar,Dynamic,1> >(&gmmX[0], cols).transpose() << "\n"; 156 } 157 #endif 158 159 // MTL4 160 #ifndef NOMTL 161 { 162 std::cout << "MTL4\t" << density*100 << "%\n"; 163 MtlSparse m1(rows,cols); 164 MtlSparseRowMajor m2(rows,cols); 165 eiToMtl(sm1, m1); 166 m2 = m1; 167 mtl::dense_vector<Scalar> x(rows, 1.0); 168 mtl::dense_vector<Scalar> b(rows, 1.0); 169 170 BENCH(x = mtl::upper_trisolve(m1,b);) 171 std::cout << " colmajor^-1 * b:\t" << timer.value() << endl; 172 // std::cerr << x << "\n"; 173 174 BENCH(x = mtl::upper_trisolve(m2,b);) 175 std::cout << " rowmajor^-1 * b:\t" << timer.value() << endl; 176 // std::cerr << x << "\n"; 177 } 178 #endif 179 180 181 std::cout << "\n\n"; 182 } 183 #endif 184 185 #if 0 186 // bench small matrices (in-place versus return bye value) 187 { 188 timer.reset(); 189 for (int _j=0; _j<10; ++_j) { 190 Matrix4f m = Matrix4f::Random(); 191 Vector4f b = Vector4f::Random(); 192 Vector4f x = Vector4f::Random(); 193 timer.start(); 194 for (int _k=0; _k<1000000; ++_k) { 195 b = m.inverseProduct(b); 196 } 197 timer.stop(); 198 } 199 std::cout << "4x4 :\t" << timer.value() << endl; 200 } 201 202 { 203 timer.reset(); 204 for (int _j=0; _j<10; ++_j) { 205 Matrix4f m = Matrix4f::Random(); 206 Vector4f b = Vector4f::Random(); 207 Vector4f x = Vector4f::Random(); 208 timer.start(); 209 for (int _k=0; _k<1000000; ++_k) { 210 m.inverseProductInPlace(x); 211 } 212 timer.stop(); 213 } 214 std::cout << "4x4 IP :\t" << timer.value() << endl; 215 } 216 #endif 217 218 return 0; 219 } 220 221
