From 3b56abbbb0e231c8683b3e9ba0567f86f9ed5f26 Mon Sep 17 00:00:00 2001 From: Francesco Rizzi Date: Fri, 1 Nov 2024 16:24:48 +0100 Subject: [PATCH 1/2] wip on gmres for eigen dense --- ...inear_eigen_iterative_matrix_free_impl.hpp | 224 ++++++++++++++++++ .../solvers_linear_solver_selector_impl.hpp | 11 + .../impl/solvers_linear_traits.hpp | 18 ++ .../solvers_linear/solvers_linear_tags.hpp | 1 + .../solvers_linear/CMakeLists.txt | 6 + .../solvers_linear/eigen_gmres.cc | 172 ++++++++++++++ .../eigen_matrix_free_sparse.cc | 146 ++++++++++++ 7 files changed, 578 insertions(+) create mode 100644 include/pressio/solvers_linear/impl/solvers_linear_eigen_iterative_matrix_free_impl.hpp create mode 100644 tests/functional_small/solvers_linear/eigen_gmres.cc create mode 100644 tests/functional_small/solvers_linear/eigen_matrix_free_sparse.cc diff --git a/include/pressio/solvers_linear/impl/solvers_linear_eigen_iterative_matrix_free_impl.hpp b/include/pressio/solvers_linear/impl/solvers_linear_eigen_iterative_matrix_free_impl.hpp new file mode 100644 index 000000000..542d2a1e3 --- /dev/null +++ b/include/pressio/solvers_linear/impl/solvers_linear_eigen_iterative_matrix_free_impl.hpp @@ -0,0 +1,224 @@ +/* +//@HEADER +// ************************************************************************ +// +// solvers_linear_eigen_iterative_matrix_free_impl.hpp +// Pressio +// Copyright 2019 +// National Technology & Engineering Solutions of Sandia, LLC (NTESS) +// +// Under the terms of Contract DE-NA0003525 with NTESS, the +// U.S. Government retains certain rights in this software. +// +// Pressio is licensed under BSD-3-Clause terms of use: +// +// Redistribution and use in source and binary forms, with or without +// modification, are permitted provided that the following conditions +// are met: +// +// 1. Redistributions of source code must retain the above copyright +// notice, this list of conditions and the following disclaimer. +// +// 2. Redistributions in binary form must reproduce the above copyright +// notice, this list of conditions and the following disclaimer in the +// documentation and/or other materials provided with the distribution. +// +// 3. Neither the name of the copyright holder nor the names of its +// contributors may be used to endorse or promote products derived +// from this software without specific prior written permission. +// +// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +// "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT +// LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS +// FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE +// COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, +// INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES +// (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR +// SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) +// HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, +// STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING +// IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +// POSSIBILITY OF SUCH DAMAGE. +// +// Questions? Contact Francesco Rizzi (fnrizzi@sandia.gov) +// +// ************************************************************************ +//@HEADER +*/ + +#ifndef PRESSIO_SOLVERS_LINEAR_IMPL_SOLVERS_LINEAR_EIGEN_ITERATIVE_MATRIX_FREE_IMPL_HPP_ +#define PRESSIO_SOLVERS_LINEAR_IMPL_SOLVERS_LINEAR_EIGEN_ITERATIVE_MATRIX_FREE_IMPL_HPP_ + +#include "solvers_linear_iterative_base.hpp" +#include +#include + +namespace pressio { namespace linearsolvers{ + +template +class OperatorWrapper; + +}} + +namespace Eigen { + namespace internal { + + template + struct traits< pressio::linearsolvers::OperatorWrapper > + : public Eigen::internal::traits< + Eigen::Matrix + > + {}; + + } +} + +namespace pressio { namespace linearsolvers{ + +template +class OperatorWrapper : + public Eigen::EigenBase< + OperatorWrapper + > +{ +public: + using Scalar = typename UserDefinedOperatorType::scalar_type; + using RealScalar = Scalar; + using StorageIndex = int; + enum { + ColsAtCompileTime = Eigen::Dynamic, + MaxColsAtCompileTime = Eigen::Dynamic, + }; + + OperatorWrapper() = default; + + OperatorWrapper(UserDefinedOperatorType const & valueIn) + : m_userOperator(&valueIn) {} + + + int rows() const { return m_userOperator->rows(); } + int cols() const { return m_userOperator->cols(); } + + template + Eigen::Product, Rhs, Eigen::AliasFreeProduct> + operator*(const Eigen::MatrixBase& x) const{ + using r_t = Eigen::Product, Rhs, Eigen::AliasFreeProduct>; + return r_t(*this, x.derived()); + } + + void replace(const UserDefinedOperatorType & opIn) { + m_userOperator = &opIn; + } + + template + void applyAndAddTo(OperandT const & operand, ResultT & out) const { + // + // compute: out += operator * operand + + // out += *m_userOperator * operand; + m_userOperator->applyAndAddTo(operand, out); + } + +private: + UserDefinedOperatorType const *m_userOperator = nullptr; +}; + +}} // end namespace pressio::linearsolvers + +namespace Eigen { + namespace internal { + + template + struct generic_product_impl< + pressio::linearsolvers::OperatorWrapper, + Rhs, DenseShape, DenseShape, GemvProduct + > + : generic_product_impl_base< + pressio::linearsolvers::OperatorWrapper, Rhs, + generic_product_impl, Rhs> + > + { + using Scalar = typename Product< + pressio::linearsolvers::OperatorWrapper, + Rhs + >::Scalar; + + template + static void scaleAndAddTo( + Dest& dst, + const pressio::linearsolvers::OperatorWrapper & lhs, + const Rhs& rhs, + const Scalar& alpha) + { + // This method should implement "dst += alpha * lhs * rhs" inplace, + // however, for iterative solvers, alpha is always equal to 1, + // so let's not bother about it. + assert(alpha==Scalar(1) && "scaling is not implemented"); + EIGEN_ONLY_USED_FOR_DEBUG(alpha); + + lhs.applyAndAddTo(rhs, dst); + } + }; + } +} + +namespace pressio { namespace linearsolvers{ namespace impl{ + +template +class EigenIterativeMatrixFree + : public IterativeBase< EigenIterativeMatrixFree> +{ + +public: + using this_type = EigenIterative; + // using matrix_type = UserDefinedOperatorType; + using scalar_type = typename UserDefinedOperatorType::scalar_type; + using solver_traits = ::pressio::linearsolvers::Traits; + using op_wrapper_t = OperatorWrapper; + using native_solver_type = typename solver_traits::template eigen_solver_type; + using base_iterative_type = IterativeBase; + using iteration_type = typename base_iterative_type::iteration_type; + + static_assert(solver_traits::eigen_enabled == true, + "the native solver must be from Eigen to use in EigenIterativeMatrixFree"); + static_assert(solver_traits::direct == false, + "The native eigen solver must be iterative to use in EigenIterativeMatrixFree"); + +public: + EigenIterativeMatrixFree() = default; + + iteration_type numIterationsExecuted() const{ + return mysolver_.iterations(); + } + + scalar_type finalError() const{ + return mysolver_.error(); + } + + void resetLinearSystem(const UserDefinedOperatorType& A) + { + mysolver_.setMaxIterations(this->maxIters_); + m_wrapper.replace(A); + mysolver_.compute(m_wrapper); + } + + template + void solve(const T& b, T & y){ + mysolver_.setMaxIterations(this->maxIters_); + y = mysolver_.solve(b); + } + + template + void solve(const UserDefinedOperatorType & A, const T& b, T & y){ + this->resetLinearSystem(A); + this->solve(b, y); + } + +private: + friend base_iterative_type; + native_solver_type mysolver_ = {}; + op_wrapper_t m_wrapper; +}; + +}}} // end namespace pressio::solvers::iterarive::impl +#endif diff --git a/include/pressio/solvers_linear/impl/solvers_linear_solver_selector_impl.hpp b/include/pressio/solvers_linear/impl/solvers_linear_solver_selector_impl.hpp index 736b8d6de..3c8a708ef 100644 --- a/include/pressio/solvers_linear/impl/solvers_linear_solver_selector_impl.hpp +++ b/include/pressio/solvers_linear/impl/solvers_linear_solver_selector_impl.hpp @@ -52,6 +52,7 @@ #ifdef PRESSIO_ENABLE_TPL_EIGEN #include "solvers_linear_eigen_direct_impl.hpp" #include "solvers_linear_eigen_iterative_impl.hpp" +#include "solvers_linear_eigen_iterative_matrix_free_impl.hpp" #endif #ifdef PRESSIO_ENABLE_TPL_KOKKOS #include "solvers_linear_kokkos_direct_geqrf_impl.hpp" @@ -69,6 +70,16 @@ struct Selector{ }; #ifdef PRESSIO_ENABLE_TPL_EIGEN +template +struct Selector< + iterative::GMRES, UserDefinedOperatorType, void + > +{ + using tag_t = iterative::GMRES; + using solver_traits = ::pressio::linearsolvers::Traits; + using type = EigenIterativeMatrixFree; +}; + template struct Selector< TagType, MatrixType, diff --git a/include/pressio/solvers_linear/impl/solvers_linear_traits.hpp b/include/pressio/solvers_linear/impl/solvers_linear_traits.hpp index 3478d589a..2d39d403c 100644 --- a/include/pressio/solvers_linear/impl/solvers_linear_traits.hpp +++ b/include/pressio/solvers_linear/impl/solvers_linear_traits.hpp @@ -57,6 +57,7 @@ #include #include #include +#include #endif namespace pressio{ namespace linearsolvers{ @@ -73,6 +74,23 @@ struct Traits { #endif }; +template <> +struct Traits<::pressio::linearsolvers::iterative::GMRES> +{ + static constexpr bool direct = false; + static constexpr bool iterative = true; + +#ifdef PRESSIO_ENABLE_TPL_EIGEN + template < + typename MatrixOrOperatorT, + typename PrecT = Eigen::IdentityPreconditioner + > + using eigen_solver_type = Eigen::GMRES; + + static constexpr bool eigen_enabled = true; +#endif +}; + template <> struct Traits<::pressio::linearsolvers::iterative::CG> { diff --git a/include/pressio/solvers_linear/solvers_linear_tags.hpp b/include/pressio/solvers_linear/solvers_linear_tags.hpp index b44032a2d..345db561f 100644 --- a/include/pressio/solvers_linear/solvers_linear_tags.hpp +++ b/include/pressio/solvers_linear/solvers_linear_tags.hpp @@ -55,6 +55,7 @@ namespace iterative{ struct CG {}; struct LSCG {}; struct Bicgstab {}; +struct GMRES{}; } namespace direct{ diff --git a/tests/functional_small/solvers_linear/CMakeLists.txt b/tests/functional_small/solvers_linear/CMakeLists.txt index 4d1ffcf4a..988d6583a 100644 --- a/tests/functional_small/solvers_linear/CMakeLists.txt +++ b/tests/functional_small/solvers_linear/CMakeLists.txt @@ -2,6 +2,12 @@ if(PRESSIO_ENABLE_TPL_EIGEN) set(SRC1 ${CMAKE_CURRENT_SOURCE_DIR}/solvers_linear_eigen.cc) add_serial_utest(${TESTING_LEVEL}_solvers_linear_eigen ${SRC1}) + + set(SRC1 ${CMAKE_CURRENT_SOURCE_DIR}/eigen_gmres.cc) + add_serial_utest(${TESTING_LEVEL}_eigen_gmres ${SRC1}) + + set(SRC1 ${CMAKE_CURRENT_SOURCE_DIR}/eigen_matrix_free_sparse.cc) + add_serial_utest(${TESTING_LEVEL}_eigen_matrix_free_sparse ${SRC1}) endif() if(PRESSIO_ENABLE_TPL_KOKKOS) diff --git a/tests/functional_small/solvers_linear/eigen_gmres.cc b/tests/functional_small/solvers_linear/eigen_gmres.cc new file mode 100644 index 000000000..997bb0bb8 --- /dev/null +++ b/tests/functional_small/solvers_linear/eigen_gmres.cc @@ -0,0 +1,172 @@ + +#include +#include "pressio/solvers_linear.hpp" +#include +#include +#include +#include +#include + +class OperatorWrapper; + +namespace Eigen { + namespace internal { + template<> + struct traits + : public Eigen::internal::traits > + {}; + } +} + +class OperatorWrapper : public Eigen::EigenBase { +public: + using Scalar = double; + using RealScalar = double; + using StorageIndex = int; + enum { + ColsAtCompileTime = Eigen::Dynamic, + MaxColsAtCompileTime = Eigen::Dynamic, + }; + + OperatorWrapper() : mp_mat(0) {} + + Index rows() const { return mp_mat->rows(); } + Index cols() const { return mp_mat->cols(); } + + template + Eigen::Product + operator*(const Eigen::MatrixBase& x) const + { + using r_t = Eigen::Product; + return r_t(*this, x.derived()); + } + + void attachMyMatrix(const Eigen::MatrixXd &mat) { + mp_mat = &mat; + } + + const Eigen::MatrixXd my_matrix() const { return *mp_mat; } + + + template + void applyAndAddTo(OperandT const & operand, ResultT & out) const { + // + // compute: out += operator * operand + + out += *mp_mat * operand; + } + +private: + const Eigen::MatrixXd *mp_mat; +}; + + +namespace Eigen { + namespace internal { + + template + struct generic_product_impl< + OperatorWrapper, Rhs, DenseShape, DenseShape, GemvProduct + > + : generic_product_impl_base< + OperatorWrapper, Rhs, generic_product_impl + > + { + typedef typename Product::Scalar Scalar; + + template + static void scaleAndAddTo( + Dest& dst, + const OperatorWrapper& lhs, + const Rhs& rhs, + const Scalar& alpha) + { + // This method should implement "dst += alpha * lhs * rhs" inplace, + // however, for iterative solvers, alpha is always equal to 1, + // so let's not bother about it. + assert(alpha==Scalar(1) && "scaling is not implemented"); + EIGEN_ONLY_USED_FOR_DEBUG(alpha); + + // // // Here we could simply call dst.noalias() += lhs.my_matrix() * rhs, + // // // but let's do something fancier (and less efficient): + // for(Index i=0; i gmres; + gmres.compute(A); + gmres.setMaxIterations(8); + x = gmres.solve(b); + std::cout << "GMRES: #iterations: " << gmres.iterations() + << ", estimated error: " << gmres.error() + << std::endl; + std::cout << "my solution \n " << x << std::endl; +} + +////////////////////////////////////////////////////////////////// +////////////////////////////////////////////////////////////////// +////////////////////////////////////////////////////////////////// +////////////////////////////////////////////////////////////////// + +struct MyOperator{ + using scalar_type = double; + + MyOperator(int n){ + Eigen::MatrixXd S = Eigen::MatrixXd::Random(n,n); + m_A = S.transpose()*S; + } + + int rows() const { return m_A.rows(); } + int cols() const { return m_A.cols(); } + + template + void applyAndAddTo(OperandT const & operand, ResultT & out) const { + out = m_A * operand; + } + +private: + Eigen::MatrixXd m_A; +}; + +TEST(solvers_linear_eigen, gmres2) +{ + srand(3451677); + + int n = 10; + MyOperator A(n); + Eigen::VectorXd b(n), x(n); + b.setRandom(); + x.setZero(); + + using operator_t = MyOperator; + using tag = pressio::linearsolvers::iterative::GMRES; + using solver_t = pressio::linearsolvers::Solver; + solver_t solver; + solver.setMaxIterations(8); + solver.solve(A, b, x); + + std::cout << "my error \n " << solver.finalError() << std::endl; + std::cout << "my solution \n " << x << std::endl; +} diff --git a/tests/functional_small/solvers_linear/eigen_matrix_free_sparse.cc b/tests/functional_small/solvers_linear/eigen_matrix_free_sparse.cc new file mode 100644 index 000000000..ef587a6e6 --- /dev/null +++ b/tests/functional_small/solvers_linear/eigen_matrix_free_sparse.cc @@ -0,0 +1,146 @@ + +#include +#include "pressio/solvers_linear.hpp" +#include +#include +#include +#include +#include + + +/* + following is copied from https://eigen.tuxfamily.org/dox/group__MatrixfreeSolverExample.html + 3.4.90 (git rev 67eeba6e720c5745abc77ae6c92ce0a44aa7b7ae) + +*/ +class MatrixReplacement; +using Eigen::SparseMatrix; + +namespace Eigen { +namespace internal { + // MatrixReplacement looks-like a SparseMatrix, so let's inherits its traits: + template<> + struct traits : public Eigen::internal::traits > + {}; +} +} + +// Example of a matrix-free wrapper from a user type to Eigen's compatible type +// For the sake of simplicity, this example simply wrap a Eigen::SparseMatrix. +class MatrixReplacement : public Eigen::EigenBase { +public: + // Required typedefs, constants, and method: + typedef double Scalar; + typedef double RealScalar; + typedef int StorageIndex; + enum { + ColsAtCompileTime = Eigen::Dynamic, + MaxColsAtCompileTime = Eigen::Dynamic, + IsRowMajor = false + }; + + Index rows() const { return mp_mat->rows(); } + Index cols() const { return mp_mat->cols(); } + + template + Eigen::Product operator*(const Eigen::MatrixBase& x) const { + return Eigen::Product(*this, x.derived()); + } + + // Custom API: + MatrixReplacement() : mp_mat(0) {} + + void attachMyMatrix(const SparseMatrix &mat) { + mp_mat = &mat; + } + const SparseMatrix my_matrix() const { return *mp_mat; } + +private: + const SparseMatrix *mp_mat; +}; + + +// Implementation of MatrixReplacement * Eigen::DenseVector though a specialization of internal::generic_product_impl: +namespace Eigen { +namespace internal { + + template + struct generic_product_impl // GEMV stands for matrix-vector + : generic_product_impl_base > + { + typedef typename Product::Scalar Scalar; + + template + static void scaleAndAddTo(Dest& dst, const MatrixReplacement& lhs, const Rhs& rhs, const Scalar& alpha) + { + // This method should implement "dst += alpha * lhs * rhs" inplace, + // however, for iterative solvers, alpha is always equal to 1, so let's not bother about it. + assert(alpha==Scalar(1) && "scaling is not implemented"); + EIGEN_ONLY_USED_FOR_DEBUG(alpha); + + // Here we could simply call dst.noalias() += lhs.my_matrix() * rhs, + // but let's do something fancier (and less efficient): + for(Index i=0; i S = S1.sparseView(); + //S = S.transpose()*S; + std::cout << S << "\n"; + + MatrixReplacement A; + A.attachMyMatrix(S); + + Eigen::VectorXd b(n), x; + b.setRandom(); + + // Solve Ax = b using various iterative solver with matrix-free version: + { + Eigen::ConjugateGradient cg; + cg.compute(A); + x = cg.solve(b); + std::cout << "CG: #iterations: " << cg.iterations() << ", estimated error: " << cg.error() << std::endl; + } + + { + Eigen::BiCGSTAB bicg; + bicg.compute(A); + x = bicg.solve(b); + std::cout << "BiCGSTAB: #iterations: " << bicg.iterations() << ", estimated error: " << bicg.error() << std::endl; + } + + { + Eigen::GMRES gmres; + gmres.compute(A); + gmres.setMaxIterations(8); + x = gmres.solve(b); + std::cout << "GMRES: #iterations: " << gmres.iterations() << ", estimated error: " << gmres.error() << std::endl; + } + + { + Eigen::DGMRES gmres; + gmres.compute(A); + x = gmres.solve(b); + std::cout << "DGMRES: #iterations: " << gmres.iterations() << ", estimated error: " << gmres.error() << std::endl; + } + + { + Eigen::MINRES minres; + minres.compute(A); + x = minres.solve(b); + std::cout << "MINRES: #iterations: " << minres.iterations() << ", estimated error: " << minres.error() << std::endl; + } +} From 43c54b3dcb26c026086d25aed4038b93766cf856 Mon Sep 17 00:00:00 2001 From: Francesco Rizzi Date: Wed, 13 Nov 2024 07:06:49 -0700 Subject: [PATCH 2/2] add newton solver with gmres --- ...inear_eigen_iterative_matrix_free_impl.hpp | 52 +++-- .../solvers_linear_solver_selector_impl.hpp | 6 +- .../solvers_nonlinear/impl/functions.hpp | 17 ++ .../solvers_nonlinear/impl/internal_tags.hpp | 1 + .../solvers_nonlinear/impl/registries.hpp | 38 ++++ .../solvers_nonlinear/impl/root_finder.cpp | 186 +++++++++++++++++- .../impl/solvers_tagbased_registry.hpp | 2 + .../solvers_create_newton.hpp | 21 ++ .../ode_steppers/adjoint_logic_1.cc | 1 + .../ode_steppers/adjoint_logic_2.cc | 1 + .../ode_bdf1_strong_condition_correctness.cc | 1 + .../ode_bdf2_strong_condition_correctness.cc | 1 + .../ode_steppers/testing_apps.hpp | 2 + .../solvers_linear/eigen_gmres.cc | 2 +- .../solvers_nonlinear/CMakeLists.txt | 4 + .../newton_matrixfree_problem1_eigen.cc | 76 +++++++ .../solvers_nonlinear/qr/fixtures.hpp | 1 + 17 files changed, 376 insertions(+), 36 deletions(-) create mode 100644 tests/functional_small/solvers_nonlinear/newton_matrixfree_problem1_eigen.cc diff --git a/include/pressio/solvers_linear/impl/solvers_linear_eigen_iterative_matrix_free_impl.hpp b/include/pressio/solvers_linear/impl/solvers_linear_eigen_iterative_matrix_free_impl.hpp index 542d2a1e3..e57d4ae78 100644 --- a/include/pressio/solvers_linear/impl/solvers_linear_eigen_iterative_matrix_free_impl.hpp +++ b/include/pressio/solvers_linear/impl/solvers_linear_eigen_iterative_matrix_free_impl.hpp @@ -55,34 +55,32 @@ namespace pressio { namespace linearsolvers{ -template +template class OperatorWrapper; }} namespace Eigen { namespace internal { - - template - struct traits< pressio::linearsolvers::OperatorWrapper > + template + struct traits< pressio::linearsolvers::OperatorWrapper > : public Eigen::internal::traits< - Eigen::Matrix + Eigen::Matrix > {}; - } } namespace pressio { namespace linearsolvers{ -template +template class OperatorWrapper : public Eigen::EigenBase< - OperatorWrapper + OperatorWrapper > { public: - using Scalar = typename UserDefinedOperatorType::scalar_type; + using Scalar = typename UserDefinedLinearOperatorType::scalar_type; using RealScalar = Scalar; using StorageIndex = int; enum { @@ -92,7 +90,7 @@ class OperatorWrapper : OperatorWrapper() = default; - OperatorWrapper(UserDefinedOperatorType const & valueIn) + OperatorWrapper(UserDefinedLinearOperatorType const & valueIn) : m_userOperator(&valueIn) {} @@ -100,27 +98,26 @@ class OperatorWrapper : int cols() const { return m_userOperator->cols(); } template - Eigen::Product, Rhs, Eigen::AliasFreeProduct> + Eigen::Product, Rhs, Eigen::AliasFreeProduct> operator*(const Eigen::MatrixBase& x) const{ - using r_t = Eigen::Product, Rhs, Eigen::AliasFreeProduct>; + using r_t = Eigen::Product< + OperatorWrapper, Rhs, Eigen::AliasFreeProduct + >; return r_t(*this, x.derived()); } - void replace(const UserDefinedOperatorType & opIn) { + void replace(const UserDefinedLinearOperatorType & opIn) { m_userOperator = &opIn; } template void applyAndAddTo(OperandT const & operand, ResultT & out) const { - // // compute: out += operator * operand - - // out += *m_userOperator * operand; m_userOperator->applyAndAddTo(operand, out); } private: - UserDefinedOperatorType const *m_userOperator = nullptr; + UserDefinedLinearOperatorType const *m_userOperator = nullptr; }; }} // end namespace pressio::linearsolvers @@ -134,8 +131,8 @@ namespace Eigen { Rhs, DenseShape, DenseShape, GemvProduct > : generic_product_impl_base< - pressio::linearsolvers::OperatorWrapper, Rhs, - generic_product_impl, Rhs> + pressio::linearsolvers::OperatorWrapper, Rhs, + generic_product_impl, Rhs> > { using Scalar = typename Product< @@ -164,17 +161,18 @@ namespace Eigen { namespace pressio { namespace linearsolvers{ namespace impl{ -template +template class EigenIterativeMatrixFree - : public IterativeBase< EigenIterativeMatrixFree> + : public IterativeBase< + EigenIterativeMatrixFree + > { public: - using this_type = EigenIterative; - // using matrix_type = UserDefinedOperatorType; - using scalar_type = typename UserDefinedOperatorType::scalar_type; + using this_type = EigenIterative; + using scalar_type = typename UserDefinedLinearOperatorType::scalar_type; using solver_traits = ::pressio::linearsolvers::Traits; - using op_wrapper_t = OperatorWrapper; + using op_wrapper_t = OperatorWrapper; using native_solver_type = typename solver_traits::template eigen_solver_type; using base_iterative_type = IterativeBase; using iteration_type = typename base_iterative_type::iteration_type; @@ -195,7 +193,7 @@ class EigenIterativeMatrixFree return mysolver_.error(); } - void resetLinearSystem(const UserDefinedOperatorType& A) + void resetLinearSystem(const UserDefinedLinearOperatorType& A) { mysolver_.setMaxIterations(this->maxIters_); m_wrapper.replace(A); @@ -209,7 +207,7 @@ class EigenIterativeMatrixFree } template - void solve(const UserDefinedOperatorType & A, const T& b, T & y){ + void solve(const UserDefinedLinearOperatorType & A, const T& b, T & y){ this->resetLinearSystem(A); this->solve(b, y); } diff --git a/include/pressio/solvers_linear/impl/solvers_linear_solver_selector_impl.hpp b/include/pressio/solvers_linear/impl/solvers_linear_solver_selector_impl.hpp index 3c8a708ef..6aa39c666 100644 --- a/include/pressio/solvers_linear/impl/solvers_linear_solver_selector_impl.hpp +++ b/include/pressio/solvers_linear/impl/solvers_linear_solver_selector_impl.hpp @@ -70,14 +70,14 @@ struct Selector{ }; #ifdef PRESSIO_ENABLE_TPL_EIGEN -template +template struct Selector< - iterative::GMRES, UserDefinedOperatorType, void + iterative::GMRES, UserDefinedLinearOperatorType, void > { using tag_t = iterative::GMRES; using solver_traits = ::pressio::linearsolvers::Traits; - using type = EigenIterativeMatrixFree; + using type = EigenIterativeMatrixFree; }; template diff --git a/include/pressio/solvers_nonlinear/impl/functions.hpp b/include/pressio/solvers_nonlinear/impl/functions.hpp index 9700fc5e2..d4250209f 100644 --- a/include/pressio/solvers_nonlinear/impl/functions.hpp +++ b/include/pressio/solvers_nonlinear/impl/functions.hpp @@ -38,6 +38,23 @@ void compute_residual(RegistryType & reg, system.residual(state, r); } +#ifdef PRESSIO_ENABLE_CXX20 +template +requires NonlinearSystem +#else +template< + class RegistryType, class SystemType, + std::enable_if_t< NonlinearSystem::value, int> = 0 + > +#endif +void compute_residual(RegistryType & reg, + const SystemType & system) +{ + const auto & state = reg.template get(); + auto & r = reg.template get(); + system.residual(state, r); +} + template void compute_residual_and_jacobian(RegistryType & reg, const SystemType & system) diff --git a/include/pressio/solvers_nonlinear/impl/internal_tags.hpp b/include/pressio/solvers_nonlinear/impl/internal_tags.hpp index 13ffe11bc..17003ec7a 100644 --- a/include/pressio/solvers_nonlinear/impl/internal_tags.hpp +++ b/include/pressio/solvers_nonlinear/impl/internal_tags.hpp @@ -15,6 +15,7 @@ struct QTransposeResidualTag{}; struct NewtonTag{}; +struct MatrixFreeNewtonTag{}; struct GaussNewtonNormalEqTag{}; struct WeightedGaussNewtonNormalEqTag{}; struct LevenbergMarquardtNormalEqTag{}; diff --git a/include/pressio/solvers_nonlinear/impl/registries.hpp b/include/pressio/solvers_nonlinear/impl/registries.hpp index 6cb972bca..5f339860b 100644 --- a/include/pressio/solvers_nonlinear/impl/registries.hpp +++ b/include/pressio/solvers_nonlinear/impl/registries.hpp @@ -61,6 +61,44 @@ class RegistryNewton GETMETHOD(6) }; +template +class RegistryMatrixFreeNewtonKrylov +{ + using state_t = typename SystemType::state_type; + using r_t = typename SystemType::residual_type; + + using Tag1 = nonlinearsolvers::CorrectionTag; + using Tag2 = nonlinearsolvers::InitialGuessTag; + using Tag3 = nonlinearsolvers::ResidualTag; + using Tag4 = nonlinearsolvers::impl::SystemTag; + + state_t d1_; + state_t d2_; + r_t d3_; + SystemType const * d4_; + +public: + using linear_solver_tag = LinearSolverTag; + + RegistryMatrixFreeNewtonKrylov(const SystemType & system) + : d1_(system.createState()), + d2_(system.createState()), + d3_(system.createResidual()), + d4_(&system){} + + template + static constexpr bool contains(){ + return (mpl::variadic::find_if_binary_pred_t::value) < 4; + } + + GETMETHOD(1) + GETMETHOD(2) + GETMETHOD(3) + GETMETHOD(4) +}; + + template class RegistryGaussNewtonNormalEqs { diff --git a/include/pressio/solvers_nonlinear/impl/root_finder.cpp b/include/pressio/solvers_nonlinear/impl/root_finder.cpp index e0dd83514..dd97aeace 100644 --- a/include/pressio/solvers_nonlinear/impl/root_finder.cpp +++ b/include/pressio/solvers_nonlinear/impl/root_finder.cpp @@ -7,14 +7,13 @@ namespace nonlinearsolvers{ namespace impl{ template< - class ProblemTag, class UserDefinedSystemType, class RegistryType, class ToleranceType, class NormDiagnosticsContainerType, class DiagnosticsLoggerType, class UpdaterType> -void root_solving_loop_impl(ProblemTag /*problemTag*/, +void root_solving_loop_impl(NewtonTag /*problemTag*/, const UserDefinedSystemType & system, RegistryType & reg, Stop stopEnumValue, @@ -86,7 +85,167 @@ void root_solving_loop_impl(ProblemTag /*problemTag*/, /* stage 5 */ try{ const auto currentObjValue = - normDiagnostics[InternalDiagnostic::residualAbsoluteRelativel2Norm].getAbsolute(); + normDiagnostics[InternalDiagnostic::residualAbsoluteRelativel2Norm].getAbsolute(); + updater(reg, objective, currentObjValue); + } + catch (::pressio::eh::LineSearchStepTooSmall const &e) { + // nicely exist the solve + PRESSIOLOG_WARN(e.what()); + break; + } + catch (::pressio::eh::LineSearchObjFunctionChangeTooSmall const &e) { + // nicely exist the solve + PRESSIOLOG_WARN(e.what()); + break; + } + } +} + + /* + this wrapper is needed because the user-defined problem needs + to know what is the current nonlinear state to use for computing + the action of the jacobian + */ +template +class Wrapper +{ + using state_type = typename UserSystemT::state_type; + UserSystemT const * m_system = nullptr; + state_type const * m_currentNonLinearState = nullptr; + mutable state_type m_auxiliary; + int m_numJacRows = {}; + int m_numJacCols = {}; + +public: + // needed this by the gmres solver + using scalar_type = scalar_trait_t; + + Wrapper(UserSystemT const * system, int numRows, int numCols) + : m_system(system), + m_auxiliary(system->createState()), + m_numJacRows(numRows), m_numJacCols(numCols) + {} + + int rows() const { return m_numJacRows; } + int cols() const { return m_numJacCols; } + + void setCurrentNonLinearState(typename UserSystemT::state_type const & newState){ + m_currentNonLinearState = &newState; + } + + template + void applyAndAddTo(OperandT const & operand, ResultT & out) const { + // we need to compute out = out + jacobian * operand + // since user computes jacobian * operand, we need to do one extra step here + // so we do: + // m_auxiliary = jacobian * operand + // out += m_auxiliary + + pressio::ops::set_zero(m_auxiliary); + m_system->applyJacobian(*m_currentNonLinearState, operand, m_auxiliary); + pressio::ops::update(out, 1, m_auxiliary, 1); + } +}; + +template< + class UserDefinedSystemType, + class RegistryType, + class ToleranceType, + class NormDiagnosticsContainerType, + class DiagnosticsLoggerType, + class UpdaterType> +void root_solving_loop_impl(MatrixFreeNewtonTag /*problemTag*/, + const UserDefinedSystemType & system, + RegistryType & reg, + Stop stopEnumValue, + ToleranceType stopTolerance, + NormDiagnosticsContainerType & normDiagnostics, + const DiagnosticsLoggerType & logger, + int maxIters, + UpdaterType && updater) +{ + + using state_type = typename UserDefinedSystemType::state_type; + + auto objective = [®, &system](const state_type & stateIn){ + auto & r = reg.template get(); + system.residual(stateIn, r); + return ::pressio::ops::norm2(r); + }; + + auto mustStop = [&normDiag = std::as_const(normDiagnostics), + stopEnumValue, maxIters, stopTolerance] + (int stepCount) + { + const Diagnostic stopDiag = stop_criterion_to_public_diagnostic(stopEnumValue); + switch (stopEnumValue){ + case Stop::AfterMaxIters: + return stepCount == maxIters; + default: + if (is_absolute_diagnostic(stopDiag)){ + return normDiag[stopDiag].getAbsolute() < stopTolerance; + }else{ + return normDiag[stopDiag].getRelative() < stopTolerance; + } + }; + }; + + using wrapper_t = Wrapper; + using linear_solver_t = pressio::linearsolvers::Solver< + typename RegistryType::extended_registry::linear_solver_tag, + wrapper_t>; + linear_solver_t linearSolver; + + const auto & currentState = reg.template get(); + const auto & r = reg.template get(); + wrapper_t wrapper(&system, + pressio::ops::extent(r, 0), + pressio::ops::extent(currentState, 0)); + + int iStep = 0; + while (++iStep <= maxIters) + { + /* stage 1 */ + try{ + compute_residual(reg, system); + } + catch (::pressio::eh::ResidualEvaluationFailureUnrecoverable const &e){ + PRESSIOLOG_CRITICAL(e.what()); + throw ::pressio::eh::NonlinearSolveFailure(); + } + catch (::pressio::eh::ResidualHasNans const &e){ + PRESSIOLOG_CRITICAL(e.what()); + throw ::pressio::eh::NonlinearSolveFailure(); + } + + /* stage 2, matrix-free newton step */ + wrapper.setCurrentNonLinearState(currentState); + auto & c = reg.template get(); + linearSolver.solve(wrapper, r, c); + + // scale by -1 for sign convention + using c_t = mpl::remove_cvref_t; + using scalar_type = typename ::pressio::Traits::scalar_type; + pressio::ops::scale(c, static_cast(-1)); + + /* stage 3 */ + std::for_each(normDiagnostics.begin(), normDiagnostics.end(), + [®, iStep](auto & v){ + [[maybe_unused]] const bool isFirstIteration = iStep==1; + compute_norm_internal_diagnostics(reg, isFirstIteration, v); + }); + logger(iStep, normDiagnostics); + + /* stage 4*/ + if (mustStop(iStep)){ + PRESSIOLOG_DEBUG("nonlinsolver: stopping"); + break; + } + + /* stage 5 */ + try{ + const auto currentObjValue = + normDiagnostics[InternalDiagnostic::residualAbsoluteRelativel2Norm].getAbsolute(); updater(reg, objective, currentObjValue); } catch (::pressio::eh::LineSearchStepTooSmall const &e) { @@ -153,8 +312,9 @@ class RootFinder : public RegistryType void setStopTolerance(NormValueType value) { stopTolerance_ = value; } void setMaxIterations(int newMax) { maxIters_ = newMax; } - template - void solve(const SystemType & system, StateType & solutionInOut) + template + std::enable_if_t< std::is_same<_Tag, NewtonTag>::value > + solve(const SystemType & system, StateType & solutionInOut) { // deep copy the initial guess ::pressio::ops::deep_copy(this->template get(), solutionInOut); @@ -184,6 +344,22 @@ class RootFinder : public RegistryType } } + template + std::enable_if_t< std::is_same<_Tag, MatrixFreeNewtonTag>::value > + solve(const SystemType & system, StateType & solutionInOut) + { + // deep copy the initial guess + ::pressio::ops::deep_copy(this->template get(), solutionInOut); + + assert(updateEnValue_ == Update::Standard); + + auto extReg = reference_capture_registry_and_extend_with< + StateTag, StateType &>(*this, solutionInOut); + root_solving_loop_impl(tag_, system, extReg, stopEnValue_, stopTolerance_, + normDiagnostics_, diagnosticsLogger_, maxIters_, + DefaultUpdater()); + } + // this method can be used when the solver is applied // to the same system used for constructing it void solve(StateType & solutionInOut) diff --git a/include/pressio/solvers_nonlinear/impl/solvers_tagbased_registry.hpp b/include/pressio/solvers_nonlinear/impl/solvers_tagbased_registry.hpp index 8ea6196c5..35814cf09 100644 --- a/include/pressio/solvers_nonlinear/impl/solvers_tagbased_registry.hpp +++ b/include/pressio/solvers_nonlinear/impl/solvers_tagbased_registry.hpp @@ -98,6 +98,8 @@ class TagBasedStaticRegistryExtension extension_registry_type newReg_; public: + using extended_registry = Extendable; + template explicit TagBasedStaticRegistryExtension(Extendable & reg, CArgs && ... cargs) : reg_(reg), newReg_(std::forward(cargs)...){} diff --git a/include/pressio/solvers_nonlinear/solvers_create_newton.hpp b/include/pressio/solvers_nonlinear/solvers_create_newton.hpp index d82b82bce..640c880f6 100644 --- a/include/pressio/solvers_nonlinear/solvers_create_newton.hpp +++ b/include/pressio/solvers_nonlinear/solvers_create_newton.hpp @@ -128,5 +128,26 @@ auto create_newton_solver(const SystemType & system, (tag{}, diagnostics, system, std::forward(linSolver)); } +namespace experimental{ +template +auto create_matrixfree_newtonkrylov_solver(const SystemType & system) +{ + + using nonlinearsolvers::Diagnostic; + const std::vector diagnostics = + {Diagnostic::residualAbsolutel2Norm, + Diagnostic::residualRelativel2Norm, + Diagnostic::correctionAbsolutel2Norm, + Diagnostic::correctionRelativel2Norm}; + + using tag = nonlinearsolvers::impl::MatrixFreeNewtonTag; + using state_t = typename SystemType::state_type; + using reg_t = nonlinearsolvers::impl::RegistryMatrixFreeNewtonKrylov; + using scalar_t = scalar_trait_t< typename SystemType::state_type >; + return nonlinearsolvers::impl::RootFinder + (tag{}, diagnostics, system); +} +} + } //end namespace pressio #endif // PRESSIO_SOLVERS_NONLINEAR_SOLVERS_CREATE_NEWTON_HPP_ diff --git a/tests/functional_small/ode_steppers/adjoint_logic_1.cc b/tests/functional_small/ode_steppers/adjoint_logic_1.cc index 6953ca931..ef297949a 100644 --- a/tests/functional_small/ode_steppers/adjoint_logic_1.cc +++ b/tests/functional_small/ode_steppers/adjoint_logic_1.cc @@ -2,6 +2,7 @@ #include "pressio/ode_steppers_implicit.hpp" #include "pressio/ode_advancers.hpp" #include "random" +#include constexpr int _num_steps = 5; constexpr int _N = 17; diff --git a/tests/functional_small/ode_steppers/adjoint_logic_2.cc b/tests/functional_small/ode_steppers/adjoint_logic_2.cc index ad811d611..3e0f78e1c 100644 --- a/tests/functional_small/ode_steppers/adjoint_logic_2.cc +++ b/tests/functional_small/ode_steppers/adjoint_logic_2.cc @@ -2,6 +2,7 @@ #include "pressio/ode_steppers_implicit.hpp" #include "pressio/ode_advancers.hpp" #include "random" +#include constexpr int _num_steps = 5; constexpr int _N = 17; diff --git a/tests/functional_small/ode_steppers/ode_bdf1_strong_condition_correctness.cc b/tests/functional_small/ode_steppers/ode_bdf1_strong_condition_correctness.cc index 3630e1efd..c8332bf9d 100644 --- a/tests/functional_small/ode_steppers/ode_bdf1_strong_condition_correctness.cc +++ b/tests/functional_small/ode_steppers/ode_bdf1_strong_condition_correctness.cc @@ -1,6 +1,7 @@ #include #include "pressio/ode_steppers_implicit.hpp" +#include struct MyApp { diff --git a/tests/functional_small/ode_steppers/ode_bdf2_strong_condition_correctness.cc b/tests/functional_small/ode_steppers/ode_bdf2_strong_condition_correctness.cc index 8609e9266..5ff61eac0 100644 --- a/tests/functional_small/ode_steppers/ode_bdf2_strong_condition_correctness.cc +++ b/tests/functional_small/ode_steppers/ode_bdf2_strong_condition_correctness.cc @@ -1,6 +1,7 @@ #include #include "pressio/ode_steppers_implicit.hpp" +#include struct MyApp { diff --git a/tests/functional_small/ode_steppers/testing_apps.hpp b/tests/functional_small/ode_steppers/testing_apps.hpp index fcaa5210c..8f1b04577 100644 --- a/tests/functional_small/ode_steppers/testing_apps.hpp +++ b/tests/functional_small/ode_steppers/testing_apps.hpp @@ -2,6 +2,8 @@ #ifndef ODE_REF_APPS_FOR_TESTING_HPP_ #define ODE_REF_APPS_FOR_TESTING_HPP_ +#include + namespace pressio{ namespace ode{ namespace testing{ struct AppEigenA diff --git a/tests/functional_small/solvers_linear/eigen_gmres.cc b/tests/functional_small/solvers_linear/eigen_gmres.cc index 997bb0bb8..2f33c0a67 100644 --- a/tests/functional_small/solvers_linear/eigen_gmres.cc +++ b/tests/functional_small/solvers_linear/eigen_gmres.cc @@ -143,7 +143,7 @@ struct MyOperator{ template void applyAndAddTo(OperandT const & operand, ResultT & out) const { - out = m_A * operand; + out += m_A * operand; } private: diff --git a/tests/functional_small/solvers_nonlinear/CMakeLists.txt b/tests/functional_small/solvers_nonlinear/CMakeLists.txt index 3b32c3da0..f4852cb9a 100644 --- a/tests/functional_small/solvers_nonlinear/CMakeLists.txt +++ b/tests/functional_small/solvers_nonlinear/CMakeLists.txt @@ -8,6 +8,10 @@ if (PRESSIO_ENABLE_TPL_EIGEN) set(SRC1 ${CMAKE_CURRENT_SOURCE_DIR}/${name}.cc) add_serial_utest(${TESTING_LEVEL}_solvers_nonlinear_${name} ${SRC1}) + set(name newton_matrixfree_problem1_eigen) + set(SRC1 ${CMAKE_CURRENT_SOURCE_DIR}/${name}.cc) + add_serial_utest(${TESTING_LEVEL}_solvers_nonlinear_${name} ${SRC1}) + set(name newton_problem2_eigen) set(SRC1 ${CMAKE_CURRENT_SOURCE_DIR}/${name}.cc) add_serial_utest(${TESTING_LEVEL}_solvers_nonlinear_${name} ${SRC1}) diff --git a/tests/functional_small/solvers_nonlinear/newton_matrixfree_problem1_eigen.cc b/tests/functional_small/solvers_nonlinear/newton_matrixfree_problem1_eigen.cc new file mode 100644 index 000000000..06e0b1500 --- /dev/null +++ b/tests/functional_small/solvers_nonlinear/newton_matrixfree_problem1_eigen.cc @@ -0,0 +1,76 @@ + +#include +#include "pressio/solvers_linear.hpp" +#include "pressio/solvers_nonlinear_newton.hpp" + +struct Problem1MatrixFree +{ + using state_type = Eigen::VectorXd; + using residual_type = state_type; + +private: + mutable Eigen::SparseMatrix J_; + +public: + state_type createState() const { + state_type a(2); + a.setZero(); + return a; + } + + residual_type createResidual() const { + residual_type a(2); + a.setZero(); + return a; + } + + template + void applyJacobian(state_type const & x, + OperandT const & operand, + ResultT & out) const + { + // here we compute J but this is just for convenience here, + // we are only interested in the jacobian action + J_.resize(2,2); + J_.coeffRef(0, 0) = 3.0*x(0)*x(0); + J_.coeffRef(0, 1) = 1.0; + J_.coeffRef(1, 0) = -1.0; + J_.coeffRef(1, 1) = 3.0*x(1)*x(1); + + out = J_ * operand; + } + + void residual(const state_type& x, + residual_type& res) const + { + res(0) = x(0)*x(0)*x(0) + x(1) - 1.0; + res(1) = -x(0) + x(1)*x(1)*x(1) + 1.0; + } +}; + +TEST(solvers_nonlinear, problem1MatrixFree) +{ + pressio::log::initialize(pressio::logto::terminal); + pressio::log::setVerbosity({pressio::log::level::debug}); + + using namespace pressio; + using problem_t = Problem1MatrixFree; + using state_t = typename problem_t::state_type; + + problem_t P; + state_t y(2); + using tag = pressio::linearsolvers::iterative::GMRES; + auto nonLinSolver = experimental::create_matrixfree_newtonkrylov_solver(P); + + y(0) = 0.001; + y(1) = 0.0001; + + nonLinSolver.solve(y); + const auto e1 = std::abs(y(0) - (1.)); + const auto e2 = std::abs(y(1) - (0.)); + std::cout << y << std::endl; + ASSERT_TRUE(e1<1e-8); + ASSERT_TRUE(e2<1e-8); + + pressio::log::finalize(); +} diff --git a/tests/functional_small/solvers_nonlinear/qr/fixtures.hpp b/tests/functional_small/solvers_nonlinear/qr/fixtures.hpp index 96ae82b2b..f5246075e 100644 --- a/tests/functional_small/solvers_nonlinear/qr/fixtures.hpp +++ b/tests/functional_small/solvers_nonlinear/qr/fixtures.hpp @@ -5,6 +5,7 @@ #include #include "pressio/solvers_nonlinear.hpp" #include "qr_r9c4_gold.hpp" +#include #ifdef PRESSIO_ENABLE_EPETRA #include "Epetra_MpiComm.h"