1 2 // g++ -I.. sparse_lu.cpp -O3 -g0 -I /usr/include/superlu/ -lsuperlu -lgfortran -DSIZE=1000 -DDENSITY=.05 && ./a.out 3 4 #define EIGEN_SUPERLU_SUPPORT 5 #define EIGEN_UMFPACK_SUPPORT 6 #include <Eigen/Sparse> 7 8 #define NOGMM 9 #define NOMTL 10 11 #ifndef SIZE 12 #define SIZE 10 13 #endif 14 15 #ifndef DENSITY 16 #define DENSITY 0.01 17 #endif 18 19 #ifndef REPEAT 20 #define REPEAT 1 21 #endif 22 23 #include "BenchSparseUtil.h" 24 25 #ifndef MINDENSITY 26 #define MINDENSITY 0.0004 27 #endif 28 29 #ifndef NBTRIES 30 #define NBTRIES 10 31 #endif 32 33 #define BENCH(X) \ 34 timer.reset(); \ 35 for (int _j=0; _j<NBTRIES; ++_j) { \ 36 timer.start(); \ 37 for (int _k=0; _k<REPEAT; ++_k) { \ 38 X \ 39 } timer.stop(); } 40 41 typedef Matrix<Scalar,Dynamic,1> VectorX; 42 43 #include <Eigen/LU> 44 45 template<int Backend> 46 void doEigen(const char* name, const EigenSparseMatrix& sm1, const VectorX& b, VectorX& x, int flags = 0) 47 { 48 std::cout << name << "..." << std::flush; 49 BenchTimer timer; timer.start(); 50 SparseLU<EigenSparseMatrix,Backend> lu(sm1, flags); 51 timer.stop(); 52 if (lu.succeeded()) 53 std::cout << ":\t" << timer.value() << endl; 54 else 55 { 56 std::cout << ":\t FAILED" << endl; 57 return; 58 } 59 60 bool ok; 61 timer.reset(); timer.start(); 62 ok = lu.solve(b,&x); 63 timer.stop(); 64 if (ok) 65 std::cout << " solve:\t" << timer.value() << endl; 66 else 67 std::cout << " solve:\t" << " FAILED" << endl; 68 69 //std::cout << x.transpose() << "\n"; 70 } 71 72 int main(int argc, char *argv[]) 73 { 74 int rows = SIZE; 75 int cols = SIZE; 76 float density = DENSITY; 77 BenchTimer timer; 78 79 VectorX b = VectorX::Random(cols); 80 VectorX x = VectorX::Random(cols); 81 82 bool densedone = false; 83 84 //for (float density = DENSITY; density>=MINDENSITY; density*=0.5) 85 // float density = 0.5; 86 { 87 EigenSparseMatrix sm1(rows, cols); 88 fillMatrix(density, rows, cols, sm1); 89 90 // dense matrices 91 #ifdef DENSEMATRIX 92 if (!densedone) 93 { 94 densedone = true; 95 std::cout << "Eigen Dense\t" << density*100 << "%\n"; 96 DenseMatrix m1(rows,cols); 97 eiToDense(sm1, m1); 98 99 BenchTimer timer; 100 timer.start(); 101 FullPivLU<DenseMatrix> lu(m1); 102 timer.stop(); 103 std::cout << "Eigen/dense:\t" << timer.value() << endl; 104 105 timer.reset(); 106 timer.start(); 107 lu.solve(b,&x); 108 timer.stop(); 109 std::cout << " solve:\t" << timer.value() << endl; 110 // std::cout << b.transpose() << "\n"; 111 // std::cout << x.transpose() << "\n"; 112 } 113 #endif 114 115 #ifdef EIGEN_UMFPACK_SUPPORT 116 x.setZero(); 117 doEigen<Eigen::UmfPack>("Eigen/UmfPack (auto)", sm1, b, x, 0); 118 #endif 119 120 #ifdef EIGEN_SUPERLU_SUPPORT 121 x.setZero(); 122 doEigen<Eigen::SuperLU>("Eigen/SuperLU (nat)", sm1, b, x, Eigen::NaturalOrdering); 123 // doEigen<Eigen::SuperLU>("Eigen/SuperLU (MD AT+A)", sm1, b, x, Eigen::MinimumDegree_AT_PLUS_A); 124 // doEigen<Eigen::SuperLU>("Eigen/SuperLU (MD ATA)", sm1, b, x, Eigen::MinimumDegree_ATA); 125 doEigen<Eigen::SuperLU>("Eigen/SuperLU (COLAMD)", sm1, b, x, Eigen::ColApproxMinimumDegree); 126 #endif 127 128 } 129 130 return 0; 131 } 132 133