Operations.hpp

#pragma once

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

namespace AMDiS {

  // ||b - A*x||
  template <class Matrix, class RowDofMap, class ColDofMap>
  auto residuum(Matrix const& A, PetscVector<ColDofMap> const& x, PetscVector<RowDofMap> const& b)
  {
    Vec r;
    VecDuplicate(b.vector(), &r);
    MatMult(A.matrix(), x.vector(), r);
    VecAXPY(r, -1.0, b.vector());

    PetscReal res;
    VecNorm(r, NORM_2, &res);
    VecDestroy(&r);
    return res;
  }

  // ||b - A*x|| / ||b||
  template <class Matrix, class RowDofMap, class ColDofMap>
  auto relResiduum(Matrix const& A, PetscVector<ColDofMap> const& x, PetscVector<RowDofMap> const& b)
  {
    PetscReal bNrm;
    VecNorm(b.vector(), NORM_2, &bNrm);
    return residuum(A,x,b) / bNrm;
  }

} // end namespace AMDiS