mirror of
https://github.com/paboyle/Grid.git
synced 2025-06-14 05:07:05 +01:00
Merge branch 'develop' of https://github.com/paboyle/Grid into develop
This commit is contained in:
Peter Boyle
@ -75,8 +75,8 @@ int main (int argc, char ** argv)
|
||||
RealD M5=1.8;
|
||||
|
||||
{
|
||||
OverlapWilsonContFracTanhFermionR Dcf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,1.0);
|
||||
HermitianLinearOperator<OverlapWilsonContFracTanhFermionR,LatticeFermion> HermIndefOp(Dcf);
|
||||
OverlapWilsonContFracTanhFermionD Dcf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,1.0);
|
||||
HermitianLinearOperator<OverlapWilsonContFracTanhFermionD,LatticeFermion> HermIndefOp(Dcf);
|
||||
|
||||
HermIndefOp.Op(src,ref);
|
||||
HermIndefOp.OpDiag(src,result);
|
||||
@ -92,8 +92,8 @@ int main (int argc, char ** argv)
|
||||
}
|
||||
|
||||
{
|
||||
OverlapWilsonPartialFractionTanhFermionR Dpf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,1.0);
|
||||
HermitianLinearOperator<OverlapWilsonPartialFractionTanhFermionR,LatticeFermion> HermIndefOp(Dpf);
|
||||
OverlapWilsonPartialFractionTanhFermionD Dpf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,1.0);
|
||||
HermitianLinearOperator<OverlapWilsonPartialFractionTanhFermionD,LatticeFermion> HermIndefOp(Dpf);
|
||||
|
||||
HermIndefOp.Op(src,ref);
|
||||
HermIndefOp.OpDiag(src,result);
|
||||
|
||||
@ -140,14 +140,14 @@ int main (int argc, char ** argv)
|
||||
|
||||
// RealD mass=0.1;
|
||||
// RealD M5=1.8;
|
||||
// DomainWallFermionR Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5);
|
||||
// DomainWallFermionD Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5);
|
||||
|
||||
// LatticeFermion src_o(FrbGrid);
|
||||
// LatticeFermion result_o(FrbGrid);
|
||||
// pickCheckerboard(Odd,src_o,src);
|
||||
// result_o=Zero();
|
||||
|
||||
// SchurDiagMooeeOperator<DomainWallFermionR,LatticeFermion> HermOpEO(Ddwf);
|
||||
// SchurDiagMooeeOperator<DomainWallFermionD,LatticeFermion> HermOpEO(Ddwf);
|
||||
// ConjugateGradient<LatticeFermion> CG(1.0e-8,10000);
|
||||
// CG(HermOpEO,src_o,result_o);
|
||||
|
||||
|
||||
@ -53,7 +53,7 @@ static int readInt(int* argc, char*** argv, std::string&& option, int defaultVal
|
||||
|
||||
static float readFloat(int* argc, char*** argv, std::string&& option, float defaultValue) {
|
||||
std::string arg;
|
||||
float ret = defaultValue;
|
||||
double ret = defaultValue;
|
||||
if(checkPresent(argc, argv, option)) {
|
||||
arg = getContent(argc, argv, option);
|
||||
GridCmdOptionFloat(arg, ret);
|
||||
|
||||
@ -76,20 +76,20 @@ int main (int argc, char ** argv)
|
||||
RealD M5 =1.8;
|
||||
|
||||
std::cout<<GridLogMessage <<"OverlapWilsonContFracTanhFermion test"<<std::endl;
|
||||
OverlapWilsonContFracTanhFermionR Dcf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,1.0);
|
||||
TestWhat<OverlapWilsonContFracTanhFermionR>(Dcf,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5);
|
||||
OverlapWilsonContFracTanhFermionD Dcf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,1.0);
|
||||
TestWhat<OverlapWilsonContFracTanhFermionD>(Dcf,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5);
|
||||
|
||||
std::cout<<GridLogMessage <<"OverlapWilsonContFracZolotarevFermion test"<<std::endl;
|
||||
OverlapWilsonContFracZolotarevFermionR Dcfz(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,0.1,6.0);
|
||||
TestWhat<OverlapWilsonContFracZolotarevFermionR>(Dcfz,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5);
|
||||
OverlapWilsonContFracZolotarevFermionD Dcfz(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,0.1,6.0);
|
||||
TestWhat<OverlapWilsonContFracZolotarevFermionD>(Dcfz,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5);
|
||||
|
||||
std::cout<<GridLogMessage <<"OverlapWilsonPartialFractionTanhFermion test"<<std::endl;
|
||||
OverlapWilsonPartialFractionTanhFermionR Dpf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,1.0);
|
||||
TestWhat<OverlapWilsonPartialFractionTanhFermionR>(Dpf,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5);
|
||||
OverlapWilsonPartialFractionTanhFermionD Dpf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,1.0);
|
||||
TestWhat<OverlapWilsonPartialFractionTanhFermionD>(Dpf,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5);
|
||||
|
||||
std::cout<<GridLogMessage <<"OverlapWilsonPartialFractionZolotarevFermion test"<<std::endl;
|
||||
OverlapWilsonPartialFractionZolotarevFermionR Dpfz(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,0.1,6.0);
|
||||
TestWhat<OverlapWilsonPartialFractionZolotarevFermionR>(Dpfz,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5);
|
||||
OverlapWilsonPartialFractionZolotarevFermionD Dpfz(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,0.1,6.0);
|
||||
TestWhat<OverlapWilsonPartialFractionZolotarevFermionD>(Dpfz,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5);
|
||||
|
||||
Grid_finalize();
|
||||
}
|
||||
|
||||
@ -90,7 +90,7 @@ int main (int argc, char ** argv)
|
||||
RealD shift = 0.1234;
|
||||
RealD M5 = 1.8;
|
||||
int pm = 1;
|
||||
DomainWallEOFAFermionR Ddwf(Umu, *FGrid, *FrbGrid, *UGrid, *UrbGrid, mq1, mq2, mq3, shift, pm, M5);
|
||||
DomainWallEOFAFermionD Ddwf(Umu, *FGrid, *FrbGrid, *UGrid, *UrbGrid, mq1, mq2, mq3, shift, pm, M5);
|
||||
|
||||
LatticeFermion src_e (FrbGrid);
|
||||
LatticeFermion src_o (FrbGrid);
|
||||
@ -216,7 +216,7 @@ int main (int argc, char ** argv)
|
||||
pickCheckerboard(Even, phi_e, phi);
|
||||
pickCheckerboard(Odd , phi_o, phi);
|
||||
|
||||
SchurDiagMooeeOperator<DomainWallEOFAFermionR,LatticeFermion> HermOpEO(Ddwf);
|
||||
SchurDiagMooeeOperator<DomainWallEOFAFermionD,LatticeFermion> HermOpEO(Ddwf);
|
||||
HermOpEO.MpcDagMpc(chi_e, dchi_e);
|
||||
HermOpEO.MpcDagMpc(chi_o, dchi_o);
|
||||
|
||||
|
||||
@ -86,7 +86,7 @@ int main (int argc, char ** argv)
|
||||
|
||||
RealD mass=0.1;
|
||||
RealD M5 =1.8;
|
||||
DomainWallFermionR Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5);
|
||||
DomainWallFermionD Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5);
|
||||
|
||||
LatticeFermion src_e (FrbGrid);
|
||||
LatticeFermion src_o (FrbGrid);
|
||||
@ -213,7 +213,7 @@ int main (int argc, char ** argv)
|
||||
pickCheckerboard(Even,phi_e,phi);
|
||||
pickCheckerboard(Odd ,phi_o,phi);
|
||||
|
||||
SchurDiagMooeeOperator<DomainWallFermionR,LatticeFermion> HermOpEO(Ddwf);
|
||||
SchurDiagMooeeOperator<DomainWallFermionD,LatticeFermion> HermOpEO(Ddwf);
|
||||
HermOpEO.MpcDagMpc(chi_e,dchi_e);
|
||||
HermOpEO.MpcDagMpc(chi_o,dchi_o);
|
||||
|
||||
|
||||
@ -29,14 +29,10 @@ Author: Peter Boyle <paboyle@ph.ed.ac.uk>
|
||||
#include <Grid/Grid.h>
|
||||
|
||||
using namespace Grid;
|
||||
;
|
||||
|
||||
int main (int argc, char ** argv)
|
||||
{
|
||||
template<typename Gimpl>
|
||||
void run(double alpha, bool do_fft_gfix){
|
||||
std::vector<int> seeds({1,2,3,4});
|
||||
|
||||
Grid_init(&argc,&argv);
|
||||
|
||||
int threads = GridThread::GetThreads();
|
||||
|
||||
Coordinate latt_size = GridDefaultLatt();
|
||||
@ -55,10 +51,7 @@ int main (int argc, char ** argv)
|
||||
FFT theFFT(&GRID);
|
||||
|
||||
std::cout<<GridLogMessage << "Grid is setup to use "<<threads<<" threads"<<std::endl;
|
||||
|
||||
std::cout<< "*****************************************************************" <<std::endl;
|
||||
std::cout<< "* Testing we can gauge fix steep descent a RGT of Unit gauge *" <<std::endl;
|
||||
std::cout<< "*****************************************************************" <<std::endl;
|
||||
std::cout<<GridLogMessage << "Using alpha=" << alpha << std::endl;
|
||||
|
||||
// int coulomb_dir = -1;
|
||||
int coulomb_dir = Nd-1;
|
||||
@ -72,81 +65,167 @@ int main (int argc, char ** argv)
|
||||
LatticeColourMatrix xform1(&GRID); // Gauge xform
|
||||
LatticeColourMatrix xform2(&GRID); // Gauge xform
|
||||
LatticeColourMatrix xform3(&GRID); // Gauge xform
|
||||
|
||||
//#########################################################################################
|
||||
|
||||
std::cout<< "*********************************************************************************************************" <<std::endl;
|
||||
std::cout<< "* Testing steepest descent fixing to Landau gauge with randomly transformed unit gauge configuration *" <<std::endl;
|
||||
std::cout<< "*********************************************************************************************************" <<std::endl;
|
||||
|
||||
SU<Nc>::ColdConfiguration(pRNG,Umu); // Unit gauge
|
||||
Uorg=Umu;
|
||||
|
||||
Real init_plaq=WilsonLoops<Gimpl>::avgPlaquette(Umu);
|
||||
std::cout << " Initial plaquette "<< init_plaq << std::endl;
|
||||
|
||||
//Apply a random gauge transformation to the unit gauge config
|
||||
Urnd=Umu;
|
||||
SU<Nc>::RandomGaugeTransform<Gimpl>(pRNG,Urnd,g);
|
||||
|
||||
SU<Nc>::RandomGaugeTransform(pRNG,Urnd,g); // Unit gauge
|
||||
|
||||
Real plaq=WilsonLoops<PeriodicGimplR>::avgPlaquette(Umu);
|
||||
std::cout << " Initial plaquette "<<plaq << std::endl;
|
||||
|
||||
Real alpha=0.1;
|
||||
|
||||
//Gauge fix the randomly transformed field
|
||||
Umu = Urnd;
|
||||
FourierAcceleratedGaugeFixer<PeriodicGimplR>::SteepestDescentGaugeFix(Umu,xform1,alpha,10000,1.0e-12, 1.0e-12,false);
|
||||
FourierAcceleratedGaugeFixer<Gimpl>::SteepestDescentGaugeFix(Umu,xform1,alpha,10000,1.0e-12, 1.0e-12,false);
|
||||
|
||||
// Check the gauge xform matrices
|
||||
Utmp=Urnd;
|
||||
SU<Nc>::GaugeTransform(Utmp,xform1);
|
||||
SU<Nc>::GaugeTransform<Gimpl>(Utmp,xform1);
|
||||
Utmp = Utmp - Umu;
|
||||
std::cout << " Norm Difference of xformed gauge "<< norm2(Utmp) << std::endl;
|
||||
std::cout << " Check the output gauge transformation matrices applied to the original field produce the xformed field "<< norm2(Utmp) << " (expect 0)" << std::endl;
|
||||
|
||||
|
||||
plaq=WilsonLoops<PeriodicGimplR>::avgPlaquette(Umu);
|
||||
std::cout << " Final plaquette "<<plaq << std::endl;
|
||||
Real plaq=WilsonLoops<Gimpl>::avgPlaquette(Umu);
|
||||
std::cout << " Final plaquette "<<plaq << " diff " << plaq - init_plaq << " (expect 0)" << std::endl;
|
||||
|
||||
Uorg = Uorg - Umu;
|
||||
std::cout << " Norm Difference "<< norm2(Uorg) << std::endl;
|
||||
std::cout << " Norm "<< norm2(Umu) << std::endl;
|
||||
std::cout << " Norm difference between a unit gauge configuration and the gauge fixed configuration "<< norm2(Uorg) << " (expect 0)" << std::endl;
|
||||
std::cout << " Norm of gauge fixed configuration "<< norm2(Umu) << std::endl;
|
||||
|
||||
//#########################################################################################
|
||||
if(do_fft_gfix){
|
||||
std::cout<< "*************************************************************************************" <<std::endl;
|
||||
std::cout<< "* Testing Fourier accelerated fixing to Landau gauge with unit gauge configuration *" <<std::endl;
|
||||
std::cout<< "*************************************************************************************" <<std::endl;
|
||||
Umu=Urnd;
|
||||
FourierAcceleratedGaugeFixer<Gimpl>::SteepestDescentGaugeFix(Umu,xform2,alpha,10000,1.0e-12, 1.0e-12,true);
|
||||
|
||||
Utmp=Urnd;
|
||||
SU<Nc>::GaugeTransform<Gimpl>(Utmp,xform2);
|
||||
Utmp = Utmp - Umu;
|
||||
std::cout << " Check the output gauge transformation matrices applied to the original field produce the xformed field "<< norm2(Utmp) << " (expect 0)" << std::endl;
|
||||
|
||||
|
||||
std::cout<< "*****************************************************************" <<std::endl;
|
||||
std::cout<< "* Testing Fourier accelerated fixing *" <<std::endl;
|
||||
std::cout<< "*****************************************************************" <<std::endl;
|
||||
Umu=Urnd;
|
||||
FourierAcceleratedGaugeFixer<PeriodicGimplR>::SteepestDescentGaugeFix(Umu,xform2,alpha,10000,1.0e-12, 1.0e-12,true);
|
||||
plaq=WilsonLoops<Gimpl>::avgPlaquette(Umu);
|
||||
std::cout << " Final plaquette "<<plaq << " diff " << plaq - init_plaq << " (expect 0)" << std::endl;
|
||||
}
|
||||
//#########################################################################################
|
||||
|
||||
Utmp=Urnd;
|
||||
SU<Nc>::GaugeTransform(Utmp,xform2);
|
||||
Utmp = Utmp - Umu;
|
||||
std::cout << " Norm Difference of xformed gauge "<< norm2(Utmp) << std::endl;
|
||||
std::cout<< "******************************************************************************************" <<std::endl;
|
||||
std::cout<< "* Testing steepest descent fixing to Landau gauge with random configuration **" <<std::endl;
|
||||
std::cout<< "******************************************************************************************" <<std::endl;
|
||||
|
||||
SU<Nc>::HotConfiguration(pRNG,Umu);
|
||||
|
||||
plaq=WilsonLoops<PeriodicGimplR>::avgPlaquette(Umu);
|
||||
std::cout << " Final plaquette "<<plaq << std::endl;
|
||||
init_plaq=WilsonLoops<Gimpl>::avgPlaquette(Umu);
|
||||
std::cout << " Initial plaquette "<< init_plaq << std::endl;
|
||||
|
||||
std::cout<< "*****************************************************************" <<std::endl;
|
||||
std::cout<< "* Testing non-unit configuration *" <<std::endl;
|
||||
std::cout<< "*****************************************************************" <<std::endl;
|
||||
FourierAcceleratedGaugeFixer<Gimpl>::SteepestDescentGaugeFix(Umu,alpha,10000,1.0e-12, 1.0e-12,false);
|
||||
|
||||
SU<Nc>::HotConfiguration(pRNG,Umu); // Unit gauge
|
||||
plaq=WilsonLoops<Gimpl>::avgPlaquette(Umu);
|
||||
std::cout << " Final plaquette "<<plaq << " diff " << plaq - init_plaq << " (expect 0)" << std::endl;
|
||||
|
||||
plaq=WilsonLoops<PeriodicGimplR>::avgPlaquette(Umu);
|
||||
std::cout << " Initial plaquette "<<plaq << std::endl;
|
||||
//#########################################################################################
|
||||
if(do_fft_gfix){
|
||||
std::cout<< "******************************************************************************************" <<std::endl;
|
||||
std::cout<< "* Testing Fourier accelerated fixing to Landau gauge with random configuration **" <<std::endl;
|
||||
std::cout<< "******************************************************************************************" <<std::endl;
|
||||
|
||||
FourierAcceleratedGaugeFixer<PeriodicGimplR>::SteepestDescentGaugeFix(Umu,alpha,10000,1.0e-12, 1.0e-12,true);
|
||||
SU<Nc>::HotConfiguration(pRNG,Umu);
|
||||
|
||||
plaq=WilsonLoops<PeriodicGimplR>::avgPlaquette(Umu);
|
||||
std::cout << " Final plaquette "<<plaq << std::endl;
|
||||
init_plaq=WilsonLoops<Gimpl>::avgPlaquette(Umu);
|
||||
std::cout << " Initial plaquette "<< init_plaq << std::endl;
|
||||
|
||||
std::cout<< "*****************************************************************" <<std::endl;
|
||||
std::cout<< "* Testing Fourier accelerated fixing to coulomb gauge *" <<std::endl;
|
||||
std::cout<< "*****************************************************************" <<std::endl;
|
||||
FourierAcceleratedGaugeFixer<Gimpl>::SteepestDescentGaugeFix(Umu,alpha,10000,1.0e-12, 1.0e-12,true);
|
||||
|
||||
plaq=WilsonLoops<Gimpl>::avgPlaquette(Umu);
|
||||
std::cout << " Final plaquette "<<plaq << " diff " << plaq - init_plaq << " (expect 0)" << std::endl;
|
||||
}
|
||||
//#########################################################################################
|
||||
|
||||
std::cout<< "*******************************************************************************************" <<std::endl;
|
||||
std::cout<< "* Testing steepest descent fixing to coulomb gauge with random configuration *" <<std::endl;
|
||||
std::cout<< "*******************************************************************************************" <<std::endl;
|
||||
|
||||
Umu=Urnd;
|
||||
SU<Nc>::HotConfiguration(pRNG,Umu); // Unit gauge
|
||||
SU<Nc>::HotConfiguration(pRNG,Umu);
|
||||
|
||||
plaq=WilsonLoops<PeriodicGimplR>::avgPlaquette(Umu);
|
||||
std::cout << " Initial plaquette "<<plaq << std::endl;
|
||||
init_plaq=WilsonLoops<Gimpl>::avgPlaquette(Umu);
|
||||
std::cout << " Initial plaquette "<< init_plaq << std::endl;
|
||||
|
||||
FourierAcceleratedGaugeFixer<PeriodicGimplR>::SteepestDescentGaugeFix(Umu,xform3,alpha,10000,1.0e-12, 1.0e-12,true,coulomb_dir);
|
||||
FourierAcceleratedGaugeFixer<Gimpl>::SteepestDescentGaugeFix(Umu,xform3,alpha,10000,1.0e-12, 1.0e-12,false,coulomb_dir);
|
||||
|
||||
std::cout << Umu<<std::endl;
|
||||
plaq=WilsonLoops<Gimpl>::avgPlaquette(Umu);
|
||||
std::cout << " Final plaquette "<<plaq << " diff " << plaq - init_plaq << " (expect 0)" << std::endl;
|
||||
|
||||
plaq=WilsonLoops<PeriodicGimplR>::avgPlaquette(Umu);
|
||||
std::cout << " Final plaquette "<<plaq << std::endl;
|
||||
|
||||
//#########################################################################################
|
||||
if(do_fft_gfix){
|
||||
std::cout<< "*******************************************************************************************" <<std::endl;
|
||||
std::cout<< "* Testing Fourier accelerated fixing to coulomb gauge with random configuration *" <<std::endl;
|
||||
std::cout<< "*******************************************************************************************" <<std::endl;
|
||||
|
||||
Umu=Urnd;
|
||||
SU<Nc>::HotConfiguration(pRNG,Umu);
|
||||
|
||||
init_plaq=WilsonLoops<Gimpl>::avgPlaquette(Umu);
|
||||
std::cout << " Initial plaquette "<< init_plaq << std::endl;
|
||||
|
||||
FourierAcceleratedGaugeFixer<Gimpl>::SteepestDescentGaugeFix(Umu,xform3,alpha,10000,1.0e-12, 1.0e-12,true,coulomb_dir);
|
||||
|
||||
plaq=WilsonLoops<Gimpl>::avgPlaquette(Umu);
|
||||
std::cout << " Final plaquette "<<plaq << " diff " << plaq - init_plaq << " (expect 0)" << std::endl;
|
||||
}
|
||||
}
|
||||
|
||||
int main (int argc, char ** argv)
|
||||
{
|
||||
Grid_init(&argc,&argv);
|
||||
|
||||
double alpha=0.1; //step size
|
||||
std::string gimpl = "periodic";
|
||||
bool do_fft_gfix = true; //test fourier transformed gfix as well as steepest descent
|
||||
for(int i=1;i<argc;i++){
|
||||
std::string sarg(argv[i]);
|
||||
if(sarg == "--gimpl"){
|
||||
assert(i<argc-1 && "--gimpl option requires an argument");
|
||||
gimpl = argv[i+1];
|
||||
if(gimpl != "periodic" && gimpl != "conjugate")
|
||||
assert(0 && "Invalid gimpl");
|
||||
if(gimpl == "conjugate")
|
||||
alpha = 0.025; //default alpha too large for CCBC
|
||||
}else if(sarg == "--no-fft-gfix"){
|
||||
std::cout << "Not doing the Fourier accelerated gauge fixing tests" << std::endl;
|
||||
do_fft_gfix = false;
|
||||
}else if(sarg == "--alpha"){
|
||||
assert(i<argc-1 && "--alpha option requires an argument");
|
||||
std::istringstream ss(argv[i+1]); ss >> alpha;
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
if(gimpl == "periodic"){
|
||||
std::cout << GridLogMessage << "Using periodic boundary condition" << std::endl;
|
||||
run<PeriodicGimplR>(alpha, do_fft_gfix);
|
||||
}else{
|
||||
std::vector<int> conjdirs = {1,1,0,0}; //test with 2 conjugate dirs and 2 not
|
||||
std::cout << GridLogMessage << "Using complex conjugate boundary conditions in dimensions ";
|
||||
for(int i=0;i<Nd;i++)
|
||||
if(conjdirs[i])
|
||||
std::cout << i << " ";
|
||||
std::cout << std::endl;
|
||||
|
||||
ConjugateGimplR::setDirections(conjdirs);
|
||||
run<ConjugateGimplR>(alpha, do_fft_gfix);
|
||||
}
|
||||
|
||||
Grid_finalize();
|
||||
}
|
||||
|
||||
@ -1,270 +0,0 @@
|
||||
/*************************************************************************************
|
||||
grid` physics library, www.github.com/paboyle/Grid
|
||||
|
||||
Source file: ./tests/Test_cshift.cc
|
||||
|
||||
Copyright (C) 2015
|
||||
|
||||
Author: Azusa Yamaguchi <ayamaguc@staffmail.ed.ac.uk>
|
||||
Author: Peter Boyle <paboyle@ph.ed.ac.uk>
|
||||
|
||||
This program is free software; you can redistribute it and/or modify
|
||||
it under the terms of the GNU General Public License as published by
|
||||
the Free Software Foundation; either version 2 of the License, or
|
||||
(at your option) any later version.
|
||||
|
||||
This program is distributed in the hope that it will be useful,
|
||||
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||
GNU General Public License for more details.
|
||||
|
||||
You should have received a copy of the GNU General Public License along
|
||||
with this program; if not, write to the Free Software Foundation, Inc.,
|
||||
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
|
||||
|
||||
See the full license in the file "LICENSE" in the top level distribution directory
|
||||
*************************************************************************************/
|
||||
/* END LEGAL */
|
||||
#include <Grid/Grid.h>
|
||||
|
||||
using namespace Grid;
|
||||
|
||||
Gamma::Algebra Gmu [] = {
|
||||
Gamma::Algebra::GammaX,
|
||||
Gamma::Algebra::GammaY,
|
||||
Gamma::Algebra::GammaZ,
|
||||
Gamma::Algebra::GammaT,
|
||||
Gamma::Algebra::Gamma5
|
||||
};
|
||||
|
||||
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;
|
||||
|
||||
Coordinate latt_size = GridDefaultLatt();
|
||||
Coordinate simd_layout = GridDefaultSimd(Nd,vComplexD::Nsimd());
|
||||
Coordinate mpi_layout = GridDefaultMpi();
|
||||
|
||||
int vol = 1;
|
||||
for(int d=0;d<latt_size.size();d++){
|
||||
vol = vol * latt_size[d];
|
||||
}
|
||||
GridCartesian GRID(latt_size,simd_layout,mpi_layout);
|
||||
GridRedBlackCartesian RBGRID(&GRID);
|
||||
|
||||
LatticeComplexD coor(&GRID);
|
||||
ComplexD ci(0.0,1.0);
|
||||
|
||||
std::vector<int> seeds({1,2,3,4});
|
||||
GridSerialRNG sRNG; sRNG.SeedFixedIntegers(seeds); // naughty seeding
|
||||
GridParallelRNG pRNG(&GRID);
|
||||
pRNG.SeedFixedIntegers(seeds);
|
||||
|
||||
LatticeGaugeFieldD Umu(&GRID);
|
||||
SU<Nc>::ColdConfiguration(pRNG,Umu); // Unit gauge
|
||||
|
||||
////////////////////////////////////////////////////
|
||||
// Wilson test
|
||||
////////////////////////////////////////////////////
|
||||
{
|
||||
LatticeFermionD src(&GRID); gaussian(pRNG,src);
|
||||
LatticeFermionD src_p(&GRID);
|
||||
LatticeFermionD tmp(&GRID);
|
||||
LatticeFermionD ref(&GRID);
|
||||
LatticeFermionD result(&GRID);
|
||||
|
||||
RealD mass=0.1;
|
||||
WilsonFermionD Dw(Umu,GRID,RBGRID,mass);
|
||||
|
||||
Dw.M(src,ref);
|
||||
std::cout << "Norm src "<<norm2(src)<<std::endl;
|
||||
std::cout << "Norm Dw x src "<<norm2(ref)<<std::endl;
|
||||
{
|
||||
FFT theFFT(&GRID);
|
||||
|
||||
////////////////
|
||||
// operator in Fourier space
|
||||
////////////////
|
||||
tmp =ref;
|
||||
theFFT.FFT_all_dim(result,tmp,FFT::forward);
|
||||
std::cout<<"FFT[ Dw x src ] "<< norm2(result)<<std::endl;
|
||||
|
||||
tmp = src;
|
||||
theFFT.FFT_all_dim(src_p,tmp,FFT::forward);
|
||||
std::cout<<"FFT[ src ] "<< norm2(src_p)<<std::endl;
|
||||
|
||||
/////////////////////////////////////////////////////////////////
|
||||
// work out the predicted FT from Fourier
|
||||
/////////////////////////////////////////////////////////////////
|
||||
auto FGrid = &GRID;
|
||||
LatticeFermionD Kinetic(FGrid); Kinetic = Zero();
|
||||
LatticeComplexD kmu(FGrid);
|
||||
LatticeInteger scoor(FGrid);
|
||||
LatticeComplexD sk (FGrid); sk = Zero();
|
||||
LatticeComplexD sk2(FGrid); sk2= Zero();
|
||||
LatticeComplexD W(FGrid); W= Zero();
|
||||
LatticeComplexD one(FGrid); one =ComplexD(1.0,0.0);
|
||||
ComplexD ci(0.0,1.0);
|
||||
|
||||
for(int mu=0;mu<Nd;mu++) {
|
||||
|
||||
RealD TwoPiL = M_PI * 2.0/ latt_size[mu];
|
||||
|
||||
LatticeCoordinate(kmu,mu);
|
||||
|
||||
kmu = TwoPiL * kmu;
|
||||
|
||||
sk2 = sk2 + 2.0*sin(kmu*0.5)*sin(kmu*0.5);
|
||||
sk = sk + sin(kmu) *sin(kmu);
|
||||
|
||||
// -1/2 Dw -> 1/2 gmu (eip - emip) = i sinp gmu
|
||||
Kinetic = Kinetic + sin(kmu)*ci*(Gamma(Gmu[mu])*src_p);
|
||||
|
||||
}
|
||||
|
||||
W = mass + sk2;
|
||||
Kinetic = Kinetic + W * src_p;
|
||||
|
||||
std::cout<<"Momentum space src "<< norm2(src_p)<<std::endl;
|
||||
std::cout<<"Momentum space Dw x src "<< norm2(Kinetic)<<std::endl;
|
||||
std::cout<<"FT[Coordinate space Dw] "<< norm2(result)<<std::endl;
|
||||
|
||||
result = result - Kinetic;
|
||||
std::cout<<"diff "<< norm2(result)<<std::endl;
|
||||
|
||||
}
|
||||
|
||||
std::cout << " =======================================" <<std::endl;
|
||||
std::cout << " Checking FourierFreePropagator x Dw = 1" <<std::endl;
|
||||
std::cout << " =======================================" <<std::endl;
|
||||
std::cout << "Dw src = " <<norm2(src)<<std::endl;
|
||||
std::cout << "Dw tmp = " <<norm2(tmp)<<std::endl;
|
||||
Dw.M(src,tmp);
|
||||
Dw.FreePropagator(tmp,ref,mass);
|
||||
|
||||
std::cout << "Dw ref = " <<norm2(ref)<<std::endl;
|
||||
|
||||
ref = ref - src;
|
||||
|
||||
std::cout << "Dw ref-src = " <<norm2(ref)<<std::endl;
|
||||
}
|
||||
|
||||
|
||||
////////////////////////////////////////////////////
|
||||
// Wilson prop
|
||||
////////////////////////////////////////////////////
|
||||
{
|
||||
std::cout<<"****************************************"<<std::endl;
|
||||
std::cout << "Wilson Mom space 4d propagator \n";
|
||||
std::cout<<"****************************************"<<std::endl;
|
||||
|
||||
LatticeFermionD src(&GRID); gaussian(pRNG,src);
|
||||
LatticeFermionD tmp(&GRID);
|
||||
LatticeFermionD ref(&GRID);
|
||||
LatticeFermionD diff(&GRID);
|
||||
|
||||
src=Zero();
|
||||
Coordinate point(4,0); // 0,0,0,0
|
||||
SpinColourVectorD ferm;
|
||||
ferm=Zero();
|
||||
ferm()(0)(0) = ComplexD(1.0);
|
||||
pokeSite(ferm,src,point);
|
||||
|
||||
RealD mass=0.1;
|
||||
WilsonFermionD Dw(Umu,GRID,RBGRID,mass);
|
||||
|
||||
// Momentum space prop
|
||||
std::cout << " Solving by FFT and Feynman rules" <<std::endl;
|
||||
Dw.FreePropagator(src,ref,mass) ;
|
||||
|
||||
Gamma G5(Gamma::Algebra::Gamma5);
|
||||
|
||||
LatticeFermionD result(&GRID);
|
||||
const int sdir=0;
|
||||
|
||||
////////////////////////////////////////////////////////////////////////
|
||||
// Conjugate gradient on normal equations system
|
||||
////////////////////////////////////////////////////////////////////////
|
||||
std::cout << " Solving by Conjugate Gradient (CGNE)" <<std::endl;
|
||||
Dw.Mdag(src,tmp);
|
||||
src=tmp;
|
||||
MdagMLinearOperator<WilsonFermionD,LatticeFermionD> HermOp(Dw);
|
||||
ConjugateGradient<LatticeFermionD> CG(1.0e-10,10000);
|
||||
CG(HermOp,src,result);
|
||||
|
||||
////////////////////////////////////////////////////////////////////////
|
||||
std::cout << " Taking difference" <<std::endl;
|
||||
std::cout << "Dw result "<<norm2(result)<<std::endl;
|
||||
std::cout << "Dw ref "<<norm2(ref)<<std::endl;
|
||||
|
||||
diff = ref - result;
|
||||
std::cout << "result - ref "<<norm2(diff)<<std::endl;
|
||||
|
||||
DumpSliceNorm("Slice Norm Solution ",result,Nd-1);
|
||||
}
|
||||
|
||||
////////////////////////////////////////////////////
|
||||
//Gauge invariance test
|
||||
////////////////////////////////////////////////////
|
||||
{
|
||||
std::cout<<"****************************************"<<std::endl;
|
||||
std::cout << "Gauge invariance test \n";
|
||||
std::cout<<"****************************************"<<std::endl;
|
||||
LatticeGaugeField U_GT(&GRID); // Gauge transformed field
|
||||
LatticeColourMatrix g(&GRID); // local Gauge xform matrix
|
||||
U_GT = Umu;
|
||||
// Make a random xform to teh gauge field
|
||||
SU<Nc>::RandomGaugeTransform(pRNG,U_GT,g); // Unit gauge
|
||||
|
||||
LatticeFermionD src(&GRID);
|
||||
LatticeFermionD tmp(&GRID);
|
||||
LatticeFermionD ref(&GRID);
|
||||
LatticeFermionD diff(&GRID);
|
||||
|
||||
// could loop over colors
|
||||
src=Zero();
|
||||
Coordinate point(4,0); // 0,0,0,0
|
||||
SpinColourVectorD ferm;
|
||||
ferm=Zero();
|
||||
ferm()(0)(0) = ComplexD(1.0);
|
||||
pokeSite(ferm,src,point);
|
||||
|
||||
RealD mass=0.1;
|
||||
WilsonFermionD Dw(U_GT,GRID,RBGRID,mass);
|
||||
|
||||
// Momentum space prop
|
||||
std::cout << " Solving by FFT and Feynman rules" <<std::endl;
|
||||
Dw.FreePropagator(src,ref,mass) ;
|
||||
|
||||
Gamma G5(Gamma::Algebra::Gamma5);
|
||||
|
||||
LatticeFermionD result(&GRID);
|
||||
const int sdir=0;
|
||||
|
||||
////////////////////////////////////////////////////////////////////////
|
||||
// Conjugate gradient on normal equations system
|
||||
////////////////////////////////////////////////////////////////////////
|
||||
std::cout << " Solving by Conjugate Gradient (CGNE)" <<std::endl;
|
||||
Dw.Mdag(src,tmp);
|
||||
src=tmp;
|
||||
MdagMLinearOperator<WilsonFermionD,LatticeFermionD> HermOp(Dw);
|
||||
ConjugateGradient<LatticeFermionD> CG(1.0e-10,10000);
|
||||
CG(HermOp,src,result);
|
||||
|
||||
////////////////////////////////////////////////////////////////////////
|
||||
std::cout << " Taking difference" <<std::endl;
|
||||
std::cout << "Dw result "<<norm2(result)<<std::endl;
|
||||
std::cout << "Dw ref "<<norm2(ref)<<std::endl;
|
||||
|
||||
diff = ref - result;
|
||||
std::cout << "result - ref "<<norm2(diff)<<std::endl;
|
||||
|
||||
DumpSliceNorm("Slice Norm Solution ",result,Nd-1);
|
||||
}
|
||||
|
||||
|
||||
Grid_finalize();
|
||||
}
|
||||
@ -228,6 +228,59 @@ void checkGammaL(const Gamma::Algebra a, GridSerialRNG &rng)
|
||||
std::cout << std::endl;
|
||||
}
|
||||
|
||||
void checkChargeConjMatrix(){
|
||||
//Check the properties of the charge conjugation matrix
|
||||
//In the Grid basis C = -\gamma^2 \gamma^4
|
||||
SpinMatrix C = testAlgebra[Gamma::Algebra::MinusGammaY] * testAlgebra[Gamma::Algebra::GammaT];
|
||||
SpinMatrix mC = -C;
|
||||
SpinMatrix one = testAlgebra[Gamma::Algebra::Identity];
|
||||
|
||||
std::cout << "Testing properties of charge conjugation matrix C = -\\gamma^2 \\gamma^4 (in Grid's basis)" << std::endl;
|
||||
|
||||
//C^T = -C
|
||||
SpinMatrix Ct = transpose(C);
|
||||
std::cout << GridLogMessage << "C^T=-C ";
|
||||
test(Ct, mC);
|
||||
std::cout << std::endl;
|
||||
|
||||
//C^\dagger = -C
|
||||
SpinMatrix Cdag = adj(C);
|
||||
std::cout << GridLogMessage << "C^dag=-C ";
|
||||
test(Cdag, mC);
|
||||
std::cout << std::endl;
|
||||
|
||||
//C^* = C
|
||||
SpinMatrix Cstar = conjugate(C);
|
||||
std::cout << GridLogMessage << "C^*=C ";
|
||||
test(Cstar, C);
|
||||
std::cout << std::endl;
|
||||
|
||||
//C^{-1} = -C
|
||||
SpinMatrix CinvC = mC * C;
|
||||
std::cout << GridLogMessage << "C^{-1}=-C ";
|
||||
test(CinvC, one);
|
||||
std::cout << std::endl;
|
||||
|
||||
// C^{-1} \gamma^\mu C = -[\gamma^\mu]^T
|
||||
Gamma::Algebra gmu_a[4] = { Gamma::Algebra::GammaX, Gamma::Algebra::GammaY, Gamma::Algebra::GammaZ, Gamma::Algebra::GammaT };
|
||||
for(int mu=0;mu<4;mu++){
|
||||
SpinMatrix gmu = testAlgebra[gmu_a[mu]];
|
||||
SpinMatrix Cinv_gmu_C = mC * gmu * C;
|
||||
SpinMatrix mgmu_T = -transpose(gmu);
|
||||
std::cout << GridLogMessage << "C^{-1} \\gamma^" << mu << " C = -[\\gamma^" << mu << "]^T ";
|
||||
test(Cinv_gmu_C, mgmu_T);
|
||||
std::cout << std::endl;
|
||||
}
|
||||
|
||||
//[C, \gamma^5] = 0
|
||||
SpinMatrix Cg5 = C * testAlgebra[Gamma::Algebra::Gamma5];
|
||||
SpinMatrix g5C = testAlgebra[Gamma::Algebra::Gamma5] * C;
|
||||
std::cout << GridLogMessage << "C \\gamma^5 = \\gamma^5 C";
|
||||
test(Cg5, g5C);
|
||||
std::cout << std::endl;
|
||||
}
|
||||
|
||||
|
||||
int main(int argc, char *argv[])
|
||||
{
|
||||
Grid_init(&argc,&argv);
|
||||
@ -270,6 +323,13 @@ int main(int argc, char *argv[])
|
||||
{
|
||||
checkGammaL(i, sRNG);
|
||||
}
|
||||
|
||||
std::cout << GridLogMessage << "======== Charge conjugation matrix check" << std::endl;
|
||||
checkChargeConjMatrix();
|
||||
std::cout << GridLogMessage << std::endl;
|
||||
|
||||
|
||||
|
||||
|
||||
Grid_finalize();
|
||||
|
||||
|
||||
@ -55,13 +55,17 @@ static_assert(same_vComplex == 1, "Dirac Operators must have same underlying SIM
|
||||
int main (int argc, char ** argv)
|
||||
{
|
||||
int nu = 0;
|
||||
|
||||
int tbc_aprd = 0; //use antiperiodic BCs in the time direction?
|
||||
|
||||
Grid_init(&argc,&argv);
|
||||
|
||||
for(int i=1;i<argc;i++){
|
||||
if(std::string(argv[i]) == "--Gparity-dir"){
|
||||
std::stringstream ss; ss << argv[i+1]; ss >> nu;
|
||||
std::cout << GridLogMessage << "Set Gparity direction to " << nu << std::endl;
|
||||
}else if(std::string(argv[i]) == "--Tbc-APRD"){
|
||||
tbc_aprd = 1;
|
||||
std::cout << GridLogMessage << "Using antiperiodic BCs in the time direction" << std::endl;
|
||||
}
|
||||
}
|
||||
|
||||
@ -155,13 +159,18 @@ int main (int argc, char ** argv)
|
||||
|
||||
//Coordinate grid for reference
|
||||
LatticeInteger xcoor_1f5(FGrid_1f);
|
||||
LatticeCoordinate(xcoor_1f5,1+nu);
|
||||
LatticeCoordinate(xcoor_1f5,1+nu); //note '1+nu'! This is because for 5D fields the s-direction is direction 0
|
||||
Replicate(src,src_1f);
|
||||
src_1f = where( xcoor_1f5 >= Integer(L), 2.0*src_1f,src_1f );
|
||||
|
||||
RealD mass=0.0;
|
||||
RealD M5=1.8;
|
||||
StandardDiracOp Ddwf(Umu_1f,*FGrid_1f,*FrbGrid_1f,*UGrid_1f,*UrbGrid_1f,mass,M5 DOP_PARAMS);
|
||||
|
||||
//Standard Dirac op
|
||||
AcceleratorVector<Complex,4> bc_std(Nd, 1.0);
|
||||
if(tbc_aprd) bc_std[Nd-1] = -1.; //antiperiodic time BC
|
||||
StandardDiracOp::ImplParams std_params(bc_std);
|
||||
StandardDiracOp Ddwf(Umu_1f,*FGrid_1f,*FrbGrid_1f,*UGrid_1f,*UrbGrid_1f,mass,M5 DOP_PARAMS, std_params);
|
||||
|
||||
StandardFermionField src_o_1f(FrbGrid_1f);
|
||||
StandardFermionField result_o_1f(FrbGrid_1f);
|
||||
@ -172,9 +181,11 @@ int main (int argc, char ** argv)
|
||||
ConjugateGradient<StandardFermionField> CG(1.0e-8,10000);
|
||||
CG(HermOpEO,src_o_1f,result_o_1f);
|
||||
|
||||
// const int nu = 3;
|
||||
//Gparity Dirac op
|
||||
std::vector<int> twists(Nd,0);
|
||||
twists[nu] = 1;
|
||||
if(tbc_aprd) twists[Nd-1] = 1;
|
||||
|
||||
GparityDiracOp::ImplParams params;
|
||||
params.twists = twists;
|
||||
GparityDiracOp GPDdwf(Umu_2f,*FGrid_2f,*FrbGrid_2f,*UGrid_2f,*UrbGrid_2f,mass,M5 DOP_PARAMS,params);
|
||||
@ -271,8 +282,11 @@ int main (int argc, char ** argv)
|
||||
std::cout << "2f cb "<<result_o_2f.Checkerboard()<<std::endl;
|
||||
std::cout << "1f cb "<<result_o_1f.Checkerboard()<<std::endl;
|
||||
|
||||
std::cout << " result norms " <<norm2(result_o_2f)<<" " <<norm2(result_o_1f)<<std::endl;
|
||||
//Compare norms
|
||||
std::cout << " result norms 2f: " <<norm2(result_o_2f)<<" 1f: " <<norm2(result_o_1f)<<std::endl;
|
||||
|
||||
|
||||
//Take the 2f solution and convert into the corresponding 1f solution (odd cb only)
|
||||
StandardFermionField res0o (FrbGrid_2f);
|
||||
StandardFermionField res1o (FrbGrid_2f);
|
||||
StandardFermionField res0 (FGrid_2f);
|
||||
@ -281,14 +295,15 @@ int main (int argc, char ** argv)
|
||||
res0=Zero();
|
||||
res1=Zero();
|
||||
|
||||
res0o = PeekIndex<0>(result_o_2f,0);
|
||||
res1o = PeekIndex<0>(result_o_2f,1);
|
||||
res0o = PeekIndex<0>(result_o_2f,0); //flavor 0, odd cb
|
||||
res1o = PeekIndex<0>(result_o_2f,1); //flavor 1, odd cb
|
||||
|
||||
std::cout << "res cb "<<res0o.Checkerboard()<<std::endl;
|
||||
std::cout << "res cb "<<res1o.Checkerboard()<<std::endl;
|
||||
|
||||
setCheckerboard(res0,res0o);
|
||||
setCheckerboard(res1,res1o);
|
||||
//poke odd onto non-cb field
|
||||
setCheckerboard(res0,res0o);
|
||||
setCheckerboard(res1,res1o);
|
||||
|
||||
StandardFermionField replica (FGrid_1f);
|
||||
StandardFermionField replica0(FGrid_1f);
|
||||
@ -296,12 +311,13 @@ int main (int argc, char ** argv)
|
||||
Replicate(res0,replica0);
|
||||
Replicate(res1,replica1);
|
||||
|
||||
//2nd half of doubled lattice has f=1
|
||||
replica = where( xcoor_1f5 >= Integer(L), replica1,replica0 );
|
||||
|
||||
replica0 = Zero();
|
||||
setCheckerboard(replica0,result_o_1f);
|
||||
|
||||
std::cout << "Norm2 solutions is " <<norm2(replica)<<" "<< norm2(replica0)<<std::endl;
|
||||
std::cout << "Norm2 solutions 1f reconstructed from 2f: " <<norm2(replica)<<" Actual 1f: "<< norm2(replica0)<<std::endl;
|
||||
|
||||
replica = replica - replica0;
|
||||
|
||||
|
||||
177
tests/core/Test_gparity_flavour.cc
Normal file
177
tests/core/Test_gparity_flavour.cc
Normal file
@ -0,0 +1,177 @@
|
||||
/*************************************************************************************
|
||||
|
||||
Grid physics library, www.github.com/paboyle/Grid
|
||||
|
||||
Source file: ./tests/Test_gparity_flavour.cc
|
||||
|
||||
Copyright (C) 2015-2017
|
||||
|
||||
Author: Christopher Kelly <ckelly@bnl.gov>
|
||||
Author: Peter Boyle <paboyle@ph.ed.ac.uk>
|
||||
|
||||
This program is free software; you can redistribute it and/or modify
|
||||
it under the terms of the GNU General Public License as published by
|
||||
the Free Software Foundation; either version 2 of the License, or
|
||||
(at your option) any later version.
|
||||
|
||||
This program is distributed in the hope that it will be useful,
|
||||
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||
GNU General Public License for more details.
|
||||
|
||||
You should have received a copy of the GNU General Public License along
|
||||
with this program; if not, write to the Free Software Foundation, Inc.,
|
||||
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
|
||||
|
||||
See the full license in the file "LICENSE" in the top level distribution directory
|
||||
*************************************************************************************/
|
||||
/* END LEGAL */
|
||||
#include <Grid/Grid.h>
|
||||
|
||||
using namespace Grid;
|
||||
|
||||
static constexpr double tolerance = 1.0e-6;
|
||||
static std::array<GparityFlavourMatrix, GparityFlavour::nSigma> testAlgebra;
|
||||
|
||||
void print(const GparityFlavourMatrix &g)
|
||||
{
|
||||
for(int i = 0; i < Ngp; i++)
|
||||
{
|
||||
std::cout << GridLogMessage << "(";
|
||||
for(int j=0;j<Ngp;j++){
|
||||
if ( abs( g(i,j)()() ) == 0 ) {
|
||||
std::cout<< " 0";
|
||||
} else if ( abs(g(i,j)()() - Complex(0,1)) == 0){
|
||||
std::cout<< " i";
|
||||
} else if ( abs(g(i,j)()() + Complex(0,1)) == 0){
|
||||
std::cout<< "-i";
|
||||
} else if ( abs(g(i,j)()() - Complex(1,0)) == 0){
|
||||
std::cout<< " 1";
|
||||
} else if ( abs(g(i,j)()() + Complex(1,0)) == 0){
|
||||
std::cout<< "-1";
|
||||
}
|
||||
std::cout<<((j == Ngp-1) ? ")" : "," );
|
||||
}
|
||||
std::cout << std::endl;
|
||||
}
|
||||
std::cout << GridLogMessage << std::endl;
|
||||
}
|
||||
|
||||
void createTestAlgebra(void)
|
||||
{
|
||||
std::array<GparityFlavourMatrix, 3> testg;
|
||||
const Complex I(0., 1.), mI(0., -1.);
|
||||
|
||||
// 0 1
|
||||
// 1 0
|
||||
testg[0] = Zero();
|
||||
testg[0](0, 1)()() = 1.;
|
||||
testg[0](1, 0)()() = 1.;
|
||||
std::cout << GridLogMessage << "test SigmaX= " << std::endl;
|
||||
print(testg[0]);
|
||||
|
||||
// 0 -i
|
||||
// i 0
|
||||
testg[1] = Zero();
|
||||
testg[1](0, 1)()() = mI;
|
||||
testg[1](1, 0)()() = I;
|
||||
std::cout << GridLogMessage << "test SigmaY= " << std::endl;
|
||||
print(testg[1]);
|
||||
|
||||
// 1 0
|
||||
// 0 -1
|
||||
testg[2] = Zero();
|
||||
testg[2](0, 0)()() = 1.0;
|
||||
testg[2](1, 1)()() = -1.0;
|
||||
std::cout << GridLogMessage << "test SigmaZ= " << std::endl;
|
||||
print(testg[2]);
|
||||
|
||||
|
||||
#define DEFINE_TEST_G(g, exp)\
|
||||
testAlgebra[GparityFlavour::Algebra::g] = exp; \
|
||||
testAlgebra[GparityFlavour::Algebra::Minus##g] = -exp;
|
||||
|
||||
DEFINE_TEST_G(SigmaX , testg[0]);
|
||||
DEFINE_TEST_G(SigmaY , testg[1]);
|
||||
DEFINE_TEST_G(SigmaZ , testg[2]);
|
||||
DEFINE_TEST_G(Identity , 1.);
|
||||
|
||||
GparityFlavourMatrix pplus;
|
||||
pplus = 1.0;
|
||||
pplus = pplus + testg[1];
|
||||
pplus = pplus * 0.5;
|
||||
|
||||
DEFINE_TEST_G(ProjPlus , pplus);
|
||||
|
||||
GparityFlavourMatrix pminus;
|
||||
pminus = 1.0;
|
||||
pminus = pminus - testg[1];
|
||||
pminus = pminus * 0.5;
|
||||
|
||||
DEFINE_TEST_G(ProjMinus , pminus);
|
||||
|
||||
#undef DEFINE_TEST_G
|
||||
}
|
||||
|
||||
template <typename Expr>
|
||||
void test(const Expr &a, const Expr &b)
|
||||
{
|
||||
if (norm2(a - b) < tolerance)
|
||||
{
|
||||
std::cout << "[OK] ";
|
||||
}
|
||||
else
|
||||
{
|
||||
std::cout << "[fail]" << std::endl;
|
||||
std::cout << GridLogError << "a= " << a << std::endl;
|
||||
std::cout << GridLogError << "is different (tolerance= " << tolerance << ") from " << std::endl;
|
||||
std::cout << GridLogError << "b= " << b << std::endl;
|
||||
exit(EXIT_FAILURE);
|
||||
}
|
||||
}
|
||||
|
||||
void checkSigma(const GparityFlavour::Algebra a, GridSerialRNG &rng)
|
||||
{
|
||||
GparityFlavourVector v;
|
||||
GparityFlavourMatrix m, &testg = testAlgebra[a];
|
||||
GparityFlavour g(a);
|
||||
|
||||
random(rng, v);
|
||||
random(rng, m);
|
||||
|
||||
std::cout << GridLogMessage << "Checking " << GparityFlavour::name[a] << ": ";
|
||||
std::cout << "vecmul ";
|
||||
test(g*v, testg*v);
|
||||
std::cout << "matlmul ";
|
||||
test(g*m, testg*m);
|
||||
std::cout << "matrmul ";
|
||||
test(m*g, m*testg);
|
||||
std::cout << std::endl;
|
||||
}
|
||||
|
||||
int main(int argc, char *argv[])
|
||||
{
|
||||
Grid_init(&argc,&argv);
|
||||
|
||||
Coordinate latt_size = GridDefaultLatt();
|
||||
Coordinate simd_layout = GridDefaultSimd(4,vComplex::Nsimd());
|
||||
Coordinate mpi_layout = GridDefaultMpi();
|
||||
|
||||
GridCartesian Grid(latt_size,simd_layout,mpi_layout);
|
||||
GridSerialRNG sRNG;
|
||||
|
||||
sRNG.SeedFixedIntegers(std::vector<int>({45,12,81,9}));
|
||||
|
||||
std::cout << GridLogMessage << "======== Test algebra" << std::endl;
|
||||
createTestAlgebra();
|
||||
std::cout << GridLogMessage << "======== Multiplication operators check" << std::endl;
|
||||
for (int i = 0; i < GparityFlavour::nSigma; ++i)
|
||||
{
|
||||
checkSigma(i, sRNG);
|
||||
}
|
||||
std::cout << GridLogMessage << std::endl;
|
||||
|
||||
Grid_finalize();
|
||||
|
||||
return EXIT_SUCCESS;
|
||||
}
|
||||
@ -52,7 +52,7 @@ int main (int argc, char ** argv)
|
||||
// pRNG.SeedFixedIntegers(seeds);
|
||||
pRNG.SeedFixedIntegers(std::vector<int>({45,12,81,9}));
|
||||
|
||||
typedef typename GparityWilsonFermionR::FermionField FermionField;
|
||||
typedef typename GparityWilsonFermionD::FermionField FermionField;
|
||||
|
||||
FermionField src (&Grid); random(pRNG,src);
|
||||
FermionField phi (&Grid); random(pRNG,phi);
|
||||
@ -80,10 +80,10 @@ int main (int argc, char ** argv)
|
||||
|
||||
RealD mass=0.1;
|
||||
|
||||
GparityWilsonFermionR::ImplParams params;
|
||||
GparityWilsonFermionD::ImplParams params;
|
||||
std::vector<int> twists(Nd,0); twists[1] = 1;
|
||||
params.twists = twists;
|
||||
GparityWilsonFermionR Dw(Umu,Grid,RBGrid,mass,params);
|
||||
GparityWilsonFermionD Dw(Umu,Grid,RBGrid,mass,params);
|
||||
|
||||
FermionField src_e (&RBGrid);
|
||||
FermionField src_o (&RBGrid);
|
||||
@ -199,7 +199,7 @@ int main (int argc, char ** argv)
|
||||
pickCheckerboard(Even,phi_e,phi);
|
||||
pickCheckerboard(Odd ,phi_o,phi);
|
||||
|
||||
SchurDiagMooeeOperator<GparityWilsonFermionR,FermionField> HermOpEO(Dw);
|
||||
SchurDiagMooeeOperator<GparityWilsonFermionD,FermionField> HermOpEO(Dw);
|
||||
HermOpEO.MpcDagMpc(chi_e,dchi_e);
|
||||
HermOpEO.MpcDagMpc(chi_o,dchi_o);
|
||||
|
||||
|
||||
@ -92,7 +92,7 @@ int main (int argc, char ** argv)
|
||||
RealD shift = 0.1234;
|
||||
RealD M5 = 1.8;
|
||||
int pm = 1;
|
||||
MobiusEOFAFermionR Ddwf(Umu, *FGrid, *FrbGrid, *UGrid, *UrbGrid, mq1, mq2, mq3, shift, pm, M5, b, c);
|
||||
MobiusEOFAFermionD Ddwf(Umu, *FGrid, *FrbGrid, *UGrid, *UrbGrid, mq1, mq2, mq3, shift, pm, M5, b, c);
|
||||
|
||||
LatticeFermion src_e (FrbGrid);
|
||||
LatticeFermion src_o (FrbGrid);
|
||||
@ -218,7 +218,7 @@ int main (int argc, char ** argv)
|
||||
pickCheckerboard(Even, phi_e, phi);
|
||||
pickCheckerboard(Odd , phi_o, phi);
|
||||
|
||||
SchurDiagMooeeOperator<MobiusEOFAFermionR,LatticeFermion> HermOpEO(Ddwf);
|
||||
SchurDiagMooeeOperator<MobiusEOFAFermionD,LatticeFermion> HermOpEO(Ddwf);
|
||||
HermOpEO.MpcDagMpc(chi_e, dchi_e);
|
||||
HermOpEO.MpcDagMpc(chi_o, dchi_o);
|
||||
|
||||
|
||||
@ -108,8 +108,8 @@ int main (int argc, char ** argv)
|
||||
omegas.push_back( std::complex<double>(0.0686324988446592,-0.0550658530827402) );
|
||||
#endif
|
||||
|
||||
MobiusFermionR Ddwf(Umu, *FGrid, *FrbGrid, *UGrid, *UrbGrid, mass, M5, 0.5,0.5);
|
||||
// DomainWallFermionR Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5);
|
||||
MobiusFermionD Ddwf(Umu, *FGrid, *FrbGrid, *UGrid, *UrbGrid, mass, M5, 0.5,0.5);
|
||||
// DomainWallFermionD Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5);
|
||||
|
||||
LatticeFermion src_e (FrbGrid);
|
||||
LatticeFermion src_o (FrbGrid);
|
||||
@ -264,7 +264,7 @@ int main (int argc, char ** argv)
|
||||
pickCheckerboard(Odd ,phi_o,phi);
|
||||
|
||||
|
||||
SchurDiagMooeeOperator<MobiusFermionR,LatticeFermion> HermOpEO(Ddwf);
|
||||
SchurDiagMooeeOperator<MobiusFermionD,LatticeFermion> HermOpEO(Ddwf);
|
||||
HermOpEO.MpcDagMpc(chi_e,dchi_e);
|
||||
HermOpEO.MpcDagMpc(chi_o,dchi_o);
|
||||
|
||||
|
||||
124
tests/core/Test_prec_change.cc
Normal file
124
tests/core/Test_prec_change.cc
Normal file
@ -0,0 +1,124 @@
|
||||
/*************************************************************************************
|
||||
|
||||
Grid physics library, www.github.com/paboyle/Grid
|
||||
|
||||
Source file: ./tests/core/Test_prec_change.cc
|
||||
|
||||
Copyright (C) 2015
|
||||
|
||||
Author: Christopher Kelly <ckelly@bnl.gov>
|
||||
Author: Peter Boyle <paboyle@ph.ed.ac.uk>
|
||||
|
||||
This program is free software; you can redistribute it and/or modify
|
||||
it under the terms of the GNU General Public License as published by
|
||||
the Free Software Foundation; either version 2 of the License, or
|
||||
(at your option) any later version.
|
||||
|
||||
This program is distributed in the hope that it will be useful,
|
||||
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||
GNU General Public License for more details.
|
||||
|
||||
You should have received a copy of the GNU General Public License along
|
||||
with this program; if not, write to the Free Software Foundation, Inc.,
|
||||
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
|
||||
|
||||
See the full license in the file "LICENSE" in the top level distribution directory
|
||||
*************************************************************************************/
|
||||
/* END LEGAL */
|
||||
#include <Grid/Grid.h>
|
||||
|
||||
using namespace std;
|
||||
using namespace Grid;
|
||||
|
||||
int main (int argc, char ** argv)
|
||||
{
|
||||
Grid_init(&argc,&argv);
|
||||
|
||||
int Ls = 12;
|
||||
Coordinate latt4 = GridDefaultLatt();
|
||||
|
||||
GridCartesian * UGridD = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplexD::Nsimd()),GridDefaultMpi());
|
||||
GridRedBlackCartesian * UrbGridD = SpaceTimeGrid::makeFourDimRedBlackGrid(UGridD);
|
||||
GridCartesian * FGridD = SpaceTimeGrid::makeFiveDimGrid(Ls,UGridD);
|
||||
GridRedBlackCartesian * FrbGridD = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGridD);
|
||||
|
||||
GridCartesian * UGridF = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplexF::Nsimd()),GridDefaultMpi());
|
||||
GridRedBlackCartesian * UrbGridF = SpaceTimeGrid::makeFourDimRedBlackGrid(UGridF);
|
||||
GridCartesian * FGridF = SpaceTimeGrid::makeFiveDimGrid(Ls,UGridF);
|
||||
GridRedBlackCartesian * FrbGridF = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGridF);
|
||||
|
||||
|
||||
std::vector<int> seeds4({1,2,3,4});
|
||||
std::vector<int> seeds5({5,6,7,8});
|
||||
|
||||
std::cout << GridLogMessage << "Initialising 5d RNG" << std::endl;
|
||||
GridParallelRNG RNG5(FGridD); RNG5.SeedFixedIntegers(seeds5);
|
||||
GridParallelRNG RNG5F(FGridF); RNG5F.SeedFixedIntegers(seeds5);
|
||||
std::cout << GridLogMessage << "Initialised RNGs" << std::endl;
|
||||
|
||||
LatticeFermionD field_d(FGridD), tmp_d(FGridD);
|
||||
random(RNG5,field_d);
|
||||
RealD norm2_field_d = norm2(field_d);
|
||||
|
||||
LatticeFermionD2 field_d2(FGridF), tmp_d2(FGridF);
|
||||
random(RNG5F,field_d2);
|
||||
RealD norm2_field_d2 = norm2(field_d2);
|
||||
|
||||
LatticeFermionF field_f(FGridF);
|
||||
|
||||
//Test original implementation
|
||||
{
|
||||
std::cout << GridLogMessage << "Testing original implementation" << std::endl;
|
||||
field_f = Zero();
|
||||
precisionChangeOrig(field_f,field_d);
|
||||
RealD Ndiff = (norm2_field_d - norm2(field_f))/norm2_field_d;
|
||||
std::cout << GridLogMessage << (fabs(Ndiff) > 1e-05 ? "!!FAIL" : "Pass") << ": relative norm2 of single and double prec fields differs by " << Ndiff << std::endl;
|
||||
tmp_d = Zero();
|
||||
precisionChangeOrig(tmp_d, field_f);
|
||||
Ndiff = norm2( LatticeFermionD(tmp_d-field_d) ) / norm2_field_d;
|
||||
std::cout << GridLogMessage << (fabs(Ndiff) > 1e-05 ? "!!FAIL" : "Pass") << ": relative norm2 of back-converted and original double prec fields differs by " << Ndiff << std::endl;
|
||||
}
|
||||
//Test new implementation with pregenerated workspace
|
||||
{
|
||||
std::cout << GridLogMessage << "Testing new implementation with pregenerated workspace" << std::endl;
|
||||
precisionChangeWorkspace wk_sp_to_dp(field_d.Grid(),field_f.Grid());
|
||||
precisionChangeWorkspace wk_dp_to_sp(field_f.Grid(),field_d.Grid());
|
||||
|
||||
field_f = Zero();
|
||||
precisionChange(field_f,field_d,wk_dp_to_sp);
|
||||
RealD Ndiff = (norm2_field_d - norm2(field_f))/norm2_field_d;
|
||||
std::cout << GridLogMessage << (fabs(Ndiff) > 1e-05 ? "!!FAIL" : "Pass") << ": relative norm2 of single and double prec fields differs by " << Ndiff << std::endl;
|
||||
tmp_d = Zero();
|
||||
precisionChange(tmp_d, field_f,wk_sp_to_dp);
|
||||
Ndiff = norm2( LatticeFermionD(tmp_d-field_d) ) / norm2_field_d;
|
||||
std::cout << GridLogMessage << (fabs(Ndiff) > 1e-05 ? "!!FAIL" : "Pass") << ": relative norm2 of back-converted and original double prec fields differs by " << Ndiff << std::endl;
|
||||
}
|
||||
//Test new implementation without pregenerated workspace
|
||||
{
|
||||
std::cout << GridLogMessage << "Testing new implementation without pregenerated workspace" << std::endl;
|
||||
field_f = Zero();
|
||||
precisionChange(field_f,field_d);
|
||||
RealD Ndiff = (norm2_field_d - norm2(field_f))/norm2_field_d;
|
||||
std::cout << GridLogMessage << (fabs(Ndiff) > 1e-05 ? "!!FAIL" : "Pass") << ": relative norm2 of single and double prec fields differs by " << Ndiff << std::endl;
|
||||
tmp_d = Zero();
|
||||
precisionChange(tmp_d, field_f);
|
||||
Ndiff = norm2( LatticeFermionD(tmp_d-field_d) ) / norm2_field_d;
|
||||
std::cout << GridLogMessage << (fabs(Ndiff) > 1e-05 ? "!!FAIL" : "Pass") << ": relative norm2 of back-converted and original double prec fields differs by " << Ndiff << std::endl;
|
||||
}
|
||||
//Test fast implementation
|
||||
{
|
||||
std::cout << GridLogMessage << "Testing fast (double2) implementation" << std::endl;
|
||||
field_f = Zero();
|
||||
precisionChangeFast(field_f,field_d2);
|
||||
RealD Ndiff = (norm2_field_d2 - norm2(field_f))/norm2_field_d2;
|
||||
std::cout << GridLogMessage << (fabs(Ndiff) > 1e-05 ? "!!FAIL" : "Pass") << ": relative norm2 of single and double prec fields differs by " << Ndiff << std::endl;
|
||||
tmp_d2 = Zero();
|
||||
precisionChangeFast(tmp_d2, field_f);
|
||||
Ndiff = norm2( LatticeFermionD2(tmp_d2-field_d2) ) / norm2_field_d2;
|
||||
std::cout << GridLogMessage << (fabs(Ndiff) > 1e-05 ? "!!FAIL" : "Pass") << ": relative norm2 of back-converted and original double prec fields differs by " << Ndiff << std::endl;
|
||||
}
|
||||
std::cout << "Done" << std::endl;
|
||||
|
||||
Grid_finalize();
|
||||
}
|
||||
@ -53,9 +53,9 @@ int main (int argc, char ** argv)
|
||||
pRNG.SeedFixedIntegers(seeds);
|
||||
// pRNG.SeedFixedIntegers(std::vector<int>({45,12,81,9});
|
||||
|
||||
typedef typename ImprovedStaggeredFermionR::FermionField FermionField;
|
||||
typedef typename ImprovedStaggeredFermionR::ComplexField ComplexField;
|
||||
typename ImprovedStaggeredFermionR::ImplParams params;
|
||||
typedef typename ImprovedStaggeredFermionD::FermionField FermionField;
|
||||
typedef typename ImprovedStaggeredFermionD::ComplexField ComplexField;
|
||||
typename ImprovedStaggeredFermionD::ImplParams params;
|
||||
|
||||
FermionField src (&Grid); random(pRNG,src);
|
||||
FermionField result(&Grid); result=Zero();
|
||||
@ -130,7 +130,7 @@ int main (int argc, char ** argv)
|
||||
// ref = ref + mass * src;
|
||||
}
|
||||
|
||||
ImprovedStaggeredFermionR Ds(Umu,Umu,Grid,RBGrid,mass,c1,c2,u0,params);
|
||||
ImprovedStaggeredFermionD Ds(Umu,Umu,Grid,RBGrid,mass,c1,c2,u0,params);
|
||||
|
||||
|
||||
std::cout<<GridLogMessage<<"=========================================================="<<std::endl;
|
||||
@ -269,7 +269,7 @@ int main (int argc, char ** argv)
|
||||
pickCheckerboard(Even,phi_e,phi);
|
||||
pickCheckerboard(Odd ,phi_o,phi);
|
||||
|
||||
SchurDiagMooeeOperator<ImprovedStaggeredFermionR,FermionField> HermOpEO(Ds);
|
||||
SchurDiagMooeeOperator<ImprovedStaggeredFermionD,FermionField> HermOpEO(Ds);
|
||||
HermOpEO.MpcDagMpc(chi_e,dchi_e);
|
||||
HermOpEO.MpcDagMpc(chi_o,dchi_o);
|
||||
|
||||
|
||||
@ -60,9 +60,9 @@ int main (int argc, char ** argv)
|
||||
pRNG4.SeedFixedIntegers(seeds);
|
||||
pRNG5.SeedFixedIntegers(seeds);
|
||||
|
||||
typedef typename ImprovedStaggeredFermion5DR::FermionField FermionField;
|
||||
typedef typename ImprovedStaggeredFermion5DR::ComplexField ComplexField;
|
||||
typename ImprovedStaggeredFermion5DR::ImplParams params;
|
||||
typedef typename ImprovedStaggeredFermion5DD::FermionField FermionField;
|
||||
typedef typename ImprovedStaggeredFermion5DD::ComplexField ComplexField;
|
||||
typename ImprovedStaggeredFermion5DD::ImplParams params;
|
||||
|
||||
FermionField src (FGrid);
|
||||
|
||||
@ -148,7 +148,7 @@ int main (int argc, char ** argv)
|
||||
}
|
||||
}
|
||||
|
||||
ImprovedStaggeredFermion5DR Ds(Umu,Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,c1,c2,u0,params);
|
||||
ImprovedStaggeredFermion5DD Ds(Umu,Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,c1,c2,u0,params);
|
||||
|
||||
std::cout<<GridLogMessage<<"=========================================================="<<std::endl;
|
||||
std::cout<<GridLogMessage<<"= Testing Dhop against cshift implementation "<<std::endl;
|
||||
@ -288,7 +288,7 @@ int main (int argc, char ** argv)
|
||||
pickCheckerboard(Even,phi_e,phi);
|
||||
pickCheckerboard(Odd ,phi_o,phi);
|
||||
|
||||
SchurDiagMooeeOperator<ImprovedStaggeredFermion5DR,FermionField> HermOpEO(Ds);
|
||||
SchurDiagMooeeOperator<ImprovedStaggeredFermion5DD,FermionField> HermOpEO(Ds);
|
||||
HermOpEO.MpcDagMpc(chi_e,dchi_e);
|
||||
HermOpEO.MpcDagMpc(chi_o,dchi_o);
|
||||
|
||||
|
||||
@ -52,9 +52,9 @@ int main (int argc, char ** argv)
|
||||
pRNG.SeedFixedIntegers(seeds);
|
||||
// pRNG.SeedFixedIntegers(std::vector<int>({45,12,81,9});
|
||||
|
||||
typedef typename NaiveStaggeredFermionR::FermionField FermionField;
|
||||
typedef typename NaiveStaggeredFermionR::ComplexField ComplexField;
|
||||
typename NaiveStaggeredFermionR::ImplParams params;
|
||||
typedef typename NaiveStaggeredFermionD::FermionField FermionField;
|
||||
typedef typename NaiveStaggeredFermionD::ComplexField ComplexField;
|
||||
typename NaiveStaggeredFermionD::ImplParams params;
|
||||
|
||||
FermionField src (&Grid); random(pRNG,src);
|
||||
FermionField result(&Grid); result=Zero();
|
||||
@ -120,7 +120,7 @@ int main (int argc, char ** argv)
|
||||
// ref = ref + mass * src;
|
||||
}
|
||||
|
||||
NaiveStaggeredFermionR Ds(Umu,Grid,RBGrid,mass,c1,u0,params);
|
||||
NaiveStaggeredFermionD Ds(Umu,Grid,RBGrid,mass,c1,u0,params);
|
||||
|
||||
|
||||
std::cout<<GridLogMessage<<"=========================================================="<<std::endl;
|
||||
@ -258,7 +258,7 @@ int main (int argc, char ** argv)
|
||||
pickCheckerboard(Even,phi_e,phi);
|
||||
pickCheckerboard(Odd ,phi_o,phi);
|
||||
|
||||
SchurDiagMooeeOperator<NaiveStaggeredFermionR,FermionField> HermOpEO(Ds);
|
||||
SchurDiagMooeeOperator<NaiveStaggeredFermionD,FermionField> HermOpEO(Ds);
|
||||
HermOpEO.MpcDagMpc(chi_e,dchi_e);
|
||||
HermOpEO.MpcDagMpc(chi_o,dchi_o);
|
||||
|
||||
|
||||
@ -52,8 +52,8 @@ int main(int argc, char **argv)
|
||||
pRNG.SeedFixedIntegers(seeds);
|
||||
// pRNG.SeedFixedIntegers(std::vector<int>({45,12,81,9});
|
||||
|
||||
typedef typename WilsonCloverFermionR::FermionField FermionField;
|
||||
typename WilsonCloverFermionR::ImplParams params;
|
||||
typedef typename WilsonCloverFermionD::FermionField FermionField;
|
||||
typename WilsonCloverFermionD::ImplParams params;
|
||||
WilsonAnisotropyCoefficients anis;
|
||||
|
||||
FermionField src(&Grid);
|
||||
@ -88,8 +88,8 @@ int main(int argc, char **argv)
|
||||
RealD csw_r = 1.0;
|
||||
RealD csw_t = 1.0;
|
||||
|
||||
WilsonCloverFermionR Dwc(Umu, Grid, RBGrid, mass, csw_r, csw_t, anis, params);
|
||||
CompactWilsonCloverFermionR Dwc_compact(Umu, Grid, RBGrid, mass, csw_r, csw_t, 1.0, anis, params);
|
||||
WilsonCloverFermionD Dwc(Umu, Grid, RBGrid, mass, csw_r, csw_t, anis, params);
|
||||
CompactWilsonCloverFermionD Dwc_compact(Umu, Grid, RBGrid, mass, csw_r, csw_t, 1.0, anis, params);
|
||||
|
||||
std::cout << GridLogMessage << "==========================================================" << std::endl;
|
||||
std::cout << GridLogMessage << "= Testing that Deo + Doe = Dunprec " << std::endl;
|
||||
@ -324,8 +324,8 @@ int main(int argc, char **argv)
|
||||
}
|
||||
/////////////////
|
||||
|
||||
WilsonCloverFermionR Dwc_prime(U_prime, Grid, RBGrid, mass, csw_r, csw_t, anis, params);
|
||||
CompactWilsonCloverFermionR Dwc_compact_prime(U_prime, Grid, RBGrid, mass, csw_r, csw_t, 1.0, anis, params);
|
||||
WilsonCloverFermionD Dwc_prime(U_prime, Grid, RBGrid, mass, csw_r, csw_t, anis, params);
|
||||
CompactWilsonCloverFermionD Dwc_compact_prime(U_prime, Grid, RBGrid, mass, csw_r, csw_t, 1.0, anis, params);
|
||||
|
||||
tmp = Omega * src;
|
||||
pickCheckerboard(Even, src_e, tmp);
|
||||
@ -377,14 +377,14 @@ int main(int argc, char **argv)
|
||||
chi = Zero();
|
||||
phi = Zero();
|
||||
|
||||
WilsonFermionR Dw(Umu, Grid, RBGrid, mass, params);
|
||||
WilsonFermionD Dw(Umu, Grid, RBGrid, mass, params);
|
||||
|
||||
Dw.M(src, result);
|
||||
Dwc.M(src, chi);
|
||||
|
||||
Dwc_prime.M(Omega * src, phi);
|
||||
|
||||
WilsonFermionR Dw_prime(U_prime, Grid, RBGrid, mass, params);
|
||||
WilsonFermionD Dw_prime(U_prime, Grid, RBGrid, mass, params);
|
||||
Dw_prime.M(Omega * src, result2);
|
||||
|
||||
err = result - adj(Omega) * result2;
|
||||
@ -411,7 +411,7 @@ int main(int argc, char **argv)
|
||||
chi = Zero();
|
||||
phi = Zero();
|
||||
err = Zero();
|
||||
WilsonCloverFermionR Dwc_csw0(Umu, Grid, RBGrid, mass, 0.0, 0.0, anis, params); // <-- Notice: csw=0
|
||||
WilsonCloverFermionD Dwc_csw0(Umu, Grid, RBGrid, mass, 0.0, 0.0, anis, params); // <-- Notice: csw=0
|
||||
|
||||
pickCheckerboard(Even, phi_e, phi);
|
||||
pickCheckerboard(Odd, phi_o, phi);
|
||||
@ -437,7 +437,7 @@ int main(int argc, char **argv)
|
||||
chi = Zero();
|
||||
phi = Zero();
|
||||
err = Zero();
|
||||
CompactWilsonCloverFermionR Dwc_compact_csw0(Umu, Grid, RBGrid, mass, 0.0, 0.0, 1.0, anis, params); // <-- Notice: csw=0
|
||||
CompactWilsonCloverFermionD Dwc_compact_csw0(Umu, Grid, RBGrid, mass, 0.0, 0.0, 1.0, anis, params); // <-- Notice: csw=0
|
||||
|
||||
pickCheckerboard(Even, phi_e, phi);
|
||||
pickCheckerboard(Odd, phi_o, phi);
|
||||
|
||||
@ -74,7 +74,7 @@ int main (int argc, char ** argv)
|
||||
SU<Nc>::HotConfiguration(RNG4,Umu);
|
||||
}
|
||||
|
||||
typename WilsonCloverFermionR::ImplParams params;
|
||||
typename WilsonCloverFermionD::ImplParams params;
|
||||
WilsonAnisotropyCoefficients anis;
|
||||
RealD mass = 0.1;
|
||||
RealD csw_r = 1.0;
|
||||
@ -83,32 +83,32 @@ int main (int argc, char ** argv)
|
||||
std::cout<<GridLogMessage <<"=================================="<<std::endl;
|
||||
std::cout<<GridLogMessage <<"WilsonFermion test"<<std::endl;
|
||||
std::cout<<GridLogMessage <<"=================================="<<std::endl;
|
||||
WilsonFermionR Dw(Umu,*UGrid,*UrbGrid,mass,params);
|
||||
TestConserved<WilsonFermionR>(Dw,Umu,UGrid,UrbGrid,&RNG4);
|
||||
WilsonFermionD Dw(Umu,*UGrid,*UrbGrid,mass,params);
|
||||
TestConserved<WilsonFermionD>(Dw,Umu,UGrid,UrbGrid,&RNG4);
|
||||
|
||||
std::cout<<GridLogMessage <<"=================================="<<std::endl;
|
||||
std::cout<<GridLogMessage <<"WilsonCloverFermion test"<<std::endl;
|
||||
std::cout<<GridLogMessage <<"=================================="<<std::endl;
|
||||
WilsonCloverFermionR Dwc(Umu, *UGrid, *UrbGrid, mass, csw_r, csw_t, anis, params);
|
||||
TestConserved<WilsonCloverFermionR>(Dwc,Umu,UGrid,UrbGrid,&RNG4);
|
||||
WilsonCloverFermionD Dwc(Umu, *UGrid, *UrbGrid, mass, csw_r, csw_t, anis, params);
|
||||
TestConserved<WilsonCloverFermionD>(Dwc,Umu,UGrid,UrbGrid,&RNG4);
|
||||
|
||||
std::cout<<GridLogMessage <<"=================================="<<std::endl;
|
||||
std::cout<<GridLogMessage <<"CompactWilsonCloverFermion test"<<std::endl;
|
||||
std::cout<<GridLogMessage <<"=================================="<<std::endl;
|
||||
CompactWilsonCloverFermionR Dwcc(Umu, *UGrid, *UrbGrid, mass, csw_r, csw_t, 1.0, anis, params);
|
||||
TestConserved<CompactWilsonCloverFermionR>(Dwcc,Umu,UGrid,UrbGrid,&RNG4);
|
||||
CompactWilsonCloverFermionD Dwcc(Umu, *UGrid, *UrbGrid, mass, csw_r, csw_t, 1.0, anis, params);
|
||||
TestConserved<CompactWilsonCloverFermionD>(Dwcc,Umu,UGrid,UrbGrid,&RNG4);
|
||||
|
||||
std::cout<<GridLogMessage <<"=================================="<<std::endl;
|
||||
std::cout<<GridLogMessage <<"WilsonExpCloverFermion test"<<std::endl;
|
||||
std::cout<<GridLogMessage <<"=================================="<<std::endl;
|
||||
WilsonExpCloverFermionR Dewc(Umu, *UGrid, *UrbGrid, mass, csw_r, csw_t, anis, params);
|
||||
TestConserved<WilsonExpCloverFermionR>(Dewc,Umu,UGrid,UrbGrid,&RNG4);
|
||||
WilsonExpCloverFermionD Dewc(Umu, *UGrid, *UrbGrid, mass, csw_r, csw_t, anis, params);
|
||||
TestConserved<WilsonExpCloverFermionD>(Dewc,Umu,UGrid,UrbGrid,&RNG4);
|
||||
|
||||
std::cout<<GridLogMessage <<"=================================="<<std::endl;
|
||||
std::cout<<GridLogMessage <<"CompactWilsonExpCloverFermion test"<<std::endl;
|
||||
std::cout<<GridLogMessage <<"=================================="<<std::endl;
|
||||
CompactWilsonExpCloverFermionR Dewcc(Umu, *UGrid, *UrbGrid, mass, csw_r, csw_t, 1.0, anis, params);
|
||||
TestConserved<CompactWilsonExpCloverFermionR>(Dewcc,Umu,UGrid,UrbGrid,&RNG4);
|
||||
CompactWilsonExpCloverFermionD Dewcc(Umu, *UGrid, *UrbGrid, mass, csw_r, csw_t, 1.0, anis, params);
|
||||
TestConserved<CompactWilsonExpCloverFermionD>(Dewcc,Umu,UGrid,UrbGrid,&RNG4);
|
||||
|
||||
Grid_finalize();
|
||||
}
|
||||
|
||||
@ -89,7 +89,7 @@ int main (int argc, char ** argv)
|
||||
|
||||
RealD mass=0.1;
|
||||
|
||||
WilsonFermionR Dw(Umu,Grid,RBGrid,mass);
|
||||
WilsonFermionD Dw(Umu,Grid,RBGrid,mass);
|
||||
|
||||
LatticeFermion src_e (&RBGrid);
|
||||
LatticeFermion src_o (&RBGrid);
|
||||
@ -205,7 +205,7 @@ int main (int argc, char ** argv)
|
||||
pickCheckerboard(Even,phi_e,phi);
|
||||
pickCheckerboard(Odd ,phi_o,phi);
|
||||
|
||||
SchurDiagMooeeOperator<WilsonFermionR,LatticeFermion> HermOpEO(Dw);
|
||||
SchurDiagMooeeOperator<WilsonFermionD,LatticeFermion> HermOpEO(Dw);
|
||||
HermOpEO.MpcDagMpc(chi_e,dchi_e);
|
||||
HermOpEO.MpcDagMpc(chi_o,dchi_o);
|
||||
|
||||
|
||||
@ -52,8 +52,8 @@ int main(int argc, char **argv)
|
||||
pRNG.SeedFixedIntegers(seeds);
|
||||
// pRNG.SeedFixedIntegers(std::vector<int>({45,12,81,9});
|
||||
|
||||
typedef typename WilsonExpCloverFermionR::FermionField FermionField;
|
||||
typename WilsonExpCloverFermionR::ImplParams params;
|
||||
typedef typename WilsonExpCloverFermionD::FermionField FermionField;
|
||||
typename WilsonExpCloverFermionD::ImplParams params;
|
||||
WilsonAnisotropyCoefficients anis;
|
||||
|
||||
FermionField src(&Grid);
|
||||
@ -88,8 +88,8 @@ int main(int argc, char **argv)
|
||||
RealD csw_r = 1.0;
|
||||
RealD csw_t = 1.0;
|
||||
|
||||
WilsonExpCloverFermionR Dwc(Umu, Grid, RBGrid, mass, csw_r, csw_t, anis, params);
|
||||
CompactWilsonExpCloverFermionR Dwc_compact(Umu, Grid, RBGrid, mass, csw_r, csw_t, 1.0, anis, params);
|
||||
WilsonExpCloverFermionD Dwc(Umu, Grid, RBGrid, mass, csw_r, csw_t, anis, params);
|
||||
CompactWilsonExpCloverFermionD Dwc_compact(Umu, Grid, RBGrid, mass, csw_r, csw_t, 1.0, anis, params);
|
||||
|
||||
std::cout << GridLogMessage << "==========================================================" << std::endl;
|
||||
std::cout << GridLogMessage << "= Testing that Deo + Doe = Dunprec " << std::endl;
|
||||
@ -324,8 +324,8 @@ int main(int argc, char **argv)
|
||||
}
|
||||
/////////////////
|
||||
|
||||
WilsonExpCloverFermionR Dwc_prime(U_prime, Grid, RBGrid, mass, csw_r, csw_t, anis, params);
|
||||
CompactWilsonExpCloverFermionR Dwc_compact_prime(U_prime, Grid, RBGrid, mass, csw_r, csw_t, 1.0, anis, params);
|
||||
WilsonExpCloverFermionD Dwc_prime(U_prime, Grid, RBGrid, mass, csw_r, csw_t, anis, params);
|
||||
CompactWilsonExpCloverFermionD Dwc_compact_prime(U_prime, Grid, RBGrid, mass, csw_r, csw_t, 1.0, anis, params);
|
||||
|
||||
tmp = Omega * src;
|
||||
pickCheckerboard(Even, src_e, tmp);
|
||||
@ -377,14 +377,14 @@ int main(int argc, char **argv)
|
||||
chi = Zero();
|
||||
phi = Zero();
|
||||
|
||||
WilsonFermionR Dw(Umu, Grid, RBGrid, mass, params);
|
||||
WilsonFermionD Dw(Umu, Grid, RBGrid, mass, params);
|
||||
|
||||
Dw.M(src, result);
|
||||
Dwc.M(src, chi);
|
||||
|
||||
Dwc_prime.M(Omega * src, phi);
|
||||
|
||||
WilsonFermionR Dw_prime(U_prime, Grid, RBGrid, mass, params);
|
||||
WilsonFermionD Dw_prime(U_prime, Grid, RBGrid, mass, params);
|
||||
Dw_prime.M(Omega * src, result2);
|
||||
|
||||
err = result - adj(Omega) * result2;
|
||||
@ -411,7 +411,7 @@ int main(int argc, char **argv)
|
||||
chi = Zero();
|
||||
phi = Zero();
|
||||
err = Zero();
|
||||
WilsonExpCloverFermionR Dwc_csw0(Umu, Grid, RBGrid, mass, 0.0, 0.0, anis, params); // <-- Notice: csw=0
|
||||
WilsonExpCloverFermionD Dwc_csw0(Umu, Grid, RBGrid, mass, 0.0, 0.0, anis, params); // <-- Notice: csw=0
|
||||
|
||||
pickCheckerboard(Even, phi_e, phi);
|
||||
pickCheckerboard(Odd, phi_o, phi);
|
||||
@ -437,7 +437,7 @@ int main(int argc, char **argv)
|
||||
chi = Zero();
|
||||
phi = Zero();
|
||||
err = Zero();
|
||||
CompactWilsonExpCloverFermionR Dwc_compact_csw0(Umu, Grid, RBGrid, mass, 0.0, 0.0, 1.0, anis, params); // <-- Notice: csw=0
|
||||
CompactWilsonExpCloverFermionD Dwc_compact_csw0(Umu, Grid, RBGrid, mass, 0.0, 0.0, 1.0, anis, params); // <-- Notice: csw=0
|
||||
|
||||
pickCheckerboard(Even, phi_e, phi);
|
||||
pickCheckerboard(Odd, phi_o, phi);
|
||||
|
||||
@ -90,7 +90,7 @@ int main (int argc, char ** argv)
|
||||
|
||||
RealD mass=0.1;
|
||||
RealD mu = 0.1;
|
||||
WilsonTMFermionR Dw(Umu,Grid,RBGrid,mass,mu);
|
||||
WilsonTMFermionD Dw(Umu,Grid,RBGrid,mass,mu);
|
||||
|
||||
LatticeFermion src_e (&RBGrid);
|
||||
LatticeFermion src_o (&RBGrid);
|
||||
@ -206,7 +206,7 @@ int main (int argc, char ** argv)
|
||||
pickCheckerboard(Even,phi_e,phi);
|
||||
pickCheckerboard(Odd ,phi_o,phi);
|
||||
|
||||
SchurDiagMooeeOperator<WilsonTMFermionR,LatticeFermion> HermOpEO(Dw);
|
||||
SchurDiagMooeeOperator<WilsonTMFermionD,LatticeFermion> HermOpEO(Dw);
|
||||
HermOpEO.MpcDagMpc(chi_e,dchi_e);
|
||||
HermOpEO.MpcDagMpc(chi_o,dchi_o);
|
||||
|
||||
|
||||
@ -123,7 +123,7 @@ int main (int argc, char ** argv)
|
||||
RealD _mass,RealD _M5,
|
||||
std::vector<ComplexD> &gamma, RealD b,RealD c,const ImplParams &p= ImplParams()) :
|
||||
*/
|
||||
ZMobiusFermionR Ddwf(Umu, *FGrid, *FrbGrid, *UGrid, *UrbGrid, mass, M5, omegas,RealD(1.),RealD(0.));
|
||||
ZMobiusFermionD Ddwf(Umu, *FGrid, *FrbGrid, *UGrid, *UrbGrid, mass, M5, omegas,RealD(1.),RealD(0.));
|
||||
|
||||
LatticeFermion src_e (FrbGrid);
|
||||
LatticeFermion src_o (FrbGrid);
|
||||
@ -278,7 +278,7 @@ int main (int argc, char ** argv)
|
||||
pickCheckerboard(Odd ,phi_o,phi);
|
||||
|
||||
|
||||
SchurDiagMooeeOperator<ZMobiusFermionR,LatticeFermion> HermOpEO(Ddwf);
|
||||
SchurDiagMooeeOperator<ZMobiusFermionD,LatticeFermion> HermOpEO(Ddwf);
|
||||
HermOpEO.MpcDagMpc(chi_e,dchi_e);
|
||||
HermOpEO.MpcDagMpc(chi_o,dchi_o);
|
||||
|
||||
|
||||
Reference in New Issue
Block a user