mirror of
https://github.com/paboyle/Grid.git
synced 2024-11-10 07:55:35 +00:00
WilsonMG: Major cleanup
This commit is contained in:
Daniel Richtmann
parent
a3e009ba54
commit
c4ce70a821
@ -64,8 +64,8 @@ public:
|
||||
std::cout << GridLogMessage << std::scientific << std::setprecision(2) << std::setw(2) << std::showpos << "vector " << i << ": "
|
||||
<< "singular value: " << lambda << ", singular vector precision: " << mu << ", norm: " << nrm << std::endl;
|
||||
}
|
||||
std::cout << GridLogMessage << std::scientific << std::setprecision(2) << std::setw(2) << std::showpos << positiveOnes << " out of "
|
||||
<< nn << " vectors were positive" << std::endl;
|
||||
std::cout << GridLogMessage << std::scientific << std::setprecision(2) << std::setw(2) << std::showpos << positiveOnes << " out of " << nn
|
||||
<< " vectors were positive" << std::endl;
|
||||
}
|
||||
};
|
||||
|
||||
@ -87,10 +87,6 @@ public:
|
||||
};
|
||||
myclass params;
|
||||
|
||||
RealD InverseApproximation(RealD x) {
|
||||
return 1.0 / x;
|
||||
}
|
||||
|
||||
template<int nbasis> struct CoarseGrids {
|
||||
public:
|
||||
std::vector<std::vector<int>> LattSizes;
|
||||
@ -186,11 +182,10 @@ template<class Field> void testLinearOperator(LinearOperatorBase<Field> &LinOp,
|
||||
ComplexD phiMPhi = innerProduct(phi, MPhi);
|
||||
ComplexD chiMdagChi = innerProduct(chi, MdagChi);
|
||||
|
||||
std::cout << GridLogMessage << " chiMPhi = " << chiMPhi << " phiMdagChi = " << phiMdagChi
|
||||
<< " difference = " << chiMPhi - conjugate(phiMdagChi) << std::endl;
|
||||
|
||||
std::cout << GridLogMessage << " phiMPhi = " << phiMPhi << " chiMdagChi = " << chiMdagChi << " <- should be real if hermitian"
|
||||
std::cout << GridLogMessage << " chiMPhi = " << chiMPhi << " phiMdagChi = " << phiMdagChi << " difference = " << chiMPhi - conjugate(phiMdagChi)
|
||||
<< std::endl;
|
||||
|
||||
std::cout << GridLogMessage << " phiMPhi = " << phiMPhi << " chiMdagChi = " << chiMdagChi << " <- should be real if hermitian" << std::endl;
|
||||
}
|
||||
|
||||
{
|
||||
@ -255,374 +250,53 @@ public:
|
||||
, _SmootherOperator(Smooth)
|
||||
, _SmootherMatrix(SmootherMatrix) {}
|
||||
|
||||
void PowerMethod(const FineField &in) {
|
||||
|
||||
FineField p1(in._grid);
|
||||
FineField p2(in._grid);
|
||||
|
||||
MdagMLinearOperator<Matrix, FineField> fMdagMOp(_FineMatrix);
|
||||
|
||||
p1 = in;
|
||||
RealD absp2;
|
||||
for(int i = 0; i < 20; i++) {
|
||||
RealD absp1 = std::sqrt(norm2(p1));
|
||||
fMdagMOp.HermOp(p1, p2); // this is the G5 herm bit
|
||||
// _FineOperator.Op(p1,p2); // this is the G5 herm bit
|
||||
RealD absp2 = std::sqrt(norm2(p2));
|
||||
if(i % 10 == 9)
|
||||
std::cout << GridLogMessage << "Power method on mdagm " << i << " " << absp2 / absp1 << std::endl;
|
||||
p1 = p2 * (1.0 / std::sqrt(absp2));
|
||||
}
|
||||
}
|
||||
|
||||
void operator()(const FineField &in, FineField &out) {
|
||||
if(params.domaindecompose) {
|
||||
operatorSAP(in, out);
|
||||
} else {
|
||||
operatorCheby(in, out);
|
||||
}
|
||||
}
|
||||
|
||||
////////////////////////////////////////////////////////////////////////
|
||||
// ADEF2: [PTM+Q] in = [1 - Q A] M in + Q in = Min + Q [ in -A Min]
|
||||
// ADEF1: [MP+Q ] in = M [1 - A Q] in + Q in
|
||||
////////////////////////////////////////////////////////////////////////
|
||||
#if 1
|
||||
void operatorADEF2(const FineField &in, FineField &out) {
|
||||
CoarseVector coarseSrc(_CoarseOperator.Grid());
|
||||
CoarseVector coarseTmp(_CoarseOperator.Grid());
|
||||
CoarseVector coarseSol(_CoarseOperator.Grid());
|
||||
coarseSol = zero;
|
||||
|
||||
CoarseVector Csrc(_CoarseOperator.Grid());
|
||||
CoarseVector Ctmp(_CoarseOperator.Grid());
|
||||
CoarseVector Csol(_CoarseOperator.Grid());
|
||||
GeneralisedMinimalResidual<CoarseVector> coarseGMRES(5.0e-2, 100, 25, false);
|
||||
GeneralisedMinimalResidual<FineField> fineGMRES(5.0e-2, 100, 25, false);
|
||||
|
||||
ConjugateGradient<CoarseVector> CG(1.0e-10, 100000);
|
||||
ConjugateGradient<FineField> fCG(3.0e-2, 1000);
|
||||
HermitianLinearOperator<CoarseOperator, CoarseVector> coarseHermOp(_CoarseOperator);
|
||||
MdagMLinearOperator<CoarseOperator, CoarseVector> coarseMdagMOp(_CoarseOperator);
|
||||
MdagMLinearOperator<Matrix, FineField> fineMdagMOp(_SmootherMatrix);
|
||||
|
||||
HermitianLinearOperator<CoarseOperator, CoarseVector> HermOp(_CoarseOperator);
|
||||
MdagMLinearOperator<CoarseOperator, CoarseVector> MdagMOp(_CoarseOperator);
|
||||
MdagMLinearOperator<Matrix, FineField> fMdagMOp(_FineMatrix);
|
||||
|
||||
FineField tmp(in._grid);
|
||||
FineField res(in._grid);
|
||||
FineField Min(in._grid);
|
||||
|
||||
// Monitor completeness of low mode space
|
||||
_Aggregates.ProjectToSubspace(Csrc, in);
|
||||
_Aggregates.PromoteFromSubspace(Csrc, out);
|
||||
std::cout << GridLogMessage << "Coarse Grid Preconditioner\nCompleteness in: " << std::sqrt(norm2(out) / norm2(in)) << std::endl;
|
||||
|
||||
// [PTM+Q] in = [1 - Q A] M in + Q in = Min + Q [ in -A Min]
|
||||
_FineOperator.Op(in, tmp); // this is the G5 herm bit
|
||||
fCG(fMdagMOp, tmp, Min); // solves MdagM = g5 M g5M
|
||||
|
||||
// Monitor completeness of low mode space
|
||||
_Aggregates.ProjectToSubspace(Csrc, Min);
|
||||
_Aggregates.PromoteFromSubspace(Csrc, out);
|
||||
std::cout << GridLogMessage << "Completeness Min: " << std::sqrt(norm2(out) / norm2(Min)) << std::endl;
|
||||
|
||||
_FineOperator.Op(Min, tmp);
|
||||
tmp = in - tmp; // in - A Min
|
||||
|
||||
Csol = zero;
|
||||
_Aggregates.ProjectToSubspace(Csrc, tmp);
|
||||
HermOp.AdjOp(Csrc, Ctmp); // Normal equations
|
||||
CG(MdagMOp, Ctmp, Csol);
|
||||
|
||||
HermOp.Op(Csol, Ctmp);
|
||||
Ctmp = Ctmp - Csrc;
|
||||
std::cout << GridLogMessage << "coarse space true residual " << std::sqrt(norm2(Ctmp) / norm2(Csrc)) << std::endl;
|
||||
_Aggregates.PromoteFromSubspace(Csol, out);
|
||||
|
||||
_FineOperator.Op(out, res);
|
||||
res = res - tmp;
|
||||
std::cout << GridLogMessage << "promoted sol residual " << std::sqrt(norm2(res) / norm2(tmp)) << std::endl;
|
||||
_Aggregates.ProjectToSubspace(Csrc, res);
|
||||
std::cout << GridLogMessage << "coarse space proj of residual " << norm2(Csrc) << std::endl;
|
||||
|
||||
out = out + Min; // additive coarse space correction
|
||||
// out = Min; // no additive coarse space correction
|
||||
|
||||
_FineOperator.Op(out, tmp);
|
||||
tmp = tmp - in; // tmp is new residual
|
||||
|
||||
std::cout << GridLogMessage << " Preconditioner in " << norm2(in) << std::endl;
|
||||
std::cout << GridLogMessage << " Preconditioner out " << norm2(out) << std::endl;
|
||||
std::cout << GridLogMessage << "preconditioner thinks residual is " << std::sqrt(norm2(tmp) / norm2(in)) << std::endl;
|
||||
}
|
||||
#endif
|
||||
// ADEF1: [MP+Q ] in = M [1 - A Q] in + Q in
|
||||
#if 1
|
||||
void operatorADEF1(const FineField &in, FineField &out) {
|
||||
|
||||
CoarseVector Csrc(_CoarseOperator.Grid());
|
||||
CoarseVector Ctmp(_CoarseOperator.Grid());
|
||||
CoarseVector Csol(_CoarseOperator.Grid());
|
||||
Csol = zero;
|
||||
|
||||
ConjugateGradient<CoarseVector> CG(1.0e-10, 100000);
|
||||
ConjugateGradient<FineField> fCG(3.0e-2, 1000);
|
||||
|
||||
HermitianLinearOperator<CoarseOperator, CoarseVector> HermOp(_CoarseOperator);
|
||||
MdagMLinearOperator<CoarseOperator, CoarseVector> MdagMOp(_CoarseOperator);
|
||||
ShiftedMdagMLinearOperator<Matrix, FineField> fMdagMOp(_FineMatrix, 0.1);
|
||||
|
||||
FineField tmp(in._grid);
|
||||
FineField res(in._grid);
|
||||
FineField Qin(in._grid);
|
||||
|
||||
// Monitor completeness of low mode space
|
||||
// _Aggregates.ProjectToSubspace (Csrc,in);
|
||||
// _Aggregates.PromoteFromSubspace(Csrc,out);
|
||||
// std::cout<<GridLogMessage<<"Coarse Grid Preconditioner\nCompleteness in: "<<std::sqrt(norm2(out)/norm2(in))<<std::endl;
|
||||
|
||||
_Aggregates.ProjectToSubspace(Csrc, in);
|
||||
HermOp.AdjOp(Csrc, Ctmp); // Normal equations
|
||||
CG(MdagMOp, Ctmp, Csol);
|
||||
_Aggregates.PromoteFromSubspace(Csol, Qin);
|
||||
|
||||
// Qin=0;
|
||||
_FineOperator.Op(Qin, tmp); // A Q in
|
||||
tmp = in - tmp; // in - A Q in
|
||||
|
||||
_FineOperator.Op(tmp, res); // this is the G5 herm bit
|
||||
fCG(fMdagMOp, res, out); // solves MdagM = g5 M g5M
|
||||
|
||||
out = out + Qin;
|
||||
|
||||
_FineOperator.Op(out, tmp);
|
||||
tmp = tmp - in; // tmp is new residual
|
||||
|
||||
std::cout << GridLogMessage << "preconditioner thinks residual is " << std::sqrt(norm2(tmp) / norm2(in)) << std::endl;
|
||||
}
|
||||
#endif
|
||||
|
||||
void SAP(const FineField &src, FineField &psi) {
|
||||
|
||||
Lattice<iScalar<vInteger>> coor(src._grid);
|
||||
Lattice<iScalar<vInteger>> subset(src._grid);
|
||||
|
||||
FineField r(src._grid);
|
||||
FineField zz(src._grid);
|
||||
zz = zero;
|
||||
FineField vec1(src._grid);
|
||||
FineField vec2(src._grid);
|
||||
|
||||
const Integer block = params.domainsize;
|
||||
|
||||
subset = zero;
|
||||
for(int mu = 0; mu < Nd; mu++) {
|
||||
LatticeCoordinate(coor, mu + 1);
|
||||
coor = div(coor, block);
|
||||
subset = subset + coor;
|
||||
}
|
||||
subset = mod(subset, (Integer)2);
|
||||
|
||||
ShiftedMdagMLinearOperator<Matrix, FineField> fMdagMOp(_SmootherMatrix, 0.0);
|
||||
Chebyshev<FineField> Cheby(params.lo, params.hi, params.order, InverseApproximation);
|
||||
|
||||
RealD resid;
|
||||
for(int i = 0; i < params.steps; i++) {
|
||||
|
||||
// Even domain residual
|
||||
_FineOperator.Op(psi, vec1); // this is the G5 herm bit
|
||||
r = src - vec1;
|
||||
resid = norm2(r) / norm2(src);
|
||||
std::cout << "SAP " << i << " resid " << resid << std::endl;
|
||||
|
||||
// Even domain solve
|
||||
r = where(subset == (Integer)0, r, zz);
|
||||
_SmootherOperator.AdjOp(r, vec1);
|
||||
Cheby(fMdagMOp, vec1, vec2); // solves MdagM = g5 M g5M
|
||||
psi = psi + vec2;
|
||||
|
||||
// Odd domain residual
|
||||
_FineOperator.Op(psi, vec1); // this is the G5 herm bit
|
||||
r = src - vec1;
|
||||
r = where(subset == (Integer)1, r, zz);
|
||||
|
||||
resid = norm2(r) / norm2(src);
|
||||
std::cout << "SAP " << i << " resid " << resid << std::endl;
|
||||
|
||||
// Odd domain solve
|
||||
_SmootherOperator.AdjOp(r, vec1);
|
||||
Cheby(fMdagMOp, vec1, vec2); // solves MdagM = g5 M g5M
|
||||
psi = psi + vec2;
|
||||
|
||||
_FineOperator.Op(psi, vec1); // this is the G5 herm bit
|
||||
r = src - vec1;
|
||||
resid = norm2(r) / norm2(src);
|
||||
std::cout << "SAP " << i << " resid " << resid << std::endl;
|
||||
}
|
||||
};
|
||||
|
||||
void SmootherTest(const FineField &in) {
|
||||
|
||||
FineField vec1(in._grid);
|
||||
FineField vec2(in._grid);
|
||||
|
||||
RealD lo[3] = {0.5, 1.0, 2.0};
|
||||
|
||||
// MdagMLinearOperator<Matrix,FineField> fMdagMOp(_FineMatrix);
|
||||
ShiftedMdagMLinearOperator<Matrix, FineField> fMdagMOp(_SmootherMatrix, 0.0);
|
||||
|
||||
RealD Ni, r;
|
||||
|
||||
Ni = norm2(in);
|
||||
|
||||
for(int ilo = 0; ilo < 3; ilo++) {
|
||||
for(int ord = 5; ord < 50; ord *= 2) {
|
||||
|
||||
_SmootherOperator.AdjOp(in, vec1);
|
||||
|
||||
Chebyshev<FineField> Cheby(lo[ilo], 70.0, ord, InverseApproximation);
|
||||
Cheby(fMdagMOp, vec1, vec2); // solves MdagM = g5 M g5M
|
||||
|
||||
_FineOperator.Op(vec2, vec1); // this is the G5 herm bit
|
||||
vec1 = in - vec1; // tmp = in - A Min
|
||||
r = norm2(vec1);
|
||||
std::cout << GridLogMessage << "Smoother resid " << std::sqrt(r / Ni) << std::endl;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
void operatorCheby(const FineField &in, FineField &out) {
|
||||
|
||||
CoarseVector Csrc(_CoarseOperator.Grid());
|
||||
CoarseVector Ctmp(_CoarseOperator.Grid());
|
||||
CoarseVector Csol(_CoarseOperator.Grid());
|
||||
Csol = zero;
|
||||
|
||||
ConjugateGradient<CoarseVector> CG(3.0e-3, 100000);
|
||||
// ConjugateGradient<FineField> fCG(3.0e-2,1000);
|
||||
|
||||
HermitianLinearOperator<CoarseOperator, CoarseVector> HermOp(_CoarseOperator);
|
||||
MdagMLinearOperator<CoarseOperator, CoarseVector> MdagMOp(_CoarseOperator);
|
||||
// MdagMLinearOperator<Matrix,FineField> fMdagMOp(_FineMatrix);
|
||||
ShiftedMdagMLinearOperator<Matrix, FineField> fMdagMOp(_SmootherMatrix, 0.0);
|
||||
|
||||
FineField vec1(in._grid);
|
||||
FineField vec2(in._grid);
|
||||
|
||||
// Chebyshev<FineField> Cheby (0.5,70.0,30,InverseApproximation);
|
||||
// Chebyshev<FineField> ChebyAccu(0.5,70.0,30,InverseApproximation);
|
||||
Chebyshev<FineField> Cheby(params.lo, params.hi, params.order, InverseApproximation);
|
||||
Chebyshev<FineField> ChebyAccu(params.lo, params.hi, params.order, InverseApproximation);
|
||||
// Cheby.JacksonSmooth();
|
||||
// ChebyAccu.JacksonSmooth();
|
||||
|
||||
// _Aggregates.ProjectToSubspace (Csrc,in);
|
||||
// _Aggregates.PromoteFromSubspace(Csrc,out);
|
||||
// std::cout<<GridLogMessage<<"Completeness: "<<std::sqrt(norm2(out)/norm2(in))<<std::endl;
|
||||
|
||||
// ofstream fout("smoother");
|
||||
// Cheby.csv(fout);
|
||||
|
||||
// V11 multigrid.
|
||||
// Use a fixed chebyshev and hope hermiticity helps.
|
||||
|
||||
// To make a working smoother for indefinite operator
|
||||
// must multiply by "Mdag" (ouch loses all low mode content)
|
||||
// and apply to poly approx of (mdagm)^-1.
|
||||
// so that we end up with an odd polynomial.
|
||||
FineField fineTmp1(in._grid);
|
||||
FineField fineTmp2(in._grid);
|
||||
|
||||
RealD Ni = norm2(in);
|
||||
|
||||
_SmootherOperator.AdjOp(in, vec1); // this is the G5 herm bit
|
||||
ChebyAccu(fMdagMOp, vec1, out); // solves MdagM = g5 M g5M
|
||||
// no pre smoothing for now
|
||||
auto preSmootherNorm = 0;
|
||||
auto preSmootherResidual = 0;
|
||||
RealD r;
|
||||
|
||||
std::cout << GridLogMessage << "Smoother norm " << norm2(out) << std::endl;
|
||||
// Project to coarse grid, solve, project back to fine grid
|
||||
_Aggregates.ProjectToSubspace(coarseSrc, in);
|
||||
coarseGMRES(coarseMdagMOp, coarseSrc, coarseSol);
|
||||
_Aggregates.PromoteFromSubspace(coarseSol, out);
|
||||
|
||||
// Update with residual for out
|
||||
_FineOperator.Op(out, vec1); // this is the G5 herm bit
|
||||
vec1 = in - vec1; // tmp = in - A Min
|
||||
|
||||
RealD r = norm2(vec1);
|
||||
|
||||
std::cout << GridLogMessage << "Smoother resid " << std::sqrt(r / Ni) << " " << r << " " << Ni << std::endl;
|
||||
|
||||
_Aggregates.ProjectToSubspace(Csrc, vec1);
|
||||
HermOp.AdjOp(Csrc, Ctmp); // Normal equations
|
||||
CG(MdagMOp, Ctmp, Csol);
|
||||
_Aggregates.PromoteFromSubspace(Csol, vec1); // Ass^{-1} [in - A Min]_s
|
||||
// Q = Q[in - A Min]
|
||||
out = out + vec1;
|
||||
|
||||
// Three preconditioner smoothing -- hermitian if C3 = C1
|
||||
// Recompute error
|
||||
_FineOperator.Op(out, vec1); // this is the G5 herm bit
|
||||
vec1 = in - vec1; // tmp = in - A Min
|
||||
r = norm2(vec1);
|
||||
_FineOperator.Op(out, fineTmp1);
|
||||
fineTmp1 = in - fineTmp1;
|
||||
r = norm2(fineTmp1);
|
||||
auto coarseResidual = std::sqrt(r / Ni);
|
||||
|
||||
std::cout << GridLogMessage << "Coarse resid " << std::sqrt(r / Ni) << std::endl;
|
||||
// Apply smoother, use GMRES for the moment
|
||||
fineGMRES(fineMdagMOp, in, out);
|
||||
|
||||
// Reapply smoother
|
||||
_SmootherOperator.Op(vec1, vec2); // this is the G5 herm bit
|
||||
ChebyAccu(fMdagMOp, vec2, vec1); // solves MdagM = g5 M g5M
|
||||
|
||||
out = out + vec1;
|
||||
vec1 = in - vec1; // tmp = in - A Min
|
||||
r = norm2(vec1);
|
||||
std::cout << GridLogMessage << "Smoother resid " << std::sqrt(r / Ni) << std::endl;
|
||||
}
|
||||
|
||||
void operatorSAP(const FineField &in, FineField &out) {
|
||||
|
||||
CoarseVector Csrc(_CoarseOperator.Grid());
|
||||
CoarseVector Ctmp(_CoarseOperator.Grid());
|
||||
CoarseVector Csol(_CoarseOperator.Grid());
|
||||
Csol = zero;
|
||||
|
||||
ConjugateGradient<CoarseVector> CG(1.0e-3, 100000);
|
||||
|
||||
HermitianLinearOperator<CoarseOperator, CoarseVector> HermOp(_CoarseOperator);
|
||||
MdagMLinearOperator<CoarseOperator, CoarseVector> MdagMOp(_CoarseOperator);
|
||||
|
||||
FineField vec1(in._grid);
|
||||
FineField vec2(in._grid);
|
||||
|
||||
_Aggregates.ProjectToSubspace(Csrc, in);
|
||||
_Aggregates.PromoteFromSubspace(Csrc, out);
|
||||
std::cout << GridLogMessage << "Completeness: " << std::sqrt(norm2(out) / norm2(in)) << std::endl;
|
||||
|
||||
// To make a working smoother for indefinite operator
|
||||
// must multiply by "Mdag" (ouch loses all low mode content)
|
||||
// and apply to poly approx of (mdagm)^-1.
|
||||
// so that we end up with an odd polynomial.
|
||||
SAP(in, out);
|
||||
|
||||
// Update with residual for out
|
||||
_FineOperator.Op(out, vec1); // this is the G5 herm bit
|
||||
vec1 = in - vec1; // tmp = in - A Min
|
||||
|
||||
RealD r = norm2(vec1);
|
||||
RealD Ni = norm2(in);
|
||||
std::cout << GridLogMessage << "SAP resid " << std::sqrt(r / Ni) << " " << r << " " << Ni << std::endl;
|
||||
|
||||
_Aggregates.ProjectToSubspace(Csrc, vec1);
|
||||
HermOp.AdjOp(Csrc, Ctmp); // Normal equations
|
||||
CG(MdagMOp, Ctmp, Csol);
|
||||
_Aggregates.PromoteFromSubspace(Csol, vec1); // Ass^{-1} [in - A Min]_s
|
||||
// Q = Q[in - A Min]
|
||||
out = out + vec1;
|
||||
|
||||
// Three preconditioner smoothing -- hermitian if C3 = C1
|
||||
// Recompute error
|
||||
_FineOperator.Op(out, vec1); // this is the G5 herm bit
|
||||
vec1 = in - vec1; // tmp = in - A Min
|
||||
r = norm2(vec1);
|
||||
_FineOperator.Op(out, fineTmp1);
|
||||
fineTmp1 = in - fineTmp1;
|
||||
r = norm2(fineTmp1);
|
||||
auto postSmootherResidual = std::sqrt(r / Ni);
|
||||
|
||||
std::cout << GridLogMessage << "Coarse resid " << std::sqrt(r / Ni) << std::endl;
|
||||
|
||||
// Reapply smoother
|
||||
SAP(vec1, vec2);
|
||||
out = out + vec2;
|
||||
|
||||
// Update with residual for out
|
||||
_FineOperator.Op(out, vec1); // this is the G5 herm bit
|
||||
vec1 = in - vec1; // tmp = in - A Min
|
||||
|
||||
r = norm2(vec1);
|
||||
Ni = norm2(in);
|
||||
std::cout << GridLogMessage << "SAP resid(post) " << std::sqrt(r / Ni) << " " << r << " " << Ni << std::endl;
|
||||
std::cout << GridLogIterative << "Input norm = " << Ni << " Pre-Smoother norm " << preSmootherNorm
|
||||
<< " Pre-Smoother residual = " << preSmootherResidual << " Coarse residual = " << coarseResidual
|
||||
<< " Post-Smoother residual = " << postSmootherResidual << std::endl;
|
||||
}
|
||||
|
||||
void runChecks(CoarseGrids<nbasis> &cGrids, int whichCoarseGrid) {
|
||||
@ -649,16 +323,16 @@ public:
|
||||
fTmps[1] = _Aggregates.subspace[i] - fTmps[0]; // v_i - P R v_i
|
||||
auto deviation = std::sqrt(norm2(fTmps[1]) / norm2(_Aggregates.subspace[i]));
|
||||
|
||||
std::cout << GridLogMessage << "Vector " << i << ": norm2(v_i) = " << norm2(_Aggregates.subspace[i])
|
||||
<< " | norm2(R v_i) = " << norm2(cTmps[0]) << " | norm2(P R v_i) = " << norm2(fTmps[0])
|
||||
<< " | relative deviation = " << deviation << std::endl;
|
||||
std::cout << GridLogMessage << "Vector " << i << ": norm2(v_i) = " << norm2(_Aggregates.subspace[i]) << " | norm2(R v_i) = " << norm2(cTmps[0])
|
||||
<< " | norm2(P R v_i) = " << norm2(fTmps[0]) << " | relative deviation = " << deviation;
|
||||
|
||||
if(deviation > tolerance) {
|
||||
std::cout << GridLogError << "Vector " << i << ": relative deviation check failed " << deviation << " > " << tolerance << std::endl;
|
||||
abort();
|
||||
std::cout << " > " << tolerance << " -> check failed" << std::endl;
|
||||
// abort();
|
||||
} else {
|
||||
std::cout << " < " << tolerance << " -> check passed" << std::endl;
|
||||
}
|
||||
}
|
||||
std::cout << GridLogMessage << "Check passed!" << std::endl;
|
||||
|
||||
std::cout << GridLogMessage << "**************************************************" << std::endl;
|
||||
std::cout << GridLogMessage << "MG correctness check: 0 == (1 - R P) v_c" << std::endl;
|
||||
@ -673,13 +347,14 @@ public:
|
||||
auto deviation = std::sqrt(norm2(cTmps[2]) / norm2(cTmps[0]));
|
||||
|
||||
std::cout << GridLogMessage << "norm2(v_c) = " << norm2(cTmps[0]) << " | norm2(R P v_c) = " << norm2(cTmps[1])
|
||||
<< " | norm2(P v_c) = " << norm2(fTmps[0]) << " | relative deviation = " << deviation << std::endl;
|
||||
<< " | norm2(P v_c) = " << norm2(fTmps[0]) << " | relative deviation = " << deviation;
|
||||
|
||||
if(deviation > tolerance) {
|
||||
std::cout << GridLogError << "relative deviation check failed " << deviation << " > " << tolerance << std::endl;
|
||||
abort();
|
||||
std::cout << " > " << tolerance << " -> check failed" << std::endl;
|
||||
// abort();
|
||||
} else {
|
||||
std::cout << " < " << tolerance << " -> check passed" << std::endl;
|
||||
}
|
||||
std::cout << GridLogMessage << "Check passed!" << std::endl;
|
||||
|
||||
std::cout << GridLogMessage << "**************************************************" << std::endl;
|
||||
std::cout << GridLogMessage << "MG correctness check: 0 == (R D P - D_c) v_c" << std::endl;
|
||||
@ -697,13 +372,14 @@ public:
|
||||
deviation = std::sqrt(norm2(cTmps[3]) / norm2(cTmps[1]));
|
||||
|
||||
std::cout << GridLogMessage << "norm2(R D P v_c) = " << norm2(cTmps[1]) << " | norm2(D_c v_c) = " << norm2(cTmps[2])
|
||||
<< " | relative deviation = " << deviation << std::endl;
|
||||
<< " | relative deviation = " << deviation;
|
||||
|
||||
if(deviation > tolerance) {
|
||||
std::cout << GridLogError << "relative deviation check failed " << deviation << " > " << tolerance << std::endl;
|
||||
abort();
|
||||
std::cout << " > " << tolerance << " -> check failed" << std::endl;
|
||||
// abort();
|
||||
} else {
|
||||
std::cout << " < " << tolerance << " -> check passed" << std::endl;
|
||||
}
|
||||
std::cout << GridLogMessage << "Check passed!" << std::endl;
|
||||
|
||||
std::cout << GridLogMessage << "**************************************************" << std::endl;
|
||||
std::cout << GridLogMessage << "MG correctness check: 0 == |(Im(v_c^dag D_c^dag D_c v_c)|" << std::endl;
|
||||
@ -714,20 +390,18 @@ public:
|
||||
MdagMOp.Op(cTmps[0], cTmps[1]); // D_c v_c
|
||||
MdagMOp.AdjOp(cTmps[1], cTmps[2]); // D_c^dag D_c v_c
|
||||
|
||||
// // alternative impl, which is better?
|
||||
// MdagMOp.HermOp(cTmps[0], cTmps[2]); // D_c^dag D_c v_c
|
||||
|
||||
auto dot = innerProduct(cTmps[0], cTmps[2]); //v_c^dag D_c^dag D_c v_c
|
||||
deviation = abs(imag(dot)) / abs(real(dot));
|
||||
|
||||
std::cout << GridLogMessage << "Re(v_c^dag D_c^dag D_c v_c) = " << real(dot) << " | Im(v_c^dag D_c^dag D_c v_c) = " << imag(dot)
|
||||
<< " | relative deviation = " << deviation << std::endl;
|
||||
<< " | relative deviation = " << deviation;
|
||||
|
||||
if(deviation > tolerance) {
|
||||
std::cout << GridLogError << "relative deviation check failed " << deviation << " > " << tolerance << std::endl;
|
||||
abort();
|
||||
std::cout << " > " << tolerance << " -> check failed" << std::endl;
|
||||
// abort();
|
||||
} else {
|
||||
std::cout << " < " << tolerance << " -> check passed" << std::endl;
|
||||
}
|
||||
std::cout << GridLogMessage << "Check passed!" << std::endl;
|
||||
}
|
||||
};
|
||||
|
||||
@ -740,11 +414,10 @@ int main(int argc, char **argv) {
|
||||
params.domaindecompose = 0;
|
||||
params.order = 30;
|
||||
params.Ls = 1;
|
||||
// params.mq = .13;
|
||||
params.mq = .5;
|
||||
params.lo = 0.5;
|
||||
params.hi = 70.0;
|
||||
params.steps = 1;
|
||||
params.mq = -0.5;
|
||||
params.lo = 0.5;
|
||||
params.hi = 70.0;
|
||||
params.steps = 1;
|
||||
|
||||
std::cout << GridLogMessage << "**************************************************" << std::endl;
|
||||
std::cout << GridLogMessage << "Params: " << std::endl;
|
||||
@ -756,7 +429,7 @@ int main(int argc, char **argv) {
|
||||
std::cout << GridLogMessage << "Set up some fine level stuff: " << std::endl;
|
||||
std::cout << GridLogMessage << "**************************************************" << std::endl;
|
||||
|
||||
GridCartesian *FGrid = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd, vComplex::Nsimd()), GridDefaultMpi());
|
||||
GridCartesian * FGrid = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd, vComplex::Nsimd()), GridDefaultMpi());
|
||||
GridRedBlackCartesian *FrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(FGrid);
|
||||
|
||||
std::vector<int> fSeeds({1, 2, 3, 4});
|
||||
@ -766,57 +439,31 @@ int main(int argc, char **argv) {
|
||||
Gamma g5(Gamma::Algebra::Gamma5);
|
||||
|
||||
// clang-format off
|
||||
LatticeFermion src(FGrid); gaussian(fPRNG, src); // src=src + g5 * src;
|
||||
LatticeFermion src(FGrid); gaussian(fPRNG, src);
|
||||
LatticeFermion result(FGrid); result = zero;
|
||||
LatticeFermion ref(FGrid); ref = zero;
|
||||
LatticeFermion tmp(FGrid);
|
||||
LatticeFermion err(FGrid);
|
||||
LatticeGaugeField Umu(FGrid); SU3::HotConfiguration(fPRNG, Umu);
|
||||
LatticeGaugeField UmuDD(FGrid);
|
||||
LatticeColourMatrix U(FGrid);
|
||||
LatticeColourMatrix zz(FGrid);
|
||||
// clang-format on
|
||||
|
||||
if(params.domaindecompose) {
|
||||
Lattice<iScalar<vInteger>> coor(FGrid);
|
||||
zz = zero;
|
||||
for(int mu = 0; mu < Nd; mu++) {
|
||||
LatticeCoordinate(coor, mu);
|
||||
U = PeekIndex<LorentzIndex>(Umu, mu);
|
||||
U = where(mod(coor, params.domainsize) == (Integer)0, zz, U);
|
||||
PokeIndex<LorentzIndex>(UmuDD, U, mu);
|
||||
}
|
||||
} else {
|
||||
UmuDD = Umu;
|
||||
}
|
||||
|
||||
RealD mass = params.mq;
|
||||
|
||||
std::cout << GridLogMessage << "**************************************************" << std::endl;
|
||||
std::cout << GridLogMessage << "Set up some coarser levels stuff: " << std::endl;
|
||||
std::cout << GridLogMessage << "**************************************************" << std::endl;
|
||||
|
||||
const int nbasis = 20; // we fix the number of test vector to the same
|
||||
const int nbasis = 20; // fix the number of test vector to the same
|
||||
// number on every level for now
|
||||
|
||||
//////////////////////////////////////////
|
||||
// toggle to run two/three level method
|
||||
//////////////////////////////////////////
|
||||
|
||||
// // two-level algorithm
|
||||
// std::vector<std::vector<int>> blockSizes({{2, 2, 2, 2}});
|
||||
// CoarseGrids<nbasis> coarseGrids(blockSizes, 1);
|
||||
// two-level algorithm
|
||||
std::vector<std::vector<int>> blockSizes({{2, 2, 2, 2}});
|
||||
CoarseGrids<nbasis> coarseGrids(blockSizes, 1);
|
||||
|
||||
// three-level algorithm
|
||||
std::vector<std::vector<int>> blockSizes({{2, 2, 2, 2}, {2, 2, 1, 1}});
|
||||
CoarseGrids<nbasis> coarseGrids(blockSizes, 2);
|
||||
|
||||
std::cout << GridLogMessage << "**************************************************" << std::endl;
|
||||
std::cout << GridLogMessage << "Building the wilson operator on the fine grid" << std::endl;
|
||||
std::cout << GridLogMessage << "**************************************************" << std::endl;
|
||||
|
||||
WilsonFermionR Dw(Umu, *FGrid, *FrbGrid, mass);
|
||||
WilsonFermionR DwDD(UmuDD, *FGrid, *FrbGrid, mass);
|
||||
// // three-level algorithm
|
||||
// std::vector<std::vector<int>> blockSizes({{2, 2, 2, 2}, {2, 2, 1, 1}});
|
||||
// CoarseGrids<nbasis> coarseGrids(blockSizes, 2);
|
||||
|
||||
std::cout << GridLogMessage << "**************************************************" << std::endl;
|
||||
std::cout << GridLogMessage << "Some typedefs" << std::endl;
|
||||
@ -844,22 +491,40 @@ int main(int argc, char **argv) {
|
||||
|
||||
static_assert(std::is_same<CoarseVector, CoarseCoarseVector>::value, "CoarseVector and CoarseCoarseVector must be of the same type");
|
||||
|
||||
std::cout << GridLogMessage << "**************************************************" << std::endl;
|
||||
std::cout << GridLogMessage << "Building the wilson operator on the fine grid" << std::endl;
|
||||
std::cout << GridLogMessage << "**************************************************" << std::endl;
|
||||
|
||||
WilsonFermionR Dw(Umu, *FGrid, *FrbGrid, mass);
|
||||
|
||||
std::cout << GridLogMessage << "**************************************************" << std::endl;
|
||||
std::cout << GridLogMessage << "Setting up linear operators" << std::endl;
|
||||
std::cout << GridLogMessage << "**************************************************" << std::endl;
|
||||
|
||||
MdagMLinearOperator<WilsonFermionR, LatticeFermion> FineMdagMOp(Dw);
|
||||
|
||||
std::cout << GridLogMessage << "**************************************************" << std::endl;
|
||||
std::cout << GridLogMessage << "Calling Aggregation class to build subspaces" << std::endl;
|
||||
std::cout << GridLogMessage << "**************************************************" << std::endl;
|
||||
|
||||
MdagMLinearOperator<WilsonFermionR, LatticeFermion> FineHermPosdefOp(Dw);
|
||||
Subspace FineAggregates(coarseGrids.Grids[0], FGrid, 0);
|
||||
Subspace FineAggregates(coarseGrids.Grids[0], FGrid, 0);
|
||||
|
||||
assert((nbasis & 0x1) == 0);
|
||||
int nb = nbasis / 2;
|
||||
std::cout << GridLogMessage << " nbasis/2 = " << nb << std::endl;
|
||||
|
||||
FineAggregates.CreateSubspace(fPRNG, FineHermPosdefOp /*, nb */); // Don't specify nb to see the orthogonalization check
|
||||
FineAggregates.CreateSubspace(fPRNG, FineMdagMOp /*, nb */); // Don't specify nb to see the orthogonalization check
|
||||
|
||||
std::cout << GridLogMessage << "**************************************************" << std::endl;
|
||||
std::cout << GridLogMessage << "Test vector analysis after initial creation of subspace" << std::endl;
|
||||
std::cout << GridLogMessage << "**************************************************" << std::endl;
|
||||
|
||||
std::cout << GridLogMessage << "Test vector analysis after initial creation of MG test vectors" << std::endl;
|
||||
FineTVA fineTVA;
|
||||
fineTVA(FineHermPosdefOp, FineAggregates.subspace, nb);
|
||||
fineTVA(FineMdagMOp, FineAggregates.subspace);
|
||||
|
||||
std::cout << GridLogMessage << "**************************************************" << std::endl;
|
||||
std::cout << GridLogMessage << "Projecting subspace to definite chirality" << std::endl;
|
||||
std::cout << GridLogMessage << "**************************************************" << std::endl;
|
||||
|
||||
for(int n = 0; n < nb; n++) {
|
||||
FineAggregates.subspace[n + nb] = g5 * FineAggregates.subspace[n];
|
||||
@ -877,14 +542,26 @@ int main(int argc, char **argv) {
|
||||
std::cout << GridLogMessage << "Building coarse representation of Dirac operator" << std::endl;
|
||||
std::cout << GridLogMessage << "**************************************************" << std::endl;
|
||||
|
||||
// using Gamma5HermitianLinearOperator corresponds to working with H = g5 * D
|
||||
Gamma5HermitianLinearOperator<WilsonFermionR, LatticeFermion> FineHermIndefOp(Dw);
|
||||
Gamma5HermitianLinearOperator<WilsonFermionR, LatticeFermion> FineHermIndefOpDD(DwDD);
|
||||
CoarseOperator Dc(*coarseGrids.Grids[0]);
|
||||
Dc.CoarsenOperator(FGrid, FineHermIndefOp, FineAggregates); // uses only linop.OpDiag & linop.OpDir
|
||||
CoarseOperator Dc(*coarseGrids.Grids[0]);
|
||||
|
||||
std::cout << GridLogMessage << "Test vector analysis after construction of D_c" << std::endl;
|
||||
fineTVA(FineHermPosdefOp, FineAggregates.subspace, nb);
|
||||
Dc.CoarsenOperator(FGrid, FineMdagMOp, FineAggregates);
|
||||
|
||||
MdagMLinearOperator<CoarseOperator, CoarseVector> CoarseMdagMOp(Dc);
|
||||
|
||||
std::cout << GridLogMessage << "**************************************************" << std::endl;
|
||||
std::cout << GridLogMessage << "Test vector analysis after construction of coarse Dirac operator" << std::endl;
|
||||
std::cout << GridLogMessage << "**************************************************" << std::endl;
|
||||
|
||||
fineTVA(FineMdagMOp, FineAggregates.subspace);
|
||||
|
||||
std::cout << GridLogMessage << "**************************************************" << std::endl;
|
||||
std::cout << GridLogMessage << "Testing the linear operators" << std::endl;
|
||||
std::cout << GridLogMessage << "**************************************************" << std::endl;
|
||||
|
||||
// clang-format off
|
||||
testLinearOperator(FineMdagMOp, FGrid, "FineMdagMOp"); std::cout << GridLogMessage << std::endl;
|
||||
testLinearOperator(CoarseMdagMOp, coarseGrids.Grids[0], "CoarseMdagMOp"); std::cout << GridLogMessage << std::endl;
|
||||
// clang-format on
|
||||
|
||||
std::cout << GridLogMessage << "**************************************************" << std::endl;
|
||||
std::cout << GridLogMessage << "Building coarse vectors" << std::endl;
|
||||
@ -895,64 +572,85 @@ int main(int argc, char **argv) {
|
||||
gaussian(coarseGrids.PRNGs[0], coarseSource);
|
||||
coarseResult = zero;
|
||||
|
||||
std::cout << GridLogMessage << "**************************************************" << std::endl;
|
||||
std::cout << GridLogMessage << "Building some coarse space solvers" << std::endl;
|
||||
std::cout << GridLogMessage << "**************************************************" << std::endl;
|
||||
|
||||
std::vector<std::unique_ptr<OperatorFunction<CoarseVector>>> dummyCoarseSolvers;
|
||||
dummyCoarseSolvers.emplace_back(new GeneralisedMinimalResidual<CoarseVector>(5.0e-2, 100, 8, false));
|
||||
dummyCoarseSolvers.emplace_back(new MinimalResidual<CoarseVector>(5.0e-2, 100, 0.8, false));
|
||||
dummyCoarseSolvers.emplace_back(new ConjugateGradient<CoarseVector>(5.0e-2, 100, false));
|
||||
|
||||
std::cout << GridLogMessage << "**************************************************" << std::endl;
|
||||
std::cout << GridLogMessage << "Testing some coarse space solvers" << std::endl;
|
||||
std::cout << GridLogMessage << "**************************************************" << std::endl;
|
||||
|
||||
MdagMLinearOperator<CoarseOperator, CoarseVector> CoarseHermPosdefOp(Dc);
|
||||
std::cout << GridLogMessage << "checking norm of coarse src " << norm2(coarseSource) << std::endl;
|
||||
|
||||
std::vector<std::unique_ptr<OperatorFunction<CoarseVector>>> coarseSolvers;
|
||||
coarseSolvers.emplace_back(new GeneralisedMinimalResidual<CoarseVector>(5.0e-2, 100, 8, false));
|
||||
coarseSolvers.emplace_back(new MinimalResidual<CoarseVector>(5.0e-2, 100, false));
|
||||
coarseSolvers.emplace_back(new ConjugateGradient<CoarseVector>(5.0e-2, 100, false));
|
||||
|
||||
for(auto const &solver : coarseSolvers) {
|
||||
for(auto const &solver : dummyCoarseSolvers) {
|
||||
coarseResult = zero;
|
||||
(*solver)(CoarseHermPosdefOp, coarseSource, coarseResult);
|
||||
(*solver)(CoarseMdagMOp, coarseSource, coarseResult);
|
||||
}
|
||||
|
||||
std::cout << GridLogMessage << "**************************************************" << std::endl;
|
||||
std::cout << GridLogMessage << "Testing the operators" << std::endl;
|
||||
std::cout << GridLogMessage << "Building a multigrid preconditioner" << std::endl;
|
||||
std::cout << GridLogMessage << "**************************************************" << std::endl;
|
||||
|
||||
std::cout << GridLogMessage << "MdagMLinearOperator<WilsonFermionR, LatticeFermion> FineHermPosdefOp(Dw);" << std::endl;
|
||||
testOperator(FineHermPosdefOp, FGrid);
|
||||
std::cout << GridLogMessage << "Gamma5HermitianLinearOperator<WilsonFermionR, LatticeFermion> FineHermIndefOp(Dw);" << std::endl;
|
||||
testOperator(FineHermIndefOp, FGrid);
|
||||
std::cout << GridLogMessage << "Gamma5HermitianLinearOperator<WilsonFermionR, LatticeFermion> FineHermIndefOpDD(DwDD);" << std::endl;
|
||||
testOperator(FineHermIndefOpDD, FGrid);
|
||||
std::cout << GridLogMessage << "MdagMLinearOperator<CoarseOperator, CoarseVector> CoarseHermPosdefOp(Dc);" << std::endl;
|
||||
testOperator(CoarseHermPosdefOp, coarseGrids.Grids[0]);
|
||||
|
||||
std::cout << GridLogMessage << "**************************************************" << std::endl;
|
||||
std::cout << GridLogMessage << "Building deflation preconditioner " << std::endl;
|
||||
std::cout << GridLogMessage << "**************************************************" << std::endl;
|
||||
|
||||
FineMGPreconditioner FineMGPrecon(FineAggregates, Dc, FineHermIndefOp, Dw, FineHermIndefOp, Dw);
|
||||
|
||||
FineMGPreconditioner FineMGPreconDD(FineAggregates, Dc, FineHermIndefOp, Dw, FineHermIndefOpDD, DwDD);
|
||||
|
||||
FineMGPreconditioner FineMGPrecon(FineAggregates, Dc, FineMdagMOp, Dw, FineMdagMOp, Dw);
|
||||
FineTrivialPreconditioner FineSimplePrecon;
|
||||
|
||||
FineMGPrecon.runChecks(coarseGrids, 0);
|
||||
|
||||
std::cout << GridLogMessage << "**************************************************" << std::endl;
|
||||
std::cout << GridLogMessage << "Building krylov subspace solvers w/ & w/o MG Preconditioner" << std::endl;
|
||||
std::cout << GridLogMessage << "**************************************************" << std::endl;
|
||||
|
||||
std::vector<std::unique_ptr<OperatorFunction<LatticeFermion>>> solvers;
|
||||
solvers.emplace_back(new FlexibleGeneralisedMinimalResidual<LatticeFermion>(1.0e-12, 4000000, FineSimplePrecon, 25, false));
|
||||
solvers.emplace_back(new FlexibleGeneralisedMinimalResidual<LatticeFermion>(1.0e-12, 100, FineMGPrecon, 25, false));
|
||||
solvers.emplace_back(new PrecGeneralisedConjugateResidual<LatticeFermion>(1.0e-12, 4000000, FineSimplePrecon, 25, 25));
|
||||
|
||||
std::cout << GridLogMessage << "**************************************************" << std::endl;
|
||||
std::cout << GridLogMessage << "Testing the (un)?preconditioned solvers" << std::endl;
|
||||
std::cout << GridLogMessage << "**************************************************" << std::endl;
|
||||
|
||||
for(auto const &solver : solvers) {
|
||||
std::cout << GridLogMessage << "checking norm of fine src " << norm2(src) << std::endl;
|
||||
result = zero;
|
||||
(*solver)(FineMdagMOp, src, result);
|
||||
std::cout << std::endl;
|
||||
}
|
||||
|
||||
#if 0
|
||||
if(coarseGrids.LattSizes.size() == 2) {
|
||||
|
||||
std::cout << GridLogMessage << "**************************************************" << std::endl;
|
||||
std::cout << GridLogMessage << "Dummy testing for building a second coarse level" << std::endl;
|
||||
std::cout << GridLogMessage << "Some testing for construction of a second coarse level" << std::endl;
|
||||
std::cout << GridLogMessage << "**************************************************" << std::endl;
|
||||
|
||||
// std::cout << GridLogMessage << "**************************************************" << std::endl;
|
||||
// std::cout << GridLogMessage << "Calling Aggregation class to build subspaces" << std::endl;
|
||||
// std::cout << GridLogMessage << "**************************************************" << std::endl;
|
||||
|
||||
SubSubSpace CoarseAggregates(coarseGrids.Grids[1], coarseGrids.Grids[0], 0);
|
||||
CoarseAggregates.CreateSubspace(coarseGrids.PRNGs[0], CoarseHermPosdefOp);
|
||||
CoarseAggregates.CreateSubspace(coarseGrids.PRNGs[0], CoarseMdagMOp);
|
||||
|
||||
// std::cout << GridLogMessage << "**************************************************" << std::endl;
|
||||
// std::cout << GridLogMessage << "Test vector analysis after initial creation of subspace" << std::endl;
|
||||
// std::cout << GridLogMessage << "**************************************************" << std::endl;
|
||||
|
||||
// // this doesn't work because this function applies g5 to a vector, which
|
||||
// // doesn't work for coarse vectors atm -> FIXME
|
||||
// CoarseTVA coarseTVA;
|
||||
// coarseTVA(CoarseHermPosdefOp, CoarseAggregates.subspace, nb);
|
||||
// coarseTVA(CoarseMdagMOp, CoarseAggregates.subspace);
|
||||
|
||||
// std::cout << GridLogMessage << "**************************************************" << std::endl;
|
||||
// std::cout << GridLogMessage << "Projecting subspace to definite chirality" << std::endl;
|
||||
// std::cout << GridLogMessage << "**************************************************" << std::endl;
|
||||
|
||||
// // cannot apply g5 to coarse vectors atm -> FIXME
|
||||
// for(int n=0;n<nb;n++){
|
||||
// CoarseAggregates.subspace[n+nb] = g5 * CoarseAggregates.subspace[n]; // multiply with g5 normally instead of G5R5 since this specific to DWF
|
||||
// CoarseAggregates.subspace[n+nb] = g5 * CoarseAggregates.subspace[n];
|
||||
// std::cout<<GridLogMessage<<n<<" subspace "<<norm2(CoarseAggregates.subspace[n+nb])<<" "<<norm2(CoarseAggregates.subspace[n]) <<std::endl;
|
||||
// }
|
||||
|
||||
@ -965,59 +663,93 @@ int main(int argc, char **argv) {
|
||||
std::cout << GridLogMessage << "vec[" << n << "] = " << norm2(CoarseAggregates.subspace[n]) << std::endl;
|
||||
}
|
||||
|
||||
// std::cout << GridLogMessage << "**************************************************" << std::endl;
|
||||
// std::cout << GridLogMessage << "Building coarse coarse representation of Dirac operator" << std::endl;
|
||||
// std::cout << GridLogMessage << "**************************************************" << std::endl;
|
||||
|
||||
CoarseCoarseOperator Dcc(*coarseGrids.Grids[1]);
|
||||
Dcc.CoarsenOperator(coarseGrids.Grids[0], CoarseHermPosdefOp, CoarseAggregates); // uses only linop.OpDiag & linop.OpDir
|
||||
|
||||
Dcc.CoarsenOperator(coarseGrids.Grids[0], CoarseMdagMOp, CoarseAggregates);
|
||||
|
||||
MdagMLinearOperator<CoarseCoarseOperator, CoarseCoarseVector> CoarseCoarseMdagMOp(Dcc);
|
||||
|
||||
// std::cout << GridLogMessage << "**************************************************" << std::endl;
|
||||
// std::cout << GridLogMessage << "Test vector analysis after construction of coarse Dirac operator" << std::endl;
|
||||
// std::cout << GridLogMessage << "**************************************************" << std::endl;
|
||||
|
||||
// // this doesn't work because this function applies g5 to a vector, which
|
||||
// // doesn't work for coarse vectors atm -> FIXME
|
||||
// std::cout << GridLogMessage << "Test vector analysis after construction of D_c_c" << std::endl;
|
||||
// coarseTVA(CoarseHermPosdefOp, CoarseAggregates.subspace, nb);
|
||||
// coarseTVA(CoarseMdagMOp, CoarseAggregates.subspace);
|
||||
|
||||
// std::cout << GridLogMessage << "**************************************************" << std::endl;
|
||||
// std::cout << GridLogMessage << "Testing the linear operators" << std::endl;
|
||||
// std::cout << GridLogMessage << "**************************************************" << std::endl;
|
||||
|
||||
// clang-format off
|
||||
testLinearOperator(CoarseMdagMOp, coarseGrids.Grids[0], "CoarseMdagMOp");
|
||||
testLinearOperator(CoarseCoarseMdagMOp, coarseGrids.Grids[1], "CoarseCoarseMdagMOp");
|
||||
// clang-format on
|
||||
|
||||
// std::cout << GridLogMessage << "**************************************************" << std::endl;
|
||||
// std::cout << GridLogMessage << "Building coarse coarse vectors" << std::endl;
|
||||
// std::cout << GridLogMessage << "**************************************************" << std::endl;
|
||||
|
||||
CoarseCoarseVector coarseCoarseSource(coarseGrids.Grids[1]);
|
||||
CoarseCoarseVector coarseCoarseResult(coarseGrids.Grids[1]);
|
||||
gaussian(coarseGrids.PRNGs[1], coarseCoarseSource);
|
||||
coarseCoarseResult = zero;
|
||||
|
||||
MdagMLinearOperator<CoarseCoarseOperator, CoarseCoarseVector> CoarseCoarseHermPosdefOp(Dcc);
|
||||
// std::cout << GridLogMessage << "**************************************************" << std::endl;
|
||||
// std::cout << GridLogMessage << "Building some coarse space solvers" << std::endl;
|
||||
// std::cout << GridLogMessage << "**************************************************" << std::endl;
|
||||
|
||||
std::vector<std::unique_ptr<OperatorFunction<CoarseCoarseVector>>> coarseCoarseSolvers;
|
||||
coarseSolvers.emplace_back(new GeneralisedMinimalResidual<CoarseCoarseVector>(5.0e-2, 100, 8, false));
|
||||
coarseSolvers.emplace_back(new MinimalResidual<CoarseCoarseVector>(5.0e-2, 100, false));
|
||||
coarseSolvers.emplace_back(new ConjugateGradient<CoarseCoarseVector>(5.0e-2, 100, false));
|
||||
std::vector<std::unique_ptr<OperatorFunction<CoarseCoarseVector>>> dummyCoarseCoarseSolvers;
|
||||
dummyCoarseCoarseSolvers.emplace_back(new GeneralisedMinimalResidual<CoarseCoarseVector>(5.0e-2, 100, 8, false));
|
||||
dummyCoarseCoarseSolvers.emplace_back(new MinimalResidual<CoarseCoarseVector>(5.0e-2, 100, false));
|
||||
dummyCoarseCoarseSolvers.emplace_back(new ConjugateGradient<CoarseCoarseVector>(5.0e-2, 100, false));
|
||||
|
||||
for(auto const &solver : coarseCoarseSolvers) {
|
||||
// std::cout << GridLogMessage << "**************************************************" << std::endl;
|
||||
// std::cout << GridLogMessage << "Testing some coarse coarse space solvers" << std::endl;
|
||||
// std::cout << GridLogMessage << "**************************************************" << std::endl;
|
||||
|
||||
std::cout << GridLogMessage << "checking norm of coarse coarse src " << norm2(coarseCoarseSource) << std::endl;
|
||||
|
||||
for(auto const &solver : dummyCoarseCoarseSolvers) {
|
||||
coarseCoarseResult = zero;
|
||||
(*solver)(CoarseCoarseHermPosdefOp, coarseCoarseSource, coarseCoarseResult);
|
||||
(*solver)(CoarseCoarseMdagMOp, coarseCoarseSource, coarseCoarseResult);
|
||||
}
|
||||
|
||||
CoarseMGPreconditioner CoarseMGPrecon(CoarseAggregates, Dcc, CoarseHermPosdefOp, Dc, CoarseHermPosdefOp, Dc);
|
||||
// std::cout << GridLogMessage << "**************************************************" << std::endl;
|
||||
// std::cout << GridLogMessage << "Building a multigrid preconditioner" << std::endl;
|
||||
// std::cout << GridLogMessage << "**************************************************" << std::endl;
|
||||
|
||||
CoarseMGPreconditioner CoarseMGPrecon(CoarseAggregates, Dcc, CoarseMdagMOp, Dc, CoarseMdagMOp, Dc);
|
||||
CoarseTrivialPreconditioner CoarseSimplePrecon;
|
||||
|
||||
CoarseMGPrecon.runChecks(coarseGrids, 1);
|
||||
|
||||
std::cout << GridLogMessage << "ARTIFICIAL ABORT" << std::endl;
|
||||
abort();
|
||||
}
|
||||
|
||||
std::cout << GridLogMessage << "**************************************************" << std::endl;
|
||||
std::cout << GridLogMessage << "Building VPGCR and FGMRES solvers w/ & w/o MG Preconditioner" << std::endl;
|
||||
std::cout << GridLogMessage << "**************************************************" << std::endl;
|
||||
|
||||
std::vector<std::unique_ptr<OperatorFunction<LatticeFermion>>> solvers;
|
||||
solvers.emplace_back(new PrecGeneralisedConjugateResidual<LatticeFermion>(1.0e-12, 100, FineMGPrecon, 8, 8));
|
||||
solvers.emplace_back(new FlexibleGeneralisedMinimalResidual<LatticeFermion>(1.0e-12, 100, FineMGPrecon, 8));
|
||||
solvers.emplace_back(new PrecGeneralisedConjugateResidual<LatticeFermion>(1.0e-12, 4000000, FineSimplePrecon, 8, 8));
|
||||
solvers.emplace_back(new FlexibleGeneralisedMinimalResidual<LatticeFermion>(1.0e-12, 4000000, FineSimplePrecon, 8));
|
||||
|
||||
std::cout << GridLogMessage << "**************************************************" << std::endl;
|
||||
std::cout << GridLogMessage << "Testing the solvers" << std::endl;
|
||||
std::cout << GridLogMessage << "**************************************************" << std::endl;
|
||||
|
||||
std::cout << GridLogMessage << "checking norm src " << norm2(src) << std::endl;
|
||||
|
||||
for(auto const &solver : solvers) {
|
||||
result = zero;
|
||||
(*solver)(FineHermIndefOp, src, result);
|
||||
// std::cout << GridLogMessage << "**************************************************" << std::endl;
|
||||
// std::cout << GridLogMessage << "Building krylov subspace solvers w/ & w/o MG Preconditioner" << std::endl;
|
||||
// std::cout << GridLogMessage << "**************************************************" << std::endl;
|
||||
|
||||
std::vector<std::unique_ptr<OperatorFunction<CoarseVector>>> solvers;
|
||||
solvers.emplace_back(new FlexibleGeneralisedMinimalResidual<CoarseVector>(1.0e-12, 4000000, CoarseSimplePrecon, 25, false));
|
||||
solvers.emplace_back(new FlexibleGeneralisedMinimalResidual<CoarseVector>(1.0e-12, 100, CoarseMGPrecon, 25, false));
|
||||
solvers.emplace_back(new PrecGeneralisedConjugateResidual<CoarseVector>(1.0e-12, 4000000, CoarseSimplePrecon, 25, 25));
|
||||
|
||||
// std::cout << GridLogMessage << "**************************************************" << std::endl;
|
||||
// std::cout << GridLogMessage << "Testing the (un)?preconditioned solvers" << std::endl;
|
||||
// std::cout << GridLogMessage << "**************************************************" << std::endl;
|
||||
|
||||
for(auto const &solver : solvers) {
|
||||
std::cout << GridLogMessage << "checking norm of fine src " << norm2(coarseSource) << std::endl;
|
||||
coarseResult = zero;
|
||||
(*solver)(CoarseMdagMOp, coarseSource, coarseResult);
|
||||
std::cout << std::endl;
|
||||
}
|
||||
|
||||
}
|
||||
#endif
|
||||
|
||||
Grid_finalize();
|
||||
}
|
||||
|
||||
Loading…
Reference in New Issue
Block a user