ground/gcs/src/libs/eigen/bench/sparse_lu.cpp

   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