mirror of
https://github.com/paboyle/Grid.git
synced 2025-06-21 17:22:03 +01:00
Compare commits
3 Commits
37884d369f
...
feature/bo
| Author | SHA1 | Date | |
|---|---|---|---|
| 6815e138b4 | |||
| e29b97b3ea | |||
| ad2b699d2b |
@ -123,7 +123,7 @@ public:
|
||||
};
|
||||
|
||||
template<class Fobj,class CComplex,int nbasis>
|
||||
class Aggregation {
|
||||
class Aggregation {
|
||||
public:
|
||||
typedef iVector<CComplex,nbasis > siteVector;
|
||||
typedef Lattice<siteVector> CoarseVector;
|
||||
@ -158,20 +158,7 @@ public:
|
||||
blockPromote(CoarseVec,FineVec,subspace);
|
||||
}
|
||||
|
||||
virtual void CreateSubspaceRandom(GridParallelRNG &RNG) {
|
||||
int nn=nbasis;
|
||||
RealD scale;
|
||||
FineField noise(FineGrid);
|
||||
for(int b=0;b<nn;b++){
|
||||
subspace[b] = Zero();
|
||||
gaussian(RNG,noise);
|
||||
scale = std::pow(norm2(noise),-0.5);
|
||||
noise=noise*scale;
|
||||
subspace[b] = noise;
|
||||
}
|
||||
}
|
||||
virtual void CreateSubspace(GridParallelRNG &RNG,LinearOperatorBase<FineField> &hermop,int nn=nbasis)
|
||||
{
|
||||
virtual void CreateSubspace(GridParallelRNG &RNG,LinearOperatorBase<FineField> &hermop,int nn=nbasis) {
|
||||
|
||||
RealD scale;
|
||||
|
||||
@ -230,11 +217,6 @@ public:
|
||||
scale = std::pow(norm2(noise),-0.5);
|
||||
noise=noise*scale;
|
||||
|
||||
std::cout << GridLogMessage<<" Chebyshev subspace pass-1 : ord "<<orderfilter<<" ["<<lo<<","<<hi<<"]"<<std::endl;
|
||||
std::cout << GridLogMessage<<" Chebyshev subspace pass-2 : nbasis"<<nn<<" min "
|
||||
<<ordermin<<" step "<<orderstep
|
||||
<<" lo"<<filterlo<<std::endl;
|
||||
|
||||
// Initial matrix element
|
||||
hermop.Op(noise,Mn); std::cout<<GridLogMessage << "noise <n|MdagM|n> "<<norm2(Mn)<<std::endl;
|
||||
|
||||
@ -308,43 +290,6 @@ public:
|
||||
}
|
||||
assert(b==nn);
|
||||
}
|
||||
virtual void CreateSubspaceChebyshev(GridParallelRNG &RNG,LinearOperatorBase<FineField> &hermop,
|
||||
int nn,
|
||||
double hi,
|
||||
double lo,
|
||||
int orderfilter
|
||||
) {
|
||||
|
||||
RealD scale;
|
||||
|
||||
FineField noise(FineGrid);
|
||||
FineField Mn(FineGrid);
|
||||
FineField tmp(FineGrid);
|
||||
|
||||
// New normalised noise
|
||||
std::cout << GridLogMessage<<" Chebyshev subspace pure noise : ord "<<orderfilter<<" ["<<lo<<","<<hi<<"]"<<std::endl;
|
||||
std::cout << GridLogMessage<<" Chebyshev subspace pure noise : nbasis"<<nn<<std::endl;
|
||||
|
||||
|
||||
for(int b =0;b<nbasis;b++)
|
||||
{
|
||||
gaussian(RNG,noise);
|
||||
scale = std::pow(norm2(noise),-0.5);
|
||||
noise=noise*scale;
|
||||
|
||||
// Initial matrix element
|
||||
hermop.Op(noise,Mn); std::cout<<GridLogMessage << "noise <n|MdagM|n> "<<norm2(Mn)<<std::endl;
|
||||
// Filter
|
||||
Chebyshev<FineField> Cheb(lo,hi,orderfilter);
|
||||
Cheb(hermop,noise,Mn);
|
||||
// normalise
|
||||
scale = std::pow(norm2(Mn),-0.5); Mn=Mn*scale;
|
||||
subspace[b] = Mn;
|
||||
hermop.Op(Mn,tmp);
|
||||
std::cout<<GridLogMessage << "filt ["<<b<<"] <n|MdagM|n> "<<norm2(tmp)<<std::endl;
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
};
|
||||
|
||||
|
||||
@ -1,609 +0,0 @@
|
||||
/*************************************************************************************
|
||||
|
||||
Grid physics library, www.github.com/paboyle/Grid
|
||||
|
||||
Source file: ./lib/algorithms/GeneralCoarsenedMatrix.h
|
||||
|
||||
Copyright (C) 2015
|
||||
|
||||
Author: Peter Boyle <pboyle@bnl.gov>
|
||||
|
||||
This program is free software; you can redistribute it and/or modify
|
||||
it under the terms of the GNU General Public License as published by
|
||||
the Free Software Foundation; either version 2 of the License, or
|
||||
(at your option) any later version.
|
||||
|
||||
This program is distributed in the hope that it will be useful,
|
||||
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||
GNU General Public License for more details.
|
||||
|
||||
You should have received a copy of the GNU General Public License along
|
||||
with this program; if not, write to the Free Software Foundation, Inc.,
|
||||
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
|
||||
|
||||
See the full license in the file "LICENSE" in the top level distribution directory
|
||||
*************************************************************************************/
|
||||
/* END LEGAL */
|
||||
#pragma once
|
||||
|
||||
#include <Grid/qcd/QCD.h> // needed for Dagger(Yes|No), Inverse(Yes|No)
|
||||
|
||||
#include <Grid/lattice/PaddedCell.h>
|
||||
#include <Grid/stencil/GeneralLocalStencil.h>
|
||||
|
||||
NAMESPACE_BEGIN(Grid);
|
||||
|
||||
// Fixme need coalesced read gpermute
|
||||
template<class vobj> void gpermute(vobj & inout,int perm){
|
||||
vobj tmp=inout;
|
||||
if (perm & 0x1 ) { permute(inout,tmp,0); tmp=inout;}
|
||||
if (perm & 0x2 ) { permute(inout,tmp,1); tmp=inout;}
|
||||
if (perm & 0x4 ) { permute(inout,tmp,2); tmp=inout;}
|
||||
if (perm & 0x8 ) { permute(inout,tmp,3); tmp=inout;}
|
||||
}
|
||||
|
||||
/////////////////////////////////////////////////////////////////
|
||||
// Reuse Aggregation class from CoarsenedMatrix for now
|
||||
// Might think about *smoothed* Aggregation
|
||||
// Equivalent of Geometry class in cartesian case
|
||||
/////////////////////////////////////////////////////////////////
|
||||
class NonLocalStencilGeometry {
|
||||
public:
|
||||
int depth;
|
||||
int hops;
|
||||
int npoint;
|
||||
std::vector<Coordinate> shifts;
|
||||
Coordinate stencil_size;
|
||||
Coordinate stencil_lo;
|
||||
Coordinate stencil_hi;
|
||||
GridCartesian *grid;
|
||||
GridCartesian *Grid() {return grid;};
|
||||
int Depth(void){return 1;}; // Ghost zone depth
|
||||
int Hops(void){return hops;}; // # of hops=> level of corner fill in in stencil
|
||||
|
||||
virtual int DimSkip(void) =0;
|
||||
|
||||
virtual ~NonLocalStencilGeometry() {};
|
||||
|
||||
int Reverse(int point)
|
||||
{
|
||||
int Nd = Grid()->Nd();
|
||||
Coordinate shft = shifts[point];
|
||||
Coordinate rev(Nd);
|
||||
for(int mu=0;mu<Nd;mu++) rev[mu]= -shft[mu];
|
||||
for(int p=0;p<npoint;p++){
|
||||
if(rev==shifts[p]){
|
||||
return p;
|
||||
}
|
||||
}
|
||||
assert(0);
|
||||
return -1;
|
||||
}
|
||||
void BuildShifts(void)
|
||||
{
|
||||
this->shifts.resize(0);
|
||||
int Nd = this->grid->Nd();
|
||||
|
||||
int dd = this->DimSkip();
|
||||
for(int s0=this->stencil_lo[dd+0];s0<=this->stencil_hi[dd+0];s0++){
|
||||
for(int s1=this->stencil_lo[dd+1];s1<=this->stencil_hi[dd+1];s1++){
|
||||
for(int s2=this->stencil_lo[dd+2];s2<=this->stencil_hi[dd+2];s2++){
|
||||
for(int s3=this->stencil_lo[dd+3];s3<=this->stencil_hi[dd+3];s3++){
|
||||
Coordinate sft(Nd,0);
|
||||
sft[dd+0] = s0;
|
||||
sft[dd+1] = s1;
|
||||
sft[dd+2] = s2;
|
||||
sft[dd+3] = s3;
|
||||
int nhops = abs(s0)+abs(s1)+abs(s2)+abs(s3);
|
||||
if(nhops<=this->hops) this->shifts.push_back(sft);
|
||||
}}}}
|
||||
this->npoint = this->shifts.size();
|
||||
std::cout << GridLogMessage << "NonLocalStencilGeometry has "<< this->npoint << " terms in stencil "<<std::endl;
|
||||
}
|
||||
|
||||
NonLocalStencilGeometry(GridCartesian *_coarse_grid,int _hops) : grid(_coarse_grid), hops(_hops)
|
||||
{
|
||||
Coordinate latt = grid->GlobalDimensions();
|
||||
stencil_size.resize(grid->Nd());
|
||||
stencil_lo.resize(grid->Nd());
|
||||
stencil_hi.resize(grid->Nd());
|
||||
for(int d=0;d<grid->Nd();d++){
|
||||
if ( latt[d] == 1 ) {
|
||||
stencil_lo[d] = 0;
|
||||
stencil_hi[d] = 0;
|
||||
stencil_size[d]= 1;
|
||||
} else if ( latt[d] == 2 ) {
|
||||
stencil_lo[d] = -1;
|
||||
stencil_hi[d] = 0;
|
||||
stencil_size[d]= 2;
|
||||
} else if ( latt[d] > 2 ) {
|
||||
stencil_lo[d] = -1;
|
||||
stencil_hi[d] = 1;
|
||||
stencil_size[d]= 3;
|
||||
}
|
||||
}
|
||||
};
|
||||
|
||||
};
|
||||
|
||||
// Need to worry about red-black now
|
||||
class NonLocalStencilGeometry4D : public NonLocalStencilGeometry {
|
||||
public:
|
||||
virtual int DimSkip(void) { return 0;};
|
||||
NonLocalStencilGeometry4D(GridCartesian *Coarse,int _hops) : NonLocalStencilGeometry(Coarse,_hops) { };
|
||||
virtual ~NonLocalStencilGeometry4D() {};
|
||||
};
|
||||
class NonLocalStencilGeometry5D : public NonLocalStencilGeometry {
|
||||
public:
|
||||
virtual int DimSkip(void) { return 1; };
|
||||
NonLocalStencilGeometry5D(GridCartesian *Coarse,int _hops) : NonLocalStencilGeometry(Coarse,_hops) { };
|
||||
virtual ~NonLocalStencilGeometry5D() {};
|
||||
};
|
||||
/*
|
||||
* Bunch of different options classes
|
||||
*/
|
||||
class NextToNextToNextToNearestStencilGeometry4D : public NonLocalStencilGeometry4D {
|
||||
public:
|
||||
NextToNextToNextToNearestStencilGeometry4D(GridCartesian *Coarse) : NonLocalStencilGeometry4D(Coarse,4)
|
||||
{
|
||||
this->BuildShifts();
|
||||
};
|
||||
};
|
||||
class NextToNextToNextToNearestStencilGeometry5D : public NonLocalStencilGeometry5D {
|
||||
public:
|
||||
NextToNextToNextToNearestStencilGeometry5D(GridCartesian *Coarse) : NonLocalStencilGeometry5D(Coarse,4)
|
||||
{
|
||||
this->BuildShifts();
|
||||
};
|
||||
};
|
||||
class NextToNearestStencilGeometry4D : public NonLocalStencilGeometry4D {
|
||||
public:
|
||||
NextToNearestStencilGeometry4D(GridCartesian *Coarse) : NonLocalStencilGeometry4D(Coarse,2)
|
||||
{
|
||||
this->BuildShifts();
|
||||
};
|
||||
};
|
||||
class NextToNearestStencilGeometry5D : public NonLocalStencilGeometry5D {
|
||||
public:
|
||||
NextToNearestStencilGeometry5D(GridCartesian *Coarse) : NonLocalStencilGeometry5D(Coarse,2)
|
||||
{
|
||||
this->BuildShifts();
|
||||
};
|
||||
};
|
||||
class NearestStencilGeometry4D : public NonLocalStencilGeometry4D {
|
||||
public:
|
||||
NearestStencilGeometry4D(GridCartesian *Coarse) : NonLocalStencilGeometry4D(Coarse,1)
|
||||
{
|
||||
this->BuildShifts();
|
||||
};
|
||||
};
|
||||
class NearestStencilGeometry5D : public NonLocalStencilGeometry5D {
|
||||
public:
|
||||
NearestStencilGeometry5D(GridCartesian *Coarse) : NonLocalStencilGeometry5D(Coarse,1)
|
||||
{
|
||||
this->BuildShifts();
|
||||
};
|
||||
};
|
||||
|
||||
// Fine Object == (per site) type of fine field
|
||||
// nbasis == number of deflation vectors
|
||||
template<class Fobj,class CComplex,int nbasis>
|
||||
class GeneralCoarsenedMatrix : public SparseMatrixBase<Lattice<iVector<CComplex,nbasis > > > {
|
||||
public:
|
||||
|
||||
typedef iVector<CComplex,nbasis > siteVector;
|
||||
typedef Lattice<iScalar<CComplex> > CoarseComplexField;
|
||||
typedef Lattice<siteVector> CoarseVector;
|
||||
typedef Lattice<iMatrix<CComplex,nbasis > > CoarseMatrix;
|
||||
typedef iMatrix<CComplex,nbasis > Cobj;
|
||||
typedef Lattice< CComplex > CoarseScalar; // used for inner products on fine field
|
||||
typedef Lattice<Fobj > FineField;
|
||||
typedef CoarseVector Field;
|
||||
////////////////////
|
||||
// Data members
|
||||
////////////////////
|
||||
int hermitian;
|
||||
GridBase * _FineGrid;
|
||||
GridCartesian * _CoarseGrid;
|
||||
NonLocalStencilGeometry &geom;
|
||||
PaddedCell Cell;
|
||||
GeneralLocalStencil Stencil;
|
||||
|
||||
std::vector<CoarseMatrix> _A;
|
||||
std::vector<CoarseMatrix> _Adag;
|
||||
|
||||
///////////////////////
|
||||
// Interface
|
||||
///////////////////////
|
||||
GridBase * Grid(void) { return _FineGrid; }; // this is all the linalg routines need to know
|
||||
GridBase * FineGrid(void) { return _FineGrid; }; // this is all the linalg routines need to know
|
||||
GridCartesian * CoarseGrid(void) { return _CoarseGrid; }; // this is all the linalg routines need to know
|
||||
|
||||
GeneralCoarsenedMatrix(NonLocalStencilGeometry &_geom,GridBase *FineGrid, GridCartesian * CoarseGrid)
|
||||
: geom(_geom),
|
||||
_FineGrid(FineGrid),
|
||||
_CoarseGrid(CoarseGrid),
|
||||
hermitian(1),
|
||||
Cell(_geom.Depth(),_CoarseGrid),
|
||||
Stencil(Cell.grids.back(),geom.shifts)
|
||||
{
|
||||
{
|
||||
int npoint = _geom.npoint;
|
||||
StencilEntry *SE;
|
||||
autoView( Stencil_v , Stencil, AcceleratorRead);
|
||||
int osites=Stencil.Grid()->oSites();
|
||||
for(int ss=0;ss<osites;ss++){
|
||||
for(int point=0;point<npoint;point++){
|
||||
auto SE = Stencil_v.GetEntry(point,ss);
|
||||
int o = SE->_offset;
|
||||
assert( o< osites);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
_A.resize(geom.npoint,CoarseGrid);
|
||||
_Adag.resize(geom.npoint,CoarseGrid);
|
||||
}
|
||||
void M (const CoarseVector &in, CoarseVector &out)
|
||||
{
|
||||
Mult(_A,in,out);
|
||||
}
|
||||
void Mdag (const CoarseVector &in, CoarseVector &out)
|
||||
{
|
||||
Mult(_Adag,in,out);
|
||||
}
|
||||
void Mult (std::vector<CoarseMatrix> &A,const CoarseVector &in, CoarseVector &out)
|
||||
{
|
||||
conformable(CoarseGrid(),in.Grid());
|
||||
conformable(in.Grid(),out.Grid());
|
||||
out.Checkerboard() = in.Checkerboard();
|
||||
CoarseVector tin=in;
|
||||
|
||||
CoarseVector pin = Cell.Exchange(tin);
|
||||
|
||||
CoarseVector pout(pin.Grid());
|
||||
|
||||
autoView( in_v , pin, AcceleratorRead);
|
||||
autoView( out_v , pout, AcceleratorWrite);
|
||||
autoView( Stencil_v , Stencil, AcceleratorRead);
|
||||
int npoint = geom.npoint;
|
||||
typedef LatticeView<Cobj> Aview;
|
||||
|
||||
Vector<Aview> AcceleratorViewContainer;
|
||||
|
||||
for(int p=0;p<npoint;p++) AcceleratorViewContainer.push_back(A[p].View(AcceleratorRead));
|
||||
|
||||
Aview *Aview_p = & AcceleratorViewContainer[0];
|
||||
|
||||
const int Nsimd = CComplex::Nsimd();
|
||||
typedef siteVector calcVector;
|
||||
typedef CComplex calcComplex;
|
||||
|
||||
int osites=pin.Grid()->oSites();
|
||||
|
||||
for(int point=0;point<npoint;point++){
|
||||
conformable(A[point],pin);
|
||||
}
|
||||
|
||||
accelerator_for(sss, osites*nbasis, 1, {
|
||||
int ss = sss/nbasis;
|
||||
int b = sss%nbasis;
|
||||
assert(ss<osites);
|
||||
calcComplex res;
|
||||
res = Zero();
|
||||
calcVector nbr;
|
||||
int ptype;
|
||||
StencilEntry *SE;
|
||||
|
||||
for(int point=0;point<npoint;point++){
|
||||
|
||||
auto SE = Stencil_v.GetEntry(point,ss);
|
||||
|
||||
int o = SE->_offset;
|
||||
|
||||
assert( o < osites);
|
||||
// gpermute etc..
|
||||
nbr = in_v[o];
|
||||
gpermute(nbr,SE->_permute);
|
||||
|
||||
for(int bb=0;bb<nbasis;bb++) {
|
||||
res = res + Aview_p[point][ss](b,bb)*nbr(bb);
|
||||
}
|
||||
}
|
||||
out_v[ss](b)=res;
|
||||
});
|
||||
|
||||
for(int p=0;p<geom.npoint;p++) AcceleratorViewContainer[p].ViewClose();
|
||||
|
||||
out = Cell.Extract(pout);
|
||||
};
|
||||
|
||||
void PopulateAdag(void)
|
||||
{
|
||||
for(int bidx=0;bidx<CoarseGrid()->gSites() ;bidx++){
|
||||
Coordinate bcoor;
|
||||
CoarseGrid()->GlobalIndexToGlobalCoor(bidx,bcoor);
|
||||
|
||||
for(int p=0;p<geom.npoint;p++){
|
||||
Coordinate scoor = bcoor;
|
||||
for(int mu=0;mu<bcoor.size();mu++){
|
||||
int L = CoarseGrid()->GlobalDimensions()[mu];
|
||||
scoor[mu] = (bcoor[mu] - geom.shifts[p][mu] + L) % L; // Modulo arithmetic
|
||||
}
|
||||
// Flip to poke/peekLocalSite and not too bad
|
||||
auto link = peekSite(_A[p],scoor);
|
||||
int pp = geom.Reverse(p);
|
||||
pokeSite(adj(link),_Adag[pp],bcoor);
|
||||
}
|
||||
}
|
||||
}
|
||||
void CoarsenOperator(LinearOperatorBase<Lattice<Fobj> > &linop,
|
||||
Aggregation<Fobj,CComplex,nbasis> & Subspace)
|
||||
{
|
||||
RealD tproj=0.0;
|
||||
RealD tpick=0.0;
|
||||
RealD tmat=0.0;
|
||||
RealD tpeek=0.0;
|
||||
std::cout << GridLogMessage<< "CoarsenMatrix "<< std::endl;
|
||||
GridBase *grid = FineGrid();
|
||||
|
||||
////////////////////////////////////////////////
|
||||
// Orthogonalise the subblocks over the basis
|
||||
////////////////////////////////////////////////
|
||||
CoarseScalar InnerProd(CoarseGrid());
|
||||
blockOrthogonalise(InnerProd,Subspace.subspace);
|
||||
|
||||
////////////////////////////////////////////////
|
||||
// Now compute the matrix elements of linop between this orthonormal
|
||||
// set of vectors.
|
||||
////////////////////////////////////////////////
|
||||
FineField bV(grid);
|
||||
FineField MbV(grid);
|
||||
FineField tmp(grid);
|
||||
CoarseVector coarseInner(CoarseGrid());
|
||||
|
||||
// Very inefficient loop of order coarse volume.
|
||||
// First pass hack
|
||||
// Could replace with a coloring scheme our phase scheme
|
||||
// as in BFM
|
||||
for(int bidx=0;bidx<CoarseGrid()->gSites() ;bidx++){
|
||||
Coordinate bcoor;
|
||||
CoarseGrid()->GlobalIndexToGlobalCoor(bidx,bcoor);
|
||||
|
||||
for(int b=0;b<nbasis;b++){
|
||||
tpick-=usecond();
|
||||
blockPick(CoarseGrid(),Subspace.subspace[b],bV,bcoor);
|
||||
tpick+=usecond();
|
||||
tmat-=usecond();
|
||||
linop.Op(bV,MbV);
|
||||
tmat+=usecond();
|
||||
tproj-=usecond();
|
||||
blockProject(coarseInner,MbV,Subspace.subspace);
|
||||
tproj+=usecond();
|
||||
|
||||
tpeek-=usecond();
|
||||
for(int p=0;p<geom.npoint;p++){
|
||||
Coordinate scoor = bcoor;
|
||||
for(int mu=0;mu<bcoor.size();mu++){
|
||||
int L = CoarseGrid()->GlobalDimensions()[mu];
|
||||
scoor[mu] = (bcoor[mu] - geom.shifts[p][mu] + L) % L; // Modulo arithmetic
|
||||
}
|
||||
// Flip to peekLocalSite
|
||||
// Flip to pokeLocalSite
|
||||
auto ip = peekSite(coarseInner,scoor);
|
||||
auto Ab = peekSite(_A[p],scoor);
|
||||
int pp = geom.Reverse(p);
|
||||
auto Adagb = peekSite(_Adag[pp],bcoor);
|
||||
for(int bb=0;bb<nbasis;bb++){
|
||||
Ab(bb,b) = ip(bb);
|
||||
Adagb(b,bb) = conjugate(ip(bb));
|
||||
}
|
||||
pokeSite(Ab,_A[p],scoor);
|
||||
pokeSite(Adagb,_Adag[pp],bcoor);
|
||||
}
|
||||
tpeek+=usecond();
|
||||
}
|
||||
}
|
||||
for(int p=0;p<geom.npoint;p++){
|
||||
Coordinate coor({0,0,0,0,0});
|
||||
auto sval = peekSite(_A[p],coor);
|
||||
}
|
||||
for(int p=0;p<geom.npoint;p++){
|
||||
_A[p] = Cell.Exchange(_A[p]);
|
||||
_Adag[p]= Cell.Exchange(_Adag[p]);
|
||||
}
|
||||
std::cout << GridLogMessage<<"CoarsenOperator pick "<<tpick<<" us"<<std::endl;
|
||||
std::cout << GridLogMessage<<"CoarsenOperator mat "<<tmat <<" us"<<std::endl;
|
||||
std::cout << GridLogMessage<<"CoarsenOperator projection "<<tproj<<" us"<<std::endl;
|
||||
std::cout << GridLogMessage<<"CoarsenOperator peek/poke "<<tpeek<<" us"<<std::endl;
|
||||
}
|
||||
|
||||
/////////////////////////////////////////////////////////////
|
||||
//
|
||||
// A) Only reduced flops option is to use a padded cell of depth 4
|
||||
// and apply MpcDagMpc in the padded cell.
|
||||
//
|
||||
// Makes for ONE application of MpcDagMpc per vector instead of 30 or 80.
|
||||
// With the effective cell size around (B+8)^4 perhaps 12^4/4^4 ratio
|
||||
// Cost is 81x more, same as stencil size.
|
||||
//
|
||||
// But: can eliminate comms and do as local dirichlet.
|
||||
//
|
||||
// Local exchange gauge field once.
|
||||
// Apply to all vectors, local only computation.
|
||||
// Must exchange ghost subcells in reverse process of PaddedCell to take inner products
|
||||
//
|
||||
// B) Can reduce cost: pad by 1, apply Deo (4^4+6^4+8^4+8^4 )/ (4x 4^4)
|
||||
// pad by 2, apply Doe
|
||||
// pad by 3, apply Deo
|
||||
// then break out 8x directions; cost is ~10x MpcDagMpc per vector
|
||||
//
|
||||
// => almost factor of 10 in setup cost, excluding data rearrangement
|
||||
//
|
||||
// Intermediates -- ignore the corner terms, leave approximate and force Hermitian
|
||||
// Intermediates -- pad by 2 and apply 1+8+24 = 33 times.
|
||||
/////////////////////////////////////////////////////////////
|
||||
|
||||
//////////////////////////////////////////////////////////
|
||||
// BFM HDCG style approach: Solve a system of equations to get Aij
|
||||
//////////////////////////////////////////////////////////
|
||||
/*
|
||||
* Here, k,l index which possible shift within the 3^Nd "ball" connected by MdagM.
|
||||
*
|
||||
* conj(phases[block]) proj[k][ block*Nvec+j ] = \sum_ball e^{i q_k . delta} < phi_{block,j} | MdagM | phi_{(block+delta),i} >
|
||||
* = \sum_ball e^{iqk.delta} A_ji
|
||||
*
|
||||
* Must invert matrix M_k,l = e^[i q_k . delta_l]
|
||||
*
|
||||
* Where q_k = delta_k . (2*M_PI/global_nb[mu])
|
||||
*/
|
||||
void CoarsenOperatorColoured(LinearOperatorBase<Lattice<Fobj> > &linop,
|
||||
Aggregation<Fobj,CComplex,nbasis> & Subspace)
|
||||
{
|
||||
std::cout << GridLogMessage<< "CoarsenMatrixColoured "<< std::endl;
|
||||
GridBase *grid = FineGrid();
|
||||
|
||||
RealD tproj=0.0;
|
||||
RealD teigen=0.0;
|
||||
RealD tmat=0.0;
|
||||
RealD tphase=0.0;
|
||||
RealD tinv=0.0;
|
||||
|
||||
/////////////////////////////////////////////////////////////
|
||||
// Orthogonalise the subblocks over the basis
|
||||
/////////////////////////////////////////////////////////////
|
||||
CoarseScalar InnerProd(CoarseGrid());
|
||||
blockOrthogonalise(InnerProd,Subspace.subspace);
|
||||
|
||||
const int npoint = geom.npoint;
|
||||
|
||||
Coordinate clatt = CoarseGrid()->GlobalDimensions();
|
||||
int Nd = CoarseGrid()->Nd();
|
||||
|
||||
/*
|
||||
* Here, k,l index which possible momentum/shift within the N-points connected by MdagM.
|
||||
* Matrix index i is mapped to this shift via
|
||||
* geom.shifts[i]
|
||||
*
|
||||
* conj(pha[block]) proj[k (which mom)][j (basis vec cpt)][block]
|
||||
* = \sum_{l in ball} e^{i q_k . delta_l} < phi_{block,j} | MdagM | phi_{(block+delta_l),i} >
|
||||
* = \sum_{l in ball} e^{iqk.delta_l} A_ji^{b.b+l}
|
||||
* = M_{kl} A_ji^{b.b+l}
|
||||
*
|
||||
* Must assemble and invert matrix M_k,l = e^[i q_k . delta_l]
|
||||
*
|
||||
* Where q_k = delta_k . (2*M_PI/global_nb[mu])
|
||||
*
|
||||
* Then A{ji}^{b,b+l} = M^{-1}_{lm} ComputeProj_{m,b,i,j}
|
||||
*/
|
||||
teigen-=usecond();
|
||||
ComplexD ci(0.0,1.0);
|
||||
Eigen::MatrixXcd Mkl = Eigen::MatrixXcd::Zero(npoint,npoint);
|
||||
Eigen::MatrixXcd invMkl = Eigen::MatrixXcd::Zero(npoint,npoint);
|
||||
for(int k=0;k<npoint;k++){ // Loop over momenta
|
||||
|
||||
for(int l=0;l<npoint;l++){ // Loop over nbr relative
|
||||
std::complex<double> phase(0.0,0.0);
|
||||
for(int mu=0;mu<Nd;mu++){
|
||||
RealD TwoPiL = M_PI * 2.0/ clatt[mu];
|
||||
phase=phase+TwoPiL*geom.shifts[k][mu]*geom.shifts[l][mu];
|
||||
}
|
||||
phase=exp(phase*ci);
|
||||
Mkl(k,l) = phase;
|
||||
}
|
||||
}
|
||||
invMkl = Mkl.inverse();
|
||||
teigen+=usecond();
|
||||
|
||||
///////////////////////////////////////////////////////////////////////
|
||||
// Now compute the matrix elements of linop between the orthonormal
|
||||
// set of vectors.
|
||||
///////////////////////////////////////////////////////////////////////
|
||||
FineField phaV(grid); // Phased block basis vector
|
||||
FineField MphaV(grid);// Matrix applied
|
||||
CoarseVector coarseInner(CoarseGrid());
|
||||
|
||||
std::vector<CoarseVector> ComputeProj(npoint,CoarseGrid());
|
||||
std::vector<CoarseVector> FT(npoint,CoarseGrid());
|
||||
for(int i=0;i<nbasis;i++){// Loop over basis vectors
|
||||
std::cout << GridLogMessage<< "CoarsenMatrixColoured vec "<<i<<"/"<<nbasis<< std::endl;
|
||||
for(int p=0;p<npoint;p++){ // Loop over momenta in npoint
|
||||
/////////////////////////////////////////////////////
|
||||
// Stick a phase on every block
|
||||
/////////////////////////////////////////////////////
|
||||
tphase-=usecond();
|
||||
CoarseComplexField coor(CoarseGrid());
|
||||
CoarseComplexField pha(CoarseGrid()); pha=Zero();
|
||||
for(int mu=0;mu<Nd;mu++){
|
||||
LatticeCoordinate(coor,mu);
|
||||
RealD TwoPiL = M_PI * 2.0/ clatt[mu];
|
||||
pha = pha + (TwoPiL * geom.shifts[p][mu]) * coor;
|
||||
}
|
||||
pha =exp(pha*ci);
|
||||
phaV=Zero();
|
||||
blockZAXPY(phaV,pha,Subspace.subspace[i],phaV);
|
||||
tphase+=usecond();
|
||||
|
||||
/////////////////////////////////////////////////////////////////////
|
||||
// Multiple phased subspace vector by matrix and project to subspace
|
||||
// Remove local bulk phase to leave relative phases
|
||||
/////////////////////////////////////////////////////////////////////
|
||||
tmat-=usecond();
|
||||
linop.Op(phaV,MphaV);
|
||||
tmat+=usecond();
|
||||
|
||||
tproj-=usecond();
|
||||
blockProject(coarseInner,MphaV,Subspace.subspace);
|
||||
coarseInner = conjugate(pha) * coarseInner;
|
||||
|
||||
ComputeProj[p] = coarseInner;
|
||||
tproj+=usecond();
|
||||
|
||||
}
|
||||
|
||||
tinv-=usecond();
|
||||
for(int k=0;k<npoint;k++){
|
||||
FT[k] = Zero();
|
||||
for(int l=0;l<npoint;l++){
|
||||
FT[k]= FT[k]+ invMkl(l,k)*ComputeProj[l];
|
||||
}
|
||||
|
||||
int osites=CoarseGrid()->oSites();
|
||||
autoView( A_v , _A[k], AcceleratorWrite);
|
||||
autoView( FT_v , FT[k], AcceleratorRead);
|
||||
accelerator_for(sss, osites, 1, {
|
||||
for(int j=0;j<nbasis;j++){
|
||||
A_v[sss](j,i) = FT_v[sss](j);
|
||||
}
|
||||
});
|
||||
}
|
||||
tinv+=usecond();
|
||||
}
|
||||
|
||||
for(int p=0;p<geom.npoint;p++){
|
||||
Coordinate coor({0,0,0,0,0});
|
||||
auto sval = peekSite(_A[p],coor);
|
||||
}
|
||||
|
||||
PopulateAdag();
|
||||
|
||||
// Need to write something to populate Adag from A
|
||||
for(int p=0;p<geom.npoint;p++){
|
||||
_A[p] = Cell.Exchange(_A[p]);
|
||||
_Adag[p]= Cell.Exchange(_Adag[p]);
|
||||
}
|
||||
std::cout << GridLogMessage<<"CoarsenOperator eigen "<<teigen<<" us"<<std::endl;
|
||||
std::cout << GridLogMessage<<"CoarsenOperator phase "<<tphase<<" us"<<std::endl;
|
||||
std::cout << GridLogMessage<<"CoarsenOperator mat "<<tmat <<" us"<<std::endl;
|
||||
std::cout << GridLogMessage<<"CoarsenOperator proj "<<tproj<<" us"<<std::endl;
|
||||
std::cout << GridLogMessage<<"CoarsenOperator inv "<<tinv<<" us"<<std::endl;
|
||||
}
|
||||
|
||||
virtual void Mdiag (const Field &in, Field &out){ assert(0);};
|
||||
virtual void Mdir (const Field &in, Field &out,int dir, int disp){assert(0);};
|
||||
virtual void MdirAll (const Field &in, std::vector<Field> &out){assert(0);};
|
||||
};
|
||||
|
||||
|
||||
NAMESPACE_END(Grid);
|
||||
@ -90,8 +90,9 @@ public:
|
||||
order=_order;
|
||||
|
||||
if(order < 2) exit(-1);
|
||||
Coeffs.resize(order,0.0);
|
||||
Coeffs[order-1] = 1.0;
|
||||
Coeffs.resize(order);
|
||||
Coeffs.assign(0.,order);
|
||||
Coeffs[order-1] = 1.;
|
||||
};
|
||||
|
||||
// PB - more efficient low pass drops high modes above the low as 1/x uses all Chebyshev's.
|
||||
|
||||
@ -33,6 +33,15 @@ Author: Peter Boyle <paboyle@ph.ed.ac.uk>
|
||||
* Script A = SolverMatrix
|
||||
* Script P = Preconditioner
|
||||
*
|
||||
* Deflation methods considered
|
||||
* -- Solve P A x = P b [ like Luscher ]
|
||||
* DEF-1 M P A x = M P b [i.e. left precon]
|
||||
* DEF-2 P^T M A x = P^T M b
|
||||
* ADEF-1 Preconditioner = M P + Q [ Q + M + M A Q]
|
||||
* ADEF-2 Preconditioner = P^T M + Q
|
||||
* BNN Preconditioner = P^T M P + Q
|
||||
* BNN2 Preconditioner = M P + P^TM +Q - M P A M
|
||||
*
|
||||
* Implement ADEF-2
|
||||
*
|
||||
* Vstart = P^Tx + Qb
|
||||
@ -40,157 +49,45 @@ Author: Peter Boyle <paboyle@ph.ed.ac.uk>
|
||||
* M2=M3=1
|
||||
* Vout = x
|
||||
*/
|
||||
NAMESPACE_BEGIN(Grid);
|
||||
|
||||
template<class Field, class CoarseField, class Aggregation>
|
||||
// abstract base
|
||||
template<class Field, class CoarseField>
|
||||
class TwoLevelFlexiblePcg : public LinearFunction<Field>
|
||||
{
|
||||
public:
|
||||
int verbose;
|
||||
RealD Tolerance;
|
||||
Integer MaxIterations;
|
||||
const int mmax = 1;
|
||||
const int mmax = 5;
|
||||
GridBase *grid;
|
||||
GridBase *coarsegrid;
|
||||
|
||||
// Fine operator, Smoother, CoarseSolver
|
||||
LinearOperatorBase<Field> &_FineLinop;
|
||||
LinearFunction<Field> &_Smoother;
|
||||
LinearFunction<CoarseField> &_CoarseSolver;
|
||||
LinearFunction<CoarseField> &_CoarseSolverPrecise;
|
||||
LinearOperatorBase<Field> *_Linop
|
||||
OperatorFunction<Field> *_Smoother,
|
||||
LinearFunction<CoarseField> *_CoarseSolver;
|
||||
|
||||
// Need something that knows how to get from Coarse to fine and back again
|
||||
// void ProjectToSubspace(CoarseVector &CoarseVec,const FineField &FineVec){
|
||||
// void PromoteFromSubspace(const CoarseVector &CoarseVec,FineField &FineVec){
|
||||
Aggregation &_Aggregates;
|
||||
// Need somthing that knows how to get from Coarse to fine and back again
|
||||
|
||||
// more most opertor functions
|
||||
TwoLevelFlexiblePcg(RealD tol,
|
||||
Integer maxit,
|
||||
LinearOperatorBase<Field> &FineLinop,
|
||||
LinearFunction<Field> &Smoother,
|
||||
LinearFunction<CoarseField> &CoarseSolver,
|
||||
LinearFunction<CoarseField> &CoarseSolverPrecise,
|
||||
Aggregation &Aggregates
|
||||
Integer maxit,
|
||||
LinearOperatorBase<Field> *Linop,
|
||||
LinearOperatorBase<Field> *SmootherLinop,
|
||||
OperatorFunction<Field> *Smoother,
|
||||
OperatorFunction<CoarseField> CoarseLinop
|
||||
) :
|
||||
Tolerance(tol),
|
||||
MaxIterations(maxit),
|
||||
_FineLinop(FineLinop),
|
||||
_Smoother(Smoother),
|
||||
_CoarseSolver(CoarseSolver),
|
||||
_CoarseSolverPrecise(CoarseSolverPrecise),
|
||||
_Aggregates(Aggregates)
|
||||
{
|
||||
coarsegrid = Aggregates.CoarseGrid;
|
||||
grid = Aggregates.FineGrid;
|
||||
_Linop(Linop),
|
||||
_PreconditionerLinop(PrecLinop),
|
||||
_Preconditioner(Preconditioner)
|
||||
{
|
||||
verbose=0;
|
||||
};
|
||||
|
||||
void Inflexible(Field &src,Field &psi)
|
||||
{
|
||||
Field resid(grid);
|
||||
RealD f;
|
||||
RealD rtzp,rtz,a,d,b;
|
||||
RealD rptzp;
|
||||
|
||||
Field x(grid);
|
||||
Field p(grid);
|
||||
Field z(grid);
|
||||
Field tmp(grid);
|
||||
Field mmp(grid);
|
||||
Field r (grid);
|
||||
Field mu (grid);
|
||||
Field rp (grid);
|
||||
|
||||
//Initial residual computation & set up
|
||||
RealD guess = norm2(psi);
|
||||
double tn;
|
||||
|
||||
//////////////////////////
|
||||
// x0 = Vstart -- possibly modify guess
|
||||
//////////////////////////
|
||||
x=Zero();
|
||||
Vstart(x,src);
|
||||
|
||||
// r0 = b -A x0
|
||||
_FineLinop.HermOp(x,mmp);
|
||||
|
||||
axpy(r, -1.0, mmp, src); // Recomputes r=src-x0
|
||||
rp=r;
|
||||
|
||||
//////////////////////////////////
|
||||
// Compute z = M1 x
|
||||
//////////////////////////////////
|
||||
PcgM1(r,z);
|
||||
rtzp =real(innerProduct(r,z));
|
||||
|
||||
///////////////////////////////////////
|
||||
// Except Def2, M2 is trivial
|
||||
///////////////////////////////////////
|
||||
p=z;
|
||||
|
||||
RealD ssq = norm2(src);
|
||||
RealD rsq = ssq*Tolerance*Tolerance;
|
||||
|
||||
std::cout<<GridLogMessage<<"HDCG: k=0 residual "<<rtzp<<" target rsq "<<rsq<<" ssq "<<ssq<<std::endl;
|
||||
|
||||
for (int k=1;k<=MaxIterations;k++){
|
||||
|
||||
rtz=rtzp;
|
||||
d= PcgM3(p,mmp);
|
||||
a = rtz/d;
|
||||
|
||||
axpy(x,a,p,x);
|
||||
RealD rn = axpy_norm(r,-a,mmp,r);
|
||||
|
||||
PcgM1(r,z);
|
||||
|
||||
rtzp =real(innerProduct(r,z));
|
||||
|
||||
int ipcg=1; // almost free inexact preconditioned CG
|
||||
if (ipcg) {
|
||||
rptzp =real(innerProduct(rp,z));
|
||||
} else {
|
||||
rptzp =0;
|
||||
}
|
||||
b = (rtzp-rptzp)/rtz;
|
||||
|
||||
PcgM2(z,mu); // ADEF-2 this is identity. Axpy possible to eliminate
|
||||
|
||||
axpy(p,b,p,mu); // mu = A r
|
||||
|
||||
RealD rrn=sqrt(rn/ssq);
|
||||
RealD rtn=sqrt(rtz/ssq);
|
||||
std::cout<<GridLogMessage<<"HDCG: Pcg k= "<<k<<" residual = "<<rrn<<std::endl;
|
||||
|
||||
if ( ipcg ) {
|
||||
axpy(rp,0.0,r,r);
|
||||
}
|
||||
|
||||
// Stopping condition
|
||||
if ( rn <= rsq ) {
|
||||
|
||||
std::cout<<GridLogMessage<<"HDCG: Pcg converged in "<<k<<" iterations"<<std::endl;;
|
||||
|
||||
_FineLinop.HermOp(x,mmp);
|
||||
axpy(tmp,-1.0,src,mmp);
|
||||
|
||||
RealD mmpnorm = sqrt(norm2(mmp));
|
||||
RealD psinorm = sqrt(norm2(x));
|
||||
RealD srcnorm = sqrt(norm2(src));
|
||||
RealD tmpnorm = sqrt(norm2(tmp));
|
||||
RealD true_residual = tmpnorm/srcnorm;
|
||||
std::cout<<GridLogMessage<<"HDCG: true residual is "<<true_residual
|
||||
<<" solution "<<psinorm<<" source "<<srcnorm<<std::endl;
|
||||
|
||||
return;
|
||||
}
|
||||
|
||||
}
|
||||
std::cout << "HDCG: Pcg not converged"<<std::endl;
|
||||
return ;
|
||||
}
|
||||
|
||||
// The Pcg routine is common to all, but the various matrices differ from derived
|
||||
// implementation to derived implmentation
|
||||
void operator() (const Field &src, Field &psi){
|
||||
void operator() (const Field &src, Field &psi){
|
||||
|
||||
psi.Checkerboard() = src.Checkerboard();
|
||||
@ -211,7 +108,7 @@ class TwoLevelFlexiblePcg : public LinearFunction<Field>
|
||||
std::vector<Field> mmp(mmax,grid);
|
||||
std::vector<RealD> pAp(mmax);
|
||||
|
||||
Field x (grid);
|
||||
Field x (grid); x = psi;
|
||||
Field z (grid);
|
||||
Field tmp(grid);
|
||||
Field r (grid);
|
||||
@ -220,23 +117,25 @@ class TwoLevelFlexiblePcg : public LinearFunction<Field>
|
||||
//////////////////////////
|
||||
// x0 = Vstart -- possibly modify guess
|
||||
//////////////////////////
|
||||
x=Zero();
|
||||
x=src;
|
||||
Vstart(x,src);
|
||||
|
||||
// r0 = b -A x0
|
||||
_FineLinop.HermOp(x,mmp[0]); // Fine operator
|
||||
HermOp(x,mmp); // Shouldn't this be something else?
|
||||
axpy (r, -1.0,mmp[0], src); // Recomputes r=src-Ax0
|
||||
|
||||
//////////////////////////////////
|
||||
// Compute z = M1 r
|
||||
// Compute z = M1 x
|
||||
//////////////////////////////////
|
||||
PcgM1(r,z);
|
||||
M1(r,z,tmp,mp,SmootherMirs);
|
||||
rtzp =real(innerProduct(r,z));
|
||||
|
||||
///////////////////////////////////////
|
||||
// Solve for Mss mu = P A z and set p = z-mu
|
||||
// Def2: p = 1 - Q Az = Pright z
|
||||
// Other algos M2 is trivial
|
||||
///////////////////////////////////////
|
||||
PcgM2(z,p[0]);
|
||||
M2(z,p[0]);
|
||||
|
||||
for (int k=0;k<=MaxIterations;k++){
|
||||
|
||||
@ -244,38 +143,26 @@ class TwoLevelFlexiblePcg : public LinearFunction<Field>
|
||||
int peri_kp = (k+1) % mmax;
|
||||
|
||||
rtz=rtzp;
|
||||
d= PcgM3(p[peri_k],mmp[peri_k]);
|
||||
d= M3(p[peri_k],mp,mmp[peri_k],tmp);
|
||||
a = rtz/d;
|
||||
|
||||
// Memorise this
|
||||
pAp[peri_k] = d;
|
||||
std::cout << GridLogMessage << " pCG d "<< d<<std::endl;
|
||||
|
||||
axpy(x,a,p[peri_k],x);
|
||||
// std::cout << GridLogMessage << " pCG x "<< norm2(x)<<std::endl;
|
||||
RealD rn = axpy_norm(r,-a,mmp[peri_k],r);
|
||||
|
||||
std::cout << GridLogMessage << " pCG rn "<< rn<<std::endl;
|
||||
|
||||
// Compute z = M x
|
||||
PcgM1(r,z);
|
||||
// std::cout << GridLogMessage << " pCG z "<< norm2(z)<<std::endl;
|
||||
M1(r,z,tmp,mp);
|
||||
|
||||
rtzp =real(innerProduct(r,z));
|
||||
std::cout << GridLogMessage << " pCG rtzp "<<rtzp<<std::endl;
|
||||
// std::cout << GridLogMessage << " pCG r "<<norm2(r)<<std::endl;
|
||||
|
||||
PcgM2(z,mu); // ADEF-2 this is identity. Axpy possible to eliminate
|
||||
M2(z,mu); // ADEF-2 this is identity. Axpy possible to eliminate
|
||||
|
||||
// std::cout << GridLogMessage << " pCG mu "<<norm2(mu)<<std::endl;
|
||||
|
||||
p[peri_kp]=mu;
|
||||
p[peri_kp]=p[peri_k];
|
||||
|
||||
// std::cout << GridLogMessage << " pCG p[peri_kp] "<<norm2(p[peri_kp])<<std::endl;
|
||||
|
||||
// Standard search direction p -> z + b p
|
||||
// Standard search direction p -> z + b p ; b =
|
||||
b = (rtzp)/rtz;
|
||||
std::cout << GridLogMessage << " pCG b "<< b<<std::endl;
|
||||
|
||||
int northog;
|
||||
// northog = (peri_kp==0)?1:peri_kp; // This is the fCG(mmax) algorithm
|
||||
@ -287,7 +174,6 @@ class TwoLevelFlexiblePcg : public LinearFunction<Field>
|
||||
RealD beta = -pbApk/pAp[peri_back];
|
||||
axpy(p[peri_kp],beta,p[peri_back],p[peri_kp]);
|
||||
}
|
||||
// std::cout << GridLogMessage << " pCG p[peri_kp] orthog "<< norm2(p[peri_kp])<<std::endl;
|
||||
|
||||
RealD rrn=sqrt(rn/ssq);
|
||||
std::cout<<GridLogMessage<<"TwoLevelfPcg: k= "<<k<<" residual = "<<rrn<<std::endl;
|
||||
@ -295,7 +181,7 @@ class TwoLevelFlexiblePcg : public LinearFunction<Field>
|
||||
// Stopping condition
|
||||
if ( rn <= rsq ) {
|
||||
|
||||
_FineLinop.HermOp(x,mmp[0]); // Shouldn't this be something else?
|
||||
HermOp(x,mmp); // Shouldn't this be something else?
|
||||
axpy(tmp,-1.0,src,mmp[0]);
|
||||
|
||||
RealD psinorm = sqrt(norm2(x));
|
||||
@ -304,7 +190,7 @@ class TwoLevelFlexiblePcg : public LinearFunction<Field>
|
||||
RealD true_residual = tmpnorm/srcnorm;
|
||||
std::cout<<GridLogMessage<<"TwoLevelfPcg: true residual is "<<true_residual<<std::endl;
|
||||
std::cout<<GridLogMessage<<"TwoLevelfPcg: target residual was"<<Tolerance<<std::endl;
|
||||
return;
|
||||
return k;
|
||||
}
|
||||
}
|
||||
// Non-convergence
|
||||
@ -313,42 +199,52 @@ class TwoLevelFlexiblePcg : public LinearFunction<Field>
|
||||
|
||||
public:
|
||||
|
||||
virtual void PcgM1(Field & in, Field & out)
|
||||
{
|
||||
// [PTM+Q] in = [1 - Q A] M in + Q in = Min + Q [ in -A Min]
|
||||
virtual void M(Field & in,Field & out,Field & tmp) {
|
||||
|
||||
}
|
||||
|
||||
virtual void M1(Field & in, Field & out) {// the smoother
|
||||
|
||||
// [PTM+Q] in = [1 - Q A] M in + Q in = Min + Q [ in -A Min]
|
||||
Field tmp(grid);
|
||||
Field Min(grid);
|
||||
CoarseField PleftProj(coarsegrid);
|
||||
CoarseField PleftMss_proj(coarsegrid);
|
||||
|
||||
_Smoother(in,Min);
|
||||
PcgM(in,Min); // Smoother call
|
||||
|
||||
_FineLinop.HermOp(Min,out);
|
||||
HermOp(Min,out);
|
||||
axpy(tmp,-1.0,out,in); // tmp = in - A Min
|
||||
|
||||
_Aggregates.ProjectToSubspace(PleftProj,tmp);
|
||||
_CoarseSolver(PleftProj,PleftMss_proj); // Ass^{-1} [in - A Min]_s
|
||||
_Aggregates.PromoteFromSubspace(PleftMss_proj,tmp);// tmp = Q[in - A Min]
|
||||
|
||||
ProjectToSubspace(tmp,PleftProj);
|
||||
ApplyInverse(PleftProj,PleftMss_proj); // Ass^{-1} [in - A Min]_s
|
||||
PromoteFromSubspace(PleftMss_proj,tmp);// tmp = Q[in - A Min]
|
||||
axpy(out,1.0,Min,tmp); // Min+tmp
|
||||
|
||||
}
|
||||
|
||||
virtual void PcgM2(const Field & in, Field & out) {
|
||||
virtual void M2(const Field & in, Field & out) {
|
||||
out=in;
|
||||
// Must override for Def2 only
|
||||
// case PcgDef2:
|
||||
// Pright(in,out);
|
||||
// break;
|
||||
}
|
||||
|
||||
virtual RealD PcgM3(const Field & p, Field & mmp){
|
||||
RealD dd;
|
||||
_FineLinop.HermOp(p,mmp);
|
||||
ComplexD dot = innerProduct(p,mmp);
|
||||
dd=real(dot);
|
||||
virtual RealD M3(const Field & p, Field & mmp){
|
||||
double d,dd;
|
||||
HermOpAndNorm(p,mmp,d,dd);
|
||||
return dd;
|
||||
// Must override for Def1 only
|
||||
// case PcgDef1:
|
||||
// d=linop_d->Mprec(p,mmp,tmp,0,1);// Dag no
|
||||
// linop_d->Mprec(mmp,mp,tmp,1);// Dag yes
|
||||
// Pleft(mp,mmp);
|
||||
// d=real(linop_d->inner(p,mmp));
|
||||
}
|
||||
|
||||
virtual void Vstart(Field & x,const Field & src)
|
||||
{
|
||||
virtual void VstartDef2(Field & xconst Field & src){
|
||||
//case PcgDef2:
|
||||
//case PcgAdef2:
|
||||
//case PcgAdef2f:
|
||||
//case PcgV11f:
|
||||
///////////////////////////////////
|
||||
// Choose x_0 such that
|
||||
// x_0 = guess + (A_ss^inv) r_s = guess + Ass_inv [src -Aguess]
|
||||
@ -362,22 +258,140 @@ class TwoLevelFlexiblePcg : public LinearFunction<Field>
|
||||
///////////////////////////////////
|
||||
Field r(grid);
|
||||
Field mmp(grid);
|
||||
CoarseField PleftProj(coarsegrid);
|
||||
CoarseField PleftMss_proj(coarsegrid);
|
||||
|
||||
_Aggregates.ProjectToSubspace(PleftProj,src);
|
||||
_CoarseSolverPrecise(PleftProj,PleftMss_proj); // Ass^{-1} r_s
|
||||
_Aggregates.PromoteFromSubspace(PleftMss_proj,x);
|
||||
|
||||
HermOp(x,mmp);
|
||||
axpy (r, -1.0, mmp, src); // r_{-1} = src - A x
|
||||
ProjectToSubspace(r,PleftProj);
|
||||
ApplyInverseCG(PleftProj,PleftMss_proj); // Ass^{-1} r_s
|
||||
PromoteFromSubspace(PleftMss_proj,mmp);
|
||||
x=x+mmp;
|
||||
|
||||
}
|
||||
|
||||
virtual void Vstart(Field & x,const Field & src){
|
||||
return;
|
||||
}
|
||||
|
||||
/////////////////////////////////////////////////////////////////////
|
||||
// Only Def1 has non-trivial Vout.
|
||||
// Only Def1 has non-trivial Vout. Override in Def1
|
||||
/////////////////////////////////////////////////////////////////////
|
||||
virtual void Vout (Field & in, Field & out,Field & src){
|
||||
out = in;
|
||||
//case PcgDef1:
|
||||
// //Qb + PT x
|
||||
// ProjectToSubspace(src,PleftProj);
|
||||
// ApplyInverse(PleftProj,PleftMss_proj); // Ass^{-1} r_s
|
||||
// PromoteFromSubspace(PleftMss_proj,tmp);
|
||||
//
|
||||
// Pright(in,out);
|
||||
//
|
||||
// linop_d->axpy(out,tmp,out,1.0);
|
||||
// break;
|
||||
}
|
||||
};
|
||||
|
||||
NAMESPACE_END(Grid);
|
||||
////////////////////////////////////////////////////////////////////////////////////////////////
|
||||
// Pright and Pleft are common to all implementations
|
||||
////////////////////////////////////////////////////////////////////////////////////////////////
|
||||
virtual void Pright(Field & in,Field & out){
|
||||
// P_R = [ 1 0 ]
|
||||
// [ -Mss^-1 Msb 0 ]
|
||||
Field in_sbar(grid);
|
||||
|
||||
ProjectToSubspace(in,PleftProj);
|
||||
PromoteFromSubspace(PleftProj,out);
|
||||
axpy(in_sbar,-1.0,out,in); // in_sbar = in - in_s
|
||||
|
||||
HermOp(in_sbar,out);
|
||||
ProjectToSubspace(out,PleftProj); // Mssbar in_sbar (project)
|
||||
|
||||
ApplyInverse (PleftProj,PleftMss_proj); // Mss^{-1} Mssbar
|
||||
PromoteFromSubspace(PleftMss_proj,out); //
|
||||
|
||||
axpy(out,-1.0,out,in_sbar); // in_sbar - Mss^{-1} Mssbar in_sbar
|
||||
}
|
||||
virtual void Pleft (Field & in,Field & out){
|
||||
// P_L = [ 1 -Mbs Mss^-1]
|
||||
// [ 0 0 ]
|
||||
Field in_sbar(grid);
|
||||
Field tmp2(grid);
|
||||
Field Mtmp(grid);
|
||||
|
||||
ProjectToSubspace(in,PleftProj);
|
||||
PromoteFromSubspace(PleftProj,out);
|
||||
axpy(in_sbar,-1.0,out,in); // in_sbar = in - in_s
|
||||
|
||||
ApplyInverse(PleftProj,PleftMss_proj); // Mss^{-1} in_s
|
||||
PromoteFromSubspace(PleftMss_proj,out);
|
||||
|
||||
HermOp(out,Mtmp);
|
||||
|
||||
ProjectToSubspace(Mtmp,PleftProj); // Msbar s Mss^{-1}
|
||||
PromoteFromSubspace(PleftProj,tmp2);
|
||||
|
||||
axpy(out,-1.0,tmp2,Mtmp);
|
||||
axpy(out,-1.0,out,in_sbar); // in_sbar - Msbars Mss^{-1} in_s
|
||||
}
|
||||
}
|
||||
|
||||
template<class Field>
|
||||
class TwoLevelFlexiblePcgADef2 : public TwoLevelFlexiblePcg<Field> {
|
||||
public:
|
||||
virtual void M(Field & in,Field & out,Field & tmp){
|
||||
|
||||
}
|
||||
virtual void M1(Field & in, Field & out,Field & tmp,Field & mp){
|
||||
|
||||
}
|
||||
virtual void M2(Field & in, Field & out){
|
||||
|
||||
}
|
||||
virtual RealD M3(Field & p, Field & mp,Field & mmp, Field & tmp){
|
||||
|
||||
}
|
||||
virtual void Vstart(Field & in, Field & src, Field & r, Field & mp, Field & mmp, Field & tmp){
|
||||
|
||||
}
|
||||
}
|
||||
/*
|
||||
template<class Field>
|
||||
class TwoLevelFlexiblePcgAD : public TwoLevelFlexiblePcg<Field> {
|
||||
public:
|
||||
virtual void M(Field & in,Field & out,Field & tmp);
|
||||
virtual void M1(Field & in, Field & out,Field & tmp,Field & mp);
|
||||
virtual void M2(Field & in, Field & out);
|
||||
virtual RealD M3(Field & p, Field & mp,Field & mmp, Field & tmp);
|
||||
virtual void Vstart(Field & in, Field & src, Field & r, Field & mp, Field & mmp, Field & tmp);
|
||||
}
|
||||
|
||||
template<class Field>
|
||||
class TwoLevelFlexiblePcgDef1 : public TwoLevelFlexiblePcg<Field> {
|
||||
public:
|
||||
virtual void M(Field & in,Field & out,Field & tmp);
|
||||
virtual void M1(Field & in, Field & out,Field & tmp,Field & mp);
|
||||
virtual void M2(Field & in, Field & out);
|
||||
virtual RealD M3(Field & p, Field & mp,Field & mmp, Field & tmp);
|
||||
virtual void Vstart(Field & in, Field & src, Field & r, Field & mp, Field & mmp, Field & tmp);
|
||||
virtual void Vout (Field & in, Field & out,Field & src,Field & tmp);
|
||||
}
|
||||
|
||||
template<class Field>
|
||||
class TwoLevelFlexiblePcgDef2 : public TwoLevelFlexiblePcg<Field> {
|
||||
public:
|
||||
virtual void M(Field & in,Field & out,Field & tmp);
|
||||
virtual void M1(Field & in, Field & out,Field & tmp,Field & mp);
|
||||
virtual void M2(Field & in, Field & out);
|
||||
virtual RealD M3(Field & p, Field & mp,Field & mmp, Field & tmp);
|
||||
virtual void Vstart(Field & in, Field & src, Field & r, Field & mp, Field & mmp, Field & tmp);
|
||||
}
|
||||
|
||||
template<class Field>
|
||||
class TwoLevelFlexiblePcgV11: public TwoLevelFlexiblePcg<Field> {
|
||||
public:
|
||||
virtual void M(Field & in,Field & out,Field & tmp);
|
||||
virtual void M1(Field & in, Field & out,Field & tmp,Field & mp);
|
||||
virtual void M2(Field & in, Field & out);
|
||||
virtual RealD M3(Field & p, Field & mp,Field & mmp, Field & tmp);
|
||||
virtual void Vstart(Field & in, Field & src, Field & r, Field & mp, Field & mmp, Field & tmp);
|
||||
}
|
||||
*/
|
||||
#endif
|
||||
|
||||
@ -183,13 +183,13 @@ public:
|
||||
<< "\tTrue residual " << true_residual
|
||||
<< "\tTarget " << Tolerance << std::endl;
|
||||
|
||||
std::cout << GridLogMessage << "Time breakdown "<<std::endl;
|
||||
std::cout << GridLogMessage << "\tElapsed " << SolverTimer.Elapsed() <<std::endl;
|
||||
std::cout << GridLogPerformance << "Time breakdown "<<std::endl;
|
||||
std::cout << GridLogPerformance << "\tMatrix " << MatrixTimer.Elapsed() <<std::endl;
|
||||
std::cout << GridLogPerformance << "\tLinalg " << LinalgTimer.Elapsed() <<std::endl;
|
||||
std::cout << GridLogPerformance << "\tInner " << InnerTimer.Elapsed() <<std::endl;
|
||||
std::cout << GridLogPerformance << "\tAxpyNorm " << AxpyNormTimer.Elapsed() <<std::endl;
|
||||
std::cout << GridLogPerformance << "\tLinearComb " << LinearCombTimer.Elapsed() <<std::endl;
|
||||
std::cout << GridLogMessage << "\tMatrix " << MatrixTimer.Elapsed() <<std::endl;
|
||||
std::cout << GridLogMessage << "\tLinalg " << LinalgTimer.Elapsed() <<std::endl;
|
||||
std::cout << GridLogMessage << "\tInner " << InnerTimer.Elapsed() <<std::endl;
|
||||
std::cout << GridLogMessage << "\tAxpyNorm " << AxpyNormTimer.Elapsed() <<std::endl;
|
||||
std::cout << GridLogMessage << "\tLinearComb " << LinearCombTimer.Elapsed() <<std::endl;
|
||||
|
||||
std::cout << GridLogDebug << "\tMobius flop rate " << DwfFlops/ usecs<< " Gflops " <<std::endl;
|
||||
|
||||
|
||||
@ -33,7 +33,7 @@ NAMESPACE_BEGIN(Grid);
|
||||
///////////////////////////////////////////////////////////////////////////////////////////////////////
|
||||
// Take a matrix and form an NE solver calling a Herm solver
|
||||
///////////////////////////////////////////////////////////////////////////////////////////////////////
|
||||
template<class Field> class NormalEquations : public LinearFunction<Field>{
|
||||
template<class Field> class NormalEquations {
|
||||
private:
|
||||
SparseMatrixBase<Field> & _Matrix;
|
||||
OperatorFunction<Field> & _HermitianSolver;
|
||||
@ -60,7 +60,7 @@ public:
|
||||
}
|
||||
};
|
||||
|
||||
template<class Field> class HPDSolver : public LinearFunction<Field> {
|
||||
template<class Field> class HPDSolver {
|
||||
private:
|
||||
LinearOperatorBase<Field> & _Matrix;
|
||||
OperatorFunction<Field> & _HermitianSolver;
|
||||
@ -84,7 +84,7 @@ public:
|
||||
};
|
||||
|
||||
|
||||
template<class Field> class MdagMSolver : public LinearFunction<Field> {
|
||||
template<class Field> class MdagMSolver {
|
||||
private:
|
||||
SparseMatrixBase<Field> & _Matrix;
|
||||
OperatorFunction<Field> & _HermitianSolver;
|
||||
|
||||
@ -20,7 +20,7 @@ template<class Field> class PowerMethod
|
||||
RealD evalMaxApprox = 0.0;
|
||||
auto src_n = src;
|
||||
auto tmp = src;
|
||||
const int _MAX_ITER_EST_ = 100;
|
||||
const int _MAX_ITER_EST_ = 50;
|
||||
|
||||
for (int i=0;i<_MAX_ITER_EST_;i++) {
|
||||
|
||||
|
||||
@ -361,14 +361,9 @@ public:
|
||||
_bernoulli.resize(_vol,std::discrete_distribution<int32_t>{1,1});
|
||||
_uid.resize(_vol,std::uniform_int_distribution<uint32_t>() );
|
||||
}
|
||||
template <class vobj,class distribution> inline void fill(Lattice<vobj> &l,std::vector<distribution> &dist)
|
||||
{
|
||||
if ( l.Grid()->_isCheckerBoarded ) {
|
||||
Lattice<vobj> tmp(_grid);
|
||||
fill(tmp,dist);
|
||||
pickCheckerboard(l.Checkerboard(),l,tmp);
|
||||
return;
|
||||
}
|
||||
|
||||
template <class vobj,class distribution> inline void fill(Lattice<vobj> &l,std::vector<distribution> &dist){
|
||||
|
||||
typedef typename vobj::scalar_object scalar_object;
|
||||
typedef typename vobj::scalar_type scalar_type;
|
||||
typedef typename vobj::vector_type vector_type;
|
||||
|
||||
@ -45,9 +45,8 @@ public:
|
||||
dims=_grid->Nd();
|
||||
AllocateGrids();
|
||||
Coordinate local =unpadded_grid->LocalDimensions();
|
||||
Coordinate procs =unpadded_grid->ProcessorGrid();
|
||||
for(int d=0;d<dims;d++){
|
||||
if ( procs[d] > 1 ) assert(local[d]>=depth);
|
||||
assert(local[d]>=depth);
|
||||
}
|
||||
}
|
||||
void DeleteGrids(void)
|
||||
@ -112,7 +111,7 @@ public:
|
||||
if(dim==0) conformable(old_grid,unpadded_grid);
|
||||
else conformable(old_grid,grids[dim-1]);
|
||||
|
||||
// std::cout << " dim "<<dim<<" local "<<local << " padding to "<<plocal<<std::endl;
|
||||
std::cout << " dim "<<dim<<" local "<<local << " padding to "<<plocal<<std::endl;
|
||||
// Middle bit
|
||||
for(int x=0;x<local[dim];x++){
|
||||
InsertSliceLocal(in,padded,x,depth+x,dim);
|
||||
|
||||
@ -124,6 +124,11 @@ public:
|
||||
RealD _b;
|
||||
RealD _c;
|
||||
|
||||
// possible boost
|
||||
std::vector<ComplexD> qmu;
|
||||
void set_qmu(std::vector<ComplexD> _qmu) { qmu=_qmu; assert(qmu.size()==Nd);};
|
||||
void addQmu(const FermionField &in, FermionField &out, int dag);
|
||||
|
||||
// Cayley form Moebius (tanh and zolotarev)
|
||||
Vector<Coeff_t> omega;
|
||||
Vector<Coeff_t> bs; // S dependent coeffs
|
||||
|
||||
@ -60,6 +60,50 @@ public:
|
||||
// virtual void Instantiatable(void)=0;
|
||||
virtual void Instantiatable(void) =0;
|
||||
|
||||
void FreePropagator(const FermionField &in,FermionField &out,RealD mass,std::vector<Complex> boundary, std::vector<double> twist)
|
||||
{
|
||||
std::cout << "Free Propagator for PartialFraction"<<std::endl;
|
||||
FermionField in_k(in.Grid());
|
||||
FermionField prop_k(in.Grid());
|
||||
|
||||
FFT theFFT((GridCartesian *) in.Grid());
|
||||
|
||||
//phase for boundary condition
|
||||
ComplexField coor(in.Grid());
|
||||
ComplexField ph(in.Grid()); ph = Zero();
|
||||
FermionField in_buf(in.Grid()); in_buf = Zero();
|
||||
typedef typename Simd::scalar_type Scalar;
|
||||
Scalar ci(0.0,1.0);
|
||||
assert(twist.size() == Nd);//check that twist is Nd
|
||||
assert(boundary.size() == Nd);//check that boundary conditions is Nd
|
||||
int shift = 0;
|
||||
for(unsigned int nu = 0; nu < Nd; nu++)
|
||||
{
|
||||
// Shift coordinate lattice index by 1 to account for 5th dimension.
|
||||
LatticeCoordinate(coor, nu + shift);
|
||||
double boundary_phase = ::acos(real(boundary[nu]));
|
||||
ph = ph + boundary_phase*coor*((1./(in.Grid()->_fdimensions[nu+shift])));
|
||||
//momenta for propagator shifted by twist+boundary
|
||||
twist[nu] = twist[nu] + boundary_phase/((2.0*M_PI));
|
||||
}
|
||||
in_buf = exp(ci*ph*(-1.0))*in;
|
||||
|
||||
theFFT.FFT_all_dim(in_k,in,FFT::forward);
|
||||
this->MomentumSpacePropagatorHw(prop_k,in_k,mass,twist);
|
||||
theFFT.FFT_all_dim(out,prop_k,FFT::backward);
|
||||
|
||||
//phase for boundary condition
|
||||
out = out * exp(ci*ph);
|
||||
};
|
||||
|
||||
virtual void FreePropagator(const FermionField &in,FermionField &out,RealD mass) {
|
||||
std::vector<double> twist(Nd,0.0); //default: periodic boundarys in all directions
|
||||
std::vector<Complex> boundary;
|
||||
for(int i=0;i<Nd;i++) boundary.push_back(1);//default: periodic boundary conditions
|
||||
FreePropagator(in,out,mass,boundary,twist);
|
||||
};
|
||||
|
||||
|
||||
// Efficient support for multigrid coarsening
|
||||
virtual void Mdir (const FermionField &in, FermionField &out,int dir,int disp);
|
||||
virtual void MdirAll(const FermionField &in, std::vector<FermionField> &out);
|
||||
|
||||
@ -39,7 +39,7 @@ class PartialFractionFermion5D : public WilsonFermion5D<Impl>
|
||||
public:
|
||||
INHERIT_IMPL_TYPES(Impl);
|
||||
|
||||
const int part_frac_chroma_convention=1;
|
||||
const int part_frac_chroma_convention=0;
|
||||
|
||||
void Meooe_internal(const FermionField &in, FermionField &out,int dag);
|
||||
void Mooee_internal(const FermionField &in, FermionField &out,int dag);
|
||||
@ -83,12 +83,63 @@ public:
|
||||
GridRedBlackCartesian &FourDimRedBlackGrid,
|
||||
RealD _mass,RealD M5,const ImplParams &p= ImplParams());
|
||||
|
||||
PartialFractionFermion5D(GaugeField &_Umu,
|
||||
GridCartesian &FiveDimGrid,
|
||||
GridRedBlackCartesian &FiveDimRedBlackGrid,
|
||||
GridCartesian &FourDimGrid,
|
||||
GridRedBlackCartesian &FourDimRedBlackGrid,
|
||||
RealD _mass,RealD M5,std::vector<RealD> &_qmu,const ImplParams &p= ImplParams());
|
||||
|
||||
void FreePropagator(const FermionField &in,FermionField &out,RealD mass,std::vector<Complex> boundary, std::vector<double> twist)
|
||||
{
|
||||
std::cout << "Free Propagator for PartialFraction"<<std::endl;
|
||||
FermionField in_k(in.Grid());
|
||||
FermionField prop_k(in.Grid());
|
||||
|
||||
FFT theFFT((GridCartesian *) in.Grid());
|
||||
|
||||
//phase for boundary condition
|
||||
ComplexField coor(in.Grid());
|
||||
ComplexField ph(in.Grid()); ph = Zero();
|
||||
FermionField in_buf(in.Grid()); in_buf = Zero();
|
||||
typedef typename Simd::scalar_type Scalar;
|
||||
Scalar ci(0.0,1.0);
|
||||
assert(twist.size() == Nd);//check that twist is Nd
|
||||
assert(boundary.size() == Nd);//check that boundary conditions is Nd
|
||||
int shift = 0;
|
||||
for(unsigned int nu = 0; nu < Nd; nu++)
|
||||
{
|
||||
// Shift coordinate lattice index by 1 to account for 5th dimension.
|
||||
LatticeCoordinate(coor, nu + shift);
|
||||
double boundary_phase = ::acos(real(boundary[nu]));
|
||||
ph = ph + boundary_phase*coor*((1./(in.Grid()->_fdimensions[nu+shift])));
|
||||
//momenta for propagator shifted by twist+boundary
|
||||
twist[nu] = twist[nu] + boundary_phase/((2.0*M_PI));
|
||||
}
|
||||
in_buf = exp(ci*ph*(-1.0))*in;
|
||||
|
||||
theFFT.FFT_all_dim(in_k,in,FFT::forward);
|
||||
this->MomentumSpacePropagatorHw(prop_k,in_k,mass,twist);
|
||||
theFFT.FFT_all_dim(out,prop_k,FFT::backward);
|
||||
|
||||
//phase for boundary condition
|
||||
out = out * exp(ci*ph);
|
||||
};
|
||||
|
||||
virtual void FreePropagator(const FermionField &in,FermionField &out,RealD mass) {
|
||||
std::vector<double> twist(Nd,0.0); //default: periodic boundarys in all directions
|
||||
std::vector<Complex> boundary;
|
||||
for(int i=0;i<Nd;i++) boundary.push_back(1);//default: periodic boundary conditions
|
||||
FreePropagator(in,out,mass,boundary,twist);
|
||||
};
|
||||
|
||||
protected:
|
||||
|
||||
virtual void SetCoefficientsTanh(Approx::zolotarev_data *zdata,RealD scale);
|
||||
virtual void SetCoefficientsZolotarev(RealD zolo_hi,Approx::zolotarev_data *zdata);
|
||||
|
||||
// Part frac
|
||||
std::vector<RealD> qmu;
|
||||
RealD mass;
|
||||
RealD dw_diag;
|
||||
RealD R;
|
||||
|
||||
@ -48,7 +48,8 @@ CayleyFermion5D<Impl>::CayleyFermion5D(GaugeField &_Umu,
|
||||
FourDimGrid,
|
||||
FourDimRedBlackGrid,_M5,p),
|
||||
mass_plus(_mass), mass_minus(_mass)
|
||||
{
|
||||
{
|
||||
// qmu defaults to zero size;
|
||||
}
|
||||
|
||||
///////////////////////////////////////////////////////////////
|
||||
@ -270,6 +271,34 @@ void CayleyFermion5D<Impl>::MeooeDag5D (const FermionField &psi, FermionField
|
||||
M5Ddag(psi,psi,Din,lower,diag,upper);
|
||||
}
|
||||
|
||||
template<class Impl>
|
||||
void CayleyFermion5D<Impl>::addQmu(const FermionField &psi,FermionField &chi, int dag)
|
||||
{
|
||||
if ( qmu.size() ) {
|
||||
|
||||
Gamma::Algebra Gmu [] = {
|
||||
Gamma::Algebra::GammaX,
|
||||
Gamma::Algebra::GammaY,
|
||||
Gamma::Algebra::GammaZ,
|
||||
Gamma::Algebra::GammaT
|
||||
};
|
||||
std::vector<ComplexD> coeff(Nd);
|
||||
ComplexD ci(0,1);
|
||||
|
||||
assert(qmu.size()==Nd);
|
||||
|
||||
for(int mu=0;mu<Nd;mu++){
|
||||
coeff[mu] = ci*qmu[mu];
|
||||
if ( dag ) coeff[mu] = conjugate(coeff[mu]);
|
||||
}
|
||||
|
||||
chi = chi + Gamma(Gmu[0])*psi*coeff[0];
|
||||
for(int mu=1;mu<Nd;mu++){
|
||||
chi = chi + Gamma(Gmu[mu])*psi*coeff[mu];
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
template<class Impl>
|
||||
void CayleyFermion5D<Impl>::M (const FermionField &psi, FermionField &chi)
|
||||
{
|
||||
@ -277,8 +306,12 @@ void CayleyFermion5D<Impl>::M (const FermionField &psi, FermionField &chi)
|
||||
|
||||
// Assemble Din
|
||||
Meooe5D(psi,Din);
|
||||
|
||||
|
||||
this->DW(Din,chi,DaggerNo);
|
||||
|
||||
// add i q_mu gamma_mu here
|
||||
addQmu(Din,chi,DaggerNo);
|
||||
|
||||
// ((b D_W + D_w hop terms +1) on s-diag
|
||||
axpby(chi,1.0,1.0,chi,psi);
|
||||
|
||||
@ -295,6 +328,9 @@ void CayleyFermion5D<Impl>::Mdag (const FermionField &psi, FermionField &chi)
|
||||
FermionField Din(psi.Grid());
|
||||
// Apply Dw
|
||||
this->DW(psi,Din,DaggerYes);
|
||||
|
||||
// add -i conj(q_mu) gamma_mu here ... if qmu is real, gammm_5 hermitian, otherwise not.
|
||||
addQmu(psi,Din,DaggerYes);
|
||||
|
||||
MeooeDag5D(Din,chi);
|
||||
|
||||
|
||||
@ -42,13 +42,13 @@ template<class Impl>
|
||||
void ContinuedFractionFermion5D<Impl>::SetCoefficientsZolotarev(RealD zolo_hi,Approx::zolotarev_data *zdata)
|
||||
{
|
||||
// How to check Ls matches??
|
||||
// std::cout<<GridLogMessage << Ls << " Ls"<<std::endl;
|
||||
// std::cout<<GridLogMessage << zdata->n << " - n"<<std::endl;
|
||||
// std::cout<<GridLogMessage << zdata->da << " -da "<<std::endl;
|
||||
// std::cout<<GridLogMessage << zdata->db << " -db"<<std::endl;
|
||||
// std::cout<<GridLogMessage << zdata->dn << " -dn"<<std::endl;
|
||||
// std::cout<<GridLogMessage << zdata->dd << " -dd"<<std::endl;
|
||||
std::cout<<GridLogMessage << zdata->n << " - n"<<std::endl;
|
||||
std::cout<<GridLogMessage << zdata->da << " -da "<<std::endl;
|
||||
std::cout<<GridLogMessage << zdata->db << " -db"<<std::endl;
|
||||
std::cout<<GridLogMessage << zdata->dn << " -dn"<<std::endl;
|
||||
std::cout<<GridLogMessage << zdata->dd << " -dd"<<std::endl;
|
||||
int Ls = this->Ls;
|
||||
std::cout<<GridLogMessage << Ls << " Ls"<<std::endl;
|
||||
assert(zdata->db==Ls);// Beta has Ls coeffs
|
||||
|
||||
R=(1+this->mass)/(1-this->mass);
|
||||
@ -320,7 +320,7 @@ ContinuedFractionFermion5D<Impl>::ContinuedFractionFermion5D(
|
||||
int Ls = this->Ls;
|
||||
conformable(solution5d.Grid(),this->FermionGrid());
|
||||
conformable(exported4d.Grid(),this->GaugeGrid());
|
||||
ExtractSlice(exported4d, solution5d, Ls-1, Ls-1);
|
||||
ExtractSlice(exported4d, solution5d, Ls-1, 0);
|
||||
}
|
||||
template<class Impl>
|
||||
void ContinuedFractionFermion5D<Impl>::ImportPhysicalFermionSource(const FermionField &input4d,FermionField &imported5d)
|
||||
@ -330,7 +330,7 @@ ContinuedFractionFermion5D<Impl>::ContinuedFractionFermion5D(
|
||||
conformable(input4d.Grid() ,this->GaugeGrid());
|
||||
FermionField tmp(this->FermionGrid());
|
||||
tmp=Zero();
|
||||
InsertSlice(input4d, tmp, Ls-1, Ls-1);
|
||||
InsertSlice(input4d, tmp, Ls-1, 0);
|
||||
tmp=Gamma(Gamma::Algebra::Gamma5)*tmp;
|
||||
this->Dminus(tmp,imported5d);
|
||||
}
|
||||
|
||||
@ -255,15 +255,76 @@ void PartialFractionFermion5D<Impl>::M_internal(const FermionField &psi, Fermi
|
||||
}
|
||||
|
||||
{
|
||||
// The 'conventional' Cayley overlap operator is
|
||||
//
|
||||
// Dov = (1+m)/2 + (1-m)/2 g5 sgn Hw
|
||||
//
|
||||
//
|
||||
// With massless limit 1/2(1+g5 sgnHw)
|
||||
//
|
||||
// Luscher shows quite neatly that 1+g5 sgn Hw has tree level propagator i qslash +O(a^2)
|
||||
//
|
||||
// However, the conventional normalisation has both a leading order factor of 2 in Zq
|
||||
// at tree level AND a mass dependent (1-m) that are convenient to absorb.
|
||||
//
|
||||
// In WilsonFermion5DImplementation.h, the tree level propagator for Hw is
|
||||
//
|
||||
// num = -i sin kmu gmu
|
||||
//
|
||||
// denom ( sqrt(sk^2 + (2shk^2 - 1)^2
|
||||
// b_k = sk2 - M5;
|
||||
//
|
||||
// w_k = sqrt(sk + b_k*b_k);
|
||||
//
|
||||
// denom= ( w_k + b_k + mass*mass) ;
|
||||
//
|
||||
// denom= one/denom;
|
||||
// out = num*denom;
|
||||
//
|
||||
// Chroma, and Grid define partial fraction via 4d operator
|
||||
//
|
||||
// Dpf = 2/(1-m) x Dov = (1+m)/(1-m) + g5 sgn Hw
|
||||
//
|
||||
// Now since:
|
||||
//
|
||||
// (1+m)/(1-m) = (1-m)/(1-m) + 2m/(1-m) = 1 + 2m/(1-m)
|
||||
//
|
||||
// This corresponds to a modified mass parameter
|
||||
//
|
||||
// It has an annoying
|
||||
//
|
||||
//
|
||||
double R=(1+this->mass)/(1-this->mass);
|
||||
//R g5 psi[Ls] + p[0] H
|
||||
ag5xpbg5y_ssp(chi,R*scale,psi,p[nblock]*scale/amax,D,Ls-1,Ls-1);
|
||||
|
||||
|
||||
for(int b=0;b<nblock;b++){
|
||||
int s = 2*b+1;
|
||||
double pp = p[nblock-1-b];
|
||||
axpby_ssp(chi,1.0,chi,-sqrt(amax*pp)*scale*sign,psi,Ls-1,s);
|
||||
}
|
||||
|
||||
if ( qmu.size() ) {
|
||||
|
||||
FermionField qslash_psi(psi.Grid());
|
||||
|
||||
Gamma::Algebra Gmu [] = {
|
||||
Gamma::Algebra::GammaX,
|
||||
Gamma::Algebra::GammaY,
|
||||
Gamma::Algebra::GammaZ,
|
||||
Gamma::Algebra::GammaT
|
||||
};
|
||||
ComplexD ci(0,1);
|
||||
assert(qmu.size()==Nd);
|
||||
qslash_psi = Gamma(Gmu[0])*psi;
|
||||
for(int mu=1;mu<Nd;mu++){
|
||||
qslash_psi = Gamma(Gmu[mu])*psi;
|
||||
}
|
||||
// RealD coeff = 1.0;
|
||||
qslash_psi = Gamma(Gamma::Algebra::Gamma5)*qslash_psi*ci ; // i g5 qslash -- 1-m factor???
|
||||
axpby_ssp(chi,1.0,chi,1.0, qslash_psi,Ls-1,Ls-1);
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
}
|
||||
@ -411,7 +472,7 @@ void PartialFractionFermion5D<Impl>::SetCoefficientsZolotarev(RealD zolo_hi,App
|
||||
int Ls = this->Ls;
|
||||
conformable(solution5d.Grid(),this->FermionGrid());
|
||||
conformable(exported4d.Grid(),this->GaugeGrid());
|
||||
ExtractSlice(exported4d, solution5d, Ls-1, Ls-1);
|
||||
ExtractSlice(exported4d, solution5d, Ls-1, 0);
|
||||
}
|
||||
template<class Impl>
|
||||
void PartialFractionFermion5D<Impl>::ImportPhysicalFermionSource(const FermionField &input4d,FermionField &imported5d)
|
||||
@ -421,7 +482,8 @@ void PartialFractionFermion5D<Impl>::SetCoefficientsZolotarev(RealD zolo_hi,App
|
||||
conformable(input4d.Grid() ,this->GaugeGrid());
|
||||
FermionField tmp(this->FermionGrid());
|
||||
tmp=Zero();
|
||||
InsertSlice(input4d, tmp, Ls-1, Ls-1);
|
||||
std::cout << " importing to slice " << Ls-1 <<std::endl;
|
||||
InsertSlice(input4d, tmp, Ls-1, 0);
|
||||
tmp=Gamma(Gamma::Algebra::Gamma5)*tmp;
|
||||
this->Dminus(tmp,imported5d);
|
||||
}
|
||||
@ -442,7 +504,7 @@ PartialFractionFermion5D<Impl>::PartialFractionFermion5D(GaugeField &_Umu,
|
||||
|
||||
{
|
||||
int Ls = this->Ls;
|
||||
|
||||
qmu.resize(0);
|
||||
assert((Ls&0x1)==1); // Odd Ls required
|
||||
int nrational=Ls-1;
|
||||
|
||||
@ -460,6 +522,22 @@ PartialFractionFermion5D<Impl>::PartialFractionFermion5D(GaugeField &_Umu,
|
||||
Approx::zolotarev_free(zdata);
|
||||
|
||||
}
|
||||
template<class Impl>
|
||||
PartialFractionFermion5D<Impl>::PartialFractionFermion5D(GaugeField &_Umu,
|
||||
GridCartesian &FiveDimGrid,
|
||||
GridRedBlackCartesian &FiveDimRedBlackGrid,
|
||||
GridCartesian &FourDimGrid,
|
||||
GridRedBlackCartesian &FourDimRedBlackGrid,
|
||||
RealD _mass,RealD M5,
|
||||
std::vector<RealD> &_qmu,
|
||||
const ImplParams &p)
|
||||
: PartialFractionFermion5D<Impl>(_Umu,
|
||||
FiveDimGrid,FiveDimRedBlackGrid,
|
||||
FourDimGrid,FourDimRedBlackGrid,
|
||||
_mass,M5,p)
|
||||
{
|
||||
qmu=_qmu;
|
||||
}
|
||||
|
||||
NAMESPACE_END(Grid);
|
||||
|
||||
|
||||
@ -46,7 +46,7 @@ class GeneralLocalStencilView {
|
||||
accelerator_inline GeneralStencilEntry * GetEntry(int point,int osite) {
|
||||
return & this->_entries_p[point+this->_npoints*osite];
|
||||
}
|
||||
void ViewClose(void){};
|
||||
|
||||
};
|
||||
////////////////////////////////////////
|
||||
// The Stencil Class itself
|
||||
@ -61,7 +61,7 @@ protected:
|
||||
public:
|
||||
GridBase *Grid(void) const { return _grid; }
|
||||
|
||||
View_type View(int mode) const {
|
||||
View_type View(void) const {
|
||||
View_type accessor(*( (View_type *) this));
|
||||
return accessor;
|
||||
}
|
||||
|
||||
@ -94,13 +94,6 @@ static constexpr int MaxDims = GRID_MAX_LATTICE_DIMENSION;
|
||||
|
||||
typedef AcceleratorVector<int,MaxDims> Coordinate;
|
||||
|
||||
template<class T,int _ndim>
|
||||
inline bool operator==(const AcceleratorVector<T,_ndim> &v,const AcceleratorVector<T,_ndim> &w)
|
||||
{
|
||||
if (v.size()!=w.size()) return false;
|
||||
for(int i=0;i<v.size();i++) if ( v[i]!=w[i] ) return false;
|
||||
return true;
|
||||
}
|
||||
template<class T,int _ndim>
|
||||
inline std::ostream & operator<<(std::ostream &os, const AcceleratorVector<T,_ndim> &v)
|
||||
{
|
||||
|
||||
@ -1,3 +1,4 @@
|
||||
BREW=/opt/local/
|
||||
CXXFLAGS=-fsanitize=address CXX=g++ ../../configure --enable-simd=NEONv8 --enable-comms=none --enable-unified=yes --prefix $HOME/QCD/GridInstall --with-lime=/Users/peterboyle/QCD/SciDAC/install/ --with-openssl=$BREW --disable-gparity --disable-fermion-reps
|
||||
CXX=mpicxx-openmpi-mp ../../configure --enable-simd=GEN --enable-comms=mpi --enable-unified=yes --prefix $HOME/QCD/GridInstall --with-lime=/Users/peterboyle/QCD/SciDAC/install/ --with-openssl=$BREW --disable-fermion-reps --disable-gparity --disable-debug
|
||||
|
||||
|
||||
|
||||
@ -1,238 +0,0 @@
|
||||
/*************************************************************************************
|
||||
|
||||
Grid physics library, www.github.com/paboyle/Grid
|
||||
|
||||
Source file: ./tests/Test_padded_cell.cc
|
||||
|
||||
Copyright (C) 2023
|
||||
|
||||
Author: Peter Boyle <paboyle@ph.ed.ac.uk>
|
||||
|
||||
This program is free software; you can redistribute it and/or modify
|
||||
it under the terms of the GNU General Public License as published by
|
||||
the Free Software Foundation; either version 2 of the License, or
|
||||
(at your option) any later version.
|
||||
|
||||
This program is distributed in the hope that it will be useful,
|
||||
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||
GNU General Public License for more details.
|
||||
|
||||
You should have received a copy of the GNU General Public License along
|
||||
with this program; if not, write to the Free Software Foundation, Inc.,
|
||||
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
|
||||
|
||||
See the full license in the file "LICENSE" in the top level distribution directory
|
||||
*************************************************************************************/
|
||||
/* END LEGAL */
|
||||
#include <Grid/Grid.h>
|
||||
#include <Grid/lattice/PaddedCell.h>
|
||||
#include <Grid/stencil/GeneralLocalStencil.h>
|
||||
#include <Grid/algorithms/GeneralCoarsenedMatrix.h>
|
||||
|
||||
#include <Grid/algorithms/iterative/PrecGeneralisedConjugateResidual.h>
|
||||
#include <Grid/algorithms/iterative/PrecGeneralisedConjugateResidualNonHermitian.h>
|
||||
#include <Grid/algorithms/iterative/BiCGSTAB.h>
|
||||
|
||||
using namespace std;
|
||||
using namespace Grid;
|
||||
|
||||
///////////////////////
|
||||
// Tells little dirac op to use MdagM as the .Op()
|
||||
///////////////////////
|
||||
template<class Field>
|
||||
class HermOpAdaptor : public LinearOperatorBase<Field>
|
||||
{
|
||||
LinearOperatorBase<Field> & wrapped;
|
||||
public:
|
||||
HermOpAdaptor(LinearOperatorBase<Field> &wrapme) : wrapped(wrapme) {};
|
||||
void OpDiag (const Field &in, Field &out) { assert(0); }
|
||||
void OpDir (const Field &in, Field &out,int dir,int disp) { assert(0); }
|
||||
void OpDirAll (const Field &in, std::vector<Field> &out){ assert(0); };
|
||||
void Op (const Field &in, Field &out){
|
||||
wrapped.HermOp(in,out);
|
||||
}
|
||||
void AdjOp (const Field &in, Field &out){
|
||||
wrapped.HermOp(in,out);
|
||||
}
|
||||
void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){ assert(0); }
|
||||
void HermOp(const Field &in, Field &out){
|
||||
wrapped.HermOp(in,out);
|
||||
}
|
||||
};
|
||||
|
||||
|
||||
int main (int argc, char ** argv)
|
||||
{
|
||||
Grid_init(&argc,&argv);
|
||||
|
||||
const int Ls=4;
|
||||
|
||||
GridCartesian * UGrid = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(),
|
||||
GridDefaultSimd(Nd,vComplex::Nsimd()),
|
||||
GridDefaultMpi());
|
||||
GridRedBlackCartesian * UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
|
||||
|
||||
GridCartesian * FGrid = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid);
|
||||
GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid);
|
||||
|
||||
// Construct a coarsened grid
|
||||
Coordinate clatt = GridDefaultLatt();
|
||||
for(int d=0;d<clatt.size();d++){
|
||||
clatt[d] = clatt[d]/2;
|
||||
}
|
||||
|
||||
GridCartesian *Coarse4d = SpaceTimeGrid::makeFourDimGrid(clatt,
|
||||
GridDefaultSimd(Nd,vComplex::Nsimd()),
|
||||
GridDefaultMpi());;
|
||||
GridCartesian *Coarse5d = SpaceTimeGrid::makeFiveDimGrid(1,Coarse4d);
|
||||
|
||||
std::vector<int> seeds4({1,2,3,4});
|
||||
std::vector<int> seeds5({5,6,7,8});
|
||||
std::vector<int> cseeds({5,6,7,8});
|
||||
GridParallelRNG RNG5(FGrid); RNG5.SeedFixedIntegers(seeds5);
|
||||
GridParallelRNG RNG4(UGrid); RNG4.SeedFixedIntegers(seeds4);
|
||||
GridParallelRNG CRNG(Coarse5d);CRNG.SeedFixedIntegers(cseeds);
|
||||
|
||||
LatticeFermion src(FGrid); random(RNG5,src);
|
||||
LatticeFermion result(FGrid); result=Zero();
|
||||
LatticeFermion ref(FGrid); ref=Zero();
|
||||
LatticeFermion tmp(FGrid);
|
||||
LatticeFermion err(FGrid);
|
||||
LatticeGaugeField Umu(UGrid);
|
||||
SU<Nc>::HotConfiguration(RNG4,Umu);
|
||||
// Umu=Zero();
|
||||
|
||||
RealD mass=0.1;
|
||||
RealD M5=1.8;
|
||||
|
||||
DomainWallFermionD Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5);
|
||||
|
||||
const int nbasis = 4;
|
||||
const int cb = 0 ;
|
||||
LatticeFermion prom(FGrid);
|
||||
|
||||
std::vector<LatticeFermion> subspace(nbasis,FGrid);
|
||||
|
||||
std::cout<<GridLogMessage<<"Calling Aggregation class" <<std::endl;
|
||||
|
||||
///////////////////////////////////////////////////////////
|
||||
// Squared operator is in HermOp
|
||||
///////////////////////////////////////////////////////////
|
||||
MdagMLinearOperator<DomainWallFermionD,LatticeFermion> HermDefOp(Ddwf);
|
||||
|
||||
///////////////////////////////////////////////////
|
||||
// Random aggregation space
|
||||
///////////////////////////////////////////////////
|
||||
std::cout<<GridLogMessage << "Building random aggregation class"<< std::endl;
|
||||
typedef Aggregation<vSpinColourVector,vTComplex,nbasis> Subspace;
|
||||
Subspace Aggregates(Coarse5d,FGrid,cb);
|
||||
Aggregates.CreateSubspaceRandom(RNG5);
|
||||
|
||||
///////////////////////////////////////////////////
|
||||
// Build little dirac op
|
||||
///////////////////////////////////////////////////
|
||||
std::cout<<GridLogMessage << "Building little Dirac operator"<< std::endl;
|
||||
|
||||
typedef GeneralCoarsenedMatrix<vSpinColourVector,vTComplex,nbasis> LittleDiracOperator;
|
||||
typedef LittleDiracOperator::CoarseVector CoarseVector;
|
||||
|
||||
NextToNearestStencilGeometry5D geom(Coarse5d);
|
||||
LittleDiracOperator LittleDiracOp(geom,FGrid,Coarse5d);
|
||||
LittleDiracOperator LittleDiracOpCol(geom,FGrid,Coarse5d);
|
||||
|
||||
HermOpAdaptor<LatticeFermionD> HOA(HermDefOp);
|
||||
|
||||
int pp=16;
|
||||
// LittleDiracOpCol.CoarsenOperator(HOA,Aggregates);
|
||||
// std::cout << "LittleDiracOp old " << LittleDiracOpCol._A[pp]<<std::endl;
|
||||
LittleDiracOp.CoarsenOperatorColoured(HOA,Aggregates);
|
||||
// std::cout << "LittleDiracOp new " << LittleDiracOp._A[pp]<<std::endl;
|
||||
|
||||
///////////////////////////////////////////////////
|
||||
// Test the operator
|
||||
///////////////////////////////////////////////////
|
||||
CoarseVector c_src (Coarse5d);
|
||||
CoarseVector c_res (Coarse5d);
|
||||
CoarseVector c_res_dag(Coarse5d);
|
||||
CoarseVector c_proj(Coarse5d);
|
||||
|
||||
subspace=Aggregates.subspace;
|
||||
|
||||
// random(CRNG,c_src);
|
||||
c_src = 1.0;
|
||||
|
||||
blockPromote(c_src,err,subspace);
|
||||
|
||||
prom=Zero();
|
||||
for(int b=0;b<nbasis;b++){
|
||||
prom=prom+subspace[b];
|
||||
}
|
||||
err=err-prom;
|
||||
std::cout<<GridLogMessage<<"Promoted back from subspace: err "<<norm2(err)<<std::endl;
|
||||
std::cout<<GridLogMessage<<"c_src "<<norm2(c_src)<<std::endl;
|
||||
std::cout<<GridLogMessage<<"prom "<<norm2(prom)<<std::endl;
|
||||
|
||||
HermDefOp.HermOp(prom,tmp);
|
||||
|
||||
blockProject(c_proj,tmp,subspace);
|
||||
std::cout<<GridLogMessage<<" Called Big Dirac Op "<<norm2(tmp)<<std::endl;
|
||||
|
||||
LittleDiracOp.M(c_src,c_res);
|
||||
LittleDiracOp.Mdag(c_src,c_res_dag);
|
||||
|
||||
std::cout<<GridLogMessage<<"Little dop : "<<norm2(c_res)<<std::endl;
|
||||
std::cout<<GridLogMessage<<"Little dop dag : "<<norm2(c_res_dag)<<std::endl;
|
||||
std::cout<<GridLogMessage<<"Big dop in subspace : "<<norm2(c_proj)<<std::endl;
|
||||
|
||||
c_proj = c_proj - c_res;
|
||||
std::cout<<GridLogMessage<<" ldop error: "<<norm2(c_proj)<<std::endl;
|
||||
|
||||
c_res_dag = c_res_dag - c_res;
|
||||
std::cout<<GridLogMessage<<"Little dopDag - dop: "<<norm2(c_res_dag)<<std::endl;
|
||||
|
||||
std::cout<<GridLogMessage << "Testing Hermiticity stochastically "<< std::endl;
|
||||
CoarseVector phi(Coarse5d);
|
||||
CoarseVector chi(Coarse5d);
|
||||
CoarseVector Aphi(Coarse5d);
|
||||
CoarseVector Achi(Coarse5d);
|
||||
|
||||
random(CRNG,phi);
|
||||
random(CRNG,chi);
|
||||
|
||||
std::cout<<GridLogMessage<<"Made randoms "<<norm2(phi)<<" " << norm2(chi)<<std::endl;
|
||||
|
||||
LittleDiracOp.M(phi,Aphi);
|
||||
|
||||
LittleDiracOp.Mdag(chi,Achi);
|
||||
|
||||
std::cout<<GridLogMessage<<"Aphi "<<norm2(Aphi)<<" A chi" << norm2(Achi)<<std::endl;
|
||||
|
||||
ComplexD pAc = innerProduct(chi,Aphi);
|
||||
ComplexD cAp = innerProduct(phi,Achi);
|
||||
ComplexD cAc = innerProduct(chi,Achi);
|
||||
ComplexD pAp = innerProduct(phi,Aphi);
|
||||
|
||||
std::cout<<GridLogMessage<< "pAc "<<pAc<<" cAp "<< cAp<< " diff "<<pAc-adj(cAp)<<std::endl;
|
||||
std::cout<<GridLogMessage<< "pAp "<<pAp<<" cAc "<< cAc<<"Should be real"<< std::endl;
|
||||
|
||||
std::cout<<GridLogMessage<<"Testing linearity"<<std::endl;
|
||||
CoarseVector PhiPlusChi(Coarse5d);
|
||||
CoarseVector APhiPlusChi(Coarse5d);
|
||||
CoarseVector linerr(Coarse5d);
|
||||
PhiPlusChi = phi+chi;
|
||||
LittleDiracOp.M(PhiPlusChi,APhiPlusChi);
|
||||
|
||||
linerr= APhiPlusChi-Aphi;
|
||||
linerr= linerr-Achi;
|
||||
std::cout<<GridLogMessage<<"**Diff "<<norm2(linerr)<<std::endl;
|
||||
|
||||
std::cout<<GridLogMessage<<std::endl;
|
||||
std::cout<<GridLogMessage<<std::endl;
|
||||
std::cout<<GridLogMessage<<"*******************************************"<<std::endl;
|
||||
std::cout<<GridLogMessage<<"*******************************************"<<std::endl;
|
||||
std::cout<<GridLogMessage<<"*******************************************"<<std::endl;
|
||||
|
||||
Grid_finalize();
|
||||
return 0;
|
||||
}
|
||||
@ -1,272 +0,0 @@
|
||||
/*************************************************************************************
|
||||
|
||||
Grid physics library, www.github.com/paboyle/Grid
|
||||
|
||||
Source file: ./tests/Test_general_coarse_hdcg.cc
|
||||
|
||||
Copyright (C) 2023
|
||||
|
||||
Author: Peter Boyle <pboyle@bnl.gov>
|
||||
|
||||
This program is free software; you can redistribute it and/or modify
|
||||
it under the terms of the GNU General Public License as published by
|
||||
the Free Software Foundation; either version 2 of the License, or
|
||||
(at your option) any later version.
|
||||
|
||||
This program is distributed in the hope that it will be useful,
|
||||
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||
GNU General Public License for more details.
|
||||
|
||||
You should have received a copy of the GNU General Public License along
|
||||
with this program; if not, write to the Free Software Foundation, Inc.,
|
||||
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
|
||||
|
||||
See the full license in the file "LICENSE" in the top level distribution directory
|
||||
*************************************************************************************/
|
||||
/* END LEGAL */
|
||||
#include <Grid/Grid.h>
|
||||
#include <Grid/lattice/PaddedCell.h>
|
||||
#include <Grid/stencil/GeneralLocalStencil.h>
|
||||
#include <Grid/algorithms/GeneralCoarsenedMatrix.h>
|
||||
#include <Grid/algorithms/iterative/AdefGeneric.h>
|
||||
|
||||
using namespace std;
|
||||
using namespace Grid;
|
||||
|
||||
template<class Field> class TestSolver : public LinearFunction<Field> {
|
||||
public:
|
||||
TestSolver() {};
|
||||
void operator() (const Field &in, Field &out){ out = Zero(); }
|
||||
};
|
||||
|
||||
|
||||
RealD InverseApproximation(RealD x){
|
||||
return 1.0/x;
|
||||
}
|
||||
|
||||
// Want Op in CoarsenOp to call MatPcDagMatPc
|
||||
template<class Field>
|
||||
class HermOpAdaptor : public LinearOperatorBase<Field>
|
||||
{
|
||||
LinearOperatorBase<Field> & wrapped;
|
||||
public:
|
||||
HermOpAdaptor(LinearOperatorBase<Field> &wrapme) : wrapped(wrapme) {};
|
||||
void Op (const Field &in, Field &out) { wrapped.HermOp(in,out); }
|
||||
void HermOp(const Field &in, Field &out) { wrapped.HermOp(in,out); }
|
||||
void AdjOp (const Field &in, Field &out){ wrapped.HermOp(in,out); }
|
||||
void OpDiag (const Field &in, Field &out) { assert(0); }
|
||||
void OpDir (const Field &in, Field &out,int dir,int disp) { assert(0); }
|
||||
void OpDirAll (const Field &in, std::vector<Field> &out) { assert(0); };
|
||||
void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){ assert(0); }
|
||||
};
|
||||
template<class Field,class Matrix> class ChebyshevSmoother : public LinearFunction<Field>
|
||||
{
|
||||
public:
|
||||
using LinearFunction<Field>::operator();
|
||||
typedef LinearOperatorBase<Field> FineOperator;
|
||||
FineOperator & _SmootherOperator;
|
||||
Chebyshev<Field> Cheby;
|
||||
ChebyshevSmoother(RealD _lo,RealD _hi,int _ord, FineOperator &SmootherOperator) :
|
||||
_SmootherOperator(SmootherOperator),
|
||||
Cheby(_lo,_hi,_ord,InverseApproximation)
|
||||
{
|
||||
std::cout << GridLogMessage<<" Chebyshev smoother order "<<_ord<<" ["<<_lo<<","<<_hi<<"]"<<std::endl;
|
||||
};
|
||||
void operator() (const Field &in, Field &out)
|
||||
{
|
||||
Field tmp(in.Grid());
|
||||
tmp = in;
|
||||
Cheby(_SmootherOperator,tmp,out);
|
||||
}
|
||||
};
|
||||
|
||||
int main (int argc, char ** argv)
|
||||
{
|
||||
Grid_init(&argc,&argv);
|
||||
|
||||
const int Ls=16;
|
||||
|
||||
GridCartesian * UGrid = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(),
|
||||
GridDefaultSimd(Nd,vComplex::Nsimd()),
|
||||
GridDefaultMpi());
|
||||
GridRedBlackCartesian * UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
|
||||
GridCartesian * FGrid = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid);
|
||||
GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid);
|
||||
|
||||
// Construct a coarsened grid with 4^4 cell
|
||||
Coordinate clatt = GridDefaultLatt();
|
||||
for(int d=0;d<clatt.size();d++){
|
||||
clatt[d] = clatt[d]/4;
|
||||
}
|
||||
GridCartesian *Coarse4d = SpaceTimeGrid::makeFourDimGrid(clatt,
|
||||
GridDefaultSimd(Nd,vComplex::Nsimd()),
|
||||
GridDefaultMpi());;
|
||||
GridCartesian *Coarse5d = SpaceTimeGrid::makeFiveDimGrid(1,Coarse4d);
|
||||
|
||||
///////////////////////// RNGs /////////////////////////////////
|
||||
std::vector<int> seeds4({1,2,3,4});
|
||||
std::vector<int> seeds5({5,6,7,8});
|
||||
std::vector<int> cseeds({5,6,7,8});
|
||||
|
||||
GridParallelRNG RNG5(FGrid); RNG5.SeedFixedIntegers(seeds5);
|
||||
GridParallelRNG RNG4(UGrid); RNG4.SeedFixedIntegers(seeds4);
|
||||
GridParallelRNG CRNG(Coarse5d);CRNG.SeedFixedIntegers(cseeds);
|
||||
|
||||
///////////////////////// Configuration /////////////////////////////////
|
||||
LatticeGaugeField Umu(UGrid);
|
||||
|
||||
FieldMetaData header;
|
||||
std::string file("ckpoint_lat.4000");
|
||||
NerscIO::readConfiguration(Umu,header,file);
|
||||
|
||||
//////////////////////// Fermion action //////////////////////////////////
|
||||
RealD mass=0.01;
|
||||
RealD M5=1.8;
|
||||
RealD b=1.5;
|
||||
RealD c=0.5;
|
||||
MobiusFermionD Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,b,c);
|
||||
|
||||
SchurDiagMooeeOperator<MobiusFermionD, LatticeFermion> HermOpEO(Ddwf);
|
||||
|
||||
typedef HermOpAdaptor<LatticeFermionD> HermFineMatrix;
|
||||
HermFineMatrix FineHermOp(HermOpEO);
|
||||
|
||||
LatticeFermion result(FrbGrid); result=Zero();
|
||||
|
||||
LatticeFermion src(FrbGrid); random(RNG5,src);
|
||||
|
||||
// Run power method on FineHermOp
|
||||
PowerMethod<LatticeFermion> PM; PM(HermOpEO,src);
|
||||
|
||||
|
||||
////////////////////////////////////////////////////////////
|
||||
///////////// Coarse basis and Little Dirac Operator ///////
|
||||
////////////////////////////////////////////////////////////
|
||||
const int nbasis = 40;
|
||||
const int cb = 0 ;
|
||||
typedef GeneralCoarsenedMatrix<vSpinColourVector,vTComplex,nbasis> LittleDiracOperator;
|
||||
typedef LittleDiracOperator::CoarseVector CoarseVector;
|
||||
|
||||
NextToNextToNextToNearestStencilGeometry5D geom(Coarse5d);
|
||||
|
||||
// Warning: This routine calls PVdagM.Op, not PVdagM.HermOp
|
||||
typedef Aggregation<vSpinColourVector,vTComplex,nbasis> Subspace;
|
||||
Subspace Aggregates(Coarse5d,FrbGrid,cb);
|
||||
Aggregates.CreateSubspaceChebyshev(RNG5,
|
||||
HermOpEO,
|
||||
nbasis,
|
||||
// 100.0,
|
||||
// 0.1, // Low pass is pretty high still -- 311 iters
|
||||
// 250.0,
|
||||
// 0.01, // subspace too low filter power wrong
|
||||
// 250.0,
|
||||
// 0.2, // slower
|
||||
95.0,
|
||||
// 0.05, // nbasis 12 - 311 -- wrong coarse inv
|
||||
// 0.05, // nbasis 12 - 154 -- right filt
|
||||
// 0.1, // nbasis 12 - 169 oops
|
||||
// 0.05, // nbasis 16 -- 127 iters
|
||||
// 0.03, // nbasis 16 -- 13-
|
||||
// 0.1, // nbasis 16 -- 142; sloppy solve
|
||||
0.1, // nbasis 24
|
||||
300);
|
||||
////////////////////////////////////////////////////////////
|
||||
// Need to check about red-black grid coarsening
|
||||
////////////////////////////////////////////////////////////
|
||||
LittleDiracOperator LittleDiracOp(geom,FrbGrid,Coarse5d);
|
||||
LittleDiracOp.CoarsenOperatorColoured(FineHermOp,Aggregates);
|
||||
|
||||
typedef HermitianLinearOperator<LittleDiracOperator,CoarseVector> HermMatrix;
|
||||
HermMatrix CoarseOp (LittleDiracOp);
|
||||
|
||||
//////////////////////////////////////////
|
||||
// Build a coarse lanczos
|
||||
//////////////////////////////////////////
|
||||
Chebyshev<CoarseVector> IRLCheby(0.02,50.0,71); // 1 iter
|
||||
FunctionHermOp<CoarseVector> IRLOpCheby(IRLCheby,CoarseOp);
|
||||
PlainHermOp<CoarseVector> IRLOp (CoarseOp);
|
||||
int Nk=64;
|
||||
int Nm=128;
|
||||
int Nstop=Nk;
|
||||
ImplicitlyRestartedLanczos<CoarseVector> IRL(IRLOpCheby,IRLOp,Nstop,Nk,Nm,1.0e-5,20);
|
||||
|
||||
int Nconv;
|
||||
std::vector<RealD> eval(Nm);
|
||||
std::vector<CoarseVector> evec(Nm,Coarse5d);
|
||||
CoarseVector c_src(Coarse5d); c_src=1.0;
|
||||
IRL.calc(eval,evec,c_src,Nconv);
|
||||
DeflatedGuesser<CoarseVector> DeflCoarseGuesser(evec,eval);
|
||||
|
||||
//////////////////////////////////////////
|
||||
// Build a coarse space solver
|
||||
//////////////////////////////////////////
|
||||
int maxit=20000;
|
||||
ConjugateGradient<CoarseVector> CG(1.0e-8,maxit,false);
|
||||
ConjugateGradient<LatticeFermionD> CGfine(1.0e-8,10000,false);
|
||||
ZeroGuesser<CoarseVector> CoarseZeroGuesser;
|
||||
|
||||
// HPDSolver<CoarseVector> HPDSolve(CoarseOp,CG,CoarseZeroGuesser);
|
||||
HPDSolver<CoarseVector> HPDSolve(CoarseOp,CG,DeflCoarseGuesser);
|
||||
|
||||
//////////////////////////////////////////
|
||||
// Build a smoother
|
||||
//////////////////////////////////////////
|
||||
// ChebyshevSmoother<LatticeFermionD,HermFineMatrix > Smoother(10.0,100.0,10,FineHermOp); //499
|
||||
// ChebyshevSmoother<LatticeFermionD,HermFineMatrix > Smoother(3.0,100.0,10,FineHermOp); //383
|
||||
// ChebyshevSmoother<LatticeFermionD,HermFineMatrix > Smoother(1.0,100.0,10,FineHermOp); //328
|
||||
// std::vector<RealD> los({0.5,1.0,3.0}); // 147/142/146 nbasis 1
|
||||
// std::vector<RealD> los({1.0,2.0}); // Nbasis 24: 88,86 iterations
|
||||
// std::vector<RealD> los({2.0,4.0}); // Nbasis 32 == 52, iters
|
||||
// std::vector<RealD> los({2.0,4.0}); // Nbasis 40 == 36,36 iters
|
||||
|
||||
//
|
||||
// Turns approx 2700 iterations into 340 fine multiplies with Nbasis 40
|
||||
// Need to measure cost of coarse space.
|
||||
//
|
||||
// -- i) Reduce coarse residual -- 0.04
|
||||
// -- ii) Lanczos on coarse space -- done
|
||||
// -- iii) Possible 1 hop project and/or preconditioning it - easy - PrecCG it and
|
||||
// use a limited stencil. Reread BFM code to check on evecs / deflation strategy with prec
|
||||
//
|
||||
std::vector<RealD> los({3.0}); // Nbasis 40 == 36,36 iters
|
||||
std::vector<int> ords({7,8,10}); // Nbasis 40 == 40,38,36 iters (320,342,396 mults)
|
||||
|
||||
// Standard CG
|
||||
// result=Zero();
|
||||
// CGfine(HermOpEO, src, result);
|
||||
|
||||
for(int l=0;l<los.size();l++){
|
||||
|
||||
RealD lo = los[l];
|
||||
|
||||
for(int o=0;o<ords.size();o++){
|
||||
|
||||
ConjugateGradient<CoarseVector> CGsloppy(4.0e-2,maxit,false);
|
||||
HPDSolver<CoarseVector> HPDSolveSloppy(CoarseOp,CGsloppy,DeflCoarseGuesser);
|
||||
|
||||
// ChebyshevSmoother<LatticeFermionD,HermFineMatrix > Smoother(lo,92,10,FineHermOp); // 36 best case
|
||||
ChebyshevSmoother<LatticeFermionD,HermFineMatrix > Smoother(lo,92,ords[o],FineHermOp); // 311
|
||||
|
||||
//////////////////////////////////////////
|
||||
// Build a HDCG solver
|
||||
//////////////////////////////////////////
|
||||
TwoLevelFlexiblePcg<LatticeFermion,CoarseVector,Subspace>
|
||||
HDCG(1.0e-8, 3000,
|
||||
FineHermOp,
|
||||
Smoother,
|
||||
HPDSolveSloppy,
|
||||
HPDSolve,
|
||||
Aggregates);
|
||||
|
||||
// result=Zero();
|
||||
// HDCG(src,result);
|
||||
|
||||
result=Zero();
|
||||
HDCG.Inflexible(src,result);
|
||||
}
|
||||
}
|
||||
|
||||
Grid_finalize();
|
||||
return 0;
|
||||
}
|
||||
@ -1,268 +0,0 @@
|
||||
/*************************************************************************************
|
||||
|
||||
Grid physics library, www.github.com/paboyle/Grid
|
||||
|
||||
Source file: ./tests/Test_padded_cell.cc
|
||||
|
||||
Copyright (C) 2023
|
||||
|
||||
Author: Peter Boyle <paboyle@ph.ed.ac.uk>
|
||||
|
||||
This program is free software; you can redistribute it and/or modify
|
||||
it under the terms of the GNU General Public License as published by
|
||||
the Free Software Foundation; either version 2 of the License, or
|
||||
(at your option) any later version.
|
||||
|
||||
This program is distributed in the hope that it will be useful,
|
||||
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||
GNU General Public License for more details.
|
||||
|
||||
You should have received a copy of the GNU General Public License along
|
||||
with this program; if not, write to the Free Software Foundation, Inc.,
|
||||
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
|
||||
|
||||
See the full license in the file "LICENSE" in the top level distribution directory
|
||||
*************************************************************************************/
|
||||
/* END LEGAL */
|
||||
#include <Grid/Grid.h>
|
||||
#include <Grid/lattice/PaddedCell.h>
|
||||
#include <Grid/stencil/GeneralLocalStencil.h>
|
||||
#include <Grid/algorithms/GeneralCoarsenedMatrix.h>
|
||||
|
||||
#include <Grid/algorithms/iterative/PrecGeneralisedConjugateResidual.h>
|
||||
#include <Grid/algorithms/iterative/PrecGeneralisedConjugateResidualNonHermitian.h>
|
||||
#include <Grid/algorithms/iterative/BiCGSTAB.h>
|
||||
|
||||
using namespace std;
|
||||
using namespace Grid;
|
||||
|
||||
template<class Field>
|
||||
class HermOpAdaptor : public LinearOperatorBase<Field>
|
||||
{
|
||||
LinearOperatorBase<Field> & wrapped;
|
||||
public:
|
||||
HermOpAdaptor(LinearOperatorBase<Field> &wrapme) : wrapped(wrapme) {};
|
||||
void OpDiag (const Field &in, Field &out) { assert(0); }
|
||||
void OpDir (const Field &in, Field &out,int dir,int disp) { assert(0); }
|
||||
void OpDirAll (const Field &in, std::vector<Field> &out){ assert(0); };
|
||||
void Op (const Field &in, Field &out){
|
||||
wrapped.HermOp(in,out);
|
||||
}
|
||||
void AdjOp (const Field &in, Field &out){
|
||||
wrapped.HermOp(in,out);
|
||||
}
|
||||
void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){ assert(0); }
|
||||
void HermOp(const Field &in, Field &out){
|
||||
wrapped.HermOp(in,out);
|
||||
}
|
||||
|
||||
};
|
||||
|
||||
template<class Matrix,class Field>
|
||||
class PVdagMLinearOperator : public LinearOperatorBase<Field> {
|
||||
Matrix &_Mat;
|
||||
Matrix &_PV;
|
||||
public:
|
||||
PVdagMLinearOperator(Matrix &Mat,Matrix &PV): _Mat(Mat),_PV(PV){};
|
||||
|
||||
void OpDiag (const Field &in, Field &out) { assert(0); }
|
||||
void OpDir (const Field &in, Field &out,int dir,int disp) { assert(0); }
|
||||
void OpDirAll (const Field &in, std::vector<Field> &out){ assert(0); };
|
||||
void Op (const Field &in, Field &out){
|
||||
Field tmp(in.Grid());
|
||||
_Mat.M(in,tmp);
|
||||
_PV.Mdag(tmp,out);
|
||||
}
|
||||
void AdjOp (const Field &in, Field &out){
|
||||
Field tmp(in.Grid());
|
||||
_PV.M(tmp,out);
|
||||
_Mat.Mdag(in,tmp);
|
||||
}
|
||||
void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){ assert(0); }
|
||||
void HermOp(const Field &in, Field &out){
|
||||
std::cout << "HermOp"<<std::endl;
|
||||
Field tmp(in.Grid());
|
||||
_Mat.M(in,tmp);
|
||||
_PV.Mdag(tmp,out);
|
||||
_PV.M(out,tmp);
|
||||
_Mat.Mdag(tmp,out);
|
||||
std::cout << "HermOp done "<<norm2(out)<<std::endl;
|
||||
|
||||
}
|
||||
};
|
||||
|
||||
template<class Field> class DumbOperator : public LinearOperatorBase<Field> {
|
||||
public:
|
||||
LatticeComplex scale;
|
||||
DumbOperator(GridBase *grid) : scale(grid)
|
||||
{
|
||||
scale = 0.0;
|
||||
LatticeComplex scalesft(grid);
|
||||
LatticeComplex scaletmp(grid);
|
||||
for(int d=0;d<4;d++){
|
||||
Lattice<iScalar<vInteger> > x(grid); LatticeCoordinate(x,d+1);
|
||||
LatticeCoordinate(scaletmp,d+1);
|
||||
scalesft = Cshift(scaletmp,d+1,1);
|
||||
scale = 100.0*scale + where( mod(x ,2)==(Integer)0, scalesft,scaletmp);
|
||||
}
|
||||
std::cout << " scale\n" << scale << std::endl;
|
||||
}
|
||||
// Support for coarsening to a multigrid
|
||||
void OpDiag (const Field &in, Field &out) {};
|
||||
void OpDir (const Field &in, Field &out,int dir,int disp){};
|
||||
void OpDirAll (const Field &in, std::vector<Field> &out) {};
|
||||
|
||||
void Op (const Field &in, Field &out){
|
||||
out = scale * in;
|
||||
}
|
||||
void AdjOp (const Field &in, Field &out){
|
||||
out = scale * in;
|
||||
}
|
||||
void HermOp(const Field &in, Field &out){
|
||||
double n1, n2;
|
||||
HermOpAndNorm(in,out,n1,n2);
|
||||
}
|
||||
void HermOpAndNorm(const Field &in, Field &out,double &n1,double &n2){
|
||||
ComplexD dot;
|
||||
|
||||
out = scale * in;
|
||||
|
||||
dot= innerProduct(in,out);
|
||||
n1=real(dot);
|
||||
|
||||
dot = innerProduct(out,out);
|
||||
n2=real(dot);
|
||||
}
|
||||
};
|
||||
|
||||
|
||||
int main (int argc, char ** argv)
|
||||
{
|
||||
Grid_init(&argc,&argv);
|
||||
|
||||
const int Ls=2;
|
||||
|
||||
GridCartesian * UGrid = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi());
|
||||
GridRedBlackCartesian * UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
|
||||
|
||||
GridCartesian * FGrid = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid);
|
||||
GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid);
|
||||
|
||||
// Construct a coarsened grid
|
||||
Coordinate clatt = GridDefaultLatt();
|
||||
for(int d=0;d<clatt.size();d++){
|
||||
clatt[d] = clatt[d]/4;
|
||||
}
|
||||
GridCartesian *Coarse4d = SpaceTimeGrid::makeFourDimGrid(clatt, GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi());;
|
||||
GridCartesian *Coarse5d = SpaceTimeGrid::makeFiveDimGrid(1,Coarse4d);
|
||||
|
||||
std::vector<int> seeds4({1,2,3,4});
|
||||
std::vector<int> seeds5({5,6,7,8});
|
||||
std::vector<int> cseeds({5,6,7,8});
|
||||
GridParallelRNG RNG5(FGrid); RNG5.SeedFixedIntegers(seeds5);
|
||||
GridParallelRNG RNG4(UGrid); RNG4.SeedFixedIntegers(seeds4);
|
||||
GridParallelRNG CRNG(Coarse5d);CRNG.SeedFixedIntegers(cseeds);
|
||||
|
||||
LatticeFermion src(FGrid); random(RNG5,src);
|
||||
LatticeFermion result(FGrid); result=Zero();
|
||||
LatticeFermion ref(FGrid); ref=Zero();
|
||||
LatticeFermion tmp(FGrid);
|
||||
LatticeFermion err(FGrid);
|
||||
LatticeGaugeField Umu(UGrid);
|
||||
|
||||
FieldMetaData header;
|
||||
std::string file("ckpoint_lat.4000");
|
||||
NerscIO::readConfiguration(Umu,header,file);
|
||||
//Umu = 1.0;
|
||||
|
||||
RealD mass=0.5;
|
||||
RealD M5=1.8;
|
||||
|
||||
DomainWallFermionD Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5);
|
||||
DomainWallFermionD Dpv(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,1.0,M5);
|
||||
|
||||
const int nbasis = 1;
|
||||
const int cb = 0 ;
|
||||
LatticeFermion prom(FGrid);
|
||||
|
||||
typedef GeneralCoarsenedMatrix<vSpinColourVector,vTComplex,nbasis> LittleDiracOperator;
|
||||
typedef LittleDiracOperator::CoarseVector CoarseVector;
|
||||
|
||||
NextToNearestStencilGeometry5D geom;
|
||||
|
||||
std::cout<<GridLogMessage<<std::endl;
|
||||
std::cout<<GridLogMessage<<"*******************************************"<<std::endl;
|
||||
std::cout<<GridLogMessage<<std::endl;
|
||||
|
||||
PVdagMLinearOperator<DomainWallFermionD,LatticeFermionD> PVdagM(Ddwf,Dpv);
|
||||
HermOpAdaptor<LatticeFermionD> HOA(PVdagM);
|
||||
|
||||
// Run power method on HOA??
|
||||
PowerMethod<LatticeFermion> PM; PM(HOA,src);
|
||||
|
||||
// Warning: This routine calls PVdagM.Op, not PVdagM.HermOp
|
||||
typedef Aggregation<vSpinColourVector,vTComplex,nbasis> Subspace;
|
||||
Subspace AggregatesPD(Coarse5d,FGrid,cb);
|
||||
AggregatesPD.CreateSubspaceChebyshev(RNG5,
|
||||
HOA,
|
||||
nbasis,
|
||||
5000.0,
|
||||
0.02,
|
||||
100,
|
||||
50,
|
||||
50,
|
||||
0.0);
|
||||
|
||||
LittleDiracOperator LittleDiracOpPV(geom,FGrid,Coarse5d);
|
||||
LittleDiracOpPV.CoarsenOperator(PVdagM,AggregatesPD);
|
||||
|
||||
std::cout<<GridLogMessage<<std::endl;
|
||||
std::cout<<GridLogMessage<<"*******************************************"<<std::endl;
|
||||
std::cout<<GridLogMessage<<std::endl;
|
||||
std::cout<<GridLogMessage<<"Testing coarsened operator "<<std::endl;
|
||||
|
||||
CoarseVector c_src (Coarse5d);
|
||||
CoarseVector c_res (Coarse5d);
|
||||
CoarseVector c_proj(Coarse5d);
|
||||
|
||||
std::vector<LatticeFermion> subspace(nbasis,FGrid);
|
||||
subspace=AggregatesPD.subspace;
|
||||
|
||||
Complex one(1.0);
|
||||
c_src = one; // 1 in every element for vector 1.
|
||||
blockPromote(c_src,err,subspace);
|
||||
|
||||
prom=Zero();
|
||||
for(int b=0;b<nbasis;b++){
|
||||
prom=prom+subspace[b];
|
||||
}
|
||||
err=err-prom;
|
||||
std::cout<<GridLogMessage<<"Promoted back from subspace: err "<<norm2(err)<<std::endl;
|
||||
std::cout<<GridLogMessage<<"c_src "<<norm2(c_src)<<std::endl;
|
||||
std::cout<<GridLogMessage<<"prom "<<norm2(prom)<<std::endl;
|
||||
|
||||
PVdagM.Op(prom,tmp);
|
||||
blockProject(c_proj,tmp,subspace);
|
||||
std::cout<<GridLogMessage<<" Called Big Dirac Op "<<norm2(tmp)<<std::endl;
|
||||
|
||||
LittleDiracOpPV.M(c_src,c_res);
|
||||
std::cout<<GridLogMessage<<" Called Little Dirac Op c_src "<< norm2(c_src) << " c_res "<< norm2(c_res) <<std::endl;
|
||||
|
||||
std::cout<<GridLogMessage<<"Little dop : "<<norm2(c_res)<<std::endl;
|
||||
// std::cout<<GridLogMessage<<" Little "<< c_res<<std::endl;
|
||||
|
||||
std::cout<<GridLogMessage<<"Big dop in subspace : "<<norm2(c_proj)<<std::endl;
|
||||
// std::cout<<GridLogMessage<<" Big "<< c_proj<<std::endl;
|
||||
c_proj = c_proj - c_res;
|
||||
std::cout<<GridLogMessage<<" ldop error: "<<norm2(c_proj)<<std::endl;
|
||||
// std::cout<<GridLogMessage<<" error "<< c_proj<<std::endl;
|
||||
|
||||
std::cout<<GridLogMessage<<std::endl;
|
||||
std::cout<<GridLogMessage<<"*******************************************"<<std::endl;
|
||||
std::cout<<GridLogMessage<<std::endl;
|
||||
std::cout<<GridLogMessage << "Done "<< std::endl;
|
||||
|
||||
Grid_finalize();
|
||||
return 0;
|
||||
}
|
||||
@ -1,7 +1,6 @@
|
||||
/*************************************************************************************
|
||||
|
||||
Grid physics library, www.github.com/paboyle/Grid
|
||||
|
||||
Source file: ./tests/qdpxx/Test_qdpxx_munprec.cc
|
||||
|
||||
Copyright (C) 2015
|
||||
@ -26,13 +25,17 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
|
||||
See the full license in the file "LICENSE" in the top level distribution directory
|
||||
*************************************************************************************/
|
||||
/* END LEGAL */
|
||||
#include <chroma.h>
|
||||
#include <actions/ferm/invert/syssolver_linop_cg_array.h>
|
||||
#include <actions/ferm/invert/syssolver_linop_aggregate.h>
|
||||
|
||||
#include <Grid/Grid.h>
|
||||
|
||||
int Ls=8;
|
||||
double M5=1.6;
|
||||
double mq=0.01;
|
||||
double zolo_lo = 0.1;
|
||||
double zolo_hi = 2.0;
|
||||
double zolo_lo = 0.01;
|
||||
double zolo_hi = 7.0;
|
||||
double mobius_scale=2.0;
|
||||
|
||||
enum ChromaAction {
|
||||
@ -55,11 +58,6 @@ enum ChromaAction {
|
||||
void calc_grid (ChromaAction action,Grid::LatticeGaugeField & lat, Grid::LatticeFermion &src, Grid::LatticeFermion &res,int dag);
|
||||
void calc_chroma (ChromaAction action,Grid::LatticeGaugeField & lat, Grid::LatticeFermion &src, Grid::LatticeFermion &res,int dag);
|
||||
|
||||
#include <chroma.h>
|
||||
#include <actions/ferm/invert/syssolver_linop_cg_array.h>
|
||||
#include <actions/ferm/invert/syssolver_linop_aggregate.h>
|
||||
|
||||
|
||||
|
||||
namespace Chroma {
|
||||
|
||||
@ -81,7 +79,7 @@ public:
|
||||
|
||||
std::vector<int> x(4);
|
||||
QDP::multi1d<int> cx(4);
|
||||
std::vector<int> gd= gr.Grid()->GlobalDimensions();
|
||||
Grid::Coordinate gd = gr.Grid()->GlobalDimensions();
|
||||
|
||||
for (x[0]=0;x[0]<gd[0];x[0]++){
|
||||
for (x[1]=0;x[1]<gd[1];x[1]++){
|
||||
@ -124,7 +122,7 @@ public:
|
||||
|
||||
std::vector<int> x(5);
|
||||
QDP::multi1d<int> cx(4);
|
||||
std::vector<int> gd= gr.Grid()->GlobalDimensions();
|
||||
Grid::Coordinate gd= gr.Grid()->GlobalDimensions();
|
||||
|
||||
for (x[0]=0;x[0]<gd[0];x[0]++){
|
||||
for (x[1]=0;x[1]<gd[1];x[1]++){
|
||||
@ -166,7 +164,7 @@ public:
|
||||
|
||||
std::vector<int> x(5);
|
||||
QDP::multi1d<int> cx(4);
|
||||
std::vector<int> gd= gr.Grid()->GlobalDimensions();
|
||||
Grid::Coordinate gd= gr.Grid()->GlobalDimensions();
|
||||
|
||||
for (x[0]=0;x[0]<gd[0];x[0]++){
|
||||
for (x[1]=0;x[1]<gd[1];x[1]++){
|
||||
@ -304,7 +302,30 @@ public:
|
||||
// param.approximation_type=COEFF_TYPE_TANH_UNSCALED;
|
||||
// param.approximation_type=COEFF_TYPE_TANH;
|
||||
param.tuning_strategy_xml=
|
||||
"<TuningStrategy><Name>OVEXT_CONSTANT_STRATEGY</Name></TuningStrategy>\n";
|
||||
"<TuningStrategy><Name>OVEXT_CONSTANT_STRATEGY</Name><TuningConstant>1.0</TuningConstant></TuningStrategy>\n";
|
||||
UnprecOvExtFermActArray S_f(cfs,param);
|
||||
Handle< FermState<T4,U,U> > fs( S_f.createState(u) );
|
||||
Handle< LinearOperatorArray<T4> > M(S_f.linOp(fs));
|
||||
return M;
|
||||
}
|
||||
if ( parms == HwPartFracTanh ) {
|
||||
if ( Ls%2 == 0 ) {
|
||||
printf("Ls is not odd\n");
|
||||
exit(-1);
|
||||
}
|
||||
UnprecOvExtFermActArrayParams param;
|
||||
param.OverMass=M5;
|
||||
param.Mass=_mq;
|
||||
param.RatPolyDeg = Ls;
|
||||
param.ApproxMin =eps_lo;
|
||||
param.ApproxMax =eps_hi;
|
||||
param.b5 =1.0;
|
||||
param.c5 =1.0;
|
||||
// param.approximation_type=COEFF_TYPE_ZOLOTAREV;
|
||||
param.approximation_type=COEFF_TYPE_TANH_UNSCALED;
|
||||
//param.approximation_type=COEFF_TYPE_TANH;
|
||||
param.tuning_strategy_xml=
|
||||
"<TuningStrategy><Name>OVEXT_CONSTANT_STRATEGY</Name><TuningConstant>1.0</TuningConstant></TuningStrategy>\n";
|
||||
UnprecOvExtFermActArray S_f(cfs,param);
|
||||
Handle< FermState<T4,U,U> > fs( S_f.createState(u) );
|
||||
Handle< LinearOperatorArray<T4> > M(S_f.linOp(fs));
|
||||
@ -316,7 +337,35 @@ public:
|
||||
param.ApproxMin=eps_lo;
|
||||
param.ApproxMax=eps_hi;
|
||||
param.approximation_type=COEFF_TYPE_ZOLOTAREV;
|
||||
param.RatPolyDeg=Ls;
|
||||
param.RatPolyDeg=Ls-1;
|
||||
// The following is why I think Chroma made some directional errors:
|
||||
param.AuxFermAct= std::string(
|
||||
"<AuxFermAct>\n"
|
||||
" <FermAct>UNPRECONDITIONED_WILSON</FermAct>\n"
|
||||
" <Mass>-1.8</Mass>\n"
|
||||
" <b5>1</b5>\n"
|
||||
" <c5>0</c5>\n"
|
||||
" <MaxCG>1000</MaxCG>\n"
|
||||
" <RsdCG>1.0e-9</RsdCG>\n"
|
||||
" <FermionBC>\n"
|
||||
" <FermBC>SIMPLE_FERMBC</FermBC>\n"
|
||||
" <boundary>1 1 1 1</boundary>\n"
|
||||
" </FermionBC> \n"
|
||||
"</AuxFermAct>"
|
||||
);
|
||||
param.AuxFermActGrp= std::string("");
|
||||
UnprecOvlapContFrac5DFermActArray S_f(fbc,param);
|
||||
Handle< FermState<T4,U,U> > fs( S_f.createState(u) );
|
||||
Handle< LinearOperatorArray<T4> > M(S_f.linOp(fs));
|
||||
return M;
|
||||
}
|
||||
if ( parms == HwContFracTanh ) {
|
||||
UnprecOvlapContFrac5DFermActParams param;
|
||||
param.Mass=_mq; // How is M5 set? Wilson mass In AuxFermAct
|
||||
param.ApproxMin=eps_lo;
|
||||
param.ApproxMax=eps_hi;
|
||||
param.approximation_type=COEFF_TYPE_TANH_UNSCALED;
|
||||
param.RatPolyDeg=Ls-1;
|
||||
// The following is why I think Chroma made some directional errors:
|
||||
param.AuxFermAct= std::string(
|
||||
"<AuxFermAct>\n"
|
||||
@ -378,7 +427,14 @@ int main (int argc,char **argv )
|
||||
* Setup QDP
|
||||
*********************************************************/
|
||||
Chroma::initialize(&argc,&argv);
|
||||
Chroma::WilsonTypeFermActs4DEnv::registerAll();
|
||||
// Chroma::WilsonTypeFermActs4DEnv::registerAll();
|
||||
Chroma::WilsonTypeFermActsEnv::registerAll();
|
||||
//bool linkageHack(void)
|
||||
//{
|
||||
// bool foo = true;
|
||||
// Inline Measurements
|
||||
// InlineAggregateEnv::registerAll();
|
||||
// GaugeInitEnv::registerAll();
|
||||
|
||||
/********************************************************
|
||||
* Setup Grid
|
||||
@ -388,26 +444,34 @@ int main (int argc,char **argv )
|
||||
Grid::GridDefaultSimd(Grid::Nd,Grid::vComplex::Nsimd()),
|
||||
Grid::GridDefaultMpi());
|
||||
|
||||
std::vector<int> gd = UGrid->GlobalDimensions();
|
||||
Grid::Coordinate gd = UGrid->GlobalDimensions();
|
||||
QDP::multi1d<int> nrow(QDP::Nd);
|
||||
for(int mu=0;mu<4;mu++) nrow[mu] = gd[mu];
|
||||
|
||||
QDP::Layout::setLattSize(nrow);
|
||||
QDP::Layout::create();
|
||||
|
||||
Grid::GridCartesian * FGrid = Grid::SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid);
|
||||
Grid::LatticeGaugeField lat(UGrid);
|
||||
Grid::LatticeFermion src(FGrid);
|
||||
Grid::LatticeFermion res_chroma(FGrid);
|
||||
Grid::LatticeFermion res_grid (FGrid);
|
||||
|
||||
std::vector<ChromaAction> ActionList({
|
||||
HtCayleyTanh, // Plain old DWF.
|
||||
HmCayleyTanh,
|
||||
HwCayleyTanh,
|
||||
HtCayleyZolo, // Plain old DWF.
|
||||
HmCayleyZolo,
|
||||
HwCayleyZolo
|
||||
HwCayleyZolo,
|
||||
HwPartFracZolo,
|
||||
HwContFracZolo,
|
||||
HwContFracTanh
|
||||
});
|
||||
std::vector<int> LsList({
|
||||
8,//HtCayleyTanh, // Plain old DWF.
|
||||
8,//HmCayleyTanh,
|
||||
8,//HwCayleyTanh,
|
||||
8,//HtCayleyZolo, // Plain old DWF.
|
||||
8,//HmCayleyZolo,
|
||||
8,//HwCayleyZolo,
|
||||
9,//HwPartFracZolo
|
||||
9, //HwContFracZolo
|
||||
9 //HwContFracTanh
|
||||
});
|
||||
std::vector<std::string> ActionName({
|
||||
"HtCayleyTanh",
|
||||
@ -415,10 +479,19 @@ int main (int argc,char **argv )
|
||||
"HwCayleyTanh",
|
||||
"HtCayleyZolo",
|
||||
"HmCayleyZolo",
|
||||
"HwCayleyZolo"
|
||||
"HwCayleyZolo",
|
||||
"HwPartFracZolo",
|
||||
"HwContFracZolo",
|
||||
"HwContFracTanh"
|
||||
});
|
||||
|
||||
for(int i=0;i<ActionList.size();i++) {
|
||||
Ls = LsList[i];
|
||||
Grid::GridCartesian * FGrid = Grid::SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid);
|
||||
Grid::LatticeGaugeField lat(UGrid);
|
||||
Grid::LatticeFermion src(FGrid);
|
||||
Grid::LatticeFermion res_chroma(FGrid);
|
||||
Grid::LatticeFermion res_grid (FGrid);
|
||||
std::cout << "*****************************"<<std::endl;
|
||||
std::cout << "Action "<<ActionName[i]<<std::endl;
|
||||
std::cout << "*****************************"<<std::endl;
|
||||
@ -439,6 +512,7 @@ int main (int argc,char **argv )
|
||||
|
||||
std::cout << "Norm of difference "<<Grid::norm2(res_chroma)<<std::endl;
|
||||
}
|
||||
delete FGrid;
|
||||
}
|
||||
|
||||
std::cout << "Finished test "<<std::endl;
|
||||
@ -502,7 +576,7 @@ void calc_grid(ChromaAction action,Grid::LatticeGaugeField & Umu, Grid::LatticeF
|
||||
Grid::gaussian(RNG5,src);
|
||||
Grid::gaussian(RNG5,res);
|
||||
|
||||
Grid::SU<Nc>::HotConfiguration(RNG4,Umu);
|
||||
Grid::SU<Grid::Nc>::HotConfiguration(RNG4,Umu);
|
||||
|
||||
/*
|
||||
Grid::LatticeColourMatrix U(UGrid);
|
||||
@ -519,7 +593,7 @@ void calc_grid(ChromaAction action,Grid::LatticeGaugeField & Umu, Grid::LatticeF
|
||||
|
||||
if ( action == HtCayleyTanh ) {
|
||||
|
||||
Grid::DomainWallFermionR Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,_mass,_M5);
|
||||
Grid::DomainWallFermionD Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,_mass,_M5);
|
||||
|
||||
std::cout << Grid::GridLogMessage <<" Calling domain wall multiply "<<std::endl;
|
||||
|
||||
@ -535,7 +609,7 @@ void calc_grid(ChromaAction action,Grid::LatticeGaugeField & Umu, Grid::LatticeF
|
||||
|
||||
Grid::Real _b = 0.5*(mobius_scale +1.0);
|
||||
Grid::Real _c = 0.5*(mobius_scale -1.0);
|
||||
Grid::MobiusZolotarevFermionR D(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,_mass,_M5,_b,_c,zolo_lo,zolo_hi);
|
||||
Grid::MobiusZolotarevFermionD D(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,_mass,_M5,_b,_c,zolo_lo,zolo_hi);
|
||||
|
||||
std::cout << Grid::GridLogMessage <<" Calling mobius zolo multiply "<<std::endl;
|
||||
|
||||
@ -549,7 +623,7 @@ void calc_grid(ChromaAction action,Grid::LatticeGaugeField & Umu, Grid::LatticeF
|
||||
|
||||
if ( action == HtCayleyZolo ) {
|
||||
|
||||
Grid::ShamirZolotarevFermionR D(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,_mass,_M5,zolo_lo,zolo_hi);
|
||||
Grid::ShamirZolotarevFermionD D(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,_mass,_M5,zolo_lo,zolo_hi);
|
||||
|
||||
std::cout << Grid::GridLogMessage <<" Calling shamir zolo multiply "<<std::endl;
|
||||
|
||||
@ -561,6 +635,60 @@ void calc_grid(ChromaAction action,Grid::LatticeGaugeField & Umu, Grid::LatticeF
|
||||
return;
|
||||
}
|
||||
|
||||
if ( action == HwPartFracTanh ) {
|
||||
|
||||
Grid::OverlapWilsonPartialFractionTanhFermionD Dov(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,_mass,_M5,1.0);
|
||||
|
||||
std::cout << Grid::GridLogMessage <<" Calling part frac tanh multiply "<<std::endl;
|
||||
|
||||
if ( dag )
|
||||
Dov.Mdag(src,res);
|
||||
else
|
||||
Dov.M(src,res);
|
||||
|
||||
return;
|
||||
}
|
||||
|
||||
if ( action == HwContFracTanh ) {
|
||||
|
||||
Grid::OverlapWilsonContFracTanhFermionD Dov(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,_mass,_M5,1.0);
|
||||
|
||||
std::cout << Grid::GridLogMessage <<" Calling cont frac tanh multiply "<<std::endl;
|
||||
|
||||
if ( dag )
|
||||
Dov.Mdag(src,res);
|
||||
else
|
||||
Dov.M(src,res);
|
||||
|
||||
return;
|
||||
}
|
||||
if ( action == HwContFracZolo ) {
|
||||
|
||||
Grid::OverlapWilsonContFracZolotarevFermionD Dov(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,_mass,_M5,zolo_lo,zolo_hi);
|
||||
|
||||
std::cout << Grid::GridLogMessage <<" Calling cont frac zolo multiply "<<std::endl;
|
||||
|
||||
if ( dag )
|
||||
Dov.Mdag(src,res);
|
||||
else
|
||||
Dov.M(src,res);
|
||||
|
||||
return;
|
||||
}
|
||||
|
||||
if ( action == HwPartFracZolo ) {
|
||||
|
||||
Grid::OverlapWilsonPartialFractionZolotarevFermionD Dov(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,_mass,_M5,zolo_lo,zolo_hi);
|
||||
std::cout << Grid::GridLogMessage <<" Calling part frac zolotarev multiply "<<std::endl;
|
||||
|
||||
if ( dag )
|
||||
Dov.Mdag(src,res);
|
||||
else
|
||||
Dov.M(src,res);
|
||||
|
||||
return;
|
||||
}
|
||||
|
||||
/*
|
||||
if ( action == HmCayleyTanh ) {
|
||||
Grid::Real _b = 0.5*(mobius_scale +1.0);
|
||||
@ -581,7 +709,7 @@ void calc_grid(ChromaAction action,Grid::LatticeGaugeField & Umu, Grid::LatticeF
|
||||
|
||||
if ( action == HmCayleyTanh ) {
|
||||
|
||||
Grid::ScaledShamirFermionR D(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,_mass,_M5,mobius_scale);
|
||||
Grid::ScaledShamirFermionD D(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,_mass,_M5,mobius_scale);
|
||||
|
||||
std::cout << Grid::GridLogMessage <<" Calling scaled shamir multiply "<<std::endl;
|
||||
|
||||
@ -595,7 +723,7 @@ void calc_grid(ChromaAction action,Grid::LatticeGaugeField & Umu, Grid::LatticeF
|
||||
|
||||
if ( action == HwCayleyTanh ) {
|
||||
|
||||
Grid::OverlapWilsonCayleyTanhFermionR D(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,_mass,_M5,1.0);
|
||||
Grid::OverlapWilsonCayleyTanhFermionD D(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,_mass,_M5,1.0);
|
||||
|
||||
if ( dag )
|
||||
D.Mdag(src,res);
|
||||
@ -607,7 +735,7 @@ void calc_grid(ChromaAction action,Grid::LatticeGaugeField & Umu, Grid::LatticeF
|
||||
|
||||
if ( action == HwCayleyZolo ) {
|
||||
|
||||
Grid::OverlapWilsonCayleyZolotarevFermionR D(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,_mass,_M5,zolo_lo,zolo_hi);
|
||||
Grid::OverlapWilsonCayleyZolotarevFermionD D(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,_mass,_M5,zolo_lo,zolo_hi);
|
||||
|
||||
if ( dag )
|
||||
D.Mdag(src,res);
|
||||
|
||||
@ -1,4 +1,4 @@
|
||||
/*************************************************************************************
|
||||
*************************************************************************************
|
||||
|
||||
Grid physics library, www.github.com/paboyle/Grid
|
||||
|
||||
@ -67,7 +67,13 @@ int main(int argc, char** argv) {
|
||||
result = Zero();
|
||||
LatticeGaugeField Umu(UGrid);
|
||||
|
||||
#if 0
|
||||
FieldMetaData header;
|
||||
std::string file("ckpoint_lat.4000");
|
||||
NerscIO::readConfiguration(Umu,header,file);
|
||||
#else
|
||||
SU<Nc>::HotConfiguration(RNG4, Umu);
|
||||
#endif
|
||||
|
||||
std::cout << GridLogMessage << "Lattice dimensions: " << GridDefaultLatt()
|
||||
<< " Ls: " << Ls << std::endl;
|
||||
|
||||
@ -54,15 +54,30 @@ int main (int argc, char ** argv)
|
||||
GridCartesian * FGrid = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid);
|
||||
GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid);
|
||||
|
||||
std::vector<ComplexD> qmu;
|
||||
qmu.push_back(ComplexD(0.1,0.0));
|
||||
qmu.push_back(ComplexD(0.0,0.0));
|
||||
qmu.push_back(ComplexD(0.0,0.0));
|
||||
qmu.push_back(ComplexD(0.0,0.01));
|
||||
|
||||
|
||||
std::vector<int> seeds4({1,2,3,4});
|
||||
std::vector<int> seeds5({5,6,7,8});
|
||||
GridParallelRNG RNG5(FGrid); RNG5.SeedFixedIntegers(seeds5);
|
||||
GridParallelRNG RNG4(UGrid); RNG4.SeedFixedIntegers(seeds4);
|
||||
|
||||
LatticeFermion tmp(FGrid);
|
||||
LatticeFermion src(FGrid); random(RNG5,src);
|
||||
LatticeFermion result(FGrid); result=Zero();
|
||||
LatticeGaugeField Umu(UGrid); SU<Nc>::HotConfiguration(RNG4,Umu);
|
||||
|
||||
LatticeGaugeField Umu(UGrid);
|
||||
#if 0
|
||||
FieldMetaData header;
|
||||
std::string file("ckpoint_lat.4000");
|
||||
NerscIO::readConfiguration(Umu,header,file);
|
||||
#else
|
||||
SU<Nc>::HotConfiguration(RNG4,Umu);
|
||||
#endif
|
||||
|
||||
std::vector<LatticeColourMatrix> U(4,UGrid);
|
||||
for(int mu=0;mu<Nd;mu++){
|
||||
U[mu] = PeekIndex<LorentzIndex>(Umu,mu);
|
||||
@ -71,8 +86,15 @@ int main (int argc, char ** argv)
|
||||
RealD mass=0.1;
|
||||
RealD M5=1.8;
|
||||
DomainWallFermionD Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5);
|
||||
Ddwf.qmu = qmu;
|
||||
|
||||
Ddwf.M(src,tmp);
|
||||
std::cout << " |M src|^2 "<<norm2(tmp)<<std::endl;
|
||||
MdagMLinearOperator<DomainWallFermionD,LatticeFermion> HermOp(Ddwf);
|
||||
HermOp.HermOp(src,tmp);
|
||||
|
||||
std::cout << " <src|MdagM| src> "<<innerProduct(src,tmp)<<std::endl;
|
||||
|
||||
ConjugateGradient<LatticeFermion> CG(1.0e-6,10000);
|
||||
CG(HermOp,src,result);
|
||||
|
||||
|
||||
Reference in New Issue
Block a user