From 2489c81562b2905384dc3bec6a30debe864f1c33 Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Thu, 7 Jul 2011 14:19:42 +0200 Subject: add new interface to SuperLU --- unsupported/Eigen/SuperLUSupport | 5 +- unsupported/Eigen/src/SparseExtra/SuperLUSupport.h | 679 +++++++++++++++------ .../Eigen/src/SparseExtra/SuperLUSupportLegacy.h | 404 ++++++++++++ unsupported/test/sparse_lu.cpp | 37 +- 4 files changed, 936 insertions(+), 189 deletions(-) create mode 100644 unsupported/Eigen/src/SparseExtra/SuperLUSupportLegacy.h diff --git a/unsupported/Eigen/SuperLUSupport b/unsupported/Eigen/SuperLUSupport index 89cb649b2..a0248dfb2 100644 --- a/unsupported/Eigen/SuperLUSupport +++ b/unsupported/Eigen/SuperLUSupport @@ -24,10 +24,11 @@ namespace Eigen { * \endcode */ -struct SuperLU {}; - #include "src/SparseExtra/SuperLUSupport.h" +struct SuperLULegacy {}; +#include "src/SparseExtra/SuperLUSupportLegacy.h" + } // namespace Eigen #include "../../Eigen/src/Core/util/ReenableStupidWarnings.h" diff --git a/unsupported/Eigen/src/SparseExtra/SuperLUSupport.h b/unsupported/Eigen/src/SparseExtra/SuperLUSupport.h index bb7312190..c4292c283 100644 --- a/unsupported/Eigen/src/SparseExtra/SuperLUSupport.h +++ b/unsupported/Eigen/src/SparseExtra/SuperLUSupport.h @@ -194,7 +194,7 @@ struct SluMatrix : SuperMatrix res.ncol = mat.cols(); } - res.Mtype = SLU_GE; + res.Mtype = SLU_GE; res.storage.nnz = mat.nonZeros(); res.storage.values = mat.derived()._valuePtr(); @@ -252,7 +252,7 @@ struct SluMatrixMapHelper > res.ncol = mat.cols(); } - res.Mtype = SLU_GE; + res.Mtype = SLU_GE; res.storage.nnz = mat.nonZeros(); res.storage.values = mat._valuePtr(); @@ -295,83 +295,146 @@ MappedSparseMatrix map_superlu(SluMatrix& sluMat) } // end namespace internal -template -class SparseLU : public SparseLU + +template +class SuperLUBase { - protected: - typedef SparseLU Base; - typedef typename Base::Scalar Scalar; - typedef typename Base::RealScalar RealScalar; + public: + typedef _MatrixType MatrixType; + typedef typename MatrixType::Scalar Scalar; + typedef typename MatrixType::RealScalar RealScalar; + typedef typename MatrixType::Index Index; typedef Matrix Vector; typedef Matrix IntRowVectorType; - typedef Matrix IntColVectorType; - typedef SparseMatrix LMatrixType; - typedef SparseMatrix UMatrixType; - using Base::m_flags; - using Base::m_status; + typedef Matrix IntColVectorType; + typedef SparseMatrix LUMatrixType; public: - SparseLU(int flags = NaturalOrdering) - : Base(flags) - { - } + SuperLUBase() {} - SparseLU(const MatrixType& matrix, int flags = NaturalOrdering) - : Base(flags) - { - compute(matrix); - } - - ~SparseLU() + ~SuperLUBase() { Destroy_SuperNode_Matrix(&m_sluL); Destroy_CompCol_Matrix(&m_sluU); } - - inline const LMatrixType& matrixL() const + + Derived& derived() { return *static_cast(this); } + const Derived& derived() const { return *static_cast(this); } + + inline Index rows() const { return m_matrix.rows(); } + inline Index cols() const { return m_matrix.cols(); } + + /** \returns a reference to the Super LU option object to configure the Super LU algorithms. */ + inline superlu_options_t& options() { return m_sluOptions; } + + /** \brief Reports whether previous computation was successful. + * + * \returns \c Success if computation was succesful, + * \c NumericalIssue if the matrix.appears to be negative. + */ + ComputationInfo info() const { - if (m_extractedDataAreDirty) extractData(); - return m_l; + eigen_assert(m_isInitialized && "Decomposition is not initialized."); + return m_info; } - inline const UMatrixType& matrixU() const + /** Computes the sparse Cholesky decomposition of \a matrix */ + void compute(const MatrixType& matrix) { - if (m_extractedDataAreDirty) extractData(); - return m_u; + derived().analyzePattern(matrix); + derived().factorize(matrix); } - - inline const IntColVectorType& permutationP() const + + /** \returns the solution x of \f$ A x = b \f$ using the current decomposition of A. + * + * \sa compute() + */ + template + inline const internal::solve_retval solve(const MatrixBase& b) const { - if (m_extractedDataAreDirty) extractData(); - return m_p; + eigen_assert(m_isInitialized && "SuperLU is not initialized."); + eigen_assert(rows()==b.rows() + && "SuperLU::solve(): invalid number of rows of the right hand side matrix b"); + return internal::solve_retval(*this, b.derived()); } - - inline const IntRowVectorType& permutationQ() const + + /** \returns the solution x of \f$ A x = b \f$ using the current decomposition of A. + * + * \sa compute() + */ +// template +// inline const internal::sparse_solve_retval solve(const SparseMatrixBase& b) const +// { +// eigen_assert(m_isInitialized && "SuperLU is not initialized."); +// eigen_assert(rows()==b.rows() +// && "SuperLU::solve(): invalid number of rows of the right hand side matrix b"); +// return internal::sparse_solve_retval(*this, b.derived()); +// } + + /** Performs a symbolic decomposition on the sparcity of \a matrix. + * + * This function is particularly useful when solving for several problems having the same structure. + * + * \sa factorize() + */ + void analyzePattern(const MatrixType& /*matrix*/) { - if (m_extractedDataAreDirty) extractData(); - return m_q; + m_isInitialized = true; + m_info = Success; + m_analysisIsOk = true; + m_factorizationIsOk = false; } - - Scalar determinant() const; - - template - bool solve(const MatrixBase &b, MatrixBase* x, const int transposed = SvNoTrans) const; - - void compute(const MatrixType& matrix); - + + template + void dumpMemory(Stream& s) + {} + protected: - + + void initFactorization(const MatrixType& a) + { + const int size = a.rows(); + m_matrix = a; + + m_sluA = internal::asSluMatrix(m_matrix); + memset(&m_sluL,0,sizeof m_sluL); + memset(&m_sluU,0,sizeof m_sluU); + + m_p.resize(size); + m_q.resize(size); + m_sluRscale.resize(size); + m_sluCscale.resize(size); + m_sluEtree.resize(size); + + // set empty B and X + m_sluB.setStorageType(SLU_DN); + m_sluB.setScalarType(); + m_sluB.Mtype = SLU_GE; + m_sluB.storage.values = 0; + m_sluB.nrow = 0; + m_sluB.ncol = 0; + m_sluB.storage.lda = size; + m_sluX = m_sluB; + + m_extractedDataAreDirty = true; + } + + void init() + { + m_info = InvalidInput; + m_isInitialized = false; + } + void extractData() const; - protected: // cached data to reduce reallocation, etc. - mutable LMatrixType m_l; - mutable UMatrixType m_u; + mutable LUMatrixType m_l; + mutable LUMatrixType m_u; mutable IntColVectorType m_p; mutable IntRowVectorType m_q; - mutable SparseMatrix m_matrix; + mutable LUMatrixType m_matrix; // copy of the factorized matrix mutable SluMatrix m_sluA; mutable SuperMatrix m_sluL, m_sluU; mutable SluMatrix m_sluB, m_sluX; @@ -381,171 +444,212 @@ class SparseLU : public SparseLU mutable std::vector m_sluRscale, m_sluCscale; mutable std::vector m_sluFerr, m_sluBerr; mutable char m_sluEqued; + + mutable ComputationInfo m_info; + bool m_isInitialized; + int m_factorizationIsOk; + int m_analysisIsOk; mutable bool m_extractedDataAreDirty; }; -template -void SparseLU::compute(const MatrixType& a) + +template +class SuperLU : public SuperLUBase<_MatrixType,SuperLU<_MatrixType> > { - const int size = a.rows(); - m_matrix = a; + public: + typedef SuperLUBase<_MatrixType,SuperLU> Base; + typedef _MatrixType MatrixType; + typedef typename Base::Scalar Scalar; + typedef typename Base::RealScalar RealScalar; + typedef typename Base::Index Index; + typedef typename Base::IntRowVectorType IntRowVectorType; + typedef typename Base::IntColVectorType IntColVectorType; + typedef typename Base::LUMatrixType LUMatrixType; + typedef TriangularView LMatrixType; + typedef TriangularView UMatrixType; - set_default_options(&m_sluOptions); - m_sluOptions.ColPerm = NATURAL; - m_sluOptions.PrintStat = NO; - m_sluOptions.ConditionNumber = NO; - m_sluOptions.Trans = NOTRANS; - // m_sluOptions.Equil = NO; + public: - switch (Base::orderingMethod()) - { - case NaturalOrdering : m_sluOptions.ColPerm = NATURAL; break; - case MinimumDegree_AT_PLUS_A : m_sluOptions.ColPerm = MMD_AT_PLUS_A; break; - case MinimumDegree_ATA : m_sluOptions.ColPerm = MMD_ATA; break; - case ColApproxMinimumDegree : m_sluOptions.ColPerm = COLAMD; break; - default: - //std::cerr << "Eigen: ordering method \"" << Base::orderingMethod() << "\" not supported by the SuperLU backend\n"; - m_sluOptions.ColPerm = NATURAL; - }; - - m_sluA = internal::asSluMatrix(m_matrix); - memset(&m_sluL,0,sizeof m_sluL); - memset(&m_sluU,0,sizeof m_sluU); - //m_sluEqued = 'B'; - int info = 0; + SuperLU() : Base() { init(); } - m_p.resize(size); - m_q.resize(size); - m_sluRscale.resize(size); - m_sluCscale.resize(size); - m_sluEtree.resize(size); + SuperLU(const MatrixType& matrix) : Base() + { + init(); + compute(matrix); + } - RealScalar recip_pivot_gross, rcond; - RealScalar ferr, berr; + ~SuperLU() + { + } + + /** Performs a symbolic decomposition on the sparcity of \a matrix. + * + * This function is particularly useful when solving for several problems having the same structure. + * + * \sa factorize() + */ + void analyzePattern(const MatrixType& matrix) + { + Base::analyzePattern(matrix); + } + + /** Performs a numeric decomposition of \a matrix + * + * The given matrix must has the same sparcity than the matrix on which the symbolic decomposition has been performed. + * + * \sa analyzePattern() + */ + void factorize(const MatrixType& matrix); + + #ifndef EIGEN_PARSED_BY_DOXYGEN + /** \internal */ + template + void _solve(const MatrixBase &b, MatrixBase &dest) const; + #endif // EIGEN_PARSED_BY_DOXYGEN + + inline const LMatrixType& matrixL() const + { + if (m_extractedDataAreDirty) this->extractData(); + return m_l; + } - // set empty B and X - m_sluB.setStorageType(SLU_DN); - m_sluB.setScalarType(); - m_sluB.Mtype = SLU_GE; - m_sluB.storage.values = 0; - m_sluB.nrow = m_sluB.ncol = 0; - m_sluB.storage.lda = size; - m_sluX = m_sluB; + inline const UMatrixType& matrixU() const + { + if (m_extractedDataAreDirty) this->extractData(); + return m_u; + } - StatInit(&m_sluStat); - if (m_flags&IncompleteFactorization) + inline const IntColVectorType& permutationP() const + { + if (m_extractedDataAreDirty) this->extractData(); + return m_p; + } + + inline const IntRowVectorType& permutationQ() const + { + if (m_extractedDataAreDirty) this->extractData(); + return m_q; + } + + Scalar determinant() const; + + protected: + + using Base::m_matrix; + using Base::m_sluOptions; + using Base::m_sluA; + using Base::m_sluB; + using Base::m_sluX; + using Base::m_p; + using Base::m_q; + using Base::m_sluEtree; + using Base::m_sluEqued; + using Base::m_sluRscale; + using Base::m_sluCscale; + using Base::m_sluL; + using Base::m_sluU; + using Base::m_sluStat; + using Base::m_sluFerr; + using Base::m_sluBerr; + using Base::m_l; + using Base::m_u; + + using Base::m_analysisIsOk; + using Base::m_factorizationIsOk; + using Base::m_extractedDataAreDirty; + using Base::m_isInitialized; + using Base::m_info; + + void init() + { + Base::init(); + + set_default_options(&this->m_sluOptions); + m_sluOptions.PrintStat = NO; + m_sluOptions.ConditionNumber = NO; + m_sluOptions.Trans = NOTRANS; + m_sluOptions.ColPerm = COLAMD; + } +}; + +template +void SuperLU::factorize(const MatrixType& a) +{ + eigen_assert(m_analysisIsOk && "You must first call analyzePattern()"); + if(!m_analysisIsOk) { - #ifdef EIGEN_SUPERLU_HAS_ILU - ilu_set_default_options(&m_sluOptions); - - // no attempt to preserve column sum - m_sluOptions.ILU_MILU = SILU; - - // only basic ILU(k) support -- no direct control over memory consumption - // better to use ILU_DropRule = DROP_BASIC | DROP_AREA - // and set ILU_FillFactor to max memory growth - m_sluOptions.ILU_DropRule = DROP_BASIC; - m_sluOptions.ILU_DropTol = Base::m_precision; - - SuperLU_gsisx(&m_sluOptions, &m_sluA, m_q.data(), m_p.data(), &m_sluEtree[0], - &m_sluEqued, &m_sluRscale[0], &m_sluCscale[0], - &m_sluL, &m_sluU, - NULL, 0, - &m_sluB, &m_sluX, - &recip_pivot_gross, &rcond, - &m_sluStat, &info, Scalar()); - #else - //std::cerr << "Incomplete factorization is only available in SuperLU v4\n"; - Base::m_succeeded = false; + m_info = InvalidInput; return; - #endif - } - else - { - SuperLU_gssvx(&m_sluOptions, &m_sluA, m_q.data(), m_p.data(), &m_sluEtree[0], - &m_sluEqued, &m_sluRscale[0], &m_sluCscale[0], - &m_sluL, &m_sluU, - NULL, 0, - &m_sluB, &m_sluX, - &recip_pivot_gross, &rcond, - &ferr, &berr, - &m_sluStat, &info, Scalar()); } + + initFactorization(a); + + int info = 0; + RealScalar recip_pivot_growth, rcond; + RealScalar ferr, berr; + + StatInit(&m_sluStat); + SuperLU_gssvx(&m_sluOptions, &m_sluA, m_q.data(), m_p.data(), &m_sluEtree[0], + &m_sluEqued, &m_sluRscale[0], &m_sluCscale[0], + &m_sluL, &m_sluU, + NULL, 0, + &m_sluB, &m_sluX, + &recip_pivot_growth, &rcond, + &ferr, &berr, + &m_sluStat, &info, Scalar()); StatFree(&m_sluStat); m_extractedDataAreDirty = true; // FIXME how to better check for errors ??? - Base::m_succeeded = (info == 0); + m_info = info == 0 ? Success : NumericalIssue; + m_factorizationIsOk = true; } template -template -bool SparseLU::solve(const MatrixBase &b, - MatrixBase *x, const int transposed) const +template +void SuperLU::_solve(const MatrixBase &b, MatrixBase& x) const { + eigen_assert(m_factorizationIsOk && "The decomposition is not in a valid state for solving, you must first call either compute() or analyzePattern()/factorize()"); + const int size = m_matrix.rows(); const int rhsCols = b.cols(); eigen_assert(size==b.rows()); - switch (transposed) { - case SvNoTrans : m_sluOptions.Trans = NOTRANS; break; - case SvTranspose : m_sluOptions.Trans = TRANS; break; - case SvAdjoint : m_sluOptions.Trans = CONJ; break; - default: - //std::cerr << "Eigen: transposition option \"" << transposed << "\" not supported by the SuperLU backend\n"; - m_sluOptions.Trans = NOTRANS; - } - + m_sluOptions.Trans = NOTRANS; m_sluOptions.Fact = FACTORED; m_sluOptions.IterRefine = NOREFINE; + m_sluFerr.resize(rhsCols); m_sluBerr.resize(rhsCols); m_sluB = SluMatrix::Map(b.const_cast_derived()); - m_sluX = SluMatrix::Map(x->derived()); + m_sluX = SluMatrix::Map(x.derived()); + + typename Rhs::PlainObject b_cpy; + if(m_sluEqued!='N') + { + b_cpy = b; + m_sluB = SluMatrix::Map(b_cpy.const_cast_derived()); + } StatInit(&m_sluStat); int info = 0; - RealScalar recip_pivot_gross, rcond; - - if (m_flags&IncompleteFactorization) - { - #ifdef EIGEN_SUPERLU_HAS_ILU - SuperLU_gsisx(&m_sluOptions, &m_sluA, m_q.data(), m_p.data(), &m_sluEtree[0], - &m_sluEqued, &m_sluRscale[0], &m_sluCscale[0], - &m_sluL, &m_sluU, - NULL, 0, - &m_sluB, &m_sluX, - &recip_pivot_gross, &rcond, - &m_sluStat, &info, Scalar()); - #else - //std::cerr << "Incomplete factorization is only available in SuperLU v4\n"; - return false; - #endif - } - else - { - SuperLU_gssvx( - &m_sluOptions, &m_sluA, - m_q.data(), m_p.data(), - &m_sluEtree[0], &m_sluEqued, - &m_sluRscale[0], &m_sluCscale[0], - &m_sluL, &m_sluU, - NULL, 0, - &m_sluB, &m_sluX, - &recip_pivot_gross, &rcond, - &m_sluFerr[0], &m_sluBerr[0], - &m_sluStat, &info, Scalar()); - } + RealScalar recip_pivot_growth, rcond; + SuperLU_gssvx(&m_sluOptions, &m_sluA, + m_q.data(), m_p.data(), + &m_sluEtree[0], &m_sluEqued, + &m_sluRscale[0], &m_sluCscale[0], + &m_sluL, &m_sluU, + NULL, 0, + &m_sluB, &m_sluX, + &recip_pivot_growth, &rcond, + &m_sluFerr[0], &m_sluBerr[0], + &m_sluStat, &info, Scalar()); StatFree(&m_sluStat); - - // reset to previous state - m_sluOptions.Trans = NOTRANS; - return info==0; + m_info = info==0 ? Success : NumericalIssue; } -// // the code of this extractData() function has been adapted from the SuperLU's Matlab support code, // // Copyright (c) 1994 by Xerox Corporation. All rights reserved. @@ -553,9 +657,10 @@ bool SparseLU::solve(const MatrixBase &b, // THIS MATERIAL IS PROVIDED AS IS, WITH ABSOLUTELY NO WARRANTY // EXPRESSED OR IMPLIED. ANY USE IS AT YOUR OWN RISK. // -template -void SparseLU::extractData() const +template +void SuperLUBase::extractData() const { + eigen_assert(m_factorizationIsOk && "The decomposition is not in a valid state for extracting factors, you must first call either compute() or analyzePattern()/factorize()"); if (m_extractedDataAreDirty) { int upper; @@ -639,13 +744,14 @@ void SparseLU::extractData() const } template -typename SparseLU::Scalar SparseLU::determinant() const +typename SuperLU::Scalar SuperLU::determinant() const { - assert((!NumTraits::IsComplex) && "This function is not implemented for complex yet"); + eigen_assert((!NumTraits::IsComplex) && "This function is not implemented for complex yet"); + eigen_assert(m_factorizationIsOk && "The decomposition is not in a valid state for computing the determinant, you must first call either compute() or analyzePattern()/factorize()"); + if (m_extractedDataAreDirty) - extractData(); + this->extractData(); - // TODO this code could be moved to the default/base backend // FIXME perhaps we have to take into account the scale factors m_sluRscale and m_sluCscale ??? Scalar det = Scalar(1); for (int j=0; j::Scalar SparseLU::dete det *= m_u._valuePtr()[lastId]; } } -// std::cout << m_sluRscale[j] << " " << m_sluCscale[j] << " \n"; } return det; } +#ifdef EIGEN_SUPERLU_HAS_ILU +template +class SuperILU : public SuperLUBase<_MatrixType,SuperILU<_MatrixType> > +{ + public: + typedef SuperLUBase<_MatrixType,SuperILU> Base; + typedef _MatrixType MatrixType; + typedef typename Base::Scalar Scalar; + typedef typename Base::RealScalar RealScalar; + typedef typename Base::Index Index; + + public: + + SuperILU() : Base() { init(); } + + SuperILU(const MatrixType& matrix) : Base() + { + init(); + compute(matrix); + } + + ~SuperILU() + { + } + + /** Performs a symbolic decomposition on the sparcity of \a matrix. + * + * This function is particularly useful when solving for several problems having the same structure. + * + * \sa factorize() + */ + void analyzePattern(const MatrixType& matrix) + { + Base::analyzePattern(matrix); + } + + /** Performs a numeric decomposition of \a matrix + * + * The given matrix must has the same sparcity than the matrix on which the symbolic decomposition has been performed. + * + * \sa analyzePattern() + */ + void factorize(const MatrixType& matrix); + + #ifndef EIGEN_PARSED_BY_DOXYGEN + /** \internal */ + template + void _solve(const MatrixBase &b, MatrixBase &dest) const; + #endif // EIGEN_PARSED_BY_DOXYGEN + + protected: + + using Base::m_matrix; + using Base::m_sluOptions; + using Base::m_sluA; + using Base::m_sluB; + using Base::m_sluX; + using Base::m_p; + using Base::m_q; + using Base::m_sluEtree; + using Base::m_sluEqued; + using Base::m_sluRscale; + using Base::m_sluCscale; + using Base::m_sluL; + using Base::m_sluU; + using Base::m_sluStat; + using Base::m_sluFerr; + using Base::m_sluBerr; + using Base::m_l; + using Base::m_u; + + using Base::m_analysisIsOk; + using Base::m_factorizationIsOk; + using Base::m_extractedDataAreDirty; + using Base::m_isInitialized; + using Base::m_info; + + void init() + { + Base::init(); + + ilu_set_default_options(&m_sluOptions); + m_sluOptions.PrintStat = NO; + m_sluOptions.ConditionNumber = NO; + m_sluOptions.Trans = NOTRANS; + m_sluOptions.ColPerm = MMD_AT_PLUS_A; + + // no attempt to preserve column sum + m_sluOptions.ILU_MILU = SILU; + // only basic ILU(k) support -- no direct control over memory consumption + // better to use ILU_DropRule = DROP_BASIC | DROP_AREA + // and set ILU_FillFactor to max memory growth + m_sluOptions.ILU_DropRule = DROP_BASIC; + m_sluOptions.ILU_DropTol = NumTraits::dummy_precision()*10; + } +}; + +template +void SuperILU::factorize(const MatrixType& a) +{ + eigen_assert(m_analysisIsOk && "You must first call analyzePattern()"); + if(!m_analysisIsOk) + { + m_info = InvalidInput; + return; + } + + this->initFactorization(a); + + int info = 0; + RealScalar recip_pivot_growth, rcond; + + StatInit(&m_sluStat); + SuperLU_gsisx(&m_sluOptions, &m_sluA, m_q.data(), m_p.data(), &m_sluEtree[0], + &m_sluEqued, &m_sluRscale[0], &m_sluCscale[0], + &m_sluL, &m_sluU, + NULL, 0, + &m_sluB, &m_sluX, + &recip_pivot_growth, &rcond, + &m_sluStat, &info, Scalar()); + StatFree(&m_sluStat); + + // FIXME how to better check for errors ??? + m_info = info == 0 ? Success : NumericalIssue; + m_factorizationIsOk = true; +} + +template +template +void SuperILU::_solve(const MatrixBase &b, MatrixBase& x) const +{ + eigen_assert(m_factorizationIsOk && "The decomposition is not in a valid state for solving, you must first call either compute() or analyzePattern()/factorize()"); + + const int size = m_matrix.rows(); + const int rhsCols = b.cols(); + eigen_assert(size==b.rows()); + + m_sluOptions.Trans = NOTRANS; + m_sluOptions.Fact = FACTORED; + m_sluOptions.IterRefine = NOREFINE; + + m_sluFerr.resize(rhsCols); + m_sluBerr.resize(rhsCols); + m_sluB = SluMatrix::Map(b.const_cast_derived()); + m_sluX = SluMatrix::Map(x.derived()); + + typename Rhs::PlainObject b_cpy; + if(m_sluEqued!='N') + { + b_cpy = b; + m_sluB = SluMatrix::Map(b_cpy.const_cast_derived()); + } + + int info = 0; + RealScalar recip_pivot_growth, rcond; + + StatInit(&m_sluStat); + SuperLU_gsisx(&m_sluOptions, &m_sluA, + m_q.data(), m_p.data(), + &m_sluEtree[0], &m_sluEqued, + &m_sluRscale[0], &m_sluCscale[0], + &m_sluL, &m_sluU, + NULL, 0, + &m_sluB, &m_sluX, + &recip_pivot_growth, &rcond, + &m_sluStat, &info, Scalar()); + StatFree(&m_sluStat); + + m_info = info==0 ? Success : NumericalIssue; +} +#endif + +namespace internal { + +template +struct solve_retval, Rhs> + : solve_retval_base, Rhs> +{ + typedef SuperLUBase<_MatrixType,Derived> Dec; + EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs) + + template void evalTo(Dest& dst) const + { + dec().derived()._solve(rhs(),dst); + } +}; + +template +struct sparse_solve_retval, Rhs> + : sparse_solve_retval_base, Rhs> +{ + typedef SuperLUBase<_MatrixType,Derived> Dec; + EIGEN_MAKE_SPARSE_SOLVE_HELPERS(Dec,Rhs) + + template void evalTo(Dest& dst) const + { + dec().derived()._solve(rhs(),dst); + } +}; + +} + #endif // EIGEN_SUPERLUSUPPORT_H diff --git a/unsupported/Eigen/src/SparseExtra/SuperLUSupportLegacy.h b/unsupported/Eigen/src/SparseExtra/SuperLUSupportLegacy.h new file mode 100644 index 000000000..e01aba4e6 --- /dev/null +++ b/unsupported/Eigen/src/SparseExtra/SuperLUSupportLegacy.h @@ -0,0 +1,404 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2008-2009 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_SUPERLUSUPPORT_LEGACY_H +#define EIGEN_SUPERLUSUPPORT_LEGACY_H + +template +class SparseLU : public SparseLU +{ + protected: + typedef SparseLU Base; + typedef typename Base::Scalar Scalar; + typedef typename Base::RealScalar RealScalar; + typedef Matrix Vector; + typedef Matrix IntRowVectorType; + typedef Matrix IntColVectorType; + typedef SparseMatrix LMatrixType; + typedef SparseMatrix UMatrixType; + using Base::m_flags; + using Base::m_status; + + public: + + SparseLU(int flags = NaturalOrdering) + : Base(flags) + { + } + + SparseLU(const MatrixType& matrix, int flags = NaturalOrdering) + : Base(flags) + { + compute(matrix); + } + + ~SparseLU() + { + Destroy_SuperNode_Matrix(&m_sluL); + Destroy_CompCol_Matrix(&m_sluU); + } + + inline const LMatrixType& matrixL() const + { + if (m_extractedDataAreDirty) extractData(); + return m_l; + } + + inline const UMatrixType& matrixU() const + { + if (m_extractedDataAreDirty) extractData(); + return m_u; + } + + inline const IntColVectorType& permutationP() const + { + if (m_extractedDataAreDirty) extractData(); + return m_p; + } + + inline const IntRowVectorType& permutationQ() const + { + if (m_extractedDataAreDirty) extractData(); + return m_q; + } + + Scalar determinant() const; + + template + bool solve(const MatrixBase &b, MatrixBase* x, const int transposed = SvNoTrans) const; + + void compute(const MatrixType& matrix); + + protected: + + void extractData() const; + + protected: + // cached data to reduce reallocation, etc. + mutable LMatrixType m_l; + mutable UMatrixType m_u; + mutable IntColVectorType m_p; + mutable IntRowVectorType m_q; + + mutable SparseMatrix m_matrix; + mutable SluMatrix m_sluA; + mutable SuperMatrix m_sluL, m_sluU; + mutable SluMatrix m_sluB, m_sluX; + mutable SuperLUStat_t m_sluStat; + mutable superlu_options_t m_sluOptions; + mutable std::vector m_sluEtree; + mutable std::vector m_sluRscale, m_sluCscale; + mutable std::vector m_sluFerr, m_sluBerr; + mutable char m_sluEqued; + mutable bool m_extractedDataAreDirty; +}; + +template +void SparseLU::compute(const MatrixType& a) +{ + const int size = a.rows(); + m_matrix = a; + + set_default_options(&m_sluOptions); + m_sluOptions.ColPerm = NATURAL; + m_sluOptions.PrintStat = NO; + m_sluOptions.ConditionNumber = NO; + m_sluOptions.Trans = NOTRANS; + // m_sluOptions.Equil = NO; + + switch (Base::orderingMethod()) + { + case NaturalOrdering : m_sluOptions.ColPerm = NATURAL; break; + case MinimumDegree_AT_PLUS_A : m_sluOptions.ColPerm = MMD_AT_PLUS_A; break; + case MinimumDegree_ATA : m_sluOptions.ColPerm = MMD_ATA; break; + case ColApproxMinimumDegree : m_sluOptions.ColPerm = COLAMD; break; + default: + //std::cerr << "Eigen: ordering method \"" << Base::orderingMethod() << "\" not supported by the SuperLU backend\n"; + m_sluOptions.ColPerm = NATURAL; + }; + + m_sluA = internal::asSluMatrix(m_matrix); + memset(&m_sluL,0,sizeof m_sluL); + memset(&m_sluU,0,sizeof m_sluU); + m_sluEqued = 'N'; + int info = 0; + + m_p.resize(size); + m_q.resize(size); + m_sluRscale.resize(size); + m_sluCscale.resize(size); + m_sluEtree.resize(size); + + RealScalar recip_pivot_gross, rcond; + RealScalar ferr, berr; + + // set empty B and X + m_sluB.setStorageType(SLU_DN); + m_sluB.setScalarType(); + m_sluB.Mtype = SLU_GE; + m_sluB.storage.values = 0; + m_sluB.nrow = m_sluB.ncol = 0; + m_sluB.storage.lda = size; + m_sluX = m_sluB; + + StatInit(&m_sluStat); + if (m_flags&IncompleteFactorization) + { + #ifdef EIGEN_SUPERLU_HAS_ILU + ilu_set_default_options(&m_sluOptions); + + // no attempt to preserve column sum + m_sluOptions.ILU_MILU = SILU; + + // only basic ILU(k) support -- no direct control over memory consumption + // better to use ILU_DropRule = DROP_BASIC | DROP_AREA + // and set ILU_FillFactor to max memory growth + m_sluOptions.ILU_DropRule = DROP_BASIC; + m_sluOptions.ILU_DropTol = Base::m_precision; + + SuperLU_gsisx(&m_sluOptions, &m_sluA, m_q.data(), m_p.data(), &m_sluEtree[0], + &m_sluEqued, &m_sluRscale[0], &m_sluCscale[0], + &m_sluL, &m_sluU, + NULL, 0, + &m_sluB, &m_sluX, + &recip_pivot_gross, &rcond, + &m_sluStat, &info, Scalar()); + #else + //std::cerr << "Incomplete factorization is only available in SuperLU v4\n"; + Base::m_succeeded = false; + return; + #endif + } + else + { + SuperLU_gssvx(&m_sluOptions, &m_sluA, m_q.data(), m_p.data(), &m_sluEtree[0], + &m_sluEqued, &m_sluRscale[0], &m_sluCscale[0], + &m_sluL, &m_sluU, + NULL, 0, + &m_sluB, &m_sluX, + &recip_pivot_gross, &rcond, + &ferr, &berr, + &m_sluStat, &info, Scalar()); + } + StatFree(&m_sluStat); + + m_extractedDataAreDirty = true; + + // FIXME how to better check for errors ??? + Base::m_succeeded = (info == 0); +} + +template +template +bool SparseLU::solve(const MatrixBase &b, + MatrixBase *x, const int transposed) const +{ + const int size = m_matrix.rows(); + const int rhsCols = b.cols(); + eigen_assert(size==b.rows()); + + switch (transposed) { + case SvNoTrans : m_sluOptions.Trans = NOTRANS; break; + case SvTranspose : m_sluOptions.Trans = TRANS; break; + case SvAdjoint : m_sluOptions.Trans = CONJ; break; + default: + //std::cerr << "Eigen: transposition option \"" << transposed << "\" not supported by the SuperLU backend\n"; + m_sluOptions.Trans = NOTRANS; + } + + m_sluOptions.Fact = FACTORED; + m_sluOptions.IterRefine = NOREFINE; + + m_sluFerr.resize(rhsCols); + m_sluBerr.resize(rhsCols); + m_sluB = SluMatrix::Map(b.const_cast_derived()); + m_sluX = SluMatrix::Map(x->derived()); + + typename BDerived::PlainObject b_cpy; + if(m_sluEqued!='N') + { + b_cpy = b; + m_sluB = SluMatrix::Map(b_cpy.const_cast_derived()); + } + + StatInit(&m_sluStat); + int info = 0; + RealScalar recip_pivot_gross, rcond; + + if (m_flags&IncompleteFactorization) + { + #ifdef EIGEN_SUPERLU_HAS_ILU + SuperLU_gsisx(&m_sluOptions, &m_sluA, m_q.data(), m_p.data(), &m_sluEtree[0], + &m_sluEqued, &m_sluRscale[0], &m_sluCscale[0], + &m_sluL, &m_sluU, + NULL, 0, + &m_sluB, &m_sluX, + &recip_pivot_gross, &rcond, + &m_sluStat, &info, Scalar()); + #else + //std::cerr << "Incomplete factorization is only available in SuperLU v4\n"; + return false; + #endif + } + else + { + SuperLU_gssvx( + &m_sluOptions, &m_sluA, + m_q.data(), m_p.data(), + &m_sluEtree[0], &m_sluEqued, + &m_sluRscale[0], &m_sluCscale[0], + &m_sluL, &m_sluU, + NULL, 0, + &m_sluB, &m_sluX, + &recip_pivot_gross, &rcond, + &m_sluFerr[0], &m_sluBerr[0], + &m_sluStat, &info, Scalar()); + } + StatFree(&m_sluStat); + + // reset to previous state + m_sluOptions.Trans = NOTRANS; + return info==0; +} + +// +// the code of this extractData() function has been adapted from the SuperLU's Matlab support code, +// +// Copyright (c) 1994 by Xerox Corporation. All rights reserved. +// +// THIS MATERIAL IS PROVIDED AS IS, WITH ABSOLUTELY NO WARRANTY +// EXPRESSED OR IMPLIED. ANY USE IS AT YOUR OWN RISK. +// +template +void SparseLU::extractData() const +{ + if (m_extractedDataAreDirty) + { + int upper; + int fsupc, istart, nsupr; + int lastl = 0, lastu = 0; + SCformat *Lstore = static_cast(m_sluL.Store); + NCformat *Ustore = static_cast(m_sluU.Store); + Scalar *SNptr; + + const int size = m_matrix.rows(); + m_l.resize(size,size); + m_l.resizeNonZeros(Lstore->nnz); + m_u.resize(size,size); + m_u.resizeNonZeros(Ustore->nnz); + + int* Lcol = m_l._outerIndexPtr(); + int* Lrow = m_l._innerIndexPtr(); + Scalar* Lval = m_l._valuePtr(); + + int* Ucol = m_u._outerIndexPtr(); + int* Urow = m_u._innerIndexPtr(); + Scalar* Uval = m_u._valuePtr(); + + Ucol[0] = 0; + Ucol[0] = 0; + + /* for each supernode */ + for (int k = 0; k <= Lstore->nsuper; ++k) + { + fsupc = L_FST_SUPC(k); + istart = L_SUB_START(fsupc); + nsupr = L_SUB_START(fsupc+1) - istart; + upper = 1; + + /* for each column in the supernode */ + for (int j = fsupc; j < L_FST_SUPC(k+1); ++j) + { + SNptr = &((Scalar*)Lstore->nzval)[L_NZ_START(j)]; + + /* Extract U */ + for (int i = U_NZ_START(j); i < U_NZ_START(j+1); ++i) + { + Uval[lastu] = ((Scalar*)Ustore->nzval)[i]; + /* Matlab doesn't like explicit zero. */ + if (Uval[lastu] != 0.0) + Urow[lastu++] = U_SUB(i); + } + for (int i = 0; i < upper; ++i) + { + /* upper triangle in the supernode */ + Uval[lastu] = SNptr[i]; + /* Matlab doesn't like explicit zero. */ + if (Uval[lastu] != 0.0) + Urow[lastu++] = L_SUB(istart+i); + } + Ucol[j+1] = lastu; + + /* Extract L */ + Lval[lastl] = 1.0; /* unit diagonal */ + Lrow[lastl++] = L_SUB(istart + upper - 1); + for (int i = upper; i < nsupr; ++i) + { + Lval[lastl] = SNptr[i]; + /* Matlab doesn't like explicit zero. */ + if (Lval[lastl] != 0.0) + Lrow[lastl++] = L_SUB(istart+i); + } + Lcol[j+1] = lastl; + + ++upper; + } /* for j ... */ + + } /* for k ... */ + + // squeeze the matrices : + m_l.resizeNonZeros(lastl); + m_u.resizeNonZeros(lastu); + + m_extractedDataAreDirty = false; + } +} + +template +typename SparseLU::Scalar SparseLU::determinant() const +{ + assert((!NumTraits::IsComplex) && "This function is not implemented for complex yet"); + if (m_extractedDataAreDirty) + extractData(); + + // TODO this code could be moved to the default/base backend + // FIXME perhaps we have to take into account the scale factors m_sluRscale and m_sluCscale ??? + Scalar det = Scalar(1); + for (int j=0; j 0) + { + int lastId = m_u._outerIndexPtr()[j+1]-1; + eigen_assert(m_u._innerIndexPtr()[lastId]<=j); + if (m_u._innerIndexPtr()[lastId]==j) + { + det *= m_u._valuePtr()[lastId]; + } + } +// std::cout << m_sluRscale[j] << " " << m_sluCscale[j] << " \n"; + } + return det; +} + +#endif // EIGEN_SUPERLUSUPPORT_LEGACY_H diff --git a/unsupported/test/sparse_lu.cpp b/unsupported/test/sparse_lu.cpp index 188d291cc..bb0dccdb5 100644 --- a/unsupported/test/sparse_lu.cpp +++ b/unsupported/test/sparse_lu.cpp @@ -76,27 +76,62 @@ template void sparse_lu(int rows, int cols) #endif #ifdef EIGEN_SUPERLU_SUPPORT + // legacy, deprecated API { x.setZero(); - SparseLU,SuperLU> slu(m2); + SparseLU,SuperLULegacy> slu(m2); if (slu.succeeded()) { + DenseVector oldb = b; if (slu.solve(b,&x)) { VERIFY(refX.isApprox(x,test_precision()) && "LU: SuperLU"); } + else + std::cerr << "super lu solving failed\n"; + VERIFY(oldb.isApprox(b) && "the rhs should not be modified!"); + // std::cerr << refDet << " == " << slu.determinant() << "\n"; if (slu.solve(b, &x, SvTranspose)) { VERIFY(b.isApprox(m2.transpose() * x, test_precision())); } + else + std::cerr << "super lu solving failed\n"; if (slu.solve(b, &x, SvAdjoint)) { VERIFY(b.isApprox(m2.adjoint() * x, test_precision())); } + else + std::cerr << "super lu solving failed\n"; if (!NumTraits::IsComplex) { VERIFY_IS_APPROX(refDet,slu.determinant()); // FIXME det is not very stable for complex } } + else + std::cerr << "super lu factorize failed\n"; + } + + // New API + { + x.setZero(); + SuperLU > slu(m2); + if (slu.info() == Success) + { + DenseVector oldb = b; + x = slu.solve(b); + VERIFY(oldb.isApprox(b) && "the rhs should not be modified!"); + if (slu.info() == Success) { + VERIFY(refX.isApprox(x,test_precision()) && "SuperLU"); + } + else + std::cerr << "super lu solving failed\n"; + + if (!NumTraits::IsComplex) { + VERIFY_IS_APPROX(refDet,slu.determinant()); // FIXME det is not very stable for complex + } + } + else + std::cerr << "super lu factorize failed\n"; } #endif -- cgit v1.2.3