Operations.hpp

#pragma once

#include <boost/numeric/mtl/operation/two_norm.hpp>

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

namespace AMDiS {

  // ||b - A*x||
  template <class T1, class T2, class T3>
  auto residuum(MTLSparseMatrix<T1> const& A, MTLVector<T2> const& x, MTLVector<T3> const& b)
  {
    auto r = b.vector();
    r -= A.matrix() * x.vector();
    return two_norm(r);
  }

  // ||b - A*x|| / ||b||
  template <class T1, class T2, class T3>
  auto relResiduum(MTLSparseMatrix<T1> const& A, MTLVector<T2> const& x, MTLVector<T3> const& b)
  {
    return residuum(A,x,b) / two_norm(b.vector());
  }

} // end namespace AMDiS