mirror of
https://github.com/paboyle/Grid.git
synced 2025-11-04 05:54:32 +00:00
Enable even-odd for CoarsenedMatrix
This commit is contained in:
Daniel Richtmann
475
tests/solver/Test_coarse_even_odd.cc
Normal file
475
tests/solver/Test_coarse_even_odd.cc
Normal file
@@ -0,0 +1,475 @@
|
||||
/*************************************************************************************
|
||||
|
||||
Grid physics library, www.github.com/paboyle/Grid
|
||||
|
||||
Source file: ./tests/solver/Test_coarse_even_odd.cc
|
||||
|
||||
Copyright (C) 2015-2020
|
||||
|
||||
Author: Daniel Richtmann <daniel.richtmann@ur.de>
|
||||
|
||||
This program is free software; you can redistribute it and/or modify
|
||||
it under the terms of the GNU General Public License as published by
|
||||
the Free Software Foundation; either version 2 of the License, or
|
||||
(at your option) any later version.
|
||||
|
||||
This program is distributed in the hope that it will be useful,
|
||||
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||
GNU General Public License for more details.
|
||||
|
||||
You should have received a copy of the GNU General Public License along
|
||||
with this program; if not, write to the Free Software Foundation, Inc.,
|
||||
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
|
||||
|
||||
See the full license in the file "LICENSE" in the top level distribution directory
|
||||
*************************************************************************************/
|
||||
/* END LEGAL */
|
||||
|
||||
#include <Grid/Grid.h>
|
||||
|
||||
using namespace Grid;
|
||||
|
||||
#ifndef NBASIS
|
||||
#define NBASIS 40
|
||||
#endif
|
||||
|
||||
// NOTE: The tests in this file are written in analogy to
|
||||
// - tests/core/Test_wilson_even_odd.cc
|
||||
// - tests/core/Test_wilson_clover.cc
|
||||
|
||||
std::vector<int> readFromCommandlineIvec(int* argc,
|
||||
char*** argv,
|
||||
std::string&& option,
|
||||
const std::vector<int>& defaultValue) {
|
||||
std::string arg;
|
||||
std::vector<int> ret(defaultValue);
|
||||
if(GridCmdOptionExists(*argv, *argv + *argc, option)) {
|
||||
arg = GridCmdOptionPayload(*argv, *argv + *argc, option);
|
||||
GridCmdOptionIntVector(arg, ret);
|
||||
}
|
||||
return ret;
|
||||
}
|
||||
|
||||
int main(int argc, char** argv) {
|
||||
Grid_init(&argc, &argv);
|
||||
|
||||
/////////////////////////////////////////////////////////////////////////////
|
||||
// Read from command line //
|
||||
/////////////////////////////////////////////////////////////////////////////
|
||||
|
||||
const int nbasis = NBASIS; static_assert((nbasis & 0x1) == 0, "");
|
||||
const int nb = nbasis/2;
|
||||
Coordinate blockSize = readFromCommandlineIvec(&argc, &argv, "--blocksize", {2, 2, 2, 2});
|
||||
|
||||
std::cout << GridLogMessage << "Compiled with nbasis = " << nbasis << " -> nb = " << nb << std::endl;
|
||||
|
||||
/////////////////////////////////////////////////////////////////////////////
|
||||
// General setup //
|
||||
/////////////////////////////////////////////////////////////////////////////
|
||||
|
||||
Coordinate clatt = GridDefaultLatt();
|
||||
for(int d=0; d<clatt.size(); d++) clatt[d] = clatt[d] / blockSize[d];
|
||||
|
||||
GridCartesian* Grid_f = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd, vComplex::Nsimd()), GridDefaultMpi());
|
||||
GridCartesian* Grid_c = SpaceTimeGrid::makeFourDimGrid(clatt, GridDefaultSimd(Nd, vComplex::Nsimd()), GridDefaultMpi());
|
||||
GridRedBlackCartesian* RBGrid_f = SpaceTimeGrid::makeFourDimRedBlackGrid(Grid_f);
|
||||
GridRedBlackCartesian* RBGrid_c = SpaceTimeGrid::makeFourDimRedBlackGrid(Grid_c);
|
||||
|
||||
std::cout << GridLogMessage << "Grid_f:" << std::endl; Grid_f->show_decomposition();
|
||||
std::cout << GridLogMessage << "Grid_c:" << std::endl; Grid_c->show_decomposition();
|
||||
std::cout << GridLogMessage << "RBGrid_f:" << std::endl; RBGrid_f->show_decomposition();
|
||||
std::cout << GridLogMessage << "RBGrid_c:" << std::endl; RBGrid_c->show_decomposition();
|
||||
|
||||
GridParallelRNG pRNG_f(Grid_f);
|
||||
GridParallelRNG pRNG_c(Grid_c);
|
||||
|
||||
std::vector<int> seeds({1, 2, 3, 4});
|
||||
|
||||
pRNG_f.SeedFixedIntegers(seeds);
|
||||
pRNG_c.SeedFixedIntegers(seeds);
|
||||
|
||||
/////////////////////////////////////////////////////////////////////////////
|
||||
// Setup of Dirac Matrix and Operator //
|
||||
/////////////////////////////////////////////////////////////////////////////
|
||||
|
||||
LatticeGaugeField Umu(Grid_f); SU3::HotConfiguration(pRNG_f, Umu);
|
||||
|
||||
RealD checkTolerance = (getPrecision<LatticeFermion>::value == 1) ? 1e-7 : 1e-15;
|
||||
|
||||
RealD mass = -0.30;
|
||||
RealD csw = 1.9192;
|
||||
|
||||
WilsonCloverFermionR Dwc(Umu, *Grid_f, *RBGrid_f, mass, csw, csw);
|
||||
MdagMLinearOperator<WilsonCloverFermionR, LatticeFermion> MdagMOp_Dwc(Dwc);
|
||||
|
||||
/////////////////////////////////////////////////////////////////////////////
|
||||
// Type definitions //
|
||||
/////////////////////////////////////////////////////////////////////////////
|
||||
|
||||
typedef Aggregation<vSpinColourVector, vTComplex, nbasis> Aggregates;
|
||||
typedef CoarsenedMatrix<vSpinColourVector, vTComplex, nbasis> CoarseDiracMatrix;
|
||||
typedef CoarseDiracMatrix::CoarseVector CoarseVector;
|
||||
|
||||
/////////////////////////////////////////////////////////////////////////////
|
||||
// Setup of Aggregation //
|
||||
/////////////////////////////////////////////////////////////////////////////
|
||||
|
||||
Aggregates Aggs(Grid_c, Grid_f, 0);
|
||||
{
|
||||
LatticeFermion tmp(Aggs.subspace[0].Grid());
|
||||
for(int n = 0; n < nb; n++) {
|
||||
gaussian(pRNG_f, Aggs.subspace[n]);
|
||||
G5C(tmp, Aggs.subspace[n]);
|
||||
axpby(Aggs.subspace[n + nb], 0.5, -0.5, Aggs.subspace[n], tmp);
|
||||
axpby(Aggs.subspace[n], 0.5, 0.5, Aggs.subspace[n], tmp);
|
||||
}
|
||||
}
|
||||
|
||||
/////////////////////////////////////////////////////////////////////////////
|
||||
// Setup of CoarsenedMatrix and Operator //
|
||||
/////////////////////////////////////////////////////////////////////////////
|
||||
|
||||
const int hermitian = 0;
|
||||
CoarseDiracMatrix Dc(*Grid_c, *RBGrid_c, hermitian);
|
||||
Dc.CoarsenOperator(Grid_f, MdagMOp_Dwc, Aggs);
|
||||
MdagMLinearOperator<CoarseDiracMatrix, CoarseVector> MdagMOp_Dc(Dc);
|
||||
|
||||
/////////////////////////////////////////////////////////////////////////////
|
||||
// Setup vectors used in all tests //
|
||||
/////////////////////////////////////////////////////////////////////////////
|
||||
|
||||
CoarseVector src(Grid_c); random(pRNG_c, src);
|
||||
CoarseVector diff(Grid_c); diff = Zero();
|
||||
|
||||
/////////////////////////////////////////////////////////////////////////////
|
||||
// Start of tests //
|
||||
/////////////////////////////////////////////////////////////////////////////
|
||||
|
||||
{
|
||||
std::cout << GridLogMessage << "===========================================================================" << std::endl;
|
||||
std::cout << GridLogMessage << "= Test Dhop + Mdiag = Munprec" << std::endl;
|
||||
std::cout << GridLogMessage << "===========================================================================" << std::endl;
|
||||
|
||||
CoarseVector phi(Grid_c); phi = Zero();
|
||||
CoarseVector chi(Grid_c); chi = Zero();
|
||||
CoarseVector res(Grid_c); res = Zero();
|
||||
CoarseVector ref(Grid_c); ref = Zero();
|
||||
|
||||
Dc.Mdiag(src, phi); std::cout << GridLogMessage << "Applied Mdiag" << std::endl;
|
||||
Dc.Dhop(src, chi, DaggerNo); std::cout << GridLogMessage << "Applied Dhop" << std::endl;
|
||||
Dc.M(src, ref); std::cout << GridLogMessage << "Applied M" << std::endl;
|
||||
|
||||
res = phi + chi;
|
||||
|
||||
diff = ref - res;
|
||||
auto absDev = norm2(diff);
|
||||
auto relDev = absDev / norm2(ref);
|
||||
std::cout << GridLogMessage << "norm2(Munprec), norm2(Dhop + Mdiag), abs. deviation, rel. deviation: "
|
||||
<< norm2(ref) << " " << norm2(res) << " " << absDev << " " << relDev << " -> check "
|
||||
<< ((relDev < checkTolerance) ? "passed" : "failed") << std::endl;
|
||||
assert(relDev <= checkTolerance);
|
||||
}
|
||||
|
||||
{
|
||||
std::cout << GridLogMessage << "===========================================================================" << std::endl;
|
||||
std::cout << GridLogMessage << "= Test Meo + Moe = Dhop" << std::endl;
|
||||
std::cout << GridLogMessage << "===========================================================================" << std::endl;
|
||||
|
||||
CoarseVector src_e(RBGrid_c); src_e = Zero();
|
||||
CoarseVector src_o(RBGrid_c); src_o = Zero();
|
||||
CoarseVector res_e(RBGrid_c); res_e = Zero();
|
||||
CoarseVector res_o(RBGrid_c); res_o = Zero();
|
||||
CoarseVector res(Grid_c); res = Zero();
|
||||
CoarseVector ref(Grid_c); ref = Zero();
|
||||
|
||||
pickCheckerboard(Even, src_e, src);
|
||||
pickCheckerboard(Odd, src_o, src);
|
||||
|
||||
Dc.Meooe(src_e, res_o); std::cout << GridLogMessage << "Applied Meo" << std::endl;
|
||||
Dc.Meooe(src_o, res_e); std::cout << GridLogMessage << "Applied Moe" << std::endl;
|
||||
Dc.Dhop(src, ref, DaggerNo); std::cout << GridLogMessage << "Applied Dhop" << std::endl;
|
||||
|
||||
setCheckerboard(res, res_o);
|
||||
setCheckerboard(res, res_e);
|
||||
|
||||
diff = ref - res;
|
||||
auto absDev = norm2(diff);
|
||||
auto relDev = absDev / norm2(ref);
|
||||
std::cout << GridLogMessage << "norm2(Dhop), norm2(Meo + Moe), abs. deviation, rel. deviation: "
|
||||
<< norm2(ref) << " " << norm2(res) << " " << absDev << " " << relDev
|
||||
<< " -> check " << ((relDev < checkTolerance) ? "passed" : "failed") << std::endl;
|
||||
assert(relDev <= checkTolerance);
|
||||
}
|
||||
|
||||
{
|
||||
std::cout << GridLogMessage << "===========================================================================" << std::endl;
|
||||
std::cout << GridLogMessage << "= Test |(Im(v^dag M^dag M v)| = 0" << std::endl;
|
||||
std::cout << GridLogMessage << "===========================================================================" << std::endl;
|
||||
|
||||
CoarseVector tmp(Grid_c); tmp = Zero();
|
||||
CoarseVector phi(Grid_c); phi = Zero();
|
||||
|
||||
Dc.M(src, tmp); std::cout << GridLogMessage << "Applied M" << std::endl;
|
||||
Dc.Mdag(tmp, phi); std::cout << GridLogMessage << "Applied Mdag" << std::endl;
|
||||
|
||||
std::cout << GridLogMessage << "src = " << norm2(src) << " tmp = " << norm2(tmp) << " phi = " << norm2(phi) << std::endl;
|
||||
|
||||
ComplexD dot = innerProduct(src, phi);
|
||||
|
||||
auto relDev = abs(imag(dot)) / abs(real(dot));
|
||||
std::cout << GridLogMessage << "Re(v^dag M^dag M v), Im(v^dag M^dag M v), rel.deviation: "
|
||||
<< real(dot) << " " << imag(dot) << " " << relDev
|
||||
<< " -> check " << ((relDev < checkTolerance) ? "passed" : "failed") << std::endl;
|
||||
assert(relDev <= checkTolerance);
|
||||
}
|
||||
|
||||
{
|
||||
std::cout << GridLogMessage << "===========================================================================" << std::endl;
|
||||
std::cout << GridLogMessage << "= Test |(Im(v^dag Mooee^dag Mooee v)| = 0 (full grid)" << std::endl;
|
||||
std::cout << GridLogMessage << "===========================================================================" << std::endl;
|
||||
|
||||
CoarseVector tmp(Grid_c); tmp = Zero();
|
||||
CoarseVector phi(Grid_c); phi = Zero();
|
||||
|
||||
Dc.Mooee(src, tmp); std::cout << GridLogMessage << "Applied Mooee" << std::endl;
|
||||
Dc.MooeeDag(tmp, phi); std::cout << GridLogMessage << "Applied MooeeDag" << std::endl;
|
||||
|
||||
ComplexD dot = innerProduct(src, phi);
|
||||
|
||||
auto relDev = abs(imag(dot)) / abs(real(dot));
|
||||
std::cout << GridLogMessage << "Re(v^dag Mooee^dag Mooee v), Im(v^dag Mooee^dag Mooee v), rel.deviation: "
|
||||
<< real(dot) << " " << imag(dot) << " " << relDev
|
||||
<< " -> check " << ((relDev < checkTolerance) ? "passed" : "failed") << std::endl;
|
||||
assert(relDev <= checkTolerance);
|
||||
}
|
||||
|
||||
{
|
||||
std::cout << GridLogMessage << "===========================================================================" << std::endl;
|
||||
std::cout << GridLogMessage << "= Test MooeeInv Mooee = 1 (full grid)" << std::endl;
|
||||
std::cout << GridLogMessage << "===========================================================================" << std::endl;
|
||||
|
||||
CoarseVector tmp(Grid_c); tmp = Zero();
|
||||
CoarseVector phi(Grid_c); phi = Zero();
|
||||
|
||||
Dc.Mooee(src, tmp); std::cout << GridLogMessage << "Applied Mooee" << std::endl;
|
||||
Dc.MooeeInv(tmp, phi); std::cout << GridLogMessage << "Applied MooeeInv" << std::endl;
|
||||
|
||||
diff = src - phi;
|
||||
auto absDev = norm2(diff);
|
||||
auto relDev = absDev / norm2(src);
|
||||
std::cout << GridLogMessage << "norm2(src), norm2(MooeeInv Mooee src), abs. deviation, rel. deviation: "
|
||||
<< norm2(src) << " " << norm2(phi) << " " << absDev << " " << relDev
|
||||
<< " -> check " << ((relDev < checkTolerance) ? "passed" : "failed") << std::endl;
|
||||
assert(relDev <= checkTolerance);
|
||||
}
|
||||
|
||||
{
|
||||
std::cout << GridLogMessage << "===========================================================================" << std::endl;
|
||||
std::cout << GridLogMessage << "= Test MeooeDagger is the dagger of Meooe by requiring" << std::endl;
|
||||
std::cout << GridLogMessage << "= < phi | Meooe | chi > * = < chi | Meooe^dag| phi>" << std::endl;
|
||||
std::cout << GridLogMessage << "===========================================================================" << std::endl;
|
||||
|
||||
// clang-format off
|
||||
CoarseVector phi(Grid_c); random(pRNG_c, phi);
|
||||
CoarseVector chi(Grid_c); random(pRNG_c, chi);
|
||||
CoarseVector chi_e(RBGrid_c); chi_e = Zero();
|
||||
CoarseVector chi_o(RBGrid_c); chi_o = Zero();
|
||||
CoarseVector dchi_e(RBGrid_c); dchi_e = Zero();
|
||||
CoarseVector dchi_o(RBGrid_c); dchi_o = Zero();
|
||||
CoarseVector phi_e(RBGrid_c); phi_e = Zero();
|
||||
CoarseVector phi_o(RBGrid_c); phi_o = Zero();
|
||||
CoarseVector dphi_e(RBGrid_c); dphi_e = Zero();
|
||||
CoarseVector dphi_o(RBGrid_c); dphi_o = Zero();
|
||||
// clang-format on
|
||||
|
||||
pickCheckerboard(Even, chi_e, chi);
|
||||
pickCheckerboard(Odd, chi_o, chi);
|
||||
pickCheckerboard(Even, phi_e, phi);
|
||||
pickCheckerboard(Odd, phi_o, phi);
|
||||
|
||||
Dc.Meooe(chi_e, dchi_o); std::cout << GridLogMessage << "Applied Meo" << std::endl;
|
||||
Dc.Meooe(chi_o, dchi_e); std::cout << GridLogMessage << "Applied Moe" << std::endl;
|
||||
Dc.MeooeDag(phi_e, dphi_o); std::cout << GridLogMessage << "Applied MeoDag" << std::endl;
|
||||
Dc.MeooeDag(phi_o, dphi_e); std::cout << GridLogMessage << "Applied MoeDag" << std::endl;
|
||||
|
||||
ComplexD phiDchi_e = innerProduct(phi_e, dchi_e);
|
||||
ComplexD phiDchi_o = innerProduct(phi_o, dchi_o);
|
||||
ComplexD chiDphi_e = innerProduct(chi_e, dphi_e);
|
||||
ComplexD chiDphi_o = innerProduct(chi_o, dphi_o);
|
||||
|
||||
std::cout << GridLogDebug << "norm dchi_e = " << norm2(dchi_e) << " norm dchi_o = " << norm2(dchi_o) << " norm dphi_e = " << norm2(dphi_e)
|
||||
<< " norm dphi_o = " << norm2(dphi_e) << std::endl;
|
||||
|
||||
std::cout << GridLogMessage << "e " << phiDchi_e << " " << chiDphi_e << std::endl;
|
||||
std::cout << GridLogMessage << "o " << phiDchi_o << " " << chiDphi_o << std::endl;
|
||||
|
||||
std::cout << GridLogMessage << "phiDchi_e - conj(chiDphi_o) " << phiDchi_e - conj(chiDphi_o) << std::endl;
|
||||
std::cout << GridLogMessage << "phiDchi_o - conj(chiDphi_e) " << phiDchi_o - conj(chiDphi_e) << std::endl;
|
||||
}
|
||||
|
||||
{
|
||||
std::cout << GridLogMessage << "===========================================================================" << std::endl;
|
||||
std::cout << GridLogMessage << "= Test MooeeInv Mooee = 1 (checkerboards separately)" << std::endl;
|
||||
std::cout << GridLogMessage << "===========================================================================" << std::endl;
|
||||
|
||||
CoarseVector chi(Grid_c); random(pRNG_c, chi);
|
||||
CoarseVector tmp(Grid_c); tmp = Zero();
|
||||
CoarseVector phi(Grid_c); phi = Zero();
|
||||
CoarseVector chi_e(RBGrid_c); chi_e = Zero();
|
||||
CoarseVector chi_o(RBGrid_c); chi_o = Zero();
|
||||
CoarseVector phi_e(RBGrid_c); phi_e = Zero();
|
||||
CoarseVector phi_o(RBGrid_c); phi_o = Zero();
|
||||
CoarseVector tmp_e(RBGrid_c); tmp_e = Zero();
|
||||
CoarseVector tmp_o(RBGrid_c); tmp_o = Zero();
|
||||
|
||||
pickCheckerboard(Even, chi_e, chi);
|
||||
pickCheckerboard(Odd, chi_o, chi);
|
||||
pickCheckerboard(Even, tmp_e, tmp);
|
||||
pickCheckerboard(Odd, tmp_o, tmp);
|
||||
|
||||
Dc.Mooee(chi_e, tmp_e); std::cout << GridLogMessage << "Applied Mee" << std::endl;
|
||||
Dc.MooeeInv(tmp_e, phi_e); std::cout << GridLogMessage << "Applied MeeInv" << std::endl;
|
||||
Dc.Mooee(chi_o, tmp_o); std::cout << GridLogMessage << "Applied Moo" << std::endl;
|
||||
Dc.MooeeInv(tmp_o, phi_o); std::cout << GridLogMessage << "Applied MooInv" << std::endl;
|
||||
|
||||
setCheckerboard(phi, phi_e);
|
||||
setCheckerboard(phi, phi_o);
|
||||
|
||||
diff = chi - phi;
|
||||
auto absDev = norm2(diff);
|
||||
auto relDev = absDev / norm2(chi);
|
||||
std::cout << GridLogMessage << "norm2(chi), norm2(MeeInv Mee chi), abs. deviation, rel. deviation: "
|
||||
<< norm2(chi) << " " << norm2(phi) << " " << absDev << " " << relDev
|
||||
<< " -> check " << ((relDev < checkTolerance) ? "passed" : "failed") << std::endl;
|
||||
assert(relDev <= checkTolerance);
|
||||
}
|
||||
|
||||
{
|
||||
std::cout << GridLogMessage << "===========================================================================" << std::endl;
|
||||
std::cout << GridLogMessage << "= Test MooeeDag MooeeInvDag = 1 (checkerboards separately)" << std::endl;
|
||||
std::cout << GridLogMessage << "===========================================================================" << std::endl;
|
||||
|
||||
CoarseVector chi(Grid_c); random(pRNG_c, chi);
|
||||
CoarseVector tmp(Grid_c); tmp = Zero();
|
||||
CoarseVector phi(Grid_c); phi = Zero();
|
||||
CoarseVector chi_e(RBGrid_c); chi_e = Zero();
|
||||
CoarseVector chi_o(RBGrid_c); chi_o = Zero();
|
||||
CoarseVector phi_e(RBGrid_c); phi_e = Zero();
|
||||
CoarseVector phi_o(RBGrid_c); phi_o = Zero();
|
||||
CoarseVector tmp_e(RBGrid_c); tmp_e = Zero();
|
||||
CoarseVector tmp_o(RBGrid_c); tmp_o = Zero();
|
||||
|
||||
pickCheckerboard(Even, chi_e, chi);
|
||||
pickCheckerboard(Odd, chi_o, chi);
|
||||
pickCheckerboard(Even, tmp_e, tmp);
|
||||
pickCheckerboard(Odd, tmp_o, tmp);
|
||||
|
||||
Dc.MooeeDag(chi_e, tmp_e); std::cout << GridLogMessage << "Applied MeeDag" << std::endl;
|
||||
Dc.MooeeInvDag(tmp_e, phi_e); std::cout << GridLogMessage << "Applied MeeInvDag" << std::endl;
|
||||
Dc.MooeeDag(chi_o, tmp_o); std::cout << GridLogMessage << "Applied MooDag" << std::endl;
|
||||
Dc.MooeeInvDag(tmp_o, phi_o); std::cout << GridLogMessage << "Applied MooInvDag" << std::endl;
|
||||
|
||||
setCheckerboard(phi, phi_e);
|
||||
setCheckerboard(phi, phi_o);
|
||||
|
||||
diff = chi - phi;
|
||||
auto absDev = norm2(diff);
|
||||
auto relDev = absDev / norm2(chi);
|
||||
std::cout << GridLogMessage << "norm2(chi), norm2(MeeDag MeeInvDag chi), abs. deviation, rel. deviation: "
|
||||
<< norm2(chi) << " " << norm2(phi) << " " << absDev << " " << relDev
|
||||
<< " -> check " << ((relDev < checkTolerance) ? "passed" : "failed") << std::endl;
|
||||
assert(relDev <= checkTolerance);
|
||||
}
|
||||
|
||||
{
|
||||
std::cout << GridLogMessage << "===========================================================================" << std::endl;
|
||||
std::cout << GridLogMessage << "= Test Meo + Moe + Moo + Mee = Munprec" << std::endl;
|
||||
std::cout << GridLogMessage << "===========================================================================" << std::endl;
|
||||
|
||||
CoarseVector chi(Grid_c); chi = Zero();
|
||||
CoarseVector phi(Grid_c); phi = Zero();
|
||||
CoarseVector ref(Grid_c); ref = Zero();
|
||||
CoarseVector src_e(RBGrid_c); src_e = Zero();
|
||||
CoarseVector src_o(RBGrid_c); src_o = Zero();
|
||||
CoarseVector phi_e(RBGrid_c); phi_e = Zero();
|
||||
CoarseVector phi_o(RBGrid_c); phi_o = Zero();
|
||||
CoarseVector chi_e(RBGrid_c); chi_e = Zero();
|
||||
CoarseVector chi_o(RBGrid_c); chi_o = Zero();
|
||||
|
||||
pickCheckerboard(Even, src_e, src);
|
||||
pickCheckerboard(Odd, src_o, src);
|
||||
pickCheckerboard(Even, phi_e, phi);
|
||||
pickCheckerboard(Odd, phi_o, phi);
|
||||
pickCheckerboard(Even, chi_e, chi);
|
||||
pickCheckerboard(Odd, chi_o, chi);
|
||||
|
||||
// M phi = (Mooee src_e + Meooe src_o , Mooee src_o + Meooe src_e)
|
||||
|
||||
Dc.M(src, ref); // Reference result from the unpreconditioned operator
|
||||
|
||||
// EO matrix
|
||||
Dc.Mooee(src_e, chi_e); std::cout << GridLogMessage << "Applied Mee" << std::endl;
|
||||
Dc.Mooee(src_o, chi_o); std::cout << GridLogMessage << "Applied Moo" << std::endl;
|
||||
Dc.Meooe(src_o, phi_e); std::cout << GridLogMessage << "Applied Moe" << std::endl;
|
||||
Dc.Meooe(src_e, phi_o); std::cout << GridLogMessage << "Applied Meo" << std::endl;
|
||||
|
||||
phi_o += chi_o;
|
||||
phi_e += chi_e;
|
||||
|
||||
setCheckerboard(phi, phi_e);
|
||||
setCheckerboard(phi, phi_o);
|
||||
|
||||
std::cout << GridLogDebug << "norm phi_e = " << norm2(phi_e) << " norm phi_o = " << norm2(phi_o) << " norm phi = " << norm2(phi) << std::endl;
|
||||
|
||||
diff = ref - phi;
|
||||
auto absDev = norm2(diff);
|
||||
auto relDev = absDev / norm2(ref);
|
||||
std::cout << GridLogMessage << "norm2(Dunprec), norm2(Deoprec), abs. deviation, rel. deviation: "
|
||||
<< norm2(ref) << " " << norm2(phi) << " " << absDev << " " << relDev
|
||||
<< " -> check " << ((relDev < checkTolerance) ? "passed" : "failed") << std::endl;
|
||||
assert(relDev <= checkTolerance);
|
||||
}
|
||||
|
||||
{
|
||||
std::cout << GridLogMessage << "===========================================================================" << std::endl;
|
||||
std::cout << GridLogMessage << "= Test MpcDagMpc is hermitian" << std::endl;
|
||||
std::cout << GridLogMessage << "===========================================================================" << std::endl;
|
||||
|
||||
CoarseVector phi(Grid_c); random(pRNG_c, phi);
|
||||
CoarseVector chi(Grid_c); random(pRNG_c, chi);
|
||||
CoarseVector chi_e(RBGrid_c); chi_e = Zero();
|
||||
CoarseVector chi_o(RBGrid_c); chi_o = Zero();
|
||||
CoarseVector dchi_e(RBGrid_c); dchi_e = Zero();
|
||||
CoarseVector dchi_o(RBGrid_c); dchi_o = Zero();
|
||||
CoarseVector phi_e(RBGrid_c); phi_e = Zero();
|
||||
CoarseVector phi_o(RBGrid_c); phi_o = Zero();
|
||||
CoarseVector dphi_e(RBGrid_c); dphi_e = Zero();
|
||||
CoarseVector dphi_o(RBGrid_c); dphi_o = Zero();
|
||||
|
||||
pickCheckerboard(Even, chi_e, chi);
|
||||
pickCheckerboard(Odd, chi_o, chi);
|
||||
pickCheckerboard(Even, phi_e, phi);
|
||||
pickCheckerboard(Odd, phi_o, phi);
|
||||
|
||||
SchurDiagMooeeOperator<CoarseDiracMatrix,CoarseVector> HermOpEO(Dc);
|
||||
|
||||
HermOpEO.MpcDagMpc(chi_e, dchi_e); std::cout << GridLogMessage << "Applied MpcDagMpc to chi_e" << std::endl;
|
||||
HermOpEO.MpcDagMpc(chi_o, dchi_o); std::cout << GridLogMessage << "Applied MpcDagMpc to chi_o" << std::endl;
|
||||
HermOpEO.MpcDagMpc(phi_e, dphi_e); std::cout << GridLogMessage << "Applied MpcDagMpc to phi_e" << std::endl;
|
||||
HermOpEO.MpcDagMpc(phi_o, dphi_o); std::cout << GridLogMessage << "Applied MpcDagMpc to phi_o" << std::endl;
|
||||
|
||||
ComplexD phiDchi_e = innerProduct(phi_e, dchi_e);
|
||||
ComplexD phiDchi_o = innerProduct(phi_o, dchi_o);
|
||||
ComplexD chiDphi_e = innerProduct(chi_e, dphi_e);
|
||||
ComplexD chiDphi_o = innerProduct(chi_o, dphi_o);
|
||||
|
||||
std::cout << GridLogMessage << "e " << phiDchi_e << " " << chiDphi_e << std::endl;
|
||||
std::cout << GridLogMessage << "o " << phiDchi_o << " " << chiDphi_o << std::endl;
|
||||
|
||||
std::cout << GridLogMessage << "phiDchi_e - conj(chiDphi_e) " << phiDchi_e - conj(chiDphi_e) << std::endl;
|
||||
std::cout << GridLogMessage << "phiDchi_o - conj(chiDphi_o) " << phiDchi_o - conj(chiDphi_o) << std::endl;
|
||||
}
|
||||
|
||||
Grid_finalize();
|
||||
}
|
||||
Reference in New Issue
Block a user