From b37036afce20e902cd5191a2a985f39b1f7e22e3 Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Mon, 7 Dec 2015 12:23:22 +0100 Subject: [PATCH] Implement wrapper for matrix-free iterative solvers --- Eigen/src/Core/util/Meta.h | 1 - Eigen/src/Core/util/StaticAssert.h | 3 +- Eigen/src/IterativeLinearSolvers/BiCGSTAB.h | 4 +- .../ConjugateGradient.h | 21 ++- .../IterativeSolverBase.h | 162 +++++++++++++++--- .../LeastSquareConjugateGradient.h | 4 +- unsupported/Eigen/IterativeSolvers | 1 + .../Eigen/src/IterativeSolvers/DGMRES.h | 16 +- .../Eigen/src/IterativeSolvers/GMRES.h | 4 +- .../Eigen/src/IterativeSolvers/MINRES.h | 24 ++- 10 files changed, 189 insertions(+), 51 deletions(-) diff --git a/Eigen/src/Core/util/Meta.h b/Eigen/src/Core/util/Meta.h index 15b80abd9..3dee2bd7c 100644 --- a/Eigen/src/Core/util/Meta.h +++ b/Eigen/src/Core/util/Meta.h @@ -237,7 +237,6 @@ protected: EIGEN_DEVICE_FUNC ~noncopyable() {} }; - /** \internal * Convenient struct to get the result type of a unary or binary functor. * diff --git a/Eigen/src/Core/util/StaticAssert.h b/Eigen/src/Core/util/StaticAssert.h index f35ddb372..108181419 100644 --- a/Eigen/src/Core/util/StaticAssert.h +++ b/Eigen/src/Core/util/StaticAssert.h @@ -95,7 +95,8 @@ IMPLICIT_CONVERSION_TO_SCALAR_IS_FOR_INNER_PRODUCT_ONLY, STORAGE_LAYOUT_DOES_NOT_MATCH, EIGEN_INTERNAL_ERROR_PLEASE_FILE_A_BUG_REPORT__INVALID_COST_VALUE, - THIS_COEFFICIENT_ACCESSOR_TAKING_ONE_ACCESS_IS_ONLY_FOR_EXPRESSIONS_ALLOWING_LINEAR_ACCESS + THIS_COEFFICIENT_ACCESSOR_TAKING_ONE_ACCESS_IS_ONLY_FOR_EXPRESSIONS_ALLOWING_LINEAR_ACCESS, + MATRIX_FREE_CONJUGATE_GRADIENT_IS_COMPATIBLE_WITH_UPPER_UNION_LOWER_MODE_ONLY }; }; diff --git a/Eigen/src/IterativeLinearSolvers/BiCGSTAB.h b/Eigen/src/IterativeLinearSolvers/BiCGSTAB.h index 4be00da47..191202138 100644 --- a/Eigen/src/IterativeLinearSolvers/BiCGSTAB.h +++ b/Eigen/src/IterativeLinearSolvers/BiCGSTAB.h @@ -156,7 +156,7 @@ template< typename _MatrixType, typename _Preconditioner> class BiCGSTAB : public IterativeSolverBase > { typedef IterativeSolverBase Base; - using Base::mp_matrix; + using Base::matrix; using Base::m_error; using Base::m_iterations; using Base::m_info; @@ -198,7 +198,7 @@ public: m_error = Base::m_tolerance; typename Dest::ColXpr xj(x,j); - if(!internal::bicgstab(mp_matrix, b.col(j), xj, Base::m_preconditioner, m_iterations, m_error)) + if(!internal::bicgstab(matrix(), b.col(j), xj, Base::m_preconditioner, m_iterations, m_error)) failed = true; } m_info = failed ? NumericalIssue diff --git a/Eigen/src/IterativeLinearSolvers/ConjugateGradient.h b/Eigen/src/IterativeLinearSolvers/ConjugateGradient.h index dbedf28fd..395daa8e4 100644 --- a/Eigen/src/IterativeLinearSolvers/ConjugateGradient.h +++ b/Eigen/src/IterativeLinearSolvers/ConjugateGradient.h @@ -149,13 +149,15 @@ struct traits > * 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 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; @@ -194,12 +196,19 @@ public: template void _solve_with_guess_impl(const Rhs& b, Dest& x) const { - typedef Ref MatRef; - typedef typename internal::conditional::IsComplex), - Transpose, MatRef const&>::type RowMajorWrapper; + 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 + typename MatrixWrapper::template ConstSelfAdjointViewReturnType::Type >::type SelfAdjointWrapper; m_iterations = Base::maxIterations(); m_error = Base::m_tolerance; @@ -210,7 +219,7 @@ public: m_error = Base::m_tolerance; typename Dest::ColXpr xj(x,j); - RowMajorWrapper row_mat(mp_matrix); + RowMajorWrapper row_mat(matrix()); internal::conjugate_gradient(SelfAdjointWrapper(row_mat), b.col(j), xj, Base::m_preconditioner, m_iterations, m_error); } diff --git a/Eigen/src/IterativeLinearSolvers/IterativeSolverBase.h b/Eigen/src/IterativeLinearSolvers/IterativeSolverBase.h index e51ff7280..3d62fef6e 100644 --- a/Eigen/src/IterativeLinearSolvers/IterativeSolverBase.h +++ b/Eigen/src/IterativeLinearSolvers/IterativeSolverBase.h @@ -12,6 +12,128 @@ namespace Eigen { +namespace internal { + +template +struct is_ref_compatible_impl +{ +private: + template + struct any_conversion + { + template any_conversion(const volatile T&); + template any_conversion(T&); + }; + struct yes {int a[1];}; + struct no {int a[2];}; + + template + static yes test(const Ref&, int); + template + static no test(any_conversion, ...); + +public: + static MatrixType ms_from; + enum { value = sizeof(test(ms_from, 0))==sizeof(yes) }; +}; + +template +struct is_ref_compatible +{ + enum { value = is_ref_compatible_impl::type>::value }; +}; + +template::value> +class generic_matrix_wrapper; + +// We have an explicit matrix at hand, compatible with Ref<> +template +class generic_matrix_wrapper +{ +public: + typedef Ref ActualMatrixType; + template struct ConstSelfAdjointViewReturnType { + typedef typename ActualMatrixType::template ConstSelfAdjointViewReturnType::Type Type; + }; + + enum { + MatrixFree = false + }; + + generic_matrix_wrapper() + : m_dummy(0,0), m_matrix(m_dummy) + {} + + template + generic_matrix_wrapper(const InputType &mat) + : m_matrix(mat) + {} + + const ActualMatrixType& matrix() const + { + return m_matrix; + } + + template + void grab(const EigenBase &mat) + { + m_matrix.~Ref(); + ::new (&m_matrix) Ref(mat.derived()); + } + + void grab(const Ref &mat) + { + if(&(mat.derived()) != &m_matrix) + { + m_matrix.~Ref(); + ::new (&m_matrix) Ref(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 +class generic_matrix_wrapper +{ +public: + typedef MatrixType ActualMatrixType; + template 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 * @@ -42,7 +164,6 @@ public: /** Default constructor. */ IterativeSolverBase() - : m_dummy(0,0), mp_matrix(m_dummy) { init(); } @@ -59,10 +180,10 @@ public: */ template explicit IterativeSolverBase(const EigenBase& A) - : mp_matrix(A.derived()) + : m_matrixWrapper(A.derived()) { init(); - compute(mp_matrix); + compute(matrix()); } ~IterativeSolverBase() {} @@ -76,7 +197,7 @@ public: Derived& analyzePattern(const EigenBase& A) { grab(A.derived()); - m_preconditioner.analyzePattern(mp_matrix); + m_preconditioner.analyzePattern(matrix()); m_isInitialized = true; m_analysisIsOk = true; m_info = m_preconditioner.info(); @@ -97,7 +218,7 @@ public: { eigen_assert(m_analysisIsOk && "You must first call analyzePattern()"); grab(A.derived()); - m_preconditioner.factorize(mp_matrix); + m_preconditioner.factorize(matrix()); m_factorizationIsOk = true; m_info = m_preconditioner.info(); return derived(); @@ -117,7 +238,7 @@ public: Derived& compute(const EigenBase& A) { grab(A.derived()); - m_preconditioner.compute(mp_matrix); + m_preconditioner.compute(matrix()); m_isInitialized = true; m_analysisIsOk = true; m_factorizationIsOk = true; @@ -126,10 +247,10 @@ public: } /** \internal */ - Index rows() const { return mp_matrix.rows(); } + Index rows() const { return matrix().rows(); } /** \internal */ - Index cols() const { return mp_matrix.cols(); } + Index cols() const { return matrix().cols(); } /** \returns the tolerance threshold used by the stopping criteria. * \sa setTolerance() @@ -159,7 +280,7 @@ public: */ Index maxIterations() const { - return (m_maxIterations<0) ? 2*mp_matrix.cols() : m_maxIterations; + return (m_maxIterations<0) ? 2*matrix().cols() : m_maxIterations; } /** Sets the max number of iterations. @@ -239,25 +360,22 @@ protected: m_maxIterations = -1; m_tolerance = NumTraits::epsilon(); } - - template - void grab(const EigenBase &A) + + typedef internal::generic_matrix_wrapper MatrixWrapper; + typedef typename MatrixWrapper::ActualMatrixType ActualMatrixType; + + const ActualMatrixType& matrix() const { - mp_matrix.~Ref(); - ::new (&mp_matrix) Ref(A.derived()); + return m_matrixWrapper.matrix(); } - void grab(const Ref &A) + template + void grab(const InputType &A) { - if(&(A.derived()) != &mp_matrix) - { - mp_matrix.~Ref(); - ::new (&mp_matrix) Ref(A); - } + m_matrixWrapper.grab(A); } - MatrixType m_dummy; - Ref mp_matrix; + MatrixWrapper m_matrixWrapper; Preconditioner m_preconditioner; Index m_maxIterations; diff --git a/Eigen/src/IterativeLinearSolvers/LeastSquareConjugateGradient.h b/Eigen/src/IterativeLinearSolvers/LeastSquareConjugateGradient.h index 1593c57b5..0aea0e099 100644 --- a/Eigen/src/IterativeLinearSolvers/LeastSquareConjugateGradient.h +++ b/Eigen/src/IterativeLinearSolvers/LeastSquareConjugateGradient.h @@ -149,7 +149,7 @@ template< typename _MatrixType, typename _Preconditioner> class LeastSquaresConjugateGradient : public IterativeSolverBase > { typedef IterativeSolverBase Base; - using Base::mp_matrix; + using Base::matrix; using Base::m_error; using Base::m_iterations; using Base::m_info; @@ -193,7 +193,7 @@ public: m_error = Base::m_tolerance; typename Dest::ColXpr xj(x,j); - internal::least_square_conjugate_gradient(mp_matrix, b.col(j), xj, Base::m_preconditioner, m_iterations, m_error); + internal::least_square_conjugate_gradient(matrix(), b.col(j), xj, Base::m_preconditioner, m_iterations, m_error); } m_isInitialized = true; diff --git a/unsupported/Eigen/IterativeSolvers b/unsupported/Eigen/IterativeSolvers index f0c017f00..31e880bdc 100644 --- a/unsupported/Eigen/IterativeSolvers +++ b/unsupported/Eigen/IterativeSolvers @@ -33,6 +33,7 @@ #include "../../Eigen/Jacobi" #include "../../Eigen/Householder" #include "src/IterativeSolvers/GMRES.h" +#include "src/IterativeSolvers/DGMRES.h" //#include "src/IterativeSolvers/SSORPreconditioner.h" #include "src/IterativeSolvers/MINRES.h" diff --git a/unsupported/Eigen/src/IterativeSolvers/DGMRES.h b/unsupported/Eigen/src/IterativeSolvers/DGMRES.h index ab82e782d..8a28fc16f 100644 --- a/unsupported/Eigen/src/IterativeSolvers/DGMRES.h +++ b/unsupported/Eigen/src/IterativeSolvers/DGMRES.h @@ -40,7 +40,6 @@ void sortWithPermutation (VectorType& vec, IndexType& perm, typename IndexType:: { eigen_assert(vec.size() == perm.size()); typedef typename IndexType::Scalar Index; - typedef typename VectorType::Scalar Scalar; bool flag; for (Index k = 0; k < ncut; k++) { @@ -101,7 +100,7 @@ template< typename _MatrixType, typename _Preconditioner> class DGMRES : public IterativeSolverBase > { typedef IterativeSolverBase Base; - using Base::mp_matrix; + using Base::matrix; using Base::m_error; using Base::m_iterations; using Base::m_info; @@ -112,6 +111,7 @@ class DGMRES : public IterativeSolverBase > typedef _MatrixType MatrixType; typedef typename MatrixType::Scalar Scalar; typedef typename MatrixType::Index Index; + typedef typename MatrixType::StorageIndex StorageIndex; typedef typename MatrixType::RealScalar RealScalar; typedef _Preconditioner Preconditioner; typedef Matrix DenseMatrix; @@ -150,7 +150,7 @@ class DGMRES : public IterativeSolverBase > m_error = Base::m_tolerance; typename Dest::ColXpr xj(x,j); - dgmres(mp_matrix, b.col(j), xj, Base::m_preconditioner); + dgmres(matrix(), b.col(j), xj, Base::m_preconditioner); } m_info = failed ? NumericalIssue : m_error <= Base::m_tolerance ? Success @@ -202,7 +202,7 @@ class DGMRES : public IterativeSolverBase > template int dgmresCycle(const MatrixType& mat, const Preconditioner& precond, Dest& x, DenseVector& r0, RealScalar& beta, const RealScalar& normRhs, int& nbIts) const; // Compute data to use for deflation - int dgmresComputeDeflationData(const MatrixType& mat, const Preconditioner& precond, const Index& it, Index& neig) const; + int dgmresComputeDeflationData(const MatrixType& mat, const Preconditioner& precond, const Index& it, StorageIndex& neig) const; // Apply deflation to a vector template int dgmresApplyDeflation(const RhsType& In, DestType& Out) const; @@ -218,7 +218,7 @@ class DGMRES : public IterativeSolverBase > mutable DenseMatrix m_MU; // matrix operator applied to m_U (for next cycles) mutable DenseMatrix m_T; /* T=U^T*M^{-1}*A*U */ mutable PartialPivLU m_luT; // LU factorization of m_T - mutable int m_neig; //Number of eigenvalues to extract at each restart + mutable StorageIndex m_neig; //Number of eigenvalues to extract at each restart mutable int m_r; // Current number of deflated eigenvalues, size of m_U mutable int m_maxNeig; // Maximum number of eigenvalues to deflate mutable RealScalar m_lambdaN; //Modulus of the largest eigenvalue of A @@ -338,7 +338,7 @@ int DGMRES<_MatrixType, _Preconditioner>::dgmresCycle(const MatrixType& mat, con beta = std::abs(g(it+1)); m_error = beta/normRhs; - std::cerr << nbIts << " Relative Residual Norm " << m_error << std::endl; + // std::cerr << nbIts << " Relative Residual Norm " << m_error << std::endl; it++; nbIts++; if (m_error < m_tolerance) @@ -416,7 +416,7 @@ inline typename DGMRES<_MatrixType, _Preconditioner>::ComplexVector DGMRES<_Matr } template< typename _MatrixType, typename _Preconditioner> -int DGMRES<_MatrixType, _Preconditioner>::dgmresComputeDeflationData(const MatrixType& mat, const Preconditioner& precond, const Index& it, Index& neig) const +int DGMRES<_MatrixType, _Preconditioner>::dgmresComputeDeflationData(const MatrixType& mat, const Preconditioner& precond, const Index& it, StorageIndex& neig) const { // First, find the Schur form of the Hessenberg matrix H typename internal::conditional::IsComplex, ComplexSchur, RealSchur >::type schurofH; @@ -426,7 +426,7 @@ int DGMRES<_MatrixType, _Preconditioner>::dgmresComputeDeflationData(const Matri schurofH.computeFromHessenberg(m_Hes.topLeftCorner(it,it), matrixQ, computeU); ComplexVector eig(it); - Matrixperm(it); + Matrixperm(it); eig = this->schurValues(schurofH); // Reorder the absolute values of Schur values diff --git a/unsupported/Eigen/src/IterativeSolvers/GMRES.h b/unsupported/Eigen/src/IterativeSolvers/GMRES.h index 2cfa60140..23bc07d61 100644 --- a/unsupported/Eigen/src/IterativeSolvers/GMRES.h +++ b/unsupported/Eigen/src/IterativeSolvers/GMRES.h @@ -257,7 +257,7 @@ template< typename _MatrixType, typename _Preconditioner> class GMRES : public IterativeSolverBase > { typedef IterativeSolverBase Base; - using Base::mp_matrix; + using Base::matrix; using Base::m_error; using Base::m_iterations; using Base::m_info; @@ -313,7 +313,7 @@ public: m_error = Base::m_tolerance; typename Dest::ColXpr xj(x,j); - if(!internal::gmres(mp_matrix, b.col(j), xj, Base::m_preconditioner, m_iterations, m_restart, m_error)) + if(!internal::gmres(matrix(), b.col(j), xj, Base::m_preconditioner, m_iterations, m_restart, m_error)) failed = true; } m_info = failed ? NumericalIssue diff --git a/unsupported/Eigen/src/IterativeSolvers/MINRES.h b/unsupported/Eigen/src/IterativeSolvers/MINRES.h index 84e491fa1..839025591 100644 --- a/unsupported/Eigen/src/IterativeSolvers/MINRES.h +++ b/unsupported/Eigen/src/IterativeSolvers/MINRES.h @@ -198,7 +198,7 @@ namespace Eigen { { typedef IterativeSolverBase Base; - using Base::mp_matrix; + using Base::matrix; using Base::m_error; using Base::m_iterations; using Base::m_info; @@ -237,21 +237,31 @@ namespace Eigen { template 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&, - SparseSelfAdjointView, UpLo> - >::type MatrixWrapperType; - + RowMajorWrapper, + typename MatrixWrapper::template ConstSelfAdjointViewReturnType::Type + >::type SelfAdjointWrapper; + m_iterations = Base::maxIterations(); m_error = Base::m_tolerance; - + RowMajorWrapper row_mat(matrix()); for(int j=0; j