diff options
Diffstat (limited to 'Eigen/src/IterativeLinearSolvers/IterativeSolverBase.h')
| -rw-r--r-- | Eigen/src/IterativeLinearSolvers/IterativeSolverBase.h | 90 |
1 files changed, 70 insertions, 20 deletions
diff --git a/Eigen/src/IterativeLinearSolvers/IterativeSolverBase.h b/Eigen/src/IterativeLinearSolvers/IterativeSolverBase.h index 7c2326eb7..28a0c5109 100644 --- a/Eigen/src/IterativeLinearSolvers/IterativeSolverBase.h +++ b/Eigen/src/IterativeLinearSolvers/IterativeSolverBase.h @@ -10,7 +10,7 @@ #ifndef EIGEN_ITERATIVE_SOLVER_BASE_H #define EIGEN_ITERATIVE_SOLVER_BASE_H -namespace Eigen { +namespace Eigen { namespace internal { @@ -145,7 +145,7 @@ 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; @@ -169,10 +169,10 @@ public: } /** 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 @@ -187,7 +187,7 @@ public: } ~IterativeSolverBase() {} - + /** Initializes the iterative solver for the sparsity pattern of the matrix \a A for further solving \c Ax=b problems. * * Currently, this function mostly calls analyzePattern on the preconditioner. In the future @@ -203,7 +203,7 @@ public: 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 calls factorize on the preconditioner. @@ -216,7 +216,7 @@ public: template<typename MatrixDerived> Derived& factorize(const EigenBase<MatrixDerived>& A) { - eigen_assert(m_analysisIsOk && "You must first call analyzePattern()"); + eigen_assert(m_analysisIsOk && "You must first call analyzePattern()"); grab(A.derived()); m_preconditioner.factorize(matrix()); m_factorizationIsOk = true; @@ -247,16 +247,16 @@ public: } /** \internal */ - Index rows() const { return matrix().rows(); } + EIGEN_CONSTEXPR Index rows() const EIGEN_NOEXCEPT { return matrix().rows(); } /** \internal */ - Index cols() const { return matrix().cols(); } + EIGEN_CONSTEXPR Index cols() const EIGEN_NOEXCEPT { return matrix().cols(); } /** \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. * * This value is used as an upper bound to the relative residual error: |Ax-b|/|b|. @@ -270,19 +270,19 @@ public: /** \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. - * It is either the value setted by setMaxIterations or, by default, + * It is either the value set by setMaxIterations or, by default, * twice the number of columns of the matrix. */ Index maxIterations() const { return (m_maxIterations<0) ? 2*matrix().cols() : m_maxIterations; } - + /** Sets the max number of iterations. * Default is twice the number of columns of the matrix. */ @@ -328,13 +328,13 @@ public: eigen_assert(m_isInitialized && "IterativeSolverBase is not initialized."); return m_info; } - + /** \internal */ template<typename Rhs, typename DestDerived> - void _solve_impl(const Rhs& b, SparseMatrixBase<DestDerived> &aDest) const + void _solve_with_guess_impl(const Rhs& b, SparseMatrixBase<DestDerived> &aDest) const { eigen_assert(rows()==b.rows()); - + Index rhsCols = b.cols(); Index size = b.rows(); DestDerived& dest(aDest.derived()); @@ -344,15 +344,65 @@ public: // 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); + ComputationInfo global_info = Success; for(Index k=0; k<rhsCols; ++k) { tb = b.col(k); - tx = derived().solve(tb); + tx = dest.col(k); + derived()._solve_vector_with_guess_impl(tb,tx); tmp.col(k) = tx.sparseView(0); + + // The call to _solve_vector_with_guess_impl updates m_info, so if it failed for a previous column + // we need to restore it to the worst value. + if(m_info==NumericalIssue) + global_info = NumericalIssue; + else if(m_info==NoConvergence) + global_info = NoConvergence; } + m_info = global_info; dest.swap(tmp); } + template<typename Rhs, typename DestDerived> + typename internal::enable_if<Rhs::ColsAtCompileTime!=1 && DestDerived::ColsAtCompileTime!=1>::type + _solve_with_guess_impl(const Rhs& b, MatrixBase<DestDerived> &aDest) const + { + eigen_assert(rows()==b.rows()); + + Index rhsCols = b.cols(); + DestDerived& dest(aDest.derived()); + ComputationInfo global_info = Success; + for(Index k=0; k<rhsCols; ++k) + { + typename DestDerived::ColXpr xk(dest,k); + typename Rhs::ConstColXpr bk(b,k); + derived()._solve_vector_with_guess_impl(bk,xk); + + // The call to _solve_vector_with_guess updates m_info, so if it failed for a previous column + // we need to restore it to the worst value. + if(m_info==NumericalIssue) + global_info = NumericalIssue; + else if(m_info==NoConvergence) + global_info = NoConvergence; + } + m_info = global_info; + } + + template<typename Rhs, typename DestDerived> + typename internal::enable_if<Rhs::ColsAtCompileTime==1 || DestDerived::ColsAtCompileTime==1>::type + _solve_with_guess_impl(const Rhs& b, MatrixBase<DestDerived> &dest) const + { + derived()._solve_vector_with_guess_impl(b,dest.derived()); + } + + /** \internal default initial guess = 0 */ + template<typename Rhs,typename Dest> + void _solve_impl(const Rhs& b, Dest& x) const + { + x.setZero(); + derived()._solve_with_guess_impl(b,x); + } + protected: void init() { @@ -370,19 +420,19 @@ protected: { return m_matrixWrapper.matrix(); } - + template<typename InputType> void grab(const InputType &A) { m_matrixWrapper.grab(A); } - + MatrixWrapper m_matrixWrapper; Preconditioner m_preconditioner; Index m_maxIterations; RealScalar m_tolerance; - + mutable RealScalar m_error; mutable Index m_iterations; mutable ComputationInfo m_info; |