diff options
author | Gael Guennebaud <g.gael@free.fr> | 2015-12-07 12:23:22 +0100 |
---|---|---|
committer | Gael Guennebaud <g.gael@free.fr> | 2015-12-07 12:23:22 +0100 |
commit | b37036afce20e902cd5191a2a985f39b1f7e22e3 (patch) | |
tree | 4c7409d679d1ecbdf55b3ec518a16264fbbb7587 | |
parent | f4ca8ad9178b5fa1b83697e1a645e55d65df5639 (diff) |
Implement wrapper for matrix-free iterative solvers
-rw-r--r-- | Eigen/src/Core/util/Meta.h | 1 | ||||
-rw-r--r-- | Eigen/src/Core/util/StaticAssert.h | 3 | ||||
-rw-r--r-- | Eigen/src/IterativeLinearSolvers/BiCGSTAB.h | 4 | ||||
-rw-r--r-- | Eigen/src/IterativeLinearSolvers/ConjugateGradient.h | 21 | ||||
-rw-r--r-- | Eigen/src/IterativeLinearSolvers/IterativeSolverBase.h | 162 | ||||
-rw-r--r-- | Eigen/src/IterativeLinearSolvers/LeastSquareConjugateGradient.h | 4 | ||||
-rw-r--r-- | unsupported/Eigen/IterativeSolvers | 1 | ||||
-rw-r--r-- | unsupported/Eigen/src/IterativeSolvers/DGMRES.h | 16 | ||||
-rw-r--r-- | unsupported/Eigen/src/IterativeSolvers/GMRES.h | 4 | ||||
-rw-r--r-- | unsupported/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<BiCGSTAB<_MatrixType,_Preconditioner> > { typedef IterativeSolverBase<BiCGSTAB> 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<ConjugateGradient<_MatrixType,_UpLo,_Preconditioner> > * 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<ConjugateGradient<_MatrixType,_UpLo,_Preconditioner> > { typedef IterativeSolverBase<ConjugateGradient> 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<typename Rhs,typename Dest> void _solve_with_guess_impl(const Rhs& b, Dest& x) const { - typedef Ref<const MatrixType> MatRef; - typedef typename internal::conditional<UpLo==(Lower|Upper) && (!MatrixType::IsRowMajor) && (!NumTraits<Scalar>::IsComplex), - Transpose<const MatRef>, MatRef const&>::type RowMajorWrapper; + typedef typename Base::MatrixWrapper MatrixWrapper; + typedef typename Base::ActualMatrixType ActualMatrixType; + enum { + TransposeInput = (!MatrixWrapper::MatrixFree) + && (UpLo==(Lower|Upper)) + && (!MatrixType::IsRowMajor) + && (!NumTraits<Scalar>::IsComplex) + }; + typedef typename internal::conditional<TransposeInput,Transpose<const ActualMatrixType>, 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<UpLo==(Lower|Upper), RowMajorWrapper, - typename MatRef::template ConstSelfAdjointViewReturnType<UpLo>::Type + typename MatrixWrapper::template ConstSelfAdjointViewReturnType<UpLo>::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<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 * @@ -42,7 +164,6 @@ public: /** Default constructor. */ IterativeSolverBase() - : m_dummy(0,0), mp_matrix(m_dummy) { init(); } @@ -59,10 +180,10 @@ public: */ template<typename MatrixDerived> explicit IterativeSolverBase(const EigenBase<MatrixDerived>& 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<MatrixDerived>& 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<MatrixDerived>& 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<Scalar>::epsilon(); } - - template<typename MatrixDerived> - void grab(const EigenBase<MatrixDerived> &A) + + typedef internal::generic_matrix_wrapper<MatrixType> MatrixWrapper; + typedef typename MatrixWrapper::ActualMatrixType ActualMatrixType; + + const ActualMatrixType& matrix() const { - mp_matrix.~Ref<const MatrixType>(); - ::new (&mp_matrix) Ref<const MatrixType>(A.derived()); + return m_matrixWrapper.matrix(); } - void grab(const Ref<const MatrixType> &A) + template<typename InputType> + void grab(const InputType &A) { - if(&(A.derived()) != &mp_matrix) - { - mp_matrix.~Ref<const MatrixType>(); - ::new (&mp_matrix) Ref<const MatrixType>(A); - } + m_matrixWrapper.grab(A); } - MatrixType m_dummy; - Ref<const MatrixType> 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<LeastSquaresConjugateGradient<_MatrixType,_Preconditioner> > { typedef IterativeSolverBase<LeastSquaresConjugateGradient> 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<DGMRES<_MatrixType,_Preconditioner> > { typedef IterativeSolverBase<DGMRES> 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<DGMRES<_MatrixType,_Preconditioner> > 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<Scalar,Dynamic,Dynamic> DenseMatrix; @@ -150,7 +150,7 @@ class DGMRES : public IterativeSolverBase<DGMRES<_MatrixType,_Preconditioner> > 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<DGMRES<_MatrixType,_Preconditioner> > template<typename Dest> 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<typename RhsType, typename DestType> int dgmresApplyDeflation(const RhsType& In, DestType& Out) const; @@ -218,7 +218,7 @@ class DGMRES : public IterativeSolverBase<DGMRES<_MatrixType,_Preconditioner> > 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<DenseMatrix> 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<NumTraits<Scalar>::IsComplex, ComplexSchur<DenseMatrix>, RealSchur<DenseMatrix> >::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); - Matrix<Index,Dynamic,1>perm(it); + Matrix<StorageIndex,Dynamic,1>perm(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<GMRES<_MatrixType,_Preconditioner> > { typedef IterativeSolverBase<GMRES> 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<MINRES> 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<typename Rhs,typename Dest> 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<Scalar>::IsComplex) + }; + typedef typename internal::conditional<TransposeInput,Transpose<const ActualMatrixType>, 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<UpLo==(Lower|Upper), - Ref<const MatrixType>&, - SparseSelfAdjointView<const Ref<const MatrixType>, UpLo> - >::type MatrixWrapperType; - + RowMajorWrapper, + typename MatrixWrapper::template ConstSelfAdjointViewReturnType<UpLo>::Type + >::type SelfAdjointWrapper; + m_iterations = Base::maxIterations(); m_error = Base::m_tolerance; - + RowMajorWrapper row_mat(matrix()); for(int j=0; j<b.cols(); ++j) { m_iterations = Base::maxIterations(); m_error = Base::m_tolerance; typename Dest::ColXpr xj(x,j); - internal::minres(MatrixWrapperType(mp_matrix), b.col(j), xj, + internal::minres(SelfAdjointWrapper(row_mat), b.col(j), xj, Base::m_preconditioner, m_iterations, m_error); } |