Operations.hpp

#pragma once

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

namespace AMDiS {

  // ||b - A*x||
  template <class T1, int O1, class T2, class T3>
  auto residuum(EigenSparseMatrix<T1,O1> const& A, EigenVector<T2> const& x, EigenVector<T3> const& b)
  {
    auto r = b.vector();
    r.noalias() -= A.matrix() * x.vector();
    return r.norm();
  }

  // ||b - A*x|| / ||b||
  template <class T1, int O1, class T2, class T3>
  auto relResiduum(EigenSparseMatrix<T1,O1> const& A, EigenVector<T2> const& x, EigenVector<T3> const& b)
  {
    return residuum(A,x,b) / b.vector().norm();
  }

} // end namespace AMDiS