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

Merge branch 'feature/ddhmc' of https://github.com/paboyle/Grid into feature/ddhmc

Browse Source
This commit is contained in:
Peter Boyle 2021-05-06 23:51:25 +02:00
commit b5b930d5bb
33 changed files with 2093 additions and 93 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

19
Grid/algorithms/LinearOperator.h
View File

@ -223,9 +223,14 @@ class SchurOperatorBase : public LinearOperatorBase<Field> {
Mpc(in,tmp); Mpc(in,tmp);
MpcDag(tmp,out); MpcDag(tmp,out);
} }
virtual void MpcMpcDag(const Field &in, Field &out) {
Field tmp(in.Grid());
tmp.Checkerboard() = in.Checkerboard();
MpcDag(in,tmp);
Mpc(tmp,out);
}
virtual void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){ virtual void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){
out.Checkerboard() = in.Checkerboard(); HermOp(in,out);
MpcDagMpc(in,out);
ComplexD dot= innerProduct(in,out); ComplexD dot= innerProduct(in,out);
n1=real(dot); n1=real(dot);
n2=norm2(out); n2=norm2(out);
@ -276,6 +281,16 @@ template<class Matrix,class Field>
axpy(out,-1.0,tmp,out); axpy(out,-1.0,tmp,out);
} }
}; };
// Mpc MpcDag system presented as the HermOp
template<class Matrix,class Field>
class SchurDiagMooeeDagOperator : public SchurDiagMooeeOperator<Matrix,Field> {
public:
virtual void HermOp(const Field &in, Field &out){
out.Checkerboard() = in.Checkerboard();
this->MpcMpcDag(in,out);
}
SchurDiagMooeeDagOperator (Matrix &Mat): SchurDiagMooeeOperator<Matrix,Field>(Mat){};
};
template<class Matrix,class Field> template<class Matrix,class Field>
class SchurDiagOneOperator : public SchurOperatorBase<Field> { class SchurDiagOneOperator : public SchurOperatorBase<Field> {
protected: protected:

154
Grid/algorithms/iterative/SchurRedBlack.h
View File

@ -40,7 +40,7 @@ Author: Peter Boyle <paboyle@ph.ed.ac.uk>
* (-MoeMee^{-1} 1 ) * (-MoeMee^{-1} 1 )
* L^{dag} = ( 1 Mee^{-dag} Moe^{dag} ) * L^{dag} = ( 1 Mee^{-dag} Moe^{dag} )
* ( 0 1 ) * ( 0 1 )
* L^{-d} = ( 1 -Mee^{-dag} Moe^{dag} ) * L^{-dag}= ( 1 -Mee^{-dag} Moe^{dag} )
* ( 0 1 ) * ( 0 1 )
* *
* U^-1 = (1 -Mee^{-1} Meo) * U^-1 = (1 -Mee^{-1} Meo)
@ -82,7 +82,8 @@ Author: Peter Boyle <paboyle@ph.ed.ac.uk>
* c) M_oo^-dag Doo^{dag} Doo Moo^-1 phi_0 = M_oo^-dag (D_oo)^dag L^{-1} eta_o * c) M_oo^-dag Doo^{dag} Doo Moo^-1 phi_0 = M_oo^-dag (D_oo)^dag L^{-1} eta_o
* eta_o' = M_oo^-dag (D_oo)^dag (eta_o - Moe Mee^{-1} eta_e) * eta_o' = M_oo^-dag (D_oo)^dag (eta_o - Moe Mee^{-1} eta_e)
* psi_o = M_oo^-1 phi_o * psi_o = M_oo^-1 phi_o
* TODO: Deflation *
*
*/ */
namespace Grid { namespace Grid {
@ -97,6 +98,7 @@ namespace Grid {
protected: protected:
typedef CheckerBoardedSparseMatrixBase<Field> Matrix; typedef CheckerBoardedSparseMatrixBase<Field> Matrix;
OperatorFunction<Field> & _HermitianRBSolver; OperatorFunction<Field> & _HermitianRBSolver;
int CBfactorise; int CBfactorise;
bool subGuess; bool subGuess;
bool useSolnAsInitGuess; // if true user-supplied solution vector is used as initial guess for solver bool useSolnAsInitGuess; // if true user-supplied solution vector is used as initial guess for solver
@ -189,13 +191,20 @@ namespace Grid {
///////////////////////////////////////////////// /////////////////////////////////////////////////
// Check unprec residual if possible // Check unprec residual if possible
///////////////////////////////////////////////// /////////////////////////////////////////////////
if ( ! subGuess ) { if ( ! subGuess ) {
_Matrix.M(out[b],resid);
if ( this->adjoint() ) _Matrix.Mdag(out[b],resid);
else _Matrix.M(out[b],resid);
resid = resid-in[b]; resid = resid-in[b];
RealD ns = norm2(in[b]); RealD ns = norm2(in[b]);
RealD nr = norm2(resid); RealD nr = norm2(resid);
std::cout<<GridLogMessage<< "SchurRedBlackBase solver true unprec resid["<<b<<"] "<<std::sqrt(nr/ns) << std::endl; std::cout<<GridLogMessage<< "SchurRedBlackBase adjoint "<< this->adjoint() << std::endl;
if ( this->adjoint() )
std::cout<<GridLogMessage<< "SchurRedBlackBase adjoint solver true unprec resid["<<b<<"] "<<std::sqrt(nr/ns) << std::endl;
else
std::cout<<GridLogMessage<< "SchurRedBlackBase solver true unprec resid["<<b<<"] "<<std::sqrt(nr/ns) << std::endl;
} else { } else {
std::cout<<GridLogMessage<< "SchurRedBlackBase Guess subtracted after solve["<<b<<"] " << std::endl; std::cout<<GridLogMessage<< "SchurRedBlackBase Guess subtracted after solve["<<b<<"] " << std::endl;
} }
@ -249,12 +258,21 @@ namespace Grid {
// Verify the unprec residual // Verify the unprec residual
if ( ! subGuess ) { if ( ! subGuess ) {
_Matrix.M(out,resid);
std::cout<<GridLogMessage<< "SchurRedBlackBase adjoint "<< this->adjoint() << std::endl;
if ( this->adjoint() ) _Matrix.Mdag(out,resid);
else _Matrix.M(out,resid);
resid = resid-in; resid = resid-in;
RealD ns = norm2(in); RealD ns = norm2(in);
RealD nr = norm2(resid); RealD nr = norm2(resid);
std::cout<<GridLogMessage << "SchurRedBlackBase solver true unprec resid "<< std::sqrt(nr/ns) << std::endl; if ( this->adjoint() )
std::cout<<GridLogMessage<< "SchurRedBlackBase adjoint solver true unprec resid "<<std::sqrt(nr/ns) << std::endl;
else
std::cout<<GridLogMessage<< "SchurRedBlackBase solver true unprec resid "<<std::sqrt(nr/ns) << std::endl;
} else { } else {
std::cout << GridLogMessage << "SchurRedBlackBase Guess subtracted after solve." << std::endl; std::cout << GridLogMessage << "SchurRedBlackBase Guess subtracted after solve." << std::endl;
} }
@ -263,6 +281,7 @@ namespace Grid {
///////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////
// Override in derived. // Override in derived.
///////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////
virtual bool adjoint(void) { return false; }
virtual void RedBlackSource (Matrix & _Matrix,const Field &src, Field &src_e,Field &src_o) =0; virtual void RedBlackSource (Matrix & _Matrix,const Field &src, Field &src_e,Field &src_o) =0;
virtual void RedBlackSolution(Matrix & _Matrix,const Field &sol_o, const Field &src_e,Field &sol) =0; virtual void RedBlackSolution(Matrix & _Matrix,const Field &sol_o, const Field &src_e,Field &sol) =0;
virtual void RedBlackSolve (Matrix & _Matrix,const Field &src_o, Field &sol_o) =0; virtual void RedBlackSolve (Matrix & _Matrix,const Field &src_o, Field &sol_o) =0;
@ -616,6 +635,127 @@ namespace Grid {
this->_HermitianRBSolver(_OpEO, src_o, sol_o); this->_HermitianRBSolver(_OpEO, src_o, sol_o);
} }
}; };
/*
* Red black Schur decomposition
*
* M = (Mee Meo) = (1 0 ) (Mee 0 ) (1 Mee^{-1} Meo)
* (Moe Moo) (Moe Mee^-1 1 ) (0 Moo-Moe Mee^-1 Meo) (0 1 )
* = L D U
*
* L^-1 = (1 0 )
* (-MoeMee^{-1} 1 )
* L^{dag} = ( 1 Mee^{-dag} Moe^{dag} )
* ( 0 1 )
*
* U^-1 = (1 -Mee^{-1} Meo)
* (0 1 )
* U^{dag} = ( 1 0)
* (Meo^dag Mee^{-dag} 1)
* U^{-dag} = ( 1 0)
* (-Meo^dag Mee^{-dag} 1)
*
*
***********************
* M^dag psi = eta
***********************
*
* Really for Mobius: (Wilson - easier to just use gamma 5 hermiticity)
*
* Mdag psi = Udag Ddag Ldag psi = eta
*
* U^{-dag} = ( 1 0)
* (-Meo^dag Mee^{-dag} 1)
*
*
* i) D^dag phi = (U^{-dag} eta)
* eta'_e = eta_e
* eta'_o = (eta_o - Meo^dag Mee^{-dag} eta_e)
*
* phi_o = D_oo^-dag eta'_o = D_oo^-dag (eta_o - Meo^dag Mee^{-dag} eta_e)
*
* phi_e = D_ee^-dag eta'_e = D_ee^-dag eta_e
*
* Solve:
*
* D_oo D_oo^dag phi_o = D_oo (eta_o - Meo^dag Mee^{-dag} eta_e)
*
* ii)
* phi = L^dag psi => psi = L^-dag phi.
*
* L^{-dag} = ( 1 -Mee^{-dag} Moe^{dag} )
* ( 0 1 )
*
* => sol_e = M_ee^-dag * ( src_e - Moe^dag phi_o )...
* => sol_o = phi_o
*/
///////////////////////////////////////////////////////////////////////////////////////////////////////
// Site diagonal has Mooee on it, but solve the Adjoint system
///////////////////////////////////////////////////////////////////////////////////////////////////////
template<class Field> class SchurRedBlackDiagMooeeDagSolve : public SchurRedBlackBase<Field> {
public:
typedef CheckerBoardedSparseMatrixBase<Field> Matrix;
virtual bool adjoint(void) { return true; }
SchurRedBlackDiagMooeeDagSolve(OperatorFunction<Field> &HermitianRBSolver,
const bool initSubGuess = false,
const bool _solnAsInitGuess = false)
: SchurRedBlackBase<Field> (HermitianRBSolver,initSubGuess,_solnAsInitGuess) {};
//////////////////////////////////////////////////////
// Override RedBlack specialisation
//////////////////////////////////////////////////////
virtual void RedBlackSource(Matrix & _Matrix,const Field &src, Field &src_e,Field &src_o)
{
GridBase *grid = _Matrix.RedBlackGrid();
GridBase *fgrid= _Matrix.Grid();
Field tmp(grid);
Field Mtmp(grid);
pickCheckerboard(Even,src_e,src);
pickCheckerboard(Odd ,src_o,src);
/////////////////////////////////////////////////////
// src_o = (source_o - Moe^dag MeeInvDag source_e)
/////////////////////////////////////////////////////
_Matrix.MooeeInvDag(src_e,tmp); assert( tmp.Checkerboard() ==Even);
_Matrix.MeooeDag (tmp,Mtmp); assert( Mtmp.Checkerboard() ==Odd);
tmp=src_o-Mtmp; assert( tmp.Checkerboard() ==Odd);
// get the right Mpc
SchurDiagMooeeOperator<Matrix,Field> _HermOpEO(_Matrix);
_HermOpEO.Mpc(tmp,src_o); assert(src_o.Checkerboard() ==Odd);
}
virtual void RedBlackSolve (Matrix & _Matrix,const Field &src_o, Field &sol_o)
{
SchurDiagMooeeDagOperator<Matrix,Field> _HermOpEO(_Matrix);
this->_HermitianRBSolver(_HermOpEO,src_o,sol_o);
};
virtual void RedBlackSolve (Matrix & _Matrix,const std::vector<Field> &src_o, std::vector<Field> &sol_o)
{
SchurDiagMooeeDagOperator<Matrix,Field> _HermOpEO(_Matrix);
this->_HermitianRBSolver(_HermOpEO,src_o,sol_o);
}
virtual void RedBlackSolution(Matrix & _Matrix,const Field &sol_o, const Field &src_e,Field &sol)
{
GridBase *grid = _Matrix.RedBlackGrid();
GridBase *fgrid= _Matrix.Grid();
Field sol_e(grid);
Field tmp(grid);
///////////////////////////////////////////////////
// sol_e = M_ee^-dag * ( src_e - Moe^dag phi_o )...
// sol_o = phi_o
///////////////////////////////////////////////////
_Matrix.MeooeDag(sol_o,tmp); assert(tmp.Checkerboard()==Even);
tmp = src_e-tmp; assert(tmp.Checkerboard()==Even);
_Matrix.MooeeInvDag(tmp,sol_e); assert(sol_e.Checkerboard()==Even);
setCheckerboard(sol,sol_e); assert( sol_e.Checkerboard() ==Even);
setCheckerboard(sol,sol_o); assert( sol_o.Checkerboard() ==Odd );
}
};
} }
#endif #endif

2
Grid/qcd/action/ActionCore.h
View File

@ -58,6 +58,8 @@ NAMESPACE_CHECK(Scalar);
//////////////////////////////////////////// ////////////////////////////////////////////
// Utility functions // Utility functions
//////////////////////////////////////////// ////////////////////////////////////////////
#include <Grid/qcd/action/momentum/MomentumFilter.h>
//#include <Grid/qcd/action/momentum/DDHMCfilter.h>
#include <Grid/qcd/utils/Metric.h> #include <Grid/qcd/utils/Metric.h>
NAMESPACE_CHECK(Metric); NAMESPACE_CHECK(Metric);
#include <Grid/qcd/utils/CovariantLaplacian.h> #include <Grid/qcd/utils/CovariantLaplacian.h>

2
Grid/qcd/action/fermion/CayleyFermion5D.h
View File

@ -60,6 +60,8 @@ public:
/////////////////////////////////////////////////////////////// ///////////////////////////////////////////////////////////////
virtual void Dminus(const FermionField &psi, FermionField &chi); virtual void Dminus(const FermionField &psi, FermionField &chi);
virtual void DminusDag(const FermionField &psi, FermionField &chi); virtual void DminusDag(const FermionField &psi, FermionField &chi);
virtual void ImportFourDimPseudoFermion(const FermionField &input,FermionField &imported);
virtual void ExportFourDimPseudoFermion(const FermionField &solution,FermionField &exported);
virtual void ExportPhysicalFermionSolution(const FermionField &solution5d,FermionField &exported4d); virtual void ExportPhysicalFermionSolution(const FermionField &solution5d,FermionField &exported4d);
virtual void ExportPhysicalFermionSource(const FermionField &solution5d, FermionField &exported4d); virtual void ExportPhysicalFermionSource(const FermionField &solution5d, FermionField &exported4d);
virtual void ImportPhysicalFermionSource(const FermionField &input4d,FermionField &imported5d); virtual void ImportPhysicalFermionSource(const FermionField &input4d,FermionField &imported5d);

10
Grid/qcd/action/fermion/FermionOperator.h
View File

@ -171,6 +171,16 @@ public:
/////////////////////////////////////////////// ///////////////////////////////////////////////
virtual void Dminus(const FermionField &psi, FermionField &chi) { chi=psi; } virtual void Dminus(const FermionField &psi, FermionField &chi) { chi=psi; }
virtual void DminusDag(const FermionField &psi, FermionField &chi) { chi=psi; } virtual void DminusDag(const FermionField &psi, FermionField &chi) { chi=psi; }
virtual void ImportFourDimPseudoFermion(const FermionField &input,FermionField &imported)
{
imported = input;
};
virtual void ExportFourDimPseudoFermion(const FermionField &solution,FermionField &exported)
{
exported=solution;
};
virtual void ImportPhysicalFermionSource(const FermionField &input,FermionField &imported) virtual void ImportPhysicalFermionSource(const FermionField &input,FermionField &imported)
{ {
imported = input; imported = input;

32
Grid/qcd/action/fermion/implementation/CayleyFermion5DImplementation.h
View File

@ -112,7 +112,6 @@ void CayleyFermion5D<Impl>::ImportUnphysicalFermion(const FermionField &input4d,
axpby_ssp_pminus(tmp, 0., tmp, 1., tmp, Ls-1, Ls-1); axpby_ssp_pminus(tmp, 0., tmp, 1., tmp, Ls-1, Ls-1);
imported5d=tmp; imported5d=tmp;
} }
template<class Impl> template<class Impl>
void CayleyFermion5D<Impl>::ImportPhysicalFermionSource(const FermionField &input4d,FermionField &imported5d) void CayleyFermion5D<Impl>::ImportPhysicalFermionSource(const FermionField &input4d,FermionField &imported5d)
{ {
@ -127,6 +126,37 @@ void CayleyFermion5D<Impl>::ImportPhysicalFermionSource(const FermionField &inpu
axpby_ssp_pminus(tmp, 0., tmp, 1., tmp, Ls-1, Ls-1); axpby_ssp_pminus(tmp, 0., tmp, 1., tmp, Ls-1, Ls-1);
Dminus(tmp,imported5d); Dminus(tmp,imported5d);
} }
////////////////////////////////////////////////////
// Added for fourD pseudofermion det estimation
////////////////////////////////////////////////////
template<class Impl>
void CayleyFermion5D<Impl>::ImportFourDimPseudoFermion(const FermionField &input4d,FermionField &imported5d)
{
int Ls = this->Ls;
FermionField tmp(this->FermionGrid());
conformable(imported5d.Grid(),this->FermionGrid());
conformable(input4d.Grid() ,this->GaugeGrid());
tmp = Zero();
InsertSlice(input4d, tmp, 0 , 0);
InsertSlice(input4d, tmp, Ls-1, 0);
axpby_ssp_pminus(tmp, 0., tmp, 1., tmp, 0, 0);
axpby_ssp_pplus (tmp, 0., tmp, 1., tmp, Ls-1, Ls-1);
imported5d=tmp;
}
template<class Impl>
void CayleyFermion5D<Impl>::ExportFourDimPseudoFermion(const FermionField &solution5d,FermionField &exported4d)
{
int Ls = this->Ls;
FermionField tmp(this->FermionGrid());
tmp = solution5d;
conformable(solution5d.Grid(),this->FermionGrid());
conformable(exported4d.Grid(),this->GaugeGrid());
axpby_ssp_pminus(tmp, 0., solution5d, 1., solution5d, 0, 0);
axpby_ssp_pplus (tmp, 1., tmp , 1., solution5d, 0, Ls-1);
ExtractSlice(exported4d, tmp, 0, 0);
}
// Dminus
template<class Impl> template<class Impl>
void CayleyFermion5D<Impl>::Dminus(const FermionField &psi, FermionField &chi) void CayleyFermion5D<Impl>::Dminus(const FermionField &psi, FermionField &chi)
{ {

79
Grid/qcd/action/momentum/DDHMCfilter.h Normal file
View File

@ -0,0 +1,79 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: ./lib/qcd/hmc/DDHMC.h
Copyright (C) 2021
Author: Peter Boyle <paboyle@ph.ed.ac.uk>
Author: Christopher Kelly
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 */
////////////////////////////////////////////////////
// DDHMC filter with sub-block size B[mu]
////////////////////////////////////////////////////
template<typename MomentaField>
struct DDHMCFilter: public MomentumFilterBase<MomentaField>
{
typedef typename MomentaField::vector_type vector_type; //SIMD-vectorized complex type
typedef typename MomentaField::scalar_type scalar_type; //scalar complex type
typedef iVector<iScalar<iScalar<vector_type> >, Nd > LorentzScalarType; //complex phase for each site/direction
typedef iScalar<iScalar<iScalar<vector_type> > > ScalarType; //complex phase for each site
typedef Lattice<LorentzScalarType> LatticeLorentzScalarType;
typedef Lattice<ScalarType> LatticeScalarType;
Coordinate Block;
// Could also zero links in plane like Luscher advocates.
DDHMCFilter(const Coordinate &_Block): Block(_Block){}
void applyFilter(MomentaField &P) const override
{
GridBase *grid = P.Grid();
LatticeScalarType mask_mu(grid);
LatticeLorentzScalarType mask(grid);
////////////////////////////////////////////////////
// Zero strictly links crossing between domains
// Luscher also zeroes links in plane of domain boundaries
// Keeping interior only. This prevents force from plaquettes
// crossing domains and keeps whole MD trajectory local.
// Step further than I want to go.
////////////////////////////////////////////////////
ScalarType zz = scalar_type(0.0);
ScalarType one= scalar_type(1.0);
LatticeInteger coor(grid);
for(int mu=0;mu<Nd;mu++) {
LatticeCoordinate(coor,mu);
mask_mu = where(mod(coor,Block[mu])==Block[mu]-1,zz,one);
PokeIndex<LorentzIndex>(mask, mask_mu, mu);
}
}
};

0
Grid/qcd/hmc/integrators/MomentumFilter.h → Grid/qcd/action/momentum/MomentumFilter.h
View File

346
Grid/qcd/action/pseudofermion/DomainDecomposedBoundary.h Normal file
View File

@ -0,0 +1,346 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: ./lib/qcd/action/pseudofermion/TwoFlavourRatio.h
Copyright (C) 2015
Author: Peter Boyle <paboyle@ph.ed.ac.uk>
Author: Peter Boyle <peterboyle@Peters-MacBook-Pro-2.local>
Author: paboyle <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 */
#pragma once
NAMESPACE_BEGIN(Grid);
///////////////////////////////////////
// Two flavour ratio
///////////////////////////////////////
template<class Impl>
class DomainDecomposedBoundary {
public:
INHERIT_IMPL_TYPES(Impl);
typedef typename GaugeField::vector_type vector_type; //SIMD-vectorized complex type
typedef typename GaugeField::scalar_type scalar_type; //scalar complex type
typedef iVector<iScalar<iScalar<vector_type> >, Nd > LorentzScalarType; //complex phase for each site/direction
typedef iScalar<iScalar<iScalar<vector_type> > > ScalarType; //complex phase for each site
typedef Lattice<LorentzScalarType> LatticeLorentzScalarType;
typedef Lattice<ScalarType> LatticeScalarType;
Coordinate Block;
DDHMCFilter Filter;
const int Omega=0;
const int OmegaBar=1;
void ProjectBoundaryBothDomains (FermionField &f,int sgn)
{
assert((sgn==1)||(sgn==-1));
Gamma::Algebra Gmu [] = {
Gamma::Algebra::GammaX,
Gamma::Algebra::GammaY,
Gamma::Algebra::GammaZ,
Gamma::Algebra::GammaT
};
GridBase *grid = f.Grid();
LatticeInteger coor(grid);
LatticeInteger face(grid);
LatticeInteger nface(grid); nface=Zero();
ComplexField zz(grid); zz=Zero();
FermionField projected(grid); projected=Zero();
FermionField sp_proj (grid);
int dims = grid->Nd();
int isDWF= (dims==Nd+1);
assert((dims==Nd)||(dims==Nd+1));
for(int mu=0;mu<Nd;mu++){
// need to worry about DWF 5th dim first
// Could extend to domain decompose in FIFTH dimension.
// With chiral projectors here
LatticeCoordinate(coor,mu+isDWF);
face = (mod(coor,Block[mu]) == 0 );
nface = nface + face;
// Lower face receives (1-gamma)/2 in normal forward hopping term
sp_proj = 0.5*(f-sgn*Gamma(Gmu[mu])*f)
projected= where(face==cb,f,projected);
face = (mod(coor,Block[mu]) == Block[mu]-1 );
nface = nface + face;
// Upper face receives (1+gamma)/2 in normal backward hopping term
sp_proj = 0.5*(f+sgn*Gamma(Gmu[mu])*f)
projected= where(face==cb,f,projected);
}
// Keep the spin projected faces where nface==1 and initial Zero() where nface==0.
projected = where(nface>1,f,projected);
}
void ProjectDomain(FermionField &f,int cb)
{
GridBase *grid = f.Grid();
ComplexField zz(grid); zz=Zero();
LatticeInteger coor(grid);
LatticeInteger domaincb(grid); domaincb=Zero();
for(int d=0;d<grid->Nd();d++){
LatticeCoordinate(coor,mu);
domaincb = domaincb + div(coor,Block[d]);
}
f = where(mod(domaincb,2)==cb,f,zz);
};
void ProjectOmegaBar (FermionField &f) {ProjectDomain(f,OmegaBar);}
void ProjectOmega (FermionField &f) {ProjectDomain(f,Omega);}
// See my notes(!).
// Notation: Following Luscher, we introduce projectors $\hPdb$ with both spinor and space structure
// projecting all spinor elements in $\Omega$ connected by $\Ddb$ to $\bar{\Omega}$,
void ProjectBoundaryBar(FermionField &f)
{
ProjectBoundaryBothDomains(f);
ProjectOmega(f);
}
// and $\hPd$ projecting all spinor elements in $\bar{\Omega}$ connected by $\Dd$ to $\Omega$.
void ProjectBoundary (FermionField &f)
{
ProjectBoundaryBothDomains(f);
ProjectOmegaBar(f);
};
void dBoundary (FermionOperator<Impl> &Op,FermionField &in,FermionField &out)
{
FermionField tmp(in);
ProjectOmegaBar(tmp);
Op.M(tmp,out);
ProjectOmega(out);
};
void dBoundaryBar (FermionOperator<Impl> &Op,FermionField &in,FermionField &out)
{
FermionField tmp(in);
ProjectOmega(tmp);
Op.M(tmp,out);
ProjectOmegaBar(out);
};
void dOmega (FermionOperator<Impl> &Op,FermionField &in,FermionField &out)
{
FermionField tmp(in);
ProjectOmega(tmp);
Op.M(tmp,out);
ProjectOmega(out);
};
void dOmegaBar (FermionOperator<Impl> &Op,FermionField &in,FermionField &out)
{
FermionField tmp(in);
ProjectOmegaBar(tmp);
Op.M(tmp,out);
ProjectOmegaBar(out);
};
void SolveOmega (FermionOperator<Impl> &Op,FermionField &in,FermionField &out){ assert(0); };
void SolveOmegaBar(FermionOperator<Impl> &Op,FermionField &in,FermionField &out){ assert(0); };
void SolveOmegaAndOmegaBar(FermionOperator<Impl> &Op,FermionField &in,FermionField &out){ assert(0); };
void dInverse (FermionOperator<Impl> &Op,FermionField &in,FermionField &out){ assert(0); };
// R = Pdbar - Pdbar DomegaInv Dd DomegabarInv Ddbar
void R(FermionOperator<Impl> &Op,FermionOperator<Impl> &OpDirichlet,FermionField &in,FermionField &out)
{
FermionField tmp1(Op.FermionGrid());
FermionField tmp2(Op.FermionGrid());
dBoundaryBar(Op,in,tmp1);
SolveOmegaBar(OpDirichlet,tmp1,tmp2); // 1/2 cost
dBoundary(Op,tmp2,tmp1);
SolveOmega(OpDirichlet,tmp1,tmp2); // 1/2 cost
out = in - tmp2 ;
ProjectBoundaryBar(out);
};
// R = Pdbar - Pdbar Dinv Ddbar
void Rinverse(FermionField &in,FermionField &out)
{
FermionField tmp1(NumOp.FermionGrid());
out = in;
ProjectBoundaryBar(out);
dInverse(out,tmp1);
ProjectBoundaryBar(tmp1);
out = out -tmp1;
};
}
template<class Impl>
class DomainDecomposedBoundaryPseudoFermionAction : public Action<typename Impl::GaugeField> {
public:
INHERIT_IMPL_TYPES(Impl);
private:
FermionOperator<Impl> & NumOp;// the basic operator
FermionOperator<Impl> & DenOp;// the basic operator
FermionOperator<Impl> & NumOpDirichlet;// the basic operator
FermionOperator<Impl> & DenOpDirichlet;// the basic operator
OperatorFunction<FermionField> &DerivativeSolver;
OperatorFunction<FermionField> &ActionSolver;
FermionField Phi; // the pseudo fermion field for this trajectory
public:
DomainBoundaryPseudoFermionAction(FermionOperator<Impl> &_NumOp,
FermionOperator<Impl> &_DenOp,
FermionOperator<Impl> &_NumOpDirichlet,
FermionOperator<Impl> &_DenOpDirichlet,
OperatorFunction<FermionField> & DS,
OperatorFunction<FermionField> & AS,
Coordinate &_Block
) : NumOp(_NumOp),
DenOp(_DenOp),
DerivativeSolver(DS),
ActionSolver(AS),
Phi(_NumOp.FermionGrid()),
Block(_Block)
// LinkFilter(Block)
{};
virtual std::string action_name(){return "DomainBoundaryPseudoFermionRatioAction";}
virtual std::string LogParameters(){
std::stringstream sstream;
sstream << GridLogMessage << "["<<action_name()<<"] Block "<<_Block << std::endl;
return sstream.str();
}
virtual void refresh(const GaugeField &U, GridParallelRNG& pRNG)
{
// P(phi) = e^{- phi^dag V (MdagM)^-1 Vdag phi}
//
// NumOp == V
// DenOp == M
//
// Take phi = Vdag^{-1} Mdag eta ; eta = Mdag^{-1} Vdag Phi
//
// P(eta) = e^{- eta^dag eta}
//
// e^{x^2/2 sig^2} => sig^2 = 0.5.
//
// So eta should be of width sig = 1/sqrt(2) and must multiply by 0.707....
//
RealD scale = std::sqrt(0.5);
FermionField eta(NumOp.FermionGrid());
FermionField tmp(NumOp.FermionGrid());
gaussian(pRNG,eta);
ProjectBoundary(eta);
NumOp.ImportGauge(U);
DenOp.ImportGauge(U);
// Note: this hard codes normal equations type solvers; alternate implementation needed for
// non-herm style solvers.
MdagMLinearOperator<FermionOperator<Impl> ,FermionField> MdagMOp(NumOp);
DenOp.Mdag(eta,Phi); // Mdag eta
tmp = Zero();
ActionSolver(MdagMOp,Phi,tmp); // (VdagV)^-1 Mdag eta = V^-1 Vdag^-1 Mdag eta
NumOp.M(tmp,Phi); // Vdag^-1 Mdag eta
Phi=Phi*scale;
};
//////////////////////////////////////////////////////
// S = phi^dag V (Mdag M)^-1 Vdag phi
//////////////////////////////////////////////////////
virtual RealD S(const GaugeField &U) {
NumOp.ImportGauge(U);
DenOp.ImportGauge(U);
FermionField X(NumOp.FermionGrid());
FermionField Y(NumOp.FermionGrid());
MdagMLinearOperator<FermionOperator<Impl> ,FermionField> MdagMOp(DenOp);
NumOp.Mdag(Phi,Y); // Y= Vdag phi
X=Zero();
ActionSolver(MdagMOp,Y,X); // X= (MdagM)^-1 Vdag phi
DenOp.M(X,Y); // Y= Mdag^-1 Vdag phi
RealD action = norm2(Y);
return action;
};
//////////////////////////////////////////////////////
// dS/du = phi^dag dV (Mdag M)^-1 V^dag phi
// - phi^dag V (Mdag M)^-1 [ Mdag dM + dMdag M ] (Mdag M)^-1 V^dag phi
// + phi^dag V (Mdag M)^-1 dV^dag phi
//////////////////////////////////////////////////////
virtual void deriv(const GaugeField &U,GaugeField & dSdU) {
NumOp.ImportGauge(U);
DenOp.ImportGauge(U);
MdagMLinearOperator<FermionOperator<Impl> ,FermionField> MdagMOp(DenOp);
FermionField X(NumOp.FermionGrid());
FermionField Y(NumOp.FermionGrid());
GaugeField force(NumOp.GaugeGrid());
//Y=Vdag phi
//X = (Mdag M)^-1 V^dag phi
//Y = (Mdag)^-1 V^dag phi
NumOp.Mdag(Phi,Y); // Y= Vdag phi
X=Zero();
DerivativeSolver(MdagMOp,Y,X); // X= (MdagM)^-1 Vdag phi
DenOp.M(X,Y); // Y= Mdag^-1 Vdag phi
// phi^dag V (Mdag M)^-1 dV^dag phi
NumOp.MDeriv(force , X, Phi, DaggerYes ); dSdU=force;
// phi^dag dV (Mdag M)^-1 V^dag phi
NumOp.MDeriv(force , Phi, X ,DaggerNo ); dSdU=dSdU+force;
// - 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
DenOp.MDeriv(force,Y,X,DaggerNo); dSdU=dSdU-force;
DenOp.MDeriv(force,X,Y,DaggerYes); dSdU=dSdU-force;
dSdU *= -1.0;
//dSdU = - Ta(dSdU);
};
};
NAMESPACE_END(Grid);

33
Grid/qcd/action/pseudofermion/TwoFlavourRatio.h
View File

@ -97,7 +97,19 @@ public:
tmp = Zero(); tmp = Zero();
ActionSolver(MdagMOp,Phi,tmp); // (VdagV)^-1 Mdag eta = V^-1 Vdag^-1 Mdag eta ActionSolver(MdagMOp,Phi,tmp); // (VdagV)^-1 Mdag eta = V^-1 Vdag^-1 Mdag eta
NumOp.M(tmp,Phi); // Vdag^-1 Mdag eta NumOp.M(tmp,Phi); // Vdag^-1 Mdag eta
#define FILTER
#ifdef FILTER
Integer OrthogDir=0;
Integer plane=0;
if ( getenv("DIR") ) OrthogDir = atoi(getenv("DIR"));
if ( getenv("COOR") ) plane = atoi(getenv("COOR"));
std::cout << " *** PseudoFermion FILTER DIR " <<OrthogDir << " plane "<<plane<<std::endl;
Lattice<iScalar<vInteger> > coor(NumOp.FermionGrid());
LatticeCoordinate(coor,OrthogDir);
tmp = Zero();
Phi = where(coor==plane,Phi,tmp);
#endif
Phi=Phi*scale; Phi=Phi*scale;
}; };
@ -165,6 +177,25 @@ public:
dSdU *= -1.0; dSdU *= -1.0;
//dSdU = - Ta(dSdU); //dSdU = - Ta(dSdU);
#ifdef FILTER
std::cout <<" In force "<<std::endl;
force = dSdU;
int mu=0;
std::cout << " FORCE mu " <<mu<<" L2 "<< norm2(force)<< " Linf " << maxLocalNorm2(force)<<std::endl;
int plane=0;
if ( getenv("COOR") ) plane = atoi(getenv("COOR"));
Lattice<iScalar<vInteger> > coor(NumOp.GaugeGrid());
LatticeCoordinate(coor,mu);
int L = NumOp.GaugeGrid()->FullDimensions()[mu];
for (Integer p=0;p<L;p++) {
force = Zero();
force = where(coor==p,dSdU,force);
std::cout << " FORCE mu " <<mu<<" PF plane "<<plane<<" T= " <<p<<" L2 "<< norm2(force)<< " Linf " << maxLocalNorm2(force)<<std::endl;
}
exit(0);
#endif
}; };
}; };

197
Grid/qcd/action/pseudofermion/TwoFlavourRatio4DPseudoFermion.h Normal file
View File

@ -0,0 +1,197 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: ./lib/qcd/action/pseudofermion/TwoFlavourRatio.h
Copyright (C) 2015
Author: Peter Boyle <paboyle@ph.ed.ac.uk>
Author: Peter Boyle <peterboyle@Peters-MacBook-Pro-2.local>
Author: paboyle <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 */
#pragma once
NAMESPACE_BEGIN(Grid);
///////////////////////////////////////
// Two flavour ratio
///////////////////////////////////////
template<class Impl>
class TwoFlavourRatio4DPseudoFermionAction : public Action<typename Impl::GaugeField> {
public:
INHERIT_IMPL_TYPES(Impl);
private:
FermionOperator<Impl> & NumOp;// the basic operator
FermionOperator<Impl> & DenOp;// the basic operator
OperatorFunction<FermionField> &DerivativeSolver;
OperatorFunction<FermionField> &ActionSolver;
FermionField phi4; // the pseudo fermion field for this trajectory
public:
TwoFlavourRatio4DPseudoFermionAction(FermionOperator<Impl> &_NumOp,
FermionOperator<Impl> &_DenOp,
OperatorFunction<FermionField> & DS,
OperatorFunction<FermionField> & AS
) : NumOp(_NumOp),
DenOp(_DenOp),
DerivativeSolver(DS),
ActionSolver(AS),
phi4(_NumOp.GaugeGrid())
{};
virtual std::string action_name(){return "TwoFlavourRatio4DPseudoFermionAction";}
virtual std::string LogParameters(){
std::stringstream sstream;
sstream << GridLogMessage << "["<<action_name()<<"] has no parameters" << std::endl;
return sstream.str();
}
virtual void refresh(const GaugeField &U, GridSerialRNG &sRNG, GridParallelRNG& pRNG) {
// P(phi) = e^{- phi^dag (V^dag M^-dag)_11 (M^-1 V)_11 phi}
//
// NumOp == V
// DenOp == M
//
// Take phi = (V^{-1} M)_11 eta ; eta = (M^{-1} V)_11 Phi
//
// P(eta) = e^{- eta^dag eta}
//
// e^{x^2/2 sig^2} => sig^2 = 0.5.
//
// So eta should be of width sig = 1/sqrt(2) and must multiply by 0.707....
//
RealD scale = std::sqrt(0.5);
FermionField eta4(NumOp.GaugeGrid());
FermionField eta5(NumOp.FermionGrid());
FermionField tmp(NumOp.FermionGrid());
FermionField phi5(NumOp.FermionGrid());
gaussian(pRNG,eta4);
NumOp.ImportFourDimPseudoFermion(eta4,eta5);
NumOp.ImportGauge(U);
DenOp.ImportGauge(U);
MdagMLinearOperator<FermionOperator<Impl> ,FermionField> MdagMOp(NumOp);
DenOp.M(eta5,phi5); // M eta
NumOp.Mdag(phi5,tmp); // Vdag M eta
phi5 = Zero();
ActionSolver(MdagMOp,tmp,phi5); // (VdagV)^-1 M eta = V^-1 Vdag^-1 Vdag M eta = V^-1 M eta
phi5=phi5*scale;
// Project to 4d
NumOp.ExportFourDimPseudoFermion(phi5,phi4);
};
//////////////////////////////////////////////////////
// S = phi^dag (V^dag M^-dag)_11 (M^-1 V)_11 phi
//////////////////////////////////////////////////////
virtual RealD S(const GaugeField &U) {
NumOp.ImportGauge(U);
DenOp.ImportGauge(U);
FermionField Y4(NumOp.GaugeGrid());
FermionField X(NumOp.FermionGrid());
FermionField Y(NumOp.FermionGrid());
FermionField phi5(NumOp.FermionGrid());
MdagMLinearOperator<FermionOperator<Impl> ,FermionField> MdagMOp(DenOp);
NumOp.ImportFourDimPseudoFermion(phi4,phi5);
NumOp.M(phi5,Y); // Y= V phi
DenOp.Mdag(Y,X); // X= Mdag V phi
Y=Zero();
ActionSolver(MdagMOp,X,Y); // Y= (MdagM)^-1 Mdag Vdag phi = M^-1 V phi
NumOp.ExportFourDimPseudoFermion(Y,Y4);
RealD action = norm2(Y4);
return action;
};
//////////////////////////////////////////////////////
// dS/du = 2 Re phi^dag (V^dag M^-dag)_11 (M^-1 d V)_11 phi
// - 2 Re phi^dag (dV^dag M^-dag)_11 (M^-1 dM M^-1 V)_11 phi
//////////////////////////////////////////////////////
virtual void deriv(const GaugeField &U,GaugeField & dSdU) {
NumOp.ImportGauge(U);
DenOp.ImportGauge(U);
MdagMLinearOperator<FermionOperator<Impl> ,FermionField> MdagMOp(DenOp);
FermionField X(NumOp.FermionGrid());
FermionField Y(NumOp.FermionGrid());
FermionField phi(NumOp.FermionGrid());
FermionField Vphi(NumOp.FermionGrid());
FermionField MinvVphi(NumOp.FermionGrid());
FermionField tmp4(NumOp.GaugeGrid());
FermionField MdagInvMinvVphi(NumOp.FermionGrid());
GaugeField force(NumOp.GaugeGrid());
//Y=V phi
//X = (Mdag V phi
//Y = (Mdag M)^-1 Mdag V phi = M^-1 V Phi
NumOp.ImportFourDimPseudoFermion(phi4,phi);
NumOp.M(phi,Vphi); // V phi
DenOp.Mdag(Vphi,X); // X= Mdag V phi
Y=Zero();
DerivativeSolver(MdagMOp,X,MinvVphi);// M^-1 V phi
// Projects onto the physical space and back
NumOp.ExportFourDimPseudoFermion(MinvVphi,tmp4);
NumOp.ImportFourDimPseudoFermion(tmp4,Y);
X=Zero();
DerivativeSolver(MdagMOp,Y,X);// X = (MdagM)^-1 proj M^-1 V phi
DenOp.M(X,MdagInvMinvVphi);
// phi^dag (Vdag Mdag^-1) (M^-1 dV) phi
NumOp.MDeriv(force ,MdagInvMinvVphi , phi, DaggerNo ); dSdU=force;
// phi^dag (dVdag Mdag^-1) (M^-1 V) phi
NumOp.MDeriv(force , phi, MdagInvMinvVphi ,DaggerYes ); dSdU=dSdU+force;
// - 2 Re phi^dag (dV^dag M^-dag)_11 (M^-1 dM M^-1 V)_11 phi
DenOp.MDeriv(force,MdagInvMinvVphi,MinvVphi,DaggerNo); dSdU=dSdU-force;
DenOp.MDeriv(force,MinvVphi,MdagInvMinvVphi,DaggerYes); dSdU=dSdU-force;
dSdU *= -1.0;
//dSdU = - Ta(dSdU);
};
};
NAMESPACE_END(Grid);

203
Grid/qcd/action/pseudofermion/TwoFlavourRatioEO4DPseudoFermion.h Normal file
View File

@ -0,0 +1,203 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: ./lib/qcd/action/pseudofermion/TwoFlavourRatio.h
Copyright (C) 2015
Author: Peter Boyle <paboyle@ph.ed.ac.uk>
Author: Peter Boyle <peterboyle@Peters-MacBook-Pro-2.local>
Author: paboyle <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 */
#pragma once
NAMESPACE_BEGIN(Grid);
///////////////////////////////////////
// Two flavour ratio
///////////////////////////////////////
template<class Impl>
class TwoFlavourRatioEO4DPseudoFermionAction : public Action<typename Impl::GaugeField> {
public:
INHERIT_IMPL_TYPES(Impl);
private:
typedef FermionOperator<Impl> FermOp;
FermionOperator<Impl> & NumOp;// the basic operator
FermionOperator<Impl> & DenOp;// the basic operator
OperatorFunction<FermionField> &DerivativeSolver;
OperatorFunction<FermionField> &DerivativeDagSolver;
OperatorFunction<FermionField> &ActionSolver;
OperatorFunction<FermionField> &HeatbathSolver;
FermionField phi4; // the pseudo fermion field for this trajectory
public:
TwoFlavourRatioEO4DPseudoFermionAction(FermionOperator<Impl> &_NumOp,
FermionOperator<Impl> &_DenOp,
OperatorFunction<FermionField> & DS,
OperatorFunction<FermionField> & AS ) :
TwoFlavourRatioEO4DPseudoFermionAction(_NumOp,_DenOp, DS,DS,AS,AS) {};
TwoFlavourRatioEO4DPseudoFermionAction(FermionOperator<Impl> &_NumOp,
FermionOperator<Impl> &_DenOp,
OperatorFunction<FermionField> & DS,
OperatorFunction<FermionField> & DDS,
OperatorFunction<FermionField> & AS,
OperatorFunction<FermionField> & HS
) : NumOp(_NumOp),
DenOp(_DenOp),
DerivativeSolver(DS),
DerivativeDagSolver(DDS),
ActionSolver(AS),
HeatbathSolver(HS),
phi4(_NumOp.GaugeGrid())
{};
virtual std::string action_name(){return "TwoFlavourRatioEO4DPseudoFermionAction";}
virtual std::string LogParameters(){
std::stringstream sstream;
sstream << GridLogMessage << "["<<action_name()<<"] has no parameters" << std::endl;
return sstream.str();
}
virtual void refresh(const GaugeField &U, GridSerialRNG &sRNG, GridParallelRNG& pRNG) {
// P(phi) = e^{- phi^dag (V^dag M^-dag)_11 (M^-1 V)_11 phi}
//
// NumOp == V
// DenOp == M
//
// Take phi = (V^{-1} M)_11 eta ; eta = (M^{-1} V)_11 Phi
//
// P(eta) = e^{- eta^dag eta}
//
// e^{x^2/2 sig^2} => sig^2 = 0.5.
//
// So eta should be of width sig = 1/sqrt(2) and must multiply by 0.707....
//
RealD scale = std::sqrt(0.5);
FermionField eta4(NumOp.GaugeGrid());
FermionField eta5(NumOp.FermionGrid());
FermionField tmp(NumOp.FermionGrid());
FermionField phi5(NumOp.FermionGrid());
gaussian(pRNG,eta4);
NumOp.ImportFourDimPseudoFermion(eta4,eta5);
NumOp.ImportGauge(U);
DenOp.ImportGauge(U);
SchurRedBlackDiagMooeeSolve<FermionField> PrecSolve(HeatbathSolver);
DenOp.M(eta5,tmp); // M eta
PrecSolve(NumOp,tmp,phi5); // phi = V^-1 M eta
phi5=phi5*scale;
std::cout << GridLogMessage << "4d pf refresh "<< norm2(phi5)<<"\n";
// Project to 4d
NumOp.ExportFourDimPseudoFermion(phi5,phi4);
};
//////////////////////////////////////////////////////
// S = phi^dag (V^dag M^-dag)_11 (M^-1 V)_11 phi
//////////////////////////////////////////////////////
virtual RealD S(const GaugeField &U) {
NumOp.ImportGauge(U);
DenOp.ImportGauge(U);
FermionField Y4(NumOp.GaugeGrid());
FermionField X(NumOp.FermionGrid());
FermionField Y(NumOp.FermionGrid());
FermionField phi5(NumOp.FermionGrid());
MdagMLinearOperator<FermionOperator<Impl> ,FermionField> MdagMOp(DenOp);
SchurRedBlackDiagMooeeSolve<FermionField> PrecSolve(ActionSolver);
NumOp.ImportFourDimPseudoFermion(phi4,phi5);
NumOp.M(phi5,X); // X= V phi
PrecSolve(DenOp,X,Y); // Y= (MdagM)^-1 Mdag Vdag phi = M^-1 V phi
NumOp.ExportFourDimPseudoFermion(Y,Y4);
RealD action = norm2(Y4);
return action;
};
//////////////////////////////////////////////////////
// dS/du = 2 Re phi^dag (V^dag M^-dag)_11 (M^-1 d V)_11 phi
// - 2 Re phi^dag (dV^dag M^-dag)_11 (M^-1 dM M^-1 V)_11 phi
//////////////////////////////////////////////////////
virtual void deriv(const GaugeField &U,GaugeField & dSdU) {
NumOp.ImportGauge(U);
DenOp.ImportGauge(U);
FermionField X(NumOp.FermionGrid());
FermionField Y(NumOp.FermionGrid());
FermionField phi(NumOp.FermionGrid());
FermionField Vphi(NumOp.FermionGrid());
FermionField MinvVphi(NumOp.FermionGrid());
FermionField tmp4(NumOp.GaugeGrid());
FermionField MdagInvMinvVphi(NumOp.FermionGrid());
GaugeField force(NumOp.GaugeGrid());
//Y=V phi
//X = (Mdag V phi
//Y = (Mdag M)^-1 Mdag V phi = M^-1 V Phi
NumOp.ImportFourDimPseudoFermion(phi4,phi);
NumOp.M(phi,Vphi); // V phi
SchurRedBlackDiagMooeeSolve<FermionField> PrecSolve(DerivativeSolver);
PrecSolve(DenOp,Vphi,MinvVphi);// M^-1 V phi
std::cout << GridLogMessage << "4d deriv solve "<< norm2(MinvVphi)<<"\n";
// Projects onto the physical space and back
NumOp.ExportFourDimPseudoFermion(MinvVphi,tmp4);
NumOp.ImportFourDimPseudoFermion(tmp4,Y);
SchurRedBlackDiagMooeeDagSolve<FermionField> PrecDagSolve(DerivativeDagSolver);
// X = proj M^-dag V phi
// Need an adjoint solve
PrecDagSolve(DenOp,Y,MdagInvMinvVphi);
std::cout << GridLogMessage << "4d deriv solve dag "<< norm2(MdagInvMinvVphi)<<"\n";
// phi^dag (Vdag Mdag^-1) (M^-1 dV) phi
NumOp.MDeriv(force ,MdagInvMinvVphi , phi, DaggerNo ); dSdU=force;
// phi^dag (dVdag Mdag^-1) (M^-1 V) phi
NumOp.MDeriv(force , phi, MdagInvMinvVphi ,DaggerYes ); dSdU=dSdU+force;
// - 2 Re phi^dag (dV^dag M^-dag)_11 (M^-1 dM M^-1 V)_11 phi
DenOp.MDeriv(force,MdagInvMinvVphi,MinvVphi,DaggerNo); dSdU=dSdU-force;
DenOp.MDeriv(force,MinvVphi,MdagInvMinvVphi,DaggerYes); dSdU=dSdU-force;
dSdU *= -1.0;
//dSdU = - Ta(dSdU);
};
};
NAMESPACE_END(Grid);

86
Grid/qcd/hmc/HMC.h
View File

@ -34,6 +34,7 @@ directory
* @brief Classes for Hybrid Monte Carlo update * @brief Classes for Hybrid Monte Carlo update
* *
* @author Guido Cossu * @author Guido Cossu
* @author Peter Boyle
*/ */
//-------------------------------------------------------------------- //--------------------------------------------------------------------
#pragma once #pragma once
@ -115,22 +116,17 @@ private:
random(sRNG, rn_test); random(sRNG, rn_test);
std::cout << GridLogHMC std::cout << GridLogHMC << "--------------------------------------------------\n";
<< "--------------------------------------------------\n"; std::cout << GridLogHMC << "exp(-dH) = " << prob << " Random = " << rn_test << "\n";
std::cout << GridLogHMC << "exp(-dH) = " << prob std::cout << GridLogHMC << "Acc. Probability = " << ((prob < 1.0) ? prob : 1.0) << "\n";
<< " Random = " << rn_test << "\n";
std::cout << GridLogHMC
<< "Acc. Probability = " << ((prob < 1.0) ? prob : 1.0) << "\n";
if ((prob > 1.0) || (rn_test <= prob)) { // accepted if ((prob > 1.0) || (rn_test <= prob)) { // accepted
std::cout << GridLogHMC << "Metropolis_test -- ACCEPTED\n"; std::cout << GridLogHMC << "Metropolis_test -- ACCEPTED\n";
std::cout << GridLogHMC std::cout << GridLogHMC << "--------------------------------------------------\n";
<< "--------------------------------------------------\n";
return true; return true;
} else { // rejected } else { // rejected
std::cout << GridLogHMC << "Metropolis_test -- REJECTED\n"; std::cout << GridLogHMC << "Metropolis_test -- REJECTED\n";
std::cout << GridLogHMC std::cout << GridLogHMC << "--------------------------------------------------\n";
<< "--------------------------------------------------\n";
return false; return false;
} }
} }
@ -139,18 +135,63 @@ private:
// Evolution // Evolution
///////////////////////////////////////////////////////// /////////////////////////////////////////////////////////
RealD evolve_hmc_step(Field &U) { RealD evolve_hmc_step(Field &U) {
TheIntegrator.refresh(U, sRNG, pRNG); // set U and initialize P and phi's
RealD H0 = TheIntegrator.S(U); // initial state action GridBase *Grid = U.Grid();
//////////////////////////////////////////////////////////////////////////////////////////////////////
// Mainly for DDHMC perform a random translation of U modulo volume
//////////////////////////////////////////////////////////////////////////////////////////////////////
std::cout << GridLogMessage << "--------------------------------------------------\n";
std::cout << GridLogMessage << "Random shifting gauge field by [";
for(int d=0;d<Grid->Nd();d++) {
int L = Grid->GlobalDimensions()[d];
RealD rn_uniform; random(sRNG, rn_uniform);
int shift = (int) (rn_uniform*L);
std::cout << shift;
if(d<Grid->Nd()-1) std::cout <<",";
else std::cout <<"]\n";
U = Cshift(U,d,shift);
}
std::cout << GridLogMessage << "--------------------------------------------------\n";
//////////////////////////////////////////////////////////////////////////////////////////////////////
// set U and initialize P and phi's
//////////////////////////////////////////////////////////////////////////////////////////////////////
std::cout << GridLogMessage << "--------------------------------------------------\n";
std::cout << GridLogMessage << "Refresh momenta and pseudofermions";
TheIntegrator.refresh(U, sRNG, pRNG);
std::cout << GridLogMessage << "--------------------------------------------------\n";
//////////////////////////////////////////////////////////////////////////////////////////////////////
// initial state action
//////////////////////////////////////////////////////////////////////////////////////////////////////
std::cout << GridLogMessage << "--------------------------------------------------\n";
std::cout << GridLogMessage << "Compute initial action";
RealD H0 = TheIntegrator.S(U);
std::cout << GridLogMessage << "--------------------------------------------------\n";
std::streamsize current_precision = std::cout.precision(); std::streamsize current_precision = std::cout.precision();
std::cout.precision(15); std::cout.precision(15);
std::cout << GridLogHMC << "Total H before trajectory = " << H0 << "\n"; std::cout << GridLogHMC << "Total H before trajectory = " << H0 << "\n";
std::cout.precision(current_precision); std::cout.precision(current_precision);
std::cout << GridLogMessage << "--------------------------------------------------\n";
std::cout << GridLogMessage << " Molecular Dynamics evolution ";
TheIntegrator.integrate(U); TheIntegrator.integrate(U);
std::cout << GridLogMessage << "--------------------------------------------------\n";
RealD H1 = TheIntegrator.S(U); // updated state action //////////////////////////////////////////////////////////////////////////////////////////////////////
// updated state action
//////////////////////////////////////////////////////////////////////////////////////////////////////
std::cout << GridLogMessage << "--------------------------------------------------\n";
std::cout << GridLogMessage << "Compute final action";
RealD H1 = TheIntegrator.S(U);
std::cout << GridLogMessage << "--------------------------------------------------\n";
/////////////////////////////////////////////////////////// ///////////////////////////////////////////////////////////
if(0){ if(0){
@ -163,17 +204,19 @@ private:
} }
/////////////////////////////////////////////////////////// ///////////////////////////////////////////////////////////
std::cout.precision(15); std::cout.precision(15);
<<<<<<< HEAD
std::cout << GridLogHMC << "Total H after trajectory = " << H1 std::cout << GridLogHMC << "Total H after trajectory = " << H1
<< " dH = " << H1 - H0 << "\n"; << " dH = " << H1 - H0 << "\n";
=======
std::cout << GridLogMessage << "--------------------------------------------------\n";
std::cout << GridLogMessage << "Total H after trajectory = " << H1 << " dH = " << H1 - H0 << "\n";
std::cout << GridLogMessage << "--------------------------------------------------\n";
>>>>>>> bd181b94819600174e91b331e24f71598b3e7bb9
std::cout.precision(current_precision); std::cout.precision(current_precision);
return (H1 - H0); return (H1 - H0);
} }
public: public:
///////////////////////////////////////// /////////////////////////////////////////
@ -195,8 +238,15 @@ public:
// Actual updates (evolve a copy Ucopy then copy back eventually) // Actual updates (evolve a copy Ucopy then copy back eventually)
unsigned int FinalTrajectory = Params.Trajectories + Params.NoMetropolisUntil + Params.StartTrajectory; unsigned int FinalTrajectory = Params.Trajectories + Params.NoMetropolisUntil + Params.StartTrajectory;
for (int traj = Params.StartTrajectory; traj < FinalTrajectory; ++traj) { for (int traj = Params.StartTrajectory; traj < FinalTrajectory; ++traj) {
<<<<<<< HEAD
std::cout << GridLogHMC << "-- # Trajectory = " << traj << "\n"; std::cout << GridLogHMC << "-- # Trajectory = " << traj << "\n";
=======
std::cout << GridLogMessage << "-- # Trajectory = " << traj << "\n";
>>>>>>> bd181b94819600174e91b331e24f71598b3e7bb9
if (traj < Params.StartTrajectory + Params.NoMetropolisUntil) { if (traj < Params.StartTrajectory + Params.NoMetropolisUntil) {
std::cout << GridLogHMC << "-- Thermalization" << std::endl; std::cout << GridLogHMC << "-- Thermalization" << std::endl;
} }
@ -216,8 +266,6 @@ public:
if (accept) if (accept)
Ucur = Ucopy; Ucur = Ucopy;
double t1=usecond(); double t1=usecond();
std::cout << GridLogHMC << "Total time for trajectory (s): " << (t1-t0)/1e6 << std::endl; std::cout << GridLogHMC << "Total time for trajectory (s): " << (t1-t0)/1e6 << std::endl;

14
Grid/qcd/hmc/integrators/Integrator.h
View File

@ -33,7 +33,6 @@ directory
#define INTEGRATOR_INCLUDED #define INTEGRATOR_INCLUDED
#include <memory> #include <memory>
#include "MomentumFilter.h"
NAMESPACE_BEGIN(Grid); NAMESPACE_BEGIN(Grid);
@ -133,17 +132,22 @@ protected:
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 << "Smearing (on/off): " << as[level].actions.at(a)->is_smeared << std::endl; std::cout << GridLogIntegrator << "Smearing (on/off): " << as[level].actions.at(a)->is_smeared << std::endl;
auto name = as[level].actions.at(a)->action_name();
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
double end_force = usecond(); double end_force = usecond();
Real force_abs = std::sqrt(norm2(force)/U.Grid()->gSites()); //average per-site norm. nb. norm2(latt) = \sum_x norm2(latt[x]) Real force_abs = std::sqrt(norm2(force)/U.Grid()->gSites()); //average per-site norm. nb. norm2(latt) = \sum_x norm2(latt[x])
Real impulse_abs = force_abs * ep * HMC_MOMENTUM_DENOMINATOR; Real impulse_abs = force_abs * ep * HMC_MOMENTUM_DENOMINATOR;
Real max_force_abs = std::sqrt(maxLocalNorm2(force)); Real force_max = std::sqrt(maxLocalNorm2(force));
Real max_impulse_abs = max_force_abs * ep * HMC_MOMENTUM_DENOMINATOR; Real impulse_max = max_force_abs * ep * HMC_MOMENTUM_DENOMINATOR;
std::cout << GridLogIntegrator<< "["<<level<<"]["<<a<<"] Force average: " << force_abs <<" "<<name<<std::endl;
std::cout << GridLogIntegrator<< "["<<level<<"]["<<a<<"] Force max : " << force_max <<" "<<name<<std::endl;
std::cout << GridLogIntegrator<< "["<<level<<"]["<<a<<"] Fdt average : " << force_abs*ep <<" "<<name<<std::endl;
std::cout << GridLogIntegrator<< "["<<level<<"]["<<a<<"] Fdt max : " << force_max*ep <<" "<<name<<std::endl;
std::cout << GridLogIntegrator << "["<<level<<"]["<<a<<"] Force average: " << force_abs << " Max force: " << max_force_abs << " Time step: " << ep << " Impulse average: " << impulse_abs << " Max impulse: " << max_impulse_abs << std::endl;
Mom -= force * ep* HMC_MOMENTUM_DENOMINATOR;; Mom -= force * ep* HMC_MOMENTUM_DENOMINATOR;;
double end_full = usecond(); double end_full = usecond();
double time_full = (end_full - start_full) / 1e3; double time_full = (end_full - start_full) / 1e3;

2
Grid/qcd/utils/CovariantCshift.h
View File

@ -182,7 +182,7 @@ namespace ConjugateBC {
GridBase *grid = Link.Grid(); GridBase *grid = Link.Grid();
int Lmu = grid->GlobalDimensions()[mu] - 1; int Lmu = grid->GlobalDimensions()[mu] - 1;
Lattice<iScalar<vInteger>> coor(grid); Lattice<iScalar<vInteger> > coor(grid);
LatticeCoordinate(coor, mu); LatticeCoordinate(coor, mu);
Lattice<gauge> tmp(grid); Lattice<gauge> tmp(grid);

20
HMC/Mobius2p1fEOFA.cc
View File

@ -188,15 +188,15 @@ int main(int argc, char **argv) {
IntegratorParameters MD; IntegratorParameters MD;
// typedef GenericHMCRunner<LeapFrog> HMCWrapper; // typedef GenericHMCRunner<LeapFrog> HMCWrapper;
// MD.name = std::string("Leap Frog"); // MD.name = std::string("Leap Frog");
typedef GenericHMCRunner<ForceGradient> HMCWrapper; // typedef GenericHMCRunner<ForceGradient> HMCWrapper;
MD.name = std::string("Force Gradient"); // MD.name = std::string("Force Gradient");
// typedef GenericHMCRunner<MinimumNorm2> HMCWrapper; typedef GenericHMCRunner<MinimumNorm2> HMCWrapper;
// MD.name = std::string("MinimumNorm2"); MD.name = std::string("MinimumNorm2");
MD.MDsteps = 6; MD.MDsteps = 12;
MD.trajL = 1.0; MD.trajL = 1.0;
HMCparameters HMCparams; HMCparameters HMCparams;
HMCparams.StartTrajectory = 590; HMCparams.StartTrajectory = 211;
HMCparams.Trajectories = 1000; HMCparams.Trajectories = 1000;
HMCparams.NoMetropolisUntil= 0; HMCparams.NoMetropolisUntil= 0;
// "[HotStart, ColdStart, TepidStart, CheckpointStart]\n"; // "[HotStart, ColdStart, TepidStart, CheckpointStart]\n";
@ -209,9 +209,9 @@ int main(int argc, char **argv) {
TheHMC.Resources.AddFourDimGrid("gauge"); // use default simd lanes decomposition TheHMC.Resources.AddFourDimGrid("gauge"); // use default simd lanes decomposition
CheckpointerParameters CPparams; CheckpointerParameters CPparams;
CPparams.config_prefix = "ckpoint_EODWF_lat"; CPparams.config_prefix = "ckpoint_EOFA_lat";
CPparams.rng_prefix = "ckpoint_EODWF_rng"; CPparams.rng_prefix = "ckpoint_EOFA_rng";
CPparams.saveInterval = 10; CPparams.saveInterval = 1;
CPparams.format = "IEEE64BIG"; CPparams.format = "IEEE64BIG";
TheHMC.Resources.LoadNerscCheckpointer(CPparams); TheHMC.Resources.LoadNerscCheckpointer(CPparams);
@ -263,7 +263,7 @@ int main(int argc, char **argv) {
FermionActionF::ImplParams ParamsF(boundary); FermionActionF::ImplParams ParamsF(boundary);
double ActionStoppingCondition = 1e-10; double ActionStoppingCondition = 1e-10;
double DerivativeStoppingCondition = 1e-6; double DerivativeStoppingCondition = 1e-8;
double MaxCGIterations = 30000; double MaxCGIterations = 30000;
//////////////////////////////////// ////////////////////////////////////

473
HMC/Mobius2p1fEOFA_4dPseudoFermion.cc Normal file
View File

@ -0,0 +1,473 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file:
Copyright (C) 2015-2016
Author: Peter Boyle <pabobyle@ph.ed.ac.uk>
Author: Guido Cossu
Author: David Murphy
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/qcd/action/pseudofermion/TwoFlavourRatioEO4DPseudoFermion.h>
#include <Grid/qcd/action/pseudofermion/TwoFlavourRatio4DPseudoFermion.h>
#ifdef GRID_DEFAULT_PRECISION_DOUBLE
#define MIXED_PRECISION
#endif
NAMESPACE_BEGIN(Grid);
/*
* Need a plan for gauge field update for mixed precision in HMC (2x speed up)
* -- Store the single prec action operator.
* -- Clone the gauge field from the operator function argument.
* -- Build the mixed precision operator dynamically from the passed operator and single prec clone.
*/
template<class FermionOperatorD, class FermionOperatorF, class SchurOperatorD, class SchurOperatorF>
class MixedPrecisionConjugateGradientOperatorFunction : public OperatorFunction<typename FermionOperatorD::FermionField> {
public:
typedef typename FermionOperatorD::FermionField FieldD;
typedef typename FermionOperatorF::FermionField FieldF;
using OperatorFunction<FieldD>::operator();
RealD Tolerance;
RealD InnerTolerance; //Initial tolerance for inner CG. Defaults to Tolerance but can be changed
Integer MaxInnerIterations;
Integer MaxOuterIterations;
GridBase* SinglePrecGrid4; //Grid for single-precision fields
GridBase* SinglePrecGrid5; //Grid for single-precision fields
RealD OuterLoopNormMult; //Stop the outer loop and move to a final double prec solve when the residual is OuterLoopNormMult * Tolerance
FermionOperatorF &FermOpF;
FermionOperatorD &FermOpD;;
SchurOperatorF &LinOpF;
SchurOperatorD &LinOpD;
Integer TotalInnerIterations; //Number of inner CG iterations
Integer TotalOuterIterations; //Number of restarts
Integer TotalFinalStepIterations; //Number of CG iterations in final patch-up step
MixedPrecisionConjugateGradientOperatorFunction(RealD tol,
Integer maxinnerit,
Integer maxouterit,
GridBase* _sp_grid4,
GridBase* _sp_grid5,
FermionOperatorF &_FermOpF,
FermionOperatorD &_FermOpD,
SchurOperatorF &_LinOpF,
SchurOperatorD &_LinOpD):
LinOpF(_LinOpF),
LinOpD(_LinOpD),
FermOpF(_FermOpF),
FermOpD(_FermOpD),
Tolerance(tol),
InnerTolerance(tol),
MaxInnerIterations(maxinnerit),
MaxOuterIterations(maxouterit),
SinglePrecGrid4(_sp_grid4),
SinglePrecGrid5(_sp_grid5),
OuterLoopNormMult(100.)
{
/* Debugging instances of objects; references are stored
std::cout << GridLogMessage << " Mixed precision CG wrapper LinOpF " <<std::hex<< &LinOpF<<std::dec <<std::endl;
std::cout << GridLogMessage << " Mixed precision CG wrapper LinOpD " <<std::hex<< &LinOpD<<std::dec <<std::endl;
std::cout << GridLogMessage << " Mixed precision CG wrapper FermOpF " <<std::hex<< &FermOpF<<std::dec <<std::endl;
std::cout << GridLogMessage << " Mixed precision CG wrapper FermOpD " <<std::hex<< &FermOpD<<std::dec <<std::endl;
*/
};
void operator()(LinearOperatorBase<FieldD> &LinOpU, const FieldD &src, FieldD &psi) {
std::cout << GridLogMessage << " Mixed precision CG wrapper operator() "<<std::endl;
SchurOperatorD * SchurOpU = static_cast<SchurOperatorD *>(&LinOpU);
// std::cout << GridLogMessage << " Mixed precision CG wrapper operator() FermOpU " <<std::hex<< &(SchurOpU->_Mat)<<std::dec <<std::endl;
// std::cout << GridLogMessage << " Mixed precision CG wrapper operator() FermOpD " <<std::hex<< &(LinOpD._Mat) <<std::dec <<std::endl;
// Assumption made in code to extract gauge field
// We could avoid storing LinopD reference alltogether ?
assert(&(SchurOpU->_Mat)==&(LinOpD._Mat));
////////////////////////////////////////////////////////////////////////////////////
// Must snarf a single precision copy of the gauge field in Linop_d argument
////////////////////////////////////////////////////////////////////////////////////
typedef typename FermionOperatorF::GaugeField GaugeFieldF;
typedef typename FermionOperatorF::GaugeLinkField GaugeLinkFieldF;
typedef typename FermionOperatorD::GaugeField GaugeFieldD;
typedef typename FermionOperatorD::GaugeLinkField GaugeLinkFieldD;
GridBase * GridPtrF = SinglePrecGrid4;
GridBase * GridPtrD = FermOpD.Umu.Grid();
GaugeFieldF U_f (GridPtrF);
GaugeLinkFieldF Umu_f(GridPtrF);
// std::cout << " Dim gauge field "<<GridPtrF->Nd()<<std::endl; // 4d
// std::cout << " Dim gauge field "<<GridPtrD->Nd()<<std::endl; // 4d
////////////////////////////////////////////////////////////////////////////////////
// Moving this to a Clone method of fermion operator would allow to duplicate the
// physics parameters and decrease gauge field copies
////////////////////////////////////////////////////////////////////////////////////
GaugeLinkFieldD Umu_d(GridPtrD);
for(int mu=0;mu<Nd*2;mu++){
Umu_d = PeekIndex<LorentzIndex>(FermOpD.Umu, mu);
precisionChange(Umu_f,Umu_d);
PokeIndex<LorentzIndex>(FermOpF.Umu, Umu_f, mu);
}
pickCheckerboard(Even,FermOpF.UmuEven,FermOpF.Umu);
pickCheckerboard(Odd ,FermOpF.UmuOdd ,FermOpF.Umu);
////////////////////////////////////////////////////////////////////////////////////
// Could test to make sure that LinOpF and LinOpD agree to single prec?
////////////////////////////////////////////////////////////////////////////////////
/*
GridBase *Fgrid = psi._grid;
FieldD tmp2(Fgrid);
FieldD tmp1(Fgrid);
LinOpU.Op(src,tmp1);
LinOpD.Op(src,tmp2);
std::cout << " Double gauge field "<< norm2(FermOpD.Umu)<<std::endl;
std::cout << " Single gauge field "<< norm2(FermOpF.Umu)<<std::endl;
std::cout << " Test of operators "<<norm2(tmp1)<<std::endl;
std::cout << " Test of operators "<<norm2(tmp2)<<std::endl;
tmp1=tmp1-tmp2;
std::cout << " Test of operators diff "<<norm2(tmp1)<<std::endl;
*/
////////////////////////////////////////////////////////////////////////////////////
// Make a mixed precision conjugate gradient
////////////////////////////////////////////////////////////////////////////////////
MixedPrecisionConjugateGradient<FieldD,FieldF> MPCG(Tolerance,MaxInnerIterations,MaxOuterIterations,SinglePrecGrid5,LinOpF,LinOpD);
std::cout << GridLogMessage << "Calling mixed precision Conjugate Gradient" <<std::endl;
MPCG(src,psi);
}
};
NAMESPACE_END(Grid);
int main(int argc, char **argv) {
using namespace Grid;
Grid_init(&argc, &argv);
int threads = GridThread::GetThreads();
// here make a routine to print all the relevant information on the run
std::cout << GridLogMessage << "Grid is setup to use " << threads << " threads" << std::endl;
// Typedefs to simplify notation
typedef WilsonImplR FermionImplPolicy;
typedef MobiusFermionR FermionAction;
typedef MobiusFermionF FermionActionF;
typedef MobiusEOFAFermionR FermionEOFAAction;
typedef MobiusEOFAFermionF FermionEOFAActionF;
typedef typename FermionAction::FermionField FermionField;
typedef typename FermionActionF::FermionField FermionFieldF;
typedef Grid::XmlReader Serialiser;
//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
IntegratorParameters MD;
// typedef GenericHMCRunner<LeapFrog> HMCWrapper;
// MD.name = std::string("Leap Frog");
// typedef GenericHMCRunner<ForceGradient> HMCWrapper;
// MD.name = std::string("Force Gradient");
typedef GenericHMCRunner<MinimumNorm2> HMCWrapper;
MD.name = std::string("MinimumNorm2");
MD.MDsteps = 12;
MD.trajL = 1.0;
HMCparameters HMCparams;
HMCparams.StartTrajectory = 211;
HMCparams.Trajectories = 1000;
HMCparams.NoMetropolisUntil= 0;
// "[HotStart, ColdStart, TepidStart, CheckpointStart]\n";
// HMCparams.StartingType =std::string("ColdStart");
HMCparams.StartingType =std::string("CheckpointStart");
HMCparams.MD = MD;
HMCWrapper TheHMC(HMCparams);
// Grid from the command line arguments --grid and --mpi
TheHMC.Resources.AddFourDimGrid("gauge"); // use default simd lanes decomposition
CheckpointerParameters CPparams;
CPparams.config_prefix = "ckpoint_EOFA4D_lat";
CPparams.rng_prefix = "ckpoint_EOFA4D_rng";
CPparams.saveInterval = 1;
CPparams.format = "IEEE64BIG";
TheHMC.Resources.LoadNerscCheckpointer(CPparams);
RNGModuleParameters RNGpar;
RNGpar.serial_seeds = "1 2 3 4 5";
RNGpar.parallel_seeds = "6 7 8 9 10";
TheHMC.Resources.SetRNGSeeds(RNGpar);
// Construct observables
// here there is too much indirection
typedef PlaquetteMod<HMCWrapper::ImplPolicy> PlaqObs;
TheHMC.Resources.AddObservable<PlaqObs>();
//////////////////////////////////////////////
const int Ls = 16;
Real beta = 2.13;
Real light_mass = 0.01;
Real strange_mass = 0.04;
Real pv_mass = 1.0;
RealD M5 = 1.8;
RealD b = 1.0;
RealD c = 0.0;
std::vector<Real> hasenbusch({ 0.1, 0.3, 0.6 });
auto GridPtr = TheHMC.Resources.GetCartesian();
auto GridRBPtr = TheHMC.Resources.GetRBCartesian();
auto FGrid = SpaceTimeGrid::makeFiveDimGrid(Ls,GridPtr);
auto FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,GridPtr);
Coordinate latt = GridDefaultLatt();
Coordinate mpi = GridDefaultMpi();
Coordinate simdF = GridDefaultSimd(Nd,vComplexF::Nsimd());
Coordinate simdD = GridDefaultSimd(Nd,vComplexD::Nsimd());
auto GridPtrF = SpaceTimeGrid::makeFourDimGrid(latt,simdF,mpi);
auto GridRBPtrF = SpaceTimeGrid::makeFourDimRedBlackGrid(GridPtrF);
auto FGridF = SpaceTimeGrid::makeFiveDimGrid(Ls,GridPtrF);
auto FrbGridF = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,GridPtrF);
IwasakiGaugeActionR GaugeAction(beta);
// temporarily need a gauge field
LatticeGaugeField U(GridPtr);
LatticeGaugeFieldF UF(GridPtrF);
// These lines are unecessary if BC are all periodic
std::vector<Complex> boundary = {1,1,1,-1};
FermionAction::ImplParams Params(boundary);
FermionActionF::ImplParams ParamsF(boundary);
double ActionStoppingCondition = 1e-10;
double DerivativeStoppingCondition = 1e-8;
double MaxCGIterations = 30000;
////////////////////////////////////
// Collect actions
////////////////////////////////////
ActionLevel<HMCWrapper::Field> Level1(1);
ActionLevel<HMCWrapper::Field> Level2(8);
////////////////////////////////////
// Strange action
////////////////////////////////////
typedef SchurDiagMooeeOperator<FermionActionF,FermionFieldF> LinearOperatorF;
typedef SchurDiagMooeeOperator<FermionAction ,FermionField > LinearOperatorD;
typedef SchurDiagMooeeDagOperator<FermionActionF,FermionFieldF> LinearOperatorDagF;
typedef SchurDiagMooeeDagOperator<FermionAction ,FermionField > LinearOperatorDagD;
typedef SchurDiagMooeeOperator<FermionEOFAActionF,FermionFieldF> LinearOperatorEOFAF;
typedef SchurDiagMooeeOperator<FermionEOFAAction ,FermionField > LinearOperatorEOFAD;
typedef MixedPrecisionConjugateGradientOperatorFunction<MobiusFermionD,MobiusFermionF,LinearOperatorD,LinearOperatorF> MxPCG;
typedef MixedPrecisionConjugateGradientOperatorFunction<MobiusFermionD,MobiusFermionF,LinearOperatorDagD,LinearOperatorDagF> MxDagPCG;
typedef MixedPrecisionConjugateGradientOperatorFunction<MobiusEOFAFermionD,MobiusEOFAFermionF,LinearOperatorEOFAD,LinearOperatorEOFAF> MxPCG_EOFA;
// DJM: setup for EOFA ratio (Mobius)
OneFlavourRationalParams OFRp;
OFRp.lo = 0.1;
OFRp.hi = 25.0;
OFRp.MaxIter = 10000;
OFRp.tolerance= 1.0e-9;
OFRp.degree = 14;
OFRp.precision= 50;
MobiusEOFAFermionR Strange_Op_L (U , *FGrid , *FrbGrid , *GridPtr , *GridRBPtr , strange_mass, strange_mass, pv_mass, 0.0, -1, M5, b, c);
MobiusEOFAFermionF Strange_Op_LF(UF, *FGridF, *FrbGridF, *GridPtrF, *GridRBPtrF, strange_mass, strange_mass, pv_mass, 0.0, -1, M5, b, c);
MobiusEOFAFermionR Strange_Op_R (U , *FGrid , *FrbGrid , *GridPtr , *GridRBPtr , pv_mass, strange_mass, pv_mass, -1.0, 1, M5, b, c);
MobiusEOFAFermionF Strange_Op_RF(UF, *FGridF, *FrbGridF, *GridPtrF, *GridRBPtrF, pv_mass, strange_mass, pv_mass, -1.0, 1, M5, b, c);
ConjugateGradient<FermionField> ActionCG(ActionStoppingCondition,MaxCGIterations);
ConjugateGradient<FermionField> DerivativeCG(DerivativeStoppingCondition,MaxCGIterations);
#ifdef MIXED_PRECISION
const int MX_inner = 1000;
// Mixed precision EOFA
LinearOperatorEOFAD Strange_LinOp_L (Strange_Op_L);
LinearOperatorEOFAD Strange_LinOp_R (Strange_Op_R);
LinearOperatorEOFAF Strange_LinOp_LF(Strange_Op_LF);
LinearOperatorEOFAF Strange_LinOp_RF(Strange_Op_RF);
MxPCG_EOFA ActionCGL(ActionStoppingCondition,
MX_inner,
MaxCGIterations,
GridPtrF,
FrbGridF,
Strange_Op_LF,Strange_Op_L,
Strange_LinOp_LF,Strange_LinOp_L);
MxPCG_EOFA DerivativeCGL(DerivativeStoppingCondition,
MX_inner,
MaxCGIterations,
GridPtrF,
FrbGridF,
Strange_Op_LF,Strange_Op_L,
Strange_LinOp_LF,Strange_LinOp_L);
MxPCG_EOFA ActionCGR(ActionStoppingCondition,
MX_inner,
MaxCGIterations,
GridPtrF,
FrbGridF,
Strange_Op_RF,Strange_Op_R,
Strange_LinOp_RF,Strange_LinOp_R);
MxPCG_EOFA DerivativeCGR(DerivativeStoppingCondition,
MX_inner,
MaxCGIterations,
GridPtrF,
FrbGridF,
Strange_Op_RF,Strange_Op_R,
Strange_LinOp_RF,Strange_LinOp_R);
ExactOneFlavourRatioPseudoFermionAction<FermionImplPolicy>
EOFA(Strange_Op_L, Strange_Op_R,
ActionCG,
ActionCGL, ActionCGR,
DerivativeCGL, DerivativeCGR,
OFRp, true);
#else
ExactOneFlavourRatioPseudoFermionAction<FermionImplPolicy>
EOFA(Strange_Op_L, Strange_Op_R,
ActionCG,
ActionCG, ActionCG,
DerivativeCG, DerivativeCG,
OFRp, true);
#endif
Level1.push_back(&EOFA);
////////////////////////////////////
// up down action
////////////////////////////////////
std::vector<Real> light_den;
std::vector<Real> light_num;
int n_hasenbusch = hasenbusch.size();
light_den.push_back(light_mass);
for(int h=0;h<n_hasenbusch;h++){
light_den.push_back(hasenbusch[h]);
light_num.push_back(hasenbusch[h]);
}
light_num.push_back(pv_mass);
//////////////////////////////////////////////////////////////
// Forced to replicate the MxPCG and DenominatorsF etc.. because
// there is no convenient way to "Clone" physics params from double op
// into single op for any operator pair.
// Same issue prevents using MxPCG in the Heatbath step
//////////////////////////////////////////////////////////////
std::vector<FermionAction *> Numerators;
std::vector<FermionAction *> Denominators;
std::vector<TwoFlavourRatioEO4DPseudoFermionAction<FermionImplPolicy> *> Quotients;
std::vector<MxPCG *> ActionMPCG;
std::vector<MxPCG *> MPCG;
std::vector<MxDagPCG *> MPCGdag;
std::vector<FermionActionF *> DenominatorsF;
std::vector<LinearOperatorD *> LinOpD;
std::vector<LinearOperatorF *> LinOpF;
std::vector<LinearOperatorDagD *> LinOpDagD;
std::vector<LinearOperatorDagF *> LinOpDagF;
for(int h=0;h<n_hasenbusch+1;h++){
std::cout << GridLogMessage << " 2f quotient Action "<< light_num[h] << " / " << light_den[h]<< std::endl;
Numerators.push_back (new FermionAction(U,*FGrid,*FrbGrid,*GridPtr,*GridRBPtr,light_num[h],M5,b,c, Params));
Denominators.push_back(new FermionAction(U,*FGrid,*FrbGrid,*GridPtr,*GridRBPtr,light_den[h],M5,b,c, Params));
#ifdef MIXED_PRECISION
////////////////////////////////////////////////////////////////////////////
// Mixed precision CG for 2f force
////////////////////////////////////////////////////////////////////////////
DenominatorsF.push_back(new FermionActionF(UF,*FGridF,*FrbGridF,*GridPtrF,*GridRBPtrF,light_den[h],M5,b,c, ParamsF));
LinOpD.push_back(new LinearOperatorD(*Denominators[h]));
LinOpF.push_back(new LinearOperatorF(*DenominatorsF[h]));
LinOpDagD.push_back(new LinearOperatorDagD(*Denominators[h]));
LinOpDagF.push_back(new LinearOperatorDagF(*DenominatorsF[h]));
MPCG.push_back(new MxPCG(DerivativeStoppingCondition,
MX_inner,
MaxCGIterations,
GridPtrF,
FrbGridF,
*DenominatorsF[h],*Denominators[h],
*LinOpF[h], *LinOpD[h]) );
MPCGdag.push_back(new MxDagPCG(DerivativeStoppingCondition,
MX_inner,
MaxCGIterations,
GridPtrF,
FrbGridF,
*DenominatorsF[h],*Denominators[h],
*LinOpDagF[h], *LinOpDagD[h]) );
ActionMPCG.push_back(new MxPCG(ActionStoppingCondition,
MX_inner,
MaxCGIterations,
GridPtrF,
FrbGridF,
*DenominatorsF[h],*Denominators[h],
*LinOpF[h], *LinOpD[h]) );
// Heatbath not mixed yet. As inverts numerators not so important as raised mass.
Quotients.push_back (new TwoFlavourRatioEO4DPseudoFermionAction<FermionImplPolicy>(*Numerators[h],*Denominators[h],*MPCG[h],*MPCGdag[h],*ActionMPCG[h],ActionCG));
#else
////////////////////////////////////////////////////////////////////////////
// Standard CG for 2f force
////////////////////////////////////////////////////////////////////////////
Quotients.push_back (new TwoFlavourRatioEO4DPseudoFermionAction<FermionImplPolicy>(*Numerators[h],*Denominators[h],DerivativeCG,ActionCG));
#endif
}
for(int h=0;h<n_hasenbusch+1;h++){
Level1.push_back(Quotients[h]);
}
/////////////////////////////////////////////////////////////
// Gauge action
/////////////////////////////////////////////////////////////
Level2.push_back(&GaugeAction);
TheHMC.TheAction.push_back(Level1);
TheHMC.TheAction.push_back(Level2);
std::cout << GridLogMessage << " Action complete "<< std::endl;
/////////////////////////////////////////////////////////////
// HMC parameters are serialisable
std::cout << GridLogMessage << " Running the HMC "<< std::endl;
TheHMC.Run(); // no smearing
Grid_finalize();
} // main

12
HMC/Mobius2p1fRHMC.cc
View File

@ -52,16 +52,16 @@ int main(int argc, char **argv) {
// MD.name = std::string("Force Gradient"); // MD.name = std::string("Force Gradient");
typedef GenericHMCRunner<MinimumNorm2> HMCWrapper; typedef GenericHMCRunner<MinimumNorm2> HMCWrapper;
MD.name = std::string("MinimumNorm2"); MD.name = std::string("MinimumNorm2");
MD.MDsteps = 20; MD.MDsteps = 12;
MD.trajL = 1.0; MD.trajL = 1.0;
HMCparameters HMCparams; HMCparameters HMCparams;
HMCparams.StartTrajectory = 0; HMCparams.StartTrajectory = 139;
HMCparams.Trajectories = 200; HMCparams.Trajectories = 200;
HMCparams.NoMetropolisUntil= 0; HMCparams.NoMetropolisUntil= 0;
// "[HotStart, ColdStart, TepidStart, CheckpointStart]\n"; // "[HotStart, ColdStart, TepidStart, CheckpointStart]\n";
HMCparams.StartingType =std::string("ColdStart"); // HMCparams.StartingType =std::string("ColdStart");
// HMCparams.StartingType =std::string("CheckpointStart"); HMCparams.StartingType =std::string("CheckpointStart");
HMCparams.MD = MD; HMCparams.MD = MD;
HMCWrapper TheHMC(HMCparams); HMCWrapper TheHMC(HMCparams);
@ -71,7 +71,7 @@ int main(int argc, char **argv) {
CheckpointerParameters CPparams; CheckpointerParameters CPparams;
CPparams.config_prefix = "ckpoint_EODWF_lat"; CPparams.config_prefix = "ckpoint_EODWF_lat";
CPparams.rng_prefix = "ckpoint_EODWF_rng"; CPparams.rng_prefix = "ckpoint_EODWF_rng";
CPparams.saveInterval = 10; CPparams.saveInterval = 1;
CPparams.format = "IEEE64BIG"; CPparams.format = "IEEE64BIG";
TheHMC.Resources.LoadNerscCheckpointer(CPparams); TheHMC.Resources.LoadNerscCheckpointer(CPparams);
@ -130,7 +130,7 @@ int main(int argc, char **argv) {
// Collect actions // Collect actions
//////////////////////////////////// ////////////////////////////////////
ActionLevel<HMCWrapper::Field> Level1(1); ActionLevel<HMCWrapper::Field> Level1(1);
ActionLevel<HMCWrapper::Field> Level2(4); ActionLevel<HMCWrapper::Field> Level2(8);
//////////////////////////////////// ////////////////////////////////////
// Strange action // Strange action

197
HMC/Mobius2p1fRHMC_4dPseudoFermion.cc Normal file
View File

@ -0,0 +1,197 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: ./tests/Test_hmc_EODWFRatio.cc
Copyright (C) 2015-2016
Author: Peter Boyle <pabobyle@ph.ed.ac.uk>
Author: Guido Cossu <guido.cossu@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>
int main(int argc, char **argv) {
using namespace Grid;
Grid_init(&argc, &argv);
int threads = GridThread::GetThreads();
// here make a routine to print all the relevant information on the run
std::cout << GridLogMessage << "Grid is setup to use " << threads << " threads" << std::endl;
// Typedefs to simplify notation
typedef WilsonImplR FermionImplPolicy;
typedef MobiusFermionR FermionAction;
typedef typename FermionAction::FermionField FermionField;
typedef Grid::XmlReader Serialiser;
//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
IntegratorParameters MD;
// typedef GenericHMCRunner<LeapFrog> HMCWrapper;
// MD.name = std::string("Leap Frog");
// typedef GenericHMCRunner<ForceGradient> HMCWrapper;
// MD.name = std::string("Force Gradient");
typedef GenericHMCRunner<MinimumNorm2> HMCWrapper;
MD.name = std::string("MinimumNorm2");
MD.MDsteps = 10;
MD.trajL = 1.0;
HMCparameters HMCparams;
HMCparams.StartTrajectory = 137;
HMCparams.Trajectories = 200;
HMCparams.NoMetropolisUntil= 0;
// "[HotStart, ColdStart, TepidStart, CheckpointStart]\n";
// HMCparams.StartingType =std::string("ColdStart");
HMCparams.StartingType =std::string("CheckpointStart");
HMCparams.MD = MD;
HMCWrapper TheHMC(HMCparams);
// Grid from the command line arguments --grid and --mpi
TheHMC.Resources.AddFourDimGrid("gauge"); // use default simd lanes decomposition
CheckpointerParameters CPparams;
CPparams.config_prefix = "ckpoint_EODWF_lat";
CPparams.rng_prefix = "ckpoint_EODWF_rng";
CPparams.saveInterval = 1;
CPparams.format = "IEEE64BIG";
TheHMC.Resources.LoadNerscCheckpointer(CPparams);
RNGModuleParameters RNGpar;
RNGpar.serial_seeds = "1 2 3 4 5";
RNGpar.parallel_seeds = "6 7 8 9 10";
TheHMC.Resources.SetRNGSeeds(RNGpar);
// Construct observables
// here there is too much indirection
typedef PlaquetteMod<HMCWrapper::ImplPolicy> PlaqObs;
TheHMC.Resources.AddObservable<PlaqObs>();
//////////////////////////////////////////////
const int Ls = 16;
Real beta = 2.13;
Real light_mass = 0.01;
Real strange_mass = 0.04;
Real pv_mass = 1.0;
RealD M5 = 1.8;
RealD b = 1.0;
RealD c = 0.0;
// FIXME:
// Same in MC and MD
// Need to mix precision too
OneFlavourRationalParams OFRp;
OFRp.lo = 4.0e-3;
OFRp.hi = 30.0;
OFRp.MaxIter = 10000;
OFRp.tolerance= 1.0e-10;
OFRp.degree = 16;
OFRp.precision= 50;
std::vector<Real> hasenbusch({ 0.1 });
auto GridPtr = TheHMC.Resources.GetCartesian();
auto GridRBPtr = TheHMC.Resources.GetRBCartesian();
auto FGrid = SpaceTimeGrid::makeFiveDimGrid(Ls,GridPtr);
auto FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,GridPtr);
IwasakiGaugeActionR GaugeAction(beta);
// temporarily need a gauge field
LatticeGaugeField U(GridPtr);
// These lines are unecessary if BC are all periodic
std::vector<Complex> boundary = {1,1,1,-1};
FermionAction::ImplParams Params(boundary);
double StoppingCondition = 1e-10;
double MaxCGIterations = 30000;
ConjugateGradient<FermionField> CG(StoppingCondition,MaxCGIterations);
////////////////////////////////////
// Collect actions
////////////////////////////////////
ActionLevel<HMCWrapper::Field> Level1(1);
ActionLevel<HMCWrapper::Field> Level2(8);
////////////////////////////////////
// Strange action
////////////////////////////////////
// FermionAction StrangeOp(U,*FGrid,*FrbGrid,*GridPtr,*GridRBPtr,light_mass,M5,b,c, Params);
// DomainWallEOFAFermionR Strange_Op_L(Umu, *FGrid, *FrbGrid, *UGrid, *UrbGrid, mf, mf, mb, shift_L, pm, M5);
// DomainWallEOFAFermionR Strange_Op_R(Umu, *FGrid, *FrbGrid, *UGrid, *UrbGrid, mb, mf, mb, shift_R, pm, M5);
// ExactOneFlavourRatioPseudoFermionAction EOFA(Strange_Op_L,Strange_Op_R,CG,ofp, false);
FermionAction StrangeOp (U,*FGrid,*FrbGrid,*GridPtr,*GridRBPtr,strange_mass,M5,b,c, Params);
FermionAction StrangePauliVillarsOp(U,*FGrid,*FrbGrid,*GridPtr,*GridRBPtr,pv_mass, M5,b,c, Params);
OneFlavourEvenOddRatioRationalPseudoFermionAction<FermionImplPolicy> StrangePseudoFermion(StrangePauliVillarsOp,StrangeOp,OFRp);
Level1.push_back(&StrangePseudoFermion);
////////////////////////////////////
// up down action
////////////////////////////////////
std::vector<Real> light_den;
std::vector<Real> light_num;
int n_hasenbusch = hasenbusch.size();
light_den.push_back(light_mass);
for(int h=0;h<n_hasenbusch;h++){
light_den.push_back(hasenbusch[h]);
light_num.push_back(hasenbusch[h]);
}
light_num.push_back(pv_mass);
std::vector<FermionAction *> Numerators;
std::vector<FermionAction *> Denominators;
std::vector<TwoFlavourRatioPseudoFermionAction<FermionImplPolicy> *> Quotients;
for(int h=0;h<n_hasenbusch+1;h++){
std::cout << GridLogMessage << " 2f quotient Action "<< light_num[h] << " / " << light_den[h]<< std::endl;
Numerators.push_back (new FermionAction(U,*FGrid,*FrbGrid,*GridPtr,*GridRBPtr,light_num[h],M5,b,c, Params));
Denominators.push_back(new FermionAction(U,*FGrid,*FrbGrid,*GridPtr,*GridRBPtr,light_den[h],M5,b,c, Params));
Quotients.push_back (new TwoFlavourRatioPseudoFermionAction<FermionImplPolicy>(*Numerators[h],*Denominators[h],CG,CG));
}
for(int h=0;h<n_hasenbusch+1;h++){
Level1.push_back(Quotients[h]);
}
/////////////////////////////////////////////////////////////
// Gauge action
/////////////////////////////////////////////////////////////
Level2.push_back(&GaugeAction);
TheHMC.TheAction.push_back(Level1);
TheHMC.TheAction.push_back(Level2);
std::cout << GridLogMessage << " Action complete "<< std::endl;
/////////////////////////////////////////////////////////////
// HMC parameters are serialisable
std::cout << GridLogMessage << " Running the HMC "<< std::endl;
TheHMC.Run(); // no smearing
Grid_finalize();
} // main

198
HMC/Mobius2p1fRHMC_4dPseudoFermionSchurSolver.cc Normal file
View File

@ -0,0 +1,198 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: ./tests/Test_hmc_EODWFRatio.cc
Copyright (C) 2015-2016
Author: Peter Boyle <pabobyle@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/qcd/action/pseudofermion/TwoFlavourRatioEO4DPseudoFermion.h>
#include <Grid/qcd/action/pseudofermion/TwoFlavourRatio4DPseudoFermion.h>
int main(int argc, char **argv) {
using namespace Grid;
Grid_init(&argc, &argv);
int threads = GridThread::GetThreads();
// here make a routine to print all the relevant information on the run
std::cout << GridLogMessage << "Grid is setup to use " << threads << " threads" << std::endl;
// Typedefs to simplify notation
typedef WilsonImplR FermionImplPolicy;
typedef MobiusFermionR FermionAction;
typedef typename FermionAction::FermionField FermionField;
typedef Grid::XmlReader Serialiser;
//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
IntegratorParameters MD;
// typedef GenericHMCRunner<LeapFrog> HMCWrapper;
// MD.name = std::string("Leap Frog");
// typedef GenericHMCRunner<ForceGradient> HMCWrapper;
// MD.name = std::string("Force Gradient");
typedef GenericHMCRunner<MinimumNorm2> HMCWrapper;
MD.name = std::string("MinimumNorm2");
MD.MDsteps = 12;
MD.trajL = 1.0;
HMCparameters HMCparams;
HMCparams.StartTrajectory = 211;
HMCparams.Trajectories = 200;
HMCparams.NoMetropolisUntil= 0;
// "[HotStart, ColdStart, TepidStart, CheckpointStart]\n";
// HMCparams.StartingType =std::string("ColdStart");
HMCparams.StartingType =std::string("CheckpointStart");
HMCparams.MD = MD;
HMCWrapper TheHMC(HMCparams);
// Grid from the command line arguments --grid and --mpi
TheHMC.Resources.AddFourDimGrid("gauge"); // use default simd lanes decomposition
CheckpointerParameters CPparams;
CPparams.config_prefix = "ckpoint_4dDWF_lat";
CPparams.rng_prefix = "ckpoint_4dDWF_rng";
CPparams.saveInterval = 1;
CPparams.format = "IEEE64BIG";
TheHMC.Resources.LoadNerscCheckpointer(CPparams);
RNGModuleParameters RNGpar;
RNGpar.serial_seeds = "1 2 3 4 5";
RNGpar.parallel_seeds = "6 7 8 9 10";
TheHMC.Resources.SetRNGSeeds(RNGpar);
// Construct observables
// here there is too much indirection
typedef PlaquetteMod<HMCWrapper::ImplPolicy> PlaqObs;
TheHMC.Resources.AddObservable<PlaqObs>();
//////////////////////////////////////////////
const int Ls = 16;
Real beta = 2.13;
Real light_mass = 0.01;
Real strange_mass = 0.04;
Real pv_mass = 1.0;
RealD M5 = 1.8;
RealD b = 1.0;
RealD c = 0.0;
// FIXME:
// Same in MC and MD
// Need to mix precision too
OneFlavourRationalParams OFRp;
OFRp.lo = 4.0e-3;
OFRp.hi = 30.0;
OFRp.MaxIter = 10000;
OFRp.tolerance= 1.0e-10;
OFRp.degree = 16;
OFRp.precision= 50;
std::vector<Real> hasenbusch({ 0.1 });
auto GridPtr = TheHMC.Resources.GetCartesian();
auto GridRBPtr = TheHMC.Resources.GetRBCartesian();
auto FGrid = SpaceTimeGrid::makeFiveDimGrid(Ls,GridPtr);
auto FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,GridPtr);
IwasakiGaugeActionR GaugeAction(beta);
// temporarily need a gauge field
LatticeGaugeField U(GridPtr);
// These lines are unecessary if BC are all periodic
std::vector<Complex> boundary = {1,1,1,-1};
FermionAction::ImplParams Params(boundary);
double StoppingCondition = 1e-10;
double MaxCGIterations = 30000;
ConjugateGradient<FermionField> CG(StoppingCondition,MaxCGIterations);
////////////////////////////////////
// Collect actions
////////////////////////////////////
ActionLevel<HMCWrapper::Field> Level1(1);
ActionLevel<HMCWrapper::Field> Level2(8);
////////////////////////////////////
// Strange action
////////////////////////////////////
// FermionAction StrangeOp(U,*FGrid,*FrbGrid,*GridPtr,*GridRBPtr,light_mass,M5,b,c, Params);
// DomainWallEOFAFermionR Strange_Op_L(Umu, *FGrid, *FrbGrid, *UGrid, *UrbGrid, mf, mf, mb, shift_L, pm, M5);
// DomainWallEOFAFermionR Strange_Op_R(Umu, *FGrid, *FrbGrid, *UGrid, *UrbGrid, mb, mf, mb, shift_R, pm, M5);
// ExactOneFlavourRatioPseudoFermionAction EOFA(Strange_Op_L,Strange_Op_R,CG,ofp, false);
FermionAction StrangeOp (U,*FGrid,*FrbGrid,*GridPtr,*GridRBPtr,strange_mass,M5,b,c, Params);
FermionAction StrangePauliVillarsOp(U,*FGrid,*FrbGrid,*GridPtr,*GridRBPtr,pv_mass, M5,b,c, Params);
OneFlavourEvenOddRatioRationalPseudoFermionAction<FermionImplPolicy> StrangePseudoFermion(StrangePauliVillarsOp,StrangeOp,OFRp);
Level1.push_back(&StrangePseudoFermion);
////////////////////////////////////
// up down action
////////////////////////////////////
std::vector<Real> light_den;
std::vector<Real> light_num;
int n_hasenbusch = hasenbusch.size();
light_den.push_back(light_mass);
for(int h=0;h<n_hasenbusch;h++){
light_den.push_back(hasenbusch[h]);
light_num.push_back(hasenbusch[h]);
}
light_num.push_back(pv_mass);
std::vector<FermionAction *> Numerators;
std::vector<FermionAction *> Denominators;
std::vector<TwoFlavourRatioEO4DPseudoFermionAction<FermionImplPolicy> *> Quotients;
for(int h=0;h<n_hasenbusch+1;h++){
std::cout << GridLogMessage << " 2f quotient Action "<< light_num[h] << " / " << light_den[h]<< std::endl;
Numerators.push_back (new FermionAction(U,*FGrid,*FrbGrid,*GridPtr,*GridRBPtr,light_num[h],M5,b,c, Params));
Denominators.push_back(new FermionAction(U,*FGrid,*FrbGrid,*GridPtr,*GridRBPtr,light_den[h],M5,b,c, Params));
Quotients.push_back (new TwoFlavourRatioEO4DPseudoFermionAction<FermionImplPolicy>(*Numerators[h],*Denominators[h],CG,CG));
}
for(int h=0;h<n_hasenbusch+1;h++){
Level1.push_back(Quotients[h]);
}
/////////////////////////////////////////////////////////////
// Gauge action
/////////////////////////////////////////////////////////////
Level2.push_back(&GaugeAction);
TheHMC.TheAction.push_back(Level1);
TheHMC.TheAction.push_back(Level2);
std::cout << GridLogMessage << " Action complete "<< std::endl;
/////////////////////////////////////////////////////////////
// HMC parameters are serialisable
std::cout << GridLogMessage << " Running the HMC "<< std::endl;
TheHMC.Run(); // no smearing
Grid_finalize();
} // main

10
TODO
View File

@ -1,5 +1,11 @@
-- comms threads issue?? --
-- Part done: Staggered kernel performance on GPU -- Comms threads issue??
-- Part done: Staggered kernel performance on GPU ; eliminate replicas
-- Antonin - Nd, Nc generic hide and make Gimpl
-- DWF 5d RB case / Shamir
-- 4D pseudofermion options
-- DDHMC
--
========================================================= =========================================================
General General

21
tests/debug/Test_cayley_even_odd.cc
View File

@ -208,6 +208,27 @@ void TestWhat(What & Ddwf,
std::cout<<GridLogMessage <<"pDce - conj(cDpo) "<< pDce-conj(cDpo) <<std::endl; std::cout<<GridLogMessage <<"pDce - conj(cDpo) "<< pDce-conj(cDpo) <<std::endl;
std::cout<<GridLogMessage <<"pDco - conj(cDpe) "<< pDco-conj(cDpe) <<std::endl; std::cout<<GridLogMessage <<"pDco - conj(cDpe) "<< pDco-conj(cDpe) <<std::endl;
std::cout<<GridLogMessage<<"=============================================================="<<std::endl;
std::cout<<GridLogMessage<<"= Test MooeeDag is the dagger of Mooee by requiring "<<std::endl;
std::cout<<GridLogMessage<<"= < phi | Deo | chi > * = < chi | Deo^dag| phi> "<<std::endl;
std::cout<<GridLogMessage<<"=============================================================="<<std::endl;
Ddwf.Mooee(chi_e,dchi_e);
Ddwf.Mooee(chi_o,dchi_o);
Ddwf.MooeeDag(phi_e,dphi_e);
Ddwf.MooeeDag(phi_o,dphi_o);
pDce = innerProduct(phi_e,dchi_e);
pDco = innerProduct(phi_o,dchi_o);
cDpe = innerProduct(chi_e,dphi_e);
cDpo = innerProduct(chi_o,dphi_o);
std::cout<<GridLogMessage <<"e "<<pDce<<" "<<cDpe <<std::endl;
std::cout<<GridLogMessage <<"o "<<pDco<<" "<<cDpo <<std::endl;
std::cout<<GridLogMessage <<"pDce - conj(cDpo) "<< pDce-conj(cDpe) <<std::endl;
std::cout<<GridLogMessage <<"pDco - conj(cDpe) "<< pDco-conj(cDpo) <<std::endl;
std::cout<<GridLogMessage<<"=============================================================="<<std::endl; std::cout<<GridLogMessage<<"=============================================================="<<std::endl;
std::cout<<GridLogMessage<<"= Test MeeInv Mee = 1 "<<std::endl; std::cout<<GridLogMessage<<"= Test MeeInv Mee = 1 "<<std::endl;
std::cout<<GridLogMessage<<"=============================================================="<<std::endl; std::cout<<GridLogMessage<<"=============================================================="<<std::endl;

4
tests/debug/Test_cayley_mres.cc
View File

@ -118,8 +118,8 @@ int main (int argc, char ** argv)
} }
else else
{ {
std::cout<<GridLogMessage <<"Using cold configuration"<<std::endl; std::cout<<GridLogMessage <<"Using hot configuration"<<std::endl;
// SU<Nc>::ColdConfiguration(Umu); // SU<Nc>::ColdConfiguration(Umu);
SU<Nc>::HotConfiguration(RNG4,Umu); SU<Nc>::HotConfiguration(RNG4,Umu);
} }

13
tests/debug/Test_heatbath_dwf_eofa.cc
View File

@ -66,12 +66,9 @@ int main(int argc, char** argv)
// Set up RNGs // Set up RNGs
std::vector<int> seeds4({1, 2, 3, 4}); std::vector<int> seeds4({1, 2, 3, 4});
std::vector<int> seeds5({5, 6, 7, 8}); std::vector<int> seeds5({5, 6, 7, 8});
GridSerialRNG sRNG; GridSerialRNG SRNG; SRNG.SeedFixedIntegers(seeds5);
GridParallelRNG RNG5(FGrid); GridParallelRNG RNG5(FGrid); RNG5.SeedFixedIntegers(seeds5);
sRNG.SeedFixedIntegers(seeds5); GridParallelRNG RNG4(UGrid); RNG4.SeedFixedIntegers(seeds4);
RNG5.SeedFixedIntegers(seeds5);
GridParallelRNG RNG4(UGrid);
RNG4.SeedFixedIntegers(seeds4);
// Random gauge field // Random gauge field
LatticeGaugeField Umu(UGrid); LatticeGaugeField Umu(UGrid);
@ -86,7 +83,7 @@ int main(int argc, char** argv)
ConjugateGradient<LatticeFermion> CG(1.0e-12, 5000); ConjugateGradient<LatticeFermion> CG(1.0e-12, 5000);
ExactOneFlavourRatioPseudoFermionAction<WilsonImplR> Meofa(Lop, Rop, CG, Params, false); ExactOneFlavourRatioPseudoFermionAction<WilsonImplR> Meofa(Lop, Rop, CG, Params, false);
Meofa.refresh(Umu,sRNG, RNG5); Meofa.refresh(Umu, SRNG, RNG5);
printf("<Phi|Meofa|Phi> = %1.15e\n", Meofa.S(Umu)); printf("<Phi|Meofa|Phi> = %1.15e\n", Meofa.S(Umu));
} }
@ -96,7 +93,7 @@ int main(int argc, char** argv)
ConjugateGradient<LatticeFermion> CG(1.0e-12, 5000); ConjugateGradient<LatticeFermion> CG(1.0e-12, 5000);
ExactOneFlavourRatioPseudoFermionAction<WilsonImplR> Meofa(Lop, Rop, CG, Params, true); ExactOneFlavourRatioPseudoFermionAction<WilsonImplR> Meofa(Lop, Rop, CG, Params, true);
Meofa.refresh(Umu,sRNG, RNG5); Meofa.refresh(Umu, SRNG, RNG5);
printf("<Phi|Meofa|Phi> = %1.15e\n", Meofa.S(Umu)); printf("<Phi|Meofa|Phi> = %1.15e\n", Meofa.S(Umu));
} }

14
tests/debug/Test_heatbath_dwf_eofa_gparity.cc
View File

@ -70,13 +70,9 @@ int main(int argc, char** argv)
// Set up RNGs // Set up RNGs
std::vector<int> seeds4({1, 2, 3, 4}); std::vector<int> seeds4({1, 2, 3, 4});
std::vector<int> seeds5({5, 6, 7, 8}); std::vector<int> seeds5({5, 6, 7, 8});
GridParallelRNG RNG5(FGrid); GridSerialRNG SRNG; SRNG.SeedFixedIntegers(seeds5);
RNG5.SeedFixedIntegers(seeds5); GridParallelRNG RNG5(FGrid); RNG5.SeedFixedIntegers(seeds5);
GridParallelRNG RNG4(UGrid); GridParallelRNG RNG4(UGrid); RNG4.SeedFixedIntegers(seeds4);
RNG4.SeedFixedIntegers(seeds4);
GridSerialRNG sRNG;
RNG4.SeedFixedIntegers(seeds4);
sRNG.SeedFixedIntegers(seeds5);
// Random gauge field // Random gauge field
LatticeGaugeField Umu(UGrid); LatticeGaugeField Umu(UGrid);
@ -93,7 +89,7 @@ int main(int argc, char** argv)
ConjugateGradient<FermionField> CG(1.0e-12, 5000); ConjugateGradient<FermionField> CG(1.0e-12, 5000);
ExactOneFlavourRatioPseudoFermionAction<FermionImplPolicy> Meofa(Lop, Rop, CG, Params, false); ExactOneFlavourRatioPseudoFermionAction<FermionImplPolicy> Meofa(Lop, Rop, CG, Params, false);
Meofa.refresh(Umu,sRNG, RNG5); Meofa.refresh(Umu, SRNG, RNG5);
printf("<Phi|Meofa|Phi> = %1.15e\n", Meofa.S(Umu)); printf("<Phi|Meofa|Phi> = %1.15e\n", Meofa.S(Umu));
} }
@ -103,7 +99,7 @@ int main(int argc, char** argv)
ConjugateGradient<FermionField> CG(1.0e-12, 5000); ConjugateGradient<FermionField> CG(1.0e-12, 5000);
ExactOneFlavourRatioPseudoFermionAction<FermionImplPolicy> Meofa(Lop, Rop, CG, Params, true); ExactOneFlavourRatioPseudoFermionAction<FermionImplPolicy> Meofa(Lop, Rop, CG, Params, true);
Meofa.refresh(Umu,sRNG, RNG5); Meofa.refresh(Umu, SRNG, RNG5);
printf("<Phi|Meofa|Phi> = %1.15e\n", Meofa.S(Umu)); printf("<Phi|Meofa|Phi> = %1.15e\n", Meofa.S(Umu));
} }

7
tests/debug/Test_heatbath_mobius_eofa.cc
View File

@ -68,10 +68,9 @@ int main(int argc, char** argv)
// Set up RNGs // Set up RNGs
std::vector<int> seeds4({1, 2, 3, 4}); std::vector<int> seeds4({1, 2, 3, 4});
std::vector<int> seeds5({5, 6, 7, 8}); std::vector<int> seeds5({5, 6, 7, 8});
GridSerialRNG sRNG; GridSerialRNG SRNG; SRNG.SeedFixedIntegers(seeds4);
GridParallelRNG RNG5(FGrid); GridParallelRNG RNG5(FGrid);
RNG5.SeedFixedIntegers(seeds5); RNG5.SeedFixedIntegers(seeds5);
sRNG.SeedFixedIntegers(seeds5);
GridParallelRNG RNG4(UGrid); GridParallelRNG RNG4(UGrid);
RNG4.SeedFixedIntegers(seeds4); RNG4.SeedFixedIntegers(seeds4);
@ -88,7 +87,7 @@ int main(int argc, char** argv)
ConjugateGradient<LatticeFermion> CG(1.0e-12, 5000); ConjugateGradient<LatticeFermion> CG(1.0e-12, 5000);
ExactOneFlavourRatioPseudoFermionAction<WilsonImplR> Meofa(Lop, Rop, CG, Params, false); ExactOneFlavourRatioPseudoFermionAction<WilsonImplR> Meofa(Lop, Rop, CG, Params, false);
Meofa.refresh(Umu, sRNG,RNG5); Meofa.refresh(Umu, SRNG, RNG5);
printf("<Phi|Meofa|Phi> = %1.15e\n", Meofa.S(Umu)); printf("<Phi|Meofa|Phi> = %1.15e\n", Meofa.S(Umu));
} }
@ -98,7 +97,7 @@ int main(int argc, char** argv)
ConjugateGradient<LatticeFermion> CG(1.0e-12, 5000); ConjugateGradient<LatticeFermion> CG(1.0e-12, 5000);
ExactOneFlavourRatioPseudoFermionAction<WilsonImplR> Meofa(Lop, Rop, CG, Params, true); ExactOneFlavourRatioPseudoFermionAction<WilsonImplR> Meofa(Lop, Rop, CG, Params, true);
Meofa.refresh(Umu, sRNG,RNG5); Meofa.refresh(Umu, SRNG, RNG5);
printf("<Phi|Meofa|Phi> = %1.15e\n", Meofa.S(Umu)); printf("<Phi|Meofa|Phi> = %1.15e\n", Meofa.S(Umu));
} }

8
tests/debug/Test_heatbath_mobius_eofa_gparity.cc
View File

@ -72,10 +72,9 @@ int main(int argc, char** argv)
// Set up RNGs // Set up RNGs
std::vector<int> seeds4({1, 2, 3, 4}); std::vector<int> seeds4({1, 2, 3, 4});
std::vector<int> seeds5({5, 6, 7, 8}); std::vector<int> seeds5({5, 6, 7, 8});
GridSerialRNG SRNG; SRNG.SeedFixedIntegers(seeds5);
GridParallelRNG RNG5(FGrid); GridParallelRNG RNG5(FGrid);
GridSerialRNG sRNG;
RNG5.SeedFixedIntegers(seeds5); RNG5.SeedFixedIntegers(seeds5);
sRNG.SeedFixedIntegers(seeds5);
GridParallelRNG RNG4(UGrid); GridParallelRNG RNG4(UGrid);
RNG4.SeedFixedIntegers(seeds4); RNG4.SeedFixedIntegers(seeds4);
@ -93,7 +92,8 @@ int main(int argc, char** argv)
ConjugateGradient<FermionField> CG(1.0e-12, 5000); ConjugateGradient<FermionField> CG(1.0e-12, 5000);
ExactOneFlavourRatioPseudoFermionAction<FermionImplPolicy> Meofa(Lop, Rop, CG, Params, false); ExactOneFlavourRatioPseudoFermionAction<FermionImplPolicy> Meofa(Lop, Rop, CG, Params, false);
Meofa.refresh(Umu, sRNG, RNG5);
Meofa.refresh(Umu, SRNG, RNG5);
printf("<Phi|Meofa|Phi> = %1.15e\n", Meofa.S(Umu)); printf("<Phi|Meofa|Phi> = %1.15e\n", Meofa.S(Umu));
} }
@ -103,7 +103,7 @@ int main(int argc, char** argv)
ConjugateGradient<FermionField> CG(1.0e-12, 5000); ConjugateGradient<FermionField> CG(1.0e-12, 5000);
ExactOneFlavourRatioPseudoFermionAction<FermionImplPolicy> Meofa(Lop, Rop, CG, Params, true); ExactOneFlavourRatioPseudoFermionAction<FermionImplPolicy> Meofa(Lop, Rop, CG, Params, true);
Meofa.refresh(Umu, sRNG, RNG5); Meofa.refresh(Umu, SRNG, RNG5);
printf("<Phi|Meofa|Phi> = %1.15e\n", Meofa.S(Umu)); printf("<Phi|Meofa|Phi> = %1.15e\n", Meofa.S(Umu));
} }

8
tests/forces/Test_gp_rect_force.cc
View File

@ -52,15 +52,15 @@ int main (int argc, char ** argv)
pRNG.SeedFixedIntegers(std::vector<int>({45,12,81,9})); pRNG.SeedFixedIntegers(std::vector<int>({45,12,81,9}));
LatticeGaugeField U(&Grid); LatticeGaugeField U(&Grid);
SU<Nc>::HotConfiguration(pRNG,U); SU<Nc>::HotConfiguration(pRNG,U);
double beta = 1.0; double beta = 1.0;
double c1 = 0.331; double c1 = 0.331;
const int nu = 1; const int nu = 1;
std::vector<int> twists(Nd,0); std::vector<int> twists(Nd,0);
twists[nu] = 1; twists[nu] = 1;
ConjugateGimplD::setDirections(twists); ConjugateGimplD::setDirections(twists);
ConjugatePlaqPlusRectangleActionR Action(beta,c1); ConjugatePlaqPlusRectangleActionR Action(beta,c1);
//ConjugateWilsonGaugeActionR Action(beta); //ConjugateWilsonGaugeActionR Action(beta);
@ -70,13 +70,13 @@ int main (int argc, char ** argv)
// get the deriv of phidag MdagM phi with respect to "U" // get the deriv of phidag MdagM phi with respect to "U"
LatticeGaugeField UdSdU(&Grid); LatticeGaugeField UdSdU(&Grid);
Action.deriv(U,UdSdU); Action.deriv(U,UdSdU);
//////////////////////////////////// ////////////////////////////////////
// Modify the gauge field a little // Modify the gauge field a little
//////////////////////////////////// ////////////////////////////////////
RealD dt = 0.01; RealD dt = 0.001;
LatticeColourMatrix mommu(&Grid); LatticeColourMatrix mommu(&Grid);
LatticeColourMatrix forcemu(&Grid); LatticeColourMatrix forcemu(&Grid);

4
tests/forces/Test_momentum_filter.cc
View File

@ -59,11 +59,13 @@ int main (int argc, char ** argv)
std::cout<<GridLogMessage << "Grid is setup to use "<<threads<<" threads"<<std::endl; std::cout<<GridLogMessage << "Grid is setup to use "<<threads<<" threads"<<std::endl;
std::vector<int> seeds({1,2,3,4}); std::vector<int> seeds({1,2,3,4});
std::vector<int> serial_seeds({5,6,7,8});
GridSerialRNG sRNG;
GridParallelRNG pRNG(&Grid); GridParallelRNG pRNG(&Grid);
GridSerialRNG sRNG; GridSerialRNG sRNG;
pRNG.SeedFixedIntegers(seeds); pRNG.SeedFixedIntegers(seeds);
sRNG.SeedFixedIntegers(seeds); sRNG.SeedFixedIntegers(serial_seeds);
typedef PeriodicGimplR Gimpl; typedef PeriodicGimplR Gimpl;
typedef WilsonGaugeAction<Gimpl> GaugeAction; typedef WilsonGaugeAction<Gimpl> GaugeAction;

4
tests/hmc/Test_hmc_GparityIwasakiGauge.cc
View File

@ -95,8 +95,10 @@ int main(int argc, char **argv) {
///////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////
// HMC parameters are serialisable // HMC parameters are serialisable
TheHMC.Parameters.MD.MDsteps = 20; TheHMC.Parameters.MD.MDsteps = 80;
TheHMC.Parameters.MD.trajL = 1.0; TheHMC.Parameters.MD.trajL = 1.0;
TheHMC.Parameters.Trajectories = 100;
TheHMC.Parameters.NoMetropolisUntil= 10;
TheHMC.ReadCommandLine(argc, argv); // these can be parameters from file TheHMC.ReadCommandLine(argc, argv); // these can be parameters from file

6
tests/hmc/Test_hmc_IwasakiGauge.cc
View File

@ -51,7 +51,7 @@ int main(int argc, char **argv) {
CheckpointerParameters CPparams; CheckpointerParameters CPparams;
CPparams.config_prefix = "ckpoint_lat"; CPparams.config_prefix = "ckpoint_lat";
CPparams.rng_prefix = "ckpoint_rng"; CPparams.rng_prefix = "ckpoint_rng";
CPparams.saveInterval = 20; CPparams.saveInterval = 1;
CPparams.format = "IEEE64BIG"; CPparams.format = "IEEE64BIG";
TheHMC.Resources.LoadNerscCheckpointer(CPparams); TheHMC.Resources.LoadNerscCheckpointer(CPparams);
@ -81,8 +81,10 @@ int main(int argc, char **argv) {
///////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////
// HMC parameters are serialisable // HMC parameters are serialisable
TheHMC.Parameters.MD.MDsteps = 20; TheHMC.Parameters.MD.MDsteps = 80;
TheHMC.Parameters.MD.trajL = 1.0; TheHMC.Parameters.MD.trajL = 1.0;
TheHMC.Parameters.Trajectories = 100;
TheHMC.Parameters.NoMetropolisUntil= 10;
TheHMC.ReadCommandLine(argc, argv); // these can be parameters from file TheHMC.ReadCommandLine(argc, argv); // these can be parameters from file
TheHMC.Run(); // no smearing TheHMC.Run(); // no smearing

4
tests/hmc/Test_hmc_ScalarActionNxN.cc
View File

@ -132,8 +132,8 @@ int main(int argc, char **argv) {
// Checkpointer definition // Checkpointer definition
CheckpointerParameters CPparams(Reader); CheckpointerParameters CPparams(Reader);
//TheHMC.Resources.LoadBinaryCheckpointer(CPparams); TheHMC.Resources.LoadBinaryCheckpointer(CPparams);
TheHMC.Resources.LoadScidacCheckpointer(CPparams, SPar); //TheHMC.Resources.LoadNerscCheckpointer(CPparams, SPar);
RNGModuleParameters RNGpar(Reader); RNGModuleParameters RNGpar(Reader);
TheHMC.Resources.SetRNGSeeds(RNGpar); TheHMC.Resources.SetRNGSeeds(RNGpar);

4
tests/hmc/Test_hmc_WG_Production.cc
View File

@ -74,10 +74,10 @@ int main(int argc, char **argv) {
// Checkpointer definition // Checkpointer definition
CheckpointerParameters CPparams(Reader); CheckpointerParameters CPparams(Reader);
//TheHMC.Resources.LoadNerscCheckpointer(CPparams); TheHMC.Resources.LoadNerscCheckpointer(CPparams);
// Store metadata in the Scidac checkpointer // Store metadata in the Scidac checkpointer
TheHMC.Resources.LoadScidacCheckpointer(CPparams, WilsonPar); //TheHMC.Resources.LoadScidacCheckpointer(CPparams, WilsonPar);
RNGModuleParameters RNGpar(Reader); RNGModuleParameters RNGpar(Reader);
TheHMC.Resources.SetRNGSeeds(RNGpar); TheHMC.Resources.SetRNGSeeds(RNGpar);