lanczos.h

Go to the documentation of this file.

#pragma once
 
#include "dg/algorithm.h"
#include "tridiaginv.h"
#include "matrixfunction.h"
 
 
namespace dg{
namespace mat{
 
template< class ContainerType >
class UniversalLanczos
{
  public:
    using value_type = get_value_type<ContainerType>; 
    UniversalLanczos(){}
    UniversalLanczos( const ContainerType& copyable, unsigned max_iterations)
    {
        m_v = m_vp = m_vm = copyable;
        m_max_iter = max_iterations;
        m_iter = max_iterations;
        //sub matrix and vector
        set_iter( max_iterations);
    }
    template<class ...Params>
    void construct( Params&& ...ps)
    {
        //construct and swap
        *this = UniversalLanczos( std::forward<Params>( ps)...);
    }
 
    void set_max( unsigned new_max) {
        m_max_iter = new_max;
        set_iter( new_max);
    }
 
    unsigned get_max() const {return m_max_iter;}
 
    void set_verbose( bool verbose){ m_verbose = verbose;}
 
    double get_bnorm() const{return m_bnorm;}
 
    unsigned get_iter() const {return m_iter;}
 
    template < class MatrixType, class ContainerType0, class ContainerType1,
             class ContainerType2, class FuncTe1>
    unsigned solve(ContainerType0& x, FuncTe1 f,
            MatrixType&& A, const ContainerType1& b,
            const ContainerType2& weights, double eps,
            double nrmb_correction = 1.,
            std::string error_norm = "universal",
            double res_fac = 1.,
            unsigned q = 1 )
    {
        tridiag( f, std::forward<MatrixType>(A), b, weights, eps,
                nrmb_correction, error_norm, res_fac, q);
        if( "residual" == error_norm)
            m_yH = f( m_TH);
        //Compute x = |b|_M V f(T) e1
        normMbVy(std::forward<MatrixType>(A), m_TH, m_yH, x, b,
                m_bnorm);
        return m_iter;
    }
 
    template< class MatrixType, class ContainerType0, class ContainerType1>
    const dg::TriDiagonal<thrust::host_vector<double>>& tridiag( MatrixType&& A, const ContainerType0& b,
            const ContainerType1& weights, double eps = 1e-4,
            double nrmb_correction = 1.,
            std::string error_norm = "compute_extreme_EV",
            double res_fac = 1.,
            unsigned q = 1 )
    {
        auto op = make_Linear_Te1( -1);
        tridiag( op, std::forward<MatrixType>(A), b, weights, eps,
                nrmb_correction, error_norm, res_fac, q);
        return m_TH;
    }
 
    template< class MatrixType, class ContainerType0,
        class ContainerType1,class ContainerType2>
    void normMbVy( MatrixType&& A,
            const dg::TriDiagonal<thrust::host_vector<double>>& T,
            const ContainerType0& y,
            ContainerType1& x,
            const ContainerType2& b, double bnorm)
    {
        dg::blas1::copy(0., x);
        if( 0 == bnorm )
        {
            return;
        }
        dg::blas1::axpby(1./bnorm, b, 0.0, m_v); //m_v[1] = b/||b||
        dg::blas1::copy(0., m_vm);
        // check if (potentially) all higher elements in y are zero
        unsigned less_iter = 0;
        for( unsigned i=0; i<y.size(); i++)
            if( y[i] != 0)
                less_iter = i+1;
        dg::blas1::axpby( y[0]*bnorm, m_v, 1., x); //Compute b= |b| V y
 
        for ( unsigned i=0; i<less_iter-1; i++)
        {
            dg::blas2::symv( std::forward<MatrixType>(A), m_v, m_vp);
            dg::blas1::axpbypgz(
                    -T.M[i]/T.P[i], m_vm,
                    -T.O[i]/T.P[i], m_v,
                               1.0/T.P[i], m_vp);
            dg::blas1::axpby( y[i+1]*bnorm, m_vp, 1., x); //Compute b= |b| V y
            m_vm.swap( m_v);
            m_v.swap( m_vp);
        }
    }
 
    template < class UnaryOp, class MatrixType,
             class ContainerType1, class ContainerType2>
    const dg::TriDiagonal<thrust::host_vector<double>>& tridiag(UnaryOp f,
            MatrixType&& A, const ContainerType1& b,
            const ContainerType2& weights, double eps,
            double nrmb_correction,
            std::string error_norm = "residual",
            double res_fac = 1.,
            unsigned q = 1 )
    {
#ifdef MPI_VERSION
        int rank;
        MPI_Comm_rank(MPI_COMM_WORLD, &rank);
#endif //MPI
        m_bnorm = nrm( weights, b);
        if( m_verbose)
        {
            DG_RANK0 std::cout << "# Norm of b  "<<m_bnorm <<"\n";
            if ( "compute_extreme_EV" == error_norm || "compute_max_EV" == error_norm)
            {
                DG_RANK0 std::cout << "# Rel EV errors: \n";
            }
            else
            {
                DG_RANK0 std::cout << "# Res factor "<<res_fac <<"\n";
                DG_RANK0 std::cout << "# Residual errors: \n";
            }
        }
        if( m_bnorm == 0)
        {
            set_iter(1);
            return m_TH;
        }
        double residual = 1e300;
        dg::blas1::axpby(1./m_bnorm, b, 0.0, m_v); //m_v[1] = x/||x||
        double betaip = 0.;
        double alphai = 0.; // For Hermitian matrices all dot products in CG are real
        double lmin = -1e300, lmax = 1e300;
        for( unsigned i=0; i<m_max_iter; i++)
        {
            m_TH.M[i] =  betaip; // -1 diagonal
            dg::blas2::symv(std::forward<MatrixType>(A), m_v, m_vp);
            dg::blas1::axpby(-betaip, m_vm, 1.0, m_vp);  // only - if i>0, therefore no if (i>0)
            alphai  = dot(m_vp, weights, m_v);
            m_TH.O[i] = alphai;
            dg::blas1::axpby(-alphai, m_v, 1.0, m_vp);
            betaip = nrm( weights, m_vp);
            if (fabs( betaip) <= 1e-12) // what this means is that x = A^{-1}b is exactly in the Krylov subspace
            {
                set_iter(i+1);
                if ( "compute_extreme_EV" == error_norm || "compute_max_EV" == error_norm)
                {
                    // If we did continue iterating after betaip = 0 the results are rather random
                    throw dg::Error( dg::Message(_ping_) << "betaip = 0 found in tridiagonalisation \n");
                }
                else if( m_verbose)
                    DG_RANK0 std::cout << "beta["<<i+1 <<"]=0 encountered\n";
                break;
            }
            m_TH.P[i] = betaip;  // +1 diagonal
 
            double xnorm = 0.;
            if( "residual" == error_norm)
            {
                residual = compute_residual_error( m_TH, i)*m_bnorm;
                xnorm = m_bnorm;
            }
            else if ( "compute_extreme_EV" == error_norm || "compute_max_EV" == error_norm)
            {
                if(i%10 == 0)
                {
                    // We use P as "subdiagonal" because it is symmetric and the first element must be on 0 index
                    // (last is ignored)
                    thrust::host_vector<double> evals( i+1), subdiagonal( i+1), Z, work;
                    for( unsigned u=0; u<i+1; u++)
                    {
                        evals[u] = m_TH.O[u];
                        subdiagonal[u] = m_TH.P[u];
                    }
                    lapack::stev('N', evals,  subdiagonal, Z, work);
                    if( m_verbose)
                        DG_RANK0 std::cout << "# lmin = " << evals[0] << " lmax = "<<evals[i]<<"\tat i = " << i << "\n";// << betaip * evecs(i,i)<<" " <<betaip * evecs(0,i)<<"\n";
                    //if( fabs(betaip * evecs(i,i)) < 0.9*eps) // CONVERGES SUPER SLOW
                    if( fabs(lmax - evals[i]) < 10*eps*evals[i])
                    {
                        if( ("compute_extreme_EV" == error_norm &&
                            //fabs( betaip * evecs(0,i)) < eps) // CONVERGES EVEN SLOWER
                            fabs( lmin - evals[0]) < 10 * (eps*fabs(lmin) + 1e-12*lmax))
                            || "compute_max_EV" == error_norm)
                        {
                            set_iter(i+1);
                            break;
                        }
                    }
                    lmin = evals[0], lmax = evals[i];
                }
            }
            else
            {
                if( i>=q &&(  (i<=10) || (i>10 && i%10 == 0) ))
                {
                    residual = compute_universal_error( m_TH, i, q, f,
                            m_yH)*m_bnorm;
                    xnorm = dg::fast_l2norm( m_yH)*m_bnorm;
                }
                else
                {
                    residual = 1e10;
                    xnorm = m_bnorm;
                }
            }
            if( m_verbose && "compute_extreme_EV" != error_norm && "compute_max_EV" != error_norm)
                DG_RANK0 std::cout << "# ||r||_W = " << residual << "\tat i = " << i << "\n";
            if (res_fac*residual< eps*(xnorm + nrmb_correction) )
            {
                set_iter(i+1);
                break;
            }
            dg::blas1::scal(m_vp, 1./betaip);
            m_vm.swap(m_v);
            m_v.swap( m_vp);
            set_iter( m_max_iter);
        }
        return m_TH;
    }
    private:
    // !! Calling this "norm" would shadow std::norm in NORM
    template<class ContainerType1, class ContainerType2>
    auto nrm( const ContainerType1& w, const ContainerType2& x)
    {
        using value_type_x = dg::get_value_type<ContainerType2>;
        constexpr bool is_complex = dg::is_scalar_v<value_type_x, dg::ComplexTag>;
        if constexpr (is_complex)
            return sqrt( blas1::vdot( dg::detail::NORM(), w, x));
        else
            return sqrt( blas2::dot( w, x));
    }
    template<class ContainerType0, class ContainerType1, class ContainerType2>
    auto dot( const ContainerType0& x0, const ContainerType1& w, const ContainerType2& x1)
    {
        using value_type_x1 = dg::get_value_type<ContainerType2>;
        constexpr bool is_complex = dg::is_scalar_v<value_type_x1, dg::ComplexTag>;
        // There is no complex_symmetric in Lanczos...
        if constexpr (not is_complex ) // or complex_mode == dg::complex_symmetric)
            return blas2::dot( x0, w, x1); // this returns value_type
        else // For Hermitian matrices all dot products in CG are real
            return blas1::vdot( dg::detail::DOT(), x0, w, x1); // this returns floating point
    }
    double compute_residual_error( const dg::TriDiagonal<thrust::host_vector<double>>& TH, unsigned iter)
    {
        double T1 = compute_Tinv_m1( TH, iter+1);
        return TH.P[iter]*fabs(T1); //Tinv_i1
    }
    template<class UnaryOp>
    double compute_universal_error( const dg::TriDiagonal<thrust::host_vector<double>>& TH, unsigned iter,
            unsigned q, UnaryOp f, HVec& yH)
    {
        unsigned new_iter = iter + 1 + q;
        set_iter( iter+1);
        dg::TriDiagonal<thrust::host_vector<double>> THtilde( new_iter);
        for( unsigned u=0; u<iter+1; u++)
        {
            THtilde.M[u] = TH.M[u];
            THtilde.O[u] = TH.O[u];
            THtilde.P[u] = TH.P[u];
        }
        for( unsigned u=1; u<=q; u++)
        {
            THtilde.M[ iter+u] = u==1 ? TH.P[iter] :
                TH.O[iter+1-u];
            THtilde.O[ iter+u] = TH.O[ iter-u];
            THtilde.P[ iter+u] = TH.M[ iter-u];
        }
        yH = f( TH);
        HVec yHtilde = f( THtilde);
        for( unsigned u=0; u<yH.size(); u++)
            yHtilde[u] -= yH[u];
        double norm = dg::fast_l2norm( yHtilde);
        return norm;
    }
 
    void set_iter( unsigned new_iter) {
        m_TH.resize(new_iter);
        m_iter = new_iter;
    }
    ContainerType  m_v, m_vp, m_vm;
    dg::TriDiagonal<thrust::host_vector<double>> m_TH;
    HVec m_yH;
    unsigned m_iter, m_max_iter;
    bool m_verbose = false;
    double m_bnorm = 0.;
};
 
 
} //namespace mat
} //namespace dg