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 !

Eigen/Sparse

@@ -50,6 +50,10 @@
   namespace Eigen { struct SluMatrix; }
 #endif
 
+#ifdef EIGEN_UMFPACK_SUPPORT
+  #include "umfpack.h"
+#endif
+
 namespace Eigen {
 
 #include "src/Sparse/SparseUtil.h"
@@ -79,6 +83,10 @@ namespace Eigen {
 # include "src/Sparse/SuperLUSupport.h"
 #endif
 
+#ifdef EIGEN_UMFPACK_SUPPORT
+# include "src/Sparse/UmfPackSupport.h"
+#endif
+
 } // namespace Eigen
 
 #endif // EIGEN_SPARSE_MODULE_H

Eigen/src/Sparse/SuperLUSupport.h

@@ -311,8 +311,8 @@ void SparseLU<MatrixType,SuperLU>::compute(const MatrixType& a)
           &m_sluStat, &info, Scalar());
   StatFree(&m_sluStat);
 
-  // FIXME how to check for errors ???
+  // FIXME how to better check for errors ???
-  Base::m_succeeded = true;
+  Base::m_succeeded = (info == 0);
 }
 
 // template<typename MatrixType>
@@ -362,6 +362,8 @@ bool SparseLU<MatrixType,SuperLU>::solve(const MatrixBase<BDerived> &b, MatrixBa
     &m_sluFerr[0], &m_sluBerr[0],
     &m_sluStat, &info, Scalar());
   StatFree(&m_sluStat);
+
+  return info==0;
 }
 
 #endif // EIGEN_SUPERLUSUPPORT_H

Eigen/src/Sparse/UmfPackSupport.h

@@ -0,0 +1,135 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra. Eigen itself is part of the KDE project.
+//
+// Copyright (C) 2008 Gael Guennebaud <g.gael@free.fr>
+//
+// Eigen is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 3 of the License, or (at your option) any later version.
+//
+// Alternatively, you can redistribute it and/or
+// modify it under the terms of the GNU General Public License as
+// published by the Free Software Foundation; either version 2 of
+// the License, or (at your option) any later version.
+//
+// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
+// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License and a copy of the GNU General Public License along with
+// Eigen. If not, see <http://www.gnu.org/licenses/>.
+
+#ifndef EIGEN_UMFPACKSUPPORT_H
+#define EIGEN_UMFPACKSUPPORT_H
+
+template<typename MatrixType>
+class SparseLU<MatrixType,UmfPack> : public SparseLU<MatrixType>
+{
+  protected:
+    typedef SparseLU<MatrixType> Base;
+    typedef typename Base::Scalar Scalar;
+    typedef typename Base::RealScalar RealScalar;
+    typedef Matrix<Scalar,Dynamic,1> Vector;
+    using Base::m_flags;
+    using Base::m_status;
+
+  public:
+
+    SparseLU(int flags = NaturalOrdering)
+      : Base(flags), m_numeric(0)
+    {
+    }
+
+    SparseLU(const MatrixType& matrix, int flags = NaturalOrdering)
+      : Base(flags), m_numeric(0)
+    {
+      compute(matrix);
+    }
+
+    ~SparseLU()
+    {
+      if (m_numeric)
+        umfpack_di_free_numeric(&m_numeric);
+    }
+
+    template<typename BDerived, typename XDerived>
+    bool solve(const MatrixBase<BDerived> &b, MatrixBase<XDerived>* x) const;
+
+    void compute(const MatrixType& matrix);
+
+  protected:
+    // cached data:
+    void* m_numeric;
+    const MatrixType* m_matrixRef;
+};
+
+template<typename MatrixType>
+void SparseLU<MatrixType,UmfPack>::compute(const MatrixType& a)
+{
+  const int size = a.rows();
+  ei_assert((MatrixType::Flags&RowMajorBit)==0);
+
+  m_matrixRef = &a;
+
+  if (m_numeric)
+    umfpack_di_free_numeric(&m_numeric);
+
+  void* symbolic;
+  int errorCode = 0;
+  errorCode = umfpack_di_symbolic(size, size, a._outerIndexPtr(), a._innerIndexPtr(), a._valuePtr(),
+                                  &symbolic, 0, 0);
+  if (errorCode==0)
+    errorCode = umfpack_di_numeric(a._outerIndexPtr(), a._innerIndexPtr(), a._valuePtr(),
+                                   symbolic, &m_numeric, 0, 0);
+
+  umfpack_di_free_symbolic(&symbolic);
+
+  Base::m_succeeded = (errorCode==0);
+}
+
+// template<typename MatrixType>
+// inline const MatrixType&
+// SparseLU<MatrixType,SuperLU>::matrixL() const
+// {
+//   ei_assert(false && "matrixL() is Not supported by the SuperLU backend");
+//   return m_matrix;
+// }
+//
+// template<typename MatrixType>
+// inline const MatrixType&
+// SparseLU<MatrixType,SuperLU>::matrixU() const
+// {
+//   ei_assert(false && "matrixU() is Not supported by the SuperLU backend");
+//   return m_matrix;
+// }
+
+template<typename MatrixType>
+template<typename BDerived,typename XDerived>
+bool SparseLU<MatrixType,UmfPack>::solve(const MatrixBase<BDerived> &b, MatrixBase<XDerived> *x) const
+{
+  //const int size = m_matrix.rows();
+  const int rhsCols = b.cols();
+//   ei_assert(size==b.rows());
+  ei_assert((BDerived::Flags&RowMajorBit)==0 && "UmfPack backend does not support non col-major rhs yet");
+  ei_assert((XDerived::Flags&RowMajorBit)==0 && "UmfPack backend does not support non col-major result yet");
+
+  int errorCode;
+  for (int j=0; j<rhsCols; ++j)
+  {
+    errorCode = umfpack_di_solve(UMFPACK_A,
+        m_matrixRef->_outerIndexPtr(), m_matrixRef->_innerIndexPtr(), m_matrixRef->_valuePtr(),
+        &x->col(j).coeffRef(0), &b.const_cast_derived().col(j).coeffRef(0), m_numeric, 0, 0);
+    if (errorCode!=0)
+      return false;
+  }
+//   errorCode = umfpack_di_solve(UMFPACK_A,
+//       m_matrixRef._outerIndexPtr(), m_matrixRef._innerIndexPtr(), m_matrixRef._valuePtr(),
+//       x->derived().data(), b.derived().data(), m_numeric, 0, 0);
+
+  return true;
+}
+
+#endif // EIGEN_UMFPACKSUPPORT_H

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
+    #ifdef EIGEN_SUPERLU_SUPPORT
-    {
+    x.setZero();
-      x.setZero();
+    doEigen<Eigen::SuperLU>("Eigen/SuperLU (nat)", sm1, b, x, Eigen::NaturalOrdering);
-      BenchTimer timer;
+//     doEigen<Eigen::SuperLU>("Eigen/SuperLU (MD AT+A)", sm1, b, x, Eigen::MinimumDegree_AT_PLUS_A);
-      timer.start();
+//     doEigen<Eigen::SuperLU>("Eigen/SuperLU (MD ATA)", sm1, b, x, Eigen::MinimumDegree_ATA);
-      SparseLU<EigenSparseMatrix,SuperLU> lu(sm1);
+    doEigen<Eigen::SuperLU>("Eigen/SuperLU (COLAMD)", sm1, b, x, Eigen::ColApproxMinimumDegree);
-      timer.stop();
+    #endif
-      std::cout << "Eigen/SuperLU:\t" << timer.value() << endl;
 
-      timer.reset();
+    #ifdef EIGEN_UMFPACK_SUPPORT
-      timer.start();
+    x.setZero();
-      lu.solve(b,&x);
+    doEigen<Eigen::UmfPack>("Eigen/UmfPack (auto)", sm1, b, x, 0);
-      timer.stop();
+    #endif
-      std::cout << "  solve:\t" << timer.value() << endl;
-
-      std::cout << x.transpose() << "\n";
-
-    }
 
   }