chebyshev.h

Go to the documentation of this file.

#ifndef _DG_CHEB_
#define _DG_CHEB_
 
#include <cmath>
 
#include "blas.h"
 
namespace dg
{
 
template< class ContainerType>
class ChebyshevIteration
{
  public:
    using container_type = ContainerType;
    using value_type = get_value_type<ContainerType>; 
    ChebyshevIteration() = default;
    ChebyshevIteration( const ContainerType& copyable):
        m_ax(copyable), m_z( m_ax), m_xm1(m_ax){}
    const ContainerType& copyable()const{ return m_ax;}
 
    void construct( const ContainerType& copyable) {
        m_xm1 = m_z = m_ax = copyable;
    }
    template< class MatrixType, class ContainerType0, class ContainerType1>
    void solve( MatrixType&& A, ContainerType0& x, const ContainerType1& b,
        value_type min_ev, value_type max_ev, unsigned num_iter, bool x_is_zero = false)
    {
        if( num_iter == 0)
            return;
        assert ( min_ev < max_ev);
        value_type theta = (min_ev+max_ev)/2., delta = (max_ev-min_ev)/2.;
        value_type rhokm1 = delta/theta, rhok=0;
        if( !x_is_zero)
        {
            dg::blas1::copy( x, m_xm1); //x_{k-1}
            dg::blas2::symv( std::forward<MatrixType>(A), x, m_ax);
            dg::blas1::axpbypgz( 1./theta, b, -1./theta, m_ax, 1., x); //x_1
        }
        else
        {
            dg::blas1::copy( 0., m_xm1); //x_{k-1}
            dg::blas1::axpby( 1./theta, b, 0., x); //x_1
        }
        for ( unsigned k=1; k<num_iter; k++)
        {
            rhok = 1./(2.*theta/delta - rhokm1);
            dg::blas2::symv( std::forward<MatrixType>(A), x, m_ax);
            dg::blas1::evaluate( m_xm1, dg::equals(), PairSum(),
                             1.+rhok*rhokm1, x,
                            -rhok*rhokm1,    m_xm1,
                             2.*rhok/delta,  b,
                            -2.*rhok/delta,  m_ax
                            );
            using std::swap;
            swap( x, m_xm1);
            rhokm1 = rhok;
        }
    }
    template< class MatrixType0, class MatrixType1, class ContainerType0, class ContainerType1>
    void solve( MatrixType0&& A, ContainerType0& x, const ContainerType1& b,
            MatrixType1&& P, value_type min_ev, value_type max_ev, unsigned num_iter,
            bool x_is_zero = false)
    {
        if( num_iter == 0)
            return;
        assert ( min_ev < max_ev);
        value_type theta = (min_ev+max_ev)/2., delta = (max_ev-min_ev)/2.;
        value_type rhokm1 = delta/theta, rhok=0;
        if( !x_is_zero)
        {
            dg::blas1::copy( x, m_xm1); //x_{k-1}
            dg::blas2::symv( std::forward<MatrixType0>(A), x, m_ax);
            dg::blas1::axpby( 1., b, -1., m_ax); //r_0
            dg::blas2::symv( std::forward<MatrixType1>(P), m_ax, m_z);
            dg::blas1::axpby( 1./theta, m_z, 1., x); //x_{k-1}
        }
        else
        {
            dg::blas2::symv( std::forward<MatrixType1>(P), b, x);
            if( num_iter == 1) return;
            dg::blas1::scal( m_xm1, 0.);
            dg::blas1::scal( x, 1./theta);
        }
        for ( unsigned k=1; k<num_iter; k++)
        {
            rhok = 1./(2.*theta/delta - rhokm1);
            dg::blas2::symv( std::forward<MatrixType0>(A), x, m_ax);
            dg::blas1::axpby( 1., b, -1., m_ax); //r_k
            dg::blas2::symv( P, m_ax, m_z);
            dg::blas1::axpbypgz(
                             1.+rhok*rhokm1, x,
                             2.*rhok/delta,  m_z,
                            -rhok*rhokm1,    m_xm1
                            );
            using std::swap;
            swap( x, m_xm1);
            rhokm1 = rhok;
        }
    }
  private:
    ContainerType m_ax, m_z, m_xm1;
};
 
template<class Matrix, class ContainerType>
struct ChebyshevPreconditioner
{
    using container_type = ContainerType;
    using value_type = get_value_type<ContainerType>; 
    ChebyshevPreconditioner( Matrix op, const ContainerType& copyable, value_type ev_min,
            value_type ev_max, unsigned degree):
        m_op(op), m_ch( copyable),
        m_ev_min(ev_min), m_ev_max(ev_max), m_degree(degree){}
 
    template<class ContainerType0, class ContainerType1>
    void symv( const ContainerType0& x, ContainerType1& y)
    {
        //m_ch.solve( m_op, y, x, m_op.precond(), m_ev_min, m_ev_max, m_degree+1, true);
        m_ch.solve( m_op, y, x, m_ev_min, m_ev_max, m_degree+1, true);
    }
    private:
    Matrix m_op;
    ChebyshevIteration<ContainerType> m_ch;
    value_type m_ev_min, m_ev_max;
    unsigned m_degree;
};
 
template<class Matrix, class ContainerType>
struct ModifiedChebyshevPreconditioner
{
    using container_type = ContainerType;
    using value_type = get_value_type<ContainerType>; 
    ModifiedChebyshevPreconditioner( Matrix op, const ContainerType& copyable, value_type ev_min,
            value_type ev_max, unsigned degree):
        m_op(op), m_ax(copyable), m_z1(m_ax), m_z2(m_ax),
        m_ev_min(ev_min), m_ev_max(ev_max), m_degree(degree){}
 
    template<class ContainerType0, class ContainerType1>
    void symv( const ContainerType0& x, ContainerType1& y)
    {
        value_type theta = (m_ev_min+m_ev_max)/2., delta = (m_ev_max-m_ev_min)/2.;
        value_type c_k = 1./sqrt(m_ev_min*m_ev_max);
        dg::blas1::axpby( c_k/2., x, 0., y);
        if( m_degree == 0) return;
        dg::blas2::symv( m_op, x, m_ax);
        dg::blas1::axpby( 1./delta, m_ax, -theta/delta, x, m_z1); //T_{k-1} x
        c_k *= (sqrt( m_ev_min/m_ev_max) - 1.)/(sqrt(m_ev_min/m_ev_max)+1);
        dg::blas1::axpby( c_k, m_z1, 1., y);
        if( m_degree == 1) return;
        dg::blas1::copy( x, m_z2); //T_{k-2} x
        for( unsigned i=1; i<m_degree; i++)
        {
            dg::blas2::symv( m_op, m_z1, m_ax);
            dg::blas1::axpby( 1./delta, m_ax, -theta/delta, m_z1, m_ax); //Z T_{k-1}
            dg::blas1::axpby( 2., m_ax, -1., m_z2, m_z2); //T_k
            c_k *= (sqrt( m_ev_min/m_ev_max) - 1.)/(sqrt(m_ev_min/m_ev_max)+1);
            dg::blas1::axpby( c_k, m_z2, 1., y);
            using std::swap;
            swap(m_z1,m_z2);
        }
    }
    private:
    Matrix m_op;
    ContainerType m_ax, m_z1, m_z2;
    value_type m_ev_min, m_ev_max;
    unsigned m_degree;
};
 
template<class Matrix, class InnerPreconditioner, class ContainerType>
struct LeastSquaresPreconditioner
{
    using container_type = ContainerType;
    using value_type = get_value_type<ContainerType>; 
    LeastSquaresPreconditioner( Matrix op, InnerPreconditioner P, const ContainerType& copyable, value_type ev_max, unsigned degree):
        m_op(op), m_p(P), m_z(copyable),
        m_ev_max( ev_max), m_degree(degree){
            m_c = coeffs(degree);
    }
    const ContainerType& copyable()const{ return m_z;}
 
    template<class ContainerType0, class ContainerType1>
    void symv( const ContainerType0& x, ContainerType1& y)
    {
        //Horner scheme
        dg::blas1::axpby( m_c[m_degree],x, 0., m_z);
        for( int i=m_degree-1; i>=0; i--)
        {
            //dg::blas1::copy( m_z, y);
            dg::blas2::symv( m_p, m_z, y);
            dg::blas2::symv( m_op, y, m_z);
            dg::blas1::axpby( m_c[i], x, +4./m_ev_max, m_z);
        }
        //dg::blas1::copy( m_z, y);
        dg::blas2::symv( m_p, m_z, y);
    }
    private:
    std::vector<value_type> coeffs( unsigned degree){
        switch( degree){
            case 0: return {1.};
            case 1: return {5., -1.};
            case 2: return { 14., -7., 1.};
            case 3: return {30., -27., 9., -1.};
            case 4: return {55., -77., 44., -11., 1.};
            case 5: return {91., -182., 156., -65., 13., -1. };
            case 6: return {140., -378., 450., -275., 90., -15., 1. };
            case 7: return {204., -714.,1122., -935., 442., -119., 17., -1.};
            case 8: return {285.,-1254., 2508., -2717., 1729., -665., 152., -19., 1.};
            case 9: return {385., -2079., 5148.,-7007.,5733.,-2940.,952.,-189.,21.,-1.};
            default:
                if (degree > 10)
                    std::cerr << "WARNING: LeastSquares Case "<<degree<<" not implemented. Taking 10 instead!\n";
                return {506., -3289., 9867.,-16445.,16744.,-10948.,4692.,-1311.,230.,-23.,1. };
        };
    }
    std::vector<value_type> m_c;
    Matrix m_op;
    InnerPreconditioner m_p;
    ContainerType m_z;
    value_type m_ev_max;
    unsigned m_degree;
};
 
template<class M, class V>
struct TensorTraits<ChebyshevPreconditioner<M,V>>
{
    using value_type      = get_value_type<V>;
    using tensor_category = SelfMadeMatrixTag;
};
template<class M, class V>
struct TensorTraits<ModifiedChebyshevPreconditioner<M,V>>
{
    using value_type      = get_value_type<V>;
    using tensor_category = SelfMadeMatrixTag;
};
template<class M, class P, class V>
struct TensorTraits<LeastSquaresPreconditioner<M,P,V>>
{
    using value_type      = get_value_type<V>;
    using tensor_category = SelfMadeMatrixTag;
};
 
 
//template<class Matrix, class Container>
//struct WrapperSpectralShift
//{
//    WrapperSpectralShift( Matrix& op, value_type ev_max):
//        m_op(op), m_ev_max(ev_max){}
//    template<class ContainerType0, class ContainerType1>
//    void symv( const ContainerType0& x, ContainerType1& y)
//    {
//        dg::blas1::axpby( m_ev_max, x, 0., y);
//        dg::blas2::symv( -1., m_op, x, 1., y);
//    }
//
//    private:
//    Matrix& m_op;
//    value_type m_ev_max;
//
//};
//template<class M, class V>
//struct TensorTraits<detail::WrapperSpectralShift<M,V>>
//{
//    using value_type      = get_value_type<V>;
//    using tensor_category = SelfMadeMatrixTag;
//};
 
} //namespace dg
 
#endif // _DG_CHEB_