mirror of
https://github.com/paboyle/Grid.git
synced 2024-11-10 07:55:35 +00:00
Merge pull request #16 from DanielRichtmann/feature/gpt-coarsenedmatrix
Enable checkerboard operations for CoarsenedMatrix
This commit is contained in:
Christoph Lehner
GitHub
commit
5534921bee
@ -31,6 +31,7 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
|
|||||||
#ifndef GRID_ALGORITHM_COARSENED_MATRIX_H
|
#ifndef GRID_ALGORITHM_COARSENED_MATRIX_H
|
||||||
#define GRID_ALGORITHM_COARSENED_MATRIX_H
|
#define GRID_ALGORITHM_COARSENED_MATRIX_H
|
||||||
|
|
||||||
|
#include <Grid/qcd/QCD.h> // needed for Dagger(Yes|No), Inverse(Yes|No)
|
||||||
|
|
||||||
NAMESPACE_BEGIN(Grid);
|
NAMESPACE_BEGIN(Grid);
|
||||||
|
|
||||||
@ -59,12 +60,14 @@ inline void blockMaskedInnerProduct(Lattice<CComplex> &CoarseInner,
|
|||||||
class Geometry {
|
class Geometry {
|
||||||
public:
|
public:
|
||||||
int npoint;
|
int npoint;
|
||||||
|
int base;
|
||||||
std::vector<int> directions ;
|
std::vector<int> directions ;
|
||||||
std::vector<int> displacements;
|
std::vector<int> displacements;
|
||||||
|
std::vector<int> points_dagger;
|
||||||
|
|
||||||
Geometry(int _d) {
|
Geometry(int _d) {
|
||||||
|
|
||||||
int base = (_d==5) ? 1:0;
|
base = (_d==5) ? 1:0;
|
||||||
|
|
||||||
// make coarse grid stencil for 4d , not 5d
|
// make coarse grid stencil for 4d , not 5d
|
||||||
if ( _d==5 ) _d=4;
|
if ( _d==5 ) _d=4;
|
||||||
@ -72,16 +75,51 @@ public:
|
|||||||
npoint = 2*_d+1;
|
npoint = 2*_d+1;
|
||||||
directions.resize(npoint);
|
directions.resize(npoint);
|
||||||
displacements.resize(npoint);
|
displacements.resize(npoint);
|
||||||
|
points_dagger.resize(npoint);
|
||||||
for(int d=0;d<_d;d++){
|
for(int d=0;d<_d;d++){
|
||||||
directions[d ] = d+base;
|
directions[d ] = d+base;
|
||||||
directions[d+_d] = d+base;
|
directions[d+_d] = d+base;
|
||||||
displacements[d ] = +1;
|
displacements[d ] = +1;
|
||||||
displacements[d+_d]= -1;
|
displacements[d+_d]= -1;
|
||||||
|
points_dagger[d ] = d+_d;
|
||||||
|
points_dagger[d+_d] = d;
|
||||||
}
|
}
|
||||||
directions [2*_d]=0;
|
directions [2*_d]=0;
|
||||||
displacements[2*_d]=0;
|
displacements[2*_d]=0;
|
||||||
|
points_dagger[2*_d]=2*_d;
|
||||||
}
|
}
|
||||||
|
|
||||||
|
int point(int dir, int disp) {
|
||||||
|
assert(disp == -1 || disp == 0 || disp == 1);
|
||||||
|
assert(base+0 <= dir && dir < base+4);
|
||||||
|
|
||||||
|
// directions faster index = new indexing
|
||||||
|
// 4d (base = 0):
|
||||||
|
// point 0 1 2 3 4 5 6 7 8
|
||||||
|
// dir 0 1 2 3 0 1 2 3 0
|
||||||
|
// disp +1 +1 +1 +1 -1 -1 -1 -1 0
|
||||||
|
// 5d (base = 1):
|
||||||
|
// point 0 1 2 3 4 5 6 7 8
|
||||||
|
// dir 1 2 3 4 1 2 3 4 0
|
||||||
|
// disp +1 +1 +1 +1 -1 -1 -1 -1 0
|
||||||
|
|
||||||
|
// displacements faster index = old indexing
|
||||||
|
// 4d (base = 0):
|
||||||
|
// point 0 1 2 3 4 5 6 7 8
|
||||||
|
// dir 0 0 1 1 2 2 3 3 0
|
||||||
|
// disp +1 -1 +1 -1 +1 -1 +1 -1 0
|
||||||
|
// 5d (base = 1):
|
||||||
|
// point 0 1 2 3 4 5 6 7 8
|
||||||
|
// dir 1 1 2 2 3 3 4 4 0
|
||||||
|
// disp +1 -1 +1 -1 +1 -1 +1 -1 0
|
||||||
|
|
||||||
|
if(dir == 0 and disp == 0)
|
||||||
|
return 8;
|
||||||
|
else // New indexing
|
||||||
|
return (1 - disp) / 2 * 4 + dir - base;
|
||||||
|
// else // Old indexing
|
||||||
|
// return (4 * (dir - base) + 1 - disp) / 2;
|
||||||
|
}
|
||||||
};
|
};
|
||||||
|
|
||||||
template<class Fobj,class CComplex,int nbasis>
|
template<class Fobj,class CComplex,int nbasis>
|
||||||
@ -258,7 +296,7 @@ public:
|
|||||||
// Fine Object == (per site) type of fine field
|
// Fine Object == (per site) type of fine field
|
||||||
// nbasis == number of deflation vectors
|
// nbasis == number of deflation vectors
|
||||||
template<class Fobj,class CComplex,int nbasis>
|
template<class Fobj,class CComplex,int nbasis>
|
||||||
class CoarsenedMatrix : public SparseMatrixBase<Lattice<iVector<CComplex,nbasis > > > {
|
class CoarsenedMatrix : public CheckerBoardedSparseMatrixBase<Lattice<iVector<CComplex,nbasis > > > {
|
||||||
public:
|
public:
|
||||||
|
|
||||||
typedef iVector<CComplex,nbasis > siteVector;
|
typedef iVector<CComplex,nbasis > siteVector;
|
||||||
@ -270,16 +308,7 @@ public:
|
|||||||
typedef Lattice<Fobj > FineField;
|
typedef Lattice<Fobj > FineField;
|
||||||
typedef CoarseVector FermionField;
|
typedef CoarseVector FermionField;
|
||||||
|
|
||||||
// enrich interface
|
// enrich interface, use default implementation as in FermionOperator ///////
|
||||||
void Dhop(CoarseVector const& in, CoarseVector& out, int dag) { assert(0); }
|
|
||||||
void DhopEO(CoarseVector const& in, CoarseVector& out, int dag) { assert(0); }
|
|
||||||
void DhopOE(CoarseVector const& in, CoarseVector& out, int dag) { assert(0); }
|
|
||||||
void Meooe(CoarseVector const& in, CoarseVector& out) { assert(0); }
|
|
||||||
void MeooeDag(CoarseVector const& in, CoarseVector& out) { assert(0); }
|
|
||||||
void Mooee(CoarseVector const& in, CoarseVector& out) { assert(0); }
|
|
||||||
void MooeeDag(CoarseVector const& in, CoarseVector& out) { assert(0); }
|
|
||||||
void MooeeInv(CoarseVector const& in, CoarseVector& out) { assert(0); }
|
|
||||||
void MooeeInvDag(CoarseVector const& in, CoarseVector& out) { assert(0); }
|
|
||||||
void Dminus(CoarseVector const& in, CoarseVector& out) { out = in; }
|
void Dminus(CoarseVector const& in, CoarseVector& out) { out = in; }
|
||||||
void DminusDag(CoarseVector const& in, CoarseVector& out) { out = in; }
|
void DminusDag(CoarseVector const& in, CoarseVector& out) { out = in; }
|
||||||
void ImportPhysicalFermionSource(CoarseVector const& input, CoarseVector& imported) { imported = input; }
|
void ImportPhysicalFermionSource(CoarseVector const& input, CoarseVector& imported) { imported = input; }
|
||||||
@ -292,21 +321,36 @@ public:
|
|||||||
////////////////////
|
////////////////////
|
||||||
Geometry geom;
|
Geometry geom;
|
||||||
GridBase * _grid;
|
GridBase * _grid;
|
||||||
|
GridBase* _cbgrid;
|
||||||
int hermitian;
|
int hermitian;
|
||||||
|
|
||||||
CartesianStencil<siteVector,siteVector,int> Stencil;
|
CartesianStencil<siteVector,siteVector,int> Stencil;
|
||||||
|
CartesianStencil<siteVector,siteVector,int> StencilEven;
|
||||||
|
CartesianStencil<siteVector,siteVector,int> StencilOdd;
|
||||||
|
|
||||||
std::vector<CoarseMatrix> A;
|
std::vector<CoarseMatrix> A;
|
||||||
|
std::vector<CoarseMatrix> Aeven;
|
||||||
|
std::vector<CoarseMatrix> Aodd;
|
||||||
|
|
||||||
|
CoarseMatrix AselfInv;
|
||||||
|
CoarseMatrix AselfInvEven;
|
||||||
|
CoarseMatrix AselfInvOdd;
|
||||||
|
|
||||||
|
Vector<RealD> dag_factor;
|
||||||
|
|
||||||
///////////////////////
|
///////////////////////
|
||||||
// Interface
|
// Interface
|
||||||
///////////////////////
|
///////////////////////
|
||||||
GridBase * Grid(void) { return _grid; }; // this is all the linalg routines need to know
|
GridBase * Grid(void) { return _grid; }; // this is all the linalg routines need to know
|
||||||
|
GridBase * RedBlackGrid() { return _cbgrid; };
|
||||||
|
|
||||||
|
int ConstEE() { return 0; }
|
||||||
|
|
||||||
void M (const CoarseVector &in, CoarseVector &out)
|
void M (const CoarseVector &in, CoarseVector &out)
|
||||||
{
|
{
|
||||||
conformable(_grid,in.Grid());
|
conformable(_grid,in.Grid());
|
||||||
conformable(in.Grid(),out.Grid());
|
conformable(in.Grid(),out.Grid());
|
||||||
|
out.Checkerboard() = in.Checkerboard();
|
||||||
|
|
||||||
SimpleCompressor<siteVector> compressor;
|
SimpleCompressor<siteVector> compressor;
|
||||||
|
|
||||||
@ -364,12 +408,72 @@ public:
|
|||||||
return M(in,out);
|
return M(in,out);
|
||||||
} else {
|
} else {
|
||||||
// corresponds to Galerkin coarsening
|
// corresponds to Galerkin coarsening
|
||||||
CoarseVector tmp(Grid());
|
return MdagNonHermitian(in, out);
|
||||||
G5C(tmp, in);
|
|
||||||
M(tmp, out);
|
|
||||||
G5C(out, out);
|
|
||||||
}
|
}
|
||||||
};
|
};
|
||||||
|
|
||||||
|
void MdagNonHermitian(const CoarseVector &in, CoarseVector &out)
|
||||||
|
{
|
||||||
|
conformable(_grid,in.Grid());
|
||||||
|
conformable(in.Grid(),out.Grid());
|
||||||
|
out.Checkerboard() = in.Checkerboard();
|
||||||
|
|
||||||
|
SimpleCompressor<siteVector> compressor;
|
||||||
|
|
||||||
|
Stencil.HaloExchange(in,compressor);
|
||||||
|
autoView( in_v , in, AcceleratorRead);
|
||||||
|
autoView( out_v , out, AcceleratorWrite);
|
||||||
|
autoView( Stencil_v , Stencil, AcceleratorRead);
|
||||||
|
auto& geom_v = geom;
|
||||||
|
typedef LatticeView<Cobj> Aview;
|
||||||
|
|
||||||
|
Vector<Aview> AcceleratorViewContainer;
|
||||||
|
|
||||||
|
for(int p=0;p<geom.npoint;p++) AcceleratorViewContainer.push_back(A[p].View(AcceleratorRead));
|
||||||
|
Aview *Aview_p = & AcceleratorViewContainer[0];
|
||||||
|
|
||||||
|
const int Nsimd = CComplex::Nsimd();
|
||||||
|
typedef decltype(coalescedRead(in_v[0])) calcVector;
|
||||||
|
typedef decltype(coalescedRead(in_v[0](0))) calcComplex;
|
||||||
|
|
||||||
|
int osites=Grid()->oSites();
|
||||||
|
|
||||||
|
Vector<int> points(geom.npoint, 0);
|
||||||
|
for(int p=0; p<geom.npoint; p++)
|
||||||
|
points[p] = geom.points_dagger[p];
|
||||||
|
|
||||||
|
RealD* dag_factor_p = &dag_factor[0];
|
||||||
|
|
||||||
|
accelerator_for(sss, Grid()->oSites()*nbasis, Nsimd, {
|
||||||
|
int ss = sss/nbasis;
|
||||||
|
int b = sss%nbasis;
|
||||||
|
calcComplex res = Zero();
|
||||||
|
calcVector nbr;
|
||||||
|
int ptype;
|
||||||
|
StencilEntry *SE;
|
||||||
|
|
||||||
|
for(int p=0;p<geom_v.npoint;p++){
|
||||||
|
int point = points[p];
|
||||||
|
|
||||||
|
SE=Stencil_v.GetEntry(ptype,point,ss);
|
||||||
|
|
||||||
|
if(SE->_is_local) {
|
||||||
|
nbr = coalescedReadPermute(in_v[SE->_offset],ptype,SE->_permute);
|
||||||
|
} else {
|
||||||
|
nbr = coalescedRead(Stencil_v.CommBuf()[SE->_offset]);
|
||||||
|
}
|
||||||
|
acceleratorSynchronise();
|
||||||
|
|
||||||
|
for(int bb=0;bb<nbasis;bb++) {
|
||||||
|
res = res + dag_factor_p[b*nbasis+bb]*coalescedRead(Aview_p[point][ss](b,bb))*nbr(bb);
|
||||||
|
}
|
||||||
|
}
|
||||||
|
coalescedWrite(out_v[ss](b),res);
|
||||||
|
});
|
||||||
|
|
||||||
|
for(int p=0;p<geom.npoint;p++) AcceleratorViewContainer[p].ViewClose();
|
||||||
|
}
|
||||||
|
|
||||||
void MdirComms(const CoarseVector &in)
|
void MdirComms(const CoarseVector &in)
|
||||||
{
|
{
|
||||||
SimpleCompressor<siteVector> compressor;
|
SimpleCompressor<siteVector> compressor;
|
||||||
@ -379,6 +483,7 @@ public:
|
|||||||
{
|
{
|
||||||
conformable(_grid,in.Grid());
|
conformable(_grid,in.Grid());
|
||||||
conformable(_grid,out.Grid());
|
conformable(_grid,out.Grid());
|
||||||
|
out.Checkerboard() = in.Checkerboard();
|
||||||
|
|
||||||
typedef LatticeView<Cobj> Aview;
|
typedef LatticeView<Cobj> Aview;
|
||||||
Vector<Aview> AcceleratorViewContainer;
|
Vector<Aview> AcceleratorViewContainer;
|
||||||
@ -434,34 +539,7 @@ public:
|
|||||||
|
|
||||||
this->MdirComms(in);
|
this->MdirComms(in);
|
||||||
|
|
||||||
int ndim = in.Grid()->Nd();
|
MdirCalc(in,out,geom.point(dir,disp));
|
||||||
|
|
||||||
//////////////
|
|
||||||
// 4D action like wilson
|
|
||||||
// 0+ => 0
|
|
||||||
// 0- => 4
|
|
||||||
// 1+ => 1
|
|
||||||
// 1- => 5
|
|
||||||
// etc..
|
|
||||||
//////////////
|
|
||||||
// 5D action like DWF
|
|
||||||
// 1+ => 0
|
|
||||||
// 1- => 4
|
|
||||||
// 2+ => 1
|
|
||||||
// 2- => 5
|
|
||||||
// etc..
|
|
||||||
auto point = [dir, disp, ndim](){
|
|
||||||
if(dir == 0 and disp == 0)
|
|
||||||
return 8;
|
|
||||||
else if ( ndim==4 ) {
|
|
||||||
return (1 - disp) / 2 * 4 + dir;
|
|
||||||
} else {
|
|
||||||
return (1 - disp) / 2 * 4 + dir - 1;
|
|
||||||
}
|
|
||||||
}();
|
|
||||||
|
|
||||||
MdirCalc(in,out,point);
|
|
||||||
|
|
||||||
};
|
};
|
||||||
|
|
||||||
void Mdiag(const CoarseVector &in, CoarseVector &out)
|
void Mdiag(const CoarseVector &in, CoarseVector &out)
|
||||||
@ -470,17 +548,269 @@ public:
|
|||||||
MdirCalc(in, out, point); // No comms
|
MdirCalc(in, out, point); // No comms
|
||||||
};
|
};
|
||||||
|
|
||||||
|
void Mooee(const CoarseVector &in, CoarseVector &out) {
|
||||||
|
MooeeInternal(in, out, DaggerNo, InverseNo);
|
||||||
|
}
|
||||||
|
|
||||||
|
void MooeeInv(const CoarseVector &in, CoarseVector &out) {
|
||||||
|
MooeeInternal(in, out, DaggerNo, InverseYes);
|
||||||
|
}
|
||||||
|
|
||||||
|
void MooeeDag(const CoarseVector &in, CoarseVector &out) {
|
||||||
|
MooeeInternal(in, out, DaggerYes, InverseNo);
|
||||||
|
}
|
||||||
|
|
||||||
|
void MooeeInvDag(const CoarseVector &in, CoarseVector &out) {
|
||||||
|
MooeeInternal(in, out, DaggerYes, InverseYes);
|
||||||
|
}
|
||||||
|
|
||||||
|
void Meooe(const CoarseVector &in, CoarseVector &out) {
|
||||||
|
if(in.Checkerboard() == Odd) {
|
||||||
|
DhopEO(in, out, DaggerNo);
|
||||||
|
} else {
|
||||||
|
DhopOE(in, out, DaggerNo);
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
void MeooeDag(const CoarseVector &in, CoarseVector &out) {
|
||||||
|
if(in.Checkerboard() == Odd) {
|
||||||
|
DhopEO(in, out, DaggerYes);
|
||||||
|
} else {
|
||||||
|
DhopOE(in, out, DaggerYes);
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
void Dhop(const CoarseVector &in, CoarseVector &out, int dag) {
|
||||||
|
conformable(in.Grid(), _grid); // verifies full grid
|
||||||
|
conformable(in.Grid(), out.Grid());
|
||||||
|
|
||||||
|
out.Checkerboard() = in.Checkerboard();
|
||||||
|
|
||||||
|
DhopInternal(Stencil, A, in, out, dag);
|
||||||
|
}
|
||||||
|
|
||||||
|
void DhopOE(const CoarseVector &in, CoarseVector &out, int dag) {
|
||||||
|
conformable(in.Grid(), _cbgrid); // verifies half grid
|
||||||
|
conformable(in.Grid(), out.Grid()); // drops the cb check
|
||||||
|
|
||||||
|
assert(in.Checkerboard() == Even);
|
||||||
|
out.Checkerboard() = Odd;
|
||||||
|
|
||||||
|
DhopInternal(StencilEven, Aodd, in, out, dag);
|
||||||
|
}
|
||||||
|
|
||||||
|
void DhopEO(const CoarseVector &in, CoarseVector &out, int dag) {
|
||||||
|
conformable(in.Grid(), _cbgrid); // verifies half grid
|
||||||
|
conformable(in.Grid(), out.Grid()); // drops the cb check
|
||||||
|
|
||||||
|
assert(in.Checkerboard() == Odd);
|
||||||
|
out.Checkerboard() = Even;
|
||||||
|
|
||||||
|
DhopInternal(StencilOdd, Aeven, in, out, dag);
|
||||||
|
}
|
||||||
|
|
||||||
|
void MooeeInternal(const CoarseVector &in, CoarseVector &out, int dag, int inv) {
|
||||||
|
out.Checkerboard() = in.Checkerboard();
|
||||||
|
assert(in.Checkerboard() == Odd || in.Checkerboard() == Even);
|
||||||
|
|
||||||
|
CoarseMatrix *Aself = nullptr;
|
||||||
|
if(in.Grid()->_isCheckerBoarded) {
|
||||||
|
if(in.Checkerboard() == Odd) {
|
||||||
|
Aself = (inv) ? &AselfInvOdd : &Aodd[geom.npoint-1];
|
||||||
|
DselfInternal(StencilOdd, *Aself, in, out, dag);
|
||||||
|
} else {
|
||||||
|
Aself = (inv) ? &AselfInvEven : &Aeven[geom.npoint-1];
|
||||||
|
DselfInternal(StencilEven, *Aself, in, out, dag);
|
||||||
|
}
|
||||||
|
} else {
|
||||||
|
Aself = (inv) ? &AselfInv : &A[geom.npoint-1];
|
||||||
|
DselfInternal(Stencil, *Aself, in, out, dag);
|
||||||
|
}
|
||||||
|
assert(Aself != nullptr);
|
||||||
|
}
|
||||||
|
|
||||||
|
void DselfInternal(CartesianStencil<siteVector,siteVector,int> &st, CoarseMatrix &a,
|
||||||
|
const CoarseVector &in, CoarseVector &out, int dag) {
|
||||||
|
int point = geom.npoint-1;
|
||||||
|
autoView( out_v, out, AcceleratorWrite);
|
||||||
|
autoView( in_v, in, AcceleratorRead);
|
||||||
|
autoView( st_v, st, AcceleratorRead);
|
||||||
|
autoView( a_v, a, AcceleratorRead);
|
||||||
|
|
||||||
|
const int Nsimd = CComplex::Nsimd();
|
||||||
|
typedef decltype(coalescedRead(in_v[0])) calcVector;
|
||||||
|
typedef decltype(coalescedRead(in_v[0](0))) calcComplex;
|
||||||
|
|
||||||
|
RealD* dag_factor_p = &dag_factor[0];
|
||||||
|
|
||||||
|
if(dag) {
|
||||||
|
accelerator_for(sss, in.Grid()->oSites()*nbasis, Nsimd, {
|
||||||
|
int ss = sss/nbasis;
|
||||||
|
int b = sss%nbasis;
|
||||||
|
calcComplex res = Zero();
|
||||||
|
calcVector nbr;
|
||||||
|
int ptype;
|
||||||
|
StencilEntry *SE;
|
||||||
|
|
||||||
|
SE=st_v.GetEntry(ptype,point,ss);
|
||||||
|
|
||||||
|
if(SE->_is_local) {
|
||||||
|
nbr = coalescedReadPermute(in_v[SE->_offset],ptype,SE->_permute);
|
||||||
|
} else {
|
||||||
|
nbr = coalescedRead(st_v.CommBuf()[SE->_offset]);
|
||||||
|
}
|
||||||
|
acceleratorSynchronise();
|
||||||
|
|
||||||
|
for(int bb=0;bb<nbasis;bb++) {
|
||||||
|
res = res + dag_factor_p[b*nbasis+bb]*coalescedRead(a_v[ss](b,bb))*nbr(bb);
|
||||||
|
}
|
||||||
|
coalescedWrite(out_v[ss](b),res);
|
||||||
|
});
|
||||||
|
} else {
|
||||||
|
accelerator_for(sss, in.Grid()->oSites()*nbasis, Nsimd, {
|
||||||
|
int ss = sss/nbasis;
|
||||||
|
int b = sss%nbasis;
|
||||||
|
calcComplex res = Zero();
|
||||||
|
calcVector nbr;
|
||||||
|
int ptype;
|
||||||
|
StencilEntry *SE;
|
||||||
|
|
||||||
|
SE=st_v.GetEntry(ptype,point,ss);
|
||||||
|
|
||||||
|
if(SE->_is_local) {
|
||||||
|
nbr = coalescedReadPermute(in_v[SE->_offset],ptype,SE->_permute);
|
||||||
|
} else {
|
||||||
|
nbr = coalescedRead(st_v.CommBuf()[SE->_offset]);
|
||||||
|
}
|
||||||
|
acceleratorSynchronise();
|
||||||
|
|
||||||
|
for(int bb=0;bb<nbasis;bb++) {
|
||||||
|
res = res + coalescedRead(a_v[ss](b,bb))*nbr(bb);
|
||||||
|
}
|
||||||
|
coalescedWrite(out_v[ss](b),res);
|
||||||
|
});
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
void DhopInternal(CartesianStencil<siteVector,siteVector,int> &st, std::vector<CoarseMatrix> &a,
|
||||||
|
const CoarseVector &in, CoarseVector &out, int dag) {
|
||||||
|
SimpleCompressor<siteVector> compressor;
|
||||||
|
|
||||||
|
st.HaloExchange(in,compressor);
|
||||||
|
autoView( in_v, in, AcceleratorRead);
|
||||||
|
autoView( out_v, out, AcceleratorWrite);
|
||||||
|
autoView( st_v , st, AcceleratorRead);
|
||||||
|
typedef LatticeView<Cobj> Aview;
|
||||||
|
|
||||||
|
// determine in what order we need the points
|
||||||
|
int npoint = geom.npoint-1;
|
||||||
|
Vector<int> points(npoint, 0);
|
||||||
|
for(int p=0; p<npoint; p++)
|
||||||
|
points[p] = (dag && !hermitian) ? geom.points_dagger[p] : p;
|
||||||
|
|
||||||
|
Vector<Aview> AcceleratorViewContainer;
|
||||||
|
for(int p=0;p<npoint;p++) AcceleratorViewContainer.push_back(a[p].View(AcceleratorRead));
|
||||||
|
Aview *Aview_p = & AcceleratorViewContainer[0];
|
||||||
|
|
||||||
|
const int Nsimd = CComplex::Nsimd();
|
||||||
|
typedef decltype(coalescedRead(in_v[0])) calcVector;
|
||||||
|
typedef decltype(coalescedRead(in_v[0](0))) calcComplex;
|
||||||
|
|
||||||
|
RealD* dag_factor_p = &dag_factor[0];
|
||||||
|
|
||||||
|
if(dag) {
|
||||||
|
accelerator_for(sss, in.Grid()->oSites()*nbasis, Nsimd, {
|
||||||
|
int ss = sss/nbasis;
|
||||||
|
int b = sss%nbasis;
|
||||||
|
calcComplex res = Zero();
|
||||||
|
calcVector nbr;
|
||||||
|
int ptype;
|
||||||
|
StencilEntry *SE;
|
||||||
|
|
||||||
|
for(int p=0;p<npoint;p++){
|
||||||
|
int point = points[p];
|
||||||
|
SE=st_v.GetEntry(ptype,point,ss);
|
||||||
|
|
||||||
|
if(SE->_is_local) {
|
||||||
|
nbr = coalescedReadPermute(in_v[SE->_offset],ptype,SE->_permute);
|
||||||
|
} else {
|
||||||
|
nbr = coalescedRead(st_v.CommBuf()[SE->_offset]);
|
||||||
|
}
|
||||||
|
acceleratorSynchronise();
|
||||||
|
|
||||||
|
for(int bb=0;bb<nbasis;bb++) {
|
||||||
|
res = res + dag_factor_p[b*nbasis+bb]*coalescedRead(Aview_p[point][ss](b,bb))*nbr(bb);
|
||||||
|
}
|
||||||
|
}
|
||||||
|
coalescedWrite(out_v[ss](b),res);
|
||||||
|
});
|
||||||
|
} else {
|
||||||
|
accelerator_for(sss, in.Grid()->oSites()*nbasis, Nsimd, {
|
||||||
|
int ss = sss/nbasis;
|
||||||
|
int b = sss%nbasis;
|
||||||
|
calcComplex res = Zero();
|
||||||
|
calcVector nbr;
|
||||||
|
int ptype;
|
||||||
|
StencilEntry *SE;
|
||||||
|
|
||||||
|
for(int p=0;p<npoint;p++){
|
||||||
|
int point = points[p];
|
||||||
|
SE=st_v.GetEntry(ptype,point,ss);
|
||||||
|
|
||||||
|
if(SE->_is_local) {
|
||||||
|
nbr = coalescedReadPermute(in_v[SE->_offset],ptype,SE->_permute);
|
||||||
|
} else {
|
||||||
|
nbr = coalescedRead(st_v.CommBuf()[SE->_offset]);
|
||||||
|
}
|
||||||
|
acceleratorSynchronise();
|
||||||
|
|
||||||
|
for(int bb=0;bb<nbasis;bb++) {
|
||||||
|
res = res + coalescedRead(Aview_p[point][ss](b,bb))*nbr(bb);
|
||||||
|
}
|
||||||
|
}
|
||||||
|
coalescedWrite(out_v[ss](b),res);
|
||||||
|
});
|
||||||
|
}
|
||||||
|
|
||||||
|
for(int p=0;p<npoint;p++) AcceleratorViewContainer[p].ViewClose();
|
||||||
|
}
|
||||||
|
|
||||||
CoarsenedMatrix(GridCartesian &CoarseGrid, int hermitian_=0) :
|
CoarsenedMatrix(GridCartesian &CoarseGrid, GridRedBlackCartesian &CoarseRBGrid, int hermitian_=0) :
|
||||||
|
|
||||||
_grid(&CoarseGrid),
|
_grid(&CoarseGrid),
|
||||||
|
_cbgrid(&CoarseRBGrid),
|
||||||
geom(CoarseGrid._ndimension),
|
geom(CoarseGrid._ndimension),
|
||||||
hermitian(hermitian_),
|
hermitian(hermitian_),
|
||||||
Stencil(&CoarseGrid,geom.npoint,Even,geom.directions,geom.displacements,0),
|
Stencil(&CoarseGrid,geom.npoint,Even,geom.directions,geom.displacements,0),
|
||||||
A(geom.npoint,&CoarseGrid)
|
StencilEven(&CoarseRBGrid,geom.npoint,Even,geom.directions,geom.displacements,0),
|
||||||
|
StencilOdd(&CoarseRBGrid,geom.npoint,Odd,geom.directions,geom.displacements,0),
|
||||||
|
A(geom.npoint,&CoarseGrid),
|
||||||
|
Aeven(geom.npoint,&CoarseRBGrid),
|
||||||
|
Aodd(geom.npoint,&CoarseRBGrid),
|
||||||
|
AselfInv(&CoarseGrid),
|
||||||
|
AselfInvEven(&CoarseRBGrid),
|
||||||
|
AselfInvOdd(&CoarseRBGrid),
|
||||||
|
dag_factor(nbasis*nbasis)
|
||||||
{
|
{
|
||||||
|
fillFactor();
|
||||||
};
|
};
|
||||||
|
|
||||||
|
void fillFactor() {
|
||||||
|
Eigen::MatrixXd dag_factor_eigen = Eigen::MatrixXd::Ones(nbasis, nbasis);
|
||||||
|
if(!hermitian) {
|
||||||
|
const int nb = nbasis/2;
|
||||||
|
dag_factor_eigen.block(0,nb,nb,nb) *= -1.0;
|
||||||
|
dag_factor_eigen.block(nb,0,nb,nb) *= -1.0;
|
||||||
|
}
|
||||||
|
|
||||||
|
// GPU readable prefactor
|
||||||
|
thread_for(i, nbasis*nbasis, {
|
||||||
|
int j = i/nbasis;
|
||||||
|
int k = i%nbasis;
|
||||||
|
dag_factor[i] = dag_factor_eigen(j, k);
|
||||||
|
});
|
||||||
|
}
|
||||||
|
|
||||||
void CoarsenOperator(GridBase *FineGrid,LinearOperatorBase<Lattice<Fobj> > &linop,
|
void CoarsenOperator(GridBase *FineGrid,LinearOperatorBase<Lattice<Fobj> > &linop,
|
||||||
Aggregation<Fobj,CComplex,nbasis> & Subspace)
|
Aggregation<Fobj,CComplex,nbasis> & Subspace)
|
||||||
{
|
{
|
||||||
@ -629,6 +959,9 @@ public:
|
|||||||
std::cout << GridLogMessage << " ForceHermitian, new code "<<std::endl;
|
std::cout << GridLogMessage << " ForceHermitian, new code "<<std::endl;
|
||||||
ForceHermitian();
|
ForceHermitian();
|
||||||
}
|
}
|
||||||
|
|
||||||
|
InvertSelfStencilLink(); std::cout << GridLogMessage << "Coarse self link inverted" << std::endl;
|
||||||
|
FillHalfCbs(); std::cout << GridLogMessage << "Coarse half checkerboards filled" << std::endl;
|
||||||
}
|
}
|
||||||
|
|
||||||
void ForceHermitian(void) {
|
void ForceHermitian(void) {
|
||||||
@ -650,6 +983,50 @@ public:
|
|||||||
}
|
}
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
|
void InvertSelfStencilLink() {
|
||||||
|
std::cout << GridLogDebug << "CoarsenedMatrix::InvertSelfStencilLink" << std::endl;
|
||||||
|
int localVolume = Grid()->lSites();
|
||||||
|
|
||||||
|
typedef typename Cobj::scalar_object scalar_object;
|
||||||
|
|
||||||
|
autoView(Aself_v, A[geom.npoint-1], CpuRead);
|
||||||
|
autoView(AselfInv_v, AselfInv, CpuWrite);
|
||||||
|
thread_for(site, localVolume, { // NOTE: Not able to bring this to GPU because of Eigen + peek/poke
|
||||||
|
Eigen::MatrixXcd selfLinkEigen = Eigen::MatrixXcd::Zero(nbasis, nbasis);
|
||||||
|
Eigen::MatrixXcd selfLinkInvEigen = Eigen::MatrixXcd::Zero(nbasis, nbasis);
|
||||||
|
|
||||||
|
scalar_object selfLink = Zero();
|
||||||
|
scalar_object selfLinkInv = Zero();
|
||||||
|
|
||||||
|
Coordinate lcoor;
|
||||||
|
|
||||||
|
Grid()->LocalIndexToLocalCoor(site, lcoor);
|
||||||
|
peekLocalSite(selfLink, Aself_v, lcoor);
|
||||||
|
|
||||||
|
for (int i = 0; i < nbasis; ++i)
|
||||||
|
for (int j = 0; j < nbasis; ++j)
|
||||||
|
selfLinkEigen(i, j) = static_cast<ComplexD>(TensorRemove(selfLink(i, j)));
|
||||||
|
|
||||||
|
selfLinkInvEigen = selfLinkEigen.inverse();
|
||||||
|
|
||||||
|
for(int i = 0; i < nbasis; ++i)
|
||||||
|
for(int j = 0; j < nbasis; ++j)
|
||||||
|
selfLinkInv(i, j) = selfLinkInvEigen(i, j);
|
||||||
|
|
||||||
|
pokeLocalSite(selfLinkInv, AselfInv_v, lcoor);
|
||||||
|
});
|
||||||
|
}
|
||||||
|
|
||||||
|
void FillHalfCbs() {
|
||||||
|
std::cout << GridLogDebug << "CoarsenedMatrix::FillHalfCbs" << std::endl;
|
||||||
|
for(int p = 0; p < geom.npoint; ++p) {
|
||||||
|
pickCheckerboard(Even, Aeven[p], A[p]);
|
||||||
|
pickCheckerboard(Odd, Aodd[p], A[p]);
|
||||||
|
}
|
||||||
|
pickCheckerboard(Even, AselfInvEven, AselfInv);
|
||||||
|
pickCheckerboard(Odd, AselfInvOdd, AselfInv);
|
||||||
|
}
|
||||||
};
|
};
|
||||||
|
|
||||||
NAMESPACE_END(Grid);
|
NAMESPACE_END(Grid);
|
||||||
|
|||||||
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();
|
||||||
|
}
|
||||||
Loading…
Reference in New Issue
Block a user