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