diff options
Diffstat (limited to 'Eigen/src/SuperLUSupport/SuperLUSupport.h')
-rw-r--r-- | Eigen/src/SuperLUSupport/SuperLUSupport.h | 105 |
1 files changed, 33 insertions, 72 deletions
diff --git a/Eigen/src/SuperLUSupport/SuperLUSupport.h b/Eigen/src/SuperLUSupport/SuperLUSupport.h index bcb355760..6de5b3dc5 100644 --- a/Eigen/src/SuperLUSupport/SuperLUSupport.h +++ b/Eigen/src/SuperLUSupport/SuperLUSupport.h @@ -1,7 +1,7 @@ // This file is part of Eigen, a lightweight C++ template library // for linear algebra. // -// Copyright (C) 2008-2011 Gael Guennebaud <gael.guennebaud@inria.fr> +// Copyright (C) 2008-2014 Gael Guennebaud <gael.guennebaud@inria.fr> // // This Source Code Form is subject to the terms of the Mozilla // Public License v. 2.0. If a copy of the MPL was not distributed @@ -288,8 +288,12 @@ MappedSparseMatrix<Scalar,Flags,Index> map_superlu(SluMatrix& sluMat) * \brief The base class for the direct and incomplete LU factorization of SuperLU */ template<typename _MatrixType, typename Derived> -class SuperLUBase : internal::noncopyable +class SuperLUBase : public SparseSolverBase<Derived> { + protected: + typedef SparseSolverBase<Derived> Base; + using Base::derived; + using Base::m_isInitialized; public: typedef _MatrixType MatrixType; typedef typename MatrixType::Scalar Scalar; @@ -309,9 +313,6 @@ class SuperLUBase : internal::noncopyable clearFactors(); } - 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(); } @@ -335,33 +336,7 @@ class SuperLUBase : internal::noncopyable derived().analyzePattern(matrix); derived().factorize(matrix); } - - /** \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 - { - 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()); - } - - /** \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<SuperLUBase, 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<SuperLUBase, 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. @@ -453,7 +428,6 @@ class SuperLUBase : internal::noncopyable mutable char m_sluEqued; mutable ComputationInfo m_info; - bool m_isInitialized; int m_factorizationIsOk; int m_analysisIsOk; mutable bool m_extractedDataAreDirty; @@ -491,10 +465,11 @@ class SuperLU : public SuperLUBase<_MatrixType,SuperLU<_MatrixType> > typedef TriangularView<LUMatrixType, Upper> UMatrixType; public: + using Base::_solve_impl; SuperLU() : Base() { init(); } - SuperLU(const MatrixType& matrix) : Base() + explicit SuperLU(const MatrixType& matrix) : Base() { init(); Base::compute(matrix); @@ -528,7 +503,7 @@ class SuperLU : public SuperLUBase<_MatrixType,SuperLU<_MatrixType> > #ifndef EIGEN_PARSED_BY_DOXYGEN /** \internal */ template<typename Rhs,typename Dest> - void _solve(const MatrixBase<Rhs> &b, MatrixBase<Dest> &dest) const; + void _solve_impl(const MatrixBase<Rhs> &b, MatrixBase<Dest> &dest) const; #endif // EIGEN_PARSED_BY_DOXYGEN inline const LMatrixType& matrixL() const @@ -637,7 +612,7 @@ void SuperLU<MatrixType>::factorize(const MatrixType& a) template<typename MatrixType> template<typename Rhs,typename Dest> -void SuperLU<MatrixType>::_solve(const MatrixBase<Rhs> &b, MatrixBase<Dest>& x) const +void SuperLU<MatrixType>::_solve_impl(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()"); @@ -652,8 +627,12 @@ void SuperLU<MatrixType>::_solve(const MatrixBase<Rhs> &b, MatrixBase<Dest>& x) m_sluFerr.resize(rhsCols); m_sluBerr.resize(rhsCols); - m_sluB = SluMatrix::Map(b.const_cast_derived()); - m_sluX = SluMatrix::Map(x.derived()); + + Ref<const Matrix<typename Rhs::Scalar,Dynamic,Dynamic,ColMajor> > b_ref(b); + Ref<const Matrix<typename Dest::Scalar,Dynamic,Dynamic,ColMajor> > x_ref(x); + + m_sluB = SluMatrix::Map(b_ref.const_cast_derived()); + m_sluX = SluMatrix::Map(x_ref.const_cast_derived()); typename Rhs::PlainObject b_cpy; if(m_sluEqued!='N') @@ -676,6 +655,10 @@ void SuperLU<MatrixType>::_solve(const MatrixBase<Rhs> &b, MatrixBase<Dest>& x) &m_sluFerr[0], &m_sluBerr[0], &m_sluStat, &info, Scalar()); StatFree(&m_sluStat); + + if(&x.coeffRef(0) != x_ref.data()) + x = x_ref; + m_info = info==0 ? Success : NumericalIssue; } @@ -828,6 +811,7 @@ class SuperILU : public SuperLUBase<_MatrixType,SuperILU<_MatrixType> > typedef typename Base::Index Index; public: + using Base::_solve_impl; SuperILU() : Base() { init(); } @@ -863,7 +847,7 @@ class SuperILU : public SuperLUBase<_MatrixType,SuperILU<_MatrixType> > #ifndef EIGEN_PARSED_BY_DOXYGEN /** \internal */ template<typename Rhs,typename Dest> - void _solve(const MatrixBase<Rhs> &b, MatrixBase<Dest> &dest) const; + void _solve_impl(const MatrixBase<Rhs> &b, MatrixBase<Dest> &dest) const; #endif // EIGEN_PARSED_BY_DOXYGEN protected: @@ -948,7 +932,7 @@ void SuperILU<MatrixType>::factorize(const MatrixType& a) template<typename MatrixType> template<typename Rhs,typename Dest> -void SuperILU<MatrixType>::_solve(const MatrixBase<Rhs> &b, MatrixBase<Dest>& x) const +void SuperILU<MatrixType>::_solve_impl(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()"); @@ -962,8 +946,12 @@ void SuperILU<MatrixType>::_solve(const MatrixBase<Rhs> &b, MatrixBase<Dest>& x) m_sluFerr.resize(rhsCols); m_sluBerr.resize(rhsCols); - m_sluB = SluMatrix::Map(b.const_cast_derived()); - m_sluX = SluMatrix::Map(x.derived()); + + Ref<const Matrix<typename Rhs::Scalar,Dynamic,Dynamic,ColMajor> > b_ref(b); + Ref<const Matrix<typename Dest::Scalar,Dynamic,Dynamic,ColMajor> > x_ref(x); + + m_sluB = SluMatrix::Map(b_ref.const_cast_derived()); + m_sluX = SluMatrix::Map(x_ref.const_cast_derived()); typename Rhs::PlainObject b_cpy; if(m_sluEqued!='N') @@ -986,41 +974,14 @@ void SuperILU<MatrixType>::_solve(const MatrixBase<Rhs> &b, MatrixBase<Dest>& x) &recip_pivot_growth, &rcond, &m_sluStat, &info, Scalar()); StatFree(&m_sluStat); + + if(&x.coeffRef(0) != x_ref.data()) + x = x_ref; 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 - { - this->defaultEvalTo(dst); - } -}; - -} // end namespace internal - } // end namespace Eigen #endif // EIGEN_SUPERLUSUPPORT_H |