diff options
author | Daniel Genrich <daniel.genrich@gmx.net> | 2012-10-15 20:29:23 +0400 |
---|---|---|
committer | Daniel Genrich <daniel.genrich@gmx.net> | 2012-10-15 20:29:23 +0400 |
commit | 827c70abd8c81089af13c4738e0586bf44e501ea (patch) | |
tree | 2485e8cadb6fcc7363322431eebd22fe846de858 /extern/Eigen3/Eigen/src/IterativeLinearSolvers/IterativeSolverBase.h | |
parent | e3ea7187cedfb92741dc6667a604a873a56670e7 (diff) |
Update to stable Eigen 3.1.1
- Fixes several bugs within the Eigen library:
http://eigen.tuxfamily.org/index.php?title=ChangeLog#Eigen_3.1.1
Diffstat (limited to 'extern/Eigen3/Eigen/src/IterativeLinearSolvers/IterativeSolverBase.h')
-rw-r--r-- | extern/Eigen3/Eigen/src/IterativeLinearSolvers/IterativeSolverBase.h | 254 |
1 files changed, 254 insertions, 0 deletions
diff --git a/extern/Eigen3/Eigen/src/IterativeLinearSolvers/IterativeSolverBase.h b/extern/Eigen3/Eigen/src/IterativeLinearSolvers/IterativeSolverBase.h new file mode 100644 index 00000000000..11706cebabd --- /dev/null +++ b/extern/Eigen3/Eigen/src/IterativeLinearSolvers/IterativeSolverBase.h @@ -0,0 +1,254 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2011 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 +// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#ifndef EIGEN_ITERATIVE_SOLVER_BASE_H +#define EIGEN_ITERATIVE_SOLVER_BASE_H + +namespace Eigen { + +/** \ingroup IterativeLinearSolvers_Module + * \brief Base class for linear iterative solvers + * + * \sa class SimplicialCholesky, DiagonalPreconditioner, IdentityPreconditioner + */ +template< typename Derived> +class IterativeSolverBase : internal::noncopyable +{ +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::RealScalar RealScalar; + +public: + + Derived& derived() { return *static_cast<Derived*>(this); } + const Derived& derived() const { return *static_cast<const Derived*>(this); } + + /** Default constructor. */ + IterativeSolverBase() + : mp_matrix(0) + { + init(); + } + + /** Initialize the solver with matrix \a A for further \c Ax=b solving. + * + * This constructor is a shortcut for the default constructor followed + * by a call to compute(). + * + * \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. + */ + IterativeSolverBase(const MatrixType& A) + { + init(); + compute(A); + } + + ~IterativeSolverBase() {} + + /** Initializes the iterative solver for the sparcity 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. + */ + Derived& analyzePattern(const MatrixType& A) + { + m_preconditioner.analyzePattern(A); + m_isInitialized = true; + m_analysisIsOk = true; + m_info = Success; + 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. + * + * \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) + { + eigen_assert(m_analysisIsOk && "You must first call analyzePattern()"); + mp_matrix = &A; + m_preconditioner.factorize(A); + m_factorizationIsOk = true; + m_info = Success; + 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. + * + * \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) + { + mp_matrix = &A; + m_preconditioner.compute(A); + m_isInitialized = true; + m_analysisIsOk = true; + m_factorizationIsOk = true; + m_info = Success; + return derived(); + } + + /** \internal */ + Index rows() const { return mp_matrix ? mp_matrix->rows() : 0; } + /** \internal */ + Index cols() const { return mp_matrix ? mp_matrix->cols() : 0; } + + /** \returns the tolerance threshold used by the stopping criteria */ + RealScalar tolerance() const { return m_tolerance; } + + /** Sets the tolerance threshold used by the stopping criteria */ + Derived& setTolerance(RealScalar tolerance) + { + m_tolerance = tolerance; + return derived(); + } + + /** \returns a read-write reference to the preconditioner for custom configuration. */ + Preconditioner& preconditioner() { return m_preconditioner; } + + /** \returns a read-only reference to the preconditioner. */ + const Preconditioner& preconditioner() const { return m_preconditioner; } + + /** \returns the max number of iterations */ + int maxIterations() const + { + return (mp_matrix && m_maxIterations<0) ? mp_matrix->cols() : m_maxIterations; + } + + /** Sets the max number of iterations */ + Derived& setMaxIterations(int maxIters) + { + m_maxIterations = maxIters; + return derived(); + } + + /** \returns the number of iterations performed during the last solve */ + int iterations() const + { + eigen_assert(m_isInitialized && "ConjugateGradient is not initialized."); + return m_iterations; + } + + /** \returns the tolerance error reached during the last solve */ + 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. + * + * \sa compute() + */ + template<typename Rhs> + inline const internal::sparse_solve_retval<IterativeSolverBase, Rhs> + solve(const SparseMatrixBase<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::sparse_solve_retval<IterativeSolverBase, Rhs>(*this, b.derived()); + } + + /** \returns Success if the iterations converged, and NoConvergence otherwise. */ + ComputationInfo info() const + { + eigen_assert(m_isInitialized && "IterativeSolverBase is not initialized."); + return m_info; + } + + /** \internal */ + template<typename Rhs, typename DestScalar, int DestOptions, typename DestIndex> + void _solve_sparse(const Rhs& b, SparseMatrix<DestScalar,DestOptions,DestIndex> &dest) const + { + eigen_assert(rows()==b.rows()); + + int rhsCols = b.cols(); + int size = b.rows(); + Eigen::Matrix<DestScalar,Dynamic,1> tb(size); + Eigen::Matrix<DestScalar,Dynamic,1> tx(size); + for(int k=0; k<rhsCols; ++k) + { + tb = b.col(k); + tx = derived().solve(tb); + dest.col(k) = tx.sparseView(0); + } + } + +protected: + void init() + { + m_isInitialized = false; + m_analysisIsOk = false; + m_factorizationIsOk = false; + m_maxIterations = -1; + m_tolerance = NumTraits<Scalar>::epsilon(); + } + const MatrixType* mp_matrix; + Preconditioner m_preconditioner; + + int m_maxIterations; + RealScalar m_tolerance; + + mutable RealScalar m_error; + mutable int 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); + } +}; + +} // end namespace internal + +} // end namespace Eigen + +#endif // EIGEN_ITERATIVE_SOLVER_BASE_H |