mirror of
https://github.com/paboyle/Grid.git
synced 2025-07-01 05:57:07 +01:00
Compare commits
3 Commits
2b43308208
...
feature/bo
| Author | SHA1 | Date | |
|---|---|---|---|
| 6815e138b4 | |||
| e29b97b3ea | |||
| ad2b699d2b |
@ -158,18 +158,6 @@ 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) {
|
||||
|
||||
RealD scale;
|
||||
|
||||
@ -1,431 +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);
|
||||
|
||||
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 npoint;
|
||||
std::vector<Coordinate> shifts;
|
||||
virtual void BuildShifts(void) { assert(0); } ;
|
||||
int Depth(void){return depth;};
|
||||
NonLocalStencilGeometry(int _depth) : depth(_depth)
|
||||
{
|
||||
};
|
||||
virtual ~NonLocalStencilGeometry() {};
|
||||
};
|
||||
// Need to worry about red-black now
|
||||
class NextToNearestStencilGeometry4D : public NonLocalStencilGeometry {
|
||||
public:
|
||||
NextToNearestStencilGeometry4D(void) : NonLocalStencilGeometry(2)
|
||||
{
|
||||
this->BuildShifts();
|
||||
};
|
||||
virtual ~NextToNearestStencilGeometry4D() {};
|
||||
virtual void BuildShifts(void)
|
||||
{
|
||||
this->shifts.resize(0);
|
||||
// Like HDCG: 81 point stencil including self connection
|
||||
this->shifts.push_back(Coordinate({0,0,0,0}));
|
||||
// +-x, +-y, +-z, +-t : 8
|
||||
for(int s=-1;s<=1;s+=2){
|
||||
this->shifts.push_back(Coordinate({s,0,0,0}));
|
||||
this->shifts.push_back(Coordinate({0,s,0,0}));
|
||||
this->shifts.push_back(Coordinate({0,0,s,0}));
|
||||
this->shifts.push_back(Coordinate({0,0,0,s}));
|
||||
}
|
||||
// +-x+-y, +-x+-z, +-x+-t, +-y+-z, +-y+-t, +-z+-t : 24
|
||||
for(int s1=-1;s1<=1;s1+=2){
|
||||
for(int s2=-1;s2<=1;s2+=2){
|
||||
this->shifts.push_back(Coordinate({s1,s2,0,0}));
|
||||
this->shifts.push_back(Coordinate({s1,0,s2,0}));
|
||||
this->shifts.push_back(Coordinate({s1,0,0,s2}));
|
||||
this->shifts.push_back(Coordinate({0,s1,s2,0}));
|
||||
this->shifts.push_back(Coordinate({0,s1,0,s2}));
|
||||
this->shifts.push_back(Coordinate({0,0,s1,s2}));
|
||||
}}
|
||||
this->npoint = this->shifts.size();
|
||||
}
|
||||
};
|
||||
// Need to worry about red-black now
|
||||
class NextToNextToNextToNearestStencilGeometry4D : public NonLocalStencilGeometry {
|
||||
public:
|
||||
NextToNextToNextToNearestStencilGeometry4D(void) : NonLocalStencilGeometry(4)
|
||||
{
|
||||
this->BuildShifts();
|
||||
};
|
||||
virtual ~NextToNextToNextToNearestStencilGeometry4D() {}
|
||||
virtual void BuildShifts(void)
|
||||
{
|
||||
this->shifts.resize(0);
|
||||
// Like HDCG: 81 point stencil including self connection
|
||||
this->shifts.push_back(Coordinate({0,0,0,0}));
|
||||
// +-x, +-y, +-z, +-t : 8
|
||||
for(int s=-1;s<=1;s+=2){
|
||||
this->shifts.push_back(Coordinate({s,0,0,0}));
|
||||
this->shifts.push_back(Coordinate({0,s,0,0}));
|
||||
this->shifts.push_back(Coordinate({0,0,s,0}));
|
||||
this->shifts.push_back(Coordinate({0,0,0,s}));
|
||||
}
|
||||
// +-x+-y, +-x+-z, +-x+-t, +-y+-z, +-y+-t, +-z+-t : 24
|
||||
for(int s1=-1;s1<=1;s1+=2){
|
||||
for(int s2=-1;s2<=1;s2+=2){
|
||||
this->shifts.push_back(Coordinate({s1,s2,0,0}));
|
||||
this->shifts.push_back(Coordinate({s1,0,s2,0}));
|
||||
this->shifts.push_back(Coordinate({s1,0,0,s2}));
|
||||
this->shifts.push_back(Coordinate({0,s1,s2,0}));
|
||||
this->shifts.push_back(Coordinate({0,s1,0,s2}));
|
||||
this->shifts.push_back(Coordinate({0,0,s1,s2}));
|
||||
}}
|
||||
// +-x+-y+-z, +-x+-y+-z, +-x+-y+-z,
|
||||
for(int s1=-1;s1<=1;s1+=2){
|
||||
for(int s2=-1;s2<=1;s2+=2){
|
||||
for(int s3=-1;s3<=1;s3+=2){
|
||||
this->shifts.push_back(Coordinate({s1,s2,s3,0})); // 8x4 = 32
|
||||
this->shifts.push_back(Coordinate({s1,s2,0,s3}));
|
||||
this->shifts.push_back(Coordinate({s1,0,s2,s3}));
|
||||
this->shifts.push_back(Coordinate({0,s1,s2,s3}));
|
||||
}}}
|
||||
for(int s1=-1;s1<=1;s1+=2){
|
||||
for(int s2=-1;s2<=1;s2+=2){
|
||||
for(int s3=-1;s3<=1;s3+=2){
|
||||
for(int s4=-1;s4<=1;s4+=2){
|
||||
this->shifts.push_back(Coordinate({s1,s2,s3,s4})); // 16
|
||||
}}}}
|
||||
this->npoint = this->shifts.size();
|
||||
}
|
||||
};
|
||||
class NextToNearestStencilGeometry5D : public NonLocalStencilGeometry {
|
||||
public:
|
||||
NextToNearestStencilGeometry5D(void) : NonLocalStencilGeometry(2)
|
||||
{
|
||||
this->BuildShifts();
|
||||
};
|
||||
virtual ~NextToNearestStencilGeometry5D() {};
|
||||
virtual void BuildShifts(void)
|
||||
{
|
||||
this->shifts.resize(0);
|
||||
// Like HDCG: 81 point stencil including self connection
|
||||
this->shifts.push_back(Coordinate({0,0,0,0,0}));
|
||||
// +-x, +-y, +-z, +-t : 8
|
||||
for(int s=-1;s<=1;s+=2){
|
||||
this->shifts.push_back(Coordinate({0,s,0,0,0}));
|
||||
this->shifts.push_back(Coordinate({0,0,s,0,0}));
|
||||
this->shifts.push_back(Coordinate({0,0,0,s,0}));
|
||||
this->shifts.push_back(Coordinate({0,0,0,0,s}));
|
||||
}
|
||||
// +-x+-y, +-x+-z, +-x+-t, +-y+-z, +-y+-t, +-z+-t : 24
|
||||
for(int s1=-1;s1<=1;s1+=2){
|
||||
for(int s2=-1;s2<=1;s2+=2){
|
||||
this->shifts.push_back(Coordinate({0,s1,s2,0,0}));
|
||||
this->shifts.push_back(Coordinate({0,s1,0,s2,0}));
|
||||
this->shifts.push_back(Coordinate({0,s1,0,0,s2}));
|
||||
this->shifts.push_back(Coordinate({0,0,s1,s2,0}));
|
||||
this->shifts.push_back(Coordinate({0,0,s1,0,s2}));
|
||||
this->shifts.push_back(Coordinate({0,0,0,s1,s2}));
|
||||
}}
|
||||
this->npoint = this->shifts.size();
|
||||
}
|
||||
};
|
||||
// Need to worry about red-black now
|
||||
class NextToNextToNextToNearestStencilGeometry5D : public NonLocalStencilGeometry {
|
||||
public:
|
||||
NextToNextToNextToNearestStencilGeometry5D(void) : NonLocalStencilGeometry(4)
|
||||
{
|
||||
this->BuildShifts();
|
||||
};
|
||||
virtual ~NextToNextToNextToNearestStencilGeometry5D() {}
|
||||
virtual void BuildShifts(void)
|
||||
{
|
||||
this->shifts.resize(0);
|
||||
// Like HDCG: 81 point stencil including self connection
|
||||
this->shifts.push_back(Coordinate({0,0,0,0,0}));
|
||||
// +-x, +-y, +-z, +-t : 8
|
||||
for(int s=-1;s<=1;s+=2){
|
||||
this->shifts.push_back(Coordinate({0,s,0,0,0}));
|
||||
this->shifts.push_back(Coordinate({0,0,s,0,0}));
|
||||
this->shifts.push_back(Coordinate({0,0,0,s,0}));
|
||||
this->shifts.push_back(Coordinate({0,0,0,0,s}));
|
||||
}
|
||||
// +-x+-y, +-x+-z, +-x+-t, +-y+-z, +-y+-t, +-z+-t : 24
|
||||
for(int s1=-1;s1<=1;s1+=2){
|
||||
for(int s2=-1;s2<=1;s2+=2){
|
||||
this->shifts.push_back(Coordinate({0,s1,s2,0,0}));
|
||||
this->shifts.push_back(Coordinate({0,s1,0,s2,0}));
|
||||
this->shifts.push_back(Coordinate({0,s1,0,0,s2}));
|
||||
this->shifts.push_back(Coordinate({0,0,s1,s2,0}));
|
||||
this->shifts.push_back(Coordinate({0,0,s1,0,s2}));
|
||||
this->shifts.push_back(Coordinate({0,0,0,s1,s2}));
|
||||
}}
|
||||
// +-x+-y+-z, +-x+-y+-z, +-x+-y+-z,
|
||||
for(int s1=-1;s1<=1;s1+=2){
|
||||
for(int s2=-1;s2<=1;s2+=2){
|
||||
for(int s3=-1;s3<=1;s3+=2){
|
||||
this->shifts.push_back(Coordinate({0,s1,s2,s3,0})); // 8x4 = 32
|
||||
this->shifts.push_back(Coordinate({0,s1,s2,0,s3}));
|
||||
this->shifts.push_back(Coordinate({0,s1,0,s2,s3}));
|
||||
this->shifts.push_back(Coordinate({0,0,s1,s2,s3}));
|
||||
}}}
|
||||
for(int s1=-1;s1<=1;s1+=2){
|
||||
for(int s2=-1;s2<=1;s2+=2){
|
||||
for(int s3=-1;s3<=1;s3+=2){
|
||||
for(int s4=-1;s4<=1;s4+=2){
|
||||
this->shifts.push_back(Coordinate({0,s1,s2,s3,s4})); // 16
|
||||
}}}}
|
||||
this->npoint = this->shifts.size();
|
||||
}
|
||||
};
|
||||
|
||||
// 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<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;
|
||||
GridCartesian * _FineGrid;
|
||||
GridCartesian * _CoarseGrid;
|
||||
NonLocalStencilGeometry &geom;
|
||||
PaddedCell Cell;
|
||||
GeneralLocalStencil Stencil;
|
||||
|
||||
std::vector<CoarseMatrix> _A;
|
||||
std::vector<CoarseMatrix> _Adag;
|
||||
|
||||
///////////////////////
|
||||
// Interface
|
||||
///////////////////////
|
||||
GridCartesian * Grid(void) { return _FineGrid; }; // this is all the linalg routines need to know
|
||||
GridCartesian * 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,GridCartesian *FineGrid, GridCartesian * CoarseGrid)
|
||||
: geom(_geom),
|
||||
_FineGrid(FineGrid),
|
||||
_CoarseGrid(CoarseGrid),
|
||||
hermitian(1),
|
||||
Cell(_geom.Depth(),_CoarseGrid),
|
||||
Stencil(Cell.grids.back(),geom.shifts)
|
||||
{
|
||||
_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;
|
||||
std::cout << "Calling Exchange"<<std::endl;
|
||||
CoarseVector pin = Cell.Exchange(tin);
|
||||
// std::cout << "Called Exchange"<<std::endl;
|
||||
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);
|
||||
}
|
||||
|
||||
// Should also exchange "A" and "Adag"
|
||||
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;
|
||||
|
||||
// FIXME -- exchange the A and the A dag
|
||||
|
||||
for(int point=0;point<npoint;point++){
|
||||
|
||||
auto SE = Stencil_v.GetEntry(point,ss);
|
||||
|
||||
int o = SE->_offset;
|
||||
|
||||
// gpermute etc..
|
||||
nbr = in_v[o];
|
||||
assert( o< osites);
|
||||
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 Test(LinearOperatorBase<Lattice<Fobj> > &linop,
|
||||
Aggregation<Fobj,CComplex,nbasis> & Subspace)
|
||||
{
|
||||
// Create a random
|
||||
GridCartesian *grid = FineGrid();
|
||||
FineField MbV(grid);
|
||||
FineField tmp(grid);
|
||||
FineField f_src(grid);
|
||||
FineField f_res(grid);
|
||||
FineField f_ref(grid);
|
||||
CoarseVector c_src(CoarseGrid());
|
||||
CoarseVector c_res(CoarseGrid());
|
||||
CoarseVector coarseInner(CoarseGrid());
|
||||
GridParallelRNG RNG(CoarseGrid()); RNG.SeedUniqueString(std::string("Coarse RNG"));
|
||||
random(RNG,c_src);
|
||||
blockPromote(c_src,f_src,Subspace.subspace);
|
||||
linop.op(f_src,f_ref);
|
||||
this->Mult (_A,c_src,c_res);
|
||||
blockPromote(c_res,f_res,Subspace.subspace);
|
||||
std::cout << " GeneralCoarsenedMatrix comparison res "<<norm2(f_res)<<std::endl;
|
||||
std::cout << " GeneralCoarsenedMatrix comparison ref "<<norm2(f_ref)<<std::endl;
|
||||
f_res = f_res - f_ref;
|
||||
std::cout << " GeneralCoarsenedMatrix comparison diff "<<norm2(f_res)<<std::endl;
|
||||
}
|
||||
void CoarsenOperator(LinearOperatorBase<Lattice<Fobj> > &linop,
|
||||
Aggregation<Fobj,CComplex,nbasis> & Subspace)
|
||||
{
|
||||
std::cout << GridLogMessage<< "CoarsenMatrix "<< std::endl;
|
||||
GridCartesian *grid = FineGrid();
|
||||
// Orthogonalise the subblocks over the basis
|
||||
CoarseScalar InnerProd(CoarseGrid());
|
||||
for(int b=0;b<nbasis;b++){
|
||||
std::cout << "subspace["<<b<<"] " <<norm2(Subspace.subspace[b])<<std::endl;
|
||||
}
|
||||
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);
|
||||
std::cout << GridLogMessage<< "CoarsenMatrix block "<< bcoor << std::endl;
|
||||
for(int b=0;b<nbasis;b++){
|
||||
blockPick(CoarseGrid(),Subspace.subspace[b],bV,bcoor);
|
||||
linop.HermOp(bV,MbV);
|
||||
blockProject(coarseInner,MbV,Subspace.subspace);
|
||||
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
|
||||
}
|
||||
auto ip = peekSite(coarseInner,scoor);
|
||||
std::cout << "A["<<b<<"]["<<p<<"]"<<scoor<<" "<<" shift "<<geom.shifts[p]<<" "<< ip <<std::endl;
|
||||
auto Ab = peekSite(_A[p],scoor);
|
||||
auto Adagb = peekSite(_Adag[p],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[p],bcoor);
|
||||
}
|
||||
}
|
||||
}
|
||||
std::cout << " Exchanging _A " <<std::endl;
|
||||
for(int p=0;p<geom.npoint;p++){
|
||||
_A[p] = Cell.Exchange(_A[p]);
|
||||
_Adag[p] = Cell.Exchange(_Adag[p]);
|
||||
}
|
||||
}
|
||||
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.
|
||||
|
||||
@ -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;
|
||||
}
|
||||
|
||||
@ -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,329 +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=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);
|
||||
SU<Nc>::ColdConfiguration(Umu);
|
||||
// auto U = peekLorentz(Umu,0);
|
||||
// Umu=Zero(); // Make operator local for now
|
||||
// pokeLorentz(Umu,U,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 = 20;
|
||||
const int cb = 0 ;
|
||||
LatticeFermion prom(FGrid);
|
||||
|
||||
std::vector<LatticeFermion> subspace(nbasis,FGrid);
|
||||
|
||||
std::cout<<GridLogMessage<<"Calling Aggregation class" <<std::endl;
|
||||
|
||||
// Possible tactics -- with zero gauge field, verify block locality of dirac op
|
||||
// Possible tactics -- with zero gauge field, take inner products
|
||||
|
||||
// Squared operator
|
||||
MdagMLinearOperator<DomainWallFermionD,LatticeFermion> HermDefOp(Ddwf);
|
||||
DumbOperator<LatticeFermion> Diagonal(FGrid);
|
||||
typedef Aggregation<vSpinColourVector,vTComplex,nbasis> Subspace;
|
||||
|
||||
Subspace Aggregates(Coarse5d,FGrid,cb);
|
||||
Aggregates.CreateSubspaceRandom(RNG5);
|
||||
|
||||
std::cout<<GridLogMessage << "Called aggregation class"<< std::endl;
|
||||
|
||||
typedef GeneralCoarsenedMatrix<vSpinColourVector,vTComplex,nbasis> LittleDiracOperator;
|
||||
typedef LittleDiracOperator::CoarseVector CoarseVector;
|
||||
|
||||
NextToNearestStencilGeometry5D geom;
|
||||
LittleDiracOperator LittleDiracOp(geom,FGrid,Coarse5d);
|
||||
LittleDiracOp.CoarsenOperator(HermDefOp,Aggregates);
|
||||
// LittleDiracOp.CoarsenOperator(Diagonal,Aggregates);
|
||||
|
||||
std::cout<<GridLogMessage<<"Coarsened operator "<<std::endl;
|
||||
|
||||
CoarseVector c_src (Coarse5d);
|
||||
CoarseVector c_res (Coarse5d);
|
||||
CoarseVector c_proj(Coarse5d);
|
||||
|
||||
subspace=Aggregates.subspace;
|
||||
|
||||
Complex one(1.0);
|
||||
c_src = one; // 1 in every element for vector 1.
|
||||
Coordinate coor(5,0);
|
||||
|
||||
std::cout << "c_src"<< c_src<< std::endl;
|
||||
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;
|
||||
|
||||
// blockPick(Coarse5d,c_src,c_src,coor);
|
||||
// blockPromote(c_src,prom,subspace);
|
||||
|
||||
// Diagonal.HermOp(prom,tmp);
|
||||
HermDefOp.HermOp(prom,tmp);
|
||||
// HermDefOp.Op(prom,tmp);
|
||||
blockProject(c_proj,tmp,subspace);
|
||||
std::cout<<GridLogMessage<<" Called Big Dirac Op "<<norm2(tmp)<<std::endl;
|
||||
|
||||
LittleDiracOp.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 << "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)<<" Adag 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;
|
||||
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
|
||||
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 << "Done "<< std::endl;
|
||||
|
||||
Grid_finalize();
|
||||
return 0;
|
||||
}
|
||||
@ -1,226 +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=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
|
||||
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);
|
||||
|
||||
FieldMetaData header;
|
||||
std::string file("ckpoint_lat.4000");
|
||||
NerscIO::readConfiguration(Umu,header,file);
|
||||
|
||||
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 = 20;
|
||||
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 << "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