portelli/Grid
1
0
Fork 0
mirror of https://github.com/paboyle/Grid.git synced 2025-04-24 12:45:56 +01:00
Code Releases Activity

Schur factored matrix

Browse Source
This commit is contained in:
Quadro 2021-06-01 13:35:38 -04:00
parent 81bbd4e4ce
commit 426d2365d1
1 changed files with 125 additions and 47 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

172
Grid/qcd/action/fermion/SchurFactoredFermionOperator.h
View File

@ -27,6 +27,8 @@ Author: Peter Boyle <paboyle@ph.ed.ac.uk>
/* END LEGAL */ /* END LEGAL */
#pragma once #pragma once
#include <Grid/qcd/utils/MixedPrecisionOperatorFunction.h>
NAMESPACE_BEGIN(Grid); NAMESPACE_BEGIN(Grid);
//////////////////////////////////////////////////////// ////////////////////////////////////////////////////////
@ -63,47 +65,91 @@ NAMESPACE_BEGIN(Grid);
// //
//////////////////////////////////////////////////////// ////////////////////////////////////////////////////////
template<class Impl>
class SchurFactoredFermionOperator : public Impl
{
INHERIT_IMPL_TYPES(Impl);
template<class ImplD,class ImplF>
class SchurFactoredFermionOperator : public ImplD
{
INHERIT_IMPL_TYPES(ImplD);
typedef typename ImplF::FermionField FermionFieldF;
typedef typename ImplD::FermionField FermionFieldD;
typedef SchurDiagMooeeOperator<FermionOperator<ImplD>,FermionFieldD> LinearOperatorD;
typedef SchurDiagMooeeOperator<FermionOperator<ImplF>,FermionFieldF> LinearOperatorF;
typedef SchurDiagMooeeDagOperator<FermionOperator<ImplD>,FermionFieldD> LinearOperatorDagD;
typedef SchurDiagMooeeDagOperator<FermionOperator<ImplF>,FermionFieldF> LinearOperatorDagF;
typedef MixedPrecisionConjugateGradientOperatorFunction<FermionOperator<ImplD>,
FermionOperator<ImplF>,
LinearOperatorD,
LinearOperatorF> MxPCG;
typedef MixedPrecisionConjugateGradientOperatorFunction<FermionOperator<ImplD>,
FermionOperator<ImplF>,
LinearOperatorDagD,
LinearOperatorDagF> MxDagPCG;
public: public:
FermionOperator<Impl> & DirichletFermOp; GridBase *FermionGrid(void) { return PeriodicFermOpD.FermionGrid(); };
FermionOperator<Impl> & FermOp; GridBase *GaugeGrid(void) { return PeriodicFermOpD.GaugeGrid(); };
OperatorFunction<FermionField> &OmegaSolver;
OperatorFunction<FermionField> &OmegaDagSolver;
OperatorFunction<FermionField> &DSolver;
OperatorFunction<FermionField> &DdagSolver;
Coordinate Block;
SchurFactoredFermionOperator(FermionOperator<Impl> & _FermOp, FermionOperator<ImplD> & DirichletFermOpD;
FermionOperator<Impl> & _DirichletFermOp, FermionOperator<ImplF> & DirichletFermOpF;
OperatorFunction<FermionField> &_OmegaSolver, FermionOperator<ImplD> & PeriodicFermOpD;
OperatorFunction<FermionField> &_OmegaDagSolver, FermionOperator<ImplF> & PeriodicFermOpF;
OperatorFunction<FermionField> &_DSolver,
OperatorFunction<FermionField> &_DdagSolver, LinearOperatorD DirichletLinOpD;
LinearOperatorF DirichletLinOpF;
LinearOperatorD PeriodicLinOpD;
LinearOperatorF PeriodicLinOpF;
LinearOperatorDagD DirichletLinOpDagD;
LinearOperatorDagF DirichletLinOpDagF;
LinearOperatorDagD PeriodicLinOpDagD;
LinearOperatorDagF PeriodicLinOpDagF;
// Can tinker with these in the pseudofermion for force vs. action solves
Integer maxinnerit;
Integer maxouterit;
RealD tol;
Coordinate Block;
Coordinate InnerBlock;
SchurFactoredFermionOperator(FermionOperator<ImplD> & _PeriodicFermOpD,
FermionOperator<ImplF> & _PeriodicFermOpF,
FermionOperator<ImplD> & _DirichletFermOpD,
FermionOperator<ImplF> & _DirichletFermOpF,
Coordinate &_Block) Coordinate &_Block)
: Block(_Block), : Block(_Block),
FermOp(_FermOp), PeriodicFermOpD(_PeriodicFermOpD),
DirichletFermOp(_DirichletFermOp), PeriodicFermOpF(_PeriodicFermOpF),
OmegaSolver(_OmegaSolver), DirichletFermOpD(_DirichletFermOpD),
OmegaDagSolver(_OmegaDagSolver), DirichletFermOpF(_DirichletFermOpF),
DSolver(_DSolver), DirichletLinOpD(DirichletFermOpD),
DdagSolver(_DdagSolver) DirichletLinOpF(DirichletFermOpF),
PeriodicLinOpD(PeriodicFermOpD),
PeriodicLinOpF(PeriodicFermOpF),
DirichletLinOpDagD(DirichletFermOpD),
DirichletLinOpDagF(DirichletFermOpF),
PeriodicLinOpDagD(PeriodicFermOpD),
PeriodicLinOpDagF(PeriodicFermOpF)
{ {
// Pass in Dirichlet FermOp because we really need two dirac operators InnerBlock = Coordinate({16,16,16,20});
// as double stored gauge fields differ and they will otherwise overwrite tol=1.0e-10;
assert(_FermOp.FermionGrid() == _DirichletFermOp.FermionGrid()); // May not be true in future if change communicator scheme maxinnerit=1000;
maxouterit=10;
assert(PeriodicFermOpD.FermionGrid() == DirichletFermOpD.FermionGrid());
assert(PeriodicFermOpF.FermionGrid() == DirichletFermOpF.FermionGrid());
}; };
enum Domain { Omega=0, OmegaBar=1 }; enum Domain { Omega=0, OmegaBar=1 };
void ImportGauge(const GaugeField &Umu) void ImportGauge(const GaugeField &Umu)
{ {
FermOp.ImportGauge(Umu); PeriodicFermOpD.ImportGauge(Umu);
DirichletFermOp.ImportGauge(Umu); DirichletFermOpD.ImportGauge(Umu);
// Single precision will update in the mixed prec CG
} }
void ProjectBoundaryBothDomains (FermionField &f,int sgn) void ProjectBoundaryBothDomains (FermionField &f,int sgn)
{ {
@ -202,56 +248,56 @@ public:
{ {
FermionField tmp(in); FermionField tmp(in);
ProjectOmegaBar(tmp); ProjectOmegaBar(tmp);
FermOp.M(tmp,out); PeriodicFermOpD.M(tmp,out);
ProjectOmega(out); ProjectOmega(out);
}; };
void dBoundaryDag (FermionField &in,FermionField &out) void dBoundaryDag (FermionField &in,FermionField &out)
{ {
FermionField tmp(in); FermionField tmp(in);
ProjectOmega(tmp); ProjectOmega(tmp);
FermOp.Mdag(tmp,out); PeriodicFermOpD.Mdag(tmp,out);
ProjectOmegaBar(out); ProjectOmegaBar(out);
}; };
void dBoundaryBar (FermionField &in,FermionField &out) void dBoundaryBar (FermionField &in,FermionField &out)
{ {
FermionField tmp(in); FermionField tmp(in);
ProjectOmega(tmp); ProjectOmega(tmp);
FermOp.M(tmp,out); PeriodicFermOpD.M(tmp,out);
ProjectOmegaBar(out); ProjectOmegaBar(out);
}; };
void dBoundaryBarDag (FermionField &in,FermionField &out) void dBoundaryBarDag (FermionField &in,FermionField &out)
{ {
FermionField tmp(in); FermionField tmp(in);
ProjectOmegaBar(tmp); ProjectOmegaBar(tmp);
FermOp.Mdag(tmp,out); PeriodicFermOpD.Mdag(tmp,out);
ProjectOmega(out); ProjectOmega(out);
}; };
void dOmega (FermionField &in,FermionField &out) void dOmega (FermionField &in,FermionField &out)
{ {
FermionField tmp(in); FermionField tmp(in);
ProjectOmega(tmp); ProjectOmega(tmp);
DirichletFermOp.M(tmp,out); DirichletFermOpD.M(tmp,out);
ProjectOmega(out); ProjectOmega(out);
}; };
void dOmegaBar (FermionField &in,FermionField &out) void dOmegaBar (FermionField &in,FermionField &out)
{ {
FermionField tmp(in); FermionField tmp(in);
ProjectOmegaBar(tmp); ProjectOmegaBar(tmp);
DirichletFermOp.M(tmp,out); DirichletFermOpD.M(tmp,out);
ProjectOmegaBar(out); ProjectOmegaBar(out);
}; };
void dOmegaDag (FermionField &in,FermionField &out) void dOmegaDag (FermionField &in,FermionField &out)
{ {
FermionField tmp(in); FermionField tmp(in);
ProjectOmega(tmp); ProjectOmega(tmp);
DirichletFermOp.Mdag(tmp,out); DirichletFermOpD.Mdag(tmp,out);
ProjectOmega(out); ProjectOmega(out);
}; };
void dOmegaBarDag (FermionField &in,FermionField &out) void dOmegaBarDag (FermionField &in,FermionField &out)
{ {
FermionField tmp(in); FermionField tmp(in);
ProjectOmegaBar(tmp); ProjectOmegaBar(tmp);
DirichletFermOp.Mdag(tmp,out); DirichletFermOpD.Mdag(tmp,out);
ProjectOmegaBar(out); ProjectOmegaBar(out);
}; };
void dOmegaInv (FermionField &in,FermionField &out) void dOmegaInv (FermionField &in,FermionField &out)
@ -284,20 +330,36 @@ public:
}; };
void dOmegaInvAndOmegaBarInv(FermionField &in,FermionField &out) void dOmegaInvAndOmegaBarInv(FermionField &in,FermionField &out)
{ {
MxPCG OmegaSolver(tol,
maxinnerit,
maxouterit,
DirichletFermOpF.FermionRedBlackGrid(),
DirichletFermOpF,
DirichletFermOpD,
DirichletLinOpF,
DirichletLinOpD);
SchurRedBlackDiagMooeeSolve<FermionField> PrecSolve(OmegaSolver); SchurRedBlackDiagMooeeSolve<FermionField> PrecSolve(OmegaSolver);
PrecSolve(DirichletFermOp,in,out); PrecSolve(DirichletFermOpD,in,out);
}; };
void dOmegaDagInvAndOmegaBarDagInv(FermionField &in,FermionField &out) void dOmegaDagInvAndOmegaBarDagInv(FermionField &in,FermionField &out)
{ {
MxDagPCG OmegaDagSolver(tol,
maxinnerit,
maxouterit,
DirichletFermOpF.FermionRedBlackGrid(),
DirichletFermOpF,
DirichletFermOpD,
DirichletLinOpDagF,
DirichletLinOpDagD);
SchurRedBlackDiagMooeeDagSolve<FermionField> PrecSolve(OmegaDagSolver); SchurRedBlackDiagMooeeDagSolve<FermionField> PrecSolve(OmegaDagSolver);
PrecSolve(DirichletFermOp,in,out); PrecSolve(DirichletFermOpD,in,out);
}; };
// Rdag = Pdbar - DdbarDag DomegabarDagInv DdDag DomegaDagInv Pdbar // Rdag = Pdbar - DdbarDag DomegabarDagInv DdDag DomegaDagInv Pdbar
void RDag(FermionField &in,FermionField &out) void RDag(FermionField &in,FermionField &out)
{ {
FermionField tmp1(FermOp.FermionGrid()); FermionField tmp1(PeriodicFermOpD.FermionGrid());
FermionField tmp2(FermOp.FermionGrid()); FermionField tmp2(PeriodicFermOpD.FermionGrid());
out = in; out = in;
ProjectBoundaryBar(out); ProjectBoundaryBar(out);
dOmegaDagInv(out,tmp1); dOmegaDagInv(out,tmp1);
@ -310,8 +372,8 @@ public:
// R = Pdbar - Pdbar DomegaInv Dd DomegabarInv Ddbar // R = Pdbar - Pdbar DomegaInv Dd DomegabarInv Ddbar
void R(FermionField &in,FermionField &out) void R(FermionField &in,FermionField &out)
{ {
FermionField tmp1(FermOp.FermionGrid()); FermionField tmp1(PeriodicFermOpD.FermionGrid());
FermionField tmp2(FermOp.FermionGrid()); FermionField tmp2(PeriodicFermOpD.FermionGrid());
out = in; out = in;
ProjectBoundaryBar(out); ProjectBoundaryBar(out);
dBoundaryBar(out,tmp1); dBoundaryBar(out,tmp1);
@ -326,7 +388,7 @@ public:
// R = Pdbar - Pdbar Dinv Ddbar // R = Pdbar - Pdbar Dinv Ddbar
void RInv(FermionField &in,FermionField &out) void RInv(FermionField &in,FermionField &out)
{ {
FermionField tmp1(FermOp.FermionGrid()); FermionField tmp1(PeriodicFermOpD.FermionGrid());
dBoundaryBar(in,out); dBoundaryBar(in,out);
Dinverse(out,tmp1); Dinverse(out,tmp1);
out =in -tmp1; out =in -tmp1;
@ -335,8 +397,8 @@ public:
// R = Pdbar - DdbarDag DinvDag Pdbar // R = Pdbar - DdbarDag DinvDag Pdbar
void RDagInv(FermionField &in,FermionField &out) void RDagInv(FermionField &in,FermionField &out)
{ {
FermionField tmp(FermOp.FermionGrid()); FermionField tmp(PeriodicFermOpD.FermionGrid());
FermionField Pin(FermOp.FermionGrid()); FermionField Pin(PeriodicFermOpD.FermionGrid());
Pin = in; ProjectBoundaryBar(Pin); Pin = in; ProjectBoundaryBar(Pin);
DinverseDag(Pin,out); DinverseDag(Pin,out);
dBoundaryBarDag(out,tmp); dBoundaryBarDag(out,tmp);
@ -345,13 +407,29 @@ public:
// Non-dirichlet inverter using red-black preconditioning // Non-dirichlet inverter using red-black preconditioning
void Dinverse(FermionField &in,FermionField &out) void Dinverse(FermionField &in,FermionField &out)
{ {
MxPCG DSolver(tol,
maxinnerit,
maxouterit,
PeriodicFermOpF.FermionRedBlackGrid(),
PeriodicFermOpF,
PeriodicFermOpD,
PeriodicLinOpF,
PeriodicLinOpD);
SchurRedBlackDiagMooeeSolve<FermionField> Solve(DSolver); SchurRedBlackDiagMooeeSolve<FermionField> Solve(DSolver);
Solve(FermOp,in,out); Solve(PeriodicFermOpD,in,out);
} }
void DinverseDag(FermionField &in,FermionField &out) void DinverseDag(FermionField &in,FermionField &out)
{ {
MxDagPCG DdagSolver(tol,
maxinnerit,
maxouterit,
PeriodicFermOpF.FermionRedBlackGrid(),
PeriodicFermOpF,
PeriodicFermOpD,
PeriodicLinOpDagF,
PeriodicLinOpDagD);
SchurRedBlackDiagMooeeDagSolve<FermionField> Solve(DdagSolver); SchurRedBlackDiagMooeeDagSolve<FermionField> Solve(DdagSolver);
Solve(FermOp,in,out); Solve(PeriodicFermOpD,in,out);
} }
}; };