From 9292be0b69cfb8938d3be1bd8a4f2e4ce49cfed5 Mon Sep 17 00:00:00 2001
From: Daniel Richtmann <daniel.richtmann@physik.uni-regensburg.de>
Date: Wed, 31 Jan 2018 14:03:30 +0100
Subject: [PATCH] =?UTF-8?q?WilsonMG:=20Add=20check=20for=20Mdiag=20+=20?=
 =?UTF-8?q?=CE=A3=20Mdir=20=3D=3D=20M?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

Need to test my implementations of CoarsenedMatrix::Mdiag &
CoarsenedMatrix::Mdir.
---
 tests/solver/Test_wilson_mg.cc | 47 ++++++++++++++++++++++++++++++++++
 1 file changed, 47 insertions(+)

diff --git a/tests/solver/Test_wilson_mg.cc b/tests/solver/Test_wilson_mg.cc
index 25d18d9a..71fe5837 100644
--- a/tests/solver/Test_wilson_mg.cc
+++ b/tests/solver/Test_wilson_mg.cc
@@ -126,6 +126,40 @@ public:
   }
 };
 
+template<class Field> void testOperator(LinearOperatorBase<Field> &LinOp, GridBase *Grid) {
+
+  std::vector<int> seeds({1, 2, 3, 4});
+  GridParallelRNG  RNG(Grid);
+  RNG.SeedFixedIntegers(seeds);
+
+  // clang-format off
+  Field src(Grid);    random(RNG, src);
+  Field result(Grid); result = zero;
+  Field ref(Grid);    ref    = zero;
+  Field tmp(Grid);
+  Field err(Grid);
+  // clang-format on
+
+  LinOp.Op(src, ref);
+
+  LinOp.OpDiag(src, result);
+  std::cout << GridLogMessage << "diag:  norm2(result) = " << norm2(result) << std::endl;
+
+  for(int d = 0; d < 4; d++) {
+    LinOp.OpDir(src, tmp, d, +1);
+    std::cout << GridLogMessage << "dir + " << d << ": norm2(tmp) = " << norm2(tmp) << std::endl;
+    result = result + tmp;
+
+    LinOp.OpDir(src, tmp, d, -1);
+    std::cout << GridLogMessage << "dir - " << d << ": norm2(tmp) = " << norm2(tmp) << std::endl;
+    result = result + tmp;
+  }
+
+  err = result - ref;
+
+  std::cout << GridLogMessage << "Error: absolute = " << norm2(err) << " relative = " << norm2(err) / norm2(ref) << std::endl;
+}
+
 // template < class Fobj, class CComplex, int coarseSpins, int nbasis, class Matrix >
 // class MultiGridPreconditioner : public LinearFunction< Lattice< Fobj > > {
 template<class Fobj, class CComplex, int nbasis, class Matrix> class MultiGridPreconditioner : public LinearFunction<Lattice<Fobj>> {
@@ -815,6 +849,19 @@ int main(int argc, char **argv) {
   coarseResult = zero;
   CoarseCG(CoarsePosDefHermOp, coarseSource, coarseResult);
 
+  std::cout << GridLogMessage << "**************************************************" << std::endl;
+  std::cout << GridLogMessage << "Testing the operators" << std::endl;
+  std::cout << GridLogMessage << "**************************************************" << std::endl;
+
+  std::cout << GridLogMessage << "MdagMLinearOperator<WilsonFermionR, LatticeFermion> FineHermOp(Dw);" << std::endl;
+  testOperator(FineHermOp, FGrid);
+  std::cout << GridLogMessage << "Gamma5HermitianLinearOperator<WilsonFermionR, LatticeFermion> FineHermIndefOp(Dw);" << std::endl;
+  testOperator(FineHermIndefOp, FGrid);
+  std::cout << GridLogMessage << "Gamma5HermitianLinearOperator<WilsonFermionR, LatticeFermion> FineHermIndefOpDD(DwDD);" << std::endl;
+  testOperator(FineHermIndefOpDD, FGrid);
+  std::cout << GridLogMessage << "MdagMLinearOperator<CoarseOperator, CoarseVector> CoarsePosDefHermOp(Dc);" << std::endl;
+  testOperator(CoarsePosDefHermOp, coarseGrids.Grids[0]);
+
   std::cout << GridLogMessage << "**************************************************" << std::endl;
   std::cout << GridLogMessage << "Building deflation preconditioner " << std::endl;
   std::cout << GridLogMessage << "**************************************************" << std::endl;