diff options
Diffstat (limited to 'unsupported/Eigen/src/SparseExtra/SuperLUSupport.h')
| -rw-r--r-- | unsupported/Eigen/src/SparseExtra/SuperLUSupport.h | 679 |
1 files changed, 493 insertions, 186 deletions
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 |