Operations.hpp

#pragma once

#include <amdis/linearalgebra/istl/MatrixBackend.hpp>
#include <amdis/linearalgebra/istl/VectorBackend.hpp>

namespace AMDiS {

  // ||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 r.two_norm();
  }

  // ||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) / b.vector().two_norm();
  }

} // end namespace AMDiS