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 41typedef Matrix<Scalar,Dynamic,1> VectorX; 42 43#include <Eigen/LU> 44 45template<int Backend> 46void 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 72int 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