mirror of
https://github.com/paboyle/Grid.git
synced 2025-06-13 20:57:06 +01:00
Some small steps towards a multigrid
This commit is contained in:
Peter Boyle
@ -1,15 +1,11 @@
|
||||
|
||||
bin_PROGRAMS = Test_GaugeAction Test_Metropolis Test_cayley_cg Test_cayley_coarsen_support Test_cayley_even_odd Test_cayley_ldop_cg Test_cayley_ldop_cr Test_cf_coarsen_support Test_cf_cr_unprec Test_contfrac_cg Test_contfrac_even_odd Test_cshift Test_cshift_red_black Test_dwf_cg_prec Test_dwf_cg_schur Test_dwf_cg_unprec Test_dwf_cr_unprec Test_dwf_even_odd Test_dwf_fpgcr Test_gamma Test_lie_generators Test_main Test_multishift_sqrt Test_nersc_io Test_quenched_update Test_remez Test_rng Test_rng_fixed Test_simd Test_stencil Test_wilson_cg_prec Test_wilson_cg_schur Test_wilson_cg_unprec Test_wilson_cr_unprec Test_wilson_even_odd
|
||||
bin_PROGRAMS = Test_GaugeAction Test_cayley_cg Test_cayley_coarsen_support Test_cayley_even_odd Test_cayley_ldop_cr Test_cf_coarsen_support Test_cf_cr_unprec Test_contfrac_cg Test_contfrac_even_odd Test_cshift Test_cshift_red_black Test_dwf_cg_prec Test_dwf_cg_schur Test_dwf_cg_unprec Test_dwf_cr_unprec Test_dwf_even_odd Test_dwf_fpgcr Test_dwf_hdcr Test_gamma Test_lie_generators Test_main Test_multishift_sqrt Test_nersc_io Test_quenched_update Test_remez Test_rng Test_rng_fixed Test_simd Test_stencil Test_wilson_cg_prec Test_wilson_cg_schur Test_wilson_cg_unprec Test_wilson_cr_unprec Test_wilson_even_odd
|
||||
|
||||
|
||||
Test_GaugeAction_SOURCES=Test_GaugeAction.cc
|
||||
Test_GaugeAction_LDADD=-lGrid
|
||||
|
||||
|
||||
Test_Metropolis_SOURCES=Test_Metropolis.cc
|
||||
Test_Metropolis_LDADD=-lGrid
|
||||
|
||||
|
||||
Test_cayley_cg_SOURCES=Test_cayley_cg.cc
|
||||
Test_cayley_cg_LDADD=-lGrid
|
||||
|
||||
@ -22,10 +18,6 @@ Test_cayley_even_odd_SOURCES=Test_cayley_even_odd.cc
|
||||
Test_cayley_even_odd_LDADD=-lGrid
|
||||
|
||||
|
||||
Test_cayley_ldop_cg_SOURCES=Test_cayley_ldop_cg.cc
|
||||
Test_cayley_ldop_cg_LDADD=-lGrid
|
||||
|
||||
|
||||
Test_cayley_ldop_cr_SOURCES=Test_cayley_ldop_cr.cc
|
||||
Test_cayley_ldop_cr_LDADD=-lGrid
|
||||
|
||||
@ -78,6 +70,10 @@ Test_dwf_fpgcr_SOURCES=Test_dwf_fpgcr.cc
|
||||
Test_dwf_fpgcr_LDADD=-lGrid
|
||||
|
||||
|
||||
Test_dwf_hdcr_SOURCES=Test_dwf_hdcr.cc
|
||||
Test_dwf_hdcr_LDADD=-lGrid
|
||||
|
||||
|
||||
Test_gamma_SOURCES=Test_gamma.cc
|
||||
Test_gamma_LDADD=-lGrid
|
||||
|
||||
|
||||
@ -1,2 +0,0 @@
|
||||
f(x) = 6.81384+(-2.34645e-06/(x+0.000228091))+(-1.51593e-05/(x+0.00112084))+(-6.89254e-05/(x+0.003496))+(-0.000288983/(x+0.00954309))+(-0.00119277/(x+0.024928))+(-0.0050183/(x+0.0646627))+(-0.0226449/(x+0.171576))+(-0.123767/(x+0.491792))+(-1.1705/(x+1.78667))+(-102.992/(x+18.4866));
|
||||
f(x) = 0.14676+(0.00952992/(x+5.40933e-05))+(0.0115952/(x+0.000559699))+(0.0161824/(x+0.00203338))+(0.0243252/(x+0.00582831))+(0.0379533/(x+0.0154649))+(0.060699/(x+0.0401156))+(0.100345/(x+0.104788))+(0.178335/(x+0.286042))+(0.381586/(x+0.892189))+(1.42625/(x+4.38422));
|
||||
@ -1,155 +0,0 @@
|
||||
#include <Grid.h>
|
||||
#include <qcd/utils/WilsonLoops.h>
|
||||
#include <qcd/utils/SUn.h>
|
||||
|
||||
using namespace std;
|
||||
using namespace Grid;
|
||||
using namespace Grid::QCD;
|
||||
|
||||
|
||||
template<class d>
|
||||
struct scal {
|
||||
d internal;
|
||||
};
|
||||
|
||||
Gamma::GammaMatrix Gmu [] = {
|
||||
Gamma::GammaX,
|
||||
Gamma::GammaY,
|
||||
Gamma::GammaZ,
|
||||
Gamma::GammaT
|
||||
};
|
||||
|
||||
double lowpass(double x)
|
||||
{
|
||||
return pow(x*x+1.0,-2);
|
||||
}
|
||||
|
||||
int main (int argc, char ** argv)
|
||||
{
|
||||
Grid_init(&argc,&argv);
|
||||
|
||||
Chebyshev<LatticeFermion> filter(-150.0, 150.0,16, lowpass);
|
||||
ofstream csv(std::string("filter.dat"),std::ios::out|std::ios::trunc);
|
||||
filter.csv(csv);
|
||||
csv.close();
|
||||
|
||||
const int Ls=8;
|
||||
|
||||
GridCartesian * UGrid = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplexF::Nsimd()),GridDefaultMpi());
|
||||
GridRedBlackCartesian * UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
|
||||
|
||||
GridCartesian * FGrid = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid);
|
||||
GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid);
|
||||
|
||||
// Construct a coarsened grid
|
||||
std::vector<int> clatt = GridDefaultLatt();
|
||||
for(int d=0;d<clatt.size();d++){
|
||||
clatt[d] = clatt[d]/2;
|
||||
}
|
||||
GridCartesian *Coarse4d = SpaceTimeGrid::makeFourDimGrid(clatt, GridDefaultSimd(Nd,vComplexF::Nsimd()),GridDefaultMpi());;
|
||||
GridCartesian *Coarse5d = SpaceTimeGrid::makeFiveDimGrid(1,Coarse4d);
|
||||
|
||||
std::vector<int> seeds4({1,2,3,4});
|
||||
std::vector<int> seeds5({5,6,7,8});
|
||||
std::vector<int> cseeds({5,6,7,8});
|
||||
GridParallelRNG RNG5(FGrid); RNG5.SeedFixedIntegers(seeds5);
|
||||
GridParallelRNG RNG4(UGrid); RNG4.SeedFixedIntegers(seeds4);
|
||||
GridParallelRNG CRNG(Coarse5d);CRNG.SeedFixedIntegers(cseeds);
|
||||
|
||||
LatticeFermion src(FGrid); gaussian(RNG5,src);
|
||||
LatticeFermion result(FGrid); result=zero;
|
||||
LatticeFermion ref(FGrid); ref=zero;
|
||||
LatticeFermion tmp(FGrid);
|
||||
LatticeFermion err(FGrid);
|
||||
LatticeGaugeField Umu(UGrid);
|
||||
|
||||
|
||||
|
||||
//gaussian(RNG4,Umu);
|
||||
//random(RNG4,Umu);
|
||||
|
||||
NerscField header;
|
||||
std::string file("./ckpoint_lat.400");
|
||||
readNerscConfiguration(Umu,header,file);
|
||||
// SU3::ColdConfiguration(RNG4,Umu);
|
||||
// SU3::TepidConfiguration(RNG4,Umu);
|
||||
// SU3::HotConfiguration(RNG4,Umu);
|
||||
// Umu=zero;
|
||||
|
||||
#if 0
|
||||
LatticeColourMatrix U(UGrid);
|
||||
for(int nn=0;nn<Nd;nn++){
|
||||
U=peekIndex<LorentzIndex>(Umu,nn);
|
||||
U=U*adj(U)-1.0;
|
||||
std::cout<<"SU3 test "<<norm2(U)<<std::endl;
|
||||
}
|
||||
#endif
|
||||
RealD mass=0.1;
|
||||
RealD M5=1.5;
|
||||
|
||||
DomainWallFermion Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5);
|
||||
Gamma5R5HermitianLinearOperator<DomainWallFermion,LatticeFermion> HermIndefOp(Ddwf);
|
||||
|
||||
const int nbasis = 8;
|
||||
|
||||
#if 0
|
||||
std::vector<LatticeFermion> subspace(nbasis,FGrid);
|
||||
LatticeFermion noise(FGrid);
|
||||
LatticeFermion ms(FGrid);
|
||||
for(int b=0;b<nbasis;b++){
|
||||
|
||||
gaussian(RNG5,noise);
|
||||
RealD scale = pow(norm2(noise),-0.5);
|
||||
noise=noise*scale;
|
||||
|
||||
HermIndefOp.Op(noise,ms); std::cout << "Noise "<<b<<" Ms "<<norm2(ms)<< " "<< norm2(noise)<<std::endl;
|
||||
|
||||
// filter(HermIndefOp,noise,subspace[b]);
|
||||
// inverse iteration
|
||||
MdagMLinearOperator<DomainWallFermion,LatticeFermion> HermDefOp(Ddwf);
|
||||
ConjugateGradient<LatticeFermion> CG(1.0e-4,10000);
|
||||
|
||||
for(int i=0;i<1;i++){
|
||||
|
||||
CG(HermDefOp,noise,subspace[b]);
|
||||
noise = subspace[b];
|
||||
|
||||
scale = pow(norm2(noise),-0.5);
|
||||
noise=noise*scale;
|
||||
HermDefOp.Op(noise,ms); std::cout << "filt "<<b<<" <u|H|u> "<<norm2(ms)<< " "<< norm2(noise)<<std::endl;
|
||||
}
|
||||
|
||||
subspace[b] = noise;
|
||||
HermIndefOp.Op(subspace[b],ms); std::cout << "Filtered "<<b<<" Ms "<<norm2(ms)<< " "<<norm2(subspace[b]) <<std::endl;
|
||||
|
||||
}
|
||||
std::cout << "Computed randoms"<< std::endl;
|
||||
#else
|
||||
std::cout<<"Calling Aggregation class" <<std::endl;
|
||||
MdagMLinearOperator<DomainWallFermion,LatticeFermion> HermDefOp(Ddwf);
|
||||
typedef Aggregation<vSpinColourVector,vTComplex,nbasis> Subspace;
|
||||
Subspace Aggregates(Coarse5d,FGrid);
|
||||
Aggregates.CreateSubspace(RNG5,HermDefOp);
|
||||
std::cout << "Called aggregation class"<< std::endl;
|
||||
#endif
|
||||
|
||||
|
||||
typedef CoarsenedMatrix<vSpinColourVector,vTComplex,nbasis> LittleDiracOperator;
|
||||
typedef LittleDiracOperator::CoarseVector CoarseVector;
|
||||
|
||||
LittleDiracOperator LittleDiracOp(*Coarse5d);
|
||||
LittleDiracOp.CoarsenOperator(FGrid,HermIndefOp,Aggregates);
|
||||
|
||||
CoarseVector c_src (Coarse5d);
|
||||
CoarseVector c_res (Coarse5d);
|
||||
gaussian(CRNG,c_src);
|
||||
c_res=zero;
|
||||
|
||||
std::cout << "Solving CG on coarse space "<< std::endl;
|
||||
MdagMLinearOperator<LittleDiracOperator,CoarseVector> PosdefLdop(LittleDiracOp);
|
||||
ConjugateGradient<CoarseVector> CG(1.0e-6,10000);
|
||||
CG(PosdefLdop,c_src,c_res);
|
||||
|
||||
std::cout << "Done "<< std::endl;
|
||||
Grid_finalize();
|
||||
}
|
||||
@ -1,38 +1,13 @@
|
||||
#include <Grid.h>
|
||||
#include <qcd/utils/WilsonLoops.h>
|
||||
#include <qcd/utils/SUn.h>
|
||||
|
||||
using namespace std;
|
||||
using namespace Grid;
|
||||
using namespace Grid::QCD;
|
||||
|
||||
|
||||
template<class d>
|
||||
struct scal {
|
||||
d internal;
|
||||
};
|
||||
|
||||
Gamma::GammaMatrix Gmu [] = {
|
||||
Gamma::GammaX,
|
||||
Gamma::GammaY,
|
||||
Gamma::GammaZ,
|
||||
Gamma::GammaT
|
||||
};
|
||||
|
||||
double lowpass(double x)
|
||||
{
|
||||
return pow(x*x+1.0,-2);
|
||||
}
|
||||
|
||||
int main (int argc, char ** argv)
|
||||
{
|
||||
Grid_init(&argc,&argv);
|
||||
|
||||
Chebyshev<LatticeFermion> filter(-150.0, 150.0,16, lowpass);
|
||||
ofstream csv(std::string("filter.dat"),std::ios::out|std::ios::trunc);
|
||||
filter.csv(csv);
|
||||
csv.close();
|
||||
|
||||
const int Ls=8;
|
||||
|
||||
GridCartesian * UGrid = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplexF::Nsimd()),GridDefaultMpi());
|
||||
@ -41,7 +16,9 @@ int main (int argc, char ** argv)
|
||||
GridCartesian * FGrid = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid);
|
||||
GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid);
|
||||
|
||||
// Construct a coarsened grid
|
||||
///////////////////////////////////////////////////
|
||||
// Construct a coarsened grid; utility for this?
|
||||
///////////////////////////////////////////////////
|
||||
std::vector<int> clatt = GridDefaultLatt();
|
||||
for(int d=0;d<clatt.size();d++){
|
||||
clatt[d] = clatt[d]/2;
|
||||
@ -63,80 +40,38 @@ int main (int argc, char ** argv)
|
||||
LatticeFermion err(FGrid);
|
||||
LatticeGaugeField Umu(UGrid);
|
||||
|
||||
|
||||
|
||||
//gaussian(RNG4,Umu);
|
||||
//random(RNG4,Umu);
|
||||
|
||||
NerscField header;
|
||||
std::string file("./ckpoint_lat.400");
|
||||
readNerscConfiguration(Umu,header,file);
|
||||
|
||||
// SU3::ColdConfiguration(RNG4,Umu);
|
||||
// SU3::TepidConfiguration(RNG4,Umu);
|
||||
// SU3::HotConfiguration(RNG4,Umu);
|
||||
// Umu=zero;
|
||||
|
||||
#if 0
|
||||
LatticeColourMatrix U(UGrid);
|
||||
for(int nn=0;nn<Nd;nn++){
|
||||
U=peekIndex<LorentzIndex>(Umu,nn);
|
||||
U=U*adj(U)-1.0;
|
||||
std::cout<<"SU3 test "<<norm2(U)<<std::endl;
|
||||
}
|
||||
#endif
|
||||
RealD mass=0.1;
|
||||
RealD M5=1.5;
|
||||
|
||||
std::cout << "**************************************************"<< std::endl;
|
||||
std::cout << "Building g5R5 hermitian DWF operator" <<std::endl;
|
||||
std::cout << "**************************************************"<< std::endl;
|
||||
DomainWallFermion Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5);
|
||||
Gamma5R5HermitianLinearOperator<DomainWallFermion,LatticeFermion> HermIndefOp(Ddwf);
|
||||
|
||||
const int nbasis = 8;
|
||||
|
||||
#if 0
|
||||
std::vector<LatticeFermion> subspace(nbasis,FGrid);
|
||||
LatticeFermion noise(FGrid);
|
||||
LatticeFermion ms(FGrid);
|
||||
for(int b=0;b<nbasis;b++){
|
||||
typedef Aggregation<vSpinColourVector,vTComplex,nbasis> Subspace;
|
||||
typedef CoarsenedMatrix<vSpinColourVector,vTComplex,nbasis> LittleDiracOperator;
|
||||
typedef LittleDiracOperator::CoarseVector CoarseVector;
|
||||
|
||||
gaussian(RNG5,noise);
|
||||
RealD scale = pow(norm2(noise),-0.5);
|
||||
noise=noise*scale;
|
||||
|
||||
HermIndefOp.Op(noise,ms); std::cout << "Noise "<<b<<" Ms "<<norm2(ms)<< " "<< norm2(noise)<<std::endl;
|
||||
|
||||
// filter(HermIndefOp,noise,subspace[b]);
|
||||
// inverse iteration
|
||||
MdagMLinearOperator<DomainWallFermion,LatticeFermion> HermDefOp(Ddwf);
|
||||
ConjugateGradient<LatticeFermion> CG(1.0e-4,10000);
|
||||
|
||||
for(int i=0;i<1;i++){
|
||||
|
||||
CG(HermDefOp,noise,subspace[b]);
|
||||
noise = subspace[b];
|
||||
|
||||
scale = pow(norm2(noise),-0.5);
|
||||
noise=noise*scale;
|
||||
HermDefOp.Op(noise,ms); std::cout << "filt "<<b<<" <u|H|u> "<<norm2(ms)<< " "<< norm2(noise)<<std::endl;
|
||||
}
|
||||
|
||||
subspace[b] = noise;
|
||||
HermIndefOp.Op(subspace[b],ms); std::cout << "Filtered "<<b<<" Ms "<<norm2(ms)<< " "<<norm2(subspace[b]) <<std::endl;
|
||||
|
||||
}
|
||||
std::cout << "Computed randoms"<< std::endl;
|
||||
#else
|
||||
std::cout<<"Calling Aggregation class" <<std::endl;
|
||||
std::cout << "**************************************************"<< std::endl;
|
||||
std::cout << "Calling Aggregation class to build subspace" <<std::endl;
|
||||
std::cout << "**************************************************"<< std::endl;
|
||||
MdagMLinearOperator<DomainWallFermion,LatticeFermion> HermDefOp(Ddwf);
|
||||
typedef Aggregation<vSpinColourVector,vTComplex,nbasis> Subspace;
|
||||
Subspace Aggregates(Coarse5d,FGrid);
|
||||
Aggregates.CreateSubspace(RNG5,HermDefOp);
|
||||
std::cout << "Called aggregation class"<< std::endl;
|
||||
#endif
|
||||
|
||||
|
||||
typedef CoarsenedMatrix<vSpinColourVector,vTComplex,nbasis> LittleDiracOperator;
|
||||
typedef LittleDiracOperator::CoarseVector CoarseVector;
|
||||
|
||||
LittleDiracOperator LittleDiracOp(*Coarse5d);
|
||||
LittleDiracOp.CoarsenOperator(FGrid,HermIndefOp,Aggregates);
|
||||
|
||||
@ -145,18 +80,22 @@ int main (int argc, char ** argv)
|
||||
gaussian(CRNG,c_src);
|
||||
c_res=zero;
|
||||
|
||||
|
||||
std::cout << "Solving CG on coarse space "<< std::endl;
|
||||
std::cout << "**************************************************"<< std::endl;
|
||||
std::cout << "Solving mdagm-CG on coarse space "<< std::endl;
|
||||
std::cout << "**************************************************"<< std::endl;
|
||||
MdagMLinearOperator<LittleDiracOperator,CoarseVector> PosdefLdop(LittleDiracOp);
|
||||
ConjugateGradient<CoarseVector> CG(1.0e-6,10000);
|
||||
CG(PosdefLdop,c_src,c_res);
|
||||
|
||||
|
||||
std::cout << "Solving MCR on coarse space "<< std::endl;
|
||||
std::cout << "**************************************************"<< std::endl;
|
||||
std::cout << "Solving indef-MCR on coarse space "<< std::endl;
|
||||
std::cout << "**************************************************"<< std::endl;
|
||||
HermitianLinearOperator<LittleDiracOperator,CoarseVector> HermIndefLdop(LittleDiracOp);
|
||||
ConjugateResidual<CoarseVector> MCR(1.0e-6,10000);
|
||||
MCR(HermIndefLdop,c_src,c_res);
|
||||
|
||||
std::cout << "**************************************************"<< std::endl;
|
||||
std::cout << "Done "<< std::endl;
|
||||
std::cout << "**************************************************"<< std::endl;
|
||||
Grid_finalize();
|
||||
}
|
||||
|
||||
188
tests/Test_dwf_hdcr.cc
Normal file
188
tests/Test_dwf_hdcr.cc
Normal file
@ -0,0 +1,188 @@
|
||||
#include <Grid.h>
|
||||
#include <algorithms/iterative/PrecGeneralisedConjugateResidual.h>
|
||||
|
||||
using namespace std;
|
||||
using namespace Grid;
|
||||
using namespace Grid::QCD;
|
||||
|
||||
template<class Fobj,class CComplex,int nbasis>
|
||||
class MultiGridPreconditioner : public LinearFunction< Lattice<Fobj> > {
|
||||
public:
|
||||
|
||||
typedef Aggregation<Fobj,CComplex,nbasis> Aggregates;
|
||||
typedef CoarsenedMatrix<Fobj,CComplex,nbasis> CoarseOperator;
|
||||
|
||||
typedef typename Aggregation<Fobj,CComplex,nbasis>::siteVector siteVector;
|
||||
typedef typename Aggregation<Fobj,CComplex,nbasis>::CoarseScalar CoarseScalar;
|
||||
typedef typename Aggregation<Fobj,CComplex,nbasis>::CoarseVector CoarseVector;
|
||||
typedef typename Aggregation<Fobj,CComplex,nbasis>::CoarseMatrix CoarseMatrix;
|
||||
typedef typename Aggregation<Fobj,CComplex,nbasis>::FineField FineField;
|
||||
typedef LinearOperatorBase<FineField> FineOperator;
|
||||
|
||||
Aggregates & _Aggregates;
|
||||
CoarseOperator & _CoarseOperator;
|
||||
FineOperator & _FineOperator;
|
||||
|
||||
// Constructor
|
||||
MultiGridPreconditioner(Aggregates &Agg, CoarseOperator &Coarse, FineOperator &Fine)
|
||||
: _Aggregates(Agg),
|
||||
_CoarseOperator(Coarse),
|
||||
_FineOperator(Fine)
|
||||
{
|
||||
}
|
||||
|
||||
void operator()(const FineField &in, FineField & out) {
|
||||
|
||||
FineField Min(in._grid);
|
||||
|
||||
CoarseVector Csrc(_CoarseOperator.Grid());
|
||||
CoarseVector Ctmp(_CoarseOperator.Grid());
|
||||
CoarseVector Csol(_CoarseOperator.Grid());
|
||||
|
||||
// Monitor completeness of low mode space
|
||||
_Aggregates.ProjectToSubspace (Csrc,in);
|
||||
_Aggregates.PromoteFromSubspace(Csrc,out);
|
||||
std::cout<<"Completeness: "<<std::sqrt(norm2(out)/norm2(in))<<std::endl;
|
||||
|
||||
ConjugateResidual<FineField> MCR(1.0e-2,1000);
|
||||
ConjugateGradient<CoarseVector> CG(1.0e-2,10000);
|
||||
|
||||
// [PTM+Q] in = [1 - Q A] M in + Q in = Min + Q [ in -A Min]
|
||||
// Smoothing step, followed by coarse grid correction
|
||||
|
||||
MCR(_FineOperator,in,Min);
|
||||
_FineOperator.Op(Min,out);
|
||||
out = in -out; // out = in - A Min
|
||||
|
||||
MdagMLinearOperator<CoarseOperator,CoarseVector> MdagMOp(_CoarseOperator);
|
||||
HermitianLinearOperator<CoarseOperator,CoarseVector> HermOp(_CoarseOperator);
|
||||
Csol=zero;
|
||||
_Aggregates.ProjectToSubspace (Csrc,out);
|
||||
HermOp.AdjOp(Csrc,Ctmp);// Normal equations
|
||||
CG(MdagMOp ,Ctmp,Csol);
|
||||
_Aggregates.PromoteFromSubspace(Csol,out);
|
||||
|
||||
out = Min + out;;
|
||||
}
|
||||
|
||||
};
|
||||
|
||||
int main (int argc, char ** argv)
|
||||
{
|
||||
Grid_init(&argc,&argv);
|
||||
|
||||
const int Ls=8;
|
||||
|
||||
GridCartesian * UGrid = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplexF::Nsimd()),GridDefaultMpi());
|
||||
GridRedBlackCartesian * UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
|
||||
|
||||
GridCartesian * FGrid = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid);
|
||||
GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid);
|
||||
|
||||
///////////////////////////////////////////////////
|
||||
// Construct a coarsened grid; utility for this?
|
||||
///////////////////////////////////////////////////
|
||||
std::vector<int> clatt = GridDefaultLatt();
|
||||
for(int d=0;d<clatt.size();d++){
|
||||
clatt[d] = clatt[d]/4;
|
||||
}
|
||||
GridCartesian *Coarse4d = SpaceTimeGrid::makeFourDimGrid(clatt, GridDefaultSimd(Nd,vComplexF::Nsimd()),GridDefaultMpi());;
|
||||
GridCartesian *Coarse5d = SpaceTimeGrid::makeFiveDimGrid(1,Coarse4d);
|
||||
|
||||
std::vector<int> seeds4({1,2,3,4});
|
||||
std::vector<int> seeds5({5,6,7,8});
|
||||
std::vector<int> cseeds({5,6,7,8});
|
||||
GridParallelRNG RNG5(FGrid); RNG5.SeedFixedIntegers(seeds5);
|
||||
GridParallelRNG RNG4(UGrid); RNG4.SeedFixedIntegers(seeds4);
|
||||
GridParallelRNG CRNG(Coarse5d);CRNG.SeedFixedIntegers(cseeds);
|
||||
|
||||
LatticeFermion src(FGrid); gaussian(RNG5,src);
|
||||
LatticeFermion result(FGrid); result=zero;
|
||||
LatticeFermion ref(FGrid); ref=zero;
|
||||
LatticeFermion tmp(FGrid);
|
||||
LatticeFermion err(FGrid);
|
||||
LatticeGaugeField Umu(UGrid);
|
||||
|
||||
NerscField header;
|
||||
std::string file("./ckpoint_lat.4000");
|
||||
readNerscConfiguration(Umu,header,file);
|
||||
|
||||
// SU3::ColdConfiguration(RNG4,Umu);
|
||||
// SU3::TepidConfiguration(RNG4,Umu);
|
||||
// SU3::HotConfiguration(RNG4,Umu);
|
||||
// Umu=zero;
|
||||
|
||||
RealD mass=0.04;
|
||||
RealD M5=1.8;
|
||||
|
||||
std::cout << "**************************************************"<< std::endl;
|
||||
std::cout << "Building g5R5 hermitian DWF operator" <<std::endl;
|
||||
std::cout << "**************************************************"<< std::endl;
|
||||
DomainWallFermion Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5);
|
||||
|
||||
const int nbasis = 4;
|
||||
|
||||
typedef Aggregation<vSpinColourVector,vTComplex,nbasis> Subspace;
|
||||
typedef CoarsenedMatrix<vSpinColourVector,vTComplex,nbasis> CoarseOperator;
|
||||
typedef CoarseOperator::CoarseVector CoarseVector;
|
||||
|
||||
std::cout << "**************************************************"<< std::endl;
|
||||
std::cout << "Calling Aggregation class to build subspace" <<std::endl;
|
||||
std::cout << "**************************************************"<< std::endl;
|
||||
MdagMLinearOperator<DomainWallFermion,LatticeFermion> HermDefOp(Ddwf);
|
||||
Subspace Aggregates(Coarse5d,FGrid);
|
||||
Aggregates.CreateSubspace(RNG5,HermDefOp);
|
||||
|
||||
std::cout << "**************************************************"<< std::endl;
|
||||
std::cout << "Building coarse representation of Indef operator" <<std::endl;
|
||||
std::cout << "**************************************************"<< std::endl;
|
||||
Gamma5R5HermitianLinearOperator<DomainWallFermion,LatticeFermion> HermIndefOp(Ddwf);
|
||||
CoarsenedMatrix<vSpinColourVector,vTComplex,nbasis> LDOp(*Coarse5d);
|
||||
LDOp.CoarsenOperator(FGrid,HermIndefOp,Aggregates);
|
||||
|
||||
std::cout << "**************************************************"<< std::endl;
|
||||
std::cout << "Testing some coarse space solvers " <<std::endl;
|
||||
std::cout << "**************************************************"<< std::endl;
|
||||
CoarseVector c_src (Coarse5d);
|
||||
CoarseVector c_res (Coarse5d);
|
||||
gaussian(CRNG,c_src);
|
||||
c_res=zero;
|
||||
|
||||
std::cout << "**************************************************"<< std::endl;
|
||||
std::cout << "Solving posdef-CG on coarse space "<< std::endl;
|
||||
std::cout << "**************************************************"<< std::endl;
|
||||
MdagMLinearOperator<CoarseOperator,CoarseVector> PosdefLdop(LDOp);
|
||||
ConjugateGradient<CoarseVector> CG(1.0e-6,10000);
|
||||
CG(PosdefLdop,c_src,c_res);
|
||||
|
||||
std::cout << "**************************************************"<< std::endl;
|
||||
std::cout << "Solving indef-MCR on coarse space "<< std::endl;
|
||||
std::cout << "**************************************************"<< std::endl;
|
||||
HermitianLinearOperator<CoarseOperator,CoarseVector> HermIndefLdop(LDOp);
|
||||
ConjugateResidual<CoarseVector> MCR(1.0e-6,10000);
|
||||
//MCR(HermIndefLdop,c_src,c_res);
|
||||
|
||||
std::cout << "**************************************************"<< std::endl;
|
||||
std::cout << "Building deflation preconditioner "<< std::endl;
|
||||
std::cout << "**************************************************"<< std::endl;
|
||||
|
||||
MultiGridPreconditioner <vSpinColourVector,vTComplex,nbasis> Precon(Aggregates, LDOp,HermIndefOp);
|
||||
|
||||
std::cout << "**************************************************"<< std::endl;
|
||||
std::cout << "Building a one level PGCR "<< std::endl;
|
||||
std::cout << "**************************************************"<< std::endl;
|
||||
TrivialPrecon<LatticeFermion> simple;
|
||||
PrecGeneralisedConjugateResidual<LatticeFermion> GCR(1.0e-6,10000,simple,8,64);
|
||||
GCR(HermIndefOp,src,result);
|
||||
|
||||
std::cout << "**************************************************"<< std::endl;
|
||||
std::cout << "Building a two level PGCR "<< std::endl;
|
||||
std::cout << "**************************************************"<< std::endl;
|
||||
PrecGeneralisedConjugateResidual<LatticeFermion> PGCR(1.0e-6,10000,Precon,8,64);
|
||||
PGCR(HermIndefOp,src,result);
|
||||
|
||||
std::cout << "**************************************************"<< std::endl;
|
||||
std::cout << "Done "<< std::endl;
|
||||
std::cout << "**************************************************"<< std::endl;
|
||||
Grid_finalize();
|
||||
}
|
||||
Reference in New Issue
Block a user