Home | History | Annotate | Download | only in bench 
      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