portelli/Grid
1
0
Fork 0
mirror of https://github.com/paboyle/Grid.git synced 2025-11-04 14:04:32 +00:00
Code Releases Activity

Sizable improvement in multigrid for unsquared.

Browse Source 
6000 matmuls CG unprec
2000 matmuls CG prec (4000 eo muls)
1050 matmuls PGCR on 16^3 x 32 x 8 m=.01

Substantial effort on timing and logging infrastructure
This commit is contained in:
Peter Boyle
2015-07-24 01:31:13 +09:00
parent 5b475d5e08
commit 5e370db6c5
67 changed files with 945 additions and 753 deletions
Download Patch File Download Diff File Expand all files Collapse all files

22
benchmarks/Benchmark_comms.cc
View File

@@ -11,15 +11,15 @@ int main (int argc, char ** argv)
std::vector<int> simd_layout = GridDefaultSimd(Nd,vComplexD::Nsimd());
std::vector<int> mpi_layout = GridDefaultMpi();
int threads = GridThread::GetThreads();
std::cout << "Grid is setup to use "<<threads<<" threads"<<std::endl;
std::cout<<GridLogMessage << "Grid is setup to use "<<threads<<" threads"<<std::endl;
int Nloop=10;
int nmu=0;
for(int mu=0;mu<4;mu++) if (mpi_layout[mu]>1) nmu++;
std::cout << "===================================================================================================="<<std::endl;
std::cout << "= Benchmarking concurrent halo exchange in "<<nmu<<" dimensions"<<std::endl;
std::cout << "===================================================================================================="<<std::endl;
std::cout << " L "<<"\t\t"<<" Ls "<<"\t\t"<<"bytes"<<"\t\t"<<"MB/s uni"<<"\t\t"<<"MB/s bidi"<<std::endl;
std::cout<<GridLogMessage << "===================================================================================================="<<std::endl;
std::cout<<GridLogMessage << "= Benchmarking concurrent halo exchange in "<<nmu<<" dimensions"<<std::endl;
std::cout<<GridLogMessage << "===================================================================================================="<<std::endl;
std::cout<<GridLogMessage << " L "<<"\t\t"<<" Ls "<<"\t\t"<<"bytes"<<"\t\t"<<"MB/s uni"<<"\t\t"<<"MB/s bidi"<<std::endl;
@@ -87,15 +87,15 @@ int main (int argc, char ** argv)
double time = stop-start; // microseconds
std::cout << lat<<"\t\t"<<Ls<<"\t\t"<<bytes<<"\t\t"<<xbytes/time<<"\t\t"<<bidibytes/time<<std::endl;
std::cout<<GridLogMessage << lat<<"\t\t"<<Ls<<"\t\t"<<bytes<<"\t\t"<<xbytes/time<<"\t\t"<<bidibytes/time<<std::endl;
}
}
std::cout << "===================================================================================================="<<std::endl;
std::cout << "= Benchmarking sequential halo exchange in "<<nmu<<" dimensions"<<std::endl;
std::cout << "===================================================================================================="<<std::endl;
std::cout << " L "<<"\t\t"<<" Ls "<<"\t\t"<<"bytes"<<"\t\t"<<"MB/s uni"<<"\t\t"<<"MB/s bidi"<<std::endl;
std::cout<<GridLogMessage << "===================================================================================================="<<std::endl;
std::cout<<GridLogMessage << "= Benchmarking sequential halo exchange in "<<nmu<<" dimensions"<<std::endl;
std::cout<<GridLogMessage << "===================================================================================================="<<std::endl;
std::cout<<GridLogMessage << " L "<<"\t\t"<<" Ls "<<"\t\t"<<"bytes"<<"\t\t"<<"MB/s uni"<<"\t\t"<<"MB/s bidi"<<std::endl;
for(int lat=4;lat<=32;lat+=2){
@@ -163,7 +163,7 @@ int main (int argc, char ** argv)
double time = stop-start;
std::cout << lat<<"\t\t"<<Ls<<"\t\t"<<bytes<<"\t\t"<<xbytes/time<<"\t\t"<<bidibytes/time<<std::endl;
std::cout<<GridLogMessage << lat<<"\t\t"<<Ls<<"\t\t"<<bytes<<"\t\t"<<xbytes/time<<"\t\t"<<bidibytes/time<<std::endl;
}
}

38
benchmarks/Benchmark_dwf.cc
View File

@@ -21,7 +21,7 @@ int main (int argc, char ** argv)
Grid_init(&argc,&argv);
int threads = GridThread::GetThreads();
std::cout << "Grid is setup to use "<<threads<<" threads"<<std::endl;
std::cout<<GridLogMessage << "Grid is setup to use "<<threads<<" threads"<<std::endl;
std::vector<int> latt4 = GridDefaultLatt();
const int Ls=8;
@@ -81,7 +81,7 @@ int main (int argc, char ** argv)
RealD M5 =1.8;
DomainWallFermion Dw(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5);
std::cout << "Calling Dw"<<std::endl;
std::cout<<GridLogMessage << "Calling Dw"<<std::endl;
int ncall=10;
double t0=usecond();
for(int i=0;i<ncall;i++){
@@ -92,12 +92,12 @@ int main (int argc, char ** argv)
double volume=Ls; for(int mu=0;mu<Nd;mu++) volume=volume*latt4[mu];
double flops=1344*volume*ncall;
std::cout << "Called Dw"<<std::endl;
std::cout << "norm result "<< norm2(result)<<std::endl;
std::cout << "norm ref "<< norm2(ref)<<std::endl;
std::cout << "mflop/s = "<< flops/(t1-t0)<<std::endl;
std::cout<<GridLogMessage << "Called Dw"<<std::endl;
std::cout<<GridLogMessage << "norm result "<< norm2(result)<<std::endl;
std::cout<<GridLogMessage << "norm ref "<< norm2(ref)<<std::endl;
std::cout<<GridLogMessage << "mflop/s = "<< flops/(t1-t0)<<std::endl;
err = ref-result;
std::cout << "norm diff "<< norm2(err)<<std::endl;
std::cout<<GridLogMessage << "norm diff "<< norm2(err)<<std::endl;
if (1)
@@ -120,11 +120,11 @@ int main (int argc, char ** argv)
ref = -0.5*ref;
}
Dw.Dhop(src,result,1);
std::cout << "Called DwDag"<<std::endl;
std::cout << "norm result "<< norm2(result)<<std::endl;
std::cout << "norm ref "<< norm2(ref)<<std::endl;
std::cout<<GridLogMessage << "Called DwDag"<<std::endl;
std::cout<<GridLogMessage << "norm result "<< norm2(result)<<std::endl;
std::cout<<GridLogMessage << "norm ref "<< norm2(ref)<<std::endl;
err = ref-result;
std::cout << "norm diff "<< norm2(err)<<std::endl;
std::cout<<GridLogMessage << "norm diff "<< norm2(err)<<std::endl;
LatticeFermion src_e (FrbGrid);
LatticeFermion src_o (FrbGrid);
@@ -133,24 +133,32 @@ int main (int argc, char ** argv)
LatticeFermion r_eo (FGrid);
std::cout << "Calling Deo and Doe"<<std::endl;
std::cout<<GridLogMessage << "Calling Deo and Doe"<<std::endl;
pickCheckerboard(Even,src_e,src);
pickCheckerboard(Odd,src_o,src);
std::cout<<GridLogMessage << "src_e"<<norm2(src_e)<<std::endl;
std::cout<<GridLogMessage << "src_o"<<norm2(src_o)<<std::endl;
Dw.DhopEO(src_o,r_e,DaggerNo);
Dw.DhopOE(src_e,r_o,DaggerNo);
Dw.Dhop (src ,result,DaggerNo);
std::cout<<GridLogMessage << "r_e"<<norm2(r_e)<<std::endl;
std::cout<<GridLogMessage << "r_o"<<norm2(r_o)<<std::endl;
std::cout<<GridLogMessage << "res"<<norm2(result)<<std::endl;
setCheckerboard(r_eo,r_o);
setCheckerboard(r_eo,r_e);
err = r_eo-result;
std::cout << "norm diff "<< norm2(err)<<std::endl;
std::cout<<GridLogMessage << "norm diff "<< norm2(err)<<std::endl;
pickCheckerboard(Even,src_e,err);
pickCheckerboard(Odd,src_o,err);
std::cout << "norm diff even "<< norm2(src_e)<<std::endl;
std::cout << "norm diff odd "<< norm2(src_o)<<std::endl;
std::cout<<GridLogMessage << "norm diff even "<< norm2(src_e)<<std::endl;
std::cout<<GridLogMessage << "norm diff odd "<< norm2(src_o)<<std::endl;
Grid_finalize();
}

48
benchmarks/Benchmark_memory_bandwidth.cc
View File

@@ -17,13 +17,13 @@ int main (int argc, char ** argv)
std::vector<int> mpi_layout = GridDefaultMpi();
int threads = GridThread::GetThreads();
std::cout << "Grid is setup to use "<<threads<<" threads"<<std::endl;
std::cout<<GridLogMessage << "Grid is setup to use "<<threads<<" threads"<<std::endl;
std::cout << "===================================================================================================="<<std::endl;
std::cout << "= Benchmarking fused AXPY bandwidth ; sizeof(Real) "<<sizeof(Real)<<std::endl;
std::cout << "===================================================================================================="<<std::endl;
std::cout << " L "<<"\t\t"<<"bytes"<<"\t\t\t"<<"GB/s"<<"\t\t"<<"Gflop/s"<<std::endl;
std::cout << "----------------------------------------------------------"<<std::endl;
std::cout<<GridLogMessage << "===================================================================================================="<<std::endl;
std::cout<<GridLogMessage << "= Benchmarking fused AXPY bandwidth ; sizeof(Real) "<<sizeof(Real)<<std::endl;
std::cout<<GridLogMessage << "===================================================================================================="<<std::endl;
std::cout<<GridLogMessage << " L "<<"\t\t"<<"bytes"<<"\t\t\t"<<"GB/s"<<"\t\t"<<"Gflop/s"<<std::endl;
std::cout<<GridLogMessage << "----------------------------------------------------------"<<std::endl;
for(int lat=4;lat<=32;lat+=4){
@@ -49,15 +49,15 @@ int main (int argc, char ** argv)
double flops=vol*Nvec*2;// mul,add
double bytes=3*vol*Nvec*sizeof(Real);
std::cout<<std::setprecision(3) << lat<<"\t\t"<<bytes<<" \t\t"<<bytes/time<<"\t\t"<<flops/time<<std::endl;
std::cout<<GridLogMessage<<std::setprecision(3) << lat<<"\t\t"<<bytes<<" \t\t"<<bytes/time<<"\t\t"<<flops/time<<std::endl;
}
std::cout << "===================================================================================================="<<std::endl;
std::cout << "= Benchmarking a*x + y bandwidth"<<std::endl;
std::cout << "===================================================================================================="<<std::endl;
std::cout << " L "<<"\t\t"<<"bytes"<<"\t\t\t"<<"GB/s"<<"\t\t"<<"Gflop/s"<<std::endl;
std::cout << "----------------------------------------------------------"<<std::endl;
std::cout<<GridLogMessage << "===================================================================================================="<<std::endl;
std::cout<<GridLogMessage << "= Benchmarking a*x + y bandwidth"<<std::endl;
std::cout<<GridLogMessage << "===================================================================================================="<<std::endl;
std::cout<<GridLogMessage << " L "<<"\t\t"<<"bytes"<<"\t\t\t"<<"GB/s"<<"\t\t"<<"Gflop/s"<<std::endl;
std::cout<<GridLogMessage << "----------------------------------------------------------"<<std::endl;
for(int lat=4;lat<=32;lat+=4){
@@ -81,14 +81,14 @@ int main (int argc, char ** argv)
double flops=vol*Nvec*2;// mul,add
double bytes=3*vol*Nvec*sizeof(Real);
std::cout<<std::setprecision(3) << lat<<"\t\t"<<bytes<<" \t\t"<<bytes/time<<"\t\t"<<flops/time<<std::endl;
std::cout<<GridLogMessage<<std::setprecision(3) << lat<<"\t\t"<<bytes<<" \t\t"<<bytes/time<<"\t\t"<<flops/time<<std::endl;
}
std::cout << "===================================================================================================="<<std::endl;
std::cout << "= Benchmarking SCALE bandwidth"<<std::endl;
std::cout << "===================================================================================================="<<std::endl;
std::cout << " L "<<"\t\t"<<"bytes"<<"\t\t\t"<<"GB/s"<<"\t\t"<<"Gflop/s"<<std::endl;
std::cout<<GridLogMessage << "===================================================================================================="<<std::endl;
std::cout<<GridLogMessage << "= Benchmarking SCALE bandwidth"<<std::endl;
std::cout<<GridLogMessage << "===================================================================================================="<<std::endl;
std::cout<<GridLogMessage << " L "<<"\t\t"<<"bytes"<<"\t\t\t"<<"GB/s"<<"\t\t"<<"Gflop/s"<<std::endl;
for(int lat=4;lat<=32;lat+=4){
@@ -114,15 +114,15 @@ int main (int argc, char ** argv)
double bytes=2*vol*Nvec*sizeof(Real);
double flops=vol*Nvec*1;// mul
std::cout <<std::setprecision(3) << lat<<"\t\t"<<bytes<<" \t\t"<<bytes/time<<"\t\t"<<flops/time<<std::endl;
std::cout<<GridLogMessage <<std::setprecision(3) << lat<<"\t\t"<<bytes<<" \t\t"<<bytes/time<<"\t\t"<<flops/time<<std::endl;
}
std::cout << "===================================================================================================="<<std::endl;
std::cout << "= Benchmarking READ bandwidth"<<std::endl;
std::cout << "===================================================================================================="<<std::endl;
std::cout << " L "<<"\t\t"<<"bytes"<<"\t\t\t"<<"GB/s"<<"\t\t"<<"Gflop/s"<<std::endl;
std::cout << "----------------------------------------------------------"<<std::endl;
std::cout<<GridLogMessage << "===================================================================================================="<<std::endl;
std::cout<<GridLogMessage << "= Benchmarking READ bandwidth"<<std::endl;
std::cout<<GridLogMessage << "===================================================================================================="<<std::endl;
std::cout<<GridLogMessage << " L "<<"\t\t"<<"bytes"<<"\t\t\t"<<"GB/s"<<"\t\t"<<"Gflop/s"<<std::endl;
std::cout<<GridLogMessage << "----------------------------------------------------------"<<std::endl;
for(int lat=4;lat<=32;lat+=4){
@@ -147,7 +147,7 @@ int main (int argc, char ** argv)
double bytes=vol*Nvec*sizeof(Real);
double flops=vol*Nvec*2;// mul,add
std::cout<<std::setprecision(3) << lat<<"\t\t"<<bytes<<" \t\t"<<bytes/time<<"\t\t"<<flops/time<<std::endl;
std::cout<<GridLogMessage<<std::setprecision(3) << lat<<"\t\t"<<bytes<<" \t\t"<<bytes/time<<"\t\t"<<flops/time<<std::endl;
}

50
benchmarks/Benchmark_su3.cc
View File

@@ -14,13 +14,13 @@ int main (int argc, char ** argv)
std::vector<int> mpi_layout = GridDefaultMpi();
int threads = GridThread::GetThreads();
std::cout << "Grid is setup to use "<<threads<<" threads"<<std::endl;
std::cout<<GridLogMessage << "Grid is setup to use "<<threads<<" threads"<<std::endl;
std::cout << "===================================================================================================="<<std::endl;
std::cout << "= Benchmarking SU3xSU3 x= x*y"<<std::endl;
std::cout << "===================================================================================================="<<std::endl;
std::cout << " L "<<"\t\t"<<"bytes"<<"\t\t\t"<<"GB/s\t\t GFlop/s"<<std::endl;
std::cout << "----------------------------------------------------------"<<std::endl;
std::cout<<GridLogMessage << "===================================================================================================="<<std::endl;
std::cout<<GridLogMessage << "= Benchmarking SU3xSU3 x= x*y"<<std::endl;
std::cout<<GridLogMessage << "===================================================================================================="<<std::endl;
std::cout<<GridLogMessage << " L "<<"\t\t"<<"bytes"<<"\t\t\t"<<"GB/s\t\t GFlop/s"<<std::endl;
std::cout<<GridLogMessage << "----------------------------------------------------------"<<std::endl;
for(int lat=2;lat<=24;lat+=2){
@@ -43,16 +43,16 @@ int main (int argc, char ** argv)
double bytes=3.0*vol*Nc*Nc*sizeof(Complex);
double footprint=2.0*vol*Nc*Nc*sizeof(Complex);
double flops=Nc*Nc*(6.0+8.0+8.0)*vol;
std::cout<<std::setprecision(3) << lat<<"\t\t"<<footprint<<" \t\t"<<bytes/time<<"\t\t" << flops/time<<std::endl;
std::cout<<GridLogMessage<<std::setprecision(3) << lat<<"\t\t"<<footprint<<" \t\t"<<bytes/time<<"\t\t" << flops/time<<std::endl;
}
std::cout << "===================================================================================================="<<std::endl;
std::cout << "= Benchmarking SU3xSU3 z= x*y"<<std::endl;
std::cout << "===================================================================================================="<<std::endl;
std::cout << " L "<<"\t\t"<<"bytes"<<"\t\t\t"<<"GB/s\t\t GFlop/s"<<std::endl;
std::cout << "----------------------------------------------------------"<<std::endl;
std::cout<<GridLogMessage << "===================================================================================================="<<std::endl;
std::cout<<GridLogMessage << "= Benchmarking SU3xSU3 z= x*y"<<std::endl;
std::cout<<GridLogMessage << "===================================================================================================="<<std::endl;
std::cout<<GridLogMessage << " L "<<"\t\t"<<"bytes"<<"\t\t\t"<<"GB/s\t\t GFlop/s"<<std::endl;
std::cout<<GridLogMessage << "----------------------------------------------------------"<<std::endl;
for(int lat=2;lat<=24;lat+=2){
@@ -75,15 +75,15 @@ int main (int argc, char ** argv)
double bytes=3*vol*Nc*Nc*sizeof(Complex);
double flops=Nc*Nc*(6+8+8)*vol;
std::cout<<std::setprecision(3) << lat<<"\t\t"<<bytes<<" \t\t"<<bytes/time<<"\t\t" << flops/time<<std::endl;
std::cout<<GridLogMessage<<std::setprecision(3) << lat<<"\t\t"<<bytes<<" \t\t"<<bytes/time<<"\t\t" << flops/time<<std::endl;
}
std::cout << "===================================================================================================="<<std::endl;
std::cout << "= Benchmarking SU3xSU3 mult(z,x,y)"<<std::endl;
std::cout << "===================================================================================================="<<std::endl;
std::cout << " L "<<"\t\t"<<"bytes"<<"\t\t\t"<<"GB/s\t\t GFlop/s"<<std::endl;
std::cout << "----------------------------------------------------------"<<std::endl;
std::cout<<GridLogMessage << "===================================================================================================="<<std::endl;
std::cout<<GridLogMessage << "= Benchmarking SU3xSU3 mult(z,x,y)"<<std::endl;
std::cout<<GridLogMessage << "===================================================================================================="<<std::endl;
std::cout<<GridLogMessage << " L "<<"\t\t"<<"bytes"<<"\t\t\t"<<"GB/s\t\t GFlop/s"<<std::endl;
std::cout<<GridLogMessage << "----------------------------------------------------------"<<std::endl;
for(int lat=2;lat<=24;lat+=2){
@@ -106,15 +106,15 @@ int main (int argc, char ** argv)
double bytes=3*vol*Nc*Nc*sizeof(Complex);
double flops=Nc*Nc*(6+8+8)*vol;
std::cout<<std::setprecision(3) << lat<<"\t\t"<<bytes<<" \t\t"<<bytes/time<<"\t\t" << flops/time<<std::endl;
std::cout<<GridLogMessage<<std::setprecision(3) << lat<<"\t\t"<<bytes<<" \t\t"<<bytes/time<<"\t\t" << flops/time<<std::endl;
}
std::cout << "===================================================================================================="<<std::endl;
std::cout << "= Benchmarking SU3xSU3 mac(z,x,y)"<<std::endl;
std::cout << "===================================================================================================="<<std::endl;
std::cout << " L "<<"\t\t"<<"bytes"<<"\t\t\t"<<"GB/s\t\t GFlop/s"<<std::endl;
std::cout << "----------------------------------------------------------"<<std::endl;
std::cout<<GridLogMessage << "===================================================================================================="<<std::endl;
std::cout<<GridLogMessage << "= Benchmarking SU3xSU3 mac(z,x,y)"<<std::endl;
std::cout<<GridLogMessage << "===================================================================================================="<<std::endl;
std::cout<<GridLogMessage << " L "<<"\t\t"<<"bytes"<<"\t\t\t"<<"GB/s\t\t GFlop/s"<<std::endl;
std::cout<<GridLogMessage << "----------------------------------------------------------"<<std::endl;
for(int lat=2;lat<=24;lat+=2){
@@ -137,7 +137,7 @@ int main (int argc, char ** argv)
double bytes=3*vol*Nc*Nc*sizeof(Complex);
double flops=Nc*Nc*(8+8+8)*vol;
std::cout<<std::setprecision(3) << lat<<"\t\t"<<bytes<<" \t\t"<<bytes/time<<"\t\t" << flops/time<<std::endl;
std::cout<<GridLogMessage<<std::setprecision(3) << lat<<"\t\t"<<bytes<<" \t\t"<<bytes/time<<"\t\t" << flops/time<<std::endl;
}

34
benchmarks/Benchmark_wilson.cc
View File

@@ -28,10 +28,10 @@ int main (int argc, char ** argv)
GridRedBlackCartesian RBGrid(latt_size,simd_layout,mpi_layout);
int threads = GridThread::GetThreads();
std::cout << "Grid is setup to use "<<threads<<" threads"<<std::endl;
std::cout << "Grid floating point word size is REALF"<< sizeof(RealF)<<std::endl;
std::cout << "Grid floating point word size is REALD"<< sizeof(RealD)<<std::endl;
std::cout << "Grid floating point word size is REAL"<< sizeof(Real)<<std::endl;
std::cout<<GridLogMessage << "Grid is setup to use "<<threads<<" threads"<<std::endl;
std::cout<<GridLogMessage << "Grid floating point word size is REALF"<< sizeof(RealF)<<std::endl;
std::cout<<GridLogMessage << "Grid floating point word size is REALD"<< sizeof(RealD)<<std::endl;
std::cout<<GridLogMessage << "Grid floating point word size is REAL"<< sizeof(Real)<<std::endl;
std::vector<int> seeds({1,2,3,4});
GridParallelRNG pRNG(&Grid);
@@ -58,8 +58,8 @@ int main (int argc, char ** argv)
for(int nn=0;nn<Nd;nn++){
random(pRNG,U[nn]);
if(0) {
if (nn==-1) { U[nn]=zero; std::cout << "zeroing gauge field in dir "<<nn<<std::endl; }
else { U[nn] = cone;std::cout << "unit gauge field in dir "<<nn<<std::endl; }
if (nn==-1) { U[nn]=zero; std::cout<<GridLogMessage << "zeroing gauge field in dir "<<nn<<std::endl; }
else { U[nn] = cone;std::cout<<GridLogMessage << "unit gauge field in dir "<<nn<<std::endl; }
}
pokeIndex<LorentzIndex>(Umu,U[nn],nn);
}
@@ -89,7 +89,7 @@ int main (int argc, char ** argv)
RealD mass=0.1;
WilsonFermion Dw(Umu,Grid,RBGrid,mass);
std::cout << "Calling Dw"<<std::endl;
std::cout<<GridLogMessage << "Calling Dw"<<std::endl;
int ncall=10000;
double t0=usecond();
for(int i=0;i<ncall;i++){
@@ -98,12 +98,12 @@ int main (int argc, char ** argv)
double t1=usecond();
double flops=1344*volume*ncall;
std::cout << "Called Dw"<<std::endl;
std::cout << "norm result "<< norm2(result)<<std::endl;
std::cout << "norm ref "<< norm2(ref)<<std::endl;
std::cout << "mflop/s = "<< flops/(t1-t0)<<std::endl;
std::cout<<GridLogMessage << "Called Dw"<<std::endl;
std::cout<<GridLogMessage << "norm result "<< norm2(result)<<std::endl;
std::cout<<GridLogMessage << "norm ref "<< norm2(ref)<<std::endl;
std::cout<<GridLogMessage << "mflop/s = "<< flops/(t1-t0)<<std::endl;
err = ref-result;
std::cout << "norm diff "<< norm2(err)<<std::endl;
std::cout<<GridLogMessage << "norm diff "<< norm2(err)<<std::endl;
// for(int ss=0;ss<10;ss++ ){
@@ -112,7 +112,7 @@ int main (int argc, char ** argv)
for(int j=0;j<Nc;j++){
ComplexF * ref_p = (ComplexF *)&ref._odata[ss]()(i)(j);
ComplexF * res_p = (ComplexF *)&result._odata[ss]()(i)(j);
std::cout << ss<< " "<<i<<" "<<j<<" "<< (*ref_p)<<" " <<(*res_p)<<std::endl;
std::cout<<GridLogMessage << ss<< " "<<i<<" "<<j<<" "<< (*ref_p)<<" " <<(*res_p)<<std::endl;
}
}
}
@@ -136,11 +136,11 @@ int main (int argc, char ** argv)
}
ref = -0.5*ref;
Dw.Dhop(src,result,1);
std::cout << "Called DwDag"<<std::endl;
std::cout << "norm result "<< norm2(result)<<std::endl;
std::cout << "norm ref "<< norm2(ref)<<std::endl;
std::cout<<GridLogMessage << "Called DwDag"<<std::endl;
std::cout<<GridLogMessage << "norm result "<< norm2(result)<<std::endl;
std::cout<<GridLogMessage << "norm ref "<< norm2(ref)<<std::endl;
err = ref-result;
std::cout << "norm diff "<< norm2(err)<<std::endl;
std::cout<<GridLogMessage << "norm diff "<< norm2(err)<<std::endl;
Grid_finalize();
}

4
gcc-bug-report/broken.cc
View File

@@ -29,12 +29,12 @@ public:
template<int N,class obj,typename std::enable_if<N==obj::NestLevel >::type * = nullptr > auto function(const obj &arg)-> obj
{
std::cout<<"Leaf "<<obj::NestLevel<<std::endl;
std::cout<<GridLogMessage<<"Leaf "<<obj::NestLevel<<std::endl;
return arg;
}
template<int N,class obj,typename std::enable_if<N!=obj::NestLevel >::type * = nullptr > auto function(const obj &arg)-> obj
{
std::cout<<"Node "<<obj::NestLevel<<std::endl;
std::cout<<GridLogMessage<<"Node "<<obj::NestLevel<<std::endl;
obj ret;
ret.internal=function<N>(arg.internal);
return ret;

3
lib/Grid.h
View File

@@ -35,6 +35,9 @@
#include <GridConfig.h>
#include <GridTime.h>
#include <GridLog.h>
////////////////////////////////////////////////////////////
// Tunable header includes
////////////////////////////////////////////////////////////

93
lib/GridInit.cc
View File

@@ -31,9 +31,11 @@ namespace Grid {
//////////////////////////////////////////////////////
static std::vector<int> Grid_default_latt;
static std::vector<int> Grid_default_mpi;
int GridThread::_threads;
const std::vector<int> &GridDefaultLatt(void) {return Grid_default_latt;};
const std::vector<int> &GridDefaultMpi(void) {return Grid_default_mpi;};
const std::vector<int> GridDefaultSimd(int dims,int nsimd)
{
std::vector<int> layout(dims);
@@ -50,10 +52,6 @@ namespace Grid {
return layout;
}
const std::vector<int> &GridDefaultLatt(void) {return Grid_default_latt;};
const std::vector<int> &GridDefaultMpi(void) {return Grid_default_mpi;};
////////////////////////////////////////////////////////////
// Command line parsing assist for stock controls
////////////////////////////////////////////////////////////
@@ -70,6 +68,23 @@ bool GridCmdOptionExists(char** begin, char** end, const std::string& option)
{
return std::find(begin, end, option) != end;
}
// Comma separated list
void GridCmdOptionCSL(std::string str,std::vector<std::string> & vec)
{
size_t pos = 0;
std::string token;
std::string delimiter(",");
vec.resize(0);
while ((pos = str.find(delimiter)) != std::string::npos) {
token = str.substr(0, pos);
vec.push_back(token);
str.erase(0, pos + delimiter.length());
}
token = str;
vec.push_back(token);
return;
}
void GridCmdOptionIntVector(std::string &str,std::vector<int> & vec)
{
@@ -84,6 +99,7 @@ void GridCmdOptionIntVector(std::string &str,std::vector<int> & vec)
return;
}
void GridParseLayout(char **argv,int argc,
std::vector<int> &latt,
std::vector<int> &mpi)
@@ -118,6 +134,7 @@ std::string GridCmdVectorIntToString(const std::vector<int> & vec){
return oss.str();
}
/////////////////////////////////////////////////////////
//
/////////////////////////////////////////////////////////
void Grid_init(int *argc,char ***argv)
{
@@ -126,15 +143,33 @@ void Grid_init(int *argc,char ***argv)
#endif
// Parse command line args.
GridLogger::StopWatch.Start();
std::string arg;
std::vector<std::string> logstreams;
std::string defaultLog("Error,Warning,Message,Performance");
GridCmdOptionCSL(defaultLog,logstreams);
GridLogConfigure(logstreams);
if( GridCmdOptionExists(*argv,*argv+*argc,"--help") ){
std::cout<<"--help : this message"<<std::endl;
std::cout<<"--debug-signals : catch sigsegv and print a blame report"<<std::endl;
std::cout<<"--debug-stdout : print stdout from EVERY node"<<std::endl;
std::cout<<"--decomposition : report on default omp,mpi and simd decomposition"<<std::endl;
std::cout<<"--mpi n.n.n.n : default MPI decomposition"<<std::endl;
std::cout<<"--omp n : default number of OMP threads"<<std::endl;
std::cout<<"--grid n.n.n.n : default Grid size"<<std::endl;
std::cout<<GridLogMessage<<"--help : this message"<<std::endl;
std::cout<<GridLogMessage<<"--debug-signals : catch sigsegv and print a blame report"<<std::endl;
std::cout<<GridLogMessage<<"--debug-stdout : print stdout from EVERY node"<<std::endl;
std::cout<<GridLogMessage<<"--decomposition : report on default omp,mpi and simd decomposition"<<std::endl;
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;
}
if( GridCmdOptionExists(*argv,*argv+*argc,"--log") ){
arg = GridCmdOptionPayload(*argv,*argv+*argc,"--log");
GridCmdOptionCSL(arg,logstreams);
GridLogConfigure(logstreams);
}
if( GridCmdOptionExists(*argv,*argv+*argc,"--debug-signals") ){
Grid_debug_handler_init();
}
@@ -152,36 +187,16 @@ void Grid_init(int *argc,char ***argv)
Grid_default_latt,
Grid_default_mpi);
if( GridCmdOptionExists(*argv,*argv+*argc,"--decomposition") ){
std::cout<<"Grid Decomposition\n";
std::cout<<"\tOpenMP threads : "<<GridThread::GetThreads()<<std::endl;
std::cout<<"\tMPI tasks : "<<GridCmdVectorIntToString(GridDefaultMpi())<<std::endl;
std::cout<<"\tvRealF : "<<sizeof(vRealF)*8 <<"bits ; " <<GridCmdVectorIntToString(GridDefaultSimd(4,vRealF::Nsimd()))<<std::endl;
std::cout<<"\tvRealD : "<<sizeof(vRealD)*8 <<"bits ; " <<GridCmdVectorIntToString(GridDefaultSimd(4,vRealD::Nsimd()))<<std::endl;
std::cout<<"\tvComplexF : "<<sizeof(vComplexF)*8 <<"bits ; " <<GridCmdVectorIntToString(GridDefaultSimd(4,vComplexF::Nsimd()))<<std::endl;
std::cout<<"\tvComplexD : "<<sizeof(vComplexD)*8 <<"bits ; " <<GridCmdVectorIntToString(GridDefaultSimd(4,vComplexD::Nsimd()))<<std::endl;
}
std::cout<<GridLogMessage<<"Grid Decomposition\n";
std::cout<<GridLogMessage<<"\tOpenMP threads : "<<GridThread::GetThreads()<<std::endl;
std::cout<<GridLogMessage<<"\tMPI tasks : "<<GridCmdVectorIntToString(GridDefaultMpi())<<std::endl;
std::cout<<GridLogMessage<<"\tvRealF : "<<sizeof(vRealF)*8 <<"bits ; " <<GridCmdVectorIntToString(GridDefaultSimd(4,vRealF::Nsimd()))<<std::endl;
std::cout<<GridLogMessage<<"\tvRealD : "<<sizeof(vRealD)*8 <<"bits ; " <<GridCmdVectorIntToString(GridDefaultSimd(4,vRealD::Nsimd()))<<std::endl;
std::cout<<GridLogMessage<<"\tvComplexF : "<<sizeof(vComplexF)*8 <<"bits ; " <<GridCmdVectorIntToString(GridDefaultSimd(4,vComplexF::Nsimd()))<<std::endl;
std::cout<<GridLogMessage<<"\tvComplexD : "<<sizeof(vComplexD)*8 <<"bits ; " <<GridCmdVectorIntToString(GridDefaultSimd(4,vComplexD::Nsimd()))<<std::endl;
}
////////////////////////////////////////////////////////////
// Verbose limiter on MPI tasks
////////////////////////////////////////////////////////////
void Grid_quiesce_nodes(void)
{
#ifdef GRID_COMMS_MPI
int me;
MPI_Comm_rank(MPI_COMM_WORLD,&me);
if ( me ) {
std::cout.setstate(std::ios::badbit);
}
#endif
}
void Grid_unquiesce_nodes(void)
{
#ifdef GRID_COMMS_MPI
std::cout.clear();
#endif
}

62
lib/GridLog.cc Normal file
View File

@@ -0,0 +1,62 @@
#include <Grid.h>
namespace Grid {
GridStopWatch GridLogger::StopWatch;
std::ostream GridLogger::devnull(0);
GridLogger GridLogError (1,"Error");
GridLogger GridLogWarning (1,"Warning");
GridLogger GridLogMessage (1,"Message");
GridLogger GridLogDebug (1,"Debug");
GridLogger GridLogPerformance(1,"Performance");
GridLogger GridLogIterative (1,"Iterative");
void GridLogConfigure(std::vector<std::string> &logstreams)
{
GridLogError.Active(0);
GridLogWarning.Active(0);
GridLogMessage.Active(0);
GridLogIterative.Active(0);
GridLogDebug.Active(0);
GridLogPerformance.Active(0);
for(int i=0;i<logstreams.size();i++){
if ( logstreams[i]== std::string("Error") ) GridLogError.Active(1);
if ( logstreams[i]== std::string("Warning") ) GridLogWarning.Active(1);
if ( logstreams[i]== std::string("Message") ) GridLogMessage.Active(1);
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);
}
}
////////////////////////////////////////////////////////////
// Verbose limiter on MPI tasks
////////////////////////////////////////////////////////////
void Grid_quiesce_nodes(void)
{
#ifdef GRID_COMMS_MPI
int me;
MPI_Comm_rank(MPI_COMM_WORLD,&me);
if ( me ) {
std::cout<<GridLogMessage.setstate(std::ios::badbit);
}
#endif
}
void Grid_unquiesce_nodes(void)
{
#ifdef GRID_COMMS_MPI
std::cout<<GridLogMessage.clear();
#endif
}
std::ostream& operator<< (std::ostream& stream, const GridTime& time)
{
stream << time.count()<<" ms";
return stream;
}
}

46
lib/GridLog.h Normal file
View File

@@ -0,0 +1,46 @@
#ifndef GRID_LOG_H
#define GRID_LOG_H
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 GridLogger {
int active;
std::string name;
public:
static GridStopWatch StopWatch;
static std::ostream devnull;
GridLogger(int on, std::string nm): active(on), name(nm) {
};
void Active(int on) {active = on;};
friend std::ostream& operator<< (std::ostream& stream, const GridLogger& log){
if ( log.active ) {
StopWatch.Stop();
GridTime now = StopWatch.Elapsed();
StopWatch.Start();
stream << "Grid : "<<log.name << " : " << now << " : ";
return stream;
} else {
return devnull;
}
}
};
void GridLogConfigure(std::vector<std::string> &logstreams);
extern GridLogger GridLogError;
extern GridLogger GridLogWarning;
extern GridLogger GridLogMessage;
extern GridLogger GridLogDebug ;
extern GridLogger GridLogPerformance;
extern GridLogger GridLogIterative ;
}
#endif

44
lib/GridTime.h Normal file
View File

@@ -0,0 +1,44 @@
#ifndef GRID_TIME_H
#define GRID_TIME_H
namespace Grid {
// Dress the output; use std::chrono
#include <chrono>
#include <ctime>
typedef std::chrono::system_clock GridClock;
typedef std::chrono::time_point<GridClock> GridTimePoint;
typedef std::chrono::milliseconds GridTime;
class GridStopWatch {
private:
bool running;
GridTimePoint start;
GridTime accumulator;
public:
GridStopWatch () {
Reset();
}
void Start(void) {
assert(running == false);
start = GridClock::now();
running = true;
}
void Stop(void) {
assert(running == true);
accumulator+= std::chrono::duration_cast<GridTime>(GridClock::now()-start);
running = false;
};
void Reset(void){
running = false;
start = GridClock::now();
accumulator = std::chrono::duration_cast<GridTime>(start-start);
}
GridTime Elapsed(void) {
assert(running == false);
return accumulator;
}
};
}
#endif

4
lib/Make.inc
View File

@@ -1,4 +1,4 @@
HFILES=./Cshift.h ./simd/Grid_avx.h ./simd/Grid_vector_types.h ./simd/Grid_sse4.h ./simd/Grid_avx512.h ./simd/Grid_empty.h ./simd/Grid_vector_unops.h ./simd/Grid_neon.h ./simd/Grid_qpx.h ./Tensors.h ./Algorithms.h ./communicator/Communicator_base.h ./lattice/Lattice_rng.h ./lattice/Lattice_reduction.h ./lattice/Lattice_transfer.h ./lattice/Lattice_unary.h ./lattice/Lattice_peekpoke.h ./lattice/Lattice_coordinate.h ./lattice/Lattice_comparison.h ./lattice/Lattice_overload.h ./lattice/Lattice_reality.h ./lattice/Lattice_local.h ./lattice/Lattice_conformable.h ./lattice/Lattice_where.h ./lattice/Lattice_comparison_utils.h ./lattice/Lattice_arith.h ./lattice/Lattice_base.h ./lattice/Lattice_ET.h ./lattice/Lattice_transpose.h ./lattice/Lattice_trace.h ./Stencil.h ./tensors/Tensor_arith_sub.h ./tensors/Tensor_exp.h ./tensors/Tensor_arith_mul.h ./tensors/Tensor_class.h ./tensors/Tensor_logical.h ./tensors/Tensor_transpose.h ./tensors/Tensor_arith_mac.h ./tensors/Tensor_arith_scalar.h ./tensors/Tensor_reality.h ./tensors/Tensor_trace.h ./tensors/Tensor_index.h ./tensors/Tensor_arith_add.h ./tensors/Tensor_outer.h ./tensors/Tensor_inner.h ./tensors/Tensor_traits.h ./tensors/Tensor_Ta.h ./tensors/Tensor_unary.h ./tensors/Tensor_determinant.h ./tensors/Tensor_arith.h ./tensors/Tensor_extract_merge.h ./Communicator.h ./Cartesian.h ./parallelIO/NerscIO.h ./qcd/QCD.h ./qcd/hmc/integrators/Integrator.h ./qcd/hmc/integrators/Integrator_algorithm.h ./qcd/hmc/HMC.h ./qcd/utils/SpaceTimeGrid.h ./qcd/utils/SUn.h ./qcd/utils/LinalgUtils.h ./qcd/utils/CovariantCshift.h ./qcd/utils/WilsonLoops.h ./qcd/action/ActionBase.h ./qcd/action/gauge/WilsonGaugeAction.h ./qcd/action/Actions.h ./qcd/action/fermion/CayleyFermion5D.h ./qcd/action/fermion/ScaledShamirFermion.h ./qcd/action/fermion/MobiusFermion.h ./qcd/action/fermion/OverlapWilsonContfracTanhFermion.h ./qcd/action/fermion/PartialFractionFermion5D.h ./qcd/action/fermion/ShamirZolotarevFermion.h ./qcd/action/fermion/FermionOperator.h ./qcd/action/fermion/WilsonFermion5D.h ./qcd/action/fermion/WilsonCompressor.h ./qcd/action/fermion/OverlapWilsonPartialFractionZolotarevFermion.h ./qcd/action/fermion/WilsonKernels.h ./qcd/action/fermion/DomainWallFermion.h ./qcd/action/fermion/OverlapWilsonPartialFractionTanhFermion.h ./qcd/action/fermion/OverlapWilsonContfracZolotarevFermion.h ./qcd/action/fermion/MobiusZolotarevFermion.h ./qcd/action/fermion/g5HermitianLinop.h ./qcd/action/fermion/OverlapWilsonCayleyTanhFermion.h ./qcd/action/fermion/WilsonFermion.h ./qcd/action/fermion/ContinuedFractionFermion5D.h ./qcd/action/fermion/OverlapWilsonCayleyZolotarevFermion.h ./qcd/spin/TwoSpinor.h ./qcd/spin/Dirac.h ./cshift/Cshift_common.h ./cshift/Cshift_none.h ./cshift/Cshift_mpi.h ./Simd.h ./GridConfig.h ./cartesian/Cartesian_base.h ./cartesian/Cartesian_red_black.h ./cartesian/Cartesian_full.h ./AlignedAllocator.h ./Lattice.h ./Old/Tensor_poke.h ./Old/Tensor_peek.h ./Threads.h ./Grid.h ./algorithms/Preconditioner.h ./algorithms/iterative/ConjugateResidual.h ./algorithms/iterative/PrecGeneralisedConjugateResidual.h ./algorithms/iterative/ConjugateGradientMultiShift.h ./algorithms/iterative/SchurRedBlack.h ./algorithms/iterative/NormalEquations.h ./algorithms/iterative/ConjugateGradient.h ./algorithms/iterative/AdefGeneric.h ./algorithms/approx/Chebyshev.h ./algorithms/approx/Zolotarev.h ./algorithms/approx/MultiShiftFunction.h ./algorithms/approx/bigfloat.h ./algorithms/approx/bigfloat_double.h ./algorithms/approx/Remez.h ./algorithms/LinearOperator.h ./algorithms/SparseMatrix.h ./algorithms/CoarsenedMatrix.h ./stencil/Lebesgue.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/BfmHDCG.h ./algorithms/iterative/ConjugateGradient.h ./algorithms/iterative/ConjugateGradientMultiShift.h ./algorithms/iterative/ConjugateResidual.h ./algorithms/iterative/GeneralisedMinimumResidual.h ./algorithms/iterative/gmres.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 ./GridConfig.h ./GridLog.h ./GridTime.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 ./Old/Tensor_peek.h ./Old/Tensor_poke.h ./parallelIO/NerscIO.h ./qcd/action/ActionBase.h ./qcd/action/Actions.h ./qcd/action/DiffAction.h ./qcd/action/fermion/CayleyFermion5D.h ./qcd/action/fermion/ContinuedFractionFermion5D.h ./qcd/action/fermion/DomainWallFermion.h ./qcd/action/fermion/FermionOperator.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/gauge/WilsonGaugeAction.h ./qcd/hmc/HMC.h ./qcd/hmc/integrators/Integrator.h ./qcd/hmc/integrators/Integrator_algorithm.h ./qcd/hmc/integrators/Integrator_base.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 ./simd/Grid_avx.h ./simd/Grid_avx512.h ./simd/Grid_empty.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
CCFILES=./qcd/hmc/integrators/Integrator.cc ./qcd/hmc/HMC.cc ./qcd/utils/SpaceTimeGrid.cc ./qcd/action/fermion/WilsonKernels.cc ./qcd/action/fermion/PartialFractionFermion5D.cc ./qcd/action/fermion/CayleyFermion5D.cc ./qcd/action/fermion/WilsonKernelsHand.cc ./qcd/action/fermion/WilsonFermion.cc ./qcd/action/fermion/ContinuedFractionFermion5D.cc ./qcd/action/fermion/WilsonFermion5D.cc ./qcd/spin/Dirac.cc ./GridInit.cc ./algorithms/approx/MultiShiftFunction.cc ./algorithms/approx/Remez.cc ./algorithms/approx/Zolotarev.cc ./stencil/Lebesgue.cc ./stencil/Stencil_common.cc
CCFILES=./algorithms/approx/MultiShiftFunction.cc ./algorithms/approx/Remez.cc ./algorithms/approx/Zolotarev.cc ./GridInit.cc ./GridLog.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/WilsonKernelsHand.cc ./qcd/hmc/HMC.cc ./qcd/hmc/integrators/Integrator.cc ./qcd/spin/Dirac.cc ./qcd/utils/SpaceTimeGrid.cc ./stencil/Lebesgue.cc ./stencil/Stencil_common.cc

51
lib/algorithms/CoarsenedMatrix.h
View File

@@ -32,10 +32,10 @@ namespace Grid {
displacements[2*_d]=0;
//// report back
std::cout<<"directions :";
std::cout<<GridLogMessage<<"directions :";
for(int d=0;d<npoint;d++) std::cout<< directions[d]<< " ";
std::cout <<std::endl;
std::cout<<"displacements :";
std::cout<<GridLogMessage<<"displacements :";
for(int d=0;d<npoint;d++) std::cout<< displacements[d]<< " ";
std::cout<<std::endl;
}
@@ -100,9 +100,9 @@ namespace Grid {
eProj._odata[ss](i)=CComplex(1.0);
}
eProj=eProj - iProj;
std::cout<<"Orthog check error "<<i<<" " << norm2(eProj)<<std::endl;
std::cout<<GridLogMessage<<"Orthog check error "<<i<<" " << norm2(eProj)<<std::endl;
}
std::cout <<"CheckOrthog done"<<std::endl;
std::cout<<GridLogMessage <<"CheckOrthog done"<<std::endl;
}
void ProjectToSubspace(CoarseVector &CoarseVec,const FineField &FineVec){
blockProject(CoarseVec,FineVec,subspace);
@@ -113,7 +113,7 @@ namespace Grid {
void CreateSubspaceRandom(GridParallelRNG &RNG){
for(int i=0;i<nbasis;i++){
random(RNG,subspace[i]);
std::cout<<" norm subspace["<<i<<"] "<<norm2(subspace[i])<<std::endl;
std::cout<<GridLogMessage<<" norm subspace["<<i<<"] "<<norm2(subspace[i])<<std::endl;
}
Orthogonalise();
}
@@ -121,7 +121,7 @@ namespace Grid {
RealD scale;
ConjugateGradient<FineField> CG(2.0e-3,10000);
ConjugateGradient<FineField> CG(1.0e-2,10000);
FineField noise(FineGrid);
FineField Mn(FineGrid);
@@ -131,7 +131,7 @@ namespace Grid {
scale = std::pow(norm2(noise),-0.5);
noise=noise*scale;
hermop.Op(noise,Mn); std::cout << "noise ["<<b<<"] <n|MdagM|n> "<<norm2(Mn)<<std::endl;
hermop.Op(noise,Mn); std::cout<<GridLogMessage << "noise ["<<b<<"] <n|MdagM|n> "<<norm2(Mn)<<std::endl;
for(int i=0;i<1;i++){
@@ -143,7 +143,7 @@ namespace Grid {
}
hermop.Op(noise,Mn); std::cout << "filtered["<<b<<"] <f|MdagM|f> "<<norm2(Mn)<<std::endl;
hermop.Op(noise,Mn); std::cout<<GridLogMessage << "filtered["<<b<<"] <f|MdagM|f> "<<norm2(Mn)<<std::endl;
subspace[b] = noise;
}
@@ -189,7 +189,7 @@ namespace Grid {
SimpleCompressor<siteVector> compressor;
Stencil.HaloExchange(in,comm_buf,compressor);
//PARALLEL_FOR_LOOP
PARALLEL_FOR_LOOP
for(int ss=0;ss<Grid()->oSites();ss++){
siteVector res = zero;
siteVector nbr;
@@ -252,10 +252,6 @@ namespace Grid {
// Orthogonalise the subblocks over the basis
blockOrthogonalise(InnerProd,Subspace.subspace);
//Subspace.Orthogonalise();
// Subspace.CheckOrthogonal();
//Subspace.Orthogonalise();
// Subspace.CheckOrthogonal();
// Compute the matrix elements of linop between this orthonormal
// set of vectors.
@@ -306,6 +302,7 @@ namespace Grid {
Subspace.ProjectToSubspace(oProj,oblock);
// blockProject(iProj,iblock,Subspace.subspace);
// blockProject(oProj,oblock,Subspace.subspace);
PARALLEL_FOR_LOOP
for(int ss=0;ss<Grid()->oSites();ss++){
for(int j=0;j<nbasis;j++){
if( disp!= 0 ) {
@@ -321,12 +318,12 @@ namespace Grid {
///////////////////////////
// test code worth preserving in if block
///////////////////////////
std::cout<< " Computed matrix elements "<< self_stencil <<std::endl;
std::cout<<GridLogMessage<< " Computed matrix elements "<< self_stencil <<std::endl;
for(int p=0;p<geom.npoint;p++){
std::cout<< "A["<<p<<"]" << std::endl;
std::cout<< A[p] << std::endl;
std::cout<<GridLogMessage<< "A["<<p<<"]" << std::endl;
std::cout<<GridLogMessage<< A[p] << std::endl;
}
std::cout<< " picking by block0 "<< self_stencil <<std::endl;
std::cout<<GridLogMessage<< " picking by block0 "<< self_stencil <<std::endl;
phi=Subspace.subspace[0];
std::vector<int> bc(FineGrid->_ndimension,0);
@@ -334,9 +331,9 @@ namespace Grid {
blockPick(Grid(),phi,tmp,bc); // Pick out a block
linop.Op(tmp,Mphi); // Apply big dop
blockProject(iProj,Mphi,Subspace.subspace); // project it and print it
std::cout<< " Computed matrix elements from block zero only "<<std::endl;
std::cout<< iProj <<std::endl;
std::cout<<"Computed Coarse Operator"<<std::endl;
std::cout<<GridLogMessage<< " Computed matrix elements from block zero only "<<std::endl;
std::cout<<GridLogMessage<< iProj <<std::endl;
std::cout<<GridLogMessage<<"Computed Coarse Operator"<<std::endl;
#endif
// ForceHermitian();
AssertHermitian();
@@ -345,9 +342,9 @@ namespace Grid {
void ForceDiagonal(void) {
std::cout<<"**************************************************"<<std::endl;
std::cout<<"**** Forcing coarse operator to be diagonal ****"<<std::endl;
std::cout<<"**************************************************"<<std::endl;
std::cout<<GridLogMessage<<"**************************************************"<<std::endl;
std::cout<<GridLogMessage<<"**** Forcing coarse operator to be diagonal ****"<<std::endl;
std::cout<<GridLogMessage<<"**************************************************"<<std::endl;
for(int p=0;p<8;p++){
A[p]=zero;
}
@@ -387,13 +384,13 @@ namespace Grid {
Diff = AA - adj(AAc);
std::cout<<"Norm diff dim "<<d<<" "<< norm2(Diff)<<std::endl;
std::cout<<"Norm dim "<<d<<" "<< norm2(AA)<<std::endl;
std::cout<<GridLogMessage<<"Norm diff dim "<<d<<" "<< norm2(Diff)<<std::endl;
std::cout<<GridLogMessage<<"Norm dim "<<d<<" "<< norm2(AA)<<std::endl;
}
Diff = A[8] - adj(A[8]);
std::cout<<"Norm diff local "<< norm2(Diff)<<std::endl;
std::cout<<"Norm local "<< norm2(A[8])<<std::endl;
std::cout<<GridLogMessage<<"Norm diff local "<< norm2(Diff)<<std::endl;
std::cout<<GridLogMessage<<"Norm local "<< norm2(A[8])<<std::endl;
}
};

2
lib/algorithms/approx/Chebyshev.h
View File

@@ -120,7 +120,7 @@ namespace Grid {
Field *Tn = &T1;
Field *Tnp = &T2;
std::cout << "Chebyshev ["<<lo<<","<<hi<<"]"<< " order "<<order <<std::endl;
std::cout<<GridLogMessage << "Chebyshev ["<<lo<<","<<hi<<"]"<< " order "<<order <<std::endl;
// Tn=T1 = (xscale M + mscale)in
double xscale = 2.0/(hi-lo);
double mscale = -(hi+lo)/(hi-lo);

1
lib/algorithms/approx/Remez.cc
View File

@@ -757,3 +757,4 @@ void AlgRemez::csv(std::ostream & os)
}
return;
}

1
lib/algorithms/approx/Remez.h
View File

@@ -28,6 +28,7 @@
remez.getIPFE(res,pole,&norm);
remez.csv(ostream &os);
*/
class AlgRemez
{
private:

6
lib/algorithms/iterative/AdefGeneric.h
View File

@@ -149,7 +149,7 @@ class TwoLevelFlexiblePcg : public LinearFunction<Field>
}
RealD rrn=sqrt(rn/ssq);
std::cout<<"TwoLevelfPcg: k= "<<k<<" residual = "<<rrn<<std::endl;
std::cout<<GridLogMessage<<"TwoLevelfPcg: k= "<<k<<" residual = "<<rrn<<std::endl;
// Stopping condition
if ( rn <= rsq ) {
@@ -161,8 +161,8 @@ class TwoLevelFlexiblePcg : public LinearFunction<Field>
RealD srcnorm = sqrt(norm2(src));
RealD tmpnorm = sqrt(norm2(tmp));
RealD true_residual = tmpnorm/srcnorm;
std::cout<<"TwoLevelfPcg: true residual is "<<true_residual<<std::endl;
std::cout<<"TwoLevelfPcg: target residual was"<<Tolerance<<std::endl;
std::cout<<GridLogMessage<<"TwoLevelfPcg: true residual is "<<true_residual<<std::endl;
std::cout<<GridLogMessage<<"TwoLevelfPcg: target residual was"<<Tolerance<<std::endl;
return k;
}
}

20
lib/algorithms/iterative/ConjugateGradient.h
View File

@@ -43,12 +43,12 @@ public:
ssq=norm2(src);
if ( verbose ) {
std::cout <<std::setprecision(4)<< "ConjugateGradient: guess "<<guess<<std::endl;
std::cout <<std::setprecision(4)<< "ConjugateGradient: src "<<ssq <<std::endl;
std::cout <<std::setprecision(4)<< "ConjugateGradient: mp "<<d <<std::endl;
std::cout <<std::setprecision(4)<< "ConjugateGradient: mmp "<<b <<std::endl;
std::cout <<std::setprecision(4)<< "ConjugateGradient: cp,r "<<cp <<std::endl;
std::cout <<std::setprecision(4)<< "ConjugateGradient: p "<<a <<std::endl;
std::cout<<GridLogMessage <<std::setprecision(4)<< "ConjugateGradient: guess "<<guess<<std::endl;
std::cout<<GridLogMessage <<std::setprecision(4)<< "ConjugateGradient: src "<<ssq <<std::endl;
std::cout<<GridLogMessage <<std::setprecision(4)<< "ConjugateGradient: mp "<<d <<std::endl;
std::cout<<GridLogMessage <<std::setprecision(4)<< "ConjugateGradient: mmp "<<b <<std::endl;
std::cout<<GridLogMessage <<std::setprecision(4)<< "ConjugateGradient: cp,r "<<cp <<std::endl;
std::cout<<GridLogMessage <<std::setprecision(4)<< "ConjugateGradient: p "<<a <<std::endl;
}
RealD rsq = Tolerance* Tolerance*ssq;
@@ -58,7 +58,7 @@ public:
return;
}
if(verbose) std::cout << std::setprecision(4)<< "ConjugateGradient: k=0 residual "<<cp<<" rsq"<<rsq<<std::endl;
if(verbose) std::cout<<GridLogMessage << std::setprecision(4)<< "ConjugateGradient: k=0 residual "<<cp<<" rsq"<<rsq<<std::endl;
int k;
for (k=1;k<=MaxIterations;k++){
@@ -80,7 +80,7 @@ public:
psi= a*p+psi;
p = p*b+r;
if (verbose) std::cout<<"ConjugateGradient: Iteration " <<k<<" residual "<<cp<< " target"<< rsq<<std::endl;
if (verbose) std::cout<<GridLogMessage<<"ConjugateGradient: Iteration " <<k<<" residual "<<cp<< " target"<< rsq<<std::endl;
// Stopping condition
if ( cp <= rsq ) {
@@ -94,14 +94,14 @@ public:
RealD resnorm = sqrt(norm2(p));
RealD true_residual = resnorm/srcnorm;
std::cout<<"ConjugateGradient: Converged on iteration " <<k
std::cout<<GridLogMessage<<"ConjugateGradient: Converged on iteration " <<k
<<" computed residual "<<sqrt(cp/ssq)
<<" true residual "<<true_residual
<<" target "<<Tolerance<<std::endl;
return;
}
}
std::cout<<"ConjugateGradient did NOT converge"<<std::endl;
std::cout<<GridLogMessage<<"ConjugateGradient did NOT converge"<<std::endl;
assert(0);
}
};

14
lib/algorithms/iterative/ConjugateGradientMultiShift.h
View File

@@ -91,7 +91,7 @@ void operator() (LinearOperatorBase<Field> &Linop, const Field &src, std::vector
cp = norm2(src);
for(int s=0;s<nshift;s++){
rsq[s] = cp * mresidual[s] * mresidual[s];
std::cout<<"ConjugateGradientMultiShift: shift "<<s
std::cout<<GridLogMessage<<"ConjugateGradientMultiShift: shift "<<s
<<" target resid "<<rsq[s]<<std::endl;
ps[s] = src;
}
@@ -109,7 +109,7 @@ void operator() (LinearOperatorBase<Field> &Linop, const Field &src, std::vector
// p and mmp is equal to d after this since
// the d computation is tricky
// qq = real(innerProduct(p,mmp));
// std::cout << "debug equal ? qq "<<qq<<" d "<< d<<std::endl;
// std::cout<<GridLogMessage << "debug equal ? qq "<<qq<<" d "<< d<<std::endl;
b = -cp /d;
@@ -214,7 +214,7 @@ void operator() (LinearOperatorBase<Field> &Linop, const Field &src, std::vector
if(css<rsq[s]){
if ( ! converged[s] )
std::cout<<"ConjugateGradientMultiShift k="<<k<<" Shift "<<s<<" has converged"<<std::endl;
std::cout<<GridLogMessage<<"ConjugateGradientMultiShift k="<<k<<" Shift "<<s<<" has converged"<<std::endl;
converged[s]=1;
} else {
all_converged=0;
@@ -225,8 +225,8 @@ void operator() (LinearOperatorBase<Field> &Linop, const Field &src, std::vector
if ( all_converged ){
std::cout<< "CGMultiShift: All shifts have converged iteration "<<k<<std::endl;
std::cout<< "CGMultiShift: Checking solutions"<<std::endl;
std::cout<<GridLogMessage<< "CGMultiShift: All shifts have converged iteration "<<k<<std::endl;
std::cout<<GridLogMessage<< "CGMultiShift: Checking solutions"<<std::endl;
// Check answers
for(int s=0; s < nshift; s++) {
@@ -235,13 +235,13 @@ void operator() (LinearOperatorBase<Field> &Linop, const Field &src, std::vector
axpy(r,-alpha[s],src,tmp);
RealD rn = norm2(r);
RealD cn = norm2(src);
std::cout<<"CGMultiShift: shift["<<s<<"] true residual "<<std::sqrt(rn/cn)<<std::endl;
std::cout<<GridLogMessage<<"CGMultiShift: shift["<<s<<"] true residual "<<std::sqrt(rn/cn)<<std::endl;
}
return;
}
}
// ugly hack
std::cout<<"CG multi shift did not converge"<<std::endl;
std::cout<<GridLogMessage<<"CG multi shift did not converge"<<std::endl;
assert(0);
}

8
lib/algorithms/iterative/ConjugateResidual.h
View File

@@ -41,7 +41,7 @@ namespace Grid {
ssq=norm2(src);
rsq=Tolerance*Tolerance*ssq;
if (verbose) std::cout<<"ConjugateResidual: iteration " <<0<<" residual "<<cp<< " target"<< rsq<<std::endl;
if (verbose) std::cout<<GridLogMessage<<"ConjugateResidual: iteration " <<0<<" residual "<<cp<< " target"<< rsq<<std::endl;
for(int k=1;k<MaxIterations;k++){
@@ -60,13 +60,13 @@ namespace Grid {
axpy(p,b,p,r);
pAAp=axpy_norm(Ap,b,Ap,Ar);
if(verbose) std::cout<<"ConjugateResidual: iteration " <<k<<" residual "<<cp<< " target"<< rsq<<std::endl;
if(verbose) std::cout<<GridLogMessage<<"ConjugateResidual: iteration " <<k<<" residual "<<cp<< " target"<< rsq<<std::endl;
if(cp<rsq) {
Linop.HermOp(psi,Ap);
axpy(r,-1.0,src,Ap);
RealD true_resid = norm2(r)/ssq;
std::cout<<"ConjugateResidual: Converged on iteration " <<k
std::cout<<GridLogMessage<<"ConjugateResidual: Converged on iteration " <<k
<< " computed residual "<<sqrt(cp/ssq)
<< " true residual "<<sqrt(true_resid)
<< " target " <<Tolerance <<std::endl;
@@ -75,7 +75,7 @@ namespace Grid {
}
std::cout<<"ConjugateResidual did NOT converge"<<std::endl;
std::cout<<GridLogMessage<<"ConjugateResidual did NOT converge"<<std::endl;
assert(0);
}
};

26
lib/algorithms/iterative/PrecGeneralisedConjugateResidual.h
View File

@@ -45,13 +45,13 @@ namespace Grid {
cp=GCRnStep(Linop,src,psi,rsq);
if ( verbose ) std::cout<<"VPGCR("<<mmax<<","<<nstep<<") "<< steps <<" steps cp = "<<cp<<std::endl;
if ( verbose ) std::cout<<GridLogMessage<<"VPGCR("<<mmax<<","<<nstep<<") "<< steps <<" steps cp = "<<cp<<std::endl;
if(cp<rsq) {
Linop.HermOp(psi,r);
axpy(r,-1.0,src,r);
RealD tr = norm2(r);
std::cout<<"PrecGeneralisedConjugateResidual: Converged on iteration " <<steps
std::cout<<GridLogMessage<<"PrecGeneralisedConjugateResidual: Converged on iteration " <<steps
<< " computed residual "<<sqrt(cp/ssq)
<< " true residual " <<sqrt(tr/ssq)
<< " target " <<Tolerance <<std::endl;
@@ -59,7 +59,7 @@ namespace Grid {
}
}
std::cout<<"Variable Preconditioned GCR did not converge"<<std::endl;
std::cout<<GridLogMessage<<"Variable Preconditioned GCR did not converge"<<std::endl;
assert(0);
}
RealD GCRnStep(LinearOperatorBase<Field> &Linop,const Field &src, Field &psi,RealD rsq){
@@ -96,21 +96,21 @@ namespace Grid {
/////////////////////
Preconditioner(r,z);
std::cout<< " Preconditioner in " << norm2(r)<<std::endl;
std::cout<< " Preconditioner out " << norm2(z)<<std::endl;
std::cout<<GridLogMessage<< " Preconditioner in " << norm2(r)<<std::endl;
std::cout<<GridLogMessage<< " Preconditioner out " << norm2(z)<<std::endl;
Linop.HermOp(z,tmp);
std::cout<< " Preconditioner Aout " << norm2(tmp)<<std::endl;
std::cout<<GridLogMessage<< " Preconditioner Aout " << norm2(tmp)<<std::endl;
ttmp=tmp;
tmp=tmp-r;
std::cout<< " Preconditioner resid " << std::sqrt(norm2(tmp)/norm2(r))<<std::endl;
std::cout<<GridLogMessage<< " Preconditioner resid " << std::sqrt(norm2(tmp)/norm2(r))<<std::endl;
/*
std::cout<<r<<std::endl;
std::cout<<z<<std::endl;
std::cout<<ttmp<<std::endl;
std::cout<<tmp<<std::endl;
std::cout<<GridLogMessage<<r<<std::endl;
std::cout<<GridLogMessage<<z<<std::endl;
std::cout<<GridLogMessage<<ttmp<<std::endl;
std::cout<<GridLogMessage<<tmp<<std::endl;
*/
Linop.HermOpAndNorm(z,Az,zAz,zAAz);
@@ -137,7 +137,7 @@ namespace Grid {
cp = axpy_norm(r,-a,q[peri_k],r);
std::cout<< " VPGCR_step resid" <<sqrt(cp/rsq)<<std::endl;
std::cout<<GridLogMessage<< " VPGCR_step resid" <<sqrt(cp/rsq)<<std::endl;
if((k==nstep-1)||(cp<rsq)){
return cp;
}
@@ -148,7 +148,7 @@ namespace Grid {
Linop.HermOp(z,tmp);
tmp=tmp-r;
std::cout<< " Preconditioner resid" <<sqrt(norm2(tmp)/norm2(r))<<std::endl;
std::cout<<GridLogMessage<< " Preconditioner resid" <<sqrt(norm2(tmp)/norm2(r))<<std::endl;
q[peri_kp]=Az;
p[peri_kp]=z;

4
lib/algorithms/iterative/SchurRedBlack.h
View File

@@ -89,7 +89,7 @@ namespace Grid {
//////////////////////////////////////////////////////////////
// Call the red-black solver
//////////////////////////////////////////////////////////////
std::cout << "SchurRedBlack solver calling the MpcDagMp solver" <<std::endl;
std::cout<<GridLogMessage << "SchurRedBlack solver calling the MpcDagMp solver" <<std::endl;
_HermitianRBSolver(_HermOpEO,src_o,sol_o); assert(sol_o.checkerboard==Odd);
///////////////////////////////////////////////////
@@ -108,7 +108,7 @@ namespace Grid {
RealD ns = norm2(in);
RealD nr = norm2(resid);
std::cout << "SchurRedBlackDiagMooee solver true unprec resid "<< std::sqrt(nr/ns) <<" nr "<< nr <<" ns "<<ns << std::endl;
std::cout<<GridLogMessage << "SchurRedBlackDiagMooee solver true unprec resid "<< std::sqrt(nr/ns) <<" nr "<< nr <<" ns "<<ns << std::endl;
}
};

12
lib/lattice/Lattice_ET.h
View File

@@ -132,18 +132,18 @@ inline void CBFromExpression(int &cb,const T1& lat) // Lattice leaf
assert(cb==lat.checkerboard);
}
cb=lat.checkerboard;
// std::cout<<"Lattice leaf cb "<<cb<<std::endl;
// std::cout<<GridLogMessage<<"Lattice leaf cb "<<cb<<std::endl;
}
template<class T1,typename std::enable_if<!is_lattice<T1>::value, T1>::type * = nullptr >
inline void CBFromExpression(int &cb,const T1& notlat) // non-lattice leaf
{
// std::cout<<"Non lattice leaf cb"<<cb<<std::endl;
// std::cout<<GridLogMessage<<"Non lattice leaf cb"<<cb<<std::endl;
}
template <typename Op, typename T1>
inline void CBFromExpression(int &cb,const LatticeUnaryExpression<Op,T1 > &expr)
{
CBFromExpression(cb,std::get<0>(expr.second));// recurse
// std::cout<<"Unary node cb "<<cb<<std::endl;
// std::cout<<GridLogMessage<<"Unary node cb "<<cb<<std::endl;
}
template <typename Op, typename T1, typename T2>
@@ -151,7 +151,7 @@ inline void CBFromExpression(int &cb,const LatticeBinaryExpression<Op,T1,T2> &ex
{
CBFromExpression(cb,std::get<0>(expr.second));// recurse
CBFromExpression(cb,std::get<1>(expr.second));
// std::cout<<"Binary node cb "<<cb<<std::endl;
// std::cout<<GridLogMessage<<"Binary node cb "<<cb<<std::endl;
}
template <typename Op, typename T1, typename T2, typename T3>
inline void CBFromExpression( int &cb,const LatticeTrinaryExpression<Op,T1,T2,T3 > &expr)
@@ -159,7 +159,7 @@ inline void CBFromExpression( int &cb,const LatticeTrinaryExpression<Op,T1,T2,T3
CBFromExpression(cb,std::get<0>(expr.second));// recurse
CBFromExpression(cb,std::get<1>(expr.second));
CBFromExpression(cb,std::get<2>(expr.second));
// std::cout<<"Trinary node cb "<<cb<<std::endl;
// std::cout<<GridLogMessage<<"Trinary node cb "<<cb<<std::endl;
}
////////////////////////////////////////////
@@ -370,7 +370,7 @@ using namespace Grid;
tmp.func(eval(0,v1),eval(0,v2));
auto var = v1+v2;
std::cout<<typeid(var).name()<<std::endl;
std::cout<<GridLogMessage<<typeid(var).name()<<std::endl;
v3=v1+v2;
v3=v1+v2+v1*v2;

2
lib/lattice/Lattice_base.h
View File

@@ -221,7 +221,7 @@ PARALLEL_FOR_LOOP
template<class robj> strong_inline Lattice<vobj> & operator = (const Lattice<robj> & r){
this->checkerboard = r.checkerboard;
conformable(*this,r);
std::cout<<"Lattice operator ="<<std::endl;
std::cout<<GridLogMessage<<"Lattice operator ="<<std::endl;
PARALLEL_FOR_LOOP
for(int ss=0;ss<_grid->oSites();ss++){
this->_odata[ss]=r._odata[ss];

2
lib/lattice/Lattice_reduction.h
View File

@@ -125,7 +125,7 @@ template<class vobj> inline void sliceSum(const Lattice<vobj> &Data,std::vector<
assert(grid!=NULL);
// FIXME
std::cout<<"WARNING ! SliceSum is unthreaded "<<grid->SumArraySize()<<" threads "<<std::endl;
std::cout<<GridLogMessage<<"WARNING ! SliceSum is unthreaded "<<grid->SumArraySize()<<" threads "<<std::endl;
const int Nd = grid->_ndimension;
const int Nsimd = grid->Nsimd();

1
lib/lattice/Lattice_rng.h
View File

@@ -234,6 +234,7 @@ namespace Grid {
std::vector<scalar_object> buf(Nsimd);
PARALLEL_FOR_LOOP
for(int ss=0;ss<osites;ss++){
for(int si=0;si<Nsimd;si++){

5
lib/lattice/Lattice_transfer.h
View File

@@ -23,7 +23,7 @@ inline void subdivides(GridBase *coarse,GridBase *fine)
template<class vobj> inline void pickCheckerboard(int cb,Lattice<vobj> &half,const Lattice<vobj> &full){
half.checkerboard = cb;
int ssh=0;
PARALLEL_FOR_LOOP
//PARALLEL_FOR_LOOP
for(int ss=0;ss<full._grid->oSites();ss++){
std::vector<int> coor;
int cbos;
@@ -40,7 +40,7 @@ PARALLEL_FOR_LOOP
template<class vobj> inline void setCheckerboard(Lattice<vobj> &full,const Lattice<vobj> &half){
int cb = half.checkerboard;
int ssh=0;
PARALLEL_FOR_LOOP
//PARALLEL_FOR_LOOP
for(int ss=0;ss<full._grid->oSites();ss++){
std::vector<int> coor;
int cbos;
@@ -158,6 +158,7 @@ template<class vobj,class CComplex>
fine_inner = localInnerProduct(fineX,fineY);
blockSum(coarse_inner,fine_inner);
PARALLEL_FOR_LOOP
for(int ss=0;ss<coarse->oSites();ss++){
CoarseInner._odata[ss] = coarse_inner._odata[ss];
}

14
lib/qcd/action/fermion/ContinuedFractionFermion5D.cc
View File

@@ -10,12 +10,12 @@ namespace Grid {
void ContinuedFractionFermion5D::SetCoefficientsZolotarev(RealD zolo_hi,Approx::zolotarev_data *zdata)
{
// How to check Ls matches??
// std::cout << Ls << " Ls"<<std::endl;
// std::cout << zdata->n << " - n"<<std::endl;
// std::cout << zdata->da << " -da "<<std::endl;
// std::cout << zdata->db << " -db"<<std::endl;
// std::cout << zdata->dn << " -dn"<<std::endl;
// std::cout << zdata->dd << " -dd"<<std::endl;
// std::cout<<GridLogMessage << Ls << " Ls"<<std::endl;
// std::cout<<GridLogMessage << zdata->n << " - n"<<std::endl;
// std::cout<<GridLogMessage << zdata->da << " -da "<<std::endl;
// std::cout<<GridLogMessage << zdata->db << " -db"<<std::endl;
// std::cout<<GridLogMessage << zdata->dn << " -dn"<<std::endl;
// std::cout<<GridLogMessage << zdata->dd << " -dd"<<std::endl;
assert(zdata->db==Ls);// Beta has Ls coeffs
@@ -55,7 +55,7 @@ namespace Grid {
See[s] = Aee[s] - 1.0/See[s-1];
}
for(int s=0;s<Ls;s++){
std::cout <<"s = "<<s<<" Beta "<<Beta[s]<<" Aee "<<Aee[s] <<" See "<<See[s] <<std::endl;
std::cout<<GridLogMessage <<"s = "<<s<<" Beta "<<Beta[s]<<" Aee "<<Aee[s] <<" See "<<See[s] <<std::endl;
}
}

2
lib/qcd/action/fermion/DomainWallFermion.h
View File

@@ -32,7 +32,7 @@ namespace Grid {
Approx::zolotarev_data *zdata = Approx::higham(eps,this->Ls);// eps is ignored for higham
assert(zdata->n==this->Ls);
std::cout << "DomainWallFermion with Ls="<<Ls<<std::endl;
std::cout<<GridLogMessage << "DomainWallFermion with Ls="<<Ls<<std::endl;
// Call base setter
this->CayleyFermion5D::SetCoefficientsTanh(zdata,1.0,0.0);

2
lib/qcd/action/fermion/MobiusFermion.h
View File

@@ -30,7 +30,7 @@ namespace Grid {
{
RealD eps = 1.0;
std::cout << "MobiusFermion (b="<<b<<",c="<<c<<") with Ls= "<<Ls<<" Tanh approx"<<std::endl;
std::cout<<GridLogMessage << "MobiusFermion (b="<<b<<",c="<<c<<") with Ls= "<<Ls<<" Tanh approx"<<std::endl;
Approx::zolotarev_data *zdata = Approx::higham(eps,this->Ls);// eps is ignored for higham
assert(zdata->n==this->Ls);

2
lib/qcd/action/fermion/MobiusZolotarevFermion.h
View File

@@ -34,7 +34,7 @@ namespace Grid {
Approx::zolotarev_data *zdata = Approx::zolotarev(eps,this->Ls,0);
assert(zdata->n==this->Ls);
std::cout << "MobiusZolotarevFermion (b="<<b<<",c="<<c<<") with Ls= "<<Ls<<" Zolotarev range ["<<lo<<","<<hi<<"]"<<std::endl;
std::cout<<GridLogMessage << "MobiusZolotarevFermion (b="<<b<<",c="<<c<<") with Ls= "<<Ls<<" Zolotarev range ["<<lo<<","<<hi<<"]"<<std::endl;
// Call base setter
this->CayleyFermion5D::SetCoefficientsZolotarev(hi,zdata,b,c);

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

@@ -260,12 +260,12 @@ namespace Grid {
void PartialFractionFermion5D::SetCoefficientsZolotarev(RealD zolo_hi,Approx::zolotarev_data *zdata){
// check on degree matching
// std::cout << Ls << " Ls"<<std::endl;
// std::cout << zdata->n << " - n"<<std::endl;
// std::cout << zdata->da << " -da "<<std::endl;
// std::cout << zdata->db << " -db"<<std::endl;
// std::cout << zdata->dn << " -dn"<<std::endl;
// std::cout << zdata->dd << " -dd"<<std::endl;
// std::cout<<GridLogMessage << Ls << " Ls"<<std::endl;
// std::cout<<GridLogMessage << zdata->n << " - n"<<std::endl;
// std::cout<<GridLogMessage << zdata->da << " -da "<<std::endl;
// std::cout<<GridLogMessage << zdata->db << " -db"<<std::endl;
// std::cout<<GridLogMessage << zdata->dn << " -dn"<<std::endl;
// std::cout<<GridLogMessage << zdata->dd << " -dd"<<std::endl;
assert(Ls == (2*zdata->da -1) );
// Part frac

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

@@ -100,7 +100,7 @@ void WilsonFermion5D::DhopDir(const LatticeFermion &in, LatticeFermion &out,int
assert(dirdisp<=7);
assert(dirdisp>=0);
//PARALLEL_FOR_LOOP
PARALLEL_FOR_LOOP
for(int ss=0;ss<Umu._grid->oSites();ss++){
for(int s=0;s<Ls;s++){
int sU=ss;

20
lib/qcd/action/fermion/WilsonKernelsHand.cc
View File

@@ -360,11 +360,11 @@ void DiracOptHandDhopSite(CartesianStencil &st,LatticeDoubledGaugeField &U,
MULT_2SPIN(Xp);
}
XP_RECON;
// std::cout << "XP_RECON"<<std::endl;
// std::cout << result_00 <<" "<<result_01 <<" "<<result_02 <<std::endl;
// std::cout << result_10 <<" "<<result_11 <<" "<<result_12 <<std::endl;
// std::cout << result_20 <<" "<<result_21 <<" "<<result_22 <<std::endl;
// std::cout << result_30 <<" "<<result_31 <<" "<<result_32 <<std::endl;
// std::cout<<GridLogMessage << "XP_RECON"<<std::endl;
// std::cout<<GridLogMessage << result_00 <<" "<<result_01 <<" "<<result_02 <<std::endl;
// std::cout<<GridLogMessage << result_10 <<" "<<result_11 <<" "<<result_12 <<std::endl;
// std::cout<<GridLogMessage << result_20 <<" "<<result_21 <<" "<<result_22 <<std::endl;
// std::cout<<GridLogMessage << result_30 <<" "<<result_31 <<" "<<result_32 <<std::endl;
// Yp
offset = st._offsets [Yp][ss];
@@ -446,11 +446,11 @@ void DiracOptHandDhopSite(CartesianStencil &st,LatticeDoubledGaugeField &U,
MULT_2SPIN(Xm);
}
XM_RECON_ACCUM;
// std::cout << "XM_RECON_ACCUM"<<std::endl;
// std::cout << result_00 <<" "<<result_01 <<" "<<result_02 <<std::endl;
// std::cout << result_10 <<" "<<result_11 <<" "<<result_12 <<std::endl;
// std::cout << result_20 <<" "<<result_21 <<" "<<result_22 <<std::endl;
// std::cout << result_30 <<" "<<result_31 <<" "<<result_32 <<std::endl;
// std::cout<<GridLogMessage << "XM_RECON_ACCUM"<<std::endl;
// std::cout<<GridLogMessage << result_00 <<" "<<result_01 <<" "<<result_02 <<std::endl;
// std::cout<<GridLogMessage << result_10 <<" "<<result_11 <<" "<<result_12 <<std::endl;
// std::cout<<GridLogMessage << result_20 <<" "<<result_21 <<" "<<result_22 <<std::endl;
// std::cout<<GridLogMessage << result_30 <<" "<<result_31 <<" "<<result_32 <<std::endl;
// Ym

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

@@ -18,7 +18,7 @@ namespace Grid{
virtual RealD S(const GaugeField &U) {
RealD plaq = WilsonLoops<MatrixField,GaugeField>::avgPlaquette(U);
std::cout << "Plaq : "<<plaq << "\n";
std::cout<<GridLogMessage << "Plaq : "<<plaq << "\n";
double vol = U._grid->gSites();
return beta*(1.0 -plaq)*(Nd*(Nd-1.0))*vol*0.5;
};

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

@@ -36,15 +36,15 @@ namespace Grid{
RealD prob = std::exp(-DeltaH);
random(sRNG,rn_test);
std::cout<< "--------------------------------------------\n";
std::cout<< "dH = "<<DeltaH << " Random = "<< rn_test
std::cout<<GridLogMessage<< "--------------------------------------------\n";
std::cout<<GridLogMessage<< "dH = "<<DeltaH << " Random = "<< rn_test
<< "\nAcc. Probability = " << ((prob<1.0)? prob: 1.0)<< " ";
if((prob >1.0) || (rn_test <= prob)){ // accepted
std::cout <<"-- ACCEPTED\n";
std::cout<<GridLogMessage <<"-- ACCEPTED\n";
return true;
} else { // rejected
std::cout <<"-- REJECTED\n";
std::cout<<GridLogMessage <<"-- REJECTED\n";
return false;
}
}
@@ -53,12 +53,12 @@ namespace Grid{
MD.init(U); // set U and initialize P and phi's
RealD H0 = MD.S(U); // initial state action
std::cout<<"Total H before = "<< H0 << "\n";
std::cout<<GridLogMessage<<"Total H before = "<< H0 << "\n";
MD.integrate(U);
RealD H1 = MD.S(U); // updated state action
std::cout<<"Total H after = "<< H1 << "\n";
std::cout<<GridLogMessage<<"Total H after = "<< H1 << "\n";
return (H1-H0);
}
@@ -81,17 +81,17 @@ namespace Grid{
// Thermalizations
for(int iter=1; iter <= Params.ThermalizationSteps; ++iter){
std::cout << "-- # Thermalization step = "<< iter << "\n";
std::cout<<GridLogMessage << "-- # Thermalization step = "<< iter << "\n";
DeltaH = evolve_step(Uin);
std::cout<< " dH = "<< DeltaH << "\n";
std::cout<<GridLogMessage<< " dH = "<< DeltaH << "\n";
}
// Actual updates (evolve a copy Ucopy then copy back eventually)
LatticeLorentzColourMatrix Ucopy(Uin._grid);
for(int iter=Params.StartingConfig;
iter < Params.Nsweeps+Params.StartingConfig; ++iter){
std::cout << "-- # Sweep = "<< iter << "\n";
std::cout<<GridLogMessage << "-- # Sweep = "<< iter << "\n";
Ucopy = Uin;
DeltaH = evolve_step(Ucopy);

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

@@ -95,7 +95,7 @@ namespace Grid{
//Initialization of momenta and actions
void init(LatticeLorentzColourMatrix& U){
std::cout<< "Integrator init\n";
std::cout<<GridLogMessage<< "Integrator init\n";
MDutils::generate_momenta(*P,pRNG);
for(int level=0; level< as.size(); ++level){
@@ -119,7 +119,7 @@ namespace Grid{
RealD H = Hsum.real();
std::cout << "H_p = "<< H << "\n";
std::cout<<GridLogMessage << "H_p = "<< H << "\n";
// Actions
for(int level=0; level<as.size(); ++level)

40
lib/qcd/hmc/integrators/Integrator_algorithm.h
View File

@@ -39,29 +39,29 @@ namespace Grid{
if(clock[level] == 0){ // initial half step
Integ->update_P(U,level,lambda*eps);
++clock[level];
for(int l=0; l<level;++l) std::cout<<" ";
std::cout<<"P "<< clock[level] <<std::endl;
for(int l=0; l<level;++l) std::cout<<GridLogMessage<<" ";
std::cout<<GridLogMessage<<"P "<< clock[level] <<std::endl;
}
if(level == fl){ // lowest level
Integ->update_U(U,0.5*eps);
for(int l=0; l<level;++l) std::cout<<" ";
std::cout<<"U "<< (clock[level]+1) <<std::endl;
for(int l=0; l<level;++l) std::cout<<GridLogMessage<<" ";
std::cout<<GridLogMessage<<"U "<< (clock[level]+1) <<std::endl;
}else{ // recursive function call
step(U,level+1,clock, Integ);
}
Integ->update_P(U,level,(1.0-2.0*lambda)*eps);
++clock[level];
for(int l=0; l<level;++l) std::cout<<" ";
std::cout<<"P "<< (clock[level]) <<std::endl;
for(int l=0; l<level;++l) std::cout<<GridLogMessage<<" ";
std::cout<<GridLogMessage<<"P "<< (clock[level]) <<std::endl;
if(level == fl){ // lowest level
Integ->update_U(U,0.5*eps);
for(int l=0; l<level;++l) std::cout<<" ";
std::cout<<"U "<< (clock[level]+1) <<std::endl;
for(int l=0; l<level;++l) std::cout<<GridLogMessage<<" ";
std::cout<<GridLogMessage<<"U "<< (clock[level]+1) <<std::endl;
}else{ // recursive function call
step(U,level+1,clock, Integ);
}
@@ -71,14 +71,14 @@ namespace Grid{
Integ->update_P(U,level,lambda*eps);
++clock[level];
for(int l=0; l<level;++l) std::cout<<" ";
std::cout<<"P "<< clock[level] <<std::endl;
for(int l=0; l<level;++l) std::cout<<GridLogMessage<<" ";
std::cout<<GridLogMessage<<"P "<< clock[level] <<std::endl;
}else{ // bulk step
Integ->update_P(U,level,lambda*2.0*eps);
clock[level]+=2;
for(int l=0; l<level;++l) std::cout<<" ";
std::cout<<"P "<< clock[level] <<std::endl;
for(int l=0; l<level;++l) std::cout<<GridLogMessage<<" ";
std::cout<<GridLogMessage<<"P "<< clock[level] <<std::endl;
}
}
@@ -112,13 +112,13 @@ namespace Grid{
if(clock[level] == 0){ // initial half step
Integ->update_P(U, level,eps/2.0);
++clock[level];
for(int l=0; l<level;++l) std::cout<<" ";
std::cout<<"P "<< 0.5*clock[level] <<std::endl;
for(int l=0; l<level;++l) std::cout<<GridLogMessage<<" ";
std::cout<<GridLogMessage<<"P "<< 0.5*clock[level] <<std::endl;
}
if(level == fl){ // lowest level
Integ->update_U(U, eps);
for(int l=0; l<level;++l) std::cout<<" ";
std::cout<<"U "<< 0.5*(clock[level]+1) <<std::endl;
for(int l=0; l<level;++l) std::cout<<GridLogMessage<<" ";
std::cout<<GridLogMessage<<"U "<< 0.5*(clock[level]+1) <<std::endl;
}else{ // recursive function call
step(U, level+1,clock, Integ);
}
@@ -126,14 +126,14 @@ namespace Grid{
Integ->update_P(U, level,eps/2.0);
++clock[level];
for(int l=0; l<level;++l) std::cout<<" ";
std::cout<<"P "<< 0.5*clock[level] <<std::endl;
for(int l=0; l<level;++l) std::cout<<GridLogMessage<<" ";
std::cout<<GridLogMessage<<"P "<< 0.5*clock[level] <<std::endl;
}else{ // bulk step
Integ->update_P(U, level,eps);
clock[level]+=2;
for(int l=0; l<level;++l) std::cout<<" ";
std::cout<<"P "<< 0.5*clock[level] <<std::endl;
for(int l=0; l<level;++l) std::cout<<GridLogMessage<<" ";
std::cout<<GridLogMessage<<"P "<< 0.5*clock[level] <<std::endl;
}
}

18
lib/qcd/hmc/integrators/Integrator_base.h
View File

@@ -101,7 +101,7 @@ namespace Grid{
//Initialization of momenta and actions
void init(LatticeLorentzColourMatrix& U,
GridParallelRNG& pRNG){
std::cout<< "Integrator init\n";
std::cout<<GridLogMessage<< "Integrator init\n";
if (!P)
P = new LatticeLorentzColourMatrix(U._grid);
MDutils::generate_momenta(*P,pRNG);
@@ -172,13 +172,13 @@ namespace Grid{
if(clock[level] == 0){ // initial half step
Integ->update_P(U, level,eps/2);
++clock[level];
for(int l=0; l<level;++l) std::cout<<" ";
std::cout<<"P "<< 0.5*clock[level] <<std::endl;
for(int l=0; l<level;++l) std::cout<<GridLogMessage<<" ";
std::cout<<GridLogMessage<<"P "<< 0.5*clock[level] <<std::endl;
}
if(level == fl){ // lowest level
Integ->update_U(U, eps);
for(int l=0; l<level;++l) std::cout<<" ";
std::cout<<"U "<< 0.5*(clock[level]+1) <<std::endl;
for(int l=0; l<level;++l) std::cout<<GridLogMessage<<" ";
std::cout<<GridLogMessage<<"U "<< 0.5*(clock[level]+1) <<std::endl;
}else{ // recursive function call
step(U, level+1,clock, Integ);
}
@@ -186,14 +186,14 @@ namespace Grid{
Integ->update_P(U, level,eps/2);
++clock[level];
for(int l=0; l<level;++l) std::cout<<" ";
std::cout<<"P "<< 0.5*clock[level] <<std::endl;
for(int l=0; l<level;++l) std::cout<<GridLogMessage<<" ";
std::cout<<GridLogMessage<<"P "<< 0.5*clock[level] <<std::endl;
}else{ // bulk step
Integ->update_P(U, level,eps);
clock[level]+=2;
for(int l=0; l<level;++l) std::cout<<" ";
std::cout<<"P "<< 0.5*clock[level] <<std::endl;
for(int l=0; l<level;++l) std::cout<<GridLogMessage<<" ";
std::cout<<GridLogMessage<<"P "<< 0.5*clock[level] <<std::endl;
}
}

28
lib/qcd/utils/SUn.h
View File

@@ -372,7 +372,7 @@ Note that in step D setting B ~ X - A and using B in place of A in step E will g
LatticeReal d(grid); d=zero;
LatticeReal alpha(grid);
// std::cout<<"xi "<<xi <<std::endl;
// std::cout<<GridLogMessage<<"xi "<<xi <<std::endl;
alpha = toReal(2.0*xi);
do {
@@ -468,11 +468,11 @@ Note that in step D setting B ~ X - A and using B in place of A in step E will g
LatticeMatrix Vcheck(grid);
Vcheck = zero;
Vcheck = where(Accepted,V*adj(V) - 1.0,Vcheck);
// std::cout << "SU3 check " <<norm2(Vcheck)<<std::endl;
// std::cout<<GridLogMessage << "SU3 check " <<norm2(Vcheck)<<std::endl;
assert(norm2(Vcheck)<1.0e-4);
// Verify the link stays in SU(3)
// std::cout <<"Checking the modified link"<<std::endl;
// std::cout<<GridLogMessage <<"Checking the modified link"<<std::endl;
Vcheck = link*adj(link) - 1.0;
assert(norm2(Vcheck)<1.0e-4);
/////////////////////////////////
@@ -483,42 +483,42 @@ Note that in step D setting B ~ X - A and using B in place of A in step E will g
for(int gen=0;gen<generators();gen++){
Matrix ta;
generator(gen,ta);
std::cout<< "Nc = "<<ncolour<<" t_"<<gen<<std::endl;
std::cout<<ta<<std::endl;
std::cout<<GridLogMessage<< "Nc = "<<ncolour<<" t_"<<gen<<std::endl;
std::cout<<GridLogMessage<<ta<<std::endl;
}
}
static void testGenerators(void){
Matrix ta;
Matrix tb;
std::cout<<"Checking trace ta tb is 0.5 delta_ab"<<std::endl;
std::cout<<GridLogMessage<<"Checking trace ta tb is 0.5 delta_ab"<<std::endl;
for(int a=0;a<generators();a++){
for(int b=0;b<generators();b++){
generator(a,ta);
generator(b,tb);
Complex tr =TensorRemove(trace(ta*tb));
std::cout<<tr<<" ";
std::cout<<GridLogMessage<<tr<<" ";
if(a==b) assert(abs(tr-Complex(0.5))<1.0e-6);
if(a!=b) assert(abs(tr)<1.0e-6);
}
std::cout<<std::endl;
std::cout<<GridLogMessage<<std::endl;
}
std::cout<<"Checking hermitian"<<std::endl;
std::cout<<GridLogMessage<<"Checking hermitian"<<std::endl;
for(int a=0;a<generators();a++){
generator(a,ta);
std::cout<<a<<" ";
std::cout<<GridLogMessage<<a<<" ";
assert(norm2(ta-adj(ta))<1.0e-6);
}
std::cout<<std::endl;
std::cout<<GridLogMessage<<std::endl;
std::cout<<"Checking traceless"<<std::endl;
std::cout<<GridLogMessage<<"Checking traceless"<<std::endl;
for(int a=0;a<generators();a++){
generator(a,ta);
Complex tr =TensorRemove(trace(ta));
std::cout<<a<<" ";
std::cout<<GridLogMessage<<a<<" ";
assert(abs(tr)<1.0e-6);
}
std::cout<<std::endl;
std::cout<<GridLogMessage<<std::endl;
}
// reunitarise??

6
lib/tensors/Tensor_extract_merge.h
View File

@@ -57,11 +57,11 @@ inline void extract(typename std::enable_if<!isGridTensor<vsimd>::value, const v
extracted[i]=buf[i*s];
for(int ii=1;ii<s;ii++){
if ( buf[i*s]!=buf[i*s+ii] ){
std::cout << " SIMD extract failure splat = "<<s<<" ii "<<ii<<" " <<Nextr<<" "<< Nsimd<<" "<<std::endl;
std::cout<<GridLogMessage << " SIMD extract failure splat = "<<s<<" ii "<<ii<<" " <<Nextr<<" "<< Nsimd<<" "<<std::endl;
for(int vv=0;vv<Nsimd;vv++) {
std::cout<< buf[vv]<<" ";
std::cout<<GridLogMessage<< buf[vv]<<" ";
}
std::cout<<std::endl;
std::cout<<GridLogMessage<<std::endl;
assert(0);
}
assert(buf[i*s]==buf[i*s+ii]);

10
tests/Make.inc
View File

@@ -1,5 +1,9 @@
bin_PROGRAMS = Test_cayley_cg Test_cayley_coarsen_support Test_cayley_even_odd Test_cayley_ldop_cr Test_cf_coarsen_support Test_cf_cr_unprec Test_contfrac_cg Test_contfrac_even_odd 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_fpgcr Test_dwf_hdcr Test_gamma Test_GaugeAction Test_hmc_WilsonGauge Test_lie_generators Test_main Test_multishift_sqrt Test_nersc_io Test_quenched_update Test_remez Test_rng Test_rng_fixed Test_simd Test_stencil Test_wilson_cg_prec Test_wilson_cg_schur Test_wilson_cg_unprec Test_wilson_cr_unprec Test_wilson_even_odd
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_contfrac_cg Test_contfrac_even_odd 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_fpgcr Test_dwf_hdcr Test_gamma Test_hmc_WilsonGauge Test_lie_generators Test_main Test_multishift_sqrt Test_nersc_io Test_quenched_update Test_remez Test_rng Test_rng_fixed Test_simd Test_stencil Test_wilson_cg_prec Test_wilson_cg_schur Test_wilson_cg_unprec Test_wilson_cr_unprec Test_wilson_even_odd
Test_GaugeAction_SOURCES=Test_GaugeAction.cc
Test_GaugeAction_LDADD=-lGrid
Test_cayley_cg_SOURCES=Test_cayley_cg.cc
@@ -74,10 +78,6 @@ Test_gamma_SOURCES=Test_gamma.cc
Test_gamma_LDADD=-lGrid
Test_GaugeAction_SOURCES=Test_GaugeAction.cc
Test_GaugeAction_LDADD=-lGrid
Test_hmc_WilsonGauge_SOURCES=Test_hmc_WilsonGauge.cc
Test_hmc_WilsonGauge_LDADD=-lGrid

16
tests/Test_GaugeAction.cc
View File

@@ -94,7 +94,7 @@ int main (int argc, char ** argv)
double vol = Fine.gSites();
Complex PlaqScale(1.0/vol/6.0/3.0);
RealD StapScale(1.0/vol/6.0/3.0);
std::cout <<"PlaqScale" << PlaqScale<<std::endl;
std::cout<<GridLogMessage <<"PlaqScale" << PlaqScale<<std::endl;
std::vector<TComplex> Plaq_T(orthosz);
sliceSum(Plaq,Plaq_T,Nd-1);
int Nt = Plaq_T.size();
@@ -105,20 +105,20 @@ int main (int argc, char ** argv)
for(int t=0;t<Nt;t++){
Plaq_T_sum = Plaq_T_sum+Plaq_T[t];
Complex Pt=TensorRemove(Plaq_T[t]);
std::cout << "sliced ["<<t<<"]" <<Pt*PlaqScale*Real(Nt) << std::endl;
std::cout<<GridLogMessage << "sliced ["<<t<<"]" <<Pt*PlaqScale*Real(Nt) << std::endl;
}
{
Complex Pt = TensorRemove(Plaq_T_sum);
std::cout << "total " <<Pt*PlaqScale<<std::endl;
std::cout<<GridLogMessage << "total " <<Pt*PlaqScale<<std::endl;
}
TComplex Tp = sum(Plaq);
Complex p = TensorRemove(Tp);
std::cout << "calculated plaquettes " <<p*PlaqScale<<std::endl;
std::cout<<GridLogMessage << "calculated plaquettes " <<p*PlaqScale<<std::endl;
RealD avg_plaq = ColourWilsonLoops::avgPlaquette(Umu);
std::cout << "NEW : calculated real plaquettes " <<avg_plaq<<std::endl;
std::cout<<GridLogMessage << "NEW : calculated real plaquettes " <<avg_plaq<<std::endl;
RealD stap_plaq=0.0;
LatticeColourMatrix stap(&Fine);
@@ -130,16 +130,16 @@ int main (int argc, char ** argv)
Complex s = TensorRemove(Ts);
stap_plaq+=real(s);
}
std::cout << "NEW : plaquette via staples"<< stap_plaq*StapScale*0.25<< std::endl;
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 << "calculated link trace " <<l*LinkTraceScale<<std::endl;
std::cout<<GridLogMessage << "calculated link trace " <<l*LinkTraceScale<<std::endl;
blockSum(cPlaq,Plaq);
TComplex TcP = sum(cPlaq);
Complex ll= TensorRemove(TcP);
std::cout << "coarsened plaquettes sum to " <<ll*PlaqScale<<std::endl;
std::cout<<GridLogMessage << "coarsened plaquettes sum to " <<ll*PlaqScale<<std::endl;
Grid_finalize();

22
tests/Test_cayley_cg.cc
View File

@@ -52,7 +52,7 @@ int main (int argc, char ** argv)
Grid_init(&argc,&argv);
int threads = GridThread::GetThreads();
std::cout << "Grid is setup to use "<<threads<<" threads"<<std::endl;
std::cout<<GridLogMessage << "Grid is setup to use "<<threads<<" threads"<<std::endl;
const int Ls=8;
GridCartesian * UGrid = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi());
@@ -71,33 +71,33 @@ int main (int argc, char ** argv)
RealD mass=0.1;
RealD M5 =1.8;
std::cout <<"DomainWallFermion test"<<std::endl;
std::cout<<GridLogMessage <<"DomainWallFermion test"<<std::endl;
DomainWallFermion Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5);
TestCGinversions<DomainWallFermion>(Ddwf,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5);
RealD b=1.5;// Scale factor b+c=2, b-c=1
RealD c=0.5;
std::cout <<"MobiusFermion test"<<std::endl;
std::cout<<GridLogMessage <<"MobiusFermion test"<<std::endl;
MobiusFermion Dmob(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,b,c);
TestCGinversions<MobiusFermion>(Dmob,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5);
std::cout <<"MobiusZolotarevFermion test"<<std::endl;
std::cout<<GridLogMessage <<"MobiusZolotarevFermion test"<<std::endl;
MobiusZolotarevFermion Dzolo(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,b,c,0.1,2.0);
TestCGinversions<MobiusZolotarevFermion>(Dzolo,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5);
std::cout <<"ScaledShamirFermion test"<<std::endl;
std::cout<<GridLogMessage <<"ScaledShamirFermion test"<<std::endl;
ScaledShamirFermion Dsham(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,2.0);
TestCGinversions<ScaledShamirFermion>(Dsham,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5);
std::cout <<"ShamirZolotarevFermion test"<<std::endl;
std::cout<<GridLogMessage <<"ShamirZolotarevFermion test"<<std::endl;
ShamirZolotarevFermion Dshamz(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,0.1,2.0);
TestCGinversions<ShamirZolotarevFermion>(Dshamz,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5);
std::cout <<"OverlapWilsonCayleyTanhFermion test"<<std::endl;
std::cout<<GridLogMessage <<"OverlapWilsonCayleyTanhFermion test"<<std::endl;
OverlapWilsonCayleyTanhFermion Dov(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,1.0);
TestCGinversions<OverlapWilsonCayleyTanhFermion>(Dov,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5);
std::cout <<"OverlapWilsonCayleyZolotarevFermion test"<<std::endl;
std::cout<<GridLogMessage <<"OverlapWilsonCayleyZolotarevFermion test"<<std::endl;
OverlapWilsonCayleyZolotarevFermion Dovz(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,0.1,2.0);
TestCGinversions<OverlapWilsonCayleyZolotarevFermion>(Dovz,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5);
@@ -111,11 +111,11 @@ void TestCGinversions(What & Ddwf,
GridParallelRNG *RNG4,
GridParallelRNG *RNG5)
{
std::cout << "Testing unpreconditioned inverter"<<std::endl;
std::cout<<GridLogMessage << "Testing unpreconditioned inverter"<<std::endl;
TestCGunprec<What>(Ddwf,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,RNG4,RNG5);
std::cout << "Testing red black preconditioned inverter"<<std::endl;
std::cout<<GridLogMessage << "Testing red black preconditioned inverter"<<std::endl;
TestCGprec<What>(Ddwf,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,RNG4,RNG5);
std::cout << "Testing red black Schur inverter"<<std::endl;
std::cout<<GridLogMessage << "Testing red black Schur inverter"<<std::endl;
TestCGschur<What>(Ddwf,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,RNG4,RNG5);
}

36
tests/Test_cayley_coarsen_support.cc
View File

@@ -57,8 +57,8 @@ int main (int argc, char ** argv)
Complex cone(1.0,0.0);
for(int nn=0;nn<Nd;nn++){
if(1) {
if (nn>2) { U[nn]=zero; std::cout << "zeroing gauge field in dir "<<nn<<std::endl; }
else { U[nn]=cone; std::cout << "unit gauge field in dir "<<nn<<std::endl; }
if (nn>2) { U[nn]=zero; std::cout<<GridLogMessage << "zeroing gauge field in dir "<<nn<<std::endl; }
else { U[nn]=cone; std::cout<<GridLogMessage << "unit gauge field in dir "<<nn<<std::endl; }
}
pokeIndex<LorentzIndex>(Umu,U[nn],nn);
}
@@ -75,19 +75,19 @@ int main (int argc, char ** argv)
for(int d=0;d<4;d++){
HermIndefOp.OpDir(src,tmp,d+1,+1); result=result+tmp;
std::cout<<"dir "<<d<<" tmp "<<norm2(tmp)<<std::endl;
std::cout<<GridLogMessage<<"dir "<<d<<" tmp "<<norm2(tmp)<<std::endl;
HermIndefOp.OpDir(src,tmp,d+1,-1); result=result+tmp;
std::cout<<"dir "<<d<<" tmp "<<norm2(tmp)<<std::endl;
std::cout<<GridLogMessage<<"dir "<<d<<" tmp "<<norm2(tmp)<<std::endl;
}
err = result-ref;
std::cout<<"Error "<<norm2(err)<<std::endl;
std::cout<<GridLogMessage<<"Error "<<norm2(err)<<std::endl;
const int nbasis = 2;
LatticeFermion prom(FGrid);
std::vector<LatticeFermion> subspace(nbasis,FGrid);
std::cout<<"Calling Aggregation class" <<std::endl;
std::cout<<GridLogMessage<<"Calling Aggregation class" <<std::endl;
MdagMLinearOperator<DomainWallFermion,LatticeFermion> HermDefOp(Ddwf);
typedef Aggregation<vSpinColourVector,vTComplex,nbasis> Subspace;
@@ -96,7 +96,7 @@ int main (int argc, char ** argv)
subspace=Aggregates.subspace;
std::cout << "Called aggregation class"<< std::endl;
std::cout<<GridLogMessage << "Called aggregation class"<< std::endl;
typedef CoarsenedMatrix<vSpinColourVector,vTComplex,nbasis> LittleDiracOperator;
typedef LittleDiracOperator::CoarseVector CoarseVector;
@@ -119,7 +119,7 @@ int main (int argc, char ** argv)
prom=prom+subspace[b];
}
err=err-prom;
std::cout<<"Promoted back from subspace err "<<norm2(err)<<std::endl;
std::cout<<GridLogMessage<<"Promoted back from subspace err "<<norm2(err)<<std::endl;
HermIndefOp.HermOp(prom,tmp);
blockProject(c_proj,tmp,subspace);
@@ -127,15 +127,15 @@ int main (int argc, char ** argv)
LittleDiracOp.M(c_src,c_res);
c_proj = c_proj - c_res;
std::cout<<"Representation of ldop within subspace "<<norm2(c_proj)<<std::endl;
std::cout<<GridLogMessage<<"Representation of ldop within subspace "<<norm2(c_proj)<<std::endl;
std::cout << "Multiplying by LittleDiracOp "<< std::endl;
std::cout<<GridLogMessage << "Multiplying by LittleDiracOp "<< std::endl;
LittleDiracOp.M(c_src,c_res);
std::cout<<"Testing hermiticity explicitly by inspecting matrix elements"<<std::endl;
std::cout<<GridLogMessage<<"Testing hermiticity explicitly by inspecting matrix elements"<<std::endl;
LittleDiracOp.AssertHermitian();
std::cout << "Testing Hermiticity stochastically "<< std::endl;
std::cout<<GridLogMessage << "Testing Hermiticity stochastically "<< std::endl;
CoarseVector phi(Coarse5d);
CoarseVector chi(Coarse5d);
CoarseVector Aphi(Coarse5d);
@@ -145,7 +145,7 @@ int main (int argc, char ** argv)
random(CRNG,chi);
std::cout<<"Made randoms"<<std::endl;
std::cout<<GridLogMessage<<"Made randoms"<<std::endl;
LittleDiracOp.M(phi,Aphi);
LittleDiracOp.Mdag(chi,Achi);
@@ -155,11 +155,11 @@ int main (int argc, char ** argv)
ComplexD cAc = innerProduct(chi,Achi);
ComplexD pAp = innerProduct(phi,Aphi);
std::cout<< "pAc "<<pAc<<" cAp "<< cAp<< " diff "<<pAc-adj(cAp)<<std::endl;
std::cout<<GridLogMessage<< "pAc "<<pAc<<" cAp "<< cAp<< " diff "<<pAc-adj(cAp)<<std::endl;
std::cout<< "pAp "<<pAp<<" cAc "<< cAc<<"Should be real"<< std::endl;
std::cout<<GridLogMessage<< "pAp "<<pAp<<" cAc "<< cAc<<"Should be real"<< std::endl;
std::cout<<"Testing linearity"<<std::endl;
std::cout<<GridLogMessage<<"Testing linearity"<<std::endl;
CoarseVector PhiPlusChi(Coarse5d);
CoarseVector APhiPlusChi(Coarse5d);
CoarseVector linerr(Coarse5d);
@@ -168,9 +168,9 @@ int main (int argc, char ** argv)
linerr= APhiPlusChi-Aphi;
linerr= linerr-Achi;
std::cout<<"**Diff "<<norm2(linerr)<<std::endl;
std::cout<<GridLogMessage<<"**Diff "<<norm2(linerr)<<std::endl;
std::cout << "Done "<< std::endl;
std::cout<<GridLogMessage << "Done "<< std::endl;
Grid_finalize();
}

84
tests/Test_cayley_even_odd.cc
View File

@@ -29,7 +29,7 @@ int main (int argc, char ** argv)
Grid_init(&argc,&argv);
int threads = GridThread::GetThreads();
std::cout << "Grid is setup to use "<<threads<<" threads"<<std::endl;
std::cout<<GridLogMessage << "Grid is setup to use "<<threads<<" threads"<<std::endl;
const int Ls=8;
GridCartesian * UGrid = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi());
@@ -48,34 +48,34 @@ int main (int argc, char ** argv)
RealD mass=0.1;
RealD M5 =1.8;
std::cout <<"DomainWallFermion test"<<std::endl;
std::cout<<GridLogMessage <<"DomainWallFermion test"<<std::endl;
DomainWallFermion Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5);
TestWhat<DomainWallFermion>(Ddwf,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5);
RealD b=1.5;// Scale factor b+c=2, b-c=1
RealD c=0.5;
std::cout <<"MobiusFermion test"<<std::endl;
std::cout<<GridLogMessage <<"MobiusFermion test"<<std::endl;
MobiusFermion Dmob(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,b,c);
TestWhat<MobiusFermion>(Dmob,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5);
std::cout <<"MobiusZolotarevFermion test"<<std::endl;
std::cout<<GridLogMessage <<"MobiusZolotarevFermion test"<<std::endl;
MobiusZolotarevFermion Dzolo(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,b,c,0.1,2.0);
TestWhat<MobiusZolotarevFermion>(Dzolo,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5);
std::cout <<"ScaledShamirFermion test"<<std::endl;
std::cout<<GridLogMessage <<"ScaledShamirFermion test"<<std::endl;
ScaledShamirFermion Dsham(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,2.0);
TestWhat<ScaledShamirFermion>(Dsham,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5);
std::cout <<"ShamirZolotarevFermion test"<<std::endl;
std::cout<<GridLogMessage <<"ShamirZolotarevFermion test"<<std::endl;
ShamirZolotarevFermion Dshamz(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,0.1,2.0);
TestWhat<ShamirZolotarevFermion>(Dshamz,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5);
std::cout <<"OverlapWilsonCayleyTanhFermion test"<<std::endl;
std::cout<<GridLogMessage <<"OverlapWilsonCayleyTanhFermion test"<<std::endl;
OverlapWilsonCayleyTanhFermion Dov(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,1.0);
TestWhat<OverlapWilsonCayleyTanhFermion>(Dov,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5);
std::cout <<"OverlapWilsonCayleyZolotarevFermion test"<<std::endl;
std::cout<<GridLogMessage <<"OverlapWilsonCayleyZolotarevFermion test"<<std::endl;
OverlapWilsonCayleyZolotarevFermion Dovz(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,0.1,2.0);
TestWhat<OverlapWilsonCayleyZolotarevFermion>(Dovz,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5);
@@ -106,40 +106,40 @@ void TestWhat(What & Ddwf,
LatticeFermion r_eo (FGrid);
LatticeFermion r_eeoo(FGrid);
std::cout<<"=========================================================="<<std::endl;
std::cout<<"= Testing that Meo + Moe + Moo + Mee = Munprec "<<std::endl;
std::cout<<"=========================================================="<<std::endl;
std::cout<<GridLogMessage<<"=========================================================="<<std::endl;
std::cout<<GridLogMessage<<"= Testing that Meo + Moe + Moo + Mee = Munprec "<<std::endl;
std::cout<<GridLogMessage<<"=========================================================="<<std::endl;
pickCheckerboard(Even,src_e,src);
pickCheckerboard(Odd,src_o,src);
Ddwf.Meooe(src_e,r_o); std::cout<<"Applied Meo"<<std::endl;
Ddwf.Meooe(src_o,r_e); std::cout<<"Applied Moe"<<std::endl;
Ddwf.Meooe(src_e,r_o); std::cout<<GridLogMessage<<"Applied Meo"<<std::endl;
Ddwf.Meooe(src_o,r_e); std::cout<<GridLogMessage<<"Applied Moe"<<std::endl;
setCheckerboard(r_eo,r_o);
setCheckerboard(r_eo,r_e);
Ddwf.Mooee(src_e,r_e); std::cout<<"Applied Mee"<<std::endl;
Ddwf.Mooee(src_o,r_o); std::cout<<"Applied Moo"<<std::endl;
Ddwf.Mooee(src_e,r_e); std::cout<<GridLogMessage<<"Applied Mee"<<std::endl;
Ddwf.Mooee(src_o,r_o); std::cout<<GridLogMessage<<"Applied Moo"<<std::endl;
setCheckerboard(r_eeoo,r_e);
setCheckerboard(r_eeoo,r_o);
r_eo=r_eo+r_eeoo;
Ddwf.M(src,ref);
// std::cout << r_eo<<std::endl;
// std::cout << ref <<std::endl;
// std::cout<<GridLogMessage << r_eo<<std::endl;
// std::cout<<GridLogMessage << ref <<std::endl;
err= ref - r_eo;
std::cout << "EO norm diff "<< norm2(err)<< " "<<norm2(ref)<< " " << norm2(r_eo) <<std::endl;
std::cout<<GridLogMessage << "EO norm diff "<< norm2(err)<< " "<<norm2(ref)<< " " << norm2(r_eo) <<std::endl;
LatticeComplex cerr(FGrid);
cerr = localInnerProduct(err,err);
// std::cout << cerr<<std::endl;
// std::cout<<GridLogMessage << cerr<<std::endl;
std::cout<<"=============================================================="<<std::endl;
std::cout<<"= Test Ddagger is the dagger of D by requiring "<<std::endl;
std::cout<<"= < phi | Deo | chi > * = < chi | Deo^dag| phi> "<<std::endl;
std::cout<<"=============================================================="<<std::endl;
std::cout<<GridLogMessage<<"=============================================================="<<std::endl;
std::cout<<GridLogMessage<<"= Test Ddagger is the dagger of D by requiring "<<std::endl;
std::cout<<GridLogMessage<<"= < phi | Deo | chi > * = < chi | Deo^dag| phi> "<<std::endl;
std::cout<<GridLogMessage<<"=============================================================="<<std::endl;
LatticeFermion chi_e (FrbGrid);
LatticeFermion chi_o (FrbGrid);
@@ -169,15 +169,15 @@ void TestWhat(What & Ddwf,
ComplexD cDpe = innerProduct(chi_e,dphi_e);
ComplexD cDpo = innerProduct(chi_o,dphi_o);
std::cout <<"e "<<pDce<<" "<<cDpe <<std::endl;
std::cout <<"o "<<pDco<<" "<<cDpo <<std::endl;
std::cout<<GridLogMessage <<"e "<<pDce<<" "<<cDpe <<std::endl;
std::cout<<GridLogMessage <<"o "<<pDco<<" "<<cDpo <<std::endl;
std::cout <<"pDce - conj(cDpo) "<< pDce-conj(cDpo) <<std::endl;
std::cout <<"pDco - conj(cDpe) "<< pDco-conj(cDpe) <<std::endl;
std::cout<<GridLogMessage <<"pDce - conj(cDpo) "<< pDce-conj(cDpo) <<std::endl;
std::cout<<GridLogMessage <<"pDco - conj(cDpe) "<< pDco-conj(cDpe) <<std::endl;
std::cout<<"=============================================================="<<std::endl;
std::cout<<"= Test MeeInv Mee = 1 "<<std::endl;
std::cout<<"=============================================================="<<std::endl;
std::cout<<GridLogMessage<<"=============================================================="<<std::endl;
std::cout<<GridLogMessage<<"= Test MeeInv Mee = 1 "<<std::endl;
std::cout<<GridLogMessage<<"=============================================================="<<std::endl;
pickCheckerboard(Even,chi_e,chi);
pickCheckerboard(Odd ,chi_o,chi);
@@ -192,11 +192,11 @@ void TestWhat(What & Ddwf,
setCheckerboard(phi,phi_o);
err = phi-chi;
std::cout << "norm diff "<< norm2(err)<< std::endl;
std::cout<<GridLogMessage << "norm diff "<< norm2(err)<< std::endl;
std::cout<<"=============================================================="<<std::endl;
std::cout<<"= Test MeeInvDag MeeDag = 1 "<<std::endl;
std::cout<<"=============================================================="<<std::endl;
std::cout<<GridLogMessage<<"=============================================================="<<std::endl;
std::cout<<GridLogMessage<<"= Test MeeInvDag MeeDag = 1 "<<std::endl;
std::cout<<GridLogMessage<<"=============================================================="<<std::endl;
pickCheckerboard(Even,chi_e,chi);
pickCheckerboard(Odd ,chi_o,chi);
@@ -211,11 +211,11 @@ void TestWhat(What & Ddwf,
setCheckerboard(phi,phi_o);
err = phi-chi;
std::cout << "norm diff "<< norm2(err)<< std::endl;
std::cout<<GridLogMessage << "norm diff "<< norm2(err)<< std::endl;
std::cout<<"=============================================================="<<std::endl;
std::cout<<"= Test DiagMoo MpcDagMpc is Hermitian "<<std::endl;
std::cout<<"=============================================================="<<std::endl;
std::cout<<GridLogMessage<<"=============================================================="<<std::endl;
std::cout<<GridLogMessage<<"= Test DiagMoo MpcDagMpc is Hermitian "<<std::endl;
std::cout<<GridLogMessage<<"=============================================================="<<std::endl;
random(*RNG5,phi);
random(*RNG5,chi);
@@ -237,10 +237,10 @@ void TestWhat(What & Ddwf,
cDpe = innerProduct(chi_e,dphi_e);
cDpo = innerProduct(chi_o,dphi_o);
std::cout <<"e "<<pDce<<" "<<cDpe <<std::endl;
std::cout <<"o "<<pDco<<" "<<cDpo <<std::endl;
std::cout<<GridLogMessage <<"e "<<pDce<<" "<<cDpe <<std::endl;
std::cout<<GridLogMessage <<"o "<<pDco<<" "<<cDpo <<std::endl;
std::cout <<"pDce - conj(cDpo) "<< pDco-conj(cDpo) <<std::endl;
std::cout <<"pDco - conj(cDpe) "<< pDce-conj(cDpe) <<std::endl;
std::cout<<GridLogMessage <<"pDce - conj(cDpo) "<< pDco-conj(cDpo) <<std::endl;
std::cout<<GridLogMessage <<"pDco - conj(cDpe) "<< pDce-conj(cDpe) <<std::endl;
}

30
tests/Test_cayley_ldop_cr.cc
View File

@@ -52,9 +52,9 @@ int main (int argc, char ** argv)
RealD mass=0.1;
RealD M5=1.5;
std::cout << "**************************************************"<< std::endl;
std::cout << "Building g5R5 hermitian DWF operator" <<std::endl;
std::cout << "**************************************************"<< std::endl;
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
std::cout<<GridLogMessage << "Building g5R5 hermitian DWF operator" <<std::endl;
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
DomainWallFermion Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5);
Gamma5R5HermitianLinearOperator<DomainWallFermion,LatticeFermion> HermIndefOp(Ddwf);
@@ -64,9 +64,9 @@ int main (int argc, char ** argv)
typedef CoarsenedMatrix<vSpinColourVector,vTComplex,nbasis> LittleDiracOperator;
typedef LittleDiracOperator::CoarseVector CoarseVector;
std::cout << "**************************************************"<< std::endl;
std::cout << "Calling Aggregation class to build subspace" <<std::endl;
std::cout << "**************************************************"<< std::endl;
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
std::cout<<GridLogMessage << "Calling Aggregation class to build subspace" <<std::endl;
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
MdagMLinearOperator<DomainWallFermion,LatticeFermion> HermDefOp(Ddwf);
Subspace Aggregates(Coarse5d,FGrid);
Aggregates.CreateSubspace(RNG5,HermDefOp);
@@ -80,22 +80,22 @@ int main (int argc, char ** argv)
gaussian(CRNG,c_src);
c_res=zero;
std::cout << "**************************************************"<< std::endl;
std::cout << "Solving mdagm-CG on coarse space "<< std::endl;
std::cout << "**************************************************"<< std::endl;
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
std::cout<<GridLogMessage << "Solving mdagm-CG on coarse space "<< std::endl;
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
MdagMLinearOperator<LittleDiracOperator,CoarseVector> PosdefLdop(LittleDiracOp);
ConjugateGradient<CoarseVector> CG(1.0e-6,10000);
CG(PosdefLdop,c_src,c_res);
std::cout << "**************************************************"<< std::endl;
std::cout << "Solving indef-MCR on coarse space "<< std::endl;
std::cout << "**************************************************"<< std::endl;
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
std::cout<<GridLogMessage << "Solving indef-MCR on coarse space "<< std::endl;
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
HermitianLinearOperator<LittleDiracOperator,CoarseVector> HermIndefLdop(LittleDiracOp);
ConjugateResidual<CoarseVector> MCR(1.0e-6,10000);
MCR(HermIndefLdop,c_src,c_res);
std::cout << "**************************************************"<< std::endl;
std::cout << "Done "<< std::endl;
std::cout << "**************************************************"<< std::endl;
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
std::cout<<GridLogMessage << "Done "<< std::endl;
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
Grid_finalize();
}

12
tests/Test_cf_coarsen_support.cc
View File

@@ -56,12 +56,12 @@ int main (int argc, char ** argv)
for(int d=0;d<4;d++){
HermIndefOp.OpDir(src,tmp,d,+1); result=result+tmp;
std::cout<<"dir "<<d<<" tmp "<<norm2(tmp)<<std::endl;
std::cout<<GridLogMessage<<"dir "<<d<<" tmp "<<norm2(tmp)<<std::endl;
HermIndefOp.OpDir(src,tmp,d,-1); result=result+tmp;
std::cout<<"dir "<<d<<" tmp "<<norm2(tmp)<<std::endl;
std::cout<<GridLogMessage<<"dir "<<d<<" tmp "<<norm2(tmp)<<std::endl;
}
err = result-ref;
std::cout<<"Error "<<norm2(err)<<std::endl;
std::cout<<GridLogMessage<<"Error "<<norm2(err)<<std::endl;
}
{
@@ -73,13 +73,13 @@ int main (int argc, char ** argv)
for(int d=0;d<4;d++){
HermIndefOp.OpDir(src,tmp,d,+1); result=result+tmp;
std::cout<<"dir "<<d<<" tmp "<<norm2(tmp)<<std::endl;
std::cout<<GridLogMessage<<"dir "<<d<<" tmp "<<norm2(tmp)<<std::endl;
HermIndefOp.OpDir(src,tmp,d,-1); result=result+tmp;
std::cout<<"dir "<<d<<" tmp "<<norm2(tmp)<<std::endl;
std::cout<<GridLogMessage<<"dir "<<d<<" tmp "<<norm2(tmp)<<std::endl;
}
err = result-ref;
std::cout<<"Error "<<norm2(err)<<std::endl;
std::cout<<GridLogMessage<<"Error "<<norm2(err)<<std::endl;
}

16
tests/Test_contfrac_cg.cc
View File

@@ -53,7 +53,7 @@ int main (int argc, char ** argv)
Grid_init(&argc,&argv);
int threads = GridThread::GetThreads();
std::cout << "Grid is setup to use "<<threads<<" threads"<<std::endl;
std::cout<<GridLogMessage << "Grid is setup to use "<<threads<<" threads"<<std::endl;
const int Ls=9;
GridCartesian * UGrid = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi());
@@ -74,20 +74,20 @@ int main (int argc, char ** argv)
RealD M5 =1.8;
std::cout <<"OverlapWilsonContFracTanhFermion test"<<std::endl;
std::cout<<GridLogMessage <<"OverlapWilsonContFracTanhFermion test"<<std::endl;
OverlapWilsonContFracTanhFermion Dcf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,1.0);
TestCGinversions<OverlapWilsonContFracTanhFermion>(Dcf,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5);
std::cout <<"OverlapWilsonContFracZolotarevFermion test"<<std::endl;
std::cout<<GridLogMessage <<"OverlapWilsonContFracZolotarevFermion test"<<std::endl;
OverlapWilsonContFracZolotarevFermion Dcfz(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,0.1,6.0);
TestCGinversions<OverlapWilsonContFracZolotarevFermion>(Dcfz,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5);
std::cout <<"OverlapWilsonPartialFractionTanhFermion test"<<std::endl;
std::cout<<GridLogMessage <<"OverlapWilsonPartialFractionTanhFermion test"<<std::endl;
OverlapWilsonPartialFractionTanhFermion Dpf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,1.0);
TestCGinversions<OverlapWilsonPartialFractionTanhFermion>(Dpf,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5);
std::cout <<"OverlapWilsonPartialFractionZolotarevFermion test"<<std::endl;
std::cout<<GridLogMessage <<"OverlapWilsonPartialFractionZolotarevFermion test"<<std::endl;
OverlapWilsonPartialFractionZolotarevFermion Dpfz(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,0.1,6.0);
TestCGinversions<OverlapWilsonPartialFractionZolotarevFermion>(Dpfz,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5);
@@ -102,11 +102,11 @@ void TestCGinversions(What & Ddwf,
GridParallelRNG *RNG4,
GridParallelRNG *RNG5)
{
std::cout << "Testing unpreconditioned inverter"<<std::endl;
std::cout<<GridLogMessage << "Testing unpreconditioned inverter"<<std::endl;
TestCGunprec<What>(Ddwf,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,RNG4,RNG5);
std::cout << "Testing red black preconditioned inverter"<<std::endl;
std::cout<<GridLogMessage << "Testing red black preconditioned inverter"<<std::endl;
TestCGprec<What>(Ddwf,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,RNG4,RNG5);
std::cout << "Testing red black Schur inverter"<<std::endl;
std::cout<<GridLogMessage << "Testing red black Schur inverter"<<std::endl;
TestCGschur<What>(Ddwf,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,RNG4,RNG5);
}

78
tests/Test_contfrac_even_odd.cc
View File

@@ -29,7 +29,7 @@ int main (int argc, char ** argv)
Grid_init(&argc,&argv);
int threads = GridThread::GetThreads();
std::cout << "Grid is setup to use "<<threads<<" threads"<<std::endl;
std::cout<<GridLogMessage << "Grid is setup to use "<<threads<<" threads"<<std::endl;
const int Ls=9;
GridCartesian * UGrid = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi());
@@ -49,19 +49,19 @@ int main (int argc, char ** argv)
RealD mass=0.1;
RealD M5 =1.8;
std::cout <<"OverlapWilsonContFracTanhFermion test"<<std::endl;
std::cout<<GridLogMessage <<"OverlapWilsonContFracTanhFermion test"<<std::endl;
OverlapWilsonContFracTanhFermion Dcf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,1.0);
TestWhat<OverlapWilsonContFracTanhFermion>(Dcf,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5);
std::cout <<"OverlapWilsonContFracZolotarevFermion test"<<std::endl;
std::cout<<GridLogMessage <<"OverlapWilsonContFracZolotarevFermion test"<<std::endl;
OverlapWilsonContFracZolotarevFermion Dcfz(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,0.1,6.0);
TestWhat<OverlapWilsonContFracZolotarevFermion>(Dcfz,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5);
std::cout <<"OverlapWilsonPartialFractionTanhFermion test"<<std::endl;
std::cout<<GridLogMessage <<"OverlapWilsonPartialFractionTanhFermion test"<<std::endl;
OverlapWilsonPartialFractionTanhFermion Dpf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,1.0);
TestWhat<OverlapWilsonPartialFractionTanhFermion>(Dpf,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5);
std::cout <<"OverlapWilsonPartialFractionZolotarevFermion test"<<std::endl;
std::cout<<GridLogMessage <<"OverlapWilsonPartialFractionZolotarevFermion test"<<std::endl;
OverlapWilsonPartialFractionZolotarevFermion Dpfz(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,0.1,6.0);
TestWhat<OverlapWilsonPartialFractionZolotarevFermion>(Dpfz,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5);
@@ -92,40 +92,40 @@ void TestWhat(What & Ddwf,
LatticeFermion r_eo (FGrid);
LatticeFermion r_eeoo(FGrid);
std::cout<<"=========================================================="<<std::endl;
std::cout<<"= Testing that Meo + Moe + Moo + Mee = Munprec "<<std::endl;
std::cout<<"=========================================================="<<std::endl;
std::cout<<GridLogMessage<<"=========================================================="<<std::endl;
std::cout<<GridLogMessage<<"= Testing that Meo + Moe + Moo + Mee = Munprec "<<std::endl;
std::cout<<GridLogMessage<<"=========================================================="<<std::endl;
pickCheckerboard(Even,src_e,src);
pickCheckerboard( Odd,src_o,src);
Ddwf.Meooe(src_e,r_o); std::cout<<"Applied Meo "<<norm2(r_o)<<std::endl;
Ddwf.Meooe(src_o,r_e); std::cout<<"Applied Moe "<<norm2(r_e)<<std::endl;
Ddwf.Meooe(src_e,r_o); std::cout<<GridLogMessage<<"Applied Meo "<<norm2(r_o)<<std::endl;
Ddwf.Meooe(src_o,r_e); std::cout<<GridLogMessage<<"Applied Moe "<<norm2(r_e)<<std::endl;
setCheckerboard(r_eo,r_o);
setCheckerboard(r_eo,r_e);
Ddwf.Mooee(src_e,r_e); std::cout<<"Applied Mee"<<norm2(r_e)<<std::endl;
Ddwf.Mooee(src_o,r_o); std::cout<<"Applied Moo"<<norm2(r_o)<<std::endl;
Ddwf.Mooee(src_e,r_e); std::cout<<GridLogMessage<<"Applied Mee"<<norm2(r_e)<<std::endl;
Ddwf.Mooee(src_o,r_o); std::cout<<GridLogMessage<<"Applied Moo"<<norm2(r_o)<<std::endl;
setCheckerboard(r_eeoo,r_e);
setCheckerboard(r_eeoo,r_o);
r_eo=r_eo+r_eeoo;
Ddwf.M(src,ref);
// std::cout << r_eo<<std::endl;
// std::cout << ref <<std::endl;
// std::cout<<GridLogMessage << r_eo<<std::endl;
// std::cout<<GridLogMessage << ref <<std::endl;
err= ref - r_eo;
std::cout << "EO norm diff "<< norm2(err)<< " "<<norm2(ref)<< " " << norm2(r_eo) <<std::endl;
std::cout<<GridLogMessage << "EO norm diff "<< norm2(err)<< " "<<norm2(ref)<< " " << norm2(r_eo) <<std::endl;
LatticeComplex cerr(FGrid);
cerr = localInnerProduct(err,err);
// std::cout << cerr<<std::endl;
// std::cout<<GridLogMessage << cerr<<std::endl;
std::cout<<"=============================================================="<<std::endl;
std::cout<<"= Test Ddagger is the dagger of D by requiring "<<std::endl;
std::cout<<"= < phi | Deo | chi > * = < chi | Deo^dag| phi> "<<std::endl;
std::cout<<"=============================================================="<<std::endl;
std::cout<<GridLogMessage<<"=============================================================="<<std::endl;
std::cout<<GridLogMessage<<"= Test Ddagger is the dagger of D by requiring "<<std::endl;
std::cout<<GridLogMessage<<"= < phi | Deo | chi > * = < chi | Deo^dag| phi> "<<std::endl;
std::cout<<GridLogMessage<<"=============================================================="<<std::endl;
LatticeFermion chi_e (FrbGrid);
LatticeFermion chi_o (FrbGrid);
@@ -155,15 +155,15 @@ void TestWhat(What & Ddwf,
ComplexD cDpe = innerProduct(chi_e,dphi_e);
ComplexD cDpo = innerProduct(chi_o,dphi_o);
std::cout <<"e "<<pDce<<" "<<cDpe <<std::endl;
std::cout <<"o "<<pDco<<" "<<cDpo <<std::endl;
std::cout<<GridLogMessage <<"e "<<pDce<<" "<<cDpe <<std::endl;
std::cout<<GridLogMessage <<"o "<<pDco<<" "<<cDpo <<std::endl;
std::cout <<"pDce - conj(cDpo) "<< pDce-conj(cDpo) <<std::endl;
std::cout <<"pDco - conj(cDpe) "<< pDco-conj(cDpe) <<std::endl;
std::cout<<GridLogMessage <<"pDce - conj(cDpo) "<< pDce-conj(cDpo) <<std::endl;
std::cout<<GridLogMessage <<"pDco - conj(cDpe) "<< pDco-conj(cDpe) <<std::endl;
std::cout<<"=============================================================="<<std::endl;
std::cout<<"= Test MeeInv Mee = 1 "<<std::endl;
std::cout<<"=============================================================="<<std::endl;
std::cout<<GridLogMessage<<"=============================================================="<<std::endl;
std::cout<<GridLogMessage<<"= Test MeeInv Mee = 1 "<<std::endl;
std::cout<<GridLogMessage<<"=============================================================="<<std::endl;
pickCheckerboard(Even,chi_e,chi);
pickCheckerboard(Odd ,chi_o,chi);
@@ -178,11 +178,11 @@ void TestWhat(What & Ddwf,
setCheckerboard(phi,phi_o);
err = phi-chi;
std::cout << "norm diff "<< norm2(err)<< std::endl;
std::cout<<GridLogMessage << "norm diff "<< norm2(err)<< std::endl;
std::cout<<"=============================================================="<<std::endl;
std::cout<<"= Test MeeInvDag MeeDag = 1 "<<std::endl;
std::cout<<"=============================================================="<<std::endl;
std::cout<<GridLogMessage<<"=============================================================="<<std::endl;
std::cout<<GridLogMessage<<"= Test MeeInvDag MeeDag = 1 "<<std::endl;
std::cout<<GridLogMessage<<"=============================================================="<<std::endl;
pickCheckerboard(Even,chi_e,chi);
pickCheckerboard(Odd ,chi_o,chi);
@@ -197,11 +197,11 @@ void TestWhat(What & Ddwf,
setCheckerboard(phi,phi_o);
err = phi-chi;
std::cout << "norm diff "<< norm2(err)<< std::endl;
std::cout<<GridLogMessage << "norm diff "<< norm2(err)<< std::endl;
std::cout<<"=============================================================="<<std::endl;
std::cout<<"= Test MpcDagMpc is Hermitian "<<std::endl;
std::cout<<"=============================================================="<<std::endl;
std::cout<<GridLogMessage<<"=============================================================="<<std::endl;
std::cout<<GridLogMessage<<"= Test MpcDagMpc is Hermitian "<<std::endl;
std::cout<<GridLogMessage<<"=============================================================="<<std::endl;
random(*RNG5,phi);
random(*RNG5,chi);
@@ -223,10 +223,10 @@ void TestWhat(What & Ddwf,
cDpe = innerProduct(chi_e,dphi_e);
cDpo = innerProduct(chi_o,dphi_o);
std::cout <<"e "<<pDce<<" "<<cDpe <<std::endl;
std::cout <<"o "<<pDco<<" "<<cDpo <<std::endl;
std::cout<<GridLogMessage <<"e "<<pDce<<" "<<cDpe <<std::endl;
std::cout<<GridLogMessage <<"o "<<pDco<<" "<<cDpo <<std::endl;
std::cout <<"pDce - conj(cDpo) "<< pDco-conj(cDpo) <<std::endl;
std::cout <<"pDco - conj(cDpe) "<< pDce-conj(cDpe) <<std::endl;
std::cout<<GridLogMessage <<"pDce - conj(cDpo) "<< pDco-conj(cDpo) <<std::endl;
std::cout<<GridLogMessage <<"pDco - conj(cDpe) "<< pDce-conj(cDpe) <<std::endl;
}

2
tests/Test_cshift.cc
View File

@@ -39,7 +39,7 @@ int main (int argc, char ** argv)
for(int dir=0;dir<4;dir++){
for(int shift=0;shift<latt_size[dir];shift++){
if ( Fine.IsBoss() )
std::cout<<"Shifting by "<<shift<<" in direction"<<dir<<std::endl;
std::cout<<GridLogMessage<<"Shifting by "<<shift<<" in direction"<<dir<<std::endl;
ShiftU = Cshift(U,dir,shift); // Shift everything

30
tests/Test_cshift_red_black.cc
View File

@@ -47,9 +47,9 @@ int main (int argc, char ** argv)
pickCheckerboard(Even,Ue,U);
pickCheckerboard(Odd,Uo,U);
// std::cout << U<<std::endl;
std::cout << "Ue " <<norm2(Ue)<<std::endl;
std::cout << "Uo " <<norm2(Uo)<<std::endl;
// std::cout<<GridLogMessage << U<<std::endl;
std::cout<<GridLogMessage << "Ue " <<norm2(Ue)<<std::endl;
std::cout<<GridLogMessage << "Uo " <<norm2(Uo)<<std::endl;
TComplex cm;
@@ -57,28 +57,28 @@ int main (int argc, char ** argv)
if ( dir!=1 ) continue;
for(int shift=0;shift<latt_size[dir];shift++){
std::cout<<"Shifting by "<<shift<<" in direction"<<dir<<std::endl;
std::cout<<GridLogMessage<<"Shifting by "<<shift<<" in direction"<<dir<<std::endl;
// std::cout<<"Even grid"<<std::endl;
// std::cout<<GridLogMessage<<"Even grid"<<std::endl;
ShiftUe = Cshift(Ue,dir,shift); // Shift everything cb by cb
// std::cout << "\tShiftUe " <<norm2(ShiftUe)<<std::endl;
// std::cout<<GridLogMessage << "\tShiftUe " <<norm2(ShiftUe)<<std::endl;
// std::cout<<"Odd grid"<<std::endl;
// std::cout<<GridLogMessage<<"Odd grid"<<std::endl;
ShiftUo = Cshift(Uo,dir,shift);
// std::cout << "\tShiftUo " <<norm2(ShiftUo)<<std::endl;
// std::cout<<GridLogMessage << "\tShiftUo " <<norm2(ShiftUo)<<std::endl;
// std::cout<<"Recombined Even/Odd grids"<<std::endl;
// std::cout<<GridLogMessage<<"Recombined Even/Odd grids"<<std::endl;
setCheckerboard(rbShiftU,ShiftUe);
setCheckerboard(rbShiftU,ShiftUo);
// std::cout << "\trbShiftU " <<norm2(rbShiftU)<<std::endl;
// std::cout<<GridLogMessage << "\trbShiftU " <<norm2(rbShiftU)<<std::endl;
// std::cout<<"Full grid shift"<<std::endl;
// std::cout<<GridLogMessage<<"Full grid shift"<<std::endl;
ShiftU = Cshift(U,dir,shift); // Shift everything
// std::cout << "\tShiftU " <<norm2(rbShiftU)<<std::endl;
// std::cout<<GridLogMessage << "\tShiftU " <<norm2(rbShiftU)<<std::endl;
std::vector<int> coor(4);
std::cout << "Checking the non-checkerboard shift"<<std::endl;
std::cout<<GridLogMessage << "Checking the non-checkerboard shift"<<std::endl;
for(coor[3]=0;coor[3]<latt_size[3];coor[3]++){
for(coor[2]=0;coor[2]<latt_size[2];coor[2]++){
for(coor[1]=0;coor[1]<latt_size[1];coor[1]++){
@@ -117,7 +117,7 @@ int main (int argc, char ** argv)
}}}}
std::cout << "Checking the checkerboard shift"<<std::endl;
std::cout<<GridLogMessage << "Checking the checkerboard shift"<<std::endl;
for(coor[3]=0;coor[3]<latt_size[3];coor[3]++){
for(coor[2]=0;coor[2]<latt_size[2];coor[2]++){
for(coor[1]=0;coor[1]<latt_size[1];coor[1]++){
@@ -153,7 +153,7 @@ int main (int argc, char ** argv)
std::cerr<<"Expect "<<index<<" " << peer[0]<<","<<peer[1]<<","<<peer[2]<<","<<peer[3]<<std::endl;
exit(-1);
} else if (0) {
std::cout<<"PASS shift "<< shift<<" in dir "<< dir
std::cout<<GridLogMessage<<"PASS shift "<< shift<<" in dir "<< dir
<<" ["<<coor[0]<<","<<coor[1]<<","<<coor[2]<<","<<coor[3]<<"] = "
<< cm()()()<<" expect "<<scm<<" "<<nrm<<std::endl;
}

70
tests/Test_dwf_even_odd.cc
View File

@@ -22,7 +22,7 @@ int main (int argc, char ** argv)
Grid_init(&argc,&argv);
int threads = GridThread::GetThreads();
std::cout << "Grid is setup to use "<<threads<<" threads"<<std::endl;
std::cout<<GridLogMessage << "Grid is setup to use "<<threads<<" threads"<<std::endl;
const int Ls=8;
@@ -67,40 +67,40 @@ int main (int argc, char ** argv)
LatticeFermion r_eo (FGrid);
LatticeFermion r_eeoo(FGrid);
std::cout<<"=========================================================="<<std::endl;
std::cout<<"= Testing that Meo + Moe + Moo + Mee = Munprec "<<std::endl;
std::cout<<"=========================================================="<<std::endl;
std::cout<<GridLogMessage<<"=========================================================="<<std::endl;
std::cout<<GridLogMessage<<"= Testing that Meo + Moe + Moo + Mee = Munprec "<<std::endl;
std::cout<<GridLogMessage<<"=========================================================="<<std::endl;
pickCheckerboard(Even,src_e,src);
pickCheckerboard(Odd,src_o,src);
Ddwf.Meooe(src_e,r_o); std::cout<<"Applied Meo"<<std::endl;
Ddwf.Meooe(src_o,r_e); std::cout<<"Applied Moe"<<std::endl;
Ddwf.Meooe(src_e,r_o); std::cout<<GridLogMessage<<"Applied Meo"<<std::endl;
Ddwf.Meooe(src_o,r_e); std::cout<<GridLogMessage<<"Applied Moe"<<std::endl;
setCheckerboard(r_eo,r_o);
setCheckerboard(r_eo,r_e);
Ddwf.Mooee(src_e,r_e); std::cout<<"Applied Mee"<<std::endl;
Ddwf.Mooee(src_o,r_o); std::cout<<"Applied Moo"<<std::endl;
Ddwf.Mooee(src_e,r_e); std::cout<<GridLogMessage<<"Applied Mee"<<std::endl;
Ddwf.Mooee(src_o,r_o); std::cout<<GridLogMessage<<"Applied Moo"<<std::endl;
setCheckerboard(r_eeoo,r_e);
setCheckerboard(r_eeoo,r_o);
r_eo=r_eo+r_eeoo;
Ddwf.M(src,ref);
// std::cout << r_eo<<std::endl;
// std::cout << ref <<std::endl;
// std::cout<<GridLogMessage << r_eo<<std::endl;
// std::cout<<GridLogMessage << ref <<std::endl;
err= ref - r_eo;
std::cout << "EO norm diff "<< norm2(err)<< " "<<norm2(ref)<< " " << norm2(r_eo) <<std::endl;
std::cout<<GridLogMessage << "EO norm diff "<< norm2(err)<< " "<<norm2(ref)<< " " << norm2(r_eo) <<std::endl;
LatticeComplex cerr(FGrid);
cerr = localInnerProduct(err,err);
// std::cout << cerr<<std::endl;
// std::cout<<GridLogMessage << cerr<<std::endl;
std::cout<<"=============================================================="<<std::endl;
std::cout<<"= Test Ddagger is the dagger of D by requiring "<<std::endl;
std::cout<<"= < phi | Deo | chi > * = < chi | Deo^dag| phi> "<<std::endl;
std::cout<<"=============================================================="<<std::endl;
std::cout<<GridLogMessage<<"=============================================================="<<std::endl;
std::cout<<GridLogMessage<<"= Test Ddagger is the dagger of D by requiring "<<std::endl;
std::cout<<GridLogMessage<<"= < phi | Deo | chi > * = < chi | Deo^dag| phi> "<<std::endl;
std::cout<<GridLogMessage<<"=============================================================="<<std::endl;
LatticeFermion chi_e (FrbGrid);
LatticeFermion chi_o (FrbGrid);
@@ -130,15 +130,15 @@ int main (int argc, char ** argv)
ComplexD cDpe = innerProduct(chi_e,dphi_e);
ComplexD cDpo = innerProduct(chi_o,dphi_o);
std::cout <<"e "<<pDce<<" "<<cDpe <<std::endl;
std::cout <<"o "<<pDco<<" "<<cDpo <<std::endl;
std::cout<<GridLogMessage <<"e "<<pDce<<" "<<cDpe <<std::endl;
std::cout<<GridLogMessage <<"o "<<pDco<<" "<<cDpo <<std::endl;
std::cout <<"pDce - conj(cDpo) "<< pDce-conj(cDpo) <<std::endl;
std::cout <<"pDco - conj(cDpe) "<< pDco-conj(cDpe) <<std::endl;
std::cout<<GridLogMessage <<"pDce - conj(cDpo) "<< pDce-conj(cDpo) <<std::endl;
std::cout<<GridLogMessage <<"pDco - conj(cDpe) "<< pDco-conj(cDpe) <<std::endl;
std::cout<<"=============================================================="<<std::endl;
std::cout<<"= Test MeeInv Mee = 1 "<<std::endl;
std::cout<<"=============================================================="<<std::endl;
std::cout<<GridLogMessage<<"=============================================================="<<std::endl;
std::cout<<GridLogMessage<<"= Test MeeInv Mee = 1 "<<std::endl;
std::cout<<GridLogMessage<<"=============================================================="<<std::endl;
pickCheckerboard(Even,chi_e,chi);
pickCheckerboard(Odd ,chi_o,chi);
@@ -153,11 +153,11 @@ int main (int argc, char ** argv)
setCheckerboard(phi,phi_o);
err = phi-chi;
std::cout << "norm diff "<< norm2(err)<< std::endl;
std::cout<<GridLogMessage << "norm diff "<< norm2(err)<< std::endl;
std::cout<<"=============================================================="<<std::endl;
std::cout<<"= Test MeeInvDag MeeDag = 1 "<<std::endl;
std::cout<<"=============================================================="<<std::endl;
std::cout<<GridLogMessage<<"=============================================================="<<std::endl;
std::cout<<GridLogMessage<<"= Test MeeInvDag MeeDag = 1 "<<std::endl;
std::cout<<GridLogMessage<<"=============================================================="<<std::endl;
pickCheckerboard(Even,chi_e,chi);
pickCheckerboard(Odd ,chi_o,chi);
@@ -172,11 +172,11 @@ int main (int argc, char ** argv)
setCheckerboard(phi,phi_o);
err = phi-chi;
std::cout << "norm diff "<< norm2(err)<< std::endl;
std::cout<<GridLogMessage << "norm diff "<< norm2(err)<< std::endl;
std::cout<<"=============================================================="<<std::endl;
std::cout<<"= Test MpcDagMpc is Hermitian "<<std::endl;
std::cout<<"=============================================================="<<std::endl;
std::cout<<GridLogMessage<<"=============================================================="<<std::endl;
std::cout<<GridLogMessage<<"= Test MpcDagMpc is Hermitian "<<std::endl;
std::cout<<GridLogMessage<<"=============================================================="<<std::endl;
random(RNG5,phi);
random(RNG5,chi);
@@ -199,11 +199,11 @@ int main (int argc, char ** argv)
cDpe = innerProduct(chi_e,dphi_e);
cDpo = innerProduct(chi_o,dphi_o);
std::cout <<"e "<<pDce<<" "<<cDpe <<std::endl;
std::cout <<"o "<<pDco<<" "<<cDpo <<std::endl;
std::cout<<GridLogMessage <<"e "<<pDce<<" "<<cDpe <<std::endl;
std::cout<<GridLogMessage <<"o "<<pDco<<" "<<cDpo <<std::endl;
std::cout <<"pDce - conj(cDpo) "<< pDco-conj(cDpo) <<std::endl;
std::cout <<"pDco - conj(cDpe) "<< pDce-conj(cDpe) <<std::endl;
std::cout<<GridLogMessage <<"pDce - conj(cDpo) "<< pDco-conj(cDpo) <<std::endl;
std::cout<<GridLogMessage <<"pDco - conj(cDpe) "<< pDce-conj(cDpe) <<std::endl;
Grid_finalize();
}

30
tests/Test_dwf_fpgcr.cc
View File

@@ -54,35 +54,35 @@ int main (int argc, char ** argv)
RealD M5=1.8;
DomainWallFermion Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5);
std::cout<<"*********************************************************"<<std::endl;
std::cout<<"* Solving with MdagM VPGCR "<<std::endl;
std::cout<<"*********************************************************"<<std::endl;
std::cout<<GridLogMessage<<"*********************************************************"<<std::endl;
std::cout<<GridLogMessage<<"* Solving with MdagM VPGCR "<<std::endl;
std::cout<<GridLogMessage<<"*********************************************************"<<std::endl;
MdagMLinearOperator<DomainWallFermion,LatticeFermion> HermOp(Ddwf);
result=zero;
PGCR(HermOp,src,result);
std::cout<<"*********************************************************"<<std::endl;
std::cout<<"* Solving with g5-VPGCR "<<std::endl;
std::cout<<"*********************************************************"<<std::endl;
std::cout<<GridLogMessage<<"*********************************************************"<<std::endl;
std::cout<<GridLogMessage<<"* Solving with g5-VPGCR "<<std::endl;
std::cout<<GridLogMessage<<"*********************************************************"<<std::endl;
Gamma5R5HermitianLinearOperator<DomainWallFermion,LatticeFermion> g5HermOp(Ddwf);
result=zero;
PGCR(g5HermOp,src,result);
std::cout<<"*********************************************************"<<std::endl;
std::cout<<"* Solving with MdagM-CR "<<std::endl;
std::cout<<"*********************************************************"<<std::endl;
std::cout<<GridLogMessage<<"*********************************************************"<<std::endl;
std::cout<<GridLogMessage<<"* Solving with MdagM-CR "<<std::endl;
std::cout<<GridLogMessage<<"*********************************************************"<<std::endl;
result=zero;
CR(HermOp,src,result);
std::cout<<"*********************************************************"<<std::endl;
std::cout<<"* Solving with g5-CR "<<std::endl;
std::cout<<"*********************************************************"<<std::endl;
std::cout<<GridLogMessage<<"*********************************************************"<<std::endl;
std::cout<<GridLogMessage<<"* Solving with g5-CR "<<std::endl;
std::cout<<GridLogMessage<<"*********************************************************"<<std::endl;
result=zero;
CR(g5HermOp,src,result);
std::cout<<"*********************************************************"<<std::endl;
std::cout<<"* Solving with MdagM-CG "<<std::endl;
std::cout<<"*********************************************************"<<std::endl;
std::cout<<GridLogMessage<<"*********************************************************"<<std::endl;
std::cout<<GridLogMessage<<"* Solving with MdagM-CG "<<std::endl;
std::cout<<GridLogMessage<<"*********************************************************"<<std::endl;
result=zero;
CG(HermOp,src,result);

159
tests/Test_dwf_hdcr.cc
View File

@@ -53,7 +53,7 @@ public:
// _FineOperator.Op(p1,p2);// this is the G5 herm bit
RealD absp2=std::sqrt(norm2(p2));
if(i%10==9)
std::cout << "Power method on mdagm "<<i<<" " << absp2/absp1<<std::endl;
std::cout<<GridLogMessage << "Power method on mdagm "<<i<<" " << absp2/absp1<<std::endl;
p1=p2*(1.0/std::sqrt(absp2));
}
}
@@ -71,7 +71,7 @@ public:
// Monitor completeness of low mode space
_Aggregates.ProjectToSubspace (Csrc,in);
_Aggregates.PromoteFromSubspace(Csrc,out);
std::cout<<"Completeness: "<<std::sqrt(norm2(out)/norm2(in))<<std::endl;
std::cout<<GridLogMessage<<"Completeness: "<<std::sqrt(norm2(out)/norm2(in))<<std::endl;
// Build some solvers
ConjugateGradient<FineField> fCG(1.0e-3,1000);
@@ -85,21 +85,21 @@ public:
MdagMLinearOperator<Matrix,FineField> MdagMOp(_Matrix);
Min=in;
std::cout<< " Preconditioner in " << norm2(in)<<std::endl;
std::cout<<GridLogMessage<< " Preconditioner in " << norm2(in)<<std::endl;
_FineOperator.AdjOp(Min,tmp);
std::cout<< " Preconditioner tmp " << norm2(in)<<std::endl;
std::cout<<GridLogMessage<< " Preconditioner tmp " << norm2(in)<<std::endl;
fCG(MdagMOp,tmp,out);
_FineOperator.Op(out,tmp);
std::cout<< " Preconditioner in " << norm2(in)<<std::endl;
std::cout<< " Preconditioner out " << norm2(out)<<std::endl;
std::cout<< " Preconditioner Aout" << norm2(tmp)<<std::endl;
std::cout<<GridLogMessage<< " Preconditioner in " << norm2(in)<<std::endl;
std::cout<<GridLogMessage<< " Preconditioner out " << norm2(out)<<std::endl;
std::cout<<GridLogMessage<< " Preconditioner Aout" << norm2(tmp)<<std::endl;
tmp = tmp - in;
std::cout<<"preconditioner thinks residual is "<<std::sqrt(norm2(tmp)/norm2(in))<<std::endl;
std::cout<<GridLogMessage<<"preconditioner thinks residual is "<<std::sqrt(norm2(tmp)/norm2(in))<<std::endl;
/*
// _FineOperator.Op(Min,out);
@@ -145,7 +145,7 @@ public:
// Monitor completeness of low mode space
_Aggregates.ProjectToSubspace (Csrc,in);
_Aggregates.PromoteFromSubspace(Csrc,out);
std::cout<<"Coarse Grid Preconditioner\nCompleteness in: "<<std::sqrt(norm2(out)/norm2(in))<<std::endl;
std::cout<<GridLogMessage<<"Coarse Grid Preconditioner\nCompleteness in: "<<std::sqrt(norm2(out)/norm2(in))<<std::endl;
// [PTM+Q] in = [1 - Q A] M in + Q in = Min + Q [ in -A Min]
_FineOperator.Op(in,tmp);// this is the G5 herm bit
@@ -154,7 +154,7 @@ public:
// Monitor completeness of low mode space
_Aggregates.ProjectToSubspace (Csrc,Min);
_Aggregates.PromoteFromSubspace(Csrc,out);
std::cout<<"Completeness Min: "<<std::sqrt(norm2(out)/norm2(Min))<<std::endl;
std::cout<<GridLogMessage<<"Completeness Min: "<<std::sqrt(norm2(out)/norm2(Min))<<std::endl;
_FineOperator.Op(Min,tmp);
tmp = in - tmp; // in - A Min
@@ -166,14 +166,14 @@ public:
HermOp.Op(Csol,Ctmp);
Ctmp=Ctmp-Csrc;
std::cout<<"coarse space true residual "<<std::sqrt(norm2(Ctmp)/norm2(Csrc))<<std::endl;
std::cout<<GridLogMessage<<"coarse space true residual "<<std::sqrt(norm2(Ctmp)/norm2(Csrc))<<std::endl;
_Aggregates.PromoteFromSubspace(Csol,out);
_FineOperator.Op(out,res);
res=res-tmp;
std::cout<<"promoted sol residual "<<std::sqrt(norm2(res)/norm2(tmp))<<std::endl;
std::cout<<GridLogMessage<<"promoted sol residual "<<std::sqrt(norm2(res)/norm2(tmp))<<std::endl;
_Aggregates.ProjectToSubspace (Csrc,res);
std::cout<<"coarse space proj of residual "<<norm2(Csrc)<<std::endl;
std::cout<<GridLogMessage<<"coarse space proj of residual "<<norm2(Csrc)<<std::endl;
out = out+Min; // additive coarse space correction
@@ -182,9 +182,9 @@ public:
_FineOperator.Op(out,tmp);
tmp=tmp-in; // tmp is new residual
std::cout<< " Preconditioner in " << norm2(in)<<std::endl;
std::cout<< " Preconditioner out " << norm2(out)<<std::endl;
std::cout<<"preconditioner thinks residual is "<<std::sqrt(norm2(tmp)/norm2(in))<<std::endl;
std::cout<<GridLogMessage<< " Preconditioner in " << norm2(in)<<std::endl;
std::cout<<GridLogMessage<< " Preconditioner out " << norm2(out)<<std::endl;
std::cout<<GridLogMessage<<"preconditioner thinks residual is "<<std::sqrt(norm2(tmp)/norm2(in))<<std::endl;
}
#endif
@@ -210,7 +210,7 @@ public:
// Monitor completeness of low mode space
// _Aggregates.ProjectToSubspace (Csrc,in);
// _Aggregates.PromoteFromSubspace(Csrc,out);
// std::cout<<"Coarse Grid Preconditioner\nCompleteness in: "<<std::sqrt(norm2(out)/norm2(in))<<std::endl;
// std::cout<<GridLogMessage<<"Coarse Grid Preconditioner\nCompleteness in: "<<std::sqrt(norm2(out)/norm2(in))<<std::endl;
_Aggregates.ProjectToSubspace (Csrc,in);
HermOp.AdjOp(Csrc,Ctmp);// Normal equations
@@ -229,7 +229,7 @@ public:
_FineOperator.Op(out,tmp);
tmp=tmp-in; // tmp is new residual
std::cout<<"preconditioner thinks residual is "<<std::sqrt(norm2(tmp)/norm2(in))<<std::endl;
std::cout<<GridLogMessage<<"preconditioner thinks residual is "<<std::sqrt(norm2(tmp)/norm2(in))<<std::endl;
}
#endif
@@ -247,8 +247,8 @@ public:
Ni = norm2(in);
for(int ilo=0;ilo<4;ilo++){
for(int ord=10;ord<60;ord+=10){
for(int ilo=0;ilo<3;ilo++){
for(int ord=5;ord<50;ord*=2){
_FineOperator.AdjOp(in,vec1);
@@ -258,7 +258,7 @@ public:
_FineOperator.Op(vec2,vec1);// this is the G5 herm bit
vec1 = in - vec1; // tmp = in - A Min
r=norm2(vec1);
std::cout << "Smoother resid "<<std::sqrt(r/Ni)<<std::endl;
std::cout<<GridLogMessage << "Smoother resid "<<std::sqrt(r/Ni)<<std::endl;
}
}
@@ -270,13 +270,13 @@ public:
CoarseVector Ctmp(_CoarseOperator.Grid());
CoarseVector Csol(_CoarseOperator.Grid()); Csol=zero;
ConjugateGradient<CoarseVector> CG(1.0e-5,100000);
ConjugateGradient<FineField> fCG(3.0e-2,1000);
ConjugateGradient<CoarseVector> CG(1.0e-3,100000);
// ConjugateGradient<FineField> fCG(3.0e-2,1000);
HermitianLinearOperator<CoarseOperator,CoarseVector> HermOp(_CoarseOperator);
MdagMLinearOperator<CoarseOperator,CoarseVector> MdagMOp(_CoarseOperator);
// MdagMLinearOperator<Matrix,FineField> fMdagMOp(_Matrix);
ShiftedMdagMLinearOperator<Matrix,FineField> fMdagMOp(_Matrix,1.0);
ShiftedMdagMLinearOperator<Matrix,FineField> fMdagMOp(_Matrix,0.0);
FineField vec1(in._grid);
FineField vec2(in._grid);
@@ -288,7 +288,7 @@ public:
_Aggregates.ProjectToSubspace (Csrc,in);
_Aggregates.PromoteFromSubspace(Csrc,out);
std::cout<<"Completeness: "<<std::sqrt(norm2(out)/norm2(in))<<std::endl;
std::cout<<GridLogMessage<<"Completeness: "<<std::sqrt(norm2(out)/norm2(in))<<std::endl;
// ofstream fout("smoother");
// Cheby.csv(fout);
@@ -306,7 +306,7 @@ public:
_FineOperator.AdjOp(in,vec1);// this is the G5 herm bit
ChebyAccu(fMdagMOp,vec1,out); // solves MdagM = g5 M g5M
std::cout << "Smoother norm "<<norm2(out)<<std::endl;
std::cout<<GridLogMessage << "Smoother norm "<<norm2(out)<<std::endl;
// Update with residual for out
_FineOperator.Op(out,vec1);// this is the G5 herm bit
@@ -314,7 +314,7 @@ public:
RealD r = norm2(vec1);
std::cout << "Smoother resid "<<std::sqrt(r/Ni)<< " " << r << " " << Ni <<std::endl;
std::cout<<GridLogMessage << "Smoother resid "<<std::sqrt(r/Ni)<< " " << r << " " << Ni <<std::endl;
_Aggregates.ProjectToSubspace (Csrc,vec1);
HermOp.AdjOp(Csrc,Ctmp);// Normal equations
@@ -329,7 +329,7 @@ public:
vec1 = in - vec1; // tmp = in - A Min
r=norm2(vec1);
std::cout << "Coarse resid "<<std::sqrt(r/Ni)<<std::endl;
std::cout<<GridLogMessage << "Coarse resid "<<std::sqrt(r/Ni)<<std::endl;
// Reapply smoother
_FineOperator.Op(vec1,vec2); // this is the G5 herm bit
@@ -339,7 +339,7 @@ public:
_FineOperator.Op(out,vec1);// this is the G5 herm bit
vec1 = in - vec1; // tmp = in - A Min
r=norm2(vec1);
std::cout << "Smoother resid "<<std::sqrt(r/Ni)<<std::endl;
std::cout<<GridLogMessage << "Smoother resid "<<std::sqrt(r/Ni)<<std::endl;
}
@@ -385,7 +385,7 @@ int main (int argc, char ** argv)
LatticeGaugeField Umu(UGrid);
NerscField header;
std::string file("./ckpoint_lat.4000");
std::string file("./ckpoint_lat.400");
readNerscConfiguration(Umu,header,file);
// SU3::ColdConfiguration(RNG4,Umu);
@@ -396,9 +396,9 @@ int main (int argc, char ** argv)
RealD mass=0.01;
RealD M5=1.8;
std::cout << "**************************************************"<< std::endl;
std::cout << "Building g5R5 hermitian DWF operator" <<std::endl;
std::cout << "**************************************************"<< std::endl;
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
std::cout<<GridLogMessage << "Building g5R5 hermitian DWF operator" <<std::endl;
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
DomainWallFermion Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5);
const int nbasis = 32;
@@ -408,22 +408,22 @@ int main (int argc, char ** argv)
typedef CoarsenedMatrix<vSpinColourVector,vTComplex,nbasis> CoarseOperator;
typedef CoarseOperator::CoarseVector CoarseVector;
std::cout << "**************************************************"<< std::endl;
std::cout << "Calling Aggregation class to build subspace" <<std::endl;
std::cout << "**************************************************"<< std::endl;
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
std::cout<<GridLogMessage << "Calling Aggregation class to build subspace" <<std::endl;
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
MdagMLinearOperator<DomainWallFermion,LatticeFermion> HermDefOp(Ddwf);
Subspace Aggregates(Coarse5d,FGrid);
// Aggregates.CreateSubspace(RNG5,HermDefOp,nbasis);
assert ( (nbasis & 0x1)==0);
int nb=nbasis/2;
std::cout << " nbasis/2 = "<<nb<<std::endl;
std::cout<<GridLogMessage << " nbasis/2 = "<<nb<<std::endl;
Aggregates.CreateSubspace(RNG5,HermDefOp,nb);
for(int n=0;n<nb;n++){
G5R5(Aggregates.subspace[n+nb],Aggregates.subspace[n]);
std::cout<<n<<" subspace "<<norm2(Aggregates.subspace[n+nb])<<" "<<norm2(Aggregates.subspace[n]) <<std::endl;
std::cout<<GridLogMessage<<n<<" subspace "<<norm2(Aggregates.subspace[n+nb])<<" "<<norm2(Aggregates.subspace[n]) <<std::endl;
}
for(int n=0;n<nbasis;n++){
std::cout << "vec["<<n<<"] = "<<norm2(Aggregates.subspace[n]) <<std::endl;
std::cout<<GridLogMessage << "vec["<<n<<"] = "<<norm2(Aggregates.subspace[n]) <<std::endl;
}
// for(int i=0;i<nbasis;i++){
@@ -432,77 +432,90 @@ int main (int argc, char ** argv)
// }
result=zero;
std::cout << "**************************************************"<< std::endl;
std::cout << "Building coarse representation of Indef operator" <<std::endl;
std::cout << "**************************************************"<< std::endl;
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
std::cout<<GridLogMessage << "Building coarse representation of Indef operator" <<std::endl;
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
Gamma5R5HermitianLinearOperator<DomainWallFermion,LatticeFermion> HermIndefOp(Ddwf);
CoarsenedMatrix<vSpinColourVector,vTComplex,nbasis> LDOp(*Coarse5d);
LDOp.CoarsenOperator(FGrid,HermIndefOp,Aggregates);
std::cout << "**************************************************"<< std::endl;
std::cout << "Testing some coarse space solvers " <<std::endl;
std::cout << "**************************************************"<< std::endl;
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
std::cout<<GridLogMessage << "Testing some coarse space solvers " <<std::endl;
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
CoarseVector c_src (Coarse5d);
CoarseVector c_res (Coarse5d);
gaussian(CRNG,c_src);
c_res=zero;
std::cout << "**************************************************"<< std::endl;
std::cout << "Solving posdef-CG on coarse space "<< std::endl;
std::cout << "**************************************************"<< std::endl;
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
std::cout<<GridLogMessage << "Solving posdef-CG on coarse space "<< std::endl;
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
MdagMLinearOperator<CoarseOperator,CoarseVector> PosdefLdop(LDOp);
ConjugateGradient<CoarseVector> CG(1.0e-6,100000);
CG(PosdefLdop,c_src,c_res);
// std::cout << "**************************************************"<< std::endl;
// std::cout << "Solving indef-MCR on coarse space "<< std::endl;
// std::cout << "**************************************************"<< std::endl;
// std::cout<<GridLogMessage << "**************************************************"<< std::endl;
// std::cout<<GridLogMessage << "Solving indef-MCR on coarse space "<< std::endl;
// std::cout<<GridLogMessage << "**************************************************"<< std::endl;
// HermitianLinearOperator<CoarseOperator,CoarseVector> HermIndefLdop(LDOp);
// ConjugateResidual<CoarseVector> MCR(1.0e-6,100000);
//MCR(HermIndefLdop,c_src,c_res);
std::cout << "**************************************************"<< std::endl;
std::cout << "Building deflation preconditioner "<< std::endl;
std::cout << "**************************************************"<< std::endl;
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
std::cout<<GridLogMessage << "Building deflation preconditioner "<< std::endl;
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
MultiGridPreconditioner <vSpinColourVector,vTComplex,nbasis,DomainWallFermion> Precon(Aggregates, LDOp,HermIndefOp,Ddwf);
TrivialPrecon<LatticeFermion> simple;
std::cout << "**************************************************"<< std::endl;
std::cout << "Testing smoother efficacy"<< std::endl;
std::cout << "**************************************************"<< std::endl;
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
std::cout<<GridLogMessage << "Testing smoother efficacy"<< std::endl;
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
Precon.SmootherTest(src);
std::cout << "**************************************************"<< std::endl;
std::cout << "Unprec CG "<< std::endl;
std::cout << "**************************************************"<< std::endl;
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
std::cout<<GridLogMessage << "Unprec CG "<< std::endl;
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
// TrivialPrecon<LatticeFermion> simple;
// ConjugateGradient<LatticeFermion> fCG(1.0e-8,100000);
// fCG(HermDefOp,src,result);
// exit(0);
std::cout << "**************************************************"<< std::endl;
std::cout << "Testing GCR on indef matrix "<< std::endl;
std::cout << "**************************************************"<< std::endl;
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
std::cout<<GridLogMessage << "Testing GCR on indef matrix "<< std::endl;
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
// PrecGeneralisedConjugateResidual<LatticeFermion> UPGCR(1.0e-8,100000,simple,8,128);
// UPGCR(HermIndefOp,src,result);
/// Get themax eval
std::cout << "**************************************************"<< std::endl;
std::cout <<" Applying power method to find spectral range "<<std::endl;
std::cout << "**************************************************"<< std::endl;
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
std::cout<<GridLogMessage <<" Applying power method to find spectral range "<<std::endl;
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
Precon.PowerMethod(src);
std::cout << "**************************************************"<< std::endl;
std::cout << "Building a two level PGCR "<< std::endl;
std::cout << "**************************************************"<< std::endl;
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
std::cout<<GridLogMessage << "Building a two level PGCR "<< std::endl;
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
PrecGeneralisedConjugateResidual<LatticeFermion> PGCR(1.0e-8,100000,Precon,8,128);
std::cout<<"checking norm src "<<norm2(src)<<std::endl;
std::cout<<GridLogMessage<<"checking norm src "<<norm2(src)<<std::endl;
PGCR(HermIndefOp,src,result);
std::cout << "**************************************************"<< std::endl;
std::cout << "Done "<< std::endl;
std::cout << "**************************************************"<< std::endl;
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
std::cout<<GridLogMessage << "Red Black Prec CG "<< std::endl;
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
SchurDiagMooeeOperator<DomainWallFermion,LatticeFermion> HermOpEO(Ddwf);
ConjugateGradient<LatticeFermion> pCG(1.0e-8,10000);
LatticeFermion src_o(FrbGrid);
LatticeFermion result_o(FrbGrid);
pickCheckerboard(Odd,src_o,src);
result_o=zero;
pCG(HermOpEO,src_o,result_o);
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
std::cout<<GridLogMessage << "Done "<< std::endl;
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
Grid_finalize();
}

30
tests/Test_gamma.cc
View File

@@ -35,19 +35,19 @@ int main (int argc, char ** argv)
SpinVector lv; random(sRNG,lv);
SpinVector rv; random(sRNG,rv);
// std::cout << " Is pod " << std::is_pod<SpinVector>::value << std::endl;
// std::cout << " Is pod double " << std::is_pod<double>::value << std::endl;
// std::cout << " Is pod ComplexF " << std::is_pod<ComplexF>::value << std::endl;
// std::cout << " Is triv double " << std::has_trivial_default_constructor<double>::value << std::endl;
// std::cout << " Is triv ComplexF " << std::has_trivial_default_constructor<ComplexF>::value << std::endl;
// std::cout << " Is pod Scalar<double> " << std::is_pod<iScalar<double> >::value << std::endl;
// std::cout << " Is pod Scalar<ComplexF> " << std::is_pod<iScalar<ComplexF> >::value << std::endl;
// std::cout << " Is pod Scalar<vComplexF> " << std::is_pod<iScalar<vComplexF> >::value << std::endl;
// std::cout << " Is pod Scalar<vComplexD> " << std::is_pod<iScalar<vComplexD> >::value << std::endl;
// std::cout << " Is pod Scalar<vRealF> " << std::is_pod<iScalar<vRealF> >::value << std::endl;
// std::cout << " Is pod Scalar<vRealD> " << std::is_pod<iScalar<vRealD> >::value << std::endl;
// std::cout << " Is triv Scalar<double> " <<std::has_trivial_default_constructor<iScalar<double> >::value << std::endl;
// std::cout << " Is triv Scalar<vComplexD> "<<std::has_trivial_default_constructor<iScalar<vComplexD> >::value << std::endl;
// std::cout<<GridLogMessage << " Is pod " << std::is_pod<SpinVector>::value << std::endl;
// std::cout<<GridLogMessage << " Is pod double " << std::is_pod<double>::value << std::endl;
// std::cout<<GridLogMessage << " Is pod ComplexF " << std::is_pod<ComplexF>::value << std::endl;
// std::cout<<GridLogMessage << " Is triv double " << std::has_trivial_default_constructor<double>::value << std::endl;
// std::cout<<GridLogMessage << " Is triv ComplexF " << std::has_trivial_default_constructor<ComplexF>::value << std::endl;
// std::cout<<GridLogMessage << " Is pod Scalar<double> " << std::is_pod<iScalar<double> >::value << std::endl;
// std::cout<<GridLogMessage << " Is pod Scalar<ComplexF> " << std::is_pod<iScalar<ComplexF> >::value << std::endl;
// std::cout<<GridLogMessage << " Is pod Scalar<vComplexF> " << std::is_pod<iScalar<vComplexF> >::value << std::endl;
// std::cout<<GridLogMessage << " Is pod Scalar<vComplexD> " << std::is_pod<iScalar<vComplexD> >::value << std::endl;
// std::cout<<GridLogMessage << " Is pod Scalar<vRealF> " << std::is_pod<iScalar<vRealF> >::value << std::endl;
// std::cout<<GridLogMessage << " Is pod Scalar<vRealD> " << std::is_pod<iScalar<vRealD> >::value << std::endl;
// std::cout<<GridLogMessage << " Is triv Scalar<double> " <<std::has_trivial_default_constructor<iScalar<double> >::value << std::endl;
// std::cout<<GridLogMessage << " Is triv Scalar<vComplexD> "<<std::has_trivial_default_constructor<iScalar<vComplexD> >::value << std::endl;
for(int a=0;a<Ns;a++){
ident()(a,a) = ComplexF(1.0);
@@ -66,8 +66,8 @@ int main (int argc, char ** argv)
for(int i=0;i<Ns;i++){
if(i==0) std::cout << list[mu];
else std::cout << list[12];
if(i==0) std::cout<<GridLogMessage << list[mu];
else std::cout<<GridLogMessage << list[12];
std::cout<<"(";
for(int j=0;j<Ns;j++){

24
tests/Test_lie_generators.cc
View File

@@ -20,27 +20,27 @@ int main (int argc, char ** argv)
GridRedBlackCartesian * rbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(grid);
std::cout<<"*********************************************"<<std::endl;
std::cout<<"* Generators for SU(2)"<<std::endl;
std::cout<<"*********************************************"<<std::endl;
std::cout<<GridLogMessage<<"*********************************************"<<std::endl;
std::cout<<GridLogMessage<<"* Generators for SU(2)"<<std::endl;
std::cout<<GridLogMessage<<"*********************************************"<<std::endl;
SU2::printGenerators();
SU2::testGenerators();
std::cout<<"*********************************************"<<std::endl;
std::cout<<"* Generators for SU(3)"<<std::endl;
std::cout<<"*********************************************"<<std::endl;
std::cout<<GridLogMessage<<"*********************************************"<<std::endl;
std::cout<<GridLogMessage<<"* Generators for SU(3)"<<std::endl;
std::cout<<GridLogMessage<<"*********************************************"<<std::endl;
SU3::printGenerators();
SU3::testGenerators();
// std::cout<<"*********************************************"<<std::endl;
// std::cout<<"* Generators for SU(4)"<<std::endl;
// std::cout<<"*********************************************"<<std::endl;
// std::cout<<GridLogMessage<<"*********************************************"<<std::endl;
// std::cout<<GridLogMessage<<"* Generators for SU(4)"<<std::endl;
// std::cout<<GridLogMessage<<"*********************************************"<<std::endl;
// SU4::printGenerators();
// SU4::testGenerators();
// std::cout<<"*********************************************"<<std::endl;
// std::cout<<"* Generators for SU(5)"<<std::endl;
// std::cout<<"*********************************************"<<std::endl;
// std::cout<<GridLogMessage<<"*********************************************"<<std::endl;
// std::cout<<GridLogMessage<<"* Generators for SU(5)"<<std::endl;
// std::cout<<GridLogMessage<<"*********************************************"<<std::endl;
// SU5::printGenerators();
// SU5::testGenerators();

42
tests/Test_main.cc
View File

@@ -213,29 +213,29 @@ int main (int argc, char ** argv)
random(SerialRNG, cm);
std::cout << cm << std::endl;
std::cout<<GridLogMessage << cm << std::endl;
cm = Ta(cm);
TComplex tracecm= trace(cm);
std::cout << cm << std::endl;
std::cout<<GridLogMessage << cm << std::endl;
cm = Exponentiate(cm, 2.0, 12);
std::cout << cm << " " << std::endl;
std::cout<<GridLogMessage << cm << " " << std::endl;
Complex det = Determinant(cm);
std::cout << "determinant: " << det << std::endl;
std::cout << "norm: " << norm2(cm) << std::endl;
std::cout<<GridLogMessage << "determinant: " << det << std::endl;
std::cout<<GridLogMessage << "norm: " << norm2(cm) << std::endl;
cm = ProjectOnGroup(cm);
std::cout << cm << " " << std::endl;
std::cout << "norm: " << norm2(cm) << std::endl;
std::cout<<GridLogMessage << cm << " " << std::endl;
std::cout<<GridLogMessage << "norm: " << norm2(cm) << std::endl;
cm = ProjectOnGroup(cm);
std::cout << cm << " " << std::endl;
std::cout << "norm: " << norm2(cm) << std::endl;
std::cout<<GridLogMessage << cm << " " << std::endl;
std::cout<<GridLogMessage << "norm: " << norm2(cm) << std::endl;
// det = Determinant(cm);
// std::cout << "determinant: " << det << std::endl;
// std::cout<<GridLogMessage << "determinant: " << det << std::endl;
// Foo = Foo+scalar; // LatticeColourMatrix+Scalar
@@ -253,19 +253,19 @@ int main (int argc, char ** argv)
random(FineRNG,cMat);
cMat = Ta(cMat);
peekSite(cm, cMat, mysite);
std::cout << cm << " " << std::endl;
std::cout<<GridLogMessage << cm << " " << std::endl;
cm = Exponentiate(cm, 1.0, 12);
std::cout << cm << " " << std::endl;
std::cout << "norm: " << norm2(cm) << std::endl;
std::cout<<GridLogMessage << cm << " " << std::endl;
std::cout<<GridLogMessage << "norm: " << norm2(cm) << std::endl;
std::cout << "norm cMmat : " << norm2(cMat) << std::endl;
std::cout<<GridLogMessage << "norm cMmat : " << norm2(cMat) << std::endl;
cMat = expMat(cMat, ComplexD(1.0, 0.0));
std::cout << "norm expMat: " << norm2(cMat) << std::endl;
std::cout<<GridLogMessage << "norm expMat: " << norm2(cMat) << std::endl;
peekSite(cm, cMat, mysite);
std::cout << cm << " " << std::endl;
std::cout << "determinant: " << Determinant(cm) << std::endl;
std::cout << "norm: " << norm2(cm) << std::endl;
std::cout<<GridLogMessage << cm << " " << std::endl;
std::cout<<GridLogMessage << "determinant: " << Determinant(cm) << std::endl;
std::cout<<GridLogMessage << "norm: " << norm2(cm) << std::endl;
// LatticeComplex trlcMat(&Fine);
@@ -434,7 +434,7 @@ int main (int argc, char ** argv)
pickCheckerboard(1,bFoo,Foo);
if ( Fine.IsBoss() ) {
std::cout << "Shifting both parities by "<< shift <<" direction "<< dir <<std::endl;
std::cout<<GridLogMessage << "Shifting both parities by "<< shift <<" direction "<< dir <<std::endl;
}
Shifted = Cshift(Foo,dir,shift); // Shift everything
@@ -545,7 +545,7 @@ int main (int argc, char ** argv)
}}
}}}}
if( Fine.IsBoss() ){
std::cout << "LatticeColorMatrix * LatticeColorMatrix nrm diff = "<<nrm<<std::endl;
std::cout<<GridLogMessage << "LatticeColorMatrix * LatticeColorMatrix nrm diff = "<<nrm<<std::endl;
}
}}
@@ -553,7 +553,7 @@ int main (int argc, char ** argv)
} // loop for omp
std::cout << sizeof(vComplexF) << std::endl;
std::cout<<GridLogMessage << sizeof(vComplexF) << std::endl;
Grid_finalize();
}

20
tests/Test_multishift_sqrt.cc
View File

@@ -82,7 +82,7 @@ int main (int argc, char ** argv)
////////////////////////////////////////
// sqrt and inverse sqrt
////////////////////////////////////////
std::cout << "Generating degree "<<degree<<" for x^(1/2)"<<std::endl;
std::cout<<GridLogMessage << "Generating degree "<<degree<<" for x^(1/2)"<<std::endl;
remez.generateApprox(degree,1,2);
MultiShiftFunction Sqrt(remez,1.0e-6,false);
@@ -120,19 +120,19 @@ int main (int argc, char ** argv)
error = reference - combined;
std::cout << " Reference "<<norm2(reference)<<std::endl;
std::cout << " combined "<<norm2(combined) <<std::endl;
std::cout << " error "<<norm2(error) <<std::endl;
std::cout<<GridLogMessage << " Reference "<<norm2(reference)<<std::endl;
std::cout<<GridLogMessage << " combined "<<norm2(combined) <<std::endl;
std::cout<<GridLogMessage << " error "<<norm2(error) <<std::endl;
MSCG(Diagonal,src,summed);
error = summed - combined;
std::cout << " summed-combined "<<norm2(error) <<std::endl;
std::cout<<GridLogMessage << " summed-combined "<<norm2(error) <<std::endl;
src=1.0;
Chebyshev<LatticeFermion> Cheby(0.1,40.0,200,InverseApproximation);
std::cout<<"Chebuy approx vector "<<std::endl;
std::cout<<GridLogMessage<<"Chebuy approx vector "<<std::endl;
Cheby(Diagonal,src,combined);
std::ofstream of("cheby");
Cheby.csv(of);
@@ -140,10 +140,10 @@ int main (int argc, char ** argv)
Diagonal.ApplyInverse(src,reference);
error = reference - combined;
std::cout << "Chebyshev inverse test "<<std::endl;
std::cout << " Reference "<<norm2(reference)<<std::endl;
std::cout << " combined "<<norm2(combined) <<std::endl;
std::cout << " error "<<norm2(error) <<std::endl;
std::cout<<GridLogMessage << "Chebyshev inverse test "<<std::endl;
std::cout<<GridLogMessage << " Reference "<<norm2(reference)<<std::endl;
std::cout<<GridLogMessage << " combined "<<norm2(combined) <<std::endl;
std::cout<<GridLogMessage << " error "<<norm2(error) <<std::endl;
Grid_finalize();
}

12
tests/Test_nersc_io.cc
View File

@@ -55,7 +55,7 @@ int main (int argc, char ** argv)
#endif
double vol = Fine.gSites();
Complex PlaqScale(1.0/vol/6.0/3.0);
std::cout <<"PlaqScale" << PlaqScale<<std::endl;
std::cout<<GridLogMessage <<"PlaqScale" << PlaqScale<<std::endl;
std::vector<TComplex> Plaq_T(orthosz);
sliceSum(Plaq,Plaq_T,Nd-1);
@@ -66,29 +66,29 @@ int main (int argc, char ** argv)
for(int t=0;t<Nt;t++){
Plaq_T_sum = Plaq_T_sum+Plaq_T[t];
Complex Pt=TensorRemove(Plaq_T[t]);
std::cout << "sliced ["<<t<<"]" <<Pt*PlaqScale*Real(Nt)<<std::endl;
std::cout<<GridLogMessage << "sliced ["<<t<<"]" <<Pt*PlaqScale*Real(Nt)<<std::endl;
}
{
Complex Pt = TensorRemove(Plaq_T_sum);
std::cout << "total " <<Pt*PlaqScale<<std::endl;
std::cout<<GridLogMessage << "total " <<Pt*PlaqScale<<std::endl;
}
TComplex Tp = sum(Plaq);
Complex p = TensorRemove(Tp);
std::cout << "calculated plaquettes " <<p*PlaqScale<<std::endl;
std::cout<<GridLogMessage << "calculated plaquettes " <<p*PlaqScale<<std::endl;
Complex LinkTraceScale(1.0/vol/4.0/3.0);
TComplex Tl = sum(LinkTrace);
Complex l = TensorRemove(Tl);
std::cout << "calculated link trace " <<l*LinkTraceScale<<std::endl;
std::cout<<GridLogMessage << "calculated link trace " <<l*LinkTraceScale<<std::endl;
blockSum(cPlaq,Plaq);
TComplex TcP = sum(cPlaq);
Complex ll= TensorRemove(TcP);
std::cout << "coarsened plaquettes sum to " <<ll*PlaqScale<<std::endl;
std::cout<<GridLogMessage << "coarsened plaquettes sum to " <<ll*PlaqScale<<std::endl;
Grid_finalize();
}

4
tests/Test_quenched_update.cc
View File

@@ -47,7 +47,7 @@ int main (int argc, char ** argv)
RealD plaq = ColourWilsonLoops::avgPlaquette(Umu);
std::cout<<"sweep "<<sweep<<" PLAQUETTE "<<plaq<<std::endl;
std::cout<<GridLogMessage<<"sweep "<<sweep<<" PLAQUETTE "<<plaq<<std::endl;
for( int cb=0;cb<2;cb++ ) {
@@ -55,7 +55,7 @@ int main (int argc, char ** argv)
mask= zero;
setCheckerboard(mask,one);
// std::cout<<mask<<std::endl;
// std::cout<<GridLogMessage<<mask<<std::endl;
for(int mu=0;mu<Nd;mu++){
// Get Link and Staple term in action; must contain Beta and

14
tests/Test_remez.cc
View File

@@ -8,7 +8,7 @@ int main (int argc, char ** argv)
{
Grid_init(&argc,&argv);
std::cout << "Testing Remez"<<std::endl;
std::cout<<GridLogMessage << "Testing Remez"<<std::endl;
double lo=0.01;
double hi=1.0;
@@ -20,13 +20,13 @@ int main (int argc, char ** argv)
// sqrt and inverse sqrt
////////////////////////////////////////
std::cout << "Generating degree "<<degree<<" for x^(1/2)"<<std::endl;
std::cout<<GridLogMessage << "Generating degree "<<degree<<" for x^(1/2)"<<std::endl;
remez.generateApprox(degree,1,2);
MultiShiftFunction Sqrt(remez,1.0,false);
MultiShiftFunction InvSqrt(remez,1.0,true);
std::cout << "Generating degree "<<degree<<" for x^(1/4)"<<std::endl;
std::cout<<GridLogMessage << "Generating degree "<<degree<<" for x^(1/4)"<<std::endl;
remez.generateApprox(degree,1,4);
MultiShiftFunction SqrtSqrt(remez,1.0,false);
MultiShiftFunction InvSqrtSqrt(remez,1.0,true);
@@ -49,10 +49,10 @@ int main (int argc, char ** argv)
double aisx =InvSqrt.approx(x);
double aissx=InvSqrtSqrt.approx(x);
std::cout << "x^(1/2) : "<<sx<<" "<<asx<<std::endl;
std::cout << "x^(1/4) : "<<ssx<<" "<<assx<<std::endl;
std::cout << "x^(-1/2): "<<isx<<" "<<aisx<<std::endl;
std::cout << "x^(-1/4): "<<issx<<" "<<aissx<<std::endl;
std::cout<<GridLogMessage << "x^(1/2) : "<<sx<<" "<<asx<<std::endl;
std::cout<<GridLogMessage << "x^(1/4) : "<<ssx<<" "<<assx<<std::endl;
std::cout<<GridLogMessage << "x^(-1/2): "<<isx<<" "<<aisx<<std::endl;
std::cout<<GridLogMessage << "x^(-1/4): "<<issx<<" "<<aissx<<std::endl;
assert(fabs(sx-asx)<1.0e-6);
assert(fabs(ssx-assx)<1.0e-6);
assert(fabs(isx-aisx)<1.0e-6);

14
tests/Test_rng.cc
View File

@@ -24,27 +24,27 @@ int main (int argc, char ** argv)
SpinMatrix rnd ;
random(sRNG,rnd);
std::cout<<"Random Spin Matrix (random_device)\n"<< rnd<<std::endl;
std::cout<<GridLogMessage<<"Random Spin Matrix (random_device)\n"<< rnd<<std::endl;
random(fsRNG,rnd);
std::cout<<"Random Spin Matrix (fixed seed)\n"<< rnd<<std::endl;
std::cout<<GridLogMessage<<"Random Spin Matrix (fixed seed)\n"<< rnd<<std::endl;
SpinVector rv;
random(sRNG,rv);
std::cout<<"Random Spin Vector (random device)\n"<< rv<<std::endl;
std::cout<<GridLogMessage<<"Random Spin Vector (random device)\n"<< rv<<std::endl;
random(fsRNG,rv);
std::cout<<"Random Spin Vector (fixed seed)\n"<< rv<<std::endl;
std::cout<<GridLogMessage<<"Random Spin Vector (fixed seed)\n"<< rv<<std::endl;
gaussian(fsRNG,rv);
std::cout<<"Gaussian Spin Vector (fixed seed)\n"<< rv<<std::endl;
std::cout<<GridLogMessage<<"Gaussian Spin Vector (fixed seed)\n"<< rv<<std::endl;
LatticeColourVector lcv(&Grid);
random(pRNG,lcv);
std::cout<<"Random Lattice Colour Vector (random device)\n"<< lcv<<std::endl;
std::cout<<GridLogMessage<<"Random Lattice Colour Vector (random device)\n"<< lcv<<std::endl;
random(fpRNG,lcv);
std::cout<<"Random Lattice Colour Vector (fixed seed)\n"<< lcv<<std::endl;
std::cout<<GridLogMessage<<"Random Lattice Colour Vector (fixed seed)\n"<< lcv<<std::endl;
Grid_finalize();
}

22
tests/Test_rng_fixed.cc
View File

@@ -20,31 +20,31 @@ int main (int argc, char ** argv)
GridParallelRNG fpRNG(&Grid); fpRNG.SeedFixedIntegers(seeds);
vComplexF tmp; random(fsRNG,tmp);
std::cout<<"Random vComplexF (fixed seed)\n"<< tmp<<std::endl;
std::cout<<"conjugate(tmp)\n"<< conjugate(tmp)<<std::endl;
std::cout<<"conjugate(tmp)*tmp\n"<< conjugate(tmp)*tmp<<std::endl;
std::cout<<"innerProduct"<< innerProduct(tmp,tmp)<<std::endl;
std::cout<<"Reduce(innerProduct)"<< Reduce(innerProduct(tmp,tmp))<<std::endl;
std::cout<<GridLogMessage<<"Random vComplexF (fixed seed)\n"<< tmp<<std::endl;
std::cout<<GridLogMessage<<"conjugate(tmp)\n"<< conjugate(tmp)<<std::endl;
std::cout<<GridLogMessage<<"conjugate(tmp)*tmp\n"<< conjugate(tmp)*tmp<<std::endl;
std::cout<<GridLogMessage<<"innerProduct"<< innerProduct(tmp,tmp)<<std::endl;
std::cout<<GridLogMessage<<"Reduce(innerProduct)"<< Reduce(innerProduct(tmp,tmp))<<std::endl;
SpinMatrix rnd ;
random(fsRNG,rnd);
std::cout<<"Random Spin Matrix (fixed seed)\n"<< rnd<<std::endl;
std::cout<<GridLogMessage<<"Random Spin Matrix (fixed seed)\n"<< rnd<<std::endl;
SpinVector rv;
random(fsRNG,rv);
std::cout<<"Random Spin Vector (fixed seed)\n"<< rv<<std::endl;
std::cout<<GridLogMessage<<"Random Spin Vector (fixed seed)\n"<< rv<<std::endl;
gaussian(fsRNG,rv);
std::cout<<"Gaussian Spin Vector (fixed seed)\n"<< rv<<std::endl;
std::cout<<GridLogMessage<<"Gaussian Spin Vector (fixed seed)\n"<< rv<<std::endl;
LatticeColourVector lcv(&Grid);
LatticeFermion src(&Grid); random(fpRNG,src);
std::cout << "src norm : " << norm2(src)<<std::endl;
std::cout << "src " << src<<std::endl;
std::cout<<GridLogMessage << "src norm : " << norm2(src)<<std::endl;
std::cout<<GridLogMessage << "src " << src<<std::endl;
random(fpRNG,lcv);
std::cout<<"Random Lattice Colour Vector (fixed seed)\n"<< lcv<<std::endl;
std::cout<<GridLogMessage<<"Random Lattice Colour Vector (fixed seed)\n"<< lcv<<std::endl;
Grid_finalize();

40
tests/Test_simd.cc
View File

@@ -113,18 +113,18 @@ void Tester(const functor &func)
}
extract<vec,scal>(v_result,result);
std::cout << " " << func.name()<<std::endl;
std::cout<<GridLogMessage << " " << func.name()<<std::endl;
int ok=0;
for(int i=0;i<Nsimd;i++){
if ( abs(reference[i]-result[i])>0){
std::cout<< "*****" << std::endl;
std::cout<< "["<<i<<"] "<< abs(reference[i]-result[i]) << " " <<reference[i]<< " " << result[i]<<std::endl;
std::cout<<GridLogMessage<< "*****" << std::endl;
std::cout<<GridLogMessage<< "["<<i<<"] "<< abs(reference[i]-result[i]) << " " <<reference[i]<< " " << result[i]<<std::endl;
ok++;
}
}
if ( ok==0 ) {
std::cout << " OK!" <<std::endl;
std::cout<<GridLogMessage << " OK!" <<std::endl;
}
assert(ok==0);
}
@@ -164,16 +164,16 @@ void ReductionTester(const functor &func)
reference+=tmp;
}
std::cout << " " << func.name()<<std::endl;
std::cout<<GridLogMessage << " " << func.name()<<std::endl;
int ok=0;
if ( abs(reference-result)/abs(reference) > 1.0e-6 ){ // rounding is possible for reduce order
std::cout<< "*****" << std::endl;
std::cout<< abs(reference-result) << " " <<reference<< " " << result<<std::endl;
std::cout<<GridLogMessage<< "*****" << std::endl;
std::cout<<GridLogMessage<< abs(reference-result) << " " <<reference<< " " << result<<std::endl;
ok++;
}
if ( ok==0 ) {
std::cout << " OK!" <<std::endl;
std::cout<<GridLogMessage << " OK!" <<std::endl;
}
assert(ok==0);
}
@@ -194,9 +194,9 @@ int main (int argc, char ** argv)
// Insist that operations on random scalars gives
// identical results to on vectors.
std::cout << "==================================="<< std::endl;
std::cout << "Testing vComplexF "<<std::endl;
std::cout << "==================================="<< std::endl;
std::cout<<GridLogMessage << "==================================="<< std::endl;
std::cout<<GridLogMessage << "Testing vComplexF "<<std::endl;
std::cout<<GridLogMessage << "==================================="<< std::endl;
Tester<ComplexF,vComplexF>(funcTimesI());
Tester<ComplexF,vComplexF>(funcTimesMinusI());
@@ -208,9 +208,9 @@ int main (int argc, char ** argv)
Tester<ComplexF,vComplexF>(funcInnerProduct());
ReductionTester<ComplexF,ComplexF,vComplexF>(funcReduce());
std::cout << "==================================="<< std::endl;
std::cout << "Testing vComplexD "<<std::endl;
std::cout << "==================================="<< std::endl;
std::cout<<GridLogMessage << "==================================="<< std::endl;
std::cout<<GridLogMessage << "Testing vComplexD "<<std::endl;
std::cout<<GridLogMessage << "==================================="<< std::endl;
Tester<ComplexD,vComplexD>(funcTimesI());
@@ -223,9 +223,9 @@ int main (int argc, char ** argv)
Tester<ComplexD,vComplexD>(funcInnerProduct());
ReductionTester<ComplexD,ComplexD,vComplexD>(funcReduce());
std::cout << "==================================="<< std::endl;
std::cout << "Testing vRealF "<<std::endl;
std::cout << "==================================="<< std::endl;
std::cout<<GridLogMessage << "==================================="<< std::endl;
std::cout<<GridLogMessage << "Testing vRealF "<<std::endl;
std::cout<<GridLogMessage << "==================================="<< std::endl;
Tester<RealF,vRealF>(funcPlus());
@@ -236,9 +236,9 @@ int main (int argc, char ** argv)
Tester<RealF,vRealF>(funcInnerProduct());
ReductionTester<RealF,RealF,vRealF>(funcReduce());
std::cout << "==================================="<< std::endl;
std::cout << "Testing vRealD "<<std::endl;
std::cout << "==================================="<< std::endl;
std::cout<<GridLogMessage << "==================================="<< std::endl;
std::cout<<GridLogMessage << "Testing vRealD "<<std::endl;
std::cout<<GridLogMessage << "==================================="<< std::endl;
Tester<RealD,vRealD>(funcPlus());
Tester<RealD,vRealD>(funcMinus());

6
tests/Test_stencil.cc
View File

@@ -34,7 +34,7 @@ int main (int argc, char ** argv)
for(int dir=0;dir<4;dir++){
for(int disp=0;disp<Fine._fdimensions[dir];disp++){
std::cout << "Using stencil to shift dim "<<dir<< " by "<<disp<<std::endl;
std::cout<<GridLogMessage << "Using stencil to shift dim "<<dir<< " by "<<disp<<std::endl;
// start to test the Cartesian npoint stencil infrastructure
int npoint=1;
std::vector<int> directions(npoint,dir);
@@ -73,7 +73,7 @@ int main (int argc, char ** argv)
Real nrmB = norm2(Bar);
Diff = Check-Bar;
Real nrm = norm2(Diff);
std::cout<<"N2diff ="<<nrm<<" "<<nrmC<<" " <<nrmB<<std::endl;
std::cout<<GridLogMessage<<"N2diff ="<<nrm<<" "<<nrmC<<" " <<nrmB<<std::endl;
Real snrmC =0;
Real snrmB =0;
@@ -110,7 +110,7 @@ int main (int argc, char ** argv)
}}}}
std::cout<<"scalar N2diff = "<<snrm<<" " <<snrmC<<" "<<snrmB<<std::endl;
std::cout<<GridLogMessage<<"scalar N2diff = "<<snrm<<" " <<snrmC<<" "<<snrmB<<std::endl;
}

64
tests/Test_wilson_even_odd.cc
View File

@@ -27,7 +27,7 @@ int main (int argc, char ** argv)
GridRedBlackCartesian RBGrid(latt_size,simd_layout,mpi_layout);
int threads = GridThread::GetThreads();
std::cout << "Grid is setup to use "<<threads<<" threads"<<std::endl;
std::cout<<GridLogMessage << "Grid is setup to use "<<threads<<" threads"<<std::endl;
std::vector<int> seeds({1,2,3,4});
@@ -55,9 +55,9 @@ int main (int argc, char ** argv)
Umu=zero;
for(int nn=0;nn<Nd;nn++){
random(pRNG,U[nn]);
std::cout<<"U[nn]"<<norm2(U[nn])<<std::endl;
std::cout<<GridLogMessage<<"U[nn]"<<norm2(U[nn])<<std::endl;
PokeIndex<LorentzIndex>(Umu,U[nn],nn);
std::cout<<"Umu"<<norm2(Umu)<<std::endl;
std::cout<<GridLogMessage<<"Umu"<<norm2(Umu)<<std::endl;
}
RealD mass=0.1;
@@ -70,30 +70,30 @@ int main (int argc, char ** argv)
LatticeFermion r_o (&RBGrid);
LatticeFermion r_eo (&Grid);
std::cout<<"=========================================================="<<std::endl;
std::cout<<"= Testing that Deo + Doe = Dunprec "<<std::endl;
std::cout<<"=========================================================="<<std::endl;
std::cout<<GridLogMessage<<"=========================================================="<<std::endl;
std::cout<<GridLogMessage<<"= Testing that Deo + Doe = Dunprec "<<std::endl;
std::cout<<GridLogMessage<<"=========================================================="<<std::endl;
pickCheckerboard(Even,src_e,src);
pickCheckerboard(Odd,src_o,src);
Dw.Meooe(src_e,r_o); std::cout<<"Applied Meo"<<std::endl;
Dw.Meooe(src_o,r_e); std::cout<<"Applied Moe"<<std::endl;
Dw.Meooe(src_e,r_o); std::cout<<GridLogMessage<<"Applied Meo"<<std::endl;
Dw.Meooe(src_o,r_e); std::cout<<GridLogMessage<<"Applied Moe"<<std::endl;
Dw.Dhop (src,ref,DaggerNo);
setCheckerboard(r_eo,r_o);
setCheckerboard(r_eo,r_e);
err= ref - r_eo;
std::cout << "EO norm diff "<< norm2(err)<< " "<<norm2(ref)<< " " << norm2(r_eo) <<std::endl;
std::cout<<GridLogMessage << "EO norm diff "<< norm2(err)<< " "<<norm2(ref)<< " " << norm2(r_eo) <<std::endl;
LatticeComplex cerr(&Grid);
cerr = localInnerProduct(err,err);
std::cout<<"=============================================================="<<std::endl;
std::cout<<"= Test Ddagger is the dagger of D by requiring "<<std::endl;
std::cout<<"= < phi | Deo | chi > * = < chi | Deo^dag| phi> "<<std::endl;
std::cout<<"=============================================================="<<std::endl;
std::cout<<GridLogMessage<<"=============================================================="<<std::endl;
std::cout<<GridLogMessage<<"= Test Ddagger is the dagger of D by requiring "<<std::endl;
std::cout<<GridLogMessage<<"= < phi | Deo | chi > * = < chi | Deo^dag| phi> "<<std::endl;
std::cout<<GridLogMessage<<"=============================================================="<<std::endl;
LatticeFermion chi_e (&RBGrid);
LatticeFermion chi_o (&RBGrid);
@@ -123,15 +123,15 @@ int main (int argc, char ** argv)
ComplexD cDpe = innerProduct(chi_e,dphi_e);
ComplexD cDpo = innerProduct(chi_o,dphi_o);
std::cout <<"e "<<pDce<<" "<<cDpe <<std::endl;
std::cout <<"o "<<pDco<<" "<<cDpo <<std::endl;
std::cout<<GridLogMessage <<"e "<<pDce<<" "<<cDpe <<std::endl;
std::cout<<GridLogMessage <<"o "<<pDco<<" "<<cDpo <<std::endl;
std::cout <<"pDce - conj(cDpo) "<< pDce-conj(cDpo) <<std::endl;
std::cout <<"pDco - conj(cDpe) "<< pDco-conj(cDpe) <<std::endl;
std::cout<<GridLogMessage <<"pDce - conj(cDpo) "<< pDce-conj(cDpo) <<std::endl;
std::cout<<GridLogMessage <<"pDco - conj(cDpe) "<< pDco-conj(cDpe) <<std::endl;
std::cout<<"=============================================================="<<std::endl;
std::cout<<"= Test MeeInv Mee = 1 "<<std::endl;
std::cout<<"=============================================================="<<std::endl;
std::cout<<GridLogMessage<<"=============================================================="<<std::endl;
std::cout<<GridLogMessage<<"= Test MeeInv Mee = 1 "<<std::endl;
std::cout<<GridLogMessage<<"=============================================================="<<std::endl;
pickCheckerboard(Even,chi_e,chi);
pickCheckerboard(Odd ,chi_o,chi);
@@ -146,11 +146,11 @@ int main (int argc, char ** argv)
setCheckerboard(phi,phi_o);
err = phi-chi;
std::cout << "norm diff "<< norm2(err)<< std::endl;
std::cout<<GridLogMessage << "norm diff "<< norm2(err)<< std::endl;
std::cout<<"=============================================================="<<std::endl;
std::cout<<"= Test MeeInvDag MeeDag = 1 "<<std::endl;
std::cout<<"=============================================================="<<std::endl;
std::cout<<GridLogMessage<<"=============================================================="<<std::endl;
std::cout<<GridLogMessage<<"= Test MeeInvDag MeeDag = 1 "<<std::endl;
std::cout<<GridLogMessage<<"=============================================================="<<std::endl;
pickCheckerboard(Even,chi_e,chi);
pickCheckerboard(Odd ,chi_o,chi);
@@ -165,11 +165,11 @@ int main (int argc, char ** argv)
setCheckerboard(phi,phi_o);
err = phi-chi;
std::cout << "norm diff "<< norm2(err)<< std::endl;
std::cout<<GridLogMessage << "norm diff "<< norm2(err)<< std::endl;
std::cout<<"=============================================================="<<std::endl;
std::cout<<"= Test MpcDagMpc is Hermitian "<<std::endl;
std::cout<<"=============================================================="<<std::endl;
std::cout<<GridLogMessage<<"=============================================================="<<std::endl;
std::cout<<GridLogMessage<<"= Test MpcDagMpc is Hermitian "<<std::endl;
std::cout<<GridLogMessage<<"=============================================================="<<std::endl;
random(pRNG,phi);
random(pRNG,chi);
@@ -191,11 +191,11 @@ int main (int argc, char ** argv)
cDpe = innerProduct(chi_e,dphi_e);
cDpo = innerProduct(chi_o,dphi_o);
std::cout <<"e "<<pDce<<" "<<cDpe <<std::endl;
std::cout <<"o "<<pDco<<" "<<cDpo <<std::endl;
std::cout<<GridLogMessage <<"e "<<pDce<<" "<<cDpe <<std::endl;
std::cout<<GridLogMessage <<"o "<<pDco<<" "<<cDpo <<std::endl;
std::cout <<"pDce - conj(cDpo) "<< pDco-conj(cDpo) <<std::endl;
std::cout <<"pDco - conj(cDpe) "<< pDce-conj(cDpe) <<std::endl;
std::cout<<GridLogMessage <<"pDce - conj(cDpo) "<< pDco-conj(cDpo) <<std::endl;
std::cout<<GridLogMessage <<"pDco - conj(cDpe) "<< pDce-conj(cDpe) <<std::endl;
Grid_finalize();
}