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