Operations.hpp

#pragma once
 
#include <amdis/linearalgebra/istl/MatrixBackend.hpp>
#include <amdis/linearalgebra/istl/VectorBackend.hpp>
#include <dune/istl/paamg/pinfo.hh>
 
namespace AMDiS {
 
  namespace Impl {
 
    template <class Vec, class Comm,
      std::enable_if_t<!std::is_same_v<Comm,Dune::Amg::SequentialInformation>, int> = 0>
    auto normISTL(Vec const& x, Comm const& comm)
    {
      return comm.norm(x);
    }
 
    template <class Vec>
    auto normISTL(Vec const& x, Dune::Amg::SequentialInformation const& comm)
    {
      return x.two_norm();
    }
 
  } // end namespace Impl
 
 
  // ||b - A*x||
  template <class T1, class C1, class T2, class T3>
  auto residuum(ISTLBCRSMatrix<T1,C1> const& A, ISTLBlockVector<T2> const& x, ISTLBlockVector<T3> const& b)
  {
    auto r = b.vector();
    A.matrix().mmv(x.vector(), r);
    return Impl::normISTL(r, A.comm().impl());
  }
 
  // ||b - A*x|| / ||b||
  template <class T1, class C1, class T2, class T3>
  auto relResiduum(ISTLBCRSMatrix<T1,C1> const& A, ISTLBlockVector<T2> const& x, ISTLBlockVector<T3> const& b)
  {
    return residuum(A,x,b) / Impl::normISTL(b.vector(), A.comm().impl());
  }
 
} // end namespace AMDiS