AMReX_KrylovMV.H

Go to the documentation of this file.

#ifndef AMREX_KRYLOV_MV_H_
#define AMREX_KRYLOV_MV_H_
 
#include <AMReX_Algebra.H>
#include <AMReX_Smoother_MV.H>
#include <functional>
#include <type_traits>
 
namespace amrex {
 
template <typename V, typename M> class GMRES;
 
template <template <typename, typename> class S, typename T>
class KrylovMV
{
public:
    using RT = T;
    using VEC = AlgVector<T>;
    using MAT = SpMatrix<T>;
    using Solver = S<VEC,KrylovMV<S,T>>;
    using GM = Solver;
    using PC = std::function<void(VEC&,VEC const&)>;
 
    KrylovMV (MAT const* a_mat);
 
    ~KrylovMV () = default;
 
    // The solver holds a pointer to this object.
    KrylovMV (const KrylovMV&) = delete;
    KrylovMV (KrylovMV&&) = delete;
    KrylovMV& operator= (const KrylovMV&) = delete;
    KrylovMV& operator= (KrylovMV&&) = delete;
 
    void setPrecond (PC a_pc) { m_pc = std::move(a_pc); }
 
    void solve (VEC& a_sol, VEC const& a_rhs, T a_tol_rel, T a_tol_abs);
 
    void setVerbose (int v) { m_solver.setVerbose(v); }
 
    Solver& getSolver () { return m_solver; }
 
    Solver& getGMRES () requires std::is_same_v<Solver, GMRES<VEC,KrylovMV<S,T>>>
    { return m_solver; }
 
    [[nodiscard]] VEC makeVecRHS () const;
 
    [[nodiscard]] VEC makeVecLHS () const;
 
    static T norm2 (VEC const& vec);
 
    static void scale (VEC& vec, T scale_factor);
 
    static T dotProduct (VEC const& vec1, VEC const& vec2);
 
    static void setToZero (VEC& lhs);
 
    static void assign (VEC& lhs, VEC const& rhs);
 
    static void increment (VEC& lhs, VEC const& rhs, T a);
 
    static void linComb (VEC& lhs, T a, VEC const& rhs_a, T b, VEC const& rhs_b);
 
    void apply (VEC& lhs, VEC& rhs) const;
 
    void precond (VEC& lhs, VEC const& rhs) const;
 
private:
    Solver m_solver;
    MAT const* m_mat = nullptr;
    PC m_pc;
};
 
template <template <typename, typename> class S, typename T>
KrylovMV<S,T>::KrylovMV (MAT const* a_mat)
    : m_mat(a_mat)
{
    m_solver.define(*this);
}
 
template <template <typename, typename> class S, typename T>
void KrylovMV<S,T>::solve (VEC& a_sol, VEC const& a_rhs, T a_tol_rel, T a_tol_abs)
{
    m_solver.solve(a_sol, a_rhs, a_tol_rel, a_tol_abs);
}
 
template <template <typename, typename> class S, typename T>
auto KrylovMV<S,T>::makeVecRHS () const -> VEC
{
    return VEC(m_mat->partition());
}
 
template <template <typename, typename> class S, typename T>
auto KrylovMV<S,T>::makeVecLHS () const -> VEC
{
    return VEC(m_mat->partition());
}
 
template <template <typename, typename> class S, typename T>
T KrylovMV<S,T>::norm2 (VEC const& vec)
{
    return vec.norm2();
}
 
template <template <typename, typename> class S, typename T>
void KrylovMV<S,T>::scale (VEC& vec, T scale_factor)
{
    vec.scaleAsync(scale_factor);
}
 
template <template <typename, typename> class S, typename T>
T KrylovMV<S,T>::dotProduct (VEC const& vec1, VEC const& vec2)
{
    return amrex::Dot(vec1,vec2);
}
 
template <template <typename, typename> class S, typename T>
void KrylovMV<S,T>::setToZero (VEC& lhs)
{
    lhs.setValAsync(0);
}
 
template <template <typename, typename> class S, typename T>
void KrylovMV<S,T>::assign (VEC& lhs, VEC const& rhs)
{
    lhs.copyAsync(rhs);
}
 
template <template <typename, typename> class S, typename T>
void KrylovMV<S,T>::increment (VEC& lhs, VEC const& rhs, T a)
{
    amrex::Axpy(lhs, a, rhs);
}
 
template <template <typename, typename> class S, typename T>
void KrylovMV<S,T>::linComb (VEC& lhs, T a, VEC const& rhs_a, T b, VEC const& rhs_b)
{
    amrex::LinComb(lhs, a, rhs_a, b, rhs_b);
}
 
template <template <typename, typename> class S, typename T>
void KrylovMV<S,T>::apply (VEC& lhs, VEC& rhs) const
{
    amrex::SpMV(lhs, *m_mat, rhs);
}
 
template <template <typename, typename> class S, typename T>
void KrylovMV<S,T>::precond (VEC& lhs, VEC const& rhs) const
{
    if (m_pc) {
        m_pc(lhs, rhs);
    } else {
        lhs.copyAsync(rhs);
    }
}
 
}
#endif