sparse module: add support for umfpack, the sparse direct LU

solver from suitesparse (as cholmod). It seems to be even faster
than SuperLU and it was much simpler to interface ! Well,
the factorization is faster, but for the solve part, SuperLU is
quite faster. On the other hand the solve part represents only a
fraction of the whole procedure. Moreover, I bench random matrices
that does not represents real cases, and I'm not sure at all
I use both libraries with their best settings !

1 files changed, 40 insertions, 20 deletions

diff --git a/bench/sparse_lu.cpp b/bench/sparse_lu.cpp
index 35e7e8a3a..d6020b552 100644
--- a/bench/sparse_lu.cpp
+++ b/bench/sparse_lu.cpp
@@ -1,9 +1,8 @@
 
 // g++ -I.. sparse_lu.cpp -O3 -g0 -I /usr/include/superlu/ -lsuperlu -lgfortran -DSIZE=1000 -DDENSITY=.05 && ./a.out
 
-// #define EIGEN_TAUCS_SUPPORT
-// #define EIGEN_CHOLMOD_SUPPORT
 #define EIGEN_SUPERLU_SUPPORT
+#define EIGEN_UMFPACK_SUPPORT
 #include <Eigen/Sparse>
 
 #define NOGMM
@@ -43,6 +42,33 @@ typedef Matrix<Scalar,Dynamic,1> VectorX;
 
 #include <Eigen/LU>
 
+template<int Backend>
+void doEigen(const char* name, const EigenSparseMatrix& sm1, const VectorX& b, VectorX& x, int flags = 0)
+{
+  std::cout << name << "..." << std::flush;
+  BenchTimer timer; timer.start();
+  SparseLU<EigenSparseMatrix,Backend> lu(sm1, flags);
+  timer.stop();
+  if (lu.succeeded())
+    std::cout << ":\t" << timer.value() << endl;
+  else
+  {
+    std::cout << ":\t FAILED" << endl;
+    return;
+  }
+
+  bool ok;
+  timer.reset(); timer.start();
+  ok = lu.solve(b,&x);
+  timer.stop();
+  if (ok)
+    std::cout << "  solve:\t" << timer.value() << endl;
+  else
+    std::cout << "  solve:\t" << " FAILED" << endl;
+
+  //std::cout << x.transpose() << "\n";
+}
+
 int main(int argc, char *argv[])
 {
   int rows = SIZE;
@@ -82,28 +108,22 @@ int main(int argc, char *argv[])
       timer.stop();
       std::cout << "  solve:\t" << timer.value() << endl;
 //       std::cout << b.transpose() << "\n";
-      std::cout << x.transpose() << "\n";
+//       std::cout << x.transpose() << "\n";
     }
     #endif
 
-    // eigen sparse matrices
-    {
-      x.setZero();
-      BenchTimer timer;
-      timer.start();
-      SparseLU<EigenSparseMatrix,SuperLU> lu(sm1);
-      timer.stop();
-      std::cout << "Eigen/SuperLU:\t" << timer.value() << endl;
-
-      timer.reset();
-      timer.start();
-      lu.solve(b,&x);
-      timer.stop();
-      std::cout << "  solve:\t" << timer.value() << endl;
-
-      std::cout << x.transpose() << "\n";
+    #ifdef EIGEN_SUPERLU_SUPPORT
+    x.setZero();
+    doEigen<Eigen::SuperLU>("Eigen/SuperLU (nat)", sm1, b, x, Eigen::NaturalOrdering);
+//     doEigen<Eigen::SuperLU>("Eigen/SuperLU (MD AT+A)", sm1, b, x, Eigen::MinimumDegree_AT_PLUS_A);
+//     doEigen<Eigen::SuperLU>("Eigen/SuperLU (MD ATA)", sm1, b, x, Eigen::MinimumDegree_ATA);
+    doEigen<Eigen::SuperLU>("Eigen/SuperLU (COLAMD)", sm1, b, x, Eigen::ColApproxMinimumDegree);
+    #endif
 
-    }
+    #ifdef EIGEN_UMFPACK_SUPPORT
+    x.setZero();
+    doEigen<Eigen::UmfPack>("Eigen/UmfPack (auto)", sm1, b, x, 0);
+    #endif
 
   }