diff options
author | 2011-07-07 14:19:42 +0200 | |
---|---|---|
committer | 2011-07-07 14:19:42 +0200 | |
commit | 2489c81562b2905384dc3bec6a30debe864f1c33 (patch) | |
tree | e791101af19ae074e0f1ee3182c3c2a3bff41306 | |
parent | c98cd5e564f75be2edc7c5d18491058d025fa796 (diff) |
add new interface to SuperLU
-rw-r--r-- | unsupported/Eigen/SuperLUSupport | 5 | ||||
-rw-r--r-- | unsupported/Eigen/src/SparseExtra/SuperLUSupport.h | 679 | ||||
-rw-r--r-- | unsupported/Eigen/src/SparseExtra/SuperLUSupportLegacy.h | 404 | ||||
-rw-r--r-- | unsupported/test/sparse_lu.cpp | 37 |
4 files changed, 936 insertions, 189 deletions
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<SparseMatrixBase<Derived> > 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<Scalar,Flags,Index> map_superlu(SluMatrix& sluMat) } // end namespace internal -template<typename MatrixType> -class SparseLU<MatrixType,SuperLU> : public SparseLU<MatrixType> + +template<typename _MatrixType, typename Derived> +class SuperLUBase { - protected: - typedef SparseLU<MatrixType> 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<Scalar,Dynamic,1> Vector; typedef Matrix<int, 1, MatrixType::ColsAtCompileTime> IntRowVectorType; - typedef Matrix<int, MatrixType::RowsAtCompileTime, 1> IntColVectorType; - typedef SparseMatrix<Scalar,Lower|UnitDiag> LMatrixType; - typedef SparseMatrix<Scalar,Upper> UMatrixType; - using Base::m_flags; - using Base::m_status; + typedef Matrix<int, MatrixType::RowsAtCompileTime, 1> IntColVectorType; + typedef SparseMatrix<Scalar> 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<Derived*>(this); } + const Derived& derived() const { return *static_cast<const Derived*>(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<typename Rhs> + inline const internal::solve_retval<SuperLUBase, Rhs> solve(const MatrixBase<Rhs>& 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<SuperLUBase, Rhs>(*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<typename Rhs> +// inline const internal::sparse_solve_retval<SuperLU, Rhs> solve(const SparseMatrixBase<Rhs>& 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<SuperLU, Rhs>(*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<typename BDerived, typename XDerived> - bool solve(const MatrixBase<BDerived> &b, MatrixBase<XDerived>* x, const int transposed = SvNoTrans) const; - - void compute(const MatrixType& matrix); - + + template<typename Stream> + 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<Scalar>(); + 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<Scalar> 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<MatrixType,SuperLU> : public SparseLU<MatrixType> mutable std::vector<RealScalar> m_sluRscale, m_sluCscale; mutable std::vector<RealScalar> 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<typename MatrixType> -void SparseLU<MatrixType,SuperLU>::compute(const MatrixType& a) + +template<typename _MatrixType> +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<LUMatrixType, Lower|UnitDiag> LMatrixType; + typedef TriangularView<LUMatrixType, Upper> 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<typename Rhs,typename Dest> + void _solve(const MatrixBase<Rhs> &b, MatrixBase<Dest> &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<Scalar>(); - 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<typename MatrixType> +void SuperLU<MatrixType>::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<typename MatrixType> -template<typename BDerived,typename XDerived> -bool SparseLU<MatrixType,SuperLU>::solve(const MatrixBase<BDerived> &b, - MatrixBase<XDerived> *x, const int transposed) const +template<typename Rhs,typename Dest> +void SuperLU<MatrixType>::_solve(const MatrixBase<Rhs> &b, MatrixBase<Dest>& 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<MatrixType,SuperLU>::solve(const MatrixBase<BDerived> &b, // THIS MATERIAL IS PROVIDED AS IS, WITH ABSOLUTELY NO WARRANTY // EXPRESSED OR IMPLIED. ANY USE IS AT YOUR OWN RISK. // -template<typename MatrixType> -void SparseLU<MatrixType,SuperLU>::extractData() const +template<typename MatrixType, typename Derived> +void SuperLUBase<MatrixType,Derived>::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<MatrixType,SuperLU>::extractData() const } template<typename MatrixType> -typename SparseLU<MatrixType,SuperLU>::Scalar SparseLU<MatrixType,SuperLU>::determinant() const +typename SuperLU<MatrixType>::Scalar SuperLU<MatrixType>::determinant() const { - assert((!NumTraits<Scalar>::IsComplex) && "This function is not implemented for complex yet"); + eigen_assert((!NumTraits<Scalar>::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<m_u.cols(); ++j) @@ -659,9 +765,210 @@ typename SparseLU<MatrixType,SuperLU>::Scalar SparseLU<MatrixType,SuperLU>::dete det *= m_u._valuePtr()[lastId]; } } -// std::cout << m_sluRscale[j] << " " << m_sluCscale[j] << " \n"; } return det; } +#ifdef EIGEN_SUPERLU_HAS_ILU +template<typename _MatrixType> +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<typename Rhs,typename Dest> + void _solve(const MatrixBase<Rhs> &b, MatrixBase<Dest> &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<Scalar>::dummy_precision()*10; + } +}; + +template<typename MatrixType> +void SuperILU<MatrixType>::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<typename MatrixType> +template<typename Rhs,typename Dest> +void SuperILU<MatrixType>::_solve(const MatrixBase<Rhs> &b, MatrixBase<Dest>& 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<typename _MatrixType, typename Derived, typename Rhs> +struct solve_retval<SuperLUBase<_MatrixType,Derived>, Rhs> + : solve_retval_base<SuperLUBase<_MatrixType,Derived>, Rhs> +{ + typedef SuperLUBase<_MatrixType,Derived> Dec; + EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs) + + template<typename Dest> void evalTo(Dest& dst) const + { + dec().derived()._solve(rhs(),dst); + } +}; + +template<typename _MatrixType, typename Derived, typename Rhs> +struct sparse_solve_retval<SuperLUBase<_MatrixType,Derived>, Rhs> + : sparse_solve_retval_base<SuperLUBase<_MatrixType,Derived>, Rhs> +{ + typedef SuperLUBase<_MatrixType,Derived> Dec; + EIGEN_MAKE_SPARSE_SOLVE_HELPERS(Dec,Rhs) + + template<typename Dest> 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 <gael.guennebaud@inria.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_SUPERLUSUPPORT_LEGACY_H +#define EIGEN_SUPERLUSUPPORT_LEGACY_H + +template<typename MatrixType> +class SparseLU<MatrixType,SuperLULegacy> : public SparseLU<MatrixType> +{ + protected: + typedef SparseLU<MatrixType> Base; + typedef typename Base::Scalar Scalar; + typedef typename Base::RealScalar RealScalar; + typedef Matrix<Scalar,Dynamic,1> Vector; + typedef Matrix<int, 1, MatrixType::ColsAtCompileTime> IntRowVectorType; + typedef Matrix<int, MatrixType::RowsAtCompileTime, 1> IntColVectorType; + typedef SparseMatrix<Scalar,Lower|UnitDiag> LMatrixType; + typedef SparseMatrix<Scalar,Upper> 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<typename BDerived, typename XDerived> + bool solve(const MatrixBase<BDerived> &b, MatrixBase<XDerived>* 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<Scalar> 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<int> m_sluEtree; + mutable std::vector<RealScalar> m_sluRscale, m_sluCscale; + mutable std::vector<RealScalar> m_sluFerr, m_sluBerr; + mutable char m_sluEqued; + mutable bool m_extractedDataAreDirty; +}; + +template<typename MatrixType> +void SparseLU<MatrixType,SuperLULegacy>::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<Scalar>(); + 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<typename MatrixType> +template<typename BDerived,typename XDerived> +bool SparseLU<MatrixType,SuperLULegacy>::solve(const MatrixBase<BDerived> &b, + MatrixBase<XDerived> *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<typename MatrixType> +void SparseLU<MatrixType,SuperLULegacy>::extractData() const +{ + if (m_extractedDataAreDirty) + { + int upper; + int fsupc, istart, nsupr; + int lastl = 0, lastu = 0; + SCformat *Lstore = static_cast<SCformat*>(m_sluL.Store); + NCformat *Ustore = static_cast<NCformat*>(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 MatrixType> +typename SparseLU<MatrixType,SuperLULegacy>::Scalar SparseLU<MatrixType,SuperLULegacy>::determinant() const +{ + assert((!NumTraits<Scalar>::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<m_u.cols(); ++j) + { + if (m_u._outerIndexPtr()[j+1]-m_u._outerIndexPtr()[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<typename Scalar> void sparse_lu(int rows, int cols) #endif #ifdef EIGEN_SUPERLU_SUPPORT + // legacy, deprecated API { x.setZero(); - SparseLU<SparseMatrix<Scalar>,SuperLU> slu(m2); + SparseLU<SparseMatrix<Scalar>,SuperLULegacy> slu(m2); if (slu.succeeded()) { + DenseVector oldb = b; if (slu.solve(b,&x)) { VERIFY(refX.isApprox(x,test_precision<Scalar>()) && "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<Scalar>())); } + else + std::cerr << "super lu solving failed\n"; if (slu.solve(b, &x, SvAdjoint)) { VERIFY(b.isApprox(m2.adjoint() * x, test_precision<Scalar>())); } + else + std::cerr << "super lu solving failed\n"; if (!NumTraits<Scalar>::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<SparseMatrix<Scalar> > 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<Scalar>()) && "SuperLU"); + } + else + std::cerr << "super lu solving failed\n"; + + if (!NumTraits<Scalar>::IsComplex) { + VERIFY_IS_APPROX(refDet,slu.determinant()); // FIXME det is not very stable for complex + } + } + else + std::cerr << "super lu factorize failed\n"; } #endif |