diff options
Diffstat (limited to 'Eigen/src/IterativeLinearSolvers/IterativeSolverBase.h')
| -rw-r--r-- | Eigen/src/IterativeLinearSolvers/IterativeSolverBase.h | 312 |
1 files changed, 226 insertions, 86 deletions
diff --git a/Eigen/src/IterativeLinearSolvers/IterativeSolverBase.h b/Eigen/src/IterativeLinearSolvers/IterativeSolverBase.h index 2036922d6..7c2326eb7 100644 --- a/Eigen/src/IterativeLinearSolvers/IterativeSolverBase.h +++ b/Eigen/src/IterativeLinearSolvers/IterativeSolverBase.h @@ -1,7 +1,7 @@ // This file is part of Eigen, a lightweight C++ template library // for linear algebra. // -// Copyright (C) 2011 Gael Guennebaud <gael.guennebaud@inria.fr> +// Copyright (C) 2011-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 @@ -12,29 +12,158 @@ namespace Eigen { +namespace internal { + +template<typename MatrixType> +struct is_ref_compatible_impl +{ +private: + template <typename T0> + struct any_conversion + { + template <typename T> any_conversion(const volatile T&); + template <typename T> any_conversion(T&); + }; + struct yes {int a[1];}; + struct no {int a[2];}; + + template<typename T> + static yes test(const Ref<const T>&, int); + template<typename T> + static no test(any_conversion<T>, ...); + +public: + static MatrixType ms_from; + enum { value = sizeof(test<MatrixType>(ms_from, 0))==sizeof(yes) }; +}; + +template<typename MatrixType> +struct is_ref_compatible +{ + enum { value = is_ref_compatible_impl<typename remove_all<MatrixType>::type>::value }; +}; + +template<typename MatrixType, bool MatrixFree = !internal::is_ref_compatible<MatrixType>::value> +class generic_matrix_wrapper; + +// We have an explicit matrix at hand, compatible with Ref<> +template<typename MatrixType> +class generic_matrix_wrapper<MatrixType,false> +{ +public: + typedef Ref<const MatrixType> ActualMatrixType; + template<int UpLo> struct ConstSelfAdjointViewReturnType { + typedef typename ActualMatrixType::template ConstSelfAdjointViewReturnType<UpLo>::Type Type; + }; + + enum { + MatrixFree = false + }; + + generic_matrix_wrapper() + : m_dummy(0,0), m_matrix(m_dummy) + {} + + template<typename InputType> + generic_matrix_wrapper(const InputType &mat) + : m_matrix(mat) + {} + + const ActualMatrixType& matrix() const + { + return m_matrix; + } + + template<typename MatrixDerived> + void grab(const EigenBase<MatrixDerived> &mat) + { + m_matrix.~Ref<const MatrixType>(); + ::new (&m_matrix) Ref<const MatrixType>(mat.derived()); + } + + void grab(const Ref<const MatrixType> &mat) + { + if(&(mat.derived()) != &m_matrix) + { + m_matrix.~Ref<const MatrixType>(); + ::new (&m_matrix) Ref<const MatrixType>(mat); + } + } + +protected: + MatrixType m_dummy; // used to default initialize the Ref<> object + ActualMatrixType m_matrix; +}; + +// MatrixType is not compatible with Ref<> -> matrix-free wrapper +template<typename MatrixType> +class generic_matrix_wrapper<MatrixType,true> +{ +public: + typedef MatrixType ActualMatrixType; + template<int UpLo> struct ConstSelfAdjointViewReturnType + { + typedef ActualMatrixType Type; + }; + + enum { + MatrixFree = true + }; + + generic_matrix_wrapper() + : mp_matrix(0) + {} + + generic_matrix_wrapper(const MatrixType &mat) + : mp_matrix(&mat) + {} + + const ActualMatrixType& matrix() const + { + return *mp_matrix; + } + + void grab(const MatrixType &mat) + { + mp_matrix = &mat; + } + +protected: + const ActualMatrixType *mp_matrix; +}; + +} + /** \ingroup IterativeLinearSolvers_Module * \brief Base class for linear iterative solvers * * \sa class SimplicialCholesky, DiagonalPreconditioner, IdentityPreconditioner */ template< typename Derived> -class IterativeSolverBase : internal::noncopyable +class IterativeSolverBase : public SparseSolverBase<Derived> { +protected: + typedef SparseSolverBase<Derived> Base; + using Base::m_isInitialized; + public: typedef typename internal::traits<Derived>::MatrixType MatrixType; typedef typename internal::traits<Derived>::Preconditioner Preconditioner; typedef typename MatrixType::Scalar Scalar; - typedef typename MatrixType::Index Index; + typedef typename MatrixType::StorageIndex StorageIndex; typedef typename MatrixType::RealScalar RealScalar; + enum { + ColsAtCompileTime = MatrixType::ColsAtCompileTime, + MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime + }; + public: - Derived& derived() { return *static_cast<Derived*>(this); } - const Derived& derived() const { return *static_cast<const Derived*>(this); } + using Base::derived; /** Default constructor. */ IterativeSolverBase() - : mp_matrix(0) { init(); } @@ -49,77 +178,90 @@ public: * this class becomes invalid. Call compute() to update it with the new * matrix A, or modify a copy of A. */ - IterativeSolverBase(const MatrixType& A) + template<typename MatrixDerived> + explicit IterativeSolverBase(const EigenBase<MatrixDerived>& A) + : m_matrixWrapper(A.derived()) { init(); - compute(A); + compute(matrix()); } ~IterativeSolverBase() {} - /** Initializes the iterative solver for the sparcity pattern of the matrix \a A for further solving \c Ax=b problems. + /** Initializes the iterative solver for the sparsity pattern of the matrix \a A for further solving \c Ax=b problems. * - * Currently, this function mostly call analyzePattern on the preconditioner. In the future - * we might, for instance, implement column reodering for faster matrix vector products. + * Currently, this function mostly calls analyzePattern on the preconditioner. In the future + * we might, for instance, implement column reordering for faster matrix vector products. */ - Derived& analyzePattern(const MatrixType& A) + template<typename MatrixDerived> + Derived& analyzePattern(const EigenBase<MatrixDerived>& A) { - m_preconditioner.analyzePattern(A); + grab(A.derived()); + m_preconditioner.analyzePattern(matrix()); m_isInitialized = true; m_analysisIsOk = true; - m_info = Success; + m_info = m_preconditioner.info(); return derived(); } /** Initializes the iterative solver with the numerical values of the matrix \a A for further solving \c Ax=b problems. * - * Currently, this function mostly call factorize on the preconditioner. + * Currently, this function mostly calls factorize on the preconditioner. * * \warning this class stores a reference to the matrix A as well as some * precomputed values that depend on it. Therefore, if \a A is changed * this class becomes invalid. Call compute() to update it with the new * matrix A, or modify a copy of A. */ - Derived& factorize(const MatrixType& A) + template<typename MatrixDerived> + Derived& factorize(const EigenBase<MatrixDerived>& A) { eigen_assert(m_analysisIsOk && "You must first call analyzePattern()"); - mp_matrix = &A; - m_preconditioner.factorize(A); + grab(A.derived()); + m_preconditioner.factorize(matrix()); m_factorizationIsOk = true; - m_info = Success; + m_info = m_preconditioner.info(); return derived(); } /** Initializes the iterative solver with the matrix \a A for further solving \c Ax=b problems. * - * Currently, this function mostly initialized/compute the preconditioner. In the future - * we might, for instance, implement column reodering for faster matrix vector products. + * Currently, this function mostly initializes/computes the preconditioner. In the future + * we might, for instance, implement column reordering for faster matrix vector products. * * \warning this class stores a reference to the matrix A as well as some * precomputed values that depend on it. Therefore, if \a A is changed * this class becomes invalid. Call compute() to update it with the new * matrix A, or modify a copy of A. */ - Derived& compute(const MatrixType& A) + template<typename MatrixDerived> + Derived& compute(const EigenBase<MatrixDerived>& A) { - mp_matrix = &A; - m_preconditioner.compute(A); + grab(A.derived()); + m_preconditioner.compute(matrix()); m_isInitialized = true; m_analysisIsOk = true; m_factorizationIsOk = true; - m_info = Success; + m_info = m_preconditioner.info(); return derived(); } /** \internal */ - Index rows() const { return mp_matrix ? mp_matrix->rows() : 0; } + Index rows() const { return matrix().rows(); } + /** \internal */ - Index cols() const { return mp_matrix ? mp_matrix->cols() : 0; } + Index cols() const { return matrix().cols(); } - /** \returns the tolerance threshold used by the stopping criteria */ + /** \returns the tolerance threshold used by the stopping criteria. + * \sa setTolerance() + */ RealScalar tolerance() const { return m_tolerance; } - /** Sets the tolerance threshold used by the stopping criteria */ + /** Sets the tolerance threshold used by the stopping criteria. + * + * This value is used as an upper bound to the relative residual error: |Ax-b|/|b|. + * The default value is the machine precision given by NumTraits<Scalar>::epsilon() + */ Derived& setTolerance(const RealScalar& tolerance) { m_tolerance = tolerance; @@ -132,58 +274,52 @@ public: /** \returns a read-only reference to the preconditioner. */ const Preconditioner& preconditioner() const { return m_preconditioner; } - /** \returns the max number of iterations */ - int maxIterations() const + /** \returns the max number of iterations. + * It is either the value setted by setMaxIterations or, by default, + * twice the number of columns of the matrix. + */ + Index maxIterations() const { - return (mp_matrix && m_maxIterations<0) ? mp_matrix->cols() : m_maxIterations; + return (m_maxIterations<0) ? 2*matrix().cols() : m_maxIterations; } - /** Sets the max number of iterations */ - Derived& setMaxIterations(int maxIters) + /** Sets the max number of iterations. + * Default is twice the number of columns of the matrix. + */ + Derived& setMaxIterations(Index maxIters) { m_maxIterations = maxIters; return derived(); } /** \returns the number of iterations performed during the last solve */ - int iterations() const + Index iterations() const { eigen_assert(m_isInitialized && "ConjugateGradient is not initialized."); return m_iterations; } - /** \returns the tolerance error reached during the last solve */ + /** \returns the tolerance error reached during the last solve. + * It is a close approximation of the true relative residual error |Ax-b|/|b|. + */ RealScalar error() const { eigen_assert(m_isInitialized && "ConjugateGradient is not initialized."); return m_error; } - /** \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<Derived, Rhs> - solve(const MatrixBase<Rhs>& b) const - { - eigen_assert(m_isInitialized && "IterativeSolverBase is not initialized."); - eigen_assert(rows()==b.rows() - && "IterativeSolverBase::solve(): invalid number of rows of the right hand side matrix b"); - return internal::solve_retval<Derived, Rhs>(derived(), b.derived()); - } - - /** \returns the solution x of \f$ A x = b \f$ using the current decomposition of A. + /** \returns the solution x of \f$ A x = b \f$ using the current decomposition of A + * and \a x0 as an initial solution. * - * \sa compute() + * \sa solve(), compute() */ - template<typename Rhs> - inline const internal::sparse_solve_retval<IterativeSolverBase, Rhs> - solve(const SparseMatrixBase<Rhs>& b) const + template<typename Rhs,typename Guess> + inline const SolveWithGuess<Derived, Rhs, Guess> + solveWithGuess(const MatrixBase<Rhs>& b, const Guess& x0) const { - eigen_assert(m_isInitialized && "IterativeSolverBase is not initialized."); - eigen_assert(rows()==b.rows() - && "IterativeSolverBase::solve(): invalid number of rows of the right hand side matrix b"); - return internal::sparse_solve_retval<IterativeSolverBase, Rhs>(*this, b.derived()); + eigen_assert(m_isInitialized && "Solver is not initialized."); + eigen_assert(derived().rows()==b.rows() && "solve(): invalid number of rows of the right hand side matrix b"); + return SolveWithGuess<Derived, Rhs, Guess>(derived(), b.derived(), x0); } /** \returns Success if the iterations converged, and NoConvergence otherwise. */ @@ -194,21 +330,27 @@ public: } /** \internal */ - template<typename Rhs, typename DestScalar, int DestOptions, typename DestIndex> - void _solve_sparse(const Rhs& b, SparseMatrix<DestScalar,DestOptions,DestIndex> &dest) const + template<typename Rhs, typename DestDerived> + void _solve_impl(const Rhs& b, SparseMatrixBase<DestDerived> &aDest) const { eigen_assert(rows()==b.rows()); - int rhsCols = b.cols(); - int size = b.rows(); + Index rhsCols = b.cols(); + Index size = b.rows(); + DestDerived& dest(aDest.derived()); + typedef typename DestDerived::Scalar DestScalar; Eigen::Matrix<DestScalar,Dynamic,1> tb(size); - Eigen::Matrix<DestScalar,Dynamic,1> tx(size); - for(int k=0; k<rhsCols; ++k) + Eigen::Matrix<DestScalar,Dynamic,1> tx(cols()); + // We do not directly fill dest because sparse expressions have to be free of aliasing issue. + // For non square least-square problems, b and dest might not have the same size whereas they might alias each-other. + typename DestDerived::PlainObject tmp(cols(),rhsCols); + for(Index k=0; k<rhsCols; ++k) { tb = b.col(k); tx = derived().solve(tb); - dest.col(k) = tx.sparseView(0); + tmp.col(k) = tx.sparseView(0); } + dest.swap(tmp); } protected: @@ -220,35 +362,33 @@ protected: m_maxIterations = -1; m_tolerance = NumTraits<Scalar>::epsilon(); } - const MatrixType* mp_matrix; + + typedef internal::generic_matrix_wrapper<MatrixType> MatrixWrapper; + typedef typename MatrixWrapper::ActualMatrixType ActualMatrixType; + + const ActualMatrixType& matrix() const + { + return m_matrixWrapper.matrix(); + } + + template<typename InputType> + void grab(const InputType &A) + { + m_matrixWrapper.grab(A); + } + + MatrixWrapper m_matrixWrapper; Preconditioner m_preconditioner; - int m_maxIterations; + Index m_maxIterations; RealScalar m_tolerance; mutable RealScalar m_error; - mutable int m_iterations; + mutable Index m_iterations; mutable ComputationInfo m_info; - mutable bool m_isInitialized, m_analysisIsOk, m_factorizationIsOk; -}; - -namespace internal { - -template<typename Derived, typename Rhs> -struct sparse_solve_retval<IterativeSolverBase<Derived>, Rhs> - : sparse_solve_retval_base<IterativeSolverBase<Derived>, Rhs> -{ - typedef IterativeSolverBase<Derived> Dec; - EIGEN_MAKE_SPARSE_SOLVE_HELPERS(Dec,Rhs) - - template<typename Dest> void evalTo(Dest& dst) const - { - dec().derived()._solve_sparse(rhs(),dst); - } + mutable bool m_analysisIsOk, m_factorizationIsOk; }; -} // end namespace internal - } // end namespace Eigen #endif // EIGEN_ITERATIVE_SOLVER_BASE_H |