From 3a231c2349e2a6d62cfd600aff74cd2709613d48 Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Sun, 19 Oct 2008 22:44:21 +0000 Subject: 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 | 8 +++ Eigen/src/Sparse/SuperLUSupport.h | 6 +- Eigen/src/Sparse/UmfPackSupport.h | 135 ++++++++++++++++++++++++++++++++++++++ bench/sparse_lu.cpp | 60 +++++++++++------ 4 files changed, 187 insertions(+), 22 deletions(-) create mode 100644 Eigen/src/Sparse/UmfPackSupport.h diff --git a/Eigen/Sparse b/Eigen/Sparse index 0b222d105..df0310f89 100644 --- a/Eigen/Sparse +++ b/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 diff --git a/Eigen/src/Sparse/SuperLUSupport.h b/Eigen/src/Sparse/SuperLUSupport.h index fd2899f06..55e5ed366 100644 --- a/Eigen/src/Sparse/SuperLUSupport.h +++ b/Eigen/src/Sparse/SuperLUSupport.h @@ -311,8 +311,8 @@ void SparseLU::compute(const MatrixType& a) &m_sluStat, &info, Scalar()); StatFree(&m_sluStat); - // FIXME how to check for errors ??? - Base::m_succeeded = true; + // FIXME how to better check for errors ??? + Base::m_succeeded = (info == 0); } // template @@ -362,6 +362,8 @@ bool SparseLU::solve(const MatrixBase &b, MatrixBa &m_sluFerr[0], &m_sluBerr[0], &m_sluStat, &info, Scalar()); StatFree(&m_sluStat); + + return info==0; } #endif // EIGEN_SUPERLUSUPPORT_H diff --git a/Eigen/src/Sparse/UmfPackSupport.h b/Eigen/src/Sparse/UmfPackSupport.h new file mode 100644 index 000000000..d497e27d6 --- /dev/null +++ b/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 +// +// 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 . + +#ifndef EIGEN_UMFPACKSUPPORT_H +#define EIGEN_UMFPACKSUPPORT_H + +template +class SparseLU : public SparseLU +{ + protected: + typedef SparseLU Base; + typedef typename Base::Scalar Scalar; + typedef typename Base::RealScalar RealScalar; + typedef Matrix 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 + bool solve(const MatrixBase &b, MatrixBase* x) const; + + void compute(const MatrixType& matrix); + + protected: + // cached data: + void* m_numeric; + const MatrixType* m_matrixRef; +}; + +template +void SparseLU::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 +// inline const MatrixType& +// SparseLU::matrixL() const +// { +// ei_assert(false && "matrixL() is Not supported by the SuperLU backend"); +// return m_matrix; +// } +// +// template +// inline const MatrixType& +// SparseLU::matrixU() const +// { +// ei_assert(false && "matrixU() is Not supported by the SuperLU backend"); +// return m_matrix; +// } + +template +template +bool SparseLU::solve(const MatrixBase &b, MatrixBase *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_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 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 #define NOGMM @@ -43,6 +42,33 @@ typedef Matrix VectorX; #include +template +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 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 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 (nat)", sm1, b, x, Eigen::NaturalOrdering); +// doEigen("Eigen/SuperLU (MD AT+A)", sm1, b, x, Eigen::MinimumDegree_AT_PLUS_A); +// doEigen("Eigen/SuperLU (MD ATA)", sm1, b, x, Eigen::MinimumDegree_ATA); + doEigen("Eigen/SuperLU (COLAMD)", sm1, b, x, Eigen::ColApproxMinimumDegree); + #endif - } + #ifdef EIGEN_UMFPACK_SUPPORT + x.setZero(); + doEigen("Eigen/UmfPack (auto)", sm1, b, x, 0); + #endif } -- cgit v1.2.3