portelli/Grid
1
0
Fork 0
mirror of https://github.com/paboyle/Grid.git synced 2024-11-09 23:45:36 +00:00
Code Releases Activity

Plaq, Rectangle, Iwasaki, Symanzik and DBW2 workign and HMC regresses to http://arxiv.org/pdf/hep-lat/0610075.pdf

Browse Source
This commit is contained in:
Azusa Yamaguchi 2015-12-28 16:38:31 +00:00
parent 0afcf1cf13
commit 78c4e862ef
21 changed files with 1382 additions and 38 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

22
TODO
View File

@ -1,5 +1,27 @@
TODO:
---------------
* Forces; the UdSdU term in gauge force term is half of what I think it should
be. This is a consequence of taking ONLY the first term in:
dSg/dt = dU/dt dSdU + dUdag/dt dSdUdag
in the fermion force.
Now, S_mom = - tr Pmu Pmu ; Pmu anti-herm
.
d Smom/dt = - 2.0 tr Pmu Pmu = - dSg/dt = - tr Pmu [Umu dSdUmu + UmuDag dSdUmuDag]
.
=> Pmu = Umu dSdUmu
Where the norm is half expected.
This means we must double the force in the Test_xxx_force routines, and is the origin of the factor of two.
This 2x is applied by hand in the fermion routines and in the Test_rect_force routine.
Policies:
* Link smearing/boundary conds; Policy class based implementation ; framework more in place

2
lib/Init.cc
View File

@ -166,7 +166,7 @@ void Grid_init(int *argc,char ***argv)
std::cout<<GridLogMessage<<"--mpi n.n.n.n : default MPI decomposition"<<std::endl;
std::cout<<GridLogMessage<<"--omp n : default number of OMP threads"<<std::endl;
std::cout<<GridLogMessage<<"--grid n.n.n.n : default Grid size"<<std::endl;
std::cout<<GridLogMessage<<"--log list : comma separted list of streams from Error,Warning,Message,Performance,Iterative,Debug"<<std::endl;
std::cout<<GridLogMessage<<"--log list : comma separted list of streams from Error,Warning,Message,Performance,Iterative,Integrator,Debug"<<std::endl;
}
if( GridCmdOptionExists(*argv,*argv+*argc,"--log") ){

8
lib/Log.cc
View File

@ -11,6 +11,7 @@ GridLogger GridLogMessage (1,"Message");
GridLogger GridLogDebug (1,"Debug");
GridLogger GridLogPerformance(1,"Performance");
GridLogger GridLogIterative (1,"Iterative");
GridLogger GridLogIntegrator (1,"Integrator");
void GridLogConfigure(std::vector<std::string> &logstreams)
{
@ -20,6 +21,7 @@ void GridLogConfigure(std::vector<std::string> &logstreams)
GridLogIterative.Active(0);
GridLogDebug.Active(0);
GridLogPerformance.Active(0);
GridLogIntegrator.Active(0);
for(int i=0;i<logstreams.size();i++){
if ( logstreams[i]== std::string("Error") ) GridLogError.Active(1);
@ -28,6 +30,7 @@ void GridLogConfigure(std::vector<std::string> &logstreams)
if ( logstreams[i]== std::string("Iterative") ) GridLogIterative.Active(1);
if ( logstreams[i]== std::string("Debug") ) GridLogDebug.Active(1);
if ( logstreams[i]== std::string("Performance") ) GridLogPerformance.Active(1);
if ( logstreams[i]== std::string("Integrator" ) ) GridLogIntegrator.Active(1);
}
}
@ -52,11 +55,6 @@ void Grid_unquiesce_nodes(void)
#endif
}
std::ostream& operator<< (std::ostream& stream, const GridTime& time)
{
// stream << time.count()<<" ms";
return stream;
}
}

2
lib/Log.h
View File

@ -4,7 +4,6 @@ namespace Grid {
// Dress the output; use std::chrono for time stamping via the StopWatch class
std::ostream& operator<< (std::ostream& stream, const GridTime& time);
class Logger {
protected:
@ -49,6 +48,7 @@ extern GridLogger GridLogMessage;
extern GridLogger GridLogDebug ;
extern GridLogger GridLogPerformance;
extern GridLogger GridLogIterative ;
extern GridLogger GridLogIntegrator ;
}
#endif

2
lib/Make.inc
View File

@ -1,4 +1,4 @@
HFILES=./algorithms/approx/bigfloat.h ./algorithms/approx/bigfloat_double.h ./algorithms/approx/Chebyshev.h ./algorithms/approx/MultiShiftFunction.h ./algorithms/approx/Remez.h ./algorithms/approx/Zolotarev.h ./algorithms/CoarsenedMatrix.h ./algorithms/iterative/AdefGeneric.h ./algorithms/iterative/ConjugateGradient.h ./algorithms/iterative/ConjugateGradientMultiShift.h ./algorithms/iterative/ConjugateResidual.h ./algorithms/iterative/DenseMatrix.h ./algorithms/iterative/EigenSort.h ./algorithms/iterative/Francis.h ./algorithms/iterative/Householder.h ./algorithms/iterative/ImplicitlyRestartedLanczos.h ./algorithms/iterative/Matrix.h ./algorithms/iterative/MatrixUtils.h ./algorithms/iterative/NormalEquations.h ./algorithms/iterative/PrecConjugateResidual.h ./algorithms/iterative/PrecGeneralisedConjugateResidual.h ./algorithms/iterative/SchurRedBlack.h ./algorithms/LinearOperator.h ./algorithms/Preconditioner.h ./algorithms/SparseMatrix.h ./Algorithms.h ./AlignedAllocator.h ./cartesian/Cartesian_base.h ./cartesian/Cartesian_full.h ./cartesian/Cartesian_red_black.h ./Cartesian.h ./communicator/Communicator_base.h ./Communicator.h ./cshift/Cshift_common.h ./cshift/Cshift_mpi.h ./cshift/Cshift_none.h ./Cshift.h ./Grid.h ./Init.h ./lattice/Lattice_arith.h ./lattice/Lattice_base.h ./lattice/Lattice_comparison.h ./lattice/Lattice_comparison_utils.h ./lattice/Lattice_conformable.h ./lattice/Lattice_coordinate.h ./lattice/Lattice_ET.h ./lattice/Lattice_local.h ./lattice/Lattice_overload.h ./lattice/Lattice_peekpoke.h ./lattice/Lattice_reality.h ./lattice/Lattice_reduction.h ./lattice/Lattice_rng.h ./lattice/Lattice_trace.h ./lattice/Lattice_transfer.h ./lattice/Lattice_transpose.h ./lattice/Lattice_unary.h ./lattice/Lattice_where.h ./Lattice.h ./Log.h ./Old/Tensor_peek.h ./Old/Tensor_poke.h ./parallelIO/BinaryIO.h ./parallelIO/NerscIO.h ./PerfCount.h ./pugixml/pugixml.h ./qcd/action/ActionBase.h ./qcd/action/ActionParams.h ./qcd/action/Actions.h ./qcd/action/fermion/CayleyFermion5D.h ./qcd/action/fermion/ContinuedFractionFermion5D.h ./qcd/action/fermion/DomainWallFermion.h ./qcd/action/fermion/FermionOperator.h ./qcd/action/fermion/FermionOperatorImpl.h ./qcd/action/fermion/g5HermitianLinop.h ./qcd/action/fermion/MobiusFermion.h ./qcd/action/fermion/MobiusZolotarevFermion.h ./qcd/action/fermion/OverlapWilsonCayleyTanhFermion.h ./qcd/action/fermion/OverlapWilsonCayleyZolotarevFermion.h ./qcd/action/fermion/OverlapWilsonContfracTanhFermion.h ./qcd/action/fermion/OverlapWilsonContfracZolotarevFermion.h ./qcd/action/fermion/OverlapWilsonPartialFractionTanhFermion.h ./qcd/action/fermion/OverlapWilsonPartialFractionZolotarevFermion.h ./qcd/action/fermion/PartialFractionFermion5D.h ./qcd/action/fermion/ScaledShamirFermion.h ./qcd/action/fermion/ShamirZolotarevFermion.h ./qcd/action/fermion/WilsonCompressor.h ./qcd/action/fermion/WilsonFermion.h ./qcd/action/fermion/WilsonFermion5D.h ./qcd/action/fermion/WilsonKernels.h ./qcd/action/fermion/WilsonTMFermion.h ./qcd/action/gauge/WilsonGaugeAction.h ./qcd/action/pseudofermion/EvenOddSchurDifferentiable.h ./qcd/action/pseudofermion/OneFlavourEvenOddRational.h ./qcd/action/pseudofermion/OneFlavourEvenOddRationalRatio.h ./qcd/action/pseudofermion/OneFlavourRational.h ./qcd/action/pseudofermion/OneFlavourRationalRatio.h ./qcd/action/pseudofermion/TwoFlavour.h ./qcd/action/pseudofermion/TwoFlavourEvenOdd.h ./qcd/action/pseudofermion/TwoFlavourEvenOddRatio.h ./qcd/action/pseudofermion/TwoFlavourRatio.h ./qcd/hmc/HMC.h ./qcd/hmc/integrators/Integrator.h ./qcd/hmc/integrators/Integrator_algorithm.h ./qcd/QCD.h ./qcd/spin/Dirac.h ./qcd/spin/TwoSpinor.h ./qcd/utils/CovariantCshift.h ./qcd/utils/LinalgUtils.h ./qcd/utils/SpaceTimeGrid.h ./qcd/utils/SUn.h ./qcd/utils/WilsonLoops.h ./serialisation/BaseIO.h ./serialisation/BinaryIO.h ./serialisation/MacroMagic.h ./serialisation/Serialisation.h ./serialisation/TextIO.h ./serialisation/XmlIO.h ./simd/Avx512Asm.h ./simd/Grid_avx.h ./simd/Grid_avx512.h ./simd/Grid_empty.h ./simd/Grid_imci.h ./simd/Grid_neon.h ./simd/Grid_qpx.h ./simd/Grid_sse4.h ./simd/Grid_vector_types.h ./simd/Grid_vector_unops.h ./Simd.h ./stencil/Lebesgue.h ./Stencil.h ./tensors/Tensor_arith.h ./tensors/Tensor_arith_add.h ./tensors/Tensor_arith_mac.h ./tensors/Tensor_arith_mul.h ./tensors/Tensor_arith_scalar.h ./tensors/Tensor_arith_sub.h ./tensors/Tensor_class.h ./tensors/Tensor_determinant.h ./tensors/Tensor_exp.h ./tensors/Tensor_extract_merge.h ./tensors/Tensor_index.h ./tensors/Tensor_inner.h ./tensors/Tensor_logical.h ./tensors/Tensor_outer.h ./tensors/Tensor_reality.h ./tensors/Tensor_Ta.h ./tensors/Tensor_trace.h ./tensors/Tensor_traits.h ./tensors/Tensor_transpose.h ./tensors/Tensor_unary.h ./Tensors.h ./Threads.h ./Timer.h
HFILES=./algorithms/approx/bigfloat.h ./algorithms/approx/bigfloat_double.h ./algorithms/approx/Chebyshev.h ./algorithms/approx/MultiShiftFunction.h ./algorithms/approx/Remez.h ./algorithms/approx/Zolotarev.h ./algorithms/CoarsenedMatrix.h ./algorithms/iterative/AdefGeneric.h ./algorithms/iterative/ConjugateGradient.h ./algorithms/iterative/ConjugateGradientMultiShift.h ./algorithms/iterative/ConjugateResidual.h ./algorithms/iterative/DenseMatrix.h ./algorithms/iterative/EigenSort.h ./algorithms/iterative/Francis.h ./algorithms/iterative/Householder.h ./algorithms/iterative/ImplicitlyRestartedLanczos.h ./algorithms/iterative/Matrix.h ./algorithms/iterative/MatrixUtils.h ./algorithms/iterative/NormalEquations.h ./algorithms/iterative/PrecConjugateResidual.h ./algorithms/iterative/PrecGeneralisedConjugateResidual.h ./algorithms/iterative/SchurRedBlack.h ./algorithms/LinearOperator.h ./algorithms/Preconditioner.h ./algorithms/SparseMatrix.h ./Algorithms.h ./AlignedAllocator.h ./cartesian/Cartesian_base.h ./cartesian/Cartesian_full.h ./cartesian/Cartesian_red_black.h ./Cartesian.h ./communicator/Communicator_base.h ./Communicator.h ./cshift/Cshift_common.h ./cshift/Cshift_mpi.h ./cshift/Cshift_none.h ./Cshift.h ./Grid.h ./Init.h ./lattice/Lattice_arith.h ./lattice/Lattice_base.h ./lattice/Lattice_comparison.h ./lattice/Lattice_comparison_utils.h ./lattice/Lattice_conformable.h ./lattice/Lattice_coordinate.h ./lattice/Lattice_ET.h ./lattice/Lattice_local.h ./lattice/Lattice_overload.h ./lattice/Lattice_peekpoke.h ./lattice/Lattice_reality.h ./lattice/Lattice_reduction.h ./lattice/Lattice_rng.h ./lattice/Lattice_trace.h ./lattice/Lattice_transfer.h ./lattice/Lattice_transpose.h ./lattice/Lattice_unary.h ./lattice/Lattice_where.h ./Lattice.h ./Log.h ./Old/Tensor_peek.h ./Old/Tensor_poke.h ./parallelIO/BinaryIO.h ./parallelIO/NerscIO.h ./PerfCount.h ./pugixml/pugixml.h ./qcd/action/ActionBase.h ./qcd/action/ActionParams.h ./qcd/action/Actions.h ./qcd/action/fermion/CayleyFermion5D.h ./qcd/action/fermion/ContinuedFractionFermion5D.h ./qcd/action/fermion/DomainWallFermion.h ./qcd/action/fermion/FermionOperator.h ./qcd/action/fermion/FermionOperatorImpl.h ./qcd/action/fermion/g5HermitianLinop.h ./qcd/action/fermion/MobiusFermion.h ./qcd/action/fermion/MobiusZolotarevFermion.h ./qcd/action/fermion/OverlapWilsonCayleyTanhFermion.h ./qcd/action/fermion/OverlapWilsonCayleyZolotarevFermion.h ./qcd/action/fermion/OverlapWilsonContfracTanhFermion.h ./qcd/action/fermion/OverlapWilsonContfracZolotarevFermion.h ./qcd/action/fermion/OverlapWilsonPartialFractionTanhFermion.h ./qcd/action/fermion/OverlapWilsonPartialFractionZolotarevFermion.h ./qcd/action/fermion/PartialFractionFermion5D.h ./qcd/action/fermion/ScaledShamirFermion.h ./qcd/action/fermion/ShamirZolotarevFermion.h ./qcd/action/fermion/WilsonCompressor.h ./qcd/action/fermion/WilsonFermion.h ./qcd/action/fermion/WilsonFermion5D.h ./qcd/action/fermion/WilsonKernels.h ./qcd/action/fermion/WilsonTMFermion.h ./qcd/action/gauge/PlaqPlusRectangleAction.h ./qcd/action/gauge/WilsonGaugeAction.h ./qcd/action/pseudofermion/EvenOddSchurDifferentiable.h ./qcd/action/pseudofermion/OneFlavourEvenOddRational.h ./qcd/action/pseudofermion/OneFlavourEvenOddRationalRatio.h ./qcd/action/pseudofermion/OneFlavourRational.h ./qcd/action/pseudofermion/OneFlavourRationalRatio.h ./qcd/action/pseudofermion/TwoFlavour.h ./qcd/action/pseudofermion/TwoFlavourEvenOdd.h ./qcd/action/pseudofermion/TwoFlavourEvenOddRatio.h ./qcd/action/pseudofermion/TwoFlavourRatio.h ./qcd/hmc/HMC.h ./qcd/hmc/HmcRunner.h ./qcd/hmc/integrators/Integrator.h ./qcd/hmc/integrators/Integrator_algorithm.h ./qcd/hmc/NerscCheckpointer.h ./qcd/QCD.h ./qcd/spin/Dirac.h ./qcd/spin/TwoSpinor.h ./qcd/utils/CovariantCshift.h ./qcd/utils/LinalgUtils.h ./qcd/utils/SpaceTimeGrid.h ./qcd/utils/SUn.h ./qcd/utils/WilsonLoops.h ./serialisation/BaseIO.h ./serialisation/BinaryIO.h ./serialisation/MacroMagic.h ./serialisation/Serialisation.h ./serialisation/TextIO.h ./serialisation/XmlIO.h ./simd/Avx512Asm.h ./simd/Grid_avx.h ./simd/Grid_avx512.h ./simd/Grid_empty.h ./simd/Grid_imci.h ./simd/Grid_neon.h ./simd/Grid_qpx.h ./simd/Grid_sse4.h ./simd/Grid_vector_types.h ./simd/Grid_vector_unops.h ./Simd.h ./stencil/Lebesgue.h ./Stencil.h ./tensors/Tensor_arith.h ./tensors/Tensor_arith_add.h ./tensors/Tensor_arith_mac.h ./tensors/Tensor_arith_mul.h ./tensors/Tensor_arith_scalar.h ./tensors/Tensor_arith_sub.h ./tensors/Tensor_class.h ./tensors/Tensor_determinant.h ./tensors/Tensor_exp.h ./tensors/Tensor_extract_merge.h ./tensors/Tensor_index.h ./tensors/Tensor_inner.h ./tensors/Tensor_logical.h ./tensors/Tensor_outer.h ./tensors/Tensor_reality.h ./tensors/Tensor_Ta.h ./tensors/Tensor_trace.h ./tensors/Tensor_traits.h ./tensors/Tensor_transpose.h ./tensors/Tensor_unary.h ./Tensors.h ./Threads.h ./Timer.h
CCFILES=./algorithms/approx/MultiShiftFunction.cc ./algorithms/approx/Remez.cc ./algorithms/approx/Zolotarev.cc ./Init.cc ./Log.cc ./PerfCount.cc ./pugixml/pugixml.cc ./qcd/action/fermion/CayleyFermion5D.cc ./qcd/action/fermion/ContinuedFractionFermion5D.cc ./qcd/action/fermion/PartialFractionFermion5D.cc ./qcd/action/fermion/WilsonFermion.cc ./qcd/action/fermion/WilsonFermion5D.cc ./qcd/action/fermion/WilsonKernels.cc ./qcd/action/fermion/WilsonKernelsAsm.cc ./qcd/action/fermion/WilsonKernelsHand.cc ./qcd/action/fermion/WilsonTMFermion.cc ./qcd/hmc/HMC.cc ./qcd/spin/Dirac.cc ./qcd/utils/SpaceTimeGrid.cc ./serialisation/BinaryIO.cc ./serialisation/TextIO.cc ./serialisation/XmlIO.cc ./stencil/Lebesgue.cc ./stencil/Stencil_common.cc

5
lib/Timer.h
View File

@ -17,6 +17,11 @@ typedef std::chrono::system_clock GridClock;
typedef std::chrono::time_point<GridClock> GridTimePoint;
typedef std::chrono::milliseconds GridTime;
inline std::ostream& operator<< (std::ostream & stream, const std::chrono::milliseconds & time)
{
stream << time.count()<<" ms";
return stream;
}
class GridStopWatch {
private:

10
lib/qcd/action/Actions.h
View File

@ -28,11 +28,21 @@
// Gauge Actions
////////////////////////////////////////////
#include <qcd/action/gauge/WilsonGaugeAction.h>
#include <qcd/action/gauge/PlaqPlusRectangleAction.h>
namespace Grid {
namespace QCD {
typedef WilsonGaugeAction<LatticeGaugeField> WilsonGaugeActionR;
typedef WilsonGaugeAction<LatticeGaugeFieldF> WilsonGaugeActionF;
typedef WilsonGaugeAction<LatticeGaugeFieldD> WilsonGaugeActionD;
typedef PlaqPlusRectangleAction<LatticeGaugeField> PlaqPlusRectangleActionR;
typedef PlaqPlusRectangleAction<LatticeGaugeFieldF> PlaqPlusRectangleActionF;
typedef PlaqPlusRectangleAction<LatticeGaugeFieldD> PlaqPlusRectangleActionD;
typedef IwasakiGaugeAction<LatticeGaugeField> IwasakiGaugeActionR;
typedef IwasakiGaugeAction<LatticeGaugeFieldF> IwasakiGaugeActionF;
typedef IwasakiGaugeAction<LatticeGaugeFieldD> IwasakiGaugeActionD;
typedef SymanzikGaugeAction<LatticeGaugeField> SymanzikGaugeActionR;
typedef SymanzikGaugeAction<LatticeGaugeFieldF> SymanzikGaugeActionF;
typedef SymanzikGaugeAction<LatticeGaugeFieldD> SymanzikGaugeActionD;
}}
////////////////////////////////////////////////////////////////////////////////////////////////////

7
lib/qcd/action/gauge/WilsonGaugeAction.h
View File

@ -22,17 +22,18 @@ namespace Grid{
virtual RealD S(const GaugeField &U) {
RealD plaq = WilsonLoops<GaugeField>::avgPlaquette(U);
// std::cout<<GridLogMessage << "Plaq : "<<plaq << "\n";
RealD vol = U._grid->gSites();
RealD action=beta*(1.0 -plaq)*(Nd*(Nd-1.0))*vol*0.5;
// std::cout << GridLogMessage << "WilsonGauge action "<<action<<std::endl;
return action;
};
virtual void deriv(const GaugeField &U,GaugeField & dSdU) {
//not optimal implementation FIXME
//extend Ta to include Lorentz indexes
//RealD factor = 0.5*beta/RealD(Nc);
RealD factor = 0.5*beta/RealD(Nc);
GaugeLinkField Umu(U._grid);
GaugeLinkField dSdU_mu(U._grid);
for (int mu=0; mu < Nd; mu++){
@ -41,7 +42,7 @@ namespace Grid{
// Staple in direction mu
WilsonLoops<GaugeField>::Staple(dSdU_mu,U,mu);
dSdU_mu = Ta(Umu*adj(dSdU_mu))*factor;
dSdU_mu = Ta(Umu*dSdU_mu)*factor;
PokeIndex<LorentzIndex>(dSdU, dSdU_mu, mu);
}
};

4
lib/qcd/hmc/HMC.h
View File

@ -53,8 +53,10 @@ namespace Grid{
std::string file; { std::ostringstream os; os << Stem <<"."<< traj; file = os.str(); }
std::ofstream of(file);
RealD plaq = WilsonLoops<GaugeField>::avgPlaquette(U);
of << plaq <<std::endl;
RealD rect = WilsonLoops<GaugeField>::avgRectangle(U);
of << plaq << " " << rect << std::endl;
std::cout<< GridLogMessage<< "Plaquette for trajectory "<< traj << " is " << plaq <<std::endl;
std::cout<< GridLogMessage<< "Rectangle for trajectory "<< traj << " is " << rect <<std::endl;
}
};

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

@ -77,7 +77,7 @@ namespace Grid{
t_P[level]+=ep;
update_P(P,U,level,ep);
std::cout<<GridLogMessage;
std::cout<<GridLogIntegrator;
for(int l=0; l<level;++l) std::cout<<" ";
std::cout<<"["<<level<<"] P " << " dt "<< ep <<" : t_P "<< t_P[level] <<std::endl;
}
@ -95,7 +95,7 @@ namespace Grid{
t_U+=ep;
int fl = levels-1;
std::cout<<GridLogMessage<<" ";
std::cout<<GridLogIntegrator<<" ";
for(int l=0; l<fl;++l) std::cout<<" ";
std::cout<<"["<<fl<<"] U " << " dt "<< ep <<" : t_U "<< t_U <<std::endl;
@ -133,7 +133,7 @@ namespace Grid{
//Initialization of momenta and actions
void refresh(GaugeField& U,GridParallelRNG &pRNG){
std::cout<<GridLogMessage<< "Integrator refresh\n";
std::cout<<GridLogIntegrator<< "Integrator refresh\n";
generate_momenta(P,pRNG);
for(int level=0; level< as.size(); ++level){
for(int actionID=0; actionID<as[level].actions.size(); ++actionID){
@ -186,7 +186,7 @@ namespace Grid{
// Check the clocks all match on all levels
for(int level=0; level<as.size(); ++level){
assert(fabs(t_U - t_P[level])<1.0e-6); // must be the same
std::cout<<GridLogMessage<<" times["<<level<<"]= "<<t_P[level]<< " " << t_U <<std::endl;
std::cout<<GridLogIntegrator<<" times["<<level<<"]= "<<t_P[level]<< " " << t_U <<std::endl;
}
// and that we indeed got to the end of the trajectory

158
lib/qcd/utils/WilsonLoops.h
View File

@ -108,14 +108,20 @@ public:
// mu
// ^
// |__ nu
// |__> nu
// __ __
// | |
// __| = __| *
//
staple += CovShiftForward(U[nu],nu,U[mu])*Cshift(adj(U[nu]),mu,+1);
// staple += CovShiftForward(U[nu],nu,U[mu])*Cshift(adj(U[nu]),mu,+1);
staple +=
Cshift(CovShiftForward (U[nu],nu,
CovShiftBackward(U[mu],mu,Cshift(adj(U[nu]),nu,-1))),mu,1);
// Unu(x) Umu(x+nu) Unu^dag(x+mu) ; left mult by Umu^dag(x) to close ring
//
// __ __
@ -123,8 +129,15 @@ public:
// |__ = | * __
//
//
tmp = CovShiftForward (U[mu],mu,U[nu]);
staple+= CovShiftBackward(U[nu],nu,tmp);
staple +=
Cshift(
CovShiftBackward(U[nu],nu,
CovShiftBackward(U[mu],mu,U[nu])),mu,1);
// tmp = CovShiftForward (U[mu],mu,U[nu]);
// staple+= CovShiftBackward(U[nu],nu,tmp);
// Unu^dag(x-nu) Umu(x-nu) Unu(x+mu-nu) ; left mult by Umu^dag(x) to close ring.
}
}
@ -133,13 +146,136 @@ public:
//////////////////////////////////////////////////////
// Similar to above for rectangle is required
//////////////////////////////////////////////////////
/*
void siteRectangle(GaugeMat &plaq,const std::vector<GaugeMat> &U, const int mu, const int nu){
RealD avgRectangle(const std::vector<GaugeMat> &U){}
RealD avgRectangle(const std::vector<GaugeMat> &U, const int mu, const int nu){}
void traceRectangle(LatticeComplex &plaq,const std::vector<GaugeMat> &U, const int mu, const int nu){}
void siteRectangle(GaugeMat &plaq,const std::vector<GaugeMat> &U, const int mu, const int nu){}
*/
static void dirRectangle(GaugeMat &rect,const std::vector<GaugeMat> &U, const int mu, const int nu)
{
rect = CovShiftForward(U[mu],mu,CovShiftForward(U[mu],mu,U[nu]))* // ->->|
adj(CovShiftForward(U[nu],nu,CovShiftForward(U[mu],mu,U[mu]))) ;
rect = rect +
CovShiftForward(U[mu],mu,CovShiftForward(U[nu],nu,U[nu]))* // ->||
adj(CovShiftForward(U[nu],nu,CovShiftForward(U[nu],nu,U[mu]))) ;
}
static void traceDirRectangle(LatticeComplex &rect, const std::vector<GaugeMat> &U, const int mu, const int nu)
{
GaugeMat sp(U[0]._grid);
dirRectangle(sp,U,mu,nu);
rect=trace(sp);
}
static void siteRectangle(LatticeComplex &Rect,const std::vector<GaugeMat> &U)
{
LatticeComplex siteRect(U[0]._grid);
Rect=zero;
for(int mu=1;mu<Nd;mu++){
for(int nu=0;nu<mu;nu++){
traceDirRectangle(siteRect,U,mu,nu);
Rect = Rect + siteRect;
}
}
}
//////////////////////////////////////////////////
// sum over all x,y,z,t and over all planes of plaquette
//////////////////////////////////////////////////
static RealD sumRectangle(const GaugeLorentz &Umu){
std::vector<GaugeMat> U(4,Umu._grid);
for(int mu=0;mu<Nd;mu++){
U[mu] = PeekIndex<LorentzIndex>(Umu,mu);
}
LatticeComplex Rect(Umu._grid);
siteRectangle(Rect,U);
TComplex Tp = sum(Rect);
Complex p = TensorRemove(Tp);
return p.real();
}
//////////////////////////////////////////////////
// average over all x,y,z,t and over all planes of plaquette
//////////////////////////////////////////////////
static RealD avgRectangle(const GaugeLorentz &Umu){
RealD sumrect = sumRectangle(Umu);
double vol = Umu._grid->gSites();
double faces = (1.0*Nd*(Nd-1)); // 2 distinct orientations summed
return sumrect/vol/faces/Nc; // Nd , Nc dependent... FIXME
}
//////////////////////////////////////////////////
// the sum over all staples on each site
//////////////////////////////////////////////////
static void RectStaple(GaugeMat &Stap,const GaugeLorentz &Umu,int mu){
GridBase *grid = Umu._grid;
std::vector<GaugeMat> U(4,grid);
for(int d=0;d<Nd;d++){
U[d] = PeekIndex<LorentzIndex>(Umu,d);
}
Stap=zero;
for(int nu=0;nu<Nd;nu++){
if ( nu!=mu) {
// __ ___
// | __ |
//
Stap+= Cshift(CovShiftForward (U[mu],mu,
CovShiftForward (U[nu],nu,
CovShiftBackward(U[mu],mu,
CovShiftBackward(U[mu],mu,
Cshift(adj(U[nu]),nu,-1))))) , mu, 1);
// __
// |__ __ |
Stap+= Cshift(CovShiftForward (U[mu],mu,
CovShiftBackward(U[nu],nu,
CovShiftBackward(U[mu],mu,
CovShiftBackward(U[mu],mu, U[nu])))) , mu, 1);
// __
// |__ __ |
Stap+= Cshift(CovShiftBackward(U[nu],nu,
CovShiftBackward(U[mu],mu,
CovShiftBackward(U[mu],mu,
CovShiftForward(U[nu],nu,U[mu])))) , mu, 1);
// __ ___
// |__ |
Stap+= Cshift(CovShiftForward (U[nu],nu,
CovShiftBackward(U[mu],mu,
CovShiftBackward(U[mu],mu,
CovShiftBackward(U[nu],nu,U[mu])))) , mu, 1);
// --
// | |
//
// | |
Stap+= Cshift(CovShiftForward(U[nu],nu,
CovShiftForward(U[nu],nu,
CovShiftBackward(U[mu],mu,
CovShiftBackward(U[nu],nu,
Cshift(adj(U[nu]),nu,-1))))) , mu, 1);
// | |
//
// | |
// --
Stap+= Cshift(CovShiftBackward(U[nu],nu,
CovShiftBackward(U[nu],nu,
CovShiftBackward(U[mu],mu,
CovShiftForward (U[nu],nu,U[nu])))) , mu, 1);
}}
}
};

18
tests/Make.inc
View File

@ -1,11 +1,15 @@
bin_PROGRAMS = Test_GaugeAction Test_cayley_cg Test_cayley_coarsen_support Test_cayley_even_odd Test_cayley_ldop_cr Test_cf_coarsen_support Test_cf_cr_unprec Test_cheby Test_contfrac_cg Test_contfrac_even_odd Test_contfrac_force Test_cshift Test_cshift_red_black Test_dwf_cg_prec Test_dwf_cg_schur Test_dwf_cg_unprec Test_dwf_cr_unprec Test_dwf_even_odd Test_dwf_force Test_dwf_fpgcr Test_dwf_gpforce Test_dwf_hdcr Test_dwf_lanczos Test_gamma Test_gparity Test_gpdwf_force Test_gpwilson_even_odd Test_hmc_EODWFRatio Test_hmc_EODWFRatio_Gparity Test_hmc_EOWilsonFermionGauge Test_hmc_EOWilsonRatio Test_hmc_WilsonFermionGauge Test_hmc_WilsonGauge Test_hmc_WilsonRatio Test_lie_generators Test_main Test_multishift_sqrt Test_nersc_io Test_partfrac_force Test_quenched_update Test_remez Test_rhmc_EOWilson1p1 Test_rhmc_EOWilsonRatio Test_rhmc_Wilson1p1 Test_rhmc_WilsonRatio Test_rng Test_rng_fixed Test_serialisation Test_simd Test_stencil Test_synthetic_lanczos Test_wilson_cg_prec Test_wilson_cg_schur Test_wilson_cg_unprec Test_wilson_cr_unprec Test_wilson_even_odd Test_wilson_force Test_wilson_force_phiMdagMphi Test_wilson_force_phiMphi Test_wilson_tm_even_odd
bin_PROGRAMS = Test_GaugeAction Test_RectPlaq Test_cayley_cg Test_cayley_coarsen_support Test_cayley_even_odd Test_cayley_ldop_cr Test_cf_coarsen_support Test_cf_cr_unprec Test_cheby Test_contfrac_cg Test_contfrac_even_odd Test_contfrac_force Test_cshift Test_cshift_red_black Test_dwf_cg_prec Test_dwf_cg_schur Test_dwf_cg_unprec Test_dwf_cr_unprec Test_dwf_even_odd Test_dwf_force Test_dwf_fpgcr Test_dwf_gpforce Test_dwf_hdcr Test_dwf_lanczos Test_gamma Test_gparity Test_gpdwf_force Test_gpwilson_even_odd Test_hmc_EODWFRatio Test_hmc_EODWFRatio_Gparity Test_hmc_EOWilsonFermionGauge Test_hmc_EOWilsonRatio Test_hmc_IwasakiGauge Test_hmc_RectGauge Test_hmc_WilsonFermionGauge Test_hmc_WilsonGauge Test_hmc_WilsonRatio Test_lie_generators Test_main Test_multishift_sqrt Test_nersc_io Test_partfrac_force Test_quenched_update Test_rect_force Test_remez Test_rhmc_EOWilson1p1 Test_rhmc_EOWilsonRatio Test_rhmc_Wilson1p1 Test_rhmc_WilsonRatio Test_rng Test_rng_fixed Test_serialisation Test_simd Test_stencil Test_synthetic_lanczos Test_wilson_cg_prec Test_wilson_cg_schur Test_wilson_cg_unprec Test_wilson_cr_unprec Test_wilson_even_odd Test_wilson_force Test_wilson_force_phiMdagMphi Test_wilson_force_phiMphi Test_wilson_tm_even_odd
Test_GaugeAction_SOURCES=Test_GaugeAction.cc
Test_GaugeAction_LDADD=-lGrid
Test_RectPlaq_SOURCES=Test_RectPlaq.cc
Test_RectPlaq_LDADD=-lGrid
Test_cayley_cg_SOURCES=Test_cayley_cg.cc
Test_cayley_cg_LDADD=-lGrid
@ -126,6 +130,14 @@ Test_hmc_EOWilsonRatio_SOURCES=Test_hmc_EOWilsonRatio.cc
Test_hmc_EOWilsonRatio_LDADD=-lGrid
Test_hmc_IwasakiGauge_SOURCES=Test_hmc_IwasakiGauge.cc
Test_hmc_IwasakiGauge_LDADD=-lGrid
Test_hmc_RectGauge_SOURCES=Test_hmc_RectGauge.cc
Test_hmc_RectGauge_LDADD=-lGrid
Test_hmc_WilsonFermionGauge_SOURCES=Test_hmc_WilsonFermionGauge.cc
Test_hmc_WilsonFermionGauge_LDADD=-lGrid
@ -162,6 +174,10 @@ Test_quenched_update_SOURCES=Test_quenched_update.cc
Test_quenched_update_LDADD=-lGrid
Test_rect_force_SOURCES=Test_rect_force.cc
Test_rect_force_LDADD=-lGrid
Test_remez_SOURCES=Test_remez.cc
Test_remez_LDADD=-lGrid

2
tests/Test_GaugeAction.cc
View File

@ -125,7 +125,7 @@ int main (int argc, char ** argv)
LatticeComplex stap_tr(&Fine);
for(int mu=0;mu<Nd;mu++){
ColourWilsonLoops::Staple(stap,Umu,mu);
stap_tr = trace(stap*adj(U[mu]));
stap_tr = trace(stap*U[mu]);
TComplex Ts = sum(stap_tr);
Complex s = TensorRemove(Ts);
stap_plaq+=real(s);

299
tests/Test_RectPlaq.cc Normal file
View File

@ -0,0 +1,299 @@
#include <Grid.h>
#include <qcd/utils/CovariantCshift.h>
#include <qcd/utils/WilsonLoops.h>
using namespace std;
using namespace Grid;
using namespace Grid::QCD;
/* For Metropolis */
class Metropolis {
public:
GridSerialRNG & sRNG;
Metropolis(GridSerialRNG & _sRNG) : sRNG(_sRNG) {};
bool AcceptReject(const RealD Delta)
{
RealD rand;
if(Delta <=0.0) return true;
random(sRNG,rand);
if(rand <= exp(-Delta))
return true;
else
return false;
}
};
void RectPlaq(const std::vector<LatticeColourMatrix> &U, LatticeComplex &RectPlaqValue )
{
RectPlaqValue=zero;
// 12 * vol loops
for(int mu=1;mu<Nd;mu++){
for(int nu=0;nu<mu;nu++){
RectPlaqValue = RectPlaqValue + trace(
CovShiftForward(U[mu],mu,CovShiftForward(U[mu],mu,U[nu]))* // ->->|
adj(CovShiftForward(U[nu],nu,CovShiftForward(U[mu],mu,U[mu]))) );
RectPlaqValue = RectPlaqValue + trace(
CovShiftForward(U[mu],mu,CovShiftForward(U[nu],nu,U[nu]))* // ->||
adj(CovShiftForward(U[nu],nu,CovShiftForward(U[nu],nu,U[mu]))) );
}
}
}
void RectPlaqDeriv(const std::vector<LatticeColourMatrix> &U, LatticeComplex &RectPlaqValue )
{
}
int main (int argc, char ** argv)
{
Grid_init(&argc,&argv);
std::vector<int> simd_layout = GridDefaultSimd(4,vComplex::Nsimd());
std::vector<int> mpi_layout = GridDefaultMpi();
std::vector<int> latt_size ({16,16,16,32});
std::vector<int> clatt_size ({4,4,4,8});
int orthodir=3;
int orthosz =latt_size[orthodir];
GridCartesian Fine(latt_size,simd_layout,mpi_layout);
GridCartesian Coarse(clatt_size,simd_layout,mpi_layout);
LatticeGaugeField Umu(&Fine);
std::vector<LatticeColourMatrix> U(4,&Fine);
NerscField header;
std::string file("./ckpoint_lat.4000");
NerscIO::readConfiguration(Umu,header,file);
for(int mu=0;mu<Nd;mu++){
U[mu] = PeekIndex<LorentzIndex>(Umu,mu);
}
// Painful ; fix syntactical niceness : to check reader
LatticeComplex LinkTrace(&Fine);
LinkTrace=zero;
for(int mu=0;mu<Nd;mu++){
LinkTrace = LinkTrace + trace(U[mu]);
}
LatticeComplex Plaq(&Fine);
LatticeComplex cPlaq(&Coarse);
Plaq = zero;
for(int mu=1;mu<Nd;mu++){
for(int nu=0;nu<mu;nu++){
Plaq = Plaq + trace(CovShiftForward(U[mu],mu,U[nu])*adj(CovShiftForward(U[nu],nu,U[mu])));
}
}
LatticeComplex RectPlaqValue(&Fine);
double vol = Fine.gSites();
Complex PlaqScale(1.0/vol/6.0/3.0);
Complex RectScale(1.0/vol/12.0/3.0);
std::cout<<GridLogMessage <<"PlaqScale" << PlaqScale<<std::endl;
std::cout<<GridLogMessage <<"RectScale" << RectScale<<std::endl;
RectPlaq(U,RectPlaqValue);
TComplex TRp = sum(RectPlaqValue);
Complex rp = TensorRemove(TRp);
std::cout<<GridLogMessage << "calculated Rect plaquettes A " <<rp*RectScale<<std::endl;
// Rect Plaq Calc Deriv
LatticeComplex RectPlaq_d(&Fine);
RectPlaq_d = zero;
LatticeColourMatrix ds_U(&Fine);
LatticeColourMatrix left_2(&Fine);
LatticeColourMatrix upper_l(&Fine);
LatticeColourMatrix upper_staple(&Fine);
LatticeColourMatrix down_staple(&Fine);
LatticeColourMatrix tmp(&Fine);
// 2x1 // Each link has 2*(Nd-1) + 4*(Nd-1) = 6(Nd-1) , 1x2 and 2x1 loops attached.
// //
// // For producing the rectangle term normalised to number of loops
// // there are Vol x Nd.(Nd-1) x 2 / 2 distinct loops total. (mu<nu, mu>nu)
// //
// // Expect scale factor to be
// //
for(int mu=0;mu<Nd;mu++){
ds_U=zero; // dS / dUmu
for(int nu=0;nu<Nd;nu++){
if ( nu != mu ) {
/*
(x) ---> ---> : U(x,mu)*U(x+mu, mu)
*/
left_2= CovShiftForward(U[mu],mu,U[mu]);
/*
upper_l = <---- <---
^
| =>tmp
(x+2mu)
Unu(x+2mu) Umudag(x+mu+nu) Umudag(x+nu)
*/
tmp=Cshift(U[nu],mu,2);
upper_l= CovShiftForward(tmp,nu,adj(left_2)); // i.e. upper_l
/*
upper_staple= <---- <--- ^
| |
V (x) (x + 2mu)
*/
// Unu(x+2mu) Umudag(x+mu+nu) Umudag(x+nu) Unudag(x)
upper_staple= upper_l*adj(U[nu]);
/*
down_staple= ^
| |
(x) <----- <---- V x + 2mu
*/
down_staple= adj(left_2*tmp)*U[nu];
/*
ds_U+= <---- <--- ^
| |
(x-mu) V-----> (x + mu)
*/
tmp=upper_staple*U[mu];
ds_U+= Cshift(tmp,mu,-1);
/*
ds_U+= (x-mu) ^----> (x + mu)
| |
<-----<--- V
*/
tmp=CovShiftBackward(U[mu],nu,down_staple);
ds_U+=Cshift(tmp,mu,-1);
/*
ds_U+= <----<---- ^
| |
(x) V -----> (x + 2mu)
*/
tmp=Cshift(U[mu],mu,1);
/*
ds_U+= (x) ^ ----> (x + 2mu)
| |
<---- <----V
*/
ds_U+=tmp*(upper_staple+down_staple);
/*****Part 2********/
/*
^
|
upper= ^
|
(x)
*/
LatticeColourMatrix up2= CovShiftForward(U[nu],nu,U[nu]);
/*
<----^
|
upper_l= ^
|
(x)
*/
// Unu(x+mu)Unu(x+mu+nu) UmuDag(x+nu+nu) lives at X
upper_l= CovShiftForward(Cshift(up2,mu,1),nu,Cshift(adj(U[mu]),nu,1));
/*
|<----^
upper_staple = V |
| ^
(x) V |
*/
ds_U+= upper_l*adj(up2);
/*
|
V
downer_l= |
(x)<----V
*/
upper_l= adj(CovShiftForward(U[mu],mu,up2)); //downer_l
/*
^ |
down_staple = | V
^ |
| V
(x)<----
down_staple= upper*upper_l;
*/
tmp= upper_l*up2;
ds_U+= Cshift(tmp,nu,-2);
//TRp = sum(RectPlaq_d);
//rp = TensorRemove(TRp);
//std::cout << GridLogMessage<< " Rect[" << " " << "] = "<< TensorRemove(TRp) <<std::endl;
}}
RectPlaq_d += trace( U[mu]*ds_U) * 0.25;
}
TRp = sum(RectPlaq_d);
rp = TensorRemove(TRp);
std::cout<<GridLogMessage << "calculated Rect plaquettes_d " <<rp*RectScale<<std::endl;
std::vector<TComplex> Plaq_T(orthosz);
sliceSum(Plaq,Plaq_T,Nd-1);
int Nt = Plaq_T.size();
TComplex Plaq_T_sum;
Plaq_T_sum=zero;
for(int t=0;t<Nt;t++){
Plaq_T_sum = Plaq_T_sum+Plaq_T[t];
Complex Pt=TensorRemove(Plaq_T[t]);
std::cout<<GridLogMessage << "sliced ["<<t<<"]" <<Pt*PlaqScale*Real(Nt) << std::endl;
}
{
Complex Pt = TensorRemove(Plaq_T_sum);
std::cout<<GridLogMessage << "total " <<Pt*PlaqScale<<std::endl;
}
TComplex Tp = sum(Plaq);
Complex p = TensorRemove(Tp);
std::cout<<GridLogMessage << "calculated plaquettes " <<p*PlaqScale<<std::endl;
RealD avg_plaq = ColourWilsonLoops::avgPlaquette(Umu);
std::cout<<GridLogMessage << "NEW : calculated real plaquettes " <<avg_plaq<<std::endl;
// Staple Plaq
RealD StapScale(1.0/vol/6.0/3.0);
RealD stap_plaq=0.0;
LatticeColourMatrix stap(&Fine);
LatticeComplex stap_tr(&Fine);
for(int mu=0;mu<Nd;mu++){
ColourWilsonLoops::Staple(stap,Umu,mu);
stap_tr = trace(U[mu]*stap);
TComplex Ts = sum(stap_tr);
Complex s = TensorRemove(Ts);
stap_plaq+=real(s);
}
std::cout<<GridLogMessage << "NEW : plaquette via staples"<< stap_plaq*StapScale*0.25<< std::endl;
Complex LinkTraceScale(1.0/vol/4.0/3.0);
TComplex Tl = sum(LinkTrace);
Complex l = TensorRemove(Tl);
std::cout<<GridLogMessage << "calculated link trace " <<l*LinkTraceScale<<std::endl;
blockSum(cPlaq,Plaq);
TComplex TcP = sum(cPlaq);
Complex ll= TensorRemove(TcP);
std::cout<<GridLogMessage << "coarsened plaquettes sum to " <<ll*PlaqScale<<std::endl;
Grid_finalize();
}

57
tests/Test_hmc_IwasakiGauge.cc Normal file
View File

@ -0,0 +1,57 @@
#include "Grid.h"
using namespace std;
using namespace Grid;
using namespace Grid::QCD;
namespace Grid {
namespace QCD {
class HmcRunner : public NerscHmcRunner {
public:
void BuildTheAction (int argc, char **argv)
{
typedef WilsonImplR ImplPolicy;
typedef WilsonFermionR FermionAction;
typedef typename FermionAction::FermionField FermionField;
UGrid = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi());
UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
FGrid = UGrid;
FrbGrid = UrbGrid;
// temporarily need a gauge field
LatticeGaugeField U(UGrid);
// Gauge action
IwasakiGaugeActionR Gaction(2.6);
//Collect actions
ActionLevel<LatticeGaugeField> Level1(1);
Level1.push_back(&Gaction);
TheAction.push_back(Level1);
Run(argc,argv);
};
};
}}
int main (int argc, char ** argv)
{
Grid_init(&argc,&argv);
int threads = GridThread::GetThreads();
std::cout<<GridLogMessage << "Grid is setup to use "<<threads<<" threads"<<std::endl;
HmcRunner TheHMC;
TheHMC.BuildTheAction(argc,argv);
}

57
tests/Test_hmc_RectGauge.cc Normal file
View File

@ -0,0 +1,57 @@
#include "Grid.h"
using namespace std;
using namespace Grid;
using namespace Grid::QCD;
namespace Grid {
namespace QCD {
class HmcRunner : public NerscHmcRunner {
public:
void BuildTheAction (int argc, char **argv)
{
typedef WilsonImplR ImplPolicy;
typedef WilsonFermionR FermionAction;
typedef typename FermionAction::FermionField FermionField;
UGrid = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi());
UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
FGrid = UGrid;
FrbGrid = UrbGrid;
// temporarily need a gauge field
LatticeGaugeField U(UGrid);
// Gauge action
PlaqPlusRectangleActionR Gaction(2.0,0.331);
//Collect actions
ActionLevel<LatticeGaugeField> Level1(1);
Level1.push_back(&Gaction);
TheAction.push_back(Level1);
Run(argc,argv);
};
};
}}
int main (int argc, char ** argv)
{
Grid_init(&argc,&argv);
int threads = GridThread::GetThreads();
std::cout<<GridLogMessage << "Grid is setup to use "<<threads<<" threads"<<std::endl;
HmcRunner TheHMC;
TheHMC.BuildTheAction(argc,argv);
}

1
tests/Test_quenched_update.cc
View File

@ -61,7 +61,6 @@ int main (int argc, char ** argv)
// Get Link and Staple term in action; must contain Beta and
// any other coeffs
ColourWilsonLoops::Staple(staple,Umu,mu);
staple = adj(staple);
link = PeekIndex<LorentzIndex>(Umu,mu);

97
tests/Test_rect_force.cc Normal file
View File

@ -0,0 +1,97 @@
#include <Grid.h>
using namespace std;
using namespace Grid;
using namespace Grid::QCD;
#define parallel_for PARALLEL_FOR_LOOP for
int main (int argc, char ** argv)
{
Grid_init(&argc,&argv);
std::vector<int> latt_size = GridDefaultLatt();
std::vector<int> simd_layout = GridDefaultSimd(Nd,vComplex::Nsimd());
std::vector<int> mpi_layout = GridDefaultMpi();
GridCartesian Grid(latt_size,simd_layout,mpi_layout);
GridRedBlackCartesian RBGrid(latt_size,simd_layout,mpi_layout);
int threads = GridThread::GetThreads();
std::cout<<GridLogMessage << "Grid is setup to use "<<threads<<" threads"<<std::endl;
std::vector<int> seeds({1,2,3,4});
GridParallelRNG pRNG(&Grid);
pRNG.SeedRandomDevice();
LatticeGaugeField U(&Grid);
SU3::HotConfiguration(pRNG,U);
double beta = 1.0;
double c1 = 0.331;
PlaqPlusRectangleActionR Action(beta,c1);
// WilsonGaugeActionR Action(beta);
ComplexD S = Action.S(U);
// get the deriv of phidag MdagM phi with respect to "U"
LatticeGaugeField UdSdU(&Grid);
Action.deriv(U,UdSdU);
////////////////////////////////////
// Modify the gauge field a little
////////////////////////////////////
RealD dt = 0.0001;
LatticeColourMatrix mommu(&Grid);
LatticeColourMatrix forcemu(&Grid);
LatticeGaugeField mom(&Grid);
LatticeGaugeField Uprime(&Grid);
for(int mu=0;mu<Nd;mu++){
SU3::GaussianLieAlgebraMatrix(pRNG, mommu); // Traceless antihermitian momentum; gaussian in lie alg
PokeIndex<LorentzIndex>(mom,mommu,mu);
// fourth order exponential approx
parallel_for(auto i=mom.begin();i<mom.end();i++){ // exp(pmu dt) * Umu
Uprime[i](mu) = U[i](mu) + mom[i](mu)*U[i](mu)*dt ;
}
}
ComplexD Sprime = Action.S(Uprime);
//////////////////////////////////////////////
// Use derivative to estimate dS
//////////////////////////////////////////////
LatticeComplex dS(&Grid); dS = zero;
for(int mu=0;mu<Nd;mu++){
auto UdSdUmu = PeekIndex<LorentzIndex>(UdSdU,mu);
mommu = PeekIndex<LorentzIndex>(mom,mu);
// Update gauge action density
// U = exp(p dt) U
// dU/dt = p U
// so dSdt = trace( dUdt dSdU) = trace( p UdSdUmu )
dS = dS - trace(mommu*UdSdUmu)*dt*2.0;
}
Complex dSpred = sum(dS);
std::cout << GridLogMessage << " S "<<S<<std::endl;
std::cout << GridLogMessage << " Sprime "<<Sprime<<std::endl;
std::cout << GridLogMessage << "dS "<<Sprime-S<<std::endl;
std::cout << GridLogMessage << "pred dS "<< dSpred <<std::endl;
std::cout<< GridLogMessage << "Done" <<std::endl;
Grid_finalize();
}

6
tests/qdpxx/Make.inc
View File

@ -1,5 +1,9 @@
bin_PROGRAMS = Test_qdpxx_munprec
bin_PROGRAMS = Test_qdpxx_loops_staples Test_qdpxx_munprec
Test_qdpxx_loops_staples_SOURCES=Test_qdpxx_loops_staples.cc
Test_qdpxx_loops_staples_LDADD=-lGrid
Test_qdpxx_munprec_SOURCES=Test_qdpxx_munprec.cc

648
tests/qdpxx/Test_qdpxx_loops_staples.cc Executable file
View File

@ -0,0 +1,648 @@
#include <Grid.h>
double calc_grid_p (Grid::QCD::LatticeGaugeField & lat);
double calc_chroma_p (Grid::QCD::LatticeGaugeField & lat);
double calc_grid_r (Grid::QCD::LatticeGaugeField & lat);
double calc_grid_IW (Grid::QCD::LatticeGaugeField & lat);
double calc_grid_r_dir (Grid::QCD::LatticeGaugeField & lat);
double calc_chroma_r (Grid::QCD::LatticeGaugeField & lat);
double calc_chroma_IW (Grid::QCD::LatticeGaugeField & lat);
void check_grid_r_staple(Grid::QCD::LatticeGaugeField & Umu);
void check_grid_p_staple(Grid::QCD::LatticeGaugeField & Umu);
const double beta=2.6;
const double c1=-0.331;
#include <chroma.h>
#include <actions/ferm/invert/syssolver_linop_cg_array.h>
#include <actions/ferm/invert/syssolver_linop_aggregate.h>
namespace Chroma {
class ChromaWrapper {
public:
typedef multi1d<LatticeColorMatrix> U;
static void ImportGauge(Grid::QCD::LatticeGaugeField & gr,
QDP::multi1d<QDP::LatticeColorMatrix> & ch)
{
Grid::QCD::LorentzColourMatrix LCM;
Grid::Complex cc;
QDP::ColorMatrix cm;
QDP::Complex c;
std::vector<int> x(4);
QDP::multi1d<int> cx(4);
std::vector<int> gd= gr._grid->GlobalDimensions();
for (x[0]=0;x[0]<gd[0];x[0]++){
for (x[1]=0;x[1]<gd[1];x[1]++){
for (x[2]=0;x[2]<gd[2];x[2]++){
for (x[3]=0;x[3]<gd[3];x[3]++){
cx[0] = x[0];
cx[1] = x[1];
cx[2] = x[2];
cx[3] = x[3];
Grid::peekSite(LCM,gr,x);
for(int mu=0;mu<4;mu++){
for(int i=0;i<3;i++){
for(int j=0;j<3;j++){
cc = LCM(mu)()(i,j);
c = QDP::cmplx(QDP::Real(real(cc)),QDP::Real(imag(cc)));
QDP::pokeColor(cm,c,i,j);
}}
QDP::pokeSite(ch[mu],cm,cx);
}
}}}}
}
static Chroma::Handle< Chroma::LinearGaugeAction > GetRBCAction ( U &u )
{
Chroma::CreatePeriodicGaugeState<U,U> CPGS;
Chroma::Handle< Chroma::CreateGaugeState<U,U> > cgs (new Chroma::CreatePeriodicGaugeState<U,U>);
Chroma::Handle< Chroma::LinearGaugeAction> action = new Chroma::RBCGaugeAct(cgs,beta,c1);
return action;
}
static Chroma::Handle< Chroma::LinearGaugeAction > GetRectangle ( U &u )
{
Chroma::CreatePeriodicGaugeState<U,U> CPGS;
Chroma::Handle< Chroma::CreateGaugeState<U,U> > cgs (new Chroma::CreatePeriodicGaugeState<U,U>);
Chroma::Handle< Chroma::LinearGaugeAction> action = new Chroma::RectGaugeAct(cgs, Real(c1));
return action;
}
static Chroma::Handle< Chroma::LinearGaugeAction > GetPlaquette ( U &u )
{
Chroma::CreatePeriodicGaugeState<U,U> CPGS;
Chroma::Handle< Chroma::CreateGaugeState<U,U> > cgs (new Chroma::CreatePeriodicGaugeState<U,U>);
Chroma::Handle< Chroma::LinearGaugeAction> action = new Chroma::WilsonGaugeAct(cgs, Real(beta));
return action;
}
};
}
int main (int argc,char **argv )
{
/********************************************************
* Setup QDP
*********************************************************/
Chroma::initialize(&argc,&argv);
Chroma::WilsonTypeFermActs4DEnv::registerAll();
/********************************************************
* Setup Grid
*********************************************************/
Grid::Grid_init(&argc,&argv);
Grid::GridCartesian * UGrid = Grid::QCD::SpaceTimeGrid::makeFourDimGrid(Grid::GridDefaultLatt(),
Grid::GridDefaultSimd(Grid::QCD::Nd,Grid::vComplex::Nsimd()),
Grid::GridDefaultMpi());
std::vector<int> gd = UGrid->GlobalDimensions();
QDP::multi1d<int> nrow(QDP::Nd);
for(int mu=0;mu<4;mu++) nrow[mu] = gd[mu];
QDP::Layout::setLattSize(nrow);
QDP::Layout::create();
Grid::QCD::LatticeGaugeField lat(UGrid);
double s_grid = calc_grid_p (lat);
double s_chroma = calc_chroma_p(lat);
// Match conventions
double vol = UGrid->gSites();
s_chroma+= beta * QDP::Nd * (QDP::Nd-1)*vol *0.5;
std::cout << " Chroma/Grid plaquette = " <<s_chroma<<" "<<s_grid<<std::endl;
std::cout << " Chroma avg plaquette = " <<1.0 - s_chroma / vol/ 6 / beta <<std::endl;
s_grid = calc_grid_r (lat);
s_chroma = calc_chroma_r(lat);
std::cout << std::setprecision(10);
std::cout<< " bare s_chroma "<<s_chroma<<std::endl;
s_chroma+= c1 * 12.0 *vol ;
std::cout<< " adjusted s_chroma "<<s_chroma<<std::endl;
std::cout<< " adjust "<< c1 * 12.0 *vol <<std::endl;
std::cout << " Chroma/Grid rectangle = " <<s_chroma<<" "<<s_grid<<std::endl;
std::cout << " Chroma avg rectangle = " <<1.0 - s_chroma / vol/ 12.0 / c1 <<std::endl;
check_grid_r_staple(lat);
check_grid_p_staple(lat);
calc_grid_r_dir(lat);
// Iwasaki case
std::cout << "Checking Iwasaki action (Grid/Chroma) " << std::endl;
s_grid = calc_grid_IW(lat);
s_chroma = calc_chroma_IW(lat);
// Adjust for the missing "1" in chroma's action.
s_chroma+= beta * (1.0-8*c1) * QDP::Nd * (QDP::Nd-1)*vol *0.5;
s_chroma+= beta * c1 * 12.0 *vol ;
std::cout << "Iwasaki action (Grid/Chroma) = " << s_grid<< " " << s_chroma <<std::endl;
Chroma::finalize();
}
double calc_chroma_IW(Grid::QCD::LatticeGaugeField & lat)
{
typedef QDP::multi1d<QDP::LatticeColorMatrix> U;
QDP::multi1d<QDP::LatticeColorMatrix> u(4);
Chroma::ChromaWrapper::ImportGauge(lat,u) ;
for(int mu=0;mu<4;mu++){
std::cout <<"Imported Gauge norm ["<<mu<<"] "<< QDP::norm2(u[mu])<<std::endl;
}
auto action =Chroma::ChromaWrapper::GetRBCAction(u);
Chroma::Handle<GaugeState<U,U> > gs(action->getCreateState()(u));
Real act = action->S(gs);
double s = toDouble(act);
return s;
}
double calc_chroma_r(Grid::QCD::LatticeGaugeField & lat)
{
typedef QDP::multi1d<QDP::LatticeColorMatrix> U;
QDP::multi1d<QDP::LatticeColorMatrix> u(4);
Chroma::ChromaWrapper::ImportGauge(lat,u) ;
for(int mu=0;mu<4;mu++){
std::cout <<"Imported Gauge norm ["<<mu<<"] "<< QDP::norm2(u[mu])<<std::endl;
}
auto action =Chroma::ChromaWrapper::GetRectangle(u);
Chroma::Handle<GaugeState<U,U> > gs(action->getCreateState()(u));
Real act = action->S(gs);
double s = toDouble(act);
return s;
}
// Conventions matching:
//
// Chroma:
//
// w_plaq is defined in MesPlq as
// w_plaq =( 2/(V*Nd*(Nd-1)*Nc)) * Sum Re Tr Plaq
//
// S = -(coeff/(Nc) Sum Re Tr Plaq
//
// S = -coeff * (V*Nd*(Nd-1)/2) w_plaq
// = -coeff * (V*Nd*(Nd-1)/2)*(2/(V*Nd*(Nd-1)*Nc))* Sum Re Tr Plaq
// = -coeff * (1/(Nc)) * Sum Re Tr Plaq
//
// Grid: has 1-plaq as usual to define Fmunu
//
// action = beta*(1.0 -plaq)*(Nd*(Nd-1.0))*vol*0.5;
// action = beta * Nd*Nd-1*vol*0.5 - beta * Nd*Nd-1*vol*0.5*plaq
//
// plaq == sumplaq * 2/(V*Nd*(Nd-1)*Nc)
double calc_chroma_p(Grid::QCD::LatticeGaugeField & lat)
{
typedef QDP::multi1d<QDP::LatticeColorMatrix> U;
QDP::multi1d<QDP::LatticeColorMatrix> u(4);
Chroma::ChromaWrapper::ImportGauge(lat,u) ;
for(int mu=0;mu<4;mu++){
std::cout <<"Imported Gauge norm ["<<mu<<"] "<< QDP::norm2(u[mu])<<std::endl;
}
auto action =Chroma::ChromaWrapper::GetPlaquette(u);
Chroma::Handle<GaugeState<U,U> > gs(action->getCreateState()(u));
Real act = action->S(gs);
double s = toDouble(act);
return s;
}
double calc_grid_p(Grid::QCD::LatticeGaugeField & Umu)
{
std::vector<int> seeds4({1,2,3,4});
Grid::GridCartesian * UGrid = (Grid::GridCartesian *) Umu._grid;
Grid::GridParallelRNG RNG4(UGrid); RNG4.SeedFixedIntegers(seeds4);
Grid::QCD::SU3::HotConfiguration(RNG4,Umu);
Grid::QCD::LatticeColourMatrix tmp(UGrid);
tmp = Grid::zero;
Grid::QCD::PokeIndex<Grid::QCD::LorentzIndex>(Umu,tmp,2);
Grid::QCD::PokeIndex<Grid::QCD::LorentzIndex>(Umu,tmp,3);
Grid::QCD::WilsonGaugeActionR Wilson(beta); // Just take beta = 1.0
return Wilson.S(Umu);
}
double calc_grid_r(Grid::QCD::LatticeGaugeField & Umu)
{
Grid::GridCartesian * UGrid = (Grid::GridCartesian *) Umu._grid;
Grid::QCD::PlaqPlusRectangleActionR Wilson(0.0,c1); // Just take beta = 0.0
return Wilson.S(Umu);
}
double calc_grid_IW(Grid::QCD::LatticeGaugeField & Umu)
{
Grid::GridCartesian * UGrid = (Grid::GridCartesian *) Umu._grid;
Grid::QCD::IwasakiGaugeActionR Iwasaki(beta);
return Iwasaki.S(Umu);
}
double calc_grid_r_dir(Grid::QCD::LatticeGaugeField & Umu)
{
Grid::GridCartesian * UGrid = (Grid::GridCartesian *) Umu._grid;
std::vector<Grid::QCD::LatticeColourMatrix> U(4,UGrid);
for(int mu=0;mu<Nd;mu++){
U[mu] = Grid::PeekIndex<Grid::QCD::LorentzIndex>(Umu,mu);
}
Grid::QCD::LatticeComplex rect(UGrid);
Grid::QCD::TComplex trect;
Grid::QCD::Complex crect;
Grid::RealD vol = UGrid->gSites();
for(int mu=0;mu<Grid::QCD::Nd;mu++){
for(int nu=0;nu<Grid::QCD::Nd;nu++){
if ( mu!=nu ) {
Grid::QCD::WilsonLoops<Grid::QCD::LatticeGaugeField>::traceDirRectangle(rect,U,mu,nu);
trect = Grid::sum(rect);
crect = Grid::TensorRemove(trect);
std::cout<< "mu/nu = "<<mu<<"/"<<nu<<" ; rect = "<<crect/vol/2.0/3.0<<std::endl;
Grid::GridStopWatch Peter;
Grid::GridStopWatch Azusa;
// Staple test
Peter.Start();
{
Grid::QCD::LatticeColourMatrix Stap(UGrid);
Grid::QCD::LatticeComplex SumTrStap(UGrid);
Grid::QCD::LatticeComplex TrStap(UGrid);
/*
* Make staple for loops centered at coor of link ; this one is ok. // |
*/
// __ ___
// | __ |
Stap =
Grid::Cshift(Grid::QCD::CovShiftForward (U[mu],mu,
Grid::QCD::CovShiftForward (U[nu],nu,
Grid::QCD::CovShiftBackward(U[mu],mu,
Grid::QCD::CovShiftBackward(U[mu],mu,
Grid::Cshift(adj(U[nu]),nu,-1))))) , mu, 1);
TrStap = Grid::trace (U[mu]*Stap);
SumTrStap = TrStap;
trect = Grid::sum(TrStap);
crect = Grid::TensorRemove(trect);
std::cout<< "mu/nu inline trace 2x1 code = "<<mu<<"/"<<nu<<" ; rect = "<<crect/vol/1.0/3.0<<std::endl;
// __
// |__ __ |
Stap = Grid::Cshift(Grid::QCD::CovShiftForward (U[mu],mu,
Grid::QCD::CovShiftBackward(U[nu],nu,
Grid::QCD::CovShiftBackward(U[mu],mu,
Grid::QCD::CovShiftBackward(U[mu],mu, U[nu])))) , mu, 1);
TrStap = Grid::trace (U[mu]*Stap);
trect = Grid::sum(TrStap);
crect = Grid::TensorRemove(trect);
std::cout<< "mu/nu inline trace 2x1 code = "<<mu<<"/"<<nu<<" ; rect = "<<crect/vol/1.0/3.0<<std::endl;
// __
// |__ __ |
Stap = Grid::Cshift(Grid::QCD::CovShiftBackward(U[nu],nu,
Grid::QCD::CovShiftBackward(U[mu],mu,
Grid::QCD::CovShiftBackward(U[mu],mu,
Grid::QCD::CovShiftForward(U[nu],nu,U[mu])))) , mu, 1);
TrStap = Grid::trace (U[mu]*Stap);
trect = Grid::sum(TrStap);
crect = Grid::TensorRemove(trect);
std::cout<< "mu/nu inline trace 2x1 code = "<<mu<<"/"<<nu<<" ; rect = "<<crect/vol/1.0/3.0<<std::endl;
// __ ___
// |__ |
Stap = Grid::Cshift(Grid::QCD::CovShiftForward (U[nu],nu,
Grid::QCD::CovShiftBackward(U[mu],mu,
Grid::QCD::CovShiftBackward(U[mu],mu,
Grid::QCD::CovShiftBackward(U[nu],nu,U[mu])))) , mu, 1);
TrStap = Grid::trace (U[mu]*Stap);
trect = Grid::sum(TrStap);
crect = Grid::TensorRemove(trect);
std::cout<< "mu/nu inline trace 2x1 code = "<<mu<<"/"<<nu<<" ; rect = "<<crect/vol/1.0/3.0<<std::endl;
// --
// | |
//
// | |
//
/*
* Make staple for loops centered at coor of link ; this one is ok. // |
*/
// Stap =
// Grid::Cshift(Grid::QCD::CovShiftForward(U[nu],nu,U[nu]),mu,1)* // ->||
// Grid::adj(Grid::QCD::CovShiftForward(U[nu],nu,Grid::QCD::CovShiftForward(U[nu],nu,U[mu]))) ;
Stap = Grid::Cshift(Grid::QCD::CovShiftForward(U[nu],nu,
Grid::QCD::CovShiftForward(U[nu],nu,
Grid::QCD::CovShiftBackward(U[mu],mu,
Grid::QCD::CovShiftBackward(U[nu],nu, Grid::Cshift(adj(U[nu]),nu,-1))))) , mu, 1);
TrStap = Grid::trace (U[mu]*Stap);
SumTrStap += TrStap;
trect = Grid::sum(TrStap);
crect = Grid::TensorRemove(trect);
std::cout<< "mu/nu inline trace 1x2 code = "<<mu<<"/"<<nu<<" ; rect = "<<crect/vol/1.0/3.0<<std::endl;
//
// | |
//
// | |
// --
Stap = Grid::Cshift(Grid::QCD::CovShiftBackward(U[nu],nu,
Grid::QCD::CovShiftBackward(U[nu],nu,
Grid::QCD::CovShiftBackward(U[mu],mu,
Grid::QCD::CovShiftForward (U[nu],nu,U[nu])))) , mu, 1);
TrStap = Grid::trace (U[mu]*Stap);
trect = Grid::sum(TrStap);
crect = Grid::TensorRemove(trect);
std::cout<< "mu/nu inline trace 1x2 code = "<<mu<<"/"<<nu<<" ; rect = "<<crect/vol/1.0/3.0<<std::endl;
trect = Grid::sum(SumTrStap);
crect = Grid::TensorRemove(trect);
std::cout<< "mu/nu inline trace 2x1+1x2 code = "<<mu<<"/"<<nu<<" ; rect = "<<crect/vol/2.0/3.0<<std::endl;
}
Peter.Stop();
Azusa.Start();
{
Grid::QCD::LatticeComplex RectPlaq_d(UGrid);
Grid::QCD::LatticeColourMatrix ds_U(UGrid);
Grid::QCD::LatticeColourMatrix left_2(UGrid);
Grid::QCD::LatticeColourMatrix upper_l(UGrid);
Grid::QCD::LatticeColourMatrix upper_staple(UGrid);
Grid::QCD::LatticeColourMatrix down_l(UGrid);
Grid::QCD::LatticeColourMatrix down_staple(UGrid);
Grid::QCD::LatticeColourMatrix tmp(UGrid);
// 2 (mu)x1(nu)
left_2= Grid::QCD::CovShiftForward(U[mu],mu,U[mu]); // Umu(x) Umu(x+mu)
tmp=Grid::Cshift(U[nu],mu,2); // Unu(x+2mu)
upper_l= Grid::QCD::CovShiftForward(tmp,nu,Grid::adj(left_2)); // Unu(x+2mu) Umu^dag(x+mu+nu) Umu^dag(x+nu)
// __ __
// = |
// Unu(x-2mu) Umudag(x-mu+nu) Umudag(x+nu) Unudag(x)
// __ __
// upper_staple = | |
// v
upper_staple= upper_l*adj(U[nu]); // Unu(x+2mu) Umu^dag(x+mu+nu) Umu^dag(x+nu) Unu^dag(x)
//
// down_staple = |__ __|
//
tmp = adj(left_2*tmp)*U[nu];
down_staple= Grid::Cshift(tmp,nu,-1); // Unu^dag((x+2mu-nu) Umu^dag(x+mu-nu) Umu^dag(x-nu) Unu(x-nu)
// __ __
// | __|
//
tmp=Grid::Cshift(U[mu],mu,1);
ds_U=tmp*(upper_staple); // Umu(x+mu) Unu(x+2mu) Umu^dag(x+mu+nu) Umu^dag(x+nu) Unu^dag(x)
RectPlaq_d = Grid::trace(U[mu]*ds_U);
trect = Grid::sum(RectPlaq_d);
crect = Grid::TensorRemove(trect);
std::cout<< "mu/nu inline AZUSA trace 2x1 code = "<<mu<<"/"<<nu<<" ; rect = "<<crect/vol/1.0/3.0<<std::endl;
// __ __
// |__ |
//
tmp=upper_staple*U[mu];
ds_U= Grid::Cshift(tmp,mu,-1);
RectPlaq_d = Grid::trace(U[mu]*ds_U);
trect = Grid::sum(RectPlaq_d);
crect = Grid::TensorRemove(trect);
std::cout<< "mu/nu inline AZUSA trace 2x1 code = "<<mu<<"/"<<nu<<" ; rect = "<<crect/vol/1.0/3.0<<std::endl;
// __
// |__ __ |
//
tmp=Grid::Cshift(U[mu],mu,1);
ds_U=tmp*(down_staple); // Umu(x+mu) Unu^dag((x+2mu-nu) Umu^dag(x+mu-nu) Umu^dag(x-nu) Unu(x-nu)
RectPlaq_d = Grid::trace(U[mu]*ds_U);
trect = Grid::sum(RectPlaq_d);
crect = Grid::TensorRemove(trect);
std::cout<< "mu/nu inline AZUSA trace 2x1 code = "<<mu<<"/"<<nu<<" ; rect = "<<crect/vol/1.0/3.0<<std::endl;
// __
// |__ __ |
//
tmp = down_staple*U[mu];
ds_U=Grid::Cshift(tmp,mu,-1);
RectPlaq_d = Grid::trace(U[mu]*ds_U);
trect = Grid::sum(RectPlaq_d);
crect = Grid::TensorRemove(trect);
std::cout<< "mu/nu inline AZUSA trace 2x1 code = "<<mu<<"/"<<nu<<" ; rect = "<<crect/vol/1.0/3.0<<std::endl;
// 1(mu) x 2 (nu) ** this was ok
// _
// | |
// | |
Grid::QCD::LatticeColourMatrix up2= Grid::QCD::CovShiftForward(U[nu],nu,U[nu]);
upper_l= Grid::QCD::CovShiftForward(Grid::Cshift(up2,mu,1),nu,Grid::Cshift(adj(U[mu]),nu,1));
ds_U= upper_l*Grid::adj(up2);
RectPlaq_d = Grid::trace(U[mu]*ds_U);
trect = Grid::sum(RectPlaq_d);
crect = Grid::TensorRemove(trect);
std::cout<< "mu/nu inline AZUSA trace 1x2 code = "<<mu<<"/"<<nu<<" ; rect = "<<crect/vol/1.0/3.0<<std::endl;
// 1(mu) x 2 (nu) ** this was ok
//
// | |
// |_|
/*
|
V
downer_l= |
(x)<----V
*/
down_l= Grid::adj(Grid::QCD::CovShiftForward(U[mu],mu,up2)); //downer_l
/*
^ |
down_staple = | V
^ |
| V
(x)<----
down_staple= upper*upper_l;
*/
tmp= down_l*up2;
ds_U= Grid::Cshift(tmp,nu,-2);
RectPlaq_d = Grid::trace(U[mu]*ds_U);
trect = Grid::sum(RectPlaq_d);
crect = Grid::TensorRemove(trect);
std::cout<< "mu/nu inline AZUSA trace 1x2 code = "<<mu<<"/"<<nu<<" ; rect = "<<crect/vol/1.0/3.0<<std::endl;
}
Azusa.Stop();
std::chrono::milliseconds pt = Peter.Elapsed();
std::chrono::milliseconds az = Azusa.Elapsed();
std::cout<< "Times ";
std::cout<<pt.count();
std::cout<<" A ";
std::cout<<az.count();
std::cout<<std::endl;
}
}}
Grid::QCD::PlaqPlusRectangleActionR Wilson(0.0,c1); // Just take beta = 0.0
return Wilson.S(Umu);
};
void check_grid_r_staple(Grid::QCD::LatticeGaugeField & Umu)
{
std::vector<int> seeds4({1,2,3,4});
Grid::GridCartesian * UGrid = (Grid::GridCartesian *) Umu._grid;
Grid::QCD::PlaqPlusRectangleActionR Wilson(0.0,c1); // Just take beta = 0.0
Grid::QCD::LatticeColourMatrix staple(UGrid);
Grid::QCD::LatticeColourMatrix link(UGrid);
Grid::QCD::LatticeComplex Traced(UGrid);
Grid::Complex Rplaq(0.0);
for(int mu=0;mu<Nd;mu++){
Grid::RealD vol = UGrid->gSites();
// for mu, nu!=mu => 12
// 6 loops contribute to staple for each orientation.
// Nc=3.
// Vol as for each site
Grid::RealD RectScale(1.0/vol/12.0/6.0/3.0);
Grid::QCD::WilsonLoops<Grid::QCD::LatticeGaugeField>::RectStaple(staple,Umu,mu);
link = Grid::QCD::PeekIndex<Grid::QCD::LorentzIndex>(Umu,mu);
Traced = Grid::trace( link*staple) * RectScale;
Grid::QCD::TComplex Tp = Grid::sum(Traced);
Grid::Complex p = Grid::TensorRemove(Tp);
std::cout<< "Rect from RectStaple "<<p<<std::endl;
Rplaq = Rplaq+ p;
}
std::cout<< "Rect from RectStaple "<<Rplaq<<std::endl;
}
void check_grid_p_staple(Grid::QCD::LatticeGaugeField & Umu)
{
std::vector<int> seeds4({1,2,3,4});
Grid::GridCartesian * UGrid = (Grid::GridCartesian *) Umu._grid;
Grid::QCD::PlaqPlusRectangleActionR Wilson(1.0,0.0); // Just take c1 = 0.0
Grid::QCD::LatticeColourMatrix staple(UGrid);
Grid::QCD::LatticeColourMatrix link(UGrid);
Grid::QCD::LatticeComplex Traced(UGrid);
Grid::Complex plaq(0.0);
for(int mu=0;mu<Nd;mu++){
Grid::RealD vol = UGrid->gSites();
// for mu, nu!=mu => 12
// 2 loops contribute to staple for each orientation.
// Nc=3.
// Vol as for each site
Grid::RealD Scale(1.0/vol/12.0/2.0/3.0);
Grid::QCD::WilsonLoops<Grid::QCD::LatticeGaugeField>::Staple(staple,Umu,mu);
link = Grid::QCD::PeekIndex<Grid::QCD::LorentzIndex>(Umu,mu);
Traced = Grid::trace( link*staple) * Scale;
Grid::QCD::TComplex Tp = Grid::sum(Traced);
Grid::Complex p = Grid::TensorRemove(Tp);
std::cout<< "Plaq from PlaqStaple "<<p<<std::endl;
plaq = plaq+ p;
}
std::cout<< "Plaq from PlaqStaple "<<plaq<<std::endl;
}

7
tests/qdpxx/Test_qdpxx_munprec.cc
View File

@ -425,13 +425,6 @@ void calc_chroma(ChromaAction action,Grid::QCD::LatticeGaugeField & lat, Grid::Q
// Chroma::HotSt(u);
Chroma::ChromaWrapper::ImportGauge(lat,u) ;
int lx = QDP::Layout::subgridLattSize()[0];
int ly = QDP::Layout::subgridLattSize()[1];
int lz = QDP::Layout::subgridLattSize()[2];
int lt = QDP::Layout::subgridLattSize()[3];
QDP::multi1d<int> procs = QDP::Layout::logicalSize();
QDP::multi1d<QDP::LatticeFermion> check(Ls);
QDP::multi1d<QDP::LatticeFermion> result(Ls);
QDP::multi1d<QDP::LatticeFermion> psi(Ls);