pcg.h

Go to the documentation of this file.

#ifndef _DG_PCG_
#define _DG_PCG_
 
#include <cmath>
 
#include "blas.h"
#include "functors.h"
#include "extrapolation.h"
#include "backend/typedefs.h"
 
#include "backend/timer.h"
 
namespace dg{
 
 
template< class ContainerType>
class PCG
{
  public:
    using container_type = ContainerType;
    using value_type = get_value_type<ContainerType>; 
    PCG() = default;
    PCG( const ContainerType& copyable, unsigned max_iterations):
        r(copyable), p(r), ap(r), max_iter(max_iterations){}
    void set_max( unsigned new_max) {max_iter = new_max;}
    unsigned get_max() const {return max_iter;}
    const ContainerType& copyable()const{ return r;}
 
    void set_verbose( bool verbose){ m_verbose = verbose;}
    void set_throw_on_fail( bool throw_on_fail){
        m_throw_on_fail = throw_on_fail;
    }
 
    template<class ...Params>
    void construct( Params&& ...ps)
    {
        //construct and swap
        *this = PCG( std::forward<Params>( ps)...);
    }
    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, int test_frequency = 1);
  private:
    ContainerType r, p, ap;
    unsigned max_iter;
    bool m_verbose = false, m_throw_on_fail = true;
};
 
template< class ContainerType>
template< class Matrix, class ContainerType0, class ContainerType1, class Preconditioner, class ContainerType2>
unsigned PCG< ContainerType>::solve( Matrix&& A, ContainerType0& x, const ContainerType1& b, Preconditioner&& P, const ContainerType2& W, value_type eps, value_type nrmb_correction, int save_on_dots )
{
    // self-adjoint: apply PCG algorithm to (P 1/W) (W A) x = (P 1/W) (W b) : P' A' x = P' b'
    // This effectively just replaces all scalar products with the weighted one
    value_type nrmb = sqrt( blas2::dot( W, b));
    value_type tol = eps*(nrmb + nrmb_correction);
#ifdef MPI_VERSION
    int rank;
    MPI_Comm_rank(MPI_COMM_WORLD, &rank);
#endif //MPI
    if( m_verbose)
    {
        DG_RANK0 std::cout << "# Norm of W b "<<nrmb <<"\n";
        DG_RANK0 std::cout << "# Residual errors: \n";
    }
    if( nrmb == 0)
    {
        blas1::copy( 0., x);
        return 0;
    }
    blas2::symv( std::forward<Matrix>(A),x,r);
    blas1::axpby( 1., b, -1., r);
    if( sqrt( blas2::dot(W,r) ) < tol) //if x happens to be the solution
        return 0;
    blas2::symv( std::forward<Preconditioner>(P), r, p );
    value_type nrmzr_old = blas2::dot( p,W,r); //and store the scalar product
    value_type alpha, nrmzr_new;
    for( unsigned i=1; i<max_iter; i++)
    {
        blas2::symv( std::forward<Matrix>(A), p, ap);
        alpha =  nrmzr_old/blas2::dot( p, W, ap);
        blas1::axpby( alpha, p, 1.,x);
        blas1::axpby( -alpha, ap, 1., r);
        if( 0 == i%save_on_dots )
        {
            if( m_verbose)
            {
                DG_RANK0 std::cout << "# Absolute r*W*r "<<sqrt( blas2::dot(W,r)) <<"\t ";
                DG_RANK0 std::cout << "#  < Critical "<<tol <<"\t ";
                DG_RANK0 std::cout << "# (Relative "<<sqrt( blas2::dot(W,r) )/nrmb << ")\n";
            }
            if( sqrt( blas2::dot(W,r)) < tol)
                return i;
        }
        blas2::symv(std::forward<Preconditioner>(P),r,ap);
        nrmzr_new = blas2::dot( ap, W, r);
        blas1::axpby(1.,ap, nrmzr_new/nrmzr_old, p );
        nrmzr_old=nrmzr_new;
    }
    if( m_throw_on_fail)
    {
        throw dg::Fail( tol, Message(_ping_)
            <<"After "<<max_iter<<" PCG iterations with rtol "<<eps<<" and atol "<<eps*nrmb_correction );
    }
    return max_iter;
}
 
} //namespace dg
 
 
 
#endif //_DG_PCG_