eve.h

Go to the documentation of this file.

#ifndef _DG_EVE_
#define _DG_EVE_
 
#include <cmath>
#include "blas.h"
#include "functors.h"
 
namespace dg
{
 
template< class ContainerType>
class EVE
{
  public:
    using container_type = ContainerType;
    using value_type = get_value_type<ContainerType>; 
    EVE() = default;
    EVE( const ContainerType& copyable, unsigned max_iter = 100):r( copyable), p( r), ap( r), m_max_iter( max_iter) {}
    void construct( const ContainerType& copyable, unsigned max_iter = 100) {
        ap = p = r = copyable;
        m_max_iter = max_iter;
    }
    void set_max( unsigned new_max) {
        m_max_iter = new_max;
    }
    unsigned get_max() const {   return m_max_iter; }
    void set_throw_on_fail( bool throw_on_fail){
        m_throw_on_fail = throw_on_fail;
    }
    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& ev_max, value_type eps_ev = 1e-12);
  private:
    ContainerType r, p, ap;
    unsigned m_max_iter;
    bool m_throw_on_fail = true;
};
 
 
template< class ContainerType>
template< class Matrix, class ContainerType0, class ContainerType1, class Preconditioner, class ContainerType2>
unsigned EVE< ContainerType>::solve( Matrix&& A, ContainerType0& x, const ContainerType1& b, Preconditioner&& P, const ContainerType2& W, value_type& ev_max, value_type eps_ev )
{
    blas2::symv( std::forward<Matrix>(A),x,r);
    blas1::axpby( 1., b, -1., r);
    blas2::symv( std::forward<Preconditioner>(P), r, p );
    value_type nrmzr_old = blas2::dot( p, W,r);
    value_type nrmzr_new, nrmAp;
    value_type alpha = 1., alpha_inv = 1., delta = 0.;
    value_type evdash, gamma = 0., lambda, omega, beta = 0.;
    value_type ev_est = 0.;
    ev_max = 0.;
    for( unsigned i=1; i<m_max_iter; i++)
    {
        lambda = delta*alpha_inv;    // EVE!
        blas2::symv( std::forward<Matrix>(A), p, ap);
        nrmAp = blas2::dot( p, W, ap);
        alpha =  nrmzr_old/nrmAp;
        alpha_inv = nrmAp/nrmzr_old; //EVE!
        lambda+= alpha_inv;          //EVE!
        blas1::axpby( alpha, p, 1.,x);
        blas1::axpby( -alpha, ap, 1., r);
        blas2::symv(std::forward<Preconditioner>(P),r,ap);
        nrmzr_new = blas2::dot( ap, W, r);
        delta = nrmzr_new /nrmzr_old;                  // EVE!
        evdash = ev_est -lambda;                       // EVE!
        omega = sqrt( evdash*evdash +4.*beta*gamma);   // EVE!
        gamma = 0.5 *(1. -evdash /omega);              // EVE!
        ev_max += omega*gamma;                         // EVE!
        beta = delta*alpha_inv*alpha_inv;              // EVE!
        if( fabs( ev_est - ev_max) < eps_ev*ev_max)
            return i;
        blas1::axpby(1.,ap, delta, p );
        nrmzr_old=nrmzr_new;
        ev_est = ev_max;
    }
    if( m_throw_on_fail)
    {
        throw dg::Fail( eps_ev, Message(_ping_)
            <<"After "<<m_max_iter<<" EVE iterations");
    }
    return m_max_iter;
}
 
} //namespace dg
#endif //_DG_EVE_