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

Apply clang-format to Wilson MG

Browse Source 
I can provide the configuration file I used if people want that.
This commit is contained in:
Daniel Richtmann 2018-01-18 12:38:28 +01:00
parent fa4eeb28c4
commit 9732519c41
No known key found for this signature in database
GPG Key ID: B33C490AF0772057
1 changed files with 360 additions and 368 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

600
tests/solver/Test_wilson_ddalphaamg.cc
View File

@ -1,4 +1,4 @@
/************************************************************************************* /*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid Grid physics library, www.github.com/paboyle/Grid
@ -6,7 +6,7 @@
Copyright (C) 2015 Copyright (C) 2015
Author: Daniel Richtmann <daniel.richtmann@ur.de> Author: Daniel Richtmann <daniel.richtmann@ur.de>
This program is free software; you can redistribute it and/or modify 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 it under the terms of the GNU General Public License as published by
@ -24,7 +24,8 @@ Author: Daniel Richtmann <daniel.richtmann@ur.de>
See the full license in the file "LICENSE" in the top level distribution directory See the full license in the file "LICENSE" in the top level distribution directory
*************************************************************************************/ *************************************************************************************/
/* END LEGAL */ /* END LEGAL */
#include <Grid/Grid.h> #include <Grid/Grid.h>
#include <Grid/algorithms/iterative/PrecGeneralisedConjugateResidual.h> #include <Grid/algorithms/iterative/PrecGeneralisedConjugateResidual.h>
//#include <algorithms/iterative/PrecConjugateResidual.h> //#include <algorithms/iterative/PrecConjugateResidual.h>
@ -33,11 +34,9 @@ using namespace std;
using namespace Grid; using namespace Grid;
using namespace Grid::QCD; using namespace Grid::QCD;
template<class Field, int nbasis> template<class Field, int nbasis> class TestVectorAnalyzer {
class TestVectorAnalyzer {
public: public:
void operator()(LinearOperatorBase<Field> &Linop, std::vector<Field> const & vectors, int nn=nbasis) void operator()(LinearOperatorBase<Field> &Linop, std::vector<Field> const &vectors, int nn = nbasis) {
{
// this function corresponds to testvector_analysis_PRECISION from the // this function corresponds to testvector_analysis_PRECISION from the
// DD-αAMG codebase // DD-αAMG codebase
@ -48,7 +47,7 @@ public:
std::cout << GridLogMessage << "Test vector analysis:" << std::endl; std::cout << GridLogMessage << "Test vector analysis:" << std::endl;
for (auto i = 0; i < nn; ++i) { for(auto i = 0; i < nn; ++i) {
Linop.Op(vectors[i], tmp[3]); Linop.Op(vectors[i], tmp[3]);
@ -66,16 +65,16 @@ public:
positiveOnes++; positiveOnes++;
std::cout << GridLogMessage << std::scientific << std::setprecision(2) << std::setw(2) << std::showpos << "vector " << i << ": " std::cout << GridLogMessage << std::scientific << std::setprecision(2) << std::setw(2) << std::showpos << "vector " << i << ": "
<< "singular value: " << lambda << "singular value: " << lambda << ", singular vector precision: " << mu << ", norm: " << nrm << std::endl;
<< ", 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;
} }
}; };
class myclass: Serializable { class myclass : Serializable {
public: public:
// clang-format off
GRID_SERIALIZABLE_CLASS_MEMBERS(myclass, GRID_SERIALIZABLE_CLASS_MEMBERS(myclass,
int, domaindecompose, int, domaindecompose,
int, domainsize, int, domainsize,
@ -86,19 +85,16 @@ public:
double, lo, double, lo,
double, hi, double, hi,
int, steps); int, steps);
// clang-format on
myclass(){}; myclass(){};
}; };
myclass params; myclass params;
RealD InverseApproximation(RealD x){ RealD InverseApproximation(RealD x) {
return 1.0/x; return 1.0 / x;
} }
template <int nbasis> template<int nbasis> struct CoarseGrids {
struct CoarseGrids
{
public: public:
// typedef Aggregation<vSpinColourVector,vTComplex,nbasis> Subspace; // typedef Aggregation<vSpinColourVector,vTComplex,nbasis> Subspace;
// typedef CoarsenedMatrix<vSpinColourVector,vTComplex,nbasis> // typedef CoarsenedMatrix<vSpinColourVector,vTComplex,nbasis>
@ -110,19 +106,20 @@ public:
std::vector<GridCartesian *> Grids; std::vector<GridCartesian *> Grids;
std::vector<GridParallelRNG> PRNGs; std::vector<GridParallelRNG> PRNGs;
CoarseGrids(std::vector<std::vector<int>> const &blockSizes,int coarsegrids = 1) CoarseGrids(std::vector<std::vector<int>> const &blockSizes, int coarsegrids = 1) {
{
assert( blockSizes.size() == coarsegrids ); assert(blockSizes.size() == coarsegrids);
std::cout << GridLogMessage << "Constructing " << coarsegrids << " CoarseGrids" << std::endl; std::cout << GridLogMessage << "Constructing " << coarsegrids << " CoarseGrids" << std::endl;
for(int cl=0; cl<coarsegrids; ++cl) { // may be a bit ugly and slow but not perf critical for(int cl = 0; cl < coarsegrids; ++cl) { // may be a bit ugly and slow but not perf critical
LattSizes.push_back({GridDefaultLatt()}); LattSizes.push_back({GridDefaultLatt()});
Seeds.push_back(std::vector<int>(LattSizes[cl].size())); Seeds.push_back(std::vector<int>(LattSizes[cl].size()));
for(int d=0; d<LattSizes[cl].size(); ++d) { for(int d = 0; d < LattSizes[cl].size(); ++d) {
LattSizes[cl][d] = LattSizes[cl][d] / blockSizes[cl][d]; LattSizes[cl][d] = LattSizes[cl][d] / blockSizes[cl][d];
Seeds[cl][d] = (cl + 1) * LattSizes[cl].size() + d + 1; // unimportant, just to get. e.g., {5, // 6, 7, 8} for first coarse level and // so on Seeds[cl][d] = (cl + 1) * LattSizes[cl].size() + d + 1;
// calculation unimportant, just to get. e.g., {5, 6, 7, 8} for first coarse level and so on
} }
Grids.push_back(SpaceTimeGrid::makeFourDimGrid(LattSizes[cl], GridDefaultSimd(Nd, vComplex::Nsimd()), GridDefaultMpi())); Grids.push_back(SpaceTimeGrid::makeFourDimGrid(LattSizes[cl], GridDefaultSimd(Nd, vComplex::Nsimd()), GridDefaultMpi()));
@ -138,150 +135,148 @@ public:
// template < class Fobj, class CComplex, int coarseSpins, int nbasis, class Matrix > // template < class Fobj, class CComplex, int coarseSpins, int nbasis, class Matrix >
// class MultiGridPreconditioner : public LinearFunction< Lattice< Fobj > > { // class MultiGridPreconditioner : public LinearFunction< Lattice< Fobj > > {
template<class Fobj,class CComplex,int nbasis, class Matrix> template<class Fobj, class CComplex, int nbasis, class Matrix> class MultiGridPreconditioner : public LinearFunction<Lattice<Fobj>> {
class MultiGridPreconditioner : public LinearFunction< Lattice<Fobj> > {
public: public:
typedef Aggregation<Fobj, CComplex, nbasis> Aggregates;
typedef CoarsenedMatrix<Fobj, CComplex, nbasis> CoarseOperator;
typedef Aggregation<Fobj,CComplex,nbasis> Aggregates; typedef typename Aggregation<Fobj, CComplex, nbasis>::siteVector siteVector;
typedef CoarsenedMatrix<Fobj,CComplex,nbasis> CoarseOperator; typedef typename Aggregation<Fobj, CComplex, nbasis>::CoarseScalar CoarseScalar;
typedef typename Aggregation<Fobj, CComplex, nbasis>::CoarseVector CoarseVector;
typedef typename Aggregation<Fobj,CComplex,nbasis>::siteVector siteVector; typedef typename Aggregation<Fobj, CComplex, nbasis>::CoarseMatrix CoarseMatrix;
typedef typename Aggregation<Fobj,CComplex,nbasis>::CoarseScalar CoarseScalar; typedef typename Aggregation<Fobj, CComplex, nbasis>::FineField FineField;
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; typedef LinearOperatorBase<FineField> FineOperator;
Aggregates & _Aggregates; Aggregates & _Aggregates;
CoarseOperator & _CoarseOperator; CoarseOperator &_CoarseOperator;
Matrix & _FineMatrix; Matrix & _FineMatrix;
FineOperator & _FineOperator; FineOperator & _FineOperator;
Matrix & _SmootherMatrix; Matrix & _SmootherMatrix;
FineOperator & _SmootherOperator; FineOperator & _SmootherOperator;
// Constructor // Constructor
MultiGridPreconditioner(Aggregates &Agg, CoarseOperator &Coarse, MultiGridPreconditioner(Aggregates & Agg,
FineOperator &Fine,Matrix &FineMatrix, CoarseOperator &Coarse,
FineOperator &Smooth,Matrix &SmootherMatrix) FineOperator & Fine,
: _Aggregates(Agg), Matrix & FineMatrix,
_CoarseOperator(Coarse), FineOperator & Smooth,
_FineOperator(Fine), Matrix & SmootherMatrix)
_FineMatrix(FineMatrix), : _Aggregates(Agg)
_SmootherOperator(Smooth), , _CoarseOperator(Coarse)
_SmootherMatrix(SmootherMatrix) , _FineOperator(Fine)
{ , _FineMatrix(FineMatrix)
} , _SmootherOperator(Smooth)
, _SmootherMatrix(SmootherMatrix) {}
void PowerMethod(const FineField &in) { void PowerMethod(const FineField &in) {
FineField p1(in._grid); FineField p1(in._grid);
FineField p2(in._grid); FineField p2(in._grid);
MdagMLinearOperator<Matrix,FineField> fMdagMOp(_FineMatrix); MdagMLinearOperator<Matrix, FineField> fMdagMOp(_FineMatrix);
p1=in; p1 = in;
RealD absp2; RealD absp2;
for(int i=0;i<20;i++){ for(int i = 0; i < 20; i++) {
RealD absp1=std::sqrt(norm2(p1)); RealD absp1 = std::sqrt(norm2(p1));
fMdagMOp.HermOp(p1,p2);// this is the G5 herm bit fMdagMOp.HermOp(p1, p2); // this is the G5 herm bit
// _FineOperator.Op(p1,p2);// this is the G5 herm bit // _FineOperator.Op(p1,p2); // this is the G5 herm bit
RealD absp2=std::sqrt(norm2(p2)); RealD absp2 = std::sqrt(norm2(p2));
if(i%10==9) if(i % 10 == 9)
std::cout<<GridLogMessage << "Power method on mdagm "<<i<<" " << absp2/absp1<<std::endl; std::cout << GridLogMessage << "Power method on mdagm " << i << " " << absp2 / absp1 << std::endl;
p1=p2*(1.0/std::sqrt(absp2)); p1 = p2 * (1.0 / std::sqrt(absp2));
} }
} }
void operator()(const FineField &in, FineField & out) { void operator()(const FineField &in, FineField &out) {
if ( params.domaindecompose ) { if(params.domaindecompose) {
operatorSAP(in,out); operatorSAP(in, out);
} else { } else {
operatorCheby(in,out); operatorCheby(in, out);
} }
} }
//////////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////////
// ADEF2: [PTM+Q] in = [1 - Q A] M in + Q in = Min + Q [ in -A Min] // 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 // ADEF1: [MP+Q ] in = M [1 - A Q] in + Q in
//////////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////////
#if 1 #if 1
void operatorADEF2(const FineField &in, FineField & out) { void operatorADEF2(const FineField &in, FineField &out) {
CoarseVector Csrc(_CoarseOperator.Grid()); CoarseVector Csrc(_CoarseOperator.Grid());
CoarseVector Ctmp(_CoarseOperator.Grid()); CoarseVector Ctmp(_CoarseOperator.Grid());
CoarseVector Csol(_CoarseOperator.Grid()); CoarseVector Csol(_CoarseOperator.Grid());
ConjugateGradient<CoarseVector> CG(1.0e-10,100000); ConjugateGradient<CoarseVector> CG(1.0e-10, 100000);
ConjugateGradient<FineField> fCG(3.0e-2,1000); ConjugateGradient<FineField> fCG(3.0e-2, 1000);
HermitianLinearOperator<CoarseOperator,CoarseVector> HermOp(_CoarseOperator); HermitianLinearOperator<CoarseOperator, CoarseVector> HermOp(_CoarseOperator);
MdagMLinearOperator<CoarseOperator,CoarseVector> MdagMOp(_CoarseOperator); MdagMLinearOperator<CoarseOperator, CoarseVector> MdagMOp(_CoarseOperator);
MdagMLinearOperator<Matrix,FineField> fMdagMOp(_FineMatrix); MdagMLinearOperator<Matrix, FineField> fMdagMOp(_FineMatrix);
FineField tmp(in._grid); FineField tmp(in._grid);
FineField res(in._grid); FineField res(in._grid);
FineField Min(in._grid); FineField Min(in._grid);
// Monitor completeness of low mode space // Monitor completeness of low mode space
_Aggregates.ProjectToSubspace (Csrc,in); _Aggregates.ProjectToSubspace(Csrc, in);
_Aggregates.PromoteFromSubspace(Csrc,out); _Aggregates.PromoteFromSubspace(Csrc, out);
std::cout<<GridLogMessage<<"Coarse Grid Preconditioner\nCompleteness in: "<<std::sqrt(norm2(out)/norm2(in))<<std::endl; 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] // [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 _FineOperator.Op(in, tmp); // this is the G5 herm bit
fCG(fMdagMOp,tmp,Min); // solves MdagM = g5 M g5M fCG(fMdagMOp, tmp, Min); // solves MdagM = g5 M g5M
// Monitor completeness of low mode space // Monitor completeness of low mode space
_Aggregates.ProjectToSubspace (Csrc,Min); _Aggregates.ProjectToSubspace(Csrc, Min);
_Aggregates.PromoteFromSubspace(Csrc,out); _Aggregates.PromoteFromSubspace(Csrc, out);
std::cout<<GridLogMessage<<"Completeness Min: "<<std::sqrt(norm2(out)/norm2(Min))<<std::endl; std::cout << GridLogMessage << "Completeness Min: " << std::sqrt(norm2(out) / norm2(Min)) << std::endl;
_FineOperator.Op(Min,tmp); _FineOperator.Op(Min, tmp);
tmp = in - tmp; // in - A Min tmp = in - tmp; // in - A Min
Csol=zero; Csol = zero;
_Aggregates.ProjectToSubspace (Csrc,tmp); _Aggregates.ProjectToSubspace(Csrc, tmp);
HermOp.AdjOp(Csrc,Ctmp);// Normal equations HermOp.AdjOp(Csrc, Ctmp); // Normal equations
CG(MdagMOp,Ctmp,Csol); CG(MdagMOp, Ctmp, Csol);
HermOp.Op(Csol,Ctmp); HermOp.Op(Csol, Ctmp);
Ctmp=Ctmp-Csrc; Ctmp = Ctmp - Csrc;
std::cout<<GridLogMessage<<"coarse space true residual "<<std::sqrt(norm2(Ctmp)/norm2(Csrc))<<std::endl; std::cout << GridLogMessage << "coarse space true residual " << std::sqrt(norm2(Ctmp) / norm2(Csrc)) << std::endl;
_Aggregates.PromoteFromSubspace(Csol,out); _Aggregates.PromoteFromSubspace(Csol, out);
_FineOperator.Op(out,res); _FineOperator.Op(out, res);
res=res-tmp; res = res - tmp;
std::cout<<GridLogMessage<<"promoted sol residual "<<std::sqrt(norm2(res)/norm2(tmp))<<std::endl; std::cout << GridLogMessage << "promoted sol residual " << std::sqrt(norm2(res) / norm2(tmp)) << std::endl;
_Aggregates.ProjectToSubspace (Csrc,res); _Aggregates.ProjectToSubspace(Csrc, res);
std::cout<<GridLogMessage<<"coarse space proj of residual "<<norm2(Csrc)<<std::endl; std::cout << GridLogMessage << "coarse space proj of residual " << norm2(Csrc) << std::endl;
out = out + Min; // additive coarse space correction
out = out+Min; // additive coarse space correction
// out = Min; // no additive coarse space correction // out = Min; // no additive coarse space correction
_FineOperator.Op(out,tmp); _FineOperator.Op(out, tmp);
tmp=tmp-in; // tmp is new residual 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;
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 #endif
// ADEF1: [MP+Q ] in =M [1 - A Q] in + Q in // ADEF1: [MP+Q ] in = M [1 - A Q] in + Q in
#if 1 #if 1
void operatorADEF1(const FineField &in, FineField & out) { void operatorADEF1(const FineField &in, FineField &out) {
CoarseVector Csrc(_CoarseOperator.Grid()); CoarseVector Csrc(_CoarseOperator.Grid());
CoarseVector Ctmp(_CoarseOperator.Grid()); CoarseVector Ctmp(_CoarseOperator.Grid());
CoarseVector Csol(_CoarseOperator.Grid()); Csol=zero; CoarseVector Csol(_CoarseOperator.Grid());
Csol = zero;
ConjugateGradient<CoarseVector> CG(1.0e-10,100000); ConjugateGradient<CoarseVector> CG(1.0e-10, 100000);
ConjugateGradient<FineField> fCG(3.0e-2,1000); ConjugateGradient<FineField> fCG(3.0e-2, 1000);
HermitianLinearOperator<CoarseOperator,CoarseVector> HermOp(_CoarseOperator); HermitianLinearOperator<CoarseOperator, CoarseVector> HermOp(_CoarseOperator);
MdagMLinearOperator<CoarseOperator,CoarseVector> MdagMOp(_CoarseOperator); MdagMLinearOperator<CoarseOperator, CoarseVector> MdagMOp(_CoarseOperator);
ShiftedMdagMLinearOperator<Matrix,FineField> fMdagMOp(_FineMatrix,0.1); ShiftedMdagMLinearOperator<Matrix, FineField> fMdagMOp(_FineMatrix, 0.1);
FineField tmp(in._grid); FineField tmp(in._grid);
FineField res(in._grid); FineField res(in._grid);
@ -292,139 +287,138 @@ public:
// _Aggregates.PromoteFromSubspace(Csrc,out); // _Aggregates.PromoteFromSubspace(Csrc,out);
// std::cout<<GridLogMessage<<"Coarse Grid Preconditioner\nCompleteness in: "<<std::sqrt(norm2(out)/norm2(in))<<std::endl; // std::cout<<GridLogMessage<<"Coarse Grid Preconditioner\nCompleteness in: "<<std::sqrt(norm2(out)/norm2(in))<<std::endl;
_Aggregates.ProjectToSubspace (Csrc,in); _Aggregates.ProjectToSubspace(Csrc, in);
HermOp.AdjOp(Csrc,Ctmp);// Normal equations HermOp.AdjOp(Csrc, Ctmp); // Normal equations
CG(MdagMOp,Ctmp,Csol); CG(MdagMOp, Ctmp, Csol);
_Aggregates.PromoteFromSubspace(Csol,Qin); _Aggregates.PromoteFromSubspace(Csol, Qin);
// Qin=0; // Qin=0;
_FineOperator.Op(Qin,tmp);// A Q in _FineOperator.Op(Qin, tmp); // A Q in
tmp = in - tmp; // in - A Q in tmp = in - tmp; // in - A Q in
_FineOperator.Op(tmp,res);// this is the G5 herm bit _FineOperator.Op(tmp, res); // this is the G5 herm bit
fCG(fMdagMOp,res,out); // solves MdagM = g5 M g5M fCG(fMdagMOp, res, out); // solves MdagM = g5 M g5M
out = out + Qin; out = out + Qin;
_FineOperator.Op(out,tmp); _FineOperator.Op(out, tmp);
tmp=tmp-in; // tmp is new residual tmp = tmp - in; // tmp is new residual
std::cout<<GridLogMessage<<"preconditioner thinks residual is "<<std::sqrt(norm2(tmp)/norm2(in))<<std::endl;
std::cout << GridLogMessage << "preconditioner thinks residual is " << std::sqrt(norm2(tmp) / norm2(in)) << std::endl;
} }
#endif #endif
void SAP (const FineField & src,FineField & psi){ void SAP(const FineField &src, FineField &psi) {
Lattice<iScalar<vInteger> > coor(src._grid); Lattice<iScalar<vInteger>> coor(src._grid);
Lattice<iScalar<vInteger> > subset(src._grid); Lattice<iScalar<vInteger>> subset(src._grid);
FineField r(src._grid); FineField r(src._grid);
FineField zz(src._grid); zz=zero; FineField zz(src._grid);
zz = zero;
FineField vec1(src._grid); FineField vec1(src._grid);
FineField vec2(src._grid); FineField vec2(src._grid);
const Integer block=params.domainsize; const Integer block = params.domainsize;
subset=zero; subset = zero;
for(int mu=0;mu<Nd;mu++){ for(int mu = 0; mu < Nd; mu++) {
LatticeCoordinate(coor,mu+1); LatticeCoordinate(coor, mu + 1);
coor = div(coor,block); coor = div(coor, block);
subset = subset+coor; subset = subset + coor;
} }
subset = mod(subset,(Integer)2); subset = mod(subset, (Integer)2);
ShiftedMdagMLinearOperator<Matrix,FineField> fMdagMOp(_SmootherMatrix,0.0); ShiftedMdagMLinearOperator<Matrix, FineField> fMdagMOp(_SmootherMatrix, 0.0);
Chebyshev<FineField> Cheby (params.lo,params.hi,params.order,InverseApproximation); Chebyshev<FineField> Cheby(params.lo, params.hi, params.order, InverseApproximation);
RealD resid; RealD resid;
for(int i=0;i<params.steps;i++){ for(int i = 0; i < params.steps; i++) {
// Even domain residual // Even domain residual
_FineOperator.Op(psi,vec1);// this is the G5 herm bit _FineOperator.Op(psi, vec1); // this is the G5 herm bit
r= src - vec1 ; r = src - vec1;
resid = norm2(r) /norm2(src); resid = norm2(r) / norm2(src);
std::cout << "SAP "<<i<<" resid "<<resid<<std::endl; std::cout << "SAP " << i << " resid " << resid << std::endl;
// Even domain solve // Even domain solve
r= where(subset==(Integer)0,r,zz); r = where(subset == (Integer)0, r, zz);
_SmootherOperator.AdjOp(r,vec1); _SmootherOperator.AdjOp(r, vec1);
Cheby(fMdagMOp,vec1,vec2); // solves MdagM = g5 M g5M Cheby(fMdagMOp, vec1, vec2); // solves MdagM = g5 M g5M
psi = psi + vec2; psi = psi + vec2;
// Odd domain residual // Odd domain residual
_FineOperator.Op(psi,vec1);// this is the G5 herm bit _FineOperator.Op(psi, vec1); // this is the G5 herm bit
r= src - vec1 ; r = src - vec1;
r= where(subset==(Integer)1,r,zz); r = where(subset == (Integer)1, r, zz);
resid = norm2(r) /norm2(src); resid = norm2(r) / norm2(src);
std::cout << "SAP "<<i<<" resid "<<resid<<std::endl; std::cout << "SAP " << i << " resid " << resid << std::endl;
// Odd domain solve // Odd domain solve
_SmootherOperator.AdjOp(r,vec1); _SmootherOperator.AdjOp(r, vec1);
Cheby(fMdagMOp,vec1,vec2); // solves MdagM = g5 M g5M Cheby(fMdagMOp, vec1, vec2); // solves MdagM = g5 M g5M
psi = psi + vec2; psi = psi + vec2;
_FineOperator.Op(psi,vec1);// this is the G5 herm bit _FineOperator.Op(psi, vec1); // this is the G5 herm bit
r= src - vec1 ; r = src - vec1;
resid = norm2(r) /norm2(src); resid = norm2(r) / norm2(src);
std::cout << "SAP "<<i<<" resid "<<resid<<std::endl; std::cout << "SAP " << i << " resid " << resid << std::endl;
} }
}; };
void SmootherTest (const FineField & in){ void SmootherTest(const FineField &in) {
FineField vec1(in._grid); FineField vec1(in._grid);
FineField vec2(in._grid); FineField vec2(in._grid);
RealD lo[3] = { 0.5, 1.0, 2.0};
RealD lo[3] = {0.5, 1.0, 2.0};
// MdagMLinearOperator<Matrix,FineField> fMdagMOp(_FineMatrix); // MdagMLinearOperator<Matrix,FineField> fMdagMOp(_FineMatrix);
ShiftedMdagMLinearOperator<Matrix,FineField> fMdagMOp(_SmootherMatrix,0.0); ShiftedMdagMLinearOperator<Matrix, FineField> fMdagMOp(_SmootherMatrix, 0.0);
RealD Ni,r; RealD Ni, r;
Ni = norm2(in); Ni = norm2(in);
for(int ilo=0;ilo<3;ilo++){ for(int ilo = 0; ilo < 3; ilo++) {
for(int ord=5;ord<50;ord*=2){ for(int ord = 5; ord < 50; ord *= 2) {
_SmootherOperator.AdjOp(in,vec1); _SmootherOperator.AdjOp(in, vec1);
Chebyshev<FineField> Cheby (lo[ilo],70.0,ord,InverseApproximation); Chebyshev<FineField> Cheby(lo[ilo], 70.0, ord, InverseApproximation);
Cheby(fMdagMOp,vec1,vec2); // solves MdagM = g5 M g5M Cheby(fMdagMOp, vec1, vec2); // solves MdagM = g5 M g5M
_FineOperator.Op(vec2,vec1);// this is the G5 herm bit _FineOperator.Op(vec2, vec1); // this is the G5 herm bit
vec1 = in - vec1; // tmp = in - A Min vec1 = in - vec1; // tmp = in - A Min
r=norm2(vec1); r = norm2(vec1);
std::cout<<GridLogMessage << "Smoother resid "<<std::sqrt(r/Ni)<<std::endl; std::cout << GridLogMessage << "Smoother resid " << std::sqrt(r / Ni) << std::endl;
} }
} }
} }
void operatorCheby(const FineField &in, FineField & out) { void operatorCheby(const FineField &in, FineField &out) {
CoarseVector Csrc(_CoarseOperator.Grid()); CoarseVector Csrc(_CoarseOperator.Grid());
CoarseVector Ctmp(_CoarseOperator.Grid()); CoarseVector Ctmp(_CoarseOperator.Grid());
CoarseVector Csol(_CoarseOperator.Grid()); Csol=zero; CoarseVector Csol(_CoarseOperator.Grid());
Csol = zero;
ConjugateGradient<CoarseVector> CG(3.0e-3,100000); ConjugateGradient<CoarseVector> CG(3.0e-3, 100000);
// ConjugateGradient<FineField> fCG(3.0e-2,1000); // ConjugateGradient<FineField> fCG(3.0e-2,1000);
HermitianLinearOperator<CoarseOperator,CoarseVector> HermOp(_CoarseOperator); HermitianLinearOperator<CoarseOperator, CoarseVector> HermOp(_CoarseOperator);
MdagMLinearOperator<CoarseOperator,CoarseVector> MdagMOp(_CoarseOperator); MdagMLinearOperator<CoarseOperator, CoarseVector> MdagMOp(_CoarseOperator);
// MdagMLinearOperator<Matrix,FineField> fMdagMOp(_FineMatrix); // MdagMLinearOperator<Matrix,FineField> fMdagMOp(_FineMatrix);
ShiftedMdagMLinearOperator<Matrix,FineField> fMdagMOp(_SmootherMatrix,0.0); ShiftedMdagMLinearOperator<Matrix, FineField> fMdagMOp(_SmootherMatrix, 0.0);
FineField vec1(in._grid); FineField vec1(in._grid);
FineField vec2(in._grid); FineField vec2(in._grid);
// Chebyshev<FineField> Cheby (0.5,70.0,30,InverseApproximation); // Chebyshev<FineField> Cheby (0.5,70.0,30,InverseApproximation);
// Chebyshev<FineField> ChebyAccu(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> Cheby(params.lo, params.hi, params.order, InverseApproximation);
Chebyshev<FineField> ChebyAccu(params.lo,params.hi,params.order,InverseApproximation); Chebyshev<FineField> ChebyAccu(params.lo, params.hi, params.order, InverseApproximation);
// Cheby.JacksonSmooth(); // Cheby.JacksonSmooth();
// ChebyAccu.JacksonSmooth(); // ChebyAccu.JacksonSmooth();
@ -445,130 +439,123 @@ public:
RealD Ni = norm2(in); RealD Ni = norm2(in);
_SmootherOperator.AdjOp(in,vec1);// this is the G5 herm bit _SmootherOperator.AdjOp(in, vec1); // this is the G5 herm bit
ChebyAccu(fMdagMOp,vec1,out); // solves MdagM = g5 M g5M ChebyAccu(fMdagMOp, vec1, out); // solves MdagM = g5 M g5M
std::cout<<GridLogMessage << "Smoother norm "<<norm2(out)<<std::endl; std::cout << GridLogMessage << "Smoother norm " << norm2(out) << std::endl;
// Update with residual for out // Update with residual for out
_FineOperator.Op(out,vec1);// this is the G5 herm bit _FineOperator.Op(out, vec1); // this is the G5 herm bit
vec1 = in - vec1; // tmp = in - A Min vec1 = in - vec1; // tmp = in - A Min
RealD r = norm2(vec1); RealD r = norm2(vec1);
std::cout<<GridLogMessage << "Smoother resid "<<std::sqrt(r/Ni)<< " " << r << " " << Ni <<std::endl; std::cout << GridLogMessage << "Smoother resid " << std::sqrt(r / Ni) << " " << r << " " << Ni << std::endl;
_Aggregates.ProjectToSubspace (Csrc,vec1); _Aggregates.ProjectToSubspace(Csrc, vec1);
HermOp.AdjOp(Csrc,Ctmp);// Normal equations HermOp.AdjOp(Csrc, Ctmp); // Normal equations
CG(MdagMOp,Ctmp,Csol); CG(MdagMOp, Ctmp, Csol);
_Aggregates.PromoteFromSubspace(Csol,vec1); // Ass^{-1} [in - A Min]_s _Aggregates.PromoteFromSubspace(Csol, vec1); // Ass^{-1} [in - A Min]_s
// Q = Q[in - A Min] // Q = Q[in - A Min]
out = out+vec1; out = out + vec1;
// Three preconditioner smoothing -- hermitian if C3 = C1 // Three preconditioner smoothing -- hermitian if C3 = C1
// Recompute error // Recompute error
_FineOperator.Op(out,vec1);// this is the G5 herm bit _FineOperator.Op(out, vec1); // this is the G5 herm bit
vec1 = in - vec1; // tmp = in - A Min vec1 = in - vec1; // tmp = in - A Min
r=norm2(vec1); r = norm2(vec1);
std::cout<<GridLogMessage << "Coarse resid "<<std::sqrt(r/Ni)<<std::endl; std::cout << GridLogMessage << "Coarse resid " << std::sqrt(r / Ni) << std::endl;
// Reapply smoother // Reapply smoother
_SmootherOperator.Op(vec1,vec2); // this is the G5 herm bit _SmootherOperator.Op(vec1, vec2); // this is the G5 herm bit
ChebyAccu(fMdagMOp,vec2,vec1); // solves MdagM = g5 M g5M ChebyAccu(fMdagMOp, vec2, vec1); // solves MdagM = g5 M g5M
out =out+vec1; out = out + vec1;
vec1 = in - vec1; // tmp = in - A Min vec1 = in - vec1; // tmp = in - A Min
r=norm2(vec1); r = norm2(vec1);
std::cout<<GridLogMessage << "Smoother resid "<<std::sqrt(r/Ni)<<std::endl; std::cout << GridLogMessage << "Smoother resid " << std::sqrt(r / Ni) << std::endl;
} }
void operatorSAP(const FineField &in, FineField & out) { void operatorSAP(const FineField &in, FineField &out) {
CoarseVector Csrc(_CoarseOperator.Grid()); CoarseVector Csrc(_CoarseOperator.Grid());
CoarseVector Ctmp(_CoarseOperator.Grid()); CoarseVector Ctmp(_CoarseOperator.Grid());
CoarseVector Csol(_CoarseOperator.Grid()); Csol=zero; CoarseVector Csol(_CoarseOperator.Grid());
Csol = zero;
ConjugateGradient<CoarseVector> CG(1.0e-3,100000); ConjugateGradient<CoarseVector> CG(1.0e-3, 100000);
HermitianLinearOperator<CoarseOperator,CoarseVector> HermOp(_CoarseOperator); HermitianLinearOperator<CoarseOperator, CoarseVector> HermOp(_CoarseOperator);
MdagMLinearOperator<CoarseOperator,CoarseVector> MdagMOp(_CoarseOperator); MdagMLinearOperator<CoarseOperator, CoarseVector> MdagMOp(_CoarseOperator);
FineField vec1(in._grid); FineField vec1(in._grid);
FineField vec2(in._grid); FineField vec2(in._grid);
_Aggregates.ProjectToSubspace (Csrc,in); _Aggregates.ProjectToSubspace(Csrc, in);
_Aggregates.PromoteFromSubspace(Csrc,out); _Aggregates.PromoteFromSubspace(Csrc, out);
std::cout<<GridLogMessage<<"Completeness: "<<std::sqrt(norm2(out)/norm2(in))<<std::endl; std::cout << GridLogMessage << "Completeness: " << std::sqrt(norm2(out) / norm2(in)) << std::endl;
// To make a working smoother for indefinite operator // To make a working smoother for indefinite operator
// must multiply by "Mdag" (ouch loses all low mode content) // must multiply by "Mdag" (ouch loses all low mode content)
// and apply to poly approx of (mdagm)^-1. // and apply to poly approx of (mdagm)^-1.
// so that we end up with an odd polynomial. // so that we end up with an odd polynomial.
SAP(in,out); SAP(in, out);
// Update with residual for out // Update with residual for out
_FineOperator.Op(out,vec1);// this is the G5 herm bit _FineOperator.Op(out, vec1); // this is the G5 herm bit
vec1 = in - vec1; // tmp = in - A Min vec1 = in - vec1; // tmp = in - A Min
RealD r = norm2(vec1); RealD r = norm2(vec1);
RealD Ni = norm2(in); RealD Ni = norm2(in);
std::cout<<GridLogMessage << "SAP resid "<<std::sqrt(r/Ni)<< " " << r << " " << Ni <<std::endl; std::cout << GridLogMessage << "SAP resid " << std::sqrt(r / Ni) << " " << r << " " << Ni << std::endl;
_Aggregates.ProjectToSubspace (Csrc,vec1); _Aggregates.ProjectToSubspace(Csrc, vec1);
HermOp.AdjOp(Csrc,Ctmp);// Normal equations HermOp.AdjOp(Csrc, Ctmp); // Normal equations
CG(MdagMOp,Ctmp,Csol); CG(MdagMOp, Ctmp, Csol);
_Aggregates.PromoteFromSubspace(Csol,vec1); // Ass^{-1} [in - A Min]_s _Aggregates.PromoteFromSubspace(Csol, vec1); // Ass^{-1} [in - A Min]_s
// Q = Q[in - A Min] // Q = Q[in - A Min]
out = out+vec1; out = out + vec1;
// Three preconditioner smoothing -- hermitian if C3 = C1 // Three preconditioner smoothing -- hermitian if C3 = C1
// Recompute error // Recompute error
_FineOperator.Op(out,vec1);// this is the G5 herm bit _FineOperator.Op(out, vec1); // this is the G5 herm bit
vec1 = in - vec1; // tmp = in - A Min vec1 = in - vec1; // tmp = in - A Min
r=norm2(vec1); r = norm2(vec1);
std::cout<<GridLogMessage << "Coarse resid "<<std::sqrt(r/Ni)<<std::endl; std::cout << GridLogMessage << "Coarse resid " << std::sqrt(r / Ni) << std::endl;
// Reapply smoother // Reapply smoother
SAP(vec1,vec2); SAP(vec1, vec2);
out =out+vec2; out = out + vec2;
// Update with residual for out // Update with residual for out
_FineOperator.Op(out,vec1);// this is the G5 herm bit _FineOperator.Op(out, vec1); // this is the G5 herm bit
vec1 = in - vec1; // tmp = in - A Min vec1 = in - vec1; // tmp = in - A Min
r = norm2(vec1); r = norm2(vec1);
Ni = norm2(in); Ni = norm2(in);
std::cout<<GridLogMessage << "SAP resid(post) "<<std::sqrt(r/Ni)<< " " << r << " " << Ni <<std::endl; std::cout << GridLogMessage << "SAP resid(post) " << std::sqrt(r / Ni) << " " << r << " " << Ni << std::endl;
} }
}; };
struct MGParams struct MGParams {
{ std::vector<std::vector<int>> blockSizes;
std::vector< std::vector< int > > blockSizes;
const int nbasis; const int nbasis;
MGParams() MGParams()
: blockSizes( { { 1, 1, 1, 2 } } ) : blockSizes({{1, 1, 1, 2}})
// : blockSizes({ {1,1,1,2}, {1,1,1,2} }) // : blockSizes({{1,1,1,2}, {1,1,1,2}})
// : blockSizes({ {1,1,1,2}, {1,1,1,2}, {1,1,1,2} }) // : blockSizes({{1,1,1,2}, {1,1,1,2}, {1,1,1,2}})
, nbasis( 20 ) , nbasis(20) {}
{
}
}; };
int main (int argc, char ** argv) int main(int argc, char **argv) {
{
Grid_init(&argc,&argv);
params.domainsize= 1; Grid_init(&argc, &argv);
params.coarsegrids= 1;
params.domainsize = 1;
params.coarsegrids = 1;
params.domaindecompose = 0; params.domaindecompose = 0;
params.order = 30; params.order = 30;
params.Ls = 1; params.Ls = 1;
@ -590,16 +577,17 @@ int main (int argc, char ** argv)
std::cout << GridLogMessage << "Set up some fine level stuff: " << std::endl; std::cout << GridLogMessage << "Set up some fine level stuff: " << std::endl;
std::cout << GridLogMessage << "**************************************************" << 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); GridRedBlackCartesian *FrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(FGrid);
std::vector<int> fSeeds( {1, 2, 3, 4} ); std::vector<int> fSeeds({1, 2, 3, 4});
GridParallelRNG fPRNG( FGrid ); GridParallelRNG fPRNG(FGrid);
fPRNG.SeedFixedIntegers( fSeeds ); fPRNG.SeedFixedIntegers(fSeeds);
Gamma g5(Gamma::Algebra::Gamma5); Gamma g5(Gamma::Algebra::Gamma5);
LatticeFermion src(FGrid); gaussian(fPRNG, src); // src=src+g5*src; // clang-format off
LatticeFermion src(FGrid); gaussian(fPRNG, src); // src=src + g5 * src;
LatticeFermion result(FGrid); result = zero; LatticeFermion result(FGrid); result = zero;
LatticeFermion ref(FGrid); ref = zero; LatticeFermion ref(FGrid); ref = zero;
LatticeFermion tmp(FGrid); LatticeFermion tmp(FGrid);
@ -608,29 +596,30 @@ int main (int argc, char ** argv)
LatticeGaugeField UmuDD(FGrid); LatticeGaugeField UmuDD(FGrid);
LatticeColourMatrix U(FGrid); LatticeColourMatrix U(FGrid);
LatticeColourMatrix zz(FGrid); LatticeColourMatrix zz(FGrid);
// clang-format on
if ( params.domaindecompose ) { if(params.domaindecompose) {
Lattice<iScalar<vInteger> > coor(FGrid); Lattice<iScalar<vInteger>> coor(FGrid);
zz=zero; zz = zero;
for(int mu=0;mu<Nd;mu++){ for(int mu = 0; mu < Nd; mu++) {
LatticeCoordinate(coor,mu); LatticeCoordinate(coor, mu);
U = PeekIndex<LorentzIndex>(Umu,mu); U = PeekIndex<LorentzIndex>(Umu, mu);
U = where(mod(coor,params.domainsize)==(Integer)0,zz,U); U = where(mod(coor, params.domainsize) == (Integer)0, zz, U);
PokeIndex<LorentzIndex>(UmuDD,U,mu); PokeIndex<LorentzIndex>(UmuDD, U, mu);
} }
} else { } else {
UmuDD = Umu; UmuDD = Umu;
} }
RealD mass=params.mq; RealD mass = params.mq;
std::cout << GridLogMessage << "**************************************************" << std::endl; std::cout << GridLogMessage << "**************************************************" << std::endl;
std::cout << GridLogMessage << "Set up some coarser levels stuff: " << std::endl; std::cout << GridLogMessage << "Set up some coarser levels stuff: " << std::endl;
std::cout << GridLogMessage << "**************************************************" << std::endl; std::cout << GridLogMessage << "**************************************************" << std::endl;
std::vector< std::vector< int > > blockSizes({ { 1, 1, 1, 2 } } ); // corresponds to two level algorithm std::vector<std::vector<int>> blockSizes({{1, 1, 1, 2}}); // corresponds to two level algorithm
// std::vector< std::vector<int> > blockSizes({ {1,1,1,2}, // // corresponds to three level algorithm // std::vector<std::vector<int>> blockSizes({{1, 1, 1, 2}, // corresponds to three level algorithm
// {1,1,1,2} }); // {1, 1, 1, 2}});
const int nbasis = 20; // we fix the number of test vector to the same const int nbasis = 20; // we fix the number of test vector to the same
// number on every level for now // number on every level for now
@ -659,25 +648,25 @@ int main (int argc, char ** argv)
// GridParallelRNG cPRNG(CGrid); cPRNG.SeedFixedIntegers(cSeeds); // GridParallelRNG cPRNG(CGrid); cPRNG.SeedFixedIntegers(cSeeds);
CoarseGrids< nbasis > cGrids( blockSizes ); CoarseGrids<nbasis> cGrids(blockSizes);
// assert(0); // assert(0);
std::cout<<GridLogMessage << "**************************************************"<< std::endl; std::cout << GridLogMessage << "**************************************************" << std::endl;
std::cout<<GridLogMessage << "Building the wilson operator on the fine grid" <<std::endl; std::cout << GridLogMessage << "Building the wilson operator on the fine grid" << std::endl;
std::cout<<GridLogMessage << "**************************************************"<< std::endl; std::cout << GridLogMessage << "**************************************************" << std::endl;
WilsonFermionR Dw(Umu,*FGrid,*FrbGrid,mass); WilsonFermionR Dw(Umu, *FGrid, *FrbGrid, mass);
WilsonFermionR DwDD(UmuDD,*FGrid,*FrbGrid,mass); WilsonFermionR DwDD(UmuDD, *FGrid, *FrbGrid, mass);
std::cout<<GridLogMessage<< "**************************************************"<< std::endl; std::cout << GridLogMessage << "**************************************************" << std::endl;
std::cout<<GridLogMessage<< "Some typedefs" <<std::endl; std::cout << GridLogMessage << "Some typedefs" << std::endl;
std::cout<<GridLogMessage<< "**************************************************"<< std::endl; std::cout << GridLogMessage << "**************************************************" << std::endl;
typedef Aggregation<vSpinColourVector,vTComplex,nbasis> Subspace; typedef Aggregation<vSpinColourVector, vTComplex, nbasis> Subspace;
typedef CoarsenedMatrix<vSpinColourVector,vTComplex,nbasis> CoarseOperator; typedef CoarsenedMatrix<vSpinColourVector, vTComplex, nbasis> CoarseOperator;
typedef CoarseOperator::CoarseVector CoarseVector; typedef CoarseOperator::CoarseVector CoarseVector;
typedef TestVectorAnalyzer<LatticeFermion,nbasis> TVA; typedef TestVectorAnalyzer<LatticeFermion, nbasis> TVA;
// typedef Aggregation<vSpinColourVector,vTComplex,1,nbasis> Subspace; // typedef Aggregation<vSpinColourVector,vTComplex,1,nbasis> Subspace;
// typedef CoarsenedMatrix<vSpinColourVector,vTComplex,1,nbasis> CoarseOperator; // typedef CoarsenedMatrix<vSpinColourVector,vTComplex,1,nbasis> CoarseOperator;
@ -688,6 +677,7 @@ int main (int argc, char ** argv)
// CoarseOperator::CoarseG5PMatrix CoarseG5PMatrix; // CoarseOperator::CoarseG5PMatrix CoarseG5PMatrix;
#if 1 #if 1
// clang-format off
std::cout << std::endl; std::cout << std::endl;
std::cout << "type_name<decltype(vTComplex{})>() = " << type_name<decltype(vTComplex{})>() << std::endl; std::cout << "type_name<decltype(vTComplex{})>() = " << type_name<decltype(vTComplex{})>() << std::endl;
std::cout << "type_name<GridTypeMapper<vTComplex>::scalar_type>() = " << type_name<GridTypeMapper<vTComplex>::scalar_type>() << std::endl; std::cout << "type_name<GridTypeMapper<vTComplex>::scalar_type>() = " << type_name<GridTypeMapper<vTComplex>::scalar_type>() << std::endl;
@ -711,22 +701,23 @@ int main (int argc, char ** argv)
std::cout << "type_name<GridTypeMapper<TComplex>::Realified>() = " << type_name<GridTypeMapper<TComplex>::Realified>() << std::endl; std::cout << "type_name<GridTypeMapper<TComplex>::Realified>() = " << type_name<GridTypeMapper<TComplex>::Realified>() << std::endl;
std::cout << "type_name<GridTypeMapper<TComplex>::DoublePrecision>() = " << type_name<GridTypeMapper<TComplex>::DoublePrecision>() << std::endl; std::cout << "type_name<GridTypeMapper<TComplex>::DoublePrecision>() = " << type_name<GridTypeMapper<TComplex>::DoublePrecision>() << std::endl;
std::cout << std::endl; std::cout << std::endl;
// clang-format on
#endif #endif
std::cout<<GridLogMessage << "**************************************************"<< std::endl; std::cout << GridLogMessage << "**************************************************" << std::endl;
std::cout<<GridLogMessage << "Calling Aggregation class to build subspaces" <<std::endl; std::cout << GridLogMessage << "Calling Aggregation class to build subspaces" << std::endl;
std::cout<<GridLogMessage << "**************************************************"<< std::endl; std::cout << GridLogMessage << "**************************************************" << std::endl;
// • TODO: need some way to run the smoother on the "test vectors" for a few // • TODO: need some way to run the smoother on the "test vectors" for a few
// times before constructing the subspace from them // times before constructing the subspace from them
// • Maybe an application for an mrhs (true mrhs, no block) smoother? // • Maybe an application for an mrhs (true mrhs, no block) smoother?
// • In WMG, the vectors are normalized but not orthogonalized, but here they // • In WMG, the vectors are normalized but not orthogonalized, but here they
// are constructed randomly and then orthogonalized (rather orthonormalized) against each other // are constructed randomly and then orthogonalized (rather orthonormalized) against each other
MdagMLinearOperator<WilsonFermionR,LatticeFermion> HermOp(Dw); MdagMLinearOperator<WilsonFermionR, LatticeFermion> HermOp(Dw);
Subspace Aggregates(cGrids.Grids[0],FGrid,0); Subspace Aggregates(cGrids.Grids[0], FGrid, 0);
assert ((nbasis & 0x1)==0); assert((nbasis & 0x1) == 0);
int nb=nbasis/2; int nb = nbasis / 2;
std::cout<<GridLogMessage << " nbasis/2 = "<<nb<<std::endl; std::cout << GridLogMessage << " nbasis/2 = " << nb << std::endl;
Aggregates.CreateSubspace(fPRNG, HermOp /*, nb */); // Don't specify nb to see the orthogonalization check Aggregates.CreateSubspace(fPRNG, HermOp /*, nb */); // Don't specify nb to see the orthogonalization check
@ -734,39 +725,42 @@ int main (int argc, char ** argv)
testVectorAnalyzer(HermOp, Aggregates.subspace, nb); testVectorAnalyzer(HermOp, Aggregates.subspace, nb);
for(int n=0;n<nb;n++){ for(int n = 0; n < nb; n++) {
Aggregates.subspace[n+nb] = g5 * Aggregates.subspace[n]; // multiply with g5 normally instead of G5R5 since this specific to DWF // multiply with g5 normally instead of G5R5 since this specific to DWF
std::cout<<GridLogMessage<<n<<" subspace "<<norm2(Aggregates.subspace[n+nb])<<" "<<norm2(Aggregates.subspace[n]) <<std::endl; Aggregates.subspace[n + nb] = g5 * Aggregates.subspace[n];
std::cout << GridLogMessage << n << " subspace " << norm2(Aggregates.subspace[n + nb]) << " " << norm2(Aggregates.subspace[n])
<< std::endl;
} }
for(int n=0;n<nbasis;n++){ for(int n = 0; n < nbasis; n++) {
std::cout<<GridLogMessage << "vec["<<n<<"] = "<<norm2(Aggregates.subspace[n]) <<std::endl; std::cout << GridLogMessage << "vec[" << n << "] = " << norm2(Aggregates.subspace[n]) << std::endl;
} }
// tva(HermOp, Aggregates.subspace); // tva(HermOp, Aggregates.subspace);
Aggregates.CheckOrthogonal(); Aggregates.CheckOrthogonal();
testVectorAnalyzer(HermOp, Aggregates.subspace); testVectorAnalyzer(HermOp, Aggregates.subspace);
result=zero; result = zero;
std::cout<<GridLogMessage << "**************************************************"<< std::endl; std::cout << GridLogMessage << "**************************************************" << std::endl;
std::cout<<GridLogMessage << "Building coarse representation of Dirac operator" <<std::endl; std::cout << GridLogMessage << "Building coarse representation of Dirac operator" << std::endl;
std::cout<<GridLogMessage << "**************************************************"<< std::endl; std::cout << GridLogMessage << "**************************************************" << std::endl;
Gamma5HermitianLinearOperator<WilsonFermionR,LatticeFermion> HermIndefOp(Dw); // this corresponds to working with H = g5 * D // using Gamma5HermitianLinearOperator corresponds to working with H = g5 * D
Gamma5HermitianLinearOperator<WilsonFermionR,LatticeFermion> HermIndefOpDD(DwDD); Gamma5HermitianLinearOperator<WilsonFermionR, LatticeFermion> HermIndefOp(Dw);
CoarsenedMatrix<vSpinColourVector,vTComplex,nbasis> CoarseOp(*cGrids.Grids[0]); Gamma5HermitianLinearOperator<WilsonFermionR, LatticeFermion> HermIndefOpDD(DwDD);
CoarsenedMatrix<vSpinColourVector, vTComplex, nbasis> CoarseOp(*cGrids.Grids[0]);
CoarseOp.CoarsenOperator(FGrid, HermIndefOp, Aggregates); // uses only linop.OpDiag & linop.OpDir CoarseOp.CoarsenOperator(FGrid, HermIndefOp, Aggregates); // uses only linop.OpDiag & linop.OpDir
std::cout << GridLogMessage << "**************************************************" << std::endl; std::cout << GridLogMessage << "**************************************************" << std::endl;
std::cout << GridLogMessage << "Building coarse vectors" << std::endl; std::cout << GridLogMessage << "Building coarse vectors" << std::endl;
std::cout << GridLogMessage << "**************************************************" << std::endl; std::cout << GridLogMessage << "**************************************************" << std::endl;
CoarseVector c_src (cGrids.Grids[0]); CoarseVector c_src(cGrids.Grids[0]);
CoarseVector c_res (cGrids.Grids[0]); CoarseVector c_res(cGrids.Grids[0]);
gaussian(cGrids.PRNGs[0],c_src); gaussian(cGrids.PRNGs[0], c_src);
c_res=zero; c_res = zero;
std::cout << "type_name<decltype(c_src)>() = " << type_name< decltype( c_src ) >() << std::endl; std::cout << "type_name<decltype(c_src)>() = " << type_name<decltype(c_src)>() << std::endl;
// c_res = g5 * c_src; // c_res = g5 * c_src;
@ -774,7 +768,7 @@ int main (int argc, char ** argv)
std::cout << GridLogMessage << "Solving posdef-MR on coarse space " << std::endl; std::cout << GridLogMessage << "Solving posdef-MR on coarse space " << std::endl;
std::cout << GridLogMessage << "**************************************************" << std::endl; std::cout << GridLogMessage << "**************************************************" << std::endl;
MdagMLinearOperator<CoarseOperator,CoarseVector> PosdefLdop(CoarseOp); MdagMLinearOperator<CoarseOperator, CoarseVector> PosdefLdop(CoarseOp);
MinimalResidual<CoarseVector> MR(5.0e-2, 100, false); MinimalResidual<CoarseVector> MR(5.0e-2, 100, false);
ConjugateGradient<CoarseVector> CG(5.0e-2, 100, false); ConjugateGradient<CoarseVector> CG(5.0e-2, 100, false);
@ -811,17 +805,15 @@ int main (int argc, char ** argv)
// ConjugateResidual<CoarseVector> MCR(1.0e-6,100000); // ConjugateResidual<CoarseVector> MCR(1.0e-6,100000);
// MCR(HermIndefLdop,c_src,c_res); // MCR(HermIndefLdop,c_src,c_res);
std::cout<<GridLogMessage << "**************************************************"<< std::endl; std::cout << GridLogMessage << "**************************************************" << std::endl;
std::cout<<GridLogMessage << "Building deflation preconditioner "<< std::endl; std::cout << GridLogMessage << "Building deflation preconditioner " << std::endl;
std::cout<<GridLogMessage << "**************************************************"<< std::endl; std::cout << GridLogMessage << "**************************************************" << std::endl;
MultiGridPreconditioner <vSpinColourVector,vTComplex,nbasis,WilsonFermionR> Precon (Aggregates, CoarseOp, MultiGridPreconditioner<vSpinColourVector, vTComplex, nbasis, WilsonFermionR> Precon(
HermIndefOp,Dw, Aggregates, CoarseOp, HermIndefOp, Dw, HermIndefOp, Dw);
HermIndefOp,Dw);
MultiGridPreconditioner <vSpinColourVector,vTComplex,nbasis,WilsonFermionR> PreconDD(Aggregates, CoarseOp, MultiGridPreconditioner<vSpinColourVector, vTComplex, nbasis, WilsonFermionR> PreconDD(
HermIndefOp,Dw, Aggregates, CoarseOp, HermIndefOp, Dw, HermIndefOpDD, DwDD);
HermIndefOpDD,DwDD);
// MultiGridPreconditioner(Aggregates &Agg, CoarseOperator &Coarse, // MultiGridPreconditioner(Aggregates &Agg, CoarseOperator &Coarse,
// FineOperator &Fine,Matrix &FineMatrix, // FineOperator &Fine,Matrix &FineMatrix,
// FineOperator &Smooth,Matrix &SmootherMatrix) // FineOperator &Smooth,Matrix &SmootherMatrix)
@ -831,8 +823,8 @@ int main (int argc, char ** argv)
std::cout << GridLogMessage << "Building two level VPGCR and FGMRES solvers" << std::endl; std::cout << GridLogMessage << "Building two level VPGCR and FGMRES solvers" << std::endl;
std::cout << GridLogMessage << "**************************************************" << std::endl; std::cout << GridLogMessage << "**************************************************" << std::endl;
PrecGeneralisedConjugateResidual<LatticeFermion> VPGCRMG(1.0e-12,100,Precon,8,8); PrecGeneralisedConjugateResidual<LatticeFermion> VPGCRMG(1.0e-12, 100, Precon, 8, 8);
FlexibleGeneralisedMinimalResidual<LatticeFermion> FGMRESMG(1.0e-12,100,Precon,8); FlexibleGeneralisedMinimalResidual<LatticeFermion> FGMRESMG(1.0e-12, 100, Precon, 8);
std::cout << GridLogMessage << "checking norm src " << norm2(src) << std::endl; std::cout << GridLogMessage << "checking norm src " << norm2(src) << std::endl;
@ -840,14 +832,14 @@ int main (int argc, char ** argv)
std::cout << GridLogMessage << "Building unpreconditioned VPGCR and FGMRES solvers" << std::endl; std::cout << GridLogMessage << "Building unpreconditioned VPGCR and FGMRES solvers" << std::endl;
std::cout << GridLogMessage << "**************************************************" << std::endl; std::cout << GridLogMessage << "**************************************************" << std::endl;
PrecGeneralisedConjugateResidual<LatticeFermion> VPGCRT(1.0e-12,4000000,Simple,8,8); PrecGeneralisedConjugateResidual<LatticeFermion> VPGCRT(1.0e-12, 4000000, Simple, 8, 8);
FlexibleGeneralisedMinimalResidual<LatticeFermion> FGMREST(1.0e-12,4000000,Simple,8); FlexibleGeneralisedMinimalResidual<LatticeFermion> FGMREST(1.0e-12, 4000000, Simple, 8);
std::cout << GridLogMessage << "**************************************************" << std::endl; std::cout << GridLogMessage << "**************************************************" << std::endl;
std::cout << GridLogMessage << "Testing the four solvers" << std::endl; std::cout << GridLogMessage << "Testing the four solvers" << std::endl;
std::cout << GridLogMessage << "**************************************************" << std::endl; std::cout << GridLogMessage << "**************************************************" << std::endl;
std::vector< OperatorFunction<LatticeFermion>*> solvers; std::vector<OperatorFunction<LatticeFermion> *> solvers;
solvers.push_back(&VPGCRMG); solvers.push_back(&VPGCRMG);
solvers.push_back(&FGMRESMG); solvers.push_back(&FGMRESMG);
solvers.push_back(&VPGCRT); solvers.push_back(&VPGCRT);
@ -855,7 +847,7 @@ int main (int argc, char ** argv)
for(auto elem : solvers) { for(auto elem : solvers) {
result = zero; result = zero;
(*elem)(HermIndefOp,src,result); (*elem)(HermIndefOp, src, result);
} }
Grid_finalize(); Grid_finalize();