From ffa47baf2be18e38719da935d087bbc83cf609c9 Mon Sep 17 00:00:00 2001 From: formaggia Date: Fri, 20 Sep 2024 14:54:46 +0200 Subject: [PATCH] Bettered fixed point algorithm and added Anderson acceleration --- .../src/FixedPointSolver/Accelerators.hpp | 258 +++++++++++++++--- .../FixedPointSolver/FixedPointIteration.hpp | 76 ++++-- .../src/FixedPointSolver/FixedPointTraits.hpp | 69 ++++- Examples/src/FixedPointSolver/README.md | 11 +- .../src/FixedPointSolver/main_FixedPoint.cpp | 70 +++-- 5 files changed, 397 insertions(+), 87 deletions(-) diff --git a/Examples/src/FixedPointSolver/Accelerators.hpp b/Examples/src/FixedPointSolver/Accelerators.hpp index 2ffb4ee13..d128df266 100644 --- a/Examples/src/FixedPointSolver/Accelerators.hpp +++ b/Examples/src/FixedPointSolver/Accelerators.hpp @@ -1,4 +1,4 @@ - /* +/* * Accelerators.hpp * * Created on: Jan 25, 2020 @@ -10,6 +10,7 @@ #include "FixedPointTraits.hpp" #include "VectorInterface.hpp" #include +#include #include #include //! \file Accelerators.hpp Different accelerators for fixed point iteration @@ -18,8 +19,9 @@ namespace apsc //! An enum that stores the different options made available enum FPAcceleratorId { - NoAcceleration, - ASecant + None, + ASecant, + Anderson }; //! The simple iterator @@ -27,43 +29,71 @@ enum FPAcceleratorId * This class is very thin since it only defines a call operator that * given a new iterate returns the application of the iteration function. This * class is called noAcceleration, since it just implements the simple fixed - * point iteration step + * point iteration step. * - * \tparam IteratorFunction. The iterator function. Note that at this level the - * iterator function is used only as a trait to extract the type used for the - * arguments. It must conform to a std::function! + * It is also the base class for all the other accelerators + * + * \tparam ARG. The iterator function. + * It must be equal to the type of the argument of the fixed point solver + * possible values are `FixedPointArgumentType::VECTOR` and + * `FixedPointArgumentType::EIGEN` * */ -template -class NoAccelerator +template class AcceleratorBase { public: using ArgumentType = typename FixedPointTraits::ArgumentType; using IterationFunction = typename FixedPointTraits::IterationFunction; - NoAccelerator(IterationFunction const &I) : M_I{I} {}; //!< Does nothing + explicit AcceleratorBase(IterationFunction const &I) + : M_I{I} {}; // constructor //! Call operator that returns the accelerated iterate - ArgumentType - operator()(ArgumentType const &x) + ArgumentType virtual + operator()(ArgumentType const &x) const { return M_I(x); } - //! Puts the accelerator to the original setting - void reset(){}; + //! Puts the accelerator back to the original setting + virtual void reset() const {}; //! Returns the identifier of the accelerator. Useful to implement different //! choices run-rime - constexpr FPAcceleratorId - getId() + virtual FPAcceleratorId + getId() const { - return NoAcceleration; + return id; } - //! The identifier is an attribute of the class. It may be tested on the type. - static constexpr FPAcceleratorId id = NoAcceleration; + virtual ~AcceleratorBase() = default; + static constexpr FPAcceleratorId id = None; protected: - //!@note I decided to store the function (composition). An alternative is - //!aggregation (useful if IterationFunction is a big object) + //!@note I decided to aggregate the iterator function (useful if + //! IterationFunction is a big object) + //! and not to compose it. This is a design choice. It make the class less + //! flexible but more efficient. Less flexible since it cannot have a default + //! constructors, since a reference must be bound to an object at construction + //! time. More efficient since the iterator function is not copied at each + //! call and can be shared with the fixed-point solver without useless + //! replications const IterationFunction &M_I; + // Solution with composition + // IteratorFunction M_I; }; +/* +An Alias for the identity accelerator (no acceleration) +*/ +/** + * @brief Alias template for a fixed-point solver without acceleration. + * + * This alias template defines a type `NoAcceleration` which is a specialization + * of `AcceleratorBase` for the given fixed-point argument type `ARG`. It + * represents a fixed-point solver that does not use any acceleration + * techniques. + * + * @tparam ARG The type of the argument for the fixed-point solver. + * possible values are `FixedPointArgumentType::VECTOR` and + * `FixedPointArgumentType::EIGEN` + */ +template +using NoAcceleration = AcceleratorBase; //! Alternate secant method /*! @@ -79,31 +109,31 @@ class NoAccelerator * 271–285. H. Fang, Y. Saad, Two classes of multisecant methods for nonlinear * accel- eration, Numerical Linear Algebra with Applications 16 (2009) 197–221. * - * \tparam IteratorFunction. The iterator function. Note that at this level the + * @tparam ARG. The iterator function. Note that at this level the * iterator function is used only as a trait to extract the type used for the * arguments. No other requirements are made * */ template -class ASecantAccelerator +class ASecantAccelerator : public AcceleratorBase { public: - using ArgumentType = typename FixedPointTraits::ArgumentType; - using IterationFunction = typename FixedPointTraits::IterationFunction; - ASecantAccelerator(IterationFunction const &I) : M_I{I} {}; //!< Does nothing + using ArgumentType = typename AcceleratorBase::ArgumentType; + using IterationFunction = typename AcceleratorBase::IterationFunction; + using AcceleratorBase::AcceleratorBase; // inherit constructor //! Call operator that returns the accelerated iterate - ArgumentType operator()(ArgumentType const &x); + ArgumentType operator()(ArgumentType const &x) const override; //! Put back the accelerator to initial setting void - reset() + reset() const override { length = 0; firstTime = true; } //! Returns the identifier of the accelerator. Useful to implement different //! choices run-rime - constexpr FPAcceleratorId - getId() + FPAcceleratorId + getId() const override { return id; } @@ -111,19 +141,120 @@ class ASecantAccelerator static constexpr FPAcceleratorId id = ASecant; protected: - const IterationFunction &M_I; - ArgumentType deltaXOld; - ArgumentType xOld; - ArgumentType phiOld; - size_t length{0}; - ArgumentType xNew; - ArgumentType deltaXNew; - bool firstTime = true; + mutable ArgumentType deltaXOld; + mutable ArgumentType xOld; + mutable ArgumentType phiOld; + mutable size_t length{0}; + mutable ArgumentType xNew; + mutable ArgumentType deltaXNew; + mutable bool firstTime = true; +}; + +/** + * @class AndersonAccelerator + * @brief Implements the Anderson acceleration method for fixed-point + * iterations. + * + * This class inherits from AcceleratorBase and provides an implementation of + * the Anderson acceleration technique, which is used to accelerate the + * convergence of fixed-point iterations. + * + * @tparam ARG The type of the argument for the fixed-point iteration. + * + * @note The Anderson acceleration requires at least two internal iterates to + * function properly. It can only used Eigen vectors as arguments. + */ +template + requires apsc::is_same_argument_type_v +class AndersonAccelerator : public AcceleratorBase +{ +public: + using ArgumentType = typename AcceleratorBase::ArgumentType; + using ReturnType = ArgumentType; + using IterationFunction = typename AcceleratorBase::IterationFunction; + // inherit base constructor + using AcceleratorBase::AcceleratorBase; + //! @brief Constructor for the Anderson acceleration with a maximum number of + //! stored vectors indicated + //! @param I The iterator function + //! @param m The maximum number of stored vectors + explicit AndersonAccelerator(IterationFunction const &I, unsigned int m) + : AcceleratorBase(I), m_max(m) + {} + + //! Set the maximum number of stored vectors + void + setM_max(int m) + { + m_max = m; + } + //! @brief Override the call operator + //! @param x The new iterate + //! @return The accelerated iterate + ArgumentType + operator()(ArgumentType const &x) const override + { + ArgumentType f_x = this->M_I(x); +#ifndef NDEBUG + if(m_max < 2u) + throw std::runtime_error("Anderson acceleration needs at least two " + "internal iterates, set m>2\n"); +#endif + // Update the history + if(m_history.size() == this->m_max) + { // eliminate last iterate + m_history.pop_front(); + f_history.pop_front(); + } + m_history.push_back(f_x); + f_history.push_back(f_x - x); + if(m_history.size() < 2u) + { // first or second time + return f_x; // do nothing, not enough data + } + auto const localSize = m_history.size(); + // Form the matrices F and G + Eigen::MatrixXd F(x.size(), localSize - 1u); + Eigen::MatrixXd G(x.size(), localSize - 1u); + for(std::size_t i = 0u; i < localSize - 1u; ++i) + { + F.col(i) = f_history[i + 1u] - f_history[i]; + G.col(i) = m_history[i + 1u] - m_history[i]; + } + + // Solve the least squares problem + Eigen::VectorXd gamma = F.colPivHouseholderQr().solve(f_history.back()); + + // Compute the new iterate + + return f_x - G * gamma; + } + + // Override the reset method + void + reset() const override + { + m_history.clear(); + f_history.clear(); + } + + FPAcceleratorId + getId() const override + { + return id; + } + + static constexpr FPAcceleratorId id = Anderson; + +private: + unsigned int m_max = 10u; + mutable std::deque m_history; + mutable std::deque f_history; }; template typename ASecantAccelerator::ArgumentType -ASecantAccelerator::operator()(ArgumentType const &x) +ASecantAccelerator::operator()(ArgumentType const &x) const { using namespace apsc::internals; if(firstTime) @@ -166,6 +297,55 @@ ASecantAccelerator::operator()(ArgumentType const &x) } return xNew; } + +/** + * @brief Factory function to create accelerators based on a string identifier. + * + * This function takes a string identifier and returns a unique pointer to an + * instance of the corresponding accelerator class. + * + * @tparam ARG The type of the argument for the fixed-point solver. + * Possible values are `VECTOR` and `EIGEN`. + * @param id The string identifier of the accelerator. + * @param I The iteration function. + * @param m Optional parameter for the Anderson accelerator indicating the + * maximum number of stored vectors. + * @return std::unique_ptr> A unique pointer to the created + * accelerator. + */ +template +std::unique_ptr> +createAccelerator(const std::string &id, + typename AcceleratorBase::IterationFunction const &I, + unsigned int m = 10) +{ + if(id == "NoAcceleration") + { + return std::make_unique>(I); + } + else if(id == "ASecant") + { + return std::make_unique>(I); + } + else if(id == "Anderson") + { + if constexpr(not apsc::is_same_argument_type_v< + ARG, FixedPointArgumentType::EIGEN>) + { + throw std::invalid_argument("Anderson acceleration requires Eigen " + "vectors as arguments"); + } + else + { + return std::make_unique>(I, m); + } + } + else + { + throw std::invalid_argument("Unknown accelerator type: " + id); + } +} + } // end namespace apsc #endif /* SRC_NONLINSYSSOLVER_ACCELERATORS_HPP_ */ diff --git a/Examples/src/FixedPointSolver/FixedPointIteration.hpp b/Examples/src/FixedPointSolver/FixedPointIteration.hpp index 47925929f..e003c2909 100644 --- a/Examples/src/FixedPointSolver/FixedPointIteration.hpp +++ b/Examples/src/FixedPointSolver/FixedPointIteration.hpp @@ -29,35 +29,73 @@ struct FixedPointOptions //! A class for (possibly) accelerated fixed point iterations /*! - * \tparam IterationFunction a system of equation that obeys the interface - * defined in FixedPointTraits::ArgumentType \tparam Accelerator A template - * template parameter that identify a procedure to accelerate convergence. + * @param ARG the type used for the iterator function arguments (EIGEN or + * VECTOR) * */ -template class Accelerator = NoAccelerator> +template class FixedPointIteration { public: //! The type of the argument and return value using ArgumentType = typename FixedPointTraits::ArgumentType; - using AcceleratorType = Accelerator; + using AcceleratorType = std::unique_ptr>; using IterationFunction = typename FixedPointTraits::IterationFunction; - //! Constructor that takes the iterator function and (optionally) the - //! convergence options - template - FixedPointIteration(IF && ifun = IterationFunction(), - FixedPointOptions opt = FixedPointOptions{}) - : phi{std::forward(ifun)}, options{opt}, accelerator(phi) + + /** + * @brief Constructs a FixedPointIteration object. + * + * @tparam IF Type of the iteration function. + * @param ifun The iteration function. + * @param opt Options for the fixed-point iteration, defaults to + * FixedPointOptions. + * @param AccelarationType Type of acceleration to use, defaults to + * "NoAcceleration". possible other values are "ASecant" and "Anderson". + * @param m Parameter for the accelerator (used only for Anderson), defaults + * to 10. + * @note I am using the trick of forward reference to avoid copies of the + * iteration function in case the argument can be moved. This is a C++11 + * feature. + */ + template + explicit FixedPointIteration(IF &&ifun, + FixedPointOptions opt = FixedPointOptions{}, + std::string AccelarationType = "NoAcceleration", + unsigned int m = 10) + : phi{std::forward(ifun)}, options{opt}, + accelerator(apsc::createAccelerator(AccelarationType, phi, m)) {} - //! Default constructor. IterationFunction is not set. Options are the default - FixedPointIteration() : accelerator(phi){}; + //! This class is not default constructible. Here I state it explicitely + FixedPointIteration() = delete; + //! This class is not copy constructible. Here I state it explicitely + /*! + * @note This is a design choice. The accelerator is a unique pointer to a + * polymprphic object that has not been designed to be copied. To make the + * class copyable we need a strong revision of the design. + */ + FixedPointIteration(FixedPointIteration const &) = delete; + //! This class is not copy assigneable. Here I state it explicitly + + FixedPointIteration &operator=(FixedPointIteration const &) = delete; + //! This class is move constructible. Here I state it explicitly + FixedPointIteration(FixedPointIteration &&) = default; + //! This class is move assignable. Here I state it explicitly + FixedPointIteration &operator=(FixedPointIteration &&) = default; + //! Getter for options FixedPointOptions getOptions() const { return options; } + /* + change acceleration method + */ + void + setAcceleration(std::string AccelarationType, unsigned int m = 10) + { + accelerator = apsc::createAccelerator(AccelarationType, phi, m); + } //! setter for options void setOptions(const FixedPointOptions &options) @@ -67,13 +105,13 @@ class FixedPointIteration //! To set, or change the iteration function template void - setIterationFunction(IterationFunction&& ifun) + setIterationFunction(IterationFunction &&ifun) { phi = std::forward(ifun); } //! Allow to access the accelerator /*! - * In case we need to set its state. + * In case we need to change its state. */ auto & getAccelerator() @@ -92,7 +130,7 @@ class FixedPointIteration * */ auto - compute(ArgumentType const &x0) + compute(ArgumentType const &x0) const { double currentDistance{std::numeric_limits::max()}; std::size_t iter{0}; @@ -100,7 +138,7 @@ class FixedPointIteration ArgumentType previous = x0; while(iter < options.maxIter && currentDistance > options.tolerance) { - current = this->accelerator(previous); + current = std::invoke(*accelerator, previous); currentDistance = std::sqrt(internals::squaredDistance(current, previous)); previous = current; @@ -112,7 +150,7 @@ class FixedPointIteration std::cout << std::endl; #endif } - this->accelerator.reset(); + accelerator->reset(); return std::make_tuple(current, iter, currentDistance, (iter < options.maxIter)); } diff --git a/Examples/src/FixedPointSolver/FixedPointTraits.hpp b/Examples/src/FixedPointSolver/FixedPointTraits.hpp index dd4287e68..ddf6de12a 100644 --- a/Examples/src/FixedPointSolver/FixedPointTraits.hpp +++ b/Examples/src/FixedPointSolver/FixedPointTraits.hpp @@ -8,17 +8,52 @@ #ifndef NONLINSYSSOLVER_FIXEDPOINTTRAITS_HPP_ #define NONLINSYSSOLVER_FIXEDPOINTTRAITS_HPP_ #include +#include #include +/** + * @file FixedPointTraits.hpp + * @brief This file contains the definition of traits and utilities for + * fixed-point iteration methods. + * + * The file defines the following: + * - An enumeration `FixedPointArgumentType` to specify the type of argument + * used in fixed-point iterations. + * - A primary template `FixedPointTraits` which is specialized for different + * `FixedPointArgumentType` values. + * - Specializations of `FixedPointTraits` for `FixedPointArgumentType::VECTOR` + * and `FixedPointArgumentType::EIGEN`. + * - A custom trait `is_same_argument_type` to check if two + * `FixedPointArgumentType` values are the same. + * - A constexpr variable template `is_same_argument_type_v` for ease of use of + * the `is_same_argument_type` trait. + * + * The `FixedPointTraits` template provides the following type aliases for each + * specialization: + * - `ArgumentType`: The type of the argument used in the iteration function. + * - `ReturnType`: The type returned by the iteration function. + * - `IterationFunction`: A `std::function` type representing the iteration + * function. + * + * The `is_same_argument_type` trait and `is_same_argument_type_v` variable + * template are used to perform compile-time checks on the equality of + * `FixedPointArgumentType` values. + */ namespace apsc { +/*! + * @brief An enumeration to specify the type of argument used in fixed-point + iterations. We can have two possible way of defining the argument for the + fixed-point iteration: + - `VECTOR`: The argument is a `std::vector`. + - `EIGEN`: The argument is an `Eigen::Matrix`. +*/ enum class FixedPointArgumentType { VECTOR, EIGEN }; -template -struct FixedPointTraits +template struct FixedPointTraits {}; template <> struct FixedPointTraits @@ -35,6 +70,36 @@ template <> struct FixedPointTraits using IterationFunction = std::function; }; +/*! +@brief A type trait to check if two `FixedPointArgumentType` values are the +same. +@tparam Value1 The first `FixedPointArgumentType` value. +@tparam Value2 The second `FixedPointArgumentType` value. +@details It is an example of template metaprogramming that extends the +`std::is_same` trait to work with `FixedPointArgumentType` values. Unfortunately +(so far) C++ deas not provide a specialization of `std::is_same` for +enumerators. +*/ +template +struct is_same_argument_type : std::false_type +{}; + +/* +@brief This specialization is activated if the two `FixedPointArgumentType` +values are the same. +*/ +template +struct is_same_argument_type : std::true_type +{}; + +/*! +@brief A helper variable template to check if two `FixedPointArgumentType` +values are the same. +*/ +template +constexpr bool is_same_argument_type_v = + is_same_argument_type::value; + } // namespace apsc #endif /* NONLINSYSSOLVER_FIXEDPOINTTRAITS_HPP_ */ diff --git a/Examples/src/FixedPointSolver/README.md b/Examples/src/FixedPointSolver/README.md index 30315203b..7ae776d9e 100644 --- a/Examples/src/FixedPointSolver/README.md +++ b/Examples/src/FixedPointSolver/README.md @@ -14,8 +14,8 @@ with the following features to use some standard algorithms that allow a parallel version, even if in this version I do not use the parallel capability. * We can add an accelerator, and some are provided in `Accelerators.hpp`. An accelerator takes as input the last iterate and returns a corrected value which hopefully - provides a better convergent sequence. -* Test of convergence based on the distance of two consecutive iterates (which is classical for fixed oint iterations). + provides a better convergent sequence. I have made the choice to implement the accelerators as a policy, so that they can be developed independently of the main class. Moreover, I decided to agglomerate the iterator function in the accelerator, to save memory in case phi is a callable object of big size. With agglomeration the accelerator and the main class share the same memory for the iterator function. +* Test of convergence based on the distance of two consecutive iterates (which is classical for fixed point iterations). * I have used an aggregate (a `struct` called `FixedPointOptions`) to store the options, which have a default value. You may change the options by passing a modified `FixedPointOption` to the the solver class. A possibility is to make this struct internal to the `FixedPointIteration` class. Indeed, it is used only in conjunction with the use of a `FixedPointIteration`. * If you want to have a more verbose output compile with @@ -24,6 +24,12 @@ to use some standard algorithms that allow a parallel version, even if in this v make VERBOSE=yes ``` +## Accelerators implemented ## +- `NoAccelerator`: the simplest accelerator, it does nothing. +- `ASecant`: It is two level Anderson acceleration, that may be considered as the multidimensional estension of the secant method for finding the zero of $f(x)=x-\phi(x). A giid reference of the technique is *V. Eyert, A Comparative Study on Methods for Convergence Acceleration of Iterative Vector Sequences, Journal of Computational Physics 124 (1996)271–285.* or *H. Fang, Y. Saad, Two classes of multisecant methods for nonlinearacceleration, Numerical Linear Algebra with Applications 16 (2009) 197–221.* +- `Anderson`. I timplements anderson acceleration. A good reference is *D. G. Anderson, Iterative procedures for nonlinear integral equations, J. Assoc. Comput. Mach. 12 (1965) 547–560.* or *D. G. Anderson, Iterative procedures for nonlinear integral equations, J. Assoc. Comput. Mach. 12 (1965) 547–560.*. The implementation is based on the algorithm described in *D. A. Smits, Accelerating the convergence of iterative vector sequences, Numerical Linear Algebra with Applications 4 (1997) 1–30.*. The algorithm is also described in *H. Fang, Y. Saad, Two classes of multisecant methods for nonlinearacceleration, Numerical Linear Algebra with Applications 16 (2009) 197–221.*. The implementation is based on the algorithm described in *D. A. Smits, Accelerating the convergence of iterative vector sequences, Numerical Linear Algebra with Applications 4 (1997) 1–30.*. The algorithm is also described in *H. Fang, Y. Saad, Two classes of multisecant methods for nonlinearacceleration, Numerical Linear Algebra with Applications 16 (2009) 197–221.*. The implementation is based on the algorithm described in *D. A. Smits, Accelerating the convergence of iterative vector sequences, Numerical Linear Algebra with Applications 4 (1997) 1–30.*. The algorithm is also described in *H. Fang, Y. Saad, Two classes of multisecant methods for nonlinearacceleration, Numerical Linear Algebra with Applications 16 (2009) 197–221.*. The implementation is based on the algorithm described in *D. A. Smits, Accelerating the convergence of iterative vector sequences, Numerical Linear Algebra with Applications 4 (1997) 1–30.*. The algorithm is also described in *H. Fang, Y. Saad, Two classes of multisecant methods for nonlinear acceleration, Numerical Linear Algebra with Applications 16 (2009) 197–221*. The implementation is based on the algorithm described in *D. A. Smits, Accelerating the convergence of iterative vector sequences, Numerical Linear Algebra with Applications 4 (1997) 1–30*. The algorithm is also described in the previously cited reference by *H. Fang, Y. Saad*. + + # What do I learn from this example? # - The use of adapters to make your code more general, in particular to enable it to interface with objects with different public interfaces @@ -32,4 +38,5 @@ make VERBOSE=yes - The use of some standard algorithms - The use of traits to centralize the declaration of types used throughout the code and to enable the selection among different possible choices. - The use of policies (the accelerators) to implement part of an algorithm with a set of classes that can be developed separately. it is an example of the Bridge design pattern (here implemented via templates). +- The use of user defined traits to check if an enumerator has a certain value. This is useful to check if we are using EIGEN vectors as argument for the iterator function, since the Anderson Accelerator works only in that case. diff --git a/Examples/src/FixedPointSolver/main_FixedPoint.cpp b/Examples/src/FixedPointSolver/main_FixedPoint.cpp index 8cdace05e..5d4d1180f 100644 --- a/Examples/src/FixedPointSolver/main_FixedPoint.cpp +++ b/Examples/src/FixedPointSolver/main_FixedPoint.cpp @@ -9,6 +9,9 @@ #include "FixedPointTraits.hpp" #include #include +#include +#include + template void print_result(T solution) @@ -33,7 +36,7 @@ main() { using namespace apsc; using FixedPointIterator = - apsc::FixedPointIteration; + apsc::FixedPointIteration; using IterationFunction = FixedPointIterator::IterationFunction; IterationFunction phi; @@ -46,32 +49,49 @@ main() // y=0 and is unstable, the other is stable // and we converge to the second one // Try to change lambda and see what happens - double lambda = 1.5; + double lambda = 1; + std::cout << "Give a lambda>=0; <1 for convergence to (0,0.7390085)\n"; + std::cin >> lambda; + std::cout << "You have chosen lambda=" << lambda << "\n"; std::pair exact; - phi = [lambda](FixedPointIterator::ArgumentType const &x) { - return std::vector{lambda * std::sin(x[0]), std::cos(x[1])}; + using ArgumentType = FixedPointIterator::ArgumentType; + phi = [lambda](ArgumentType const &x) -> ArgumentType { + return FixedPointIterator::ArgumentType{{lambda * std::sin(x[0])}, + {std::cos(x[1])}}; }; - FixedPointIterator iterate{phi}; - FixedPointIterator::ArgumentType startingPoint{5.0, 7.0}; - std::cout << "*** WITH BASIC METHOD:\n"; + std::set accelerators{"NoAcceleration", "ASecant", "Anderson"}; + std::string accelerator; + bool ok = false; + unsigned int m{0}; + while(not ok) + { + std::cout << "Which acceleration do you want? (NoAcceleration, ASecant " + "or Anderson)" + << std::endl; + std::cin >> accelerator; + if(accelerators.find(accelerator) != accelerators.end()) + { + ok = true; + if(accelerator == "Anderson") + { + std::cout << "Give the number of stored vectors for Anderson\n"; + std::cin >> m; + std::cout << "You have chosen " << m << " vectors\n"; + } + } + else + { + std::cout << "Wrong accelerator name, retry\n"; + } + } + apsc::FixedPointOptions options; + FixedPointIterator iterate{phi, options, accelerator, m}; + + ArgumentType startingPoint{ + {{5.0}, {7.0}}}; // strange way of construction a Eigen vector + + std::cout << "*** Result:\n"; print_result(iterate.compute(startingPoint)); - // Now with acceleration and Eigen - using FixedPointIterator2 = - apsc::FixedPointIteration; - using IterationFunction2 = FixedPointIterator2::IterationFunction; - IterationFunction2 phi2; - phi2 = [lambda](FixedPointIterator2::ArgumentType const &x) { - FixedPointIterator2::ArgumentType res(2); - res[0] = lambda * std::sin(x[0]); - res[1] = std::cos(x[1]); - return res; - }; - FixedPointIterator2 iterate2{phi2}; - FixedPointIterator2::ArgumentType startingPoint2(2); - startingPoint2[0] = 5.0; - startingPoint2[1] = 7.0; - std::cout << "*** WITH SECANT ACCELERATION:\n"; - print_result(iterate2.compute(startingPoint2)); + return 0; }