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