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

Debugged the threaded version. Cleaning up

Browse Source
This commit is contained in:
Guido Cossu 2016-12-07 04:40:25 +00:00
parent b812d5e39c
commit 143c70e29f
8 changed files with 243 additions and 285 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

36
lib/cartesian/Cartesian_base.h
View File

@ -122,6 +122,12 @@ public:
} }
inline void InOutCoorToLocalCoor (std::vector<int> &ocoor, std::vector<int> &icoor, std::vector<int> &lcoor) {
lcoor.resize(_ndimension);
for (int d = 0; d < _ndimension; d++)
lcoor[d] = ocoor[d] + _rdimensions[d] * icoor[d];
}
////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////
// SIMD lane addressing // SIMD lane addressing
////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////
@ -129,6 +135,7 @@ public:
{ {
Lexicographic::CoorFromIndex(coor,lane,_simd_layout); Lexicographic::CoorFromIndex(coor,lane,_simd_layout);
} }
inline int PermuteDim(int dimension){ inline int PermuteDim(int dimension){
return _simd_layout[dimension]>1; return _simd_layout[dimension]>1;
} }
@ -175,20 +182,22 @@ public:
inline const std::vector<int> &VirtualLocalDimensions(void) { return _ldimensions;}; inline const std::vector<int> &VirtualLocalDimensions(void) { return _ldimensions;};
//////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////
// Print decomposition // Utility to print the full decomposition details
//////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////
void show_decomposition(){ void show_decomposition(){
std::cout << GridLogMessage << "Full Dimensions : " << _fdimensions << std::endl; std::cout << GridLogMessage << "Full Dimensions : " << _fdimensions << std::endl;
std::cout << GridLogMessage << "Global Dimensions : " << _gdimensions << std::endl; std::cout << GridLogMessage << "Global Dimensions : " << _gdimensions << std::endl;
std::cout << GridLogMessage << "Local Dimensions : " << _ldimensions << std::endl; std::cout << GridLogMessage << "Local Dimensions : " << _ldimensions << std::endl;
std::cout << GridLogMessage << "Reduced Dimensions : " << _rdimensions << std::endl; std::cout << GridLogMessage << "Reduced Dimensions : " << _rdimensions << std::endl;
std::cout << GridLogMessage << "iSites : " << _isites << std::endl; std::cout << GridLogMessage << "Outer strides : " << _ostride << std::endl;
std::cout << GridLogMessage << "oSites : " << _osites << std::endl; std::cout << GridLogMessage << "Inner strides : " << _istride << std::endl;
std::cout << GridLogMessage << "lSites : " << lSites() << std::endl; std::cout << GridLogMessage << "iSites : " << _isites << std::endl;
std::cout << GridLogMessage << "gSites : " << gSites() << std::endl; std::cout << GridLogMessage << "oSites : " << _osites << std::endl;
std::cout << GridLogMessage << "Nd : " << _ndimension << std::endl; std::cout << GridLogMessage << "lSites : " << lSites() << std::endl;
} std::cout << GridLogMessage << "gSites : " << gSites() << std::endl;
std::cout << GridLogMessage << "Nd : " << _ndimension << std::endl;
}
//////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////
// Global addressing // Global addressing
@ -199,6 +208,9 @@ public:
void LocalIndexToLocalCoor(int lidx,std::vector<int> &lcoor){ void LocalIndexToLocalCoor(int lidx,std::vector<int> &lcoor){
Lexicographic::CoorFromIndex(lcoor,lidx,_ldimensions); Lexicographic::CoorFromIndex(lcoor,lidx,_ldimensions);
} }
void GlobalCoorToGlobalIndex(const std::vector<int> & gcoor,int & gidx){ void GlobalCoorToGlobalIndex(const std::vector<int> & gcoor,int & gidx){
gidx=0; gidx=0;
int mult=1; int mult=1;

52
lib/qcd/action/fermion/FermionOperatorImpl.h
View File

@ -292,35 +292,43 @@ class DomainWallVec5dImpl : public PeriodicGaugeImpl< GaugeImplTypes< S,Nrepres
} }
inline void InsertForce5D(GaugeField &mat, FermionField &Btilde, FermionField &Atilde, int mu) { inline void InsertForce5D(GaugeField &mat, FermionField &Btilde, FermionField &Atilde, int mu) {
int LLs = Btilde._grid->_rdimensions[0];
conformable(Atilde._grid,Btilde._grid);
GaugeLinkField tmp(mat._grid);
tmp = zero;
typedef decltype(traceIndex<SpinIndex>(outerProduct(Btilde[0], Atilde[0]))) result_type; typedef decltype(traceIndex<SpinIndex>(outerProduct(Btilde[0], Atilde[0]))) result_type;
std::vector<typename result_type::scalar_object> v_scalar_object(Btilde._grid->Nsimd()); unsigned int LLs = Btilde._grid->_rdimensions[0];
conformable(Atilde._grid,Btilde._grid);
GridBase* grid = mat._grid;
GridBase* Bgrid = Btilde._grid;
unsigned int dimU = grid->Nd();
unsigned int dimF = Bgrid->Nd();
GaugeLinkField tmp(grid);
tmp = zero;
// FIXME
// Current implementation works, thread safe, probably suboptimal
PARALLEL_FOR_LOOP PARALLEL_FOR_LOOP
for (int sss = 0; sss < tmp._grid->oSites(); sss++) { for (int so = 0; so < grid->oSites(); so++) {
std::vector<int> ocoor; std::vector<typename result_type::scalar_object> vres(Bgrid->Nsimd());
tmp._grid->oCoorFromOindex(ocoor,sss); std::vector<int> ocoor; grid->oCoorFromOindex(ocoor,so);
for (int si = 0; si < tmp._grid->iSites(); si++){ for (int si = 0; si < tmp._grid->iSites(); si++){
typename result_type::scalar_object scalar_object; typename result_type::scalar_object scalar_object; scalar_object = zero;
scalar_object = zero; std::vector<int> local_coor;
std::vector<int> local_coor(tmp._grid->Nd()); std::vector<int> icoor; grid->iCoorFromIindex(icoor,si);
std::vector<int> icoor; grid->InOutCoorToLocalCoor(ocoor, icoor, local_coor);
tmp._grid->iCoorFromIindex(icoor,si); //for (int i = 0; i < dimU; i++) local_coor[i] = ocoor[i] + grid->_rdimensions[i]*icoor[i];
for (int i = 0; i < tmp._grid->Nd(); i++) local_coor[i] = ocoor[i] + tmp._grid->_rdimensions[i]*icoor[i]; //std::cout << "so: " << so << " si: "<< si << " local_coor: " << local_coor << std::endl;
for (int s = 0; s < LLs; s++) { for (int s = 0; s < LLs; s++) {
std::vector<int> slocal_coor(Btilde._grid->Nd()); std::vector<int> slocal_coor(dimF);
slocal_coor[0] = s; slocal_coor[0] = s;
for (int s4d = 1; s4d< Btilde._grid->Nd(); s4d++) slocal_coor[s4d] = local_coor[s4d-1]; for (int s4d = 1; s4d< dimF; s4d++) slocal_coor[s4d] = local_coor[s4d-1];
int sF = Btilde._grid->oIndexReduced(slocal_coor); int sF = Bgrid->oIndexReduced(slocal_coor);
assert(sF < Btilde._grid->oSites()); assert(sF < Bgrid->oSites());
extract(traceIndex<SpinIndex>(outerProduct(Btilde[sF], Atilde[sF])), v_scalar_object);
for (int sv = 0; sv < v_scalar_object.size(); sv++) scalar_object += v_scalar_object[sv]; // sum across the 5d dimension extract(traceIndex<SpinIndex>(outerProduct(Btilde[sF], Atilde[sF])), vres);
// sum across the 5d dimension
for (auto v : vres) scalar_object += v;
} }
tmp._odata[sss].putlane(scalar_object, si); tmp._odata[so].putlane(scalar_object, si);
} }
} }
PokeIndex<LorentzIndex>(mat, tmp, mu); PokeIndex<LorentzIndex>(mat, tmp, mu);

12
lib/qcd/action/fermion/WilsonFermion.cc
View File

@ -222,8 +222,7 @@ void WilsonFermion<Impl>::DerivInternal(StencilImpl &st, DoubledGaugeField &U,
//////////////////////// ////////////////////////
PARALLEL_FOR_LOOP PARALLEL_FOR_LOOP
for (int sss = 0; sss < B._grid->oSites(); sss++) { for (int sss = 0; sss < B._grid->oSites(); sss++) {
Kernels::DiracOptDhopDir(st, U, st.CommBuf(), sss, sss, B, Btilde, mu, Kernels::DiracOptDhopDir(st, U, st.CommBuf(), sss, sss, B, Btilde, mu, gamma);
gamma);
} }
////////////////////////////////////////////////// //////////////////////////////////////////////////
@ -234,8 +233,7 @@ void WilsonFermion<Impl>::DerivInternal(StencilImpl &st, DoubledGaugeField &U,
} }
template <class Impl> template <class Impl>
void WilsonFermion<Impl>::DhopDeriv(GaugeField &mat, const FermionField &U, void WilsonFermion<Impl>::DhopDeriv(GaugeField &mat, const FermionField &U, const FermionField &V, int dag) {
const FermionField &V, int dag) {
conformable(U._grid, _grid); conformable(U._grid, _grid);
conformable(U._grid, V._grid); conformable(U._grid, V._grid);
conformable(U._grid, mat._grid); conformable(U._grid, mat._grid);
@ -246,8 +244,7 @@ void WilsonFermion<Impl>::DhopDeriv(GaugeField &mat, const FermionField &U,
} }
template <class Impl> template <class Impl>
void WilsonFermion<Impl>::DhopDerivOE(GaugeField &mat, const FermionField &U, void WilsonFermion<Impl>::DhopDerivOE(GaugeField &mat, const FermionField &U, const FermionField &V, int dag) {
const FermionField &V, int dag) {
conformable(U._grid, _cbgrid); conformable(U._grid, _cbgrid);
conformable(U._grid, V._grid); conformable(U._grid, V._grid);
conformable(U._grid, mat._grid); conformable(U._grid, mat._grid);
@ -260,8 +257,7 @@ void WilsonFermion<Impl>::DhopDerivOE(GaugeField &mat, const FermionField &U,
} }
template <class Impl> template <class Impl>
void WilsonFermion<Impl>::DhopDerivEO(GaugeField &mat, const FermionField &U, void WilsonFermion<Impl>::DhopDerivEO(GaugeField &mat, const FermionField &U, const FermionField &V, int dag) {
const FermionField &V, int dag) {
conformable(U._grid, _cbgrid); conformable(U._grid, _cbgrid);
conformable(U._grid, V._grid); conformable(U._grid, V._grid);
conformable(U._grid, mat._grid); conformable(U._grid, mat._grid);

70
lib/qcd/action/fermion/WilsonFermion5D.cc
View File

@ -11,6 +11,7 @@ Author: Peter Boyle <pabobyle@ph.ed.ac.uk>
Author: Peter Boyle <paboyle@ph.ed.ac.uk> Author: Peter Boyle <paboyle@ph.ed.ac.uk>
Author: Peter Boyle <peterboyle@Peters-MacBook-Pro-2.local> Author: Peter Boyle <peterboyle@Peters-MacBook-Pro-2.local>
Author: paboyle <paboyle@ph.ed.ac.uk> Author: paboyle <paboyle@ph.ed.ac.uk>
Author: Guido Cossu <guido.cossu@ed.ac.uk>
This program is free software; you can redistribute it and/or modify This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by it under the terms of the GNU General Public License as published by
@ -132,48 +133,6 @@ WilsonFermion5D<Impl>::WilsonFermion5D(GaugeField &_Umu,
// Allocate the required comms buffer // Allocate the required comms buffer
ImportGauge(_Umu); ImportGauge(_Umu);
} }
/*
template<class Impl>
WilsonFermion5D<Impl>::WilsonFermion5D(int simd,GaugeField &_Umu,
GridCartesian &FiveDimGrid,
GridRedBlackCartesian &FiveDimRedBlackGrid,
GridCartesian &FourDimGrid,
RealD _M5,const ImplParams &p) :
{
int nsimd = Simd::Nsimd();
// some assertions
assert(FiveDimGrid._ndimension==5);
assert(FiveDimRedBlackGrid._ndimension==5);
assert(FiveDimRedBlackGrid._checker_dim==0); // Checkerboard the s-direction
assert(FourDimGrid._ndimension==4);
// Dimension zero of the five-d is the Ls direction
Ls=FiveDimGrid._fdimensions[0];
assert(FiveDimGrid._processors[0] ==1);
assert(FiveDimGrid._simd_layout[0] ==nsimd);
assert(FiveDimRedBlackGrid._fdimensions[0]==Ls);
assert(FiveDimRedBlackGrid._processors[0] ==1);
assert(FiveDimRedBlackGrid._simd_layout[0]==nsimd);
// Other dimensions must match the decomposition of the four-D fields
for(int d=0;d<4;d++){
assert(FiveDimRedBlackGrid._fdimensions[d+1]==FourDimGrid._fdimensions[d]);
assert(FiveDimRedBlackGrid._processors[d+1] ==FourDimGrid._processors[d]);
assert(FourDimGrid._simd_layout[d]=1);
assert(FiveDimRedBlackGrid._simd_layout[d+1]==1);
assert(FiveDimGrid._fdimensions[d+1] ==FourDimGrid._fdimensions[d]);
assert(FiveDimGrid._processors[d+1] ==FourDimGrid._processors[d]);
assert(FiveDimGrid._simd_layout[d+1] ==FourDimGrid._simd_layout[d]);
}
{
}
}
*/
template<class Impl> template<class Impl>
void WilsonFermion5D<Impl>::Report(void) void WilsonFermion5D<Impl>::Report(void)
@ -346,13 +305,14 @@ void WilsonFermion5D<Impl>::DerivInternal(StencilImpl & st,
template<class Impl> template<class Impl>
void WilsonFermion5D<Impl>::DhopDeriv(GaugeField &mat, void WilsonFermion5D<Impl>::DhopDeriv(GaugeField &mat,
const FermionField &A, const FermionField &A,
const FermionField &B, const FermionField &B,
int dag) int dag)
{ {
conformable(A._grid,FermionGrid()); conformable(A._grid,FermionGrid());
conformable(A._grid,B._grid); conformable(A._grid,B._grid);
//conformable(GaugeGrid(),mat._grid);
//conformable(GaugeGrid(),mat._grid);// this is not general! leaving as a comment
mat.checkerboard = A.checkerboard; mat.checkerboard = A.checkerboard;
@ -361,12 +321,11 @@ void WilsonFermion5D<Impl>::DhopDeriv(GaugeField &mat,
template<class Impl> template<class Impl>
void WilsonFermion5D<Impl>::DhopDerivEO(GaugeField &mat, void WilsonFermion5D<Impl>::DhopDerivEO(GaugeField &mat,
const FermionField &A, const FermionField &A,
const FermionField &B, const FermionField &B,
int dag) int dag)
{ {
conformable(A._grid,FermionRedBlackGrid()); conformable(A._grid,FermionRedBlackGrid());
//conformable(GaugeRedBlackGrid(),mat._grid);
conformable(A._grid,B._grid); conformable(A._grid,B._grid);
assert(B.checkerboard==Odd); assert(B.checkerboard==Odd);
@ -379,12 +338,11 @@ void WilsonFermion5D<Impl>::DhopDerivEO(GaugeField &mat,
template<class Impl> template<class Impl>
void WilsonFermion5D<Impl>::DhopDerivOE(GaugeField &mat, void WilsonFermion5D<Impl>::DhopDerivOE(GaugeField &mat,
const FermionField &A, const FermionField &A,
const FermionField &B, const FermionField &B,
int dag) int dag)
{ {
conformable(A._grid,FermionRedBlackGrid()); conformable(A._grid,FermionRedBlackGrid());
//conformable(GaugeRedBlackGrid(),mat._grid);
conformable(A._grid,B._grid); conformable(A._grid,B._grid);
assert(B.checkerboard==Even); assert(B.checkerboard==Even);
@ -396,8 +354,8 @@ void WilsonFermion5D<Impl>::DhopDerivOE(GaugeField &mat,
template<class Impl> template<class Impl>
void WilsonFermion5D<Impl>::DhopInternal(StencilImpl & st, LebesgueOrder &lo, void WilsonFermion5D<Impl>::DhopInternal(StencilImpl & st, LebesgueOrder &lo,
DoubledGaugeField & U, DoubledGaugeField & U,
const FermionField &in, FermionField &out,int dag) const FermionField &in, FermionField &out,int dag)
{ {
// assert((dag==DaggerNo) ||(dag==DaggerYes)); // assert((dag==DaggerNo) ||(dag==DaggerYes));
Compressor compressor(dag); Compressor compressor(dag);

144
lib/qcd/action/pseudofermion/EvenOddSchurDifferentiable.h
View File

@ -7,6 +7,7 @@
Copyright (C) 2015 Copyright (C) 2015
Author: Peter Boyle <paboyle@ph.ed.ac.uk> Author: Peter Boyle <paboyle@ph.ed.ac.uk>
Author: Guido Cossu <guido.cossu@ed.ac.uk>
This program is free software; you can redistribute it and/or modify This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by it under the terms of the GNU General Public License as published by
@ -45,106 +46,97 @@ namespace Grid{
public: public:
INHERIT_IMPL_TYPES(Impl); INHERIT_IMPL_TYPES(Impl);
typedef FermionOperator<Impl> Matrix; typedef FermionOperator<Impl> Matrix;
SchurDifferentiableOperator (Matrix &Mat) : SchurDiagMooeeOperator<Matrix,FermionField>(Mat) {}; SchurDifferentiableOperator (Matrix &Mat) : SchurDiagMooeeOperator<Matrix,FermionField>(Mat) {};
void MpcDeriv(GaugeField &Force,const FermionField &U,const FermionField &V) { void MpcDeriv(GaugeField &Force,const FermionField &U,const FermionField &V) {
GridBase *fgrid = this->_Mat.FermionGrid(); GridBase *fgrid = this->_Mat.FermionGrid();
GridBase *fcbgrid = this->_Mat.FermionRedBlackGrid(); GridBase *fcbgrid = this->_Mat.FermionRedBlackGrid();
GridBase *ugrid = this->_Mat.GaugeGrid();
GridBase *ucbgrid = this->_Mat.GaugeRedBlackGrid();
Real coeff = 1.0; FermionField tmp1(fcbgrid);
FermionField tmp2(fcbgrid);
FermionField tmp1(fcbgrid); conformable(fcbgrid,U._grid);
FermionField tmp2(fcbgrid); conformable(fcbgrid,V._grid);
conformable(fcbgrid,U._grid); // Assert the checkerboard?? or code for either
conformable(fcbgrid,V._grid); assert(U.checkerboard==Odd);
assert(V.checkerboard==U.checkerboard);
// Assert the checkerboard?? or code for either // NOTE Guido: WE DO NOT WANT TO USE THE ucbgrid GRID FOR THE FORCE
assert(U.checkerboard==Odd); // it is not conformable with the HMC force field
assert(V.checkerboard==U.checkerboard); // INHERIT FROM THE Force field instead
// NOTE Guido: WE DO NOT WANT TO USE THIS GRID FOR THE FORCE
// INHERIT FROM THE Force field
//GaugeField ForceO(ucbgrid);
//GaugeField ForceE(ucbgrid);
GridRedBlackCartesian* forcecb = new GridRedBlackCartesian(Force._grid); GridRedBlackCartesian* forcecb = new GridRedBlackCartesian(Force._grid);
GaugeField ForceO(forcecb); GaugeField ForceO(forcecb);
GaugeField ForceE(forcecb); GaugeField ForceE(forcecb);
// X^dag Der_oe MeeInv Meo Y // X^dag Der_oe MeeInv Meo Y
// Use Mooee as nontrivial but gauge field indept // Use Mooee as nontrivial but gauge field indept
this->_Mat.Meooe (V,tmp1); // odd->even -- implicit -0.5 factor to be applied this->_Mat.Meooe (V,tmp1); // odd->even -- implicit -0.5 factor to be applied
this->_Mat.MooeeInv(tmp1,tmp2); // even->even this->_Mat.MooeeInv(tmp1,tmp2); // even->even
this->_Mat.MoeDeriv(ForceO,U,tmp2,DaggerNo); this->_Mat.MoeDeriv(ForceO,U,tmp2,DaggerNo);
// Accumulate X^dag M_oe MeeInv Der_eo Y // Accumulate X^dag M_oe MeeInv Der_eo Y
this->_Mat.MeooeDag (U,tmp1); // even->odd -- implicit -0.5 factor to be applied this->_Mat.MeooeDag (U,tmp1); // even->odd -- implicit -0.5 factor to be applied
this->_Mat.MooeeInvDag(tmp1,tmp2); // even->even this->_Mat.MooeeInvDag(tmp1,tmp2); // even->even
this->_Mat.MeoDeriv(ForceE,tmp2,V,DaggerNo); this->_Mat.MeoDeriv(ForceE,tmp2,V,DaggerNo);
assert(ForceE.checkerboard==Even); assert(ForceE.checkerboard==Even);
assert(ForceO.checkerboard==Odd); assert(ForceO.checkerboard==Odd);
setCheckerboard(Force,ForceE); setCheckerboard(Force,ForceE);
setCheckerboard(Force,ForceO); setCheckerboard(Force,ForceO);
Force=-Force; Force=-Force;
}
delete forcecb;
}
void MpcDagDeriv(GaugeField &Force,const FermionField &U,const FermionField &V) { void MpcDagDeriv(GaugeField &Force,const FermionField &U,const FermionField &V) {
GridBase *fgrid = this->_Mat.FermionGrid(); GridBase *fgrid = this->_Mat.FermionGrid();
GridBase *fcbgrid = this->_Mat.FermionRedBlackGrid(); GridBase *fcbgrid = this->_Mat.FermionRedBlackGrid();
GridBase *ugrid = this->_Mat.GaugeGrid();
GridBase *ucbgrid = this->_Mat.GaugeRedBlackGrid();
Real coeff = 1.0; FermionField tmp1(fcbgrid);
FermionField tmp2(fcbgrid);
FermionField tmp1(fcbgrid); conformable(fcbgrid,U._grid);
FermionField tmp2(fcbgrid); conformable(fcbgrid,V._grid);
conformable(fcbgrid,U._grid); // Assert the checkerboard?? or code for either
conformable(fcbgrid,V._grid); assert(V.checkerboard==Odd);
assert(V.checkerboard==V.checkerboard);
// Assert the checkerboard?? or code for either // NOTE Guido: WE DO NOT WANT TO USE THE ucbgrid GRID FOR THE FORCE
assert(V.checkerboard==Odd); // it is not conformable with the HMC force field
assert(V.checkerboard==V.checkerboard); // INHERIT FROM THE Force field instead
GridRedBlackCartesian* forcecb = new GridRedBlackCartesian(Force._grid);
// NOTE Guido: WE DO NOT WANT TO USE THIS GRID FOR THE FORCE
// INHERIT FROM THE Force field
//GaugeField ForceO(ucbgrid);
//GaugeField ForceE(ucbgrid);
GridRedBlackCartesian* forcecb = new GridRedBlackCartesian(Force._grid);
GaugeField ForceO(forcecb); GaugeField ForceO(forcecb);
GaugeField ForceE(forcecb); GaugeField ForceE(forcecb);
// X^dag Der_oe MeeInv Meo Y // X^dag Der_oe MeeInv Meo Y
// Use Mooee as nontrivial but gauge field indept // Use Mooee as nontrivial but gauge field indept
this->_Mat.MeooeDag (V,tmp1); // odd->even -- implicit -0.5 factor to be applied this->_Mat.MeooeDag (V,tmp1); // odd->even -- implicit -0.5 factor to be applied
this->_Mat.MooeeInvDag(tmp1,tmp2); // even->even this->_Mat.MooeeInvDag(tmp1,tmp2); // even->even
this->_Mat.MoeDeriv(ForceO,U,tmp2,DaggerYes); this->_Mat.MoeDeriv(ForceO,U,tmp2,DaggerYes);
// Accumulate X^dag M_oe MeeInv Der_eo Y // Accumulate X^dag M_oe MeeInv Der_eo Y
this->_Mat.Meooe (U,tmp1); // even->odd -- implicit -0.5 factor to be applied this->_Mat.Meooe (U,tmp1); // even->odd -- implicit -0.5 factor to be applied
this->_Mat.MooeeInv(tmp1,tmp2); // even->even this->_Mat.MooeeInv(tmp1,tmp2); // even->even
this->_Mat.MeoDeriv(ForceE,tmp2,V,DaggerYes); this->_Mat.MeoDeriv(ForceE,tmp2,V,DaggerYes);
assert(ForceE.checkerboard==Even); assert(ForceE.checkerboard==Even);
assert(ForceO.checkerboard==Odd); assert(ForceO.checkerboard==Odd);
setCheckerboard(Force,ForceE); setCheckerboard(Force,ForceE);
setCheckerboard(Force,ForceO); setCheckerboard(Force,ForceO);
Force=-Force; Force=-Force;
delete forcecb;
} }
}; };

188
lib/qcd/action/pseudofermion/TwoFlavourEvenOddRatio.h
View File

@ -52,68 +52,68 @@ namespace Grid{
public: public:
TwoFlavourEvenOddRatioPseudoFermionAction(FermionOperator<Impl> &_NumOp, TwoFlavourEvenOddRatioPseudoFermionAction(FermionOperator<Impl> &_NumOp,
FermionOperator<Impl> &_DenOp, FermionOperator<Impl> &_DenOp,
OperatorFunction<FermionField> & DS, OperatorFunction<FermionField> & DS,
OperatorFunction<FermionField> & AS) : OperatorFunction<FermionField> & AS) :
NumOp(_NumOp), NumOp(_NumOp),
DenOp(_DenOp), DenOp(_DenOp),
DerivativeSolver(DS), DerivativeSolver(DS),
ActionSolver(AS), ActionSolver(AS),
PhiEven(_NumOp.FermionRedBlackGrid()), PhiEven(_NumOp.FermionRedBlackGrid()),
PhiOdd(_NumOp.FermionRedBlackGrid()) PhiOdd(_NumOp.FermionRedBlackGrid())
{ {
conformable(_NumOp.FermionGrid(), _DenOp.FermionGrid()); conformable(_NumOp.FermionGrid(), _DenOp.FermionGrid());
conformable(_NumOp.FermionRedBlackGrid(), _DenOp.FermionRedBlackGrid()); conformable(_NumOp.FermionRedBlackGrid(), _DenOp.FermionRedBlackGrid());
conformable(_NumOp.GaugeGrid(), _DenOp.GaugeGrid()); conformable(_NumOp.GaugeGrid(), _DenOp.GaugeGrid());
conformable(_NumOp.GaugeRedBlackGrid(), _DenOp.GaugeRedBlackGrid()); conformable(_NumOp.GaugeRedBlackGrid(), _DenOp.GaugeRedBlackGrid());
}; };
virtual std::string action_name(){return "TwoFlavourEvenOddRatioPseudoFermionAction";} virtual std::string action_name(){return "TwoFlavourEvenOddRatioPseudoFermionAction";}
virtual void refresh(const GaugeField &U, GridParallelRNG& pRNG) { virtual void refresh(const GaugeField &U, GridParallelRNG& pRNG) {
// P(phi) = e^{- phi^dag Vpc (MpcdagMpc)^-1 Vpcdag phi} // P(phi) = e^{- phi^dag Vpc (MpcdagMpc)^-1 Vpcdag phi}
// //
// NumOp == V // NumOp == V
// DenOp == M // DenOp == M
// //
// Take phi_o = Vpcdag^{-1} Mpcdag eta_o ; eta_o = Mpcdag^{-1} Vpcdag Phi // Take phi_o = Vpcdag^{-1} Mpcdag eta_o ; eta_o = Mpcdag^{-1} Vpcdag Phi
// //
// P(eta_o) = e^{- eta_o^dag eta_o} // P(eta_o) = e^{- eta_o^dag eta_o}
// //
// e^{x^2/2 sig^2} => sig^2 = 0.5. // e^{x^2/2 sig^2} => sig^2 = 0.5.
// //
RealD scale = std::sqrt(0.5); RealD scale = std::sqrt(0.5);
FermionField eta (NumOp.FermionGrid()); FermionField eta (NumOp.FermionGrid());
FermionField etaOdd (NumOp.FermionRedBlackGrid()); FermionField etaOdd (NumOp.FermionRedBlackGrid());
FermionField etaEven(NumOp.FermionRedBlackGrid()); FermionField etaEven(NumOp.FermionRedBlackGrid());
FermionField tmp (NumOp.FermionRedBlackGrid()); FermionField tmp (NumOp.FermionRedBlackGrid());
gaussian(pRNG,eta); gaussian(pRNG,eta);
pickCheckerboard(Even,etaEven,eta); pickCheckerboard(Even,etaEven,eta);
pickCheckerboard(Odd,etaOdd,eta); pickCheckerboard(Odd,etaOdd,eta);
NumOp.ImportGauge(U); NumOp.ImportGauge(U);
DenOp.ImportGauge(U); DenOp.ImportGauge(U);
SchurDifferentiableOperator<Impl> Mpc(DenOp); SchurDifferentiableOperator<Impl> Mpc(DenOp);
SchurDifferentiableOperator<Impl> Vpc(NumOp); SchurDifferentiableOperator<Impl> Vpc(NumOp);
// Odd det factors // Odd det factors
Mpc.MpcDag(etaOdd,PhiOdd); Mpc.MpcDag(etaOdd,PhiOdd);
tmp=zero; tmp=zero;
ActionSolver(Vpc,PhiOdd,tmp); ActionSolver(Vpc,PhiOdd,tmp);
Vpc.Mpc(tmp,PhiOdd); Vpc.Mpc(tmp,PhiOdd);
// Even det factors // Even det factors
DenOp.MooeeDag(etaEven,tmp); DenOp.MooeeDag(etaEven,tmp);
NumOp.MooeeInvDag(tmp,PhiEven); NumOp.MooeeInvDag(tmp,PhiEven);
PhiOdd =PhiOdd*scale; PhiOdd =PhiOdd*scale;
PhiEven=PhiEven*scale; PhiEven=PhiEven*scale;
}; };
////////////////////////////////////////////////////// //////////////////////////////////////////////////////
@ -121,33 +121,33 @@ namespace Grid{
////////////////////////////////////////////////////// //////////////////////////////////////////////////////
virtual RealD S(const GaugeField &U) { virtual RealD S(const GaugeField &U) {
NumOp.ImportGauge(U); NumOp.ImportGauge(U);
DenOp.ImportGauge(U); DenOp.ImportGauge(U);
SchurDifferentiableOperator<Impl> Mpc(DenOp); SchurDifferentiableOperator<Impl> Mpc(DenOp);
SchurDifferentiableOperator<Impl> Vpc(NumOp); SchurDifferentiableOperator<Impl> Vpc(NumOp);
FermionField X(NumOp.FermionRedBlackGrid()); FermionField X(NumOp.FermionRedBlackGrid());
FermionField Y(NumOp.FermionRedBlackGrid()); FermionField Y(NumOp.FermionRedBlackGrid());
Vpc.MpcDag(PhiOdd,Y); // Y= Vdag phi Vpc.MpcDag(PhiOdd,Y); // Y= Vdag phi
X=zero; X=zero;
ActionSolver(Mpc,Y,X); // X= (MdagM)^-1 Vdag phi ActionSolver(Mpc,Y,X); // X= (MdagM)^-1 Vdag phi
//Mpc.Mpc(X,Y); // Y= Mdag^-1 Vdag phi //Mpc.Mpc(X,Y); // Y= Mdag^-1 Vdag phi
// Multiply by Ydag // Multiply by Ydag
RealD action = real(innerProduct(Y,X)); RealD action = real(innerProduct(Y,X));
//RealD action = norm2(Y); //RealD action = norm2(Y);
// The EE factorised block; normally can replace with zero if det is constant (gauge field indept) // The EE factorised block; normally can replace with zero if det is constant (gauge field indept)
// Only really clover term that creates this. Leave the EE portion as a future to do to make most // Only really clover term that creates this. Leave the EE portion as a future to do to make most
// rapid progresss on DWF for now. // rapid progresss on DWF for now.
// //
NumOp.MooeeDag(PhiEven,X); NumOp.MooeeDag(PhiEven,X);
DenOp.MooeeInvDag(X,Y); DenOp.MooeeInvDag(X,Y);
action = action + norm2(Y); action = action + norm2(Y);
return action; return action;
}; };
////////////////////////////////////////////////////// //////////////////////////////////////////////////////
@ -157,46 +157,44 @@ namespace Grid{
////////////////////////////////////////////////////// //////////////////////////////////////////////////////
virtual void deriv(const GaugeField &U,GaugeField & dSdU) { virtual void deriv(const GaugeField &U,GaugeField & dSdU) {
NumOp.ImportGauge(U); NumOp.ImportGauge(U);
DenOp.ImportGauge(U); DenOp.ImportGauge(U);
SchurDifferentiableOperator<Impl> Mpc(DenOp); SchurDifferentiableOperator<Impl> Mpc(DenOp);
SchurDifferentiableOperator<Impl> Vpc(NumOp); SchurDifferentiableOperator<Impl> Vpc(NumOp);
FermionField X(NumOp.FermionRedBlackGrid()); FermionField X(NumOp.FermionRedBlackGrid());
FermionField Y(NumOp.FermionRedBlackGrid()); FermionField Y(NumOp.FermionRedBlackGrid());
//GaugeField force(NumOp.GaugeGrid()); // This assignment is necessary to be compliant with the HMC grids
GaugeField force(dSdU._grid); GaugeField force(dSdU._grid);
conformable(force._grid, dSdU._grid);
//Y=Vdag phi //Y=Vdag phi
//X = (Mdag M)^-1 V^dag phi //X = (Mdag M)^-1 V^dag phi
//Y = (Mdag)^-1 V^dag phi //Y = (Mdag)^-1 V^dag phi
Vpc.MpcDag(PhiOdd,Y); // Y= Vdag phi Vpc.MpcDag(PhiOdd,Y); // Y= Vdag phi
X=zero; X=zero;
DerivativeSolver(Mpc,Y,X); // X= (MdagM)^-1 Vdag phi DerivativeSolver(Mpc,Y,X); // X= (MdagM)^-1 Vdag phi
Mpc.Mpc(X,Y); // Y= Mdag^-1 Vdag phi Mpc.Mpc(X,Y); // Y= Mdag^-1 Vdag phi
// phi^dag V (Mdag M)^-1 dV^dag phi // phi^dag V (Mdag M)^-1 dV^dag phi
Vpc.MpcDagDeriv(force , X, PhiOdd ); dSdU=force; Vpc.MpcDagDeriv(force , X, PhiOdd ); dSdU = force;
// phi^dag dV (Mdag M)^-1 V^dag phi // phi^dag dV (Mdag M)^-1 V^dag phi
Vpc.MpcDeriv(force , PhiOdd, X ); dSdU=dSdU+force; Vpc.MpcDeriv(force , PhiOdd, X ); dSdU = dSdU+force;
// - phi^dag V (Mdag M)^-1 Mdag dM (Mdag M)^-1 V^dag phi // - phi^dag V (Mdag M)^-1 Mdag dM (Mdag M)^-1 V^dag phi
// - phi^dag V (Mdag M)^-1 dMdag M (Mdag M)^-1 V^dag phi // - phi^dag V (Mdag M)^-1 dMdag M (Mdag M)^-1 V^dag phi
Mpc.MpcDeriv(force,Y,X); dSdU=dSdU-force; Mpc.MpcDeriv(force,Y,X); dSdU = dSdU-force;
Mpc.MpcDagDeriv(force,X,Y); dSdU=dSdU-force; Mpc.MpcDagDeriv(force,X,Y); dSdU = dSdU-force;
// FIXME No force contribution from EvenEven assumed here // FIXME No force contribution from EvenEven assumed here
// Needs a fix for clover. // Needs a fix for clover.
assert(NumOp.ConstEE() == 1); assert(NumOp.ConstEE() == 1);
assert(DenOp.ConstEE() == 1); assert(DenOp.ConstEE() == 1);
//dSdU = -Ta(dSdU); dSdU = -dSdU;
dSdU = -dSdU;
}; };
}; };
} }

25
lib/qcd/hmc/integrators/Integrator.h
View File

@ -57,7 +57,7 @@ struct IntegratorParameters {
void print_parameters() { void print_parameters() {
std::cout << GridLogMessage << "[Integrator] Trajectory length : " << trajL << std::endl; std::cout << GridLogMessage << "[Integrator] Trajectory length : " << trajL << std::endl;
std::cout << GridLogMessage << "[Integrator] Number of MD steps : " << MDsteps << std::endl; std::cout << GridLogMessage << "[Integrator] Number of MD steps : " << MDsteps << std::endl;
std::cout << GridLogMessage << "[Integrator] Step size : " << stepsize << std::endl; std::cout << GridLogMessage << "[Integrator] Step size : " << stepsize << std::endl;
} }
@ -99,11 +99,9 @@ class Integrator {
FieldType forceR(U._grid); FieldType forceR(U._grid);
// Implement smearing only for the fundamental representation now // Implement smearing only for the fundamental representation now
repr_set.at(a)->deriv(Rep.U, forceR); repr_set.at(a)->deriv(Rep.U, forceR);
GF force = GF force = Rep.RtoFundamentalProject(forceR); // Ta for the fundamental rep
Rep.RtoFundamentalProject(forceR); // Ta for the fundamental rep Real force_abs = std::sqrt(norm2(force)/(U._grid->gSites()));
Real force_abs = std::sqrt(norm2(force)); std::cout << GridLogIntegrator << "Hirep Force average: " << force_abs << std::endl;
std::cout << GridLogIntegrator << "Hirep Force average: "
<< force_abs / (U._grid->gSites()) << std::endl;
Mom -= force * ep ; Mom -= force * ep ;
} }
} }
@ -118,15 +116,11 @@ class Integrator {
Field& Us = Smearer.get_U(as[level].actions.at(a)->is_smeared); Field& Us = Smearer.get_U(as[level].actions.at(a)->is_smeared);
as[level].actions.at(a)->deriv(Us, force); // deriv should NOT include Ta as[level].actions.at(a)->deriv(Us, force); // deriv should NOT include Ta
std::cout << GridLogIntegrator std::cout << GridLogIntegrator << "Smearing (on/off): " << as[level].actions.at(a)->is_smeared << std::endl;
<< "Smearing (on/off): " << as[level].actions.at(a)->is_smeared
<< std::endl;
if (as[level].actions.at(a)->is_smeared) Smearer.smeared_force(force); if (as[level].actions.at(a)->is_smeared) Smearer.smeared_force(force);
force = FieldImplementation::projectForce(force); // Ta for gauge fields force = FieldImplementation::projectForce(force); // Ta for gauge fields
Real force_abs = std::sqrt(norm2(force)); Real force_abs = std::sqrt(norm2(force)/U._grid->gSites());
std::cout << GridLogIntegrator std::cout << GridLogIntegrator << "Force average: " << force_abs << std::endl;
<< "Force average: " << force_abs / (U._grid->gSites())
<< std::endl;
Mom -= force * ep; Mom -= force * ep;
} }
@ -139,10 +133,9 @@ class Integrator {
t_U += ep; t_U += ep;
int fl = levels - 1; int fl = levels - 1;
std::cout << GridLogIntegrator << " " std::cout << GridLogIntegrator << " " << "[" << fl << "] U " << " dt " << ep << " : t_U " << t_U << std::endl;
<< "[" << fl << "] U "
<< " dt " << ep << " : t_U " << t_U << std::endl;
} }
void update_U(MomentaField& Mom, Field& U, double ep) { void update_U(MomentaField& Mom, Field& U, double ep) {
// exponential of Mom*U in the gauge fields case // exponential of Mom*U in the gauge fields case
FieldImplementation::update_field(Mom, U, ep); FieldImplementation::update_field(Mom, U, ep);

1
tests/hmc/Test_hmc_EODWFRatio_Binary.cc
View File

@ -85,6 +85,7 @@ public:
FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid); FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid);
*/ */
// temporarily need a gauge field // temporarily need a gauge field
LatticeGaugeField U(UGrid); LatticeGaugeField U(UGrid);