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