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