From 35f7829af10c61e33dd2e2a7a015058e11a11ea0 Mon Sep 17 00:00:00 2001 From: Stanislaw Halik Date: Sat, 25 Mar 2017 14:17:07 +0100 Subject: update --- .../src/IterativeLinearSolvers/ConjugateGradient.h | 137 ++++++++++----------- 1 file changed, 62 insertions(+), 75 deletions(-) (limited to 'eigen/Eigen/src/IterativeLinearSolvers/ConjugateGradient.h') diff --git a/eigen/Eigen/src/IterativeLinearSolvers/ConjugateGradient.h b/eigen/Eigen/src/IterativeLinearSolvers/ConjugateGradient.h index 7dd4010..395daa8 100644 --- a/eigen/Eigen/src/IterativeLinearSolvers/ConjugateGradient.h +++ b/eigen/Eigen/src/IterativeLinearSolvers/ConjugateGradient.h @@ -1,7 +1,7 @@ // This file is part of Eigen, a lightweight C++ template library // for linear algebra. // -// Copyright (C) 2011 Gael Guennebaud +// Copyright (C) 2011-2014 Gael Guennebaud // // 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 @@ -26,7 +26,7 @@ namespace internal { template EIGEN_DONT_INLINE void conjugate_gradient(const MatrixType& mat, const Rhs& rhs, Dest& x, - const Preconditioner& precond, int& iters, + const Preconditioner& precond, Index& iters, typename Dest::RealScalar& tol_error) { using std::sqrt; @@ -36,9 +36,9 @@ void conjugate_gradient(const MatrixType& mat, const Rhs& rhs, Dest& x, typedef Matrix VectorType; RealScalar tol = tol_error; - int maxIters = iters; + Index maxIters = iters; - int n = mat.cols(); + Index n = mat.cols(); VectorType residual = rhs - mat * x; //initial residual @@ -60,29 +60,29 @@ void conjugate_gradient(const MatrixType& mat, const Rhs& rhs, Dest& x, } VectorType p(n); - p = precond.solve(residual); //initial search direction + p = precond.solve(residual); // initial search direction VectorType z(n), tmp(n); RealScalar absNew = numext::real(residual.dot(p)); // the square of the absolute value of r scaled by invM - int i = 0; + Index i = 0; while(i < maxIters) { - tmp.noalias() = mat * p; // the bottleneck of the algorithm + tmp.noalias() = mat * p; // the bottleneck of the algorithm - Scalar alpha = absNew / p.dot(tmp); // the amount we travel on dir - x += alpha * p; // update solution - residual -= alpha * tmp; // update residue + Scalar alpha = absNew / p.dot(tmp); // the amount we travel on dir + x += alpha * p; // update solution + residual -= alpha * tmp; // update residual residualNorm2 = residual.squaredNorm(); if(residualNorm2 < threshold) break; - z = precond.solve(residual); // approximately solve for "A z = residual" + z = precond.solve(residual); // approximately solve for "A z = residual" RealScalar absOld = absNew; absNew = numext::real(residual.dot(z)); // update the absolute value of r - RealScalar beta = absNew / absOld; // calculate the Gram-Schmidt value used to create the new search direction - p = z + beta * p; // update search direction + RealScalar beta = absNew / absOld; // calculate the Gram-Schmidt value used to create the new search direction + p = z + beta * p; // update search direction i++; } tol_error = sqrt(residualNorm2 / rhsNorm2); @@ -107,47 +107,57 @@ struct traits > } /** \ingroup IterativeLinearSolvers_Module - * \brief A conjugate gradient solver for sparse self-adjoint problems + * \brief A conjugate gradient solver for sparse (or dense) self-adjoint problems * - * This class allows to solve for A.x = b sparse linear problems using a conjugate gradient algorithm. - * The sparse matrix A must be selfadjoint. The vectors x and b can be either dense or sparse. + * This class allows to solve for A.x = b linear problems using an iterative conjugate gradient algorithm. + * The matrix A must be selfadjoint. The matrix A and the vectors x and b can be either dense or sparse. * * \tparam _MatrixType the type of the matrix A, can be a dense or a sparse matrix. * \tparam _UpLo the triangular part that will be used for the computations. It can be Lower, - * Upper, or Lower|Upper in which the full matrix entries will be considered. Default is Lower. + * \c Upper, or \c Lower|Upper in which the full matrix entries will be considered. + * Default is \c Lower, best performance is \c Lower|Upper. * \tparam _Preconditioner the type of the preconditioner. Default is DiagonalPreconditioner * + * \implsparsesolverconcept + * * The maximal number of iterations and tolerance value can be controlled via the setMaxIterations() * and setTolerance() methods. The defaults are the size of the problem for the maximal number of iterations * and NumTraits::epsilon() for the tolerance. * + * The tolerance corresponds to the relative residual error: |Ax-b|/|b| + * + * \b Performance: Even though the default value of \c _UpLo is \c Lower, significantly higher performance is + * achieved when using a complete matrix and \b Lower|Upper as the \a _UpLo template parameter. Moreover, in this + * case multi-threading can be exploited if the user code is compiled with OpenMP enabled. + * See \ref TopicMultiThreading for details. + * * This class can be used as the direct solver classes. Here is a typical usage example: - * \code - * int n = 10000; - * VectorXd x(n), b(n); - * SparseMatrix A(n,n); - * // fill A and b - * ConjugateGradient > cg; - * cg.compute(A); - * x = cg.solve(b); - * std::cout << "#iterations: " << cg.iterations() << std::endl; - * std::cout << "estimated error: " << cg.error() << std::endl; - * // update b, and solve again - * x = cg.solve(b); - * \endcode + \code + int n = 10000; + VectorXd x(n), b(n); + SparseMatrix A(n,n); + // fill A and b + ConjugateGradient, Lower|Upper> cg; + cg.compute(A); + x = cg.solve(b); + std::cout << "#iterations: " << cg.iterations() << std::endl; + std::cout << "estimated error: " << cg.error() << std::endl; + // update b, and solve again + x = cg.solve(b); + \endcode * * By default the iterations start with x=0 as an initial guess of the solution. * One can control the start using the solveWithGuess() method. * * ConjugateGradient can also be used in a matrix-free context, see the following \link MatrixfreeSolverExample example \endlink. * - * \sa class SimplicialCholesky, DiagonalPreconditioner, IdentityPreconditioner + * \sa class LeastSquaresConjugateGradient, class SimplicialCholesky, DiagonalPreconditioner, IdentityPreconditioner */ template< typename _MatrixType, int _UpLo, typename _Preconditioner> class ConjugateGradient : public IterativeSolverBase > { typedef IterativeSolverBase Base; - using Base::mp_matrix; + using Base::matrix; using Base::m_error; using Base::m_iterations; using Base::m_info; @@ -155,7 +165,6 @@ class ConjugateGradient : public IterativeSolverBase& A) : Base(A.derived()) {} ~ConjugateGradient() {} - - /** \returns the solution x of \f$ A x = b \f$ using the current decomposition of A - * \a x0 as an initial solution. - * - * \sa compute() - */ - template - inline const internal::solve_retval_with_guess - solveWithGuess(const MatrixBase& b, const Guess& x0) const - { - eigen_assert(m_isInitialized && "ConjugateGradient is not initialized."); - eigen_assert(Base::rows()==b.rows() - && "ConjugateGradient::solve(): invalid number of rows of the right hand side matrix b"); - return internal::solve_retval_with_guess - (*this, b.derived(), x0); - } /** \internal */ template - void _solveWithGuess(const Rhs& b, Dest& x) const + void _solve_with_guess_impl(const Rhs& b, Dest& x) const { + typedef typename Base::MatrixWrapper MatrixWrapper; + typedef typename Base::ActualMatrixType ActualMatrixType; + enum { + TransposeInput = (!MatrixWrapper::MatrixFree) + && (UpLo==(Lower|Upper)) + && (!MatrixType::IsRowMajor) + && (!NumTraits::IsComplex) + }; + typedef typename internal::conditional, ActualMatrixType const&>::type RowMajorWrapper; + EIGEN_STATIC_ASSERT(EIGEN_IMPLIES(MatrixWrapper::MatrixFree,UpLo==(Lower|Upper)),MATRIX_FREE_CONJUGATE_GRADIENT_IS_COMPATIBLE_WITH_UPPER_UNION_LOWER_MODE_ONLY); typedef typename internal::conditional - >::type MatrixWrapperType; + RowMajorWrapper, + typename MatrixWrapper::template ConstSelfAdjointViewReturnType::Type + >::type SelfAdjointWrapper; m_iterations = Base::maxIterations(); m_error = Base::m_tolerance; - for(int j=0; j - void _solve(const Rhs& b, Dest& x) const + void _solve_impl(const MatrixBase& b, Dest& x) const { x.setZero(); - _solveWithGuess(b,x); + _solve_with_guess_impl(b.derived(),x); } protected: }; - -namespace internal { - -template -struct solve_retval, Rhs> - : solve_retval_base, Rhs> -{ - typedef ConjugateGradient<_MatrixType,_UpLo,_Preconditioner> Dec; - EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs) - - template void evalTo(Dest& dst) const - { - dec()._solve(rhs(),dst); - } -}; - -} // end namespace internal - } // end namespace Eigen #endif // EIGEN_CONJUGATE_GRADIENT_H -- cgit v1.2.3