lgmres.h

Go to the documentation of this file.

#ifndef _DG_LGMRES_
#define _DG_LGMRES_
 
#include <iostream>
#include <cstring>
#include <cmath>
#include <algorithm>
 
#include "blas.h"
#include "functors.h"
 
namespace dg{
 
template< class ContainerType>
class LGMRES
{
  public:
    using container_type = ContainerType;
    using value_type = dg::get_value_type<ContainerType>; 
    LGMRES() = default;
    LGMRES( const ContainerType& copyable, unsigned max_inner, unsigned max_outer, unsigned max_restarts):
        m_tmp(copyable),
        m_dx(copyable),
        m_residual( copyable),
        m_maxRestarts( max_restarts),
        m_inner_m( max_inner),
        m_outer_k( max_outer),
        m_krylovDimension( max_inner+max_outer)
    {
        if( m_inner_m < m_outer_k)
            std::cerr << "WARNING (LGMRES): max_inner is smaller than the restart dimension max_outer. Did you swap the constructor parameters?\n";
        //Declare Hessenberg matrix
        m_H.assign( m_krylovDimension+1, std::vector<value_type>( m_krylovDimension, 0));
        m_HH = m_H; //copy of H to be stored unaltered
        //Declare givens rotation matrix
        m_givens.assign( m_krylovDimension+1, {0,0});
        //Declare s that minimizes the residual:
        m_s.assign(m_krylovDimension+1,0);
        // m+k+1 orthogonal basis vectors:
        // k augmented pairs
        m_outer_w.assign(m_outer_k,copyable);
        m_outer_Az.assign(m_outer_k,copyable);
        m_V.assign(m_krylovDimension+1,copyable);
    }
    template<class ...Params>
    void construct( Params&& ...ps)
    {
        //construct and swap
        *this = LGMRES( std::forward<Params>( ps)...);
    }
    void set_max( unsigned new_Restarts) {m_maxRestarts = new_Restarts;}
    unsigned get_max() const {return m_maxRestarts;}
    void set_throw_on_fail( bool throw_on_fail){
        m_throw_on_fail = throw_on_fail;
    }
    const ContainerType& copyable()const{ return m_tmp;}
 
    template< class MatrixType0, class ContainerType0, class ContainerType1, class MatrixType1, class ContainerType2 >
    unsigned solve( MatrixType0&& A, ContainerType0& x, const ContainerType1& b, MatrixType1&& P, const ContainerType2& W, value_type eps = 1e-12, value_type nrmb_correction = 1);
 
    bool converged() const{
        return m_converged;
    }
 
  private:
    template <class Preconditioner, class ContainerType0>
    void Update(Preconditioner&& P, ContainerType &dx, ContainerType0 &x,
            unsigned dimension, const std::vector<std::vector<value_type>> &H,
            std::vector<value_type> &s, const std::vector<const ContainerType*> &W);
    std::vector<std::array<value_type,2>> m_givens;
    std::vector<std::vector<value_type>> m_H, m_HH;
    ContainerType m_tmp, m_dx, m_residual;
    std::vector<ContainerType> m_V, m_outer_w, m_outer_Az;
    std::vector<value_type> m_s;
    unsigned m_maxRestarts, m_inner_m, m_outer_k, m_krylovDimension;
    bool m_converged = true, m_throw_on_fail = true;
};
 
template< class ContainerType>
template < class Preconditioner, class ContainerType0>
void LGMRES<ContainerType>::Update(Preconditioner&& P, ContainerType &dx,
        ContainerType0 &x,
        unsigned dimension, const std::vector<std::vector<value_type>> &H,
        std::vector<value_type> &s, const std::vector<const ContainerType*> &W)
{
    // Solve for the coefficients, i.e. solve for c in
    // H*c=s, but we do it in place.
    for (int lupe = dimension; lupe >= 0; --lupe)
    {
        s[lupe] = s[lupe]/H[lupe][lupe];
        for (int innerLupe = lupe - 1; innerLupe >= 0; --innerLupe)
        {
            // Subtract off the parts from the upper diagonal of the matrix.
            s[innerLupe] =  DG_FMA( -s[lupe],H[innerLupe][lupe], s[innerLupe]);
        }
    }
 
    // Finally update the approximation. W_m*s
    dg::blas2::gemv( dg::asDenseMatrix( W, dimension+1), std::vector<value_type>( s.begin(), s.begin()+dimension+1), dx);
    // right preconditioner
    dg::blas2::gemv( std::forward<Preconditioner>(P), dx, m_tmp);
    dg::blas1::axpby(1.,m_tmp,1.,x);
}
 
template< class ContainerType>
template< class Matrix, class ContainerType0, class ContainerType1, class Preconditioner, class ContainerType2>
unsigned LGMRES< ContainerType>::solve( Matrix&& A, ContainerType0& x, const ContainerType1& b, Preconditioner&& P, const ContainerType2& S, value_type eps, value_type nrmb_correction)
{
    // Improvements over old implementation:
    // - Use right preconditioned system such that residual norm is available in minimization
    // - do not compute Az explicitly but save on iterations
    // - first cycle equivalent to GMRES(m+k)
    // - use weights for orthogonalization (works because in Saad book 6.29 and 6.30 are also true if V_m is unitary in the S scalar product, the Hessenberg matrix is still formed in the regular 2-norm, just define J(y) with S-norm in 6.26 and form V_m with a Gram-Schmidt process in the W-norm)
    value_type nrmb = sqrt( blas2::dot( S, b));
    value_type tol = eps*(nrmb + nrmb_correction);
    m_converged = true;
    if( nrmb == 0)
    {
        blas1::copy( 0., x);
        return 0;
    }
 
    unsigned restartCycle = 0;
    unsigned counter = 0;
    value_type rho = 1.;
    // DO NOT HOLD THESE AS PRIVATE!! MAKES BUG IN COPY!!
    std::vector<ContainerType const*> m_W, m_Vptr;
    m_W.assign(m_krylovDimension,nullptr);
    m_Vptr.assign(m_krylovDimension+1,nullptr);
    for( unsigned i=0; i<m_krylovDimension+1; i++)
        m_Vptr[i] = &m_V[i];
    do
    {
        dg::blas2::gemv(std::forward<Matrix>(A),x,m_residual);
        dg::blas1::axpby(1.,b,-1.,m_residual);
        rho = sqrt(dg::blas2::dot(S,m_residual));
        counter ++;
        if( rho < tol) //if x happens to be the solution
            return counter;
        // The first vector in the Krylov subspace is the normalized residual.
        dg::blas1::axpby(1.0/rho,m_residual,0.,m_V[0]);
 
        m_s[0] = rho;
        for(unsigned lupe=1;lupe<=m_krylovDimension;++lupe)
            m_s[lupe] = 0.0;
 
        // Go through and generate the pre-determined number of vectors for the Krylov subspace.
        for( unsigned iteration=0;iteration<m_krylovDimension;++iteration)
        {
            unsigned outer_w_count = std::min(restartCycle,m_outer_k);
            if(iteration < m_krylovDimension-outer_w_count){
                m_W[iteration] = &m_V[iteration];
                dg::blas2::gemv(std::forward<Preconditioner>(P),*m_W[iteration],m_tmp);
                dg::blas2::gemv(std::forward<Matrix>(A),m_tmp,m_V[iteration+1]);
                counter++;
            } else if( iteration < m_krylovDimension){ // size of W
                unsigned w_idx = iteration - (m_krylovDimension - outer_w_count);
                m_W[iteration] = &m_outer_w[w_idx];
                dg::blas1::copy( m_outer_Az[w_idx], m_V[iteration+1]);
            }
 
            // Get the next entry in the vectors that form the basis for the Krylov subspace.
            // Arnoldi modified Gram-Schmidt orthogonalization
            for(unsigned row=0;row<=iteration;++row)
            {
                m_HH[row][iteration] = m_H[row][iteration]
                    = dg::blas2::dot(m_V[iteration+1],S,m_V[row]);
                dg::blas1::axpby(-m_H[row][iteration],m_V[row],1.,m_V[iteration+1]);
 
            }
            m_HH[iteration+1][iteration] = m_H[iteration+1][iteration]
                = sqrt(dg::blas2::dot(m_V[iteration+1],S,m_V[iteration+1]));
            dg::blas1::scal(m_V[iteration+1],1.0/m_H[iteration+1][iteration]);
 
            // Now solve the least squares problem
            // using Givens Rotations transforming H into
            // an upper triangular matrix (see Saad Chapter 6.5.3)
            // corresponding to QR-decomposition of H
 
            // First apply previous rotations to the current matrix.
            value_type tmp = 0;
            for (unsigned row = 0; row < iteration; row++)
            {
                tmp = m_givens[row][0]*m_H[row][iteration] + // c_row
                    m_givens[row][1]*m_H[row+1][iteration];  // s_row
                m_H[row+1][iteration] = -m_givens[row][1]*m_H[row][iteration]
                    + m_givens[row][0]*m_H[row+1][iteration];
                m_H[row][iteration]  = tmp;
            }
 
            // Figure out the next Givens rotation.
            if(m_H[iteration+1][iteration] == 0.0)
            {
                // It is already upper triangular. Just leave it be....
                m_givens[iteration][0] = 1.0; // c_i
                m_givens[iteration][1] = 0.0; // s_i
            }
            else if (fabs(m_H[iteration+1][iteration]) > fabs(m_H[iteration][iteration]))
            {
                // The off diagonal entry has a larger
                // magnitude. Use the ratio of the
                // diagonal entry over the off diagonal.
                tmp = m_H[iteration][iteration]/m_H[iteration+1][iteration];
                m_givens[iteration][1] = 1.0/sqrt(1.0+tmp*tmp);
                m_givens[iteration][0] = tmp*m_givens[iteration][1];
            }
            else
            {
                // The off diagonal entry has a smaller
                // magnitude. Use the ratio of the off
                // diagonal entry to the diagonal entry.
                tmp = m_H[iteration+1][iteration]/m_H[iteration][iteration];
                m_givens[iteration][0] = 1.0/sqrt(1.0+tmp*tmp);
                m_givens[iteration][1] = tmp*m_givens[iteration][0];
            }
            // Apply the new Givens rotation on the new entry in the upper Hessenberg matrix.
            tmp = m_givens[iteration][0]*m_H[iteration][iteration] +
                  m_givens[iteration][1]*m_H[iteration+1][iteration];
            m_H[iteration+1][iteration] = -m_givens[iteration][1]*m_H[iteration][iteration] +
                  m_givens[iteration][0]*m_H[iteration+1][iteration]; // zero
            m_H[iteration][iteration] = tmp;
            // Finally apply the new Givens rotation on the s vector
            tmp = m_givens[iteration][0]*m_s[iteration] + m_givens[iteration][1]*m_s[iteration+1];
            m_s[iteration+1] = -m_givens[iteration][1]*m_s[iteration] + m_givens[iteration][1]*m_s[iteration+1];
            m_s[iteration] = tmp;
 
            rho = fabs(m_s[iteration+1]);
            if( rho < tol)
            {
                Update(std::forward<Preconditioner>(P),m_dx,x,iteration,m_H,m_s,m_W);
                return counter;
            }
        }
        Update(std::forward<Preconditioner>(P),m_dx,x,m_krylovDimension-1,m_H,m_s,m_W);
        if( m_outer_k > 1)
        {
            std::rotate(m_outer_w.rbegin(),m_outer_w.rbegin()+1,m_outer_w.rend());
            std::rotate(m_outer_Az.rbegin(),m_outer_Az.rbegin()+1,m_outer_Az.rend());
        }
        if( m_outer_k > 0)
        {
            dg::blas1::copy(m_dx,m_outer_w[0]);
            // compute A P dx
            std::vector<value_type> coeffs( m_krylovDimension+1, 0.);
            for( unsigned i=0; i<m_krylovDimension+1; i++)
            {
                coeffs[i] = 0.;
                for( unsigned k=0; k<m_krylovDimension; k++)
                    coeffs[i] = DG_FMA( m_HH[i][k],m_s[k], coeffs[i]);
            }
            dg::blas2::gemv( dg::asDenseMatrix( m_Vptr), coeffs, m_outer_Az[0]);
        }
 
        restartCycle ++;
    // Go through the requisite number of restarts.
    } while( (restartCycle < m_maxRestarts) && (rho > tol));
    if( rho > tol)
    {
        if( m_throw_on_fail)
        {
            throw dg::Fail( eps, Message(_ping_)
                <<"After "<<counter<<" LGMRES iterations");
        }
        m_converged = false;
    }
    return counter;
}
}//namespace dg
#endif