<commit_before>#ifndef _DG_LGMRES_
#define _DG_LGMRES_

#include <iostream>
#include <cstring>
#include <cmath>
#include <algorithm>

#include "blas.h"
#include "functors.h"
/*!@file
 * LGMRES class
 *
 * @author <NAME>
 */

namespace dg{

/**
* @brief Functor class for the preconditioned LGMRES method to solve
* \f[ Ax=b\f]
*
* @ingroup invert
*
* @note GMRES is a method for solving non-symmetrical linear systems.
* LGMRES is a modification of restarted GMRES that aims to improve convergence.
* This implementation is adapted from:
* https://github.com/KellyBlack/GMRES/blob/master/GMRES.h
* with LGMRES elements from
* https://github.com/haranjackson/NewtonKrylov/blob/master/lgmres.cpp
* A paper can be found at
* https://www.cs.colorado.edu/~jessup/SUBPAGES/PS/lgmres.pdf
*
*/
template< class ContainerType>
class LGMRES
{
  public:
    using container_type = ContainerType;
    using value_type = dg::get_value_type<ContainerType>; //!< value type of the ContainerType class
    ///@brief Allocate nothing, Call \c construct method before usage
    LGMRES(){}
    ///@copydoc construct()
    LGMRES( const ContainerType& copyable, unsigned max_outer, unsigned max_inner, unsigned Restarts){
        construct(copyable, max_outer, max_inner, Restarts);
    }
    ///@brief Set the number of restarts
    ///@param new_Restarts New maximum number of restarts
    void set_max( unsigned new_Restarts) {numberRestarts = new_Restarts;}
    ///@brief Get the current maximum number of iterations
    ///@return the current maximum
    unsigned get_numberRestarts() const {return numberRestarts;}
    /**
     * @brief Allocate memory for the preconditioned LGMRES method
     *
     * @param copyable A ContainerType must be copy-constructible from this
     * @param max_inner Maximum number of vectors to be saved in gmres. Usually 30 seems to be a decent number.
     * @param max_outer Maximum number of solutions saved for restart. Usually 3-10 seems to be a good number.
     * @param Restarts Maximum number of restarts. This can be set high just in case. Like e.g. gridsize/max_outer.
     */
    void construct(const ContainerType& copyable, unsigned max_outer, unsigned max_inner, unsigned Restarts){
        outer_k = max_outer;
        inner_m = max_inner;
        numberRestarts = Restarts;
        krylovDimension = inner_m + outer_k;
        //Declare Hessenberg matrix
        for(unsigned i = 0; i < krylovDimension+1; i++){
            H.push_back(std::vector<value_type>());
            for(unsigned j = 0; j < krylovDimension; j++){
                H[i].push_back(0);
            }
        }
        //Declare givens rotation matrix
        for(unsigned i = 0; i < krylovDimension+1; i++){
            givens.push_back(std::vector<value_type>());
            for(unsigned j = 0; j < 2; j++){
                givens[i].push_back(0);
            }
        }
        //Declare s that minimizes the residual... something like that.
        //s(krylovDimension+1);
        s.assign(krylovDimension+1,0);

        //The residual which will be used to calculate the solution.
        V.assign(krylovDimension+1,copyable);
        W.assign(krylovDimension,copyable);
        //In principle we don't need this many... but just to be on board with the algorithm
        outer_v.assign(outer_k,copyable);
        z = copyable;
        dx = copyable;
        residual = copyable;
    }

    /**
     * @brief Solve \f$ Ax = b\f$ using a preconditioned LGMRES method
     *
     * The iteration stops if \f$ ||Ax||_S < \epsilon( ||b||_S + C) \f$ where \f$C\f$ is
     * the absolute error in units of \f$ \epsilon\f$ and \f$ S \f$ defines a square norm
     * @param A A matrix
     * @param x Contains an initial value on input and the solution on output.
     * @param b The right hand side vector. x and b may be the same vector.
     * @param P The preconditioner to be used
     * @param S (Inverse) Weights used to compute the norm for the error condition
     * @param eps The relative error to be respected
     * @param nrmb_correction the absolute error \c C in units of \c eps to be respected
     *
     * @return Number of iterations used to achieve desired precision
     * @copydoc hide_matrix<commit_msg>BugFix: outer_v in lgmres has wrong size<commit_after>#ifndef _DG_LGMRES_
#define _DG_LGMRES_

#include <iostream>
#include <cstring>
#include <cmath>
#include <algorithm>

#include "blas.h"
#include "functors.h"
/*!@file
 * LGMRES class
 *
 * @author <NAME>
 */

namespace dg{

/**
* @brief Functor class for the preconditioned LGMRES method to solve
* \f[ Ax=b\f]
*
* @ingroup invert
*
* @note GMRES is a method for solving non-symmetrical linear systems.
* LGMRES is a modification of restarted GMRES that aims to improve convergence.
* This implementation is adapted from:
* https://github.com/KellyBlack/GMRES/blob/master/GMRES.h
* with LGMRES elements from
* https://github.com/haranjackson/NewtonKrylov/blob/master/lgmres.cpp
* A paper can be found at
* https://www.cs.colorado.edu/~jessup/SUBPAGES/PS/lgmres.pdf
*
*/
template< class ContainerType>
class LGMRES
{
  public:
    using container_type = ContainerType;
    using value_type = dg::get_value_type<ContainerType>; //!< value type of the ContainerType class
    ///@brief Allocate nothing, Call \c construct method before usage
    LGMRES(){}
    ///@copydoc construct()
    LGMRES( const ContainerType& copyable, unsigned max_outer, unsigned max_inner, unsigned Restarts){
        construct(copyable, max_outer, max_inner, Restarts);
    }
    ///@brief Set the number of restarts
    ///@param new_Restarts New maximum number of restarts
    void set_max( unsigned new_Restarts) {numberRestarts = new_Restarts;}
    ///@brief Get the current maximum number of iterations
    ///@return the current maximum
    unsigned get_numberRestarts() const {return numberRestarts;}
    /**
     * @brief Allocate memory for the preconditioned LGMRES method
     *
     * @param copyable A ContainerType must be copy-constructible from this
     * @param max_inner Maximum number of vectors to be saved in gmres. Usually 30 seems to be a decent number.
     * @param max_outer Maximum number of solutions saved for restart. Usually 3-10 seems to be a good number.
     * @param Restarts Maximum number of restarts. This can be set high just in case. Like e.g. gridsize/max_outer.
     */
    void construct(const ContainerType& copyable, unsigned max_outer, unsigned max_inner, unsigned Restarts){
        outer_k = max_outer;
        inner_m = max_inner;
        numberRestarts = Restarts;
        krylovDimension = inner_m + outer_k;
        //Declare Hessenberg matrix
        for(unsigned i = 0; i < krylovDimension+1; i++){
            H.push_back(std::vector<value_type>());
            for(unsigned j = 0; j < krylovDimension; j++){
                H[i].push_back(0);
            }
        }
        //Declare givens rotation matrix
        for(unsigned i = 0; i < krylovDimension+1; i++){
            givens.push_back(std::vector<value_type>());
            for(unsigned j = 0; j < 2; j++){
                givens[i].push_back(0);
            }
        }
        //Declare s that minimizes the residual... something like that.
        //s(krylovDimension+1);
        s.assign(krylovDimension+1,0);

        //The residual which will be used to calculate the solution.
        V.assign(krylovDimension+1,copyable);
        W.assign(krylovDimension,copyable);
        //In principle we don't need this many... but just to be on board with the algorithm
        outer_v.assign(outer_k+1,copyable);
        z = copyable;
        dx = copyable;
        residual = copyable;
    }

    /**
     * @brief Solve \f$ Ax = b\f$ using a preconditioned LGMRES method
     *
     * The iteration stops if \f$ ||Ax||_S < \epsilon( ||b||_S + C) \f$ where \f$C\f$ is
     * the absolute error in units of \f$ \epsilon\f$ and \f$ S \f$ defines a square norm
     * @param A A matrix
     * @param x Contains an initial value on input and the solution on output.
     * @param b The right hand side vector. x and b may be the same vector.
     * @param P The preconditioner to be used
     * @param S (Inverse) Weights used to compute the norm for the error condition
     * @param eps The relative error to be respected
     * @param nrmb_correction the absolute error \c C in units of \c eps to be respected
     *
     * @return Number of iterations used to achieve desired precision
     * @copydoc hide_matrix
{'fraction_non_alphanumeric': 0.09798610293251678, 'fraction_numerical': 0.004357555058297021, 'mean_word_length': 2.816629213483146, 'pattern_counts': {'":': 0, '<': 33, '<?xml version=': 0, '>': 21, 'https://': 6, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '1596566', 'n_tokens_mistral': 2609, 'n_tokens_neox': 2562, 'n_words': 909}