portelli/Grid
1
0
Fork 0
mirror of https://github.com/paboyle/Grid.git synced 2025-08-05 14:07:12 +01:00
Code Releases Activity

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

Browse Source
This commit is contained in:
Chulwoo Jung
2021-09-03 17:38:10 -04:00
parent 43943008bf b06526bc1e
commit 23b9c6b5f5
280 changed files with 12382 additions and 6738 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
tests/IO/Test_ildg_io.cc
View File

@@ -69,7 +69,7 @@ int main (int argc, char ** argv)
std::vector<LatticeColourMatrix> U(4,&Fine);
SU3::HotConfiguration(pRNGa,Umu);
SU<Nc>::HotConfiguration(pRNGa,Umu);
FieldMetaData header;

2
tests/IO/Test_nersc_io.cc
View File

@@ -84,7 +84,7 @@ int main (int argc, char ** argv)
std::vector<LatticeColourMatrix> U(4,&Fine);
SU3::HotConfiguration(pRNGa,Umu);
SU<Nc>::HotConfiguration(pRNGa,Umu);
FieldMetaData header;
std::string file("./ckpoint_lat.4000");

2
tests/Test_cayley_even_odd_vec.cc
View File

@@ -80,7 +80,7 @@ int main (int argc, char ** argv)
GridParallelRNG sRNG5(sFGrid); sRNG5.SeedFixedIntegers(seeds5);
LatticeGaugeField Umu(UGrid);
SU3::HotConfiguration(RNG4,Umu);
SU<Nc>::HotConfiguration(RNG4,Umu);
RealD mass=0.1;
RealD M5 =1.8;

2
tests/Test_compressed_lanczos_hot_start.cc
View File

@@ -202,7 +202,7 @@ int main (int argc, char ** argv) {
std::vector<int> seeds4({1,2,3,4});
GridParallelRNG RNG4(UGrid); RNG4.SeedFixedIntegers(seeds4);
LatticeGaugeField Umu(UGrid);
SU3::HotConfiguration(RNG4,Umu);
SU<Nc>::HotConfiguration(RNG4,Umu);
// FieldMetaData header;
// NerscIO::readConfiguration(Umu,header,Params.config);

2
tests/Test_dwf_mixedcg_prec.cc
View File

@@ -71,7 +71,7 @@ int main (int argc, char ** argv)
LatticeGaugeFieldD Umu(UGrid);
LatticeGaugeFieldF Umu_f(UGrid_f);
SU3::HotConfiguration(RNG4,Umu);
SU<Nc>::HotConfiguration(RNG4,Umu);
precisionChange(Umu_f,Umu);

2
tests/Test_dwf_mixedcg_prec_halfcomms.cc
View File

@@ -69,7 +69,7 @@ int main (int argc, char ** argv)
LatticeGaugeFieldD Umu(UGrid);
LatticeGaugeFieldF Umu_f(UGrid_f);
SU3::HotConfiguration(RNG4,Umu);
SU<Nc>::HotConfiguration(RNG4,Umu);
precisionChange(Umu_f,Umu);

2
tests/core/Test_cf_coarsen_support.cc
View File

@@ -64,7 +64,7 @@ int main (int argc, char ** argv)
LatticeFermion ref(FGrid); ref=Zero();
LatticeFermion tmp(FGrid);
LatticeFermion err(FGrid);
LatticeGaugeField Umu(UGrid); SU3::HotConfiguration(RNG4,Umu);
LatticeGaugeField Umu(UGrid); SU<Nc>::HotConfiguration(RNG4,Umu);
std::vector<LatticeColourMatrix> U(4,UGrid);
for(int mu=0;mu<Nd;mu++){

2
tests/core/Test_checker.cc
View File

@@ -131,7 +131,7 @@ int main (int argc, char ** argv)
// LatticeFermion result(FGrid); result=Zero();
// LatticeGaugeField Umu(UGrid);
// SU3::HotConfiguration(RNG4,Umu);
// SU<Nc>::HotConfiguration(RNG4,Umu);
// std::vector<LatticeColourMatrix> U(4,UGrid);
// for(int mu=0;mu<Nd;mu++){

2
tests/core/Test_contfrac_even_odd.cc
View File

@@ -69,7 +69,7 @@ int main (int argc, char ** argv)
GridParallelRNG RNG5(FGrid); RNG5.SeedFixedIntegers(seeds5);
GridParallelRNG RNG4(UGrid); RNG4.SeedFixedIntegers(seeds4);
LatticeGaugeField Umu(UGrid); SU3::HotConfiguration(RNG4,Umu);
LatticeGaugeField Umu(UGrid); SU<Nc>::HotConfiguration(RNG4,Umu);
std::vector<LatticeColourMatrix> U(4,UGrid);
RealD mass=0.1;

2
tests/core/Test_dwf_eofa_even_odd.cc
View File

@@ -73,7 +73,7 @@ int main (int argc, char ** argv)
LatticeFermion ref (FGrid); ref = Zero();
LatticeFermion tmp (FGrid); tmp = Zero();
LatticeFermion err (FGrid); err = Zero();
LatticeGaugeField Umu (UGrid); SU3::HotConfiguration(RNG4, Umu);
LatticeGaugeField Umu (UGrid); SU<Nc>::HotConfiguration(RNG4, Umu);
std::vector<LatticeColourMatrix> U(4,UGrid);
// Only one non-zero (y)

2
tests/core/Test_dwf_even_odd.cc
View File

@@ -72,7 +72,7 @@ int main (int argc, char ** argv)
LatticeFermion ref(FGrid); ref=Zero();
LatticeFermion tmp(FGrid); tmp=Zero();
LatticeFermion err(FGrid); tmp=Zero();
LatticeGaugeField Umu(UGrid); SU3::HotConfiguration(RNG4,Umu);
LatticeGaugeField Umu(UGrid); SU<Nc>::HotConfiguration(RNG4,Umu);
std::vector<LatticeColourMatrix> U(4,UGrid);
// Only one non-zero (y)

2
tests/core/Test_fft.cc
View File

@@ -138,7 +138,7 @@ int main (int argc, char ** argv)
LatticeGaugeFieldD Umu(&GRID);
SU3::ColdConfiguration(pRNG,Umu); // Unit gauge
SU<Nc>::ColdConfiguration(pRNG,Umu); // Unit gauge
// Umu=Zero();
////////////////////////////////////////////////////
// Wilson test

8
tests/core/Test_fft_gfix.cc
View File

@@ -73,11 +73,11 @@ int main (int argc, char ** argv)
LatticeColourMatrix xform2(&GRID); // Gauge xform
LatticeColourMatrix xform3(&GRID); // Gauge xform
SU3::ColdConfiguration(pRNG,Umu); // Unit gauge
SU<Nc>::ColdConfiguration(pRNG,Umu); // Unit gauge
Uorg=Umu;
Urnd=Umu;
SU3::RandomGaugeTransform(pRNG,Urnd,g); // Unit gauge
SU<Nc>::RandomGaugeTransform(pRNG,Urnd,g); // Unit gauge
Real plaq=WilsonLoops<PeriodicGimplR>::avgPlaquette(Umu);
std::cout << " Initial plaquette "<<plaq << std::endl;
@@ -121,7 +121,7 @@ int main (int argc, char ** argv)
std::cout<< "* Testing non-unit configuration *" <<std::endl;
std::cout<< "*****************************************************************" <<std::endl;
SU3::HotConfiguration(pRNG,Umu); // Unit gauge
SU<Nc>::HotConfiguration(pRNG,Umu); // Unit gauge
plaq=WilsonLoops<PeriodicGimplR>::avgPlaquette(Umu);
std::cout << " Initial plaquette "<<plaq << std::endl;
@@ -136,7 +136,7 @@ int main (int argc, char ** argv)
std::cout<< "*****************************************************************" <<std::endl;
Umu=Urnd;
SU3::HotConfiguration(pRNG,Umu); // Unit gauge
SU<Nc>::HotConfiguration(pRNG,Umu); // Unit gauge
plaq=WilsonLoops<PeriodicGimplR>::avgPlaquette(Umu);
std::cout << " Initial plaquette "<<plaq << std::endl;

2
tests/core/Test_gparity.cc
View File

@@ -114,7 +114,7 @@ int main (int argc, char ** argv)
GridParallelRNG RNG4_2f(UGrid_2f); RNG4_2f.SeedFixedIntegers(seeds4);
GparityGaugeField Umu_2f(UGrid_2f);
SU3::HotConfiguration(RNG4_2f,Umu_2f);
SU<Nc>::HotConfiguration(RNG4_2f,Umu_2f);
StandardFermionField src (FGrid_2f);
StandardFermionField tmpsrc(FGrid_2f);

2
tests/core/Test_gpwilson_even_odd.cc
View File

@@ -61,7 +61,7 @@ int main (int argc, char ** argv)
FermionField ref(&Grid); ref=Zero();
FermionField tmp(&Grid); tmp=Zero();
FermionField err(&Grid); tmp=Zero();
LatticeGaugeField Umu(&Grid); SU3::HotConfiguration(pRNG,Umu);
LatticeGaugeField Umu(&Grid); SU<Nc>::HotConfiguration(pRNG,Umu);
std::vector<LatticeColourMatrix> U(4,&Grid);
double volume=1;

70
tests/core/Test_lie_generators.cc
View File

@@ -66,7 +66,7 @@ int main(int argc, char** argv) {
std::cout << GridLogMessage << "*********************************************"
<< std::endl;
std::cout << GridLogMessage << "* Generators for SU(3)" << std::endl;
std::cout << GridLogMessage << "* Generators for SU(Nc" << std::endl;
std::cout << GridLogMessage << "*********************************************"
<< std::endl;
SU3::printGenerators();
@@ -114,8 +114,8 @@ int main(int argc, char** argv) {
LatticeGaugeField U(grid), V(grid);
SU<Nc>::HotConfiguration<LatticeGaugeField>(gridRNG, U);
SU<Nc>::HotConfiguration<LatticeGaugeField>(gridRNG, V);
SU3::HotConfiguration<LatticeGaugeField>(gridRNG, U);
SU3::HotConfiguration<LatticeGaugeField>(gridRNG, V);
// Adjoint representation
// Test group structure
@@ -123,8 +123,8 @@ int main(int argc, char** argv) {
LatticeGaugeField UV(grid);
UV = Zero();
for (int mu = 0; mu < Nd; mu++) {
SU<Nc>::LatticeMatrix Umu = peekLorentz(U,mu);
SU<Nc>::LatticeMatrix Vmu = peekLorentz(V,mu);
SU3::LatticeMatrix Umu = peekLorentz(U,mu);
SU3::LatticeMatrix Vmu = peekLorentz(V,mu);
pokeLorentz(UV,Umu*Vmu, mu);
}
@@ -151,16 +151,16 @@ int main(int argc, char** argv) {
// Check correspondence of algebra and group transformations
// Create a random vector
SU<Nc>::LatticeAlgebraVector h_adj(grid);
SU3::LatticeAlgebraVector h_adj(grid);
typename AdjointRep<Nc>::LatticeMatrix Ar(grid);
random(gridRNG,h_adj);
h_adj = real(h_adj);
SU_Adjoint<Nc>::AdjointLieAlgebraMatrix(h_adj,Ar);
// Re-extract h_adj
SU<Nc>::LatticeAlgebraVector h_adj2(grid);
SU3::LatticeAlgebraVector h_adj2(grid);
SU_Adjoint<Nc>::projectOnAlgebra(h_adj2, Ar);
SU<Nc>::LatticeAlgebraVector h_diff = h_adj - h_adj2;
SU3::LatticeAlgebraVector h_diff = h_adj - h_adj2;
std::cout << GridLogMessage << "Projections structure check vector difference (Adjoint representation) : " << norm2(h_diff) << std::endl;
// Exponentiate
@@ -183,14 +183,14 @@ int main(int argc, char** argv) {
// Construct the fundamental matrix in the group
SU<Nc>::LatticeMatrix Af(grid);
SU<Nc>::FundamentalLieAlgebraMatrix(h_adj,Af);
SU<Nc>::LatticeMatrix Ufund(grid);
SU3::LatticeMatrix Af(grid);
SU3::FundamentalLieAlgebraMatrix(h_adj,Af);
SU3::LatticeMatrix Ufund(grid);
Ufund = expMat(Af, 1.0, 16);
// Check unitarity
SU<Nc>::LatticeMatrix uno_f(grid);
SU3::LatticeMatrix uno_f(grid);
uno_f = 1.0;
SU<Nc>::LatticeMatrix UnitCheck(grid);
SU3::LatticeMatrix UnitCheck(grid);
UnitCheck = Ufund * adj(Ufund) - uno_f;
std::cout << GridLogMessage << "unitarity check 1: " << norm2(UnitCheck)
<< std::endl;
@@ -311,14 +311,14 @@ int main(int argc, char** argv) {
// Test group structure
// (U_f * V_f)_r = U_r * V_r
LatticeGaugeField U2(grid), V2(grid);
SU<Nc>::HotConfiguration<LatticeGaugeField>(gridRNG, U2);
SU<Nc>::HotConfiguration<LatticeGaugeField>(gridRNG, V2);
SU3::HotConfiguration<LatticeGaugeField>(gridRNG, U2);
SU3::HotConfiguration<LatticeGaugeField>(gridRNG, V2);
LatticeGaugeField UV2(grid);
UV2 = Zero();
for (int mu = 0; mu < Nd; mu++) {
SU<Nc>::LatticeMatrix Umu2 = peekLorentz(U2,mu);
SU<Nc>::LatticeMatrix Vmu2 = peekLorentz(V2,mu);
SU3::LatticeMatrix Umu2 = peekLorentz(U2,mu);
SU3::LatticeMatrix Vmu2 = peekLorentz(V2,mu);
pokeLorentz(UV2,Umu2*Vmu2, mu);
}
@@ -345,16 +345,16 @@ int main(int argc, char** argv) {
// Check correspondence of algebra and group transformations
// Create a random vector
SU<Nc>::LatticeAlgebraVector h_sym(grid);
SU3::LatticeAlgebraVector h_sym(grid);
typename TwoIndexRep< Nc, Symmetric>::LatticeMatrix Ar_sym(grid);
random(gridRNG,h_sym);
h_sym = real(h_sym);
SU_TwoIndex<Nc,Symmetric>::TwoIndexLieAlgebraMatrix(h_sym,Ar_sym);
// Re-extract h_sym
SU<Nc>::LatticeAlgebraVector h_sym2(grid);
SU3::LatticeAlgebraVector h_sym2(grid);
SU_TwoIndex< Nc, Symmetric>::projectOnAlgebra(h_sym2, Ar_sym);
SU<Nc>::LatticeAlgebraVector h_diff_sym = h_sym - h_sym2;
SU3::LatticeAlgebraVector h_diff_sym = h_sym - h_sym2;
std::cout << GridLogMessage << "Projections structure check vector difference (Two Index Symmetric): " << norm2(h_diff_sym) << std::endl;
@@ -379,11 +379,11 @@ int main(int argc, char** argv) {
// Construct the fundamental matrix in the group
SU<Nc>::LatticeMatrix Af_sym(grid);
SU<Nc>::FundamentalLieAlgebraMatrix(h_sym,Af_sym);
SU<Nc>::LatticeMatrix Ufund2(grid);
SU3::LatticeMatrix Af_sym(grid);
SU3::FundamentalLieAlgebraMatrix(h_sym,Af_sym);
SU3::LatticeMatrix Ufund2(grid);
Ufund2 = expMat(Af_sym, 1.0, 16);
SU<Nc>::LatticeMatrix UnitCheck2(grid);
SU3::LatticeMatrix UnitCheck2(grid);
UnitCheck2 = Ufund2 * adj(Ufund2) - uno_f;
std::cout << GridLogMessage << "unitarity check 1: " << norm2(UnitCheck2)
<< std::endl;
@@ -421,14 +421,14 @@ int main(int argc, char** argv) {
// Test group structure
// (U_f * V_f)_r = U_r * V_r
LatticeGaugeField U2A(grid), V2A(grid);
SU<Nc>::HotConfiguration<LatticeGaugeField>(gridRNG, U2A);
SU<Nc>::HotConfiguration<LatticeGaugeField>(gridRNG, V2A);
SU3::HotConfiguration<LatticeGaugeField>(gridRNG, U2A);
SU3::HotConfiguration<LatticeGaugeField>(gridRNG, V2A);
LatticeGaugeField UV2A(grid);
UV2A = Zero();
for (int mu = 0; mu < Nd; mu++) {
SU<Nc>::LatticeMatrix Umu2A = peekLorentz(U2,mu);
SU<Nc>::LatticeMatrix Vmu2A = peekLorentz(V2,mu);
SU3::LatticeMatrix Umu2A = peekLorentz(U2,mu);
SU3::LatticeMatrix Vmu2A = peekLorentz(V2,mu);
pokeLorentz(UV2A,Umu2A*Vmu2A, mu);
}
@@ -455,16 +455,16 @@ int main(int argc, char** argv) {
// Check correspondence of algebra and group transformations
// Create a random vector
SU<Nc>::LatticeAlgebraVector h_Asym(grid);
SU3::LatticeAlgebraVector h_Asym(grid);
typename TwoIndexRep< Nc, AntiSymmetric>::LatticeMatrix Ar_Asym(grid);
random(gridRNG,h_Asym);
h_Asym = real(h_Asym);
SU_TwoIndex< Nc, AntiSymmetric>::TwoIndexLieAlgebraMatrix(h_Asym,Ar_Asym);
// Re-extract h_sym
SU<Nc>::LatticeAlgebraVector h_Asym2(grid);
SU3::LatticeAlgebraVector h_Asym2(grid);
SU_TwoIndex< Nc, AntiSymmetric>::projectOnAlgebra(h_Asym2, Ar_Asym);
SU<Nc>::LatticeAlgebraVector h_diff_Asym = h_Asym - h_Asym2;
SU3::LatticeAlgebraVector h_diff_Asym = h_Asym - h_Asym2;
std::cout << GridLogMessage << "Projections structure check vector difference (Two Index anti-Symmetric): " << norm2(h_diff_Asym) << std::endl;
@@ -489,11 +489,11 @@ int main(int argc, char** argv) {
// Construct the fundamental matrix in the group
SU<Nc>::LatticeMatrix Af_Asym(grid);
SU<Nc>::FundamentalLieAlgebraMatrix(h_Asym,Af_Asym);
SU<Nc>::LatticeMatrix Ufund2A(grid);
SU3::LatticeMatrix Af_Asym(grid);
SU3::FundamentalLieAlgebraMatrix(h_Asym,Af_Asym);
SU3::LatticeMatrix Ufund2A(grid);
Ufund2A = expMat(Af_Asym, 1.0, 16);
SU<Nc>::LatticeMatrix UnitCheck2A(grid);
SU3::LatticeMatrix UnitCheck2A(grid);
UnitCheck2A = Ufund2A * adj(Ufund2A) - uno_f;
std::cout << GridLogMessage << "unitarity check 1: " << norm2(UnitCheck2A)
<< std::endl;

19
tests/core/Test_main.cc
View File

@@ -231,6 +231,20 @@ int main(int argc, char **argv) {
scalar = localInnerProduct(cVec, cVec);
scalar = localNorm2(cVec);
std::cout << "Testing maxLocalNorm2" <<std::endl;
LatticeComplex rand_scalar(&Fine);
random(FineRNG, rand_scalar); //uniform [0,1]
for(Integer gsite=0;gsite<Fine.gSites();gsite++){ //check on every site independently
scalar = rand_scalar;
TComplex big(10.0);
Coordinate coor;
Fine.GlobalIndexToGlobalCoor(gsite,coor);
pokeSite(big,scalar,coor);
RealD Linfty = maxLocalNorm2(scalar);
assert(Linfty == 100.0);
}
// -=,+=,*=,()
// add,+,sub,-,mult,mac,*
// adj,conjugate
@@ -444,7 +458,7 @@ int main(int argc, char **argv) {
// Lattice 12x12 GEMM
scFooBar = scFoo * scBar;
// Benchmark some simple operations LatticeSU3 * Lattice SU3.
// Benchmark some simple operations LatticeSU<Nc> * Lattice SU<Nc>.
double t0, t1, flops;
double bytes;
int ncall = 5000;
@@ -549,7 +563,8 @@ int main(int argc, char **argv) {
std::vector<int> shiftcoor = coor;
shiftcoor[dir] = (shiftcoor[dir] + shift + latt_size[dir]) %
(latt_size[dir] / mpi_layout[dir]);
(latt_size[dir]);
// (latt_size[dir] / mpi_layout[dir]);
std::vector<int> rl(4);
for (int dd = 0; dd < 4; dd++) {

2
tests/core/Test_mobius_eofa_even_odd.cc
View File

@@ -73,7 +73,7 @@ int main (int argc, char ** argv)
LatticeFermion ref (FGrid); ref = Zero();
LatticeFermion tmp (FGrid); tmp = Zero();
LatticeFermion err (FGrid); err = Zero();
LatticeGaugeField Umu (UGrid); SU3::HotConfiguration(RNG4, Umu);
LatticeGaugeField Umu (UGrid); SU<Nc>::HotConfiguration(RNG4, Umu);
std::vector<LatticeColourMatrix> U(4,UGrid);
// Only one non-zero (y)

6
tests/core/Test_quenched_update.cc
View File

@@ -55,7 +55,7 @@ int main (int argc, char ** argv)
GridParallelRNG pRNG(grid); pRNG.SeedFixedIntegers(pseeds);
GridSerialRNG sRNG; sRNG.SeedFixedIntegers(sseeds);
// SU3 colour operatoions
// SU<Nc> colour operatoions
LatticeColourMatrix link(grid);
LatticeColourMatrix staple(grid);
@@ -87,10 +87,10 @@ int main (int argc, char ** argv)
link = PeekIndex<LorentzIndex>(Umu,mu);
for( int subgroup=0;subgroup<SU3::su2subgroups();subgroup++ ) {
for( int subgroup=0;subgroup<SU<Nc>::su2subgroups();subgroup++ ) {
// update Even checkerboard
SU3::SubGroupHeatBath(sRNG,pRNG,beta,link,staple,subgroup,20,mask);
SU<Nc>::SubGroupHeatBath(sRNG,pRNG,beta,link,staple,subgroup,20,mask);
}

145
tests/core/Test_reunitarise.cc Normal file
View File

@@ -0,0 +1,145 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: ./tests/Test_quenched_update.cc
Copyright (C) 2015
Author: Azusa Yamaguchi <ayamaguc@staffmail.ed.ac.uk>
Author: Peter Boyle <paboyle@ph.ed.ac.uk>
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License along
with this program; if not, write to the Free Software Foundation, Inc.,
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
See the full license in the file "LICENSE" in the top level distribution directory
*************************************************************************************/
/* END LEGAL */
#include <Grid/Grid.h>
using namespace std;
using namespace Grid;
;
int main (int argc, char ** argv)
{
Grid_init(&argc,&argv);
std::vector<int> latt({8,8,8,8});
GridCartesian * grid = SpaceTimeGrid::makeFourDimGrid(latt,
GridDefaultSimd(Nd,vComplexD::Nsimd()),
GridDefaultMpi());
GridCartesian * gridF = SpaceTimeGrid::makeFourDimGrid(latt,
GridDefaultSimd(Nd,vComplexF::Nsimd()),
GridDefaultMpi());
///////////////////////////////
// Configuration of known size
///////////////////////////////
LatticeColourMatrixD ident(grid);
LatticeColourMatrixD U(grid);
LatticeColourMatrixD UU(grid);
LatticeColourMatrixD tmp(grid);
LatticeColourMatrixD org(grid);
LatticeColourMatrixF UF(gridF);
LatticeGaugeField Umu(grid);
ident =1.0;
// RNG set up for test
std::vector<int> pseeds({1,2,3,4,5}); // once I caught a fish alive
std::vector<int> sseeds({6,7,8,9,10});// then i let it go again
GridParallelRNG pRNG(grid); pRNG.SeedFixedIntegers(pseeds);
GridSerialRNG sRNG; sRNG.SeedFixedIntegers(sseeds);
SU<Nc>::HotConfiguration(pRNG,Umu);
U = PeekIndex<LorentzIndex>(Umu,0);
org=U;
tmp= U*adj(U) - ident ;
RealD Def1 = norm2( tmp );
std::cout << " Defect1 "<<Def1<<std::endl;
tmp = U - org;
std::cout << "Diff1 "<<norm2(tmp)<<std::endl;
precisionChange(UF,U);
precisionChange(U,UF);
tmp= U*adj(U) - ident ;
RealD Def2 = norm2( tmp );
std::cout << " Defect2 "<<Def2<<std::endl;
tmp = U - org;
std::cout << "Diff2 "<<norm2(tmp)<<std::endl;
U = ProjectOnGroup(U);
tmp= U*adj(U) - ident ;
RealD Def3 = norm2( tmp);
std::cout << " Defect3 "<<Def3<<std::endl;
tmp = U - org;
std::cout << "Diff3 "<<norm2(tmp)<<std::endl;
LatticeComplexD detU(grid);
LatticeComplexD detUU(grid);
detU= Determinant(U) ;
detU=detU-1.0;
std::cout << "Determinant defect before screw up " << norm2(detU)<<std::endl;
std::cout << " Screwing up determinant " << std::endl;
RealD theta = 0.2;
ComplexD phase(cos(theta),sin(theta));
for(int i=0;i<Nc;i++){
auto element = PeekIndex<ColourIndex>(U,Nc-1,i);
element = element * phase;
PokeIndex<ColourIndex>(U,element,Nc-1,i);
}
U=U*0.1;
UU=U;
detU= Determinant(U) ;
detU=detU-1.0;
std::cout << "Determinant defect before projection " <<norm2(detU)<<std::endl;
tmp = U*adj(U) - ident;
std::cout << "Unitarity check before projection " << norm2(tmp)<<std::endl;
ProjectSU3(U);
detU= Determinant(U) ;
detU= detU -1.0;
std::cout << "Determinant ProjectSU3 defect " <<norm2(detU)<<std::endl;
tmp = U*adj(U) - ident;
std::cout << "Unitarity check after projection " << norm2(tmp)<<std::endl;
ProjectSUn(UU);
detUU= Determinant(UU);
detUU= detUU -1.0;
std::cout << "Determinant ProjectSUn defect " <<norm2(detUU)<<std::endl;
tmp = UU*adj(UU) - ident;
std::cout << "Unitarity check after projection " << norm2(tmp)<<std::endl;
Grid_finalize();
}

2
tests/core/Test_staggered.cc
View File

@@ -64,7 +64,7 @@ int main (int argc, char ** argv)
FermionField err(&Grid); tmp=Zero();
FermionField phi (&Grid); random(pRNG,phi);
FermionField chi (&Grid); random(pRNG,chi);
LatticeGaugeField Umu(&Grid); SU3::HotConfiguration(pRNG,Umu);
LatticeGaugeField Umu(&Grid); SU<Nc>::HotConfiguration(pRNG,Umu);
std::vector<LatticeColourMatrix> U(4,&Grid);

2
tests/core/Test_staggered5D.cc
View File

@@ -75,7 +75,7 @@ int main (int argc, char ** argv)
FermionField phi (FGrid); random(pRNG5,phi);
FermionField chi (FGrid); random(pRNG5,chi);
LatticeGaugeField Umu(UGrid); SU3::ColdConfiguration(pRNG4,Umu);
LatticeGaugeField Umu(UGrid); SU<Nc>::ColdConfiguration(pRNG4,Umu);
LatticeGaugeField Umua(UGrid); Umua=Umu;
double volume=Ls;

2
tests/core/Test_staggered5Dvec.cc
View File

@@ -84,7 +84,7 @@ int main (int argc, char ** argv)
FermionField chi (FGrid); random(pRNG5,chi);
LatticeGaugeField Umu(UGrid);
SU3::HotConfiguration(pRNG4,Umu);
SU<Nc>::HotConfiguration(pRNG4,Umu);
/*
for(int mu=1;mu<4;mu++){

2
tests/core/Test_staggered5DvecF.cc
View File

@@ -83,7 +83,7 @@ int main (int argc, char ** argv)
FermionField chi (FGrid); random(pRNG5,chi);
LatticeGaugeFieldF Umu(UGrid);
SU3::HotConfiguration(pRNG4,Umu);
SU<Nc>::HotConfiguration(pRNG4,Umu);
/*
for(int mu=1;mu<4;mu++){

2
tests/core/Test_staggered_naive.cc
View File

@@ -64,7 +64,7 @@ int main (int argc, char ** argv)
FermionField err(&Grid); tmp=Zero();
FermionField phi (&Grid); random(pRNG,phi);
FermionField chi (&Grid); random(pRNG,chi);
LatticeGaugeField Umu(&Grid); SU3::HotConfiguration(pRNG,Umu);
LatticeGaugeField Umu(&Grid); SU<Nc>::HotConfiguration(pRNG,Umu);
std::vector<LatticeColourMatrix> U(4,&Grid);

106
tests/core/Test_unary.cc Normal file
View File

@@ -0,0 +1,106 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: ./tests/Test_quenched_update.cc
Copyright (C) 2015
Author: Azusa Yamaguchi <ayamaguc@staffmail.ed.ac.uk>
Author: Peter Boyle <paboyle@ph.ed.ac.uk>
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License along
with this program; if not, write to the Free Software Foundation, Inc.,
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
See the full license in the file "LICENSE" in the top level distribution directory
*************************************************************************************/
/* END LEGAL */
#include <Grid/Grid.h>
using namespace std;
using namespace Grid;
;
int main (int argc, char ** argv)
{
Grid_init(&argc,&argv);
std::vector<int> latt({8,8,8,8});
GridCartesian * grid = SpaceTimeGrid::makeFourDimGrid(latt,
GridDefaultSimd(Nd,vComplexD::Nsimd()),
GridDefaultMpi());
GridCartesian * gridF = SpaceTimeGrid::makeFourDimGrid(latt,
GridDefaultSimd(Nd,vComplexF::Nsimd()),
GridDefaultMpi());
///////////////////////////////
// Configuration of known size
///////////////////////////////
LatticeColourMatrixD ident(grid);
LatticeColourMatrixD U(grid);
LatticeColourMatrixD tmp(grid);
LatticeColourMatrixD org(grid);
LatticeColourMatrixF UF(gridF);
LatticeGaugeField Umu(grid);
ident =1.0;
// RNG set up for test
std::vector<int> pseeds({1,2,3,4,5}); // once I caught a fish alive
std::vector<int> sseeds({6,7,8,9,10});// then i let it go again
GridParallelRNG pRNG(grid); pRNG.SeedFixedIntegers(pseeds);
GridSerialRNG sRNG; sRNG.SeedFixedIntegers(sseeds);
SU<Nc>::HotConfiguration(pRNG,Umu);
U = PeekIndex<LorentzIndex>(Umu,0);
org=U;
tmp= U*adj(U) - ident ;
RealD Def1 = norm2( tmp );
std::cout << " Defect1 "<<Def1<<std::endl;
tmp = U - org;
std::cout << "Diff1 "<<norm2(tmp)<<std::endl;
precisionChange(UF,U);
precisionChange(U,UF);
tmp= U*adj(U) - ident ;
RealD Def2 = norm2( tmp );
std::cout << " Defect2 "<<Def2<<std::endl;
tmp = U - org;
std::cout << "Diff2 "<<norm2(tmp)<<std::endl;
U = ProjectOnGroup(U);
tmp= U*adj(U) - ident ;
RealD Def3 = norm2( tmp);
std::cout << " Defect3 "<<Def3<<std::endl;
tmp = U - org;
std::cout << "Diff3 "<<norm2(tmp)<<std::endl;
Grid_finalize();
}

8
tests/core/Test_where.cc
View File

@@ -40,9 +40,9 @@ int main (int argc, char ** argv)
int N=16;
std::vector<int> latt_size ({N,4,4});
std::vector<int> simd_layout({vComplexD::Nsimd(),1,1});
std::vector<int> mpi_layout ({1,1,1});
std::vector<int> latt_size ({N,N,N,N});
std::vector<int> simd_layout({vComplexD::Nsimd(),1,1,1});
std::vector<int> mpi_layout ({1,1,1,1});
int vol = 1;
int nd = latt_size.size();
@@ -69,7 +69,7 @@ int main (int argc, char ** argv)
for(int t=0;t<latt_size[mu];t++){
LatticeCoordinate(coor,mu);
sl=where(coor==Integer(t),rn,zz);
std::cout <<GridLogMessage<< " sl " << sl<<std::endl;
// std::cout <<GridLogMessage<< " sl " << sl<<std::endl;
std::cout <<GridLogMessage<<" slice "<<t<<" " << norm2(sl)<<std::endl;
ns=ns+norm2(sl);
}

143
tests/core/Test_where_extended.cc Normal file
View File

@@ -0,0 +1,143 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: ./tests/Test_poisson_fft.cc
Copyright (C) 2015
Author: Azusa Yamaguchi <ayamaguc@staffmail.ed.ac.uk>
Author: Peter Boyle <paboyle@ph.ed.ac.uk>
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License along
with this program; if not, write to the Free Software Foundation, Inc.,
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
See the full license in the file "LICENSE" in the top level distribution directory
*************************************************************************************/
/* END LEGAL */
#include <Grid/Grid.h>
using namespace Grid;
;
int main (int argc, char ** argv)
{
Grid_init(&argc,&argv);
int threads = GridThread::GetThreads();
std::cout<<GridLogMessage << "Grid is setup to use "<<threads<<" threads"<<std::endl;
int N=16;
std::vector<int> latt_size ({N,4,4});
std::vector<int> simd_layout({vComplexD::Nsimd(),1,1});
std::vector<int> mpi_layout ({1,1,1});
int vol = 1;
int nd = latt_size.size();
for(int d=0;d<nd;d++){
vol = vol * latt_size[d];
}
GridCartesian GRID(latt_size,simd_layout,mpi_layout);
GridParallelRNG RNG(&GRID);
RNG.SeedFixedIntegers(std::vector<int>({45,12,81,9}));
std::cout<<GridLogMessage<<"=============================================================="<<std::endl;
std::cout<<GridLogMessage<<"== LatticeComplex =="<<std::endl;
std::cout<<GridLogMessage<<"=============================================================="<<std::endl;
{
LatticeComplexD zz(&GRID);
LatticeInteger coor(&GRID);
LatticeComplexD rn(&GRID);
LatticeComplexD sl(&GRID);
zz = ComplexD(0.0,0.0);
gaussian(RNG,rn);
RealD nn=norm2(rn);
for(int mu=0;mu<nd;mu++){
RealD ns=0.0;
for(int t=0;t<latt_size[mu];t++){
LatticeCoordinate(coor,mu);
sl=where(coor==Integer(t),rn,zz);
std::cout <<GridLogMessage<<" slice "<<t<<" " << norm2(sl)<<std::endl;
ns=ns+norm2(sl);
}
std::cout <<GridLogMessage <<" sliceNorm" <<mu<<" "<< nn <<" "<<ns<<" err " << nn-ns<<std::endl;
assert(abs(nn-ns) < 1.0e-10);
}
}
std::cout<<GridLogMessage<<"=============================================================="<<std::endl;
std::cout<<GridLogMessage<<"== LatticeFermion =="<<std::endl;
std::cout<<GridLogMessage<<"=============================================================="<<std::endl;
{
LatticeFermionD zz(&GRID);
LatticeInteger coor(&GRID);
LatticeFermionD rn(&GRID);
LatticeFermionD sl(&GRID);
zz = ComplexD(0.0,0.0);
gaussian(RNG,rn);
RealD nn=norm2(rn);
for(int mu=0;mu<nd;mu++){
RealD ns=0.0;
for(int t=0;t<latt_size[mu];t++){
LatticeCoordinate(coor,mu);
sl=where(coor==Integer(t),rn,zz);
std::cout <<GridLogMessage<<" slice "<<t<<" " << norm2(sl)<<std::endl;
ns=ns+norm2(sl);
}
std::cout <<GridLogMessage <<" sliceNorm" <<mu<<" "<< nn <<" "<<ns<<" err " << nn-ns<<std::endl;
assert(abs(nn-ns) < 1.0e-10);
}
}
std::cout<<GridLogMessage<<"=============================================================="<<std::endl;
std::cout<<GridLogMessage<<"== LatticePropagator =="<<std::endl;
std::cout<<GridLogMessage<<"=============================================================="<<std::endl;
{
LatticePropagatorD zz(&GRID);
LatticeInteger coor(&GRID);
LatticePropagatorD rn(&GRID);
LatticePropagatorD sl(&GRID);
zz = ComplexD(0.0,0.0);
gaussian(RNG,rn);
RealD nn=norm2(rn);
for(int mu=0;mu<nd;mu++){
RealD ns=0.0;
for(int t=0;t<latt_size[mu];t++){
LatticeCoordinate(coor,mu);
sl=where(coor==Integer(t),rn,zz);
std::cout <<GridLogMessage<<" slice "<<t<<" " << norm2(sl)<<std::endl;
ns=ns+norm2(sl);
}
std::cout <<GridLogMessage <<" sliceNorm" <<mu<<" "<< nn <<" "<<ns<<" err " << nn-ns<<std::endl;
assert(abs(nn-ns) < 1.0e-10);
}
}
Grid_finalize();
}

2
tests/core/Test_wilson_clover.cc
View File

@@ -74,7 +74,7 @@ int main(int argc, char **argv)
FermionField chi(&Grid);
random(pRNG, chi);
LatticeGaugeField Umu(&Grid);
SU3::HotConfiguration(pRNG, Umu);
SU<Nc>::HotConfiguration(pRNG, Umu);
std::vector<LatticeColourMatrix> U(4, &Grid);
double volume = 1;

2
tests/core/Test_wilson_even_odd.cc
View File

@@ -70,7 +70,7 @@ int main (int argc, char ** argv)
LatticeFermion tmp(&Grid); tmp=Zero();
LatticeFermion err(&Grid); tmp=Zero();
LatticeGaugeField Umu(&Grid);
SU3::HotConfiguration(pRNG,Umu);
SU<Nc>::HotConfiguration(pRNG,Umu);
std::vector<LatticeColourMatrix> U(4,&Grid);
double volume=1;

2
tests/core/Test_wilson_twisted_mass_even_odd.cc
View File

@@ -71,7 +71,7 @@ int main (int argc, char ** argv)
LatticeFermion ref(&Grid); ref=Zero();
LatticeFermion tmp(&Grid); tmp=Zero();
LatticeFermion err(&Grid); tmp=Zero();
LatticeGaugeField Umu(&Grid); SU3::HotConfiguration(pRNG,Umu);
LatticeGaugeField Umu(&Grid); SU<Nc>::HotConfiguration(pRNG,Umu);
std::vector<LatticeColourMatrix> U(4,&Grid);
double volume=1;

2
tests/debug/Test_cayley_cg.cc
View File

@@ -116,7 +116,7 @@ int main (int argc, char ** argv)
LatticeGaugeField Umu(UGrid);
LatticeGaugeFieldF UmuF(UGridF);
SU3::HotConfiguration(RNG4,Umu);
SU<Nc>::HotConfiguration(RNG4,Umu);
precisionChange(UmuF,Umu);
std::vector<LatticeColourMatrix> U(4,UGrid);

2
tests/debug/Test_cayley_coarsen_support.cc
View File

@@ -77,7 +77,7 @@ int main (int argc, char ** argv)
LatticeFermion ref(FGrid); ref=Zero();
LatticeFermion tmp(FGrid);
LatticeFermion err(FGrid);
LatticeGaugeField Umu(UGrid); SU3::HotConfiguration(RNG4,Umu);
LatticeGaugeField Umu(UGrid); SU<Nc>::HotConfiguration(RNG4,Umu);
#if 0
std::vector<LatticeColourMatrix> U(4,UGrid);

2
tests/debug/Test_cayley_even_odd.cc
View File

@@ -70,7 +70,7 @@ int main (int argc, char ** argv)
GridParallelRNG RNG5(FGrid); RNG5.SeedFixedIntegers(seeds5);
GridParallelRNG RNG4(UGrid); RNG4.SeedFixedIntegers(seeds4);
LatticeGaugeField Umu(UGrid); SU3::HotConfiguration(RNG4,Umu);
LatticeGaugeField Umu(UGrid); SU<Nc>::HotConfiguration(RNG4,Umu);
std::vector<LatticeColourMatrix> U(4,UGrid);
RealD mass=0.1;

6
tests/debug/Test_cayley_ldop_cr.cc
View File

@@ -71,9 +71,9 @@ int main (int argc, char ** argv)
std::string file("./ckpoint_lat.400");
NerscIO::readConfiguration(Umu,header,file);
// SU3::ColdConfiguration(RNG4,Umu);
// SU3::TepidConfiguration(RNG4,Umu);
// SU3::HotConfiguration(RNG4,Umu);
// SU<Nc>::ColdConfiguration(RNG4,Umu);
// SU<Nc>::TepidConfiguration(RNG4,Umu);
// SU<Nc>::HotConfiguration(RNG4,Umu);
// Umu=Zero();
RealD mass=0.1;

96
tests/debug/Test_cayley_mres.cc
View File

@@ -33,13 +33,14 @@ using namespace Grid;
template<class What>
void TestConserved(What & Ddwf, What & Ddwfrev,
void TestConserved(What & Ddwf,
LatticeGaugeField &Umu,
GridCartesian * FGrid, GridRedBlackCartesian * FrbGrid,
GridCartesian * UGrid, GridRedBlackCartesian * UrbGrid,
RealD mass, RealD M5,
GridParallelRNG *RNG4,
GridParallelRNG *RNG5);
GridParallelRNG *RNG5,
What *Ddwfrev=nullptr);
Gamma::Algebra Gmu [] = {
Gamma::Algebra::GammaX,
@@ -102,14 +103,25 @@ int main (int argc, char ** argv)
GridRedBlackCartesian * FrbGridF = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGridF);
std::vector<int> seeds4({1,2,3,4});
std::vector<int> seeds5({5,6,7,8});
GridParallelRNG RNG5(FGrid); RNG5.SeedFixedIntegers(seeds5);
GridParallelRNG RNG4(UGrid); RNG4.SeedFixedIntegers(seeds4);
GridParallelRNG RNG4(UGrid);
std::vector<int> seeds4({1,2,3,4}); RNG4.SeedFixedIntegers(seeds4);
//const std::string seeds4{ "test-gauge-3000" }; RNG4.SeedUniqueString( seeds4 );
LatticeGaugeField Umu(UGrid);
SU3::ColdConfiguration(Umu);
// SU3::HotConfiguration(RNG4,Umu);
if( argc > 1 && argv[1][0] != '-' )
{
std::cout<<GridLogMessage <<"Loading configuration from "<<argv[1]<<std::endl;
FieldMetaData header;
NerscIO::readConfiguration(Umu, header, argv[1]);
}
else
{
std::cout<<GridLogMessage <<"Using hot configuration"<<std::endl;
// SU<Nc>::ColdConfiguration(Umu);
SU<Nc>::HotConfiguration(RNG4,Umu);
}
RealD mass=0.3;
RealD M5 =1.0;
@@ -117,7 +129,7 @@ int main (int argc, char ** argv)
std::cout<<GridLogMessage <<"DomainWallFermion test"<<std::endl;
std::cout<<GridLogMessage <<"======================"<<std::endl;
DomainWallFermionR Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5);
TestConserved<DomainWallFermionR>(Ddwf,Ddwf,Umu,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5);
TestConserved<DomainWallFermionR>(Ddwf,Umu,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5);
RealD b=1.5;// Scale factor b+c=2, b-c=1
RealD c=0.5;
@@ -127,13 +139,13 @@ int main (int argc, char ** argv)
std::cout<<GridLogMessage <<"MobiusFermion test"<<std::endl;
std::cout<<GridLogMessage <<"======================"<<std::endl;
MobiusFermionR Dmob(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,b,c);
TestConserved<MobiusFermionR>(Dmob,Dmob,Umu,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5);
TestConserved<MobiusFermionR>(Dmob,Umu,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5);
std::cout<<GridLogMessage <<"======================"<<std::endl;
std::cout<<GridLogMessage <<"ScaledShamirFermion test"<<std::endl;
std::cout<<GridLogMessage <<"======================"<<std::endl;
ScaledShamirFermionR Dsham(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,2.0);
TestConserved<ScaledShamirFermionR>(Dsham,Dsham,Umu,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5);
TestConserved<ScaledShamirFermionR>(Dsham,Umu,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5);
std::cout<<GridLogMessage <<"======================"<<std::endl;
std::cout<<GridLogMessage <<"ZMobiusFermion test"<<std::endl;
@@ -142,8 +154,7 @@ int main (int argc, char ** argv)
// for(int s=0;s<Ls;s++) omegasrev[s]=omegas[Ls-1-s];
ZMobiusFermionR ZDmob(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,omegas,b,c);
ZMobiusFermionR ZDmobrev(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,omegasrev,b,c);
TestConserved<ZMobiusFermionR>(ZDmob,ZDmobrev,Umu,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5);
TestConserved<ZMobiusFermionR>(ZDmob,Umu,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5,&ZDmobrev);
Grid_finalize();
}
@@ -151,22 +162,17 @@ int main (int argc, char ** argv)
template<class Action>
void TestConserved(Action & Ddwf,
Action & Ddwfrev,
void TestConserved(Action & Ddwf,
LatticeGaugeField &Umu,
GridCartesian * FGrid, GridRedBlackCartesian * FrbGrid,
GridCartesian * UGrid, GridRedBlackCartesian * UrbGrid,
RealD mass, RealD M5,
GridParallelRNG *RNG4,
GridParallelRNG *RNG5)
GridParallelRNG *RNG5,
Action * Ddwfrev)
{
int Ls=Ddwf.Ls;
LatticePropagator phys_src(UGrid);
std::vector<LatticeColourMatrix> U(4,UGrid);
LatticePropagator seqsrc(FGrid);
LatticePropagator phys_src(UGrid);
LatticePropagator seqsrc(FGrid);
LatticePropagator prop5(FGrid);
LatticePropagator prop5rev(FGrid);
LatticePropagator prop4(UGrid);
@@ -184,9 +190,9 @@ void TestConserved(Action & Ddwf,
phys_src=Zero();
pokeSite(kronecker,phys_src,coor);
MdagMLinearOperator<Action,LatticeFermion> HermOp(Ddwf);
MdagMLinearOperator<Action,LatticeFermion> HermOprev(Ddwfrev);
ConjugateGradient<LatticeFermion> CG(1.0e-16,100000);
SchurRedBlackDiagTwoSolve<LatticeFermion> schur(CG);
ZeroGuesser<LatticeFermion> zpg;
for(int s=0;s<Nd;s++){
for(int c=0;c<Nc;c++){
LatticeFermion src4 (UGrid);
@@ -196,20 +202,20 @@ void TestConserved(Action & Ddwf,
Ddwf.ImportPhysicalFermionSource(src4,src5);
LatticeFermion result5(FGrid); result5=Zero();
// CGNE
LatticeFermion Mdagsrc5 (FGrid);
Ddwf.Mdag(src5,Mdagsrc5);
CG(HermOp,Mdagsrc5,result5);
schur(Ddwf,src5,result5,zpg);
std::cout<<GridLogMessage<<"spin "<<s<<" color "<<c<<" norm2(sourc5d) "<<norm2(src5)
<<" norm2(result5d) "<<norm2(result5)<<std::endl;
FermToProp<Action>(prop5,result5,s,c);
LatticeFermion result4(UGrid);
Ddwf.ExportPhysicalFermionSolution(result5,result4);
FermToProp<Action>(prop4,result4,s,c);
Ddwfrev.ImportPhysicalFermionSource(src4,src5);
Ddwfrev.Mdag(src5,Mdagsrc5);
CG(HermOprev,Mdagsrc5,result5);
if( Ddwfrev ) {
Ddwfrev->ImportPhysicalFermionSource(src4,src5);
result5 = Zero();
schur(*Ddwfrev,src5,result5,zpg);
}
FermToProp<Action>(prop5rev,result5,s,c);
}
}
@@ -241,11 +247,7 @@ void TestConserved(Action & Ddwf,
PropToFerm<Action>(src5,seqsrc,s,c);
LatticeFermion result5(FGrid); result5=Zero();
// CGNE
LatticeFermion Mdagsrc5 (FGrid);
Ddwf.Mdag(src5,Mdagsrc5);
CG(HermOp,Mdagsrc5,result5);
schur(Ddwf,src5,result5,zpg);
LatticeFermion result4(UGrid);
Ddwf.ExportPhysicalFermionSolution(result5,result4);
@@ -266,10 +268,10 @@ void TestConserved(Action & Ddwf,
Ddwf.ContractConservedCurrent(prop5rev,prop5,Vector_mu,phys_src,Current::Vector,Tdir);
Ddwf.ContractJ5q(prop5,PJ5q);
PA = trace(g5*Axial_mu);
SV = trace(Vector_mu);
VV = trace(gT*Vector_mu);
PP = trace(adj(prop4)*prop4);
PA = trace(g5*Axial_mu); // Pseudoscalar-Axial conserved current
SV = trace(Vector_mu); // Scalar-Vector conserved current
VV = trace(gT*Vector_mu); // (local) Vector-Vector conserved current
PP = trace(adj(prop4)*prop4); // Pseudoscalar density
// Spatial sum
sliceSum(PA,sumPA,Tdir);
@@ -278,15 +280,17 @@ void TestConserved(Action & Ddwf,
sliceSum(PP,sumPP,Tdir);
sliceSum(PJ5q,sumPJ5q,Tdir);
int Nt=sumPA.size();
const int Nt{static_cast<int>(sumPA.size())};
std::cout<<GridLogMessage<<"Vector Ward identity by timeslice (~ 0)"<<std::endl;
for(int t=0;t<Nt;t++){
std::cout <<" SV "<<real(TensorRemove(sumSV[t]));
std::cout <<" VV "<<real(TensorRemove(sumVV[t]))<<std::endl;
std::cout<<GridLogMessage <<" t "<<t<<" SV "<<real(TensorRemove(sumSV[t]))<<" VV "<<real(TensorRemove(sumVV[t]))<<std::endl;
}
std::cout<<GridLogMessage<<"Axial Ward identity by timeslice (defect ~ 0)"<<std::endl;
for(int t=0;t<Nt;t++){
std::cout <<" PAc "<<real(TensorRemove(sumPA[t]));
std::cout <<" PJ5q "<<real(TensorRemove(sumPJ5q[t]));
std::cout <<" Ward Identity defect " <<real(TensorRemove(sumPA[t]-sumPA[(t-1+Nt)%Nt] - 2.0*(Ddwf.mass*sumPP[t] + sumPJ5q[t]) ))<<"\n";
const RealD DmuPAmu{real(TensorRemove(sumPA[t]-sumPA[(t-1+Nt)%Nt]))};
std::cout<<GridLogMessage<<" t "<<t<<" DmuPAmu "<<DmuPAmu
<<" PP "<<real(TensorRemove(sumPP[t]))<<" PJ5q "<<real(TensorRemove(sumPJ5q[t]))
<<" Ward Identity defect " <<(DmuPAmu - 2.*real(TensorRemove(Ddwf.mass*sumPP[t] + sumPJ5q[t])))<<std::endl;
}
///////////////////////////////

8
tests/debug/Test_heatbath_dwf_eofa.cc
View File

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

9
tests/debug/Test_heatbath_dwf_eofa_gparity.cc
View File

@@ -74,10 +74,13 @@ int main(int argc, char** argv)
RNG5.SeedFixedIntegers(seeds5);
GridParallelRNG RNG4(UGrid);
RNG4.SeedFixedIntegers(seeds4);
GridSerialRNG sRNG;
RNG4.SeedFixedIntegers(seeds4);
sRNG.SeedFixedIntegers(seeds5);
// Random gauge field
LatticeGaugeField Umu(UGrid);
SU3::HotConfiguration(RNG4, Umu);
SU<Nc>::HotConfiguration(RNG4, Umu);
// GparityDomainWallFermionR::ImplParams params;
FermionAction::ImplParams params;
@@ -90,7 +93,7 @@ int main(int argc, char** argv)
ConjugateGradient<FermionField> CG(1.0e-12, 5000);
ExactOneFlavourRatioPseudoFermionAction<FermionImplPolicy> Meofa(Lop, Rop, CG, Params, false);
Meofa.refresh(Umu, RNG5);
Meofa.refresh(Umu,sRNG, RNG5);
printf("<Phi|Meofa|Phi> = %1.15e\n", Meofa.S(Umu));
}
@@ -100,7 +103,7 @@ int main(int argc, char** argv)
ConjugateGradient<FermionField> CG(1.0e-12, 5000);
ExactOneFlavourRatioPseudoFermionAction<FermionImplPolicy> Meofa(Lop, Rop, CG, Params, true);
Meofa.refresh(Umu, RNG5);
Meofa.refresh(Umu,sRNG, RNG5);
printf("<Phi|Meofa|Phi> = %1.15e\n", Meofa.S(Umu));
}

8
tests/debug/Test_heatbath_mobius_eofa.cc
View File

@@ -68,14 +68,16 @@ int main(int argc, char** argv)
// Set up RNGs
std::vector<int> seeds4({1, 2, 3, 4});
std::vector<int> seeds5({5, 6, 7, 8});
GridSerialRNG sRNG;
GridParallelRNG RNG5(FGrid);
RNG5.SeedFixedIntegers(seeds5);
sRNG.SeedFixedIntegers(seeds5);
GridParallelRNG RNG4(UGrid);
RNG4.SeedFixedIntegers(seeds4);
// Random gauge field
LatticeGaugeField Umu(UGrid);
SU3::HotConfiguration(RNG4, Umu);
SU<Nc>::HotConfiguration(RNG4, Umu);
MobiusEOFAFermionR Lop(Umu, *FGrid, *FrbGrid, *UGrid, *UrbGrid, mf, mf, mpv, 0.0, -1, M5, b, c);
MobiusEOFAFermionR Rop(Umu, *FGrid, *FrbGrid, *UGrid, *UrbGrid, mpv, mf, mpv, -1.0, 1, M5, b, c);
@@ -86,7 +88,7 @@ int main(int argc, char** argv)
ConjugateGradient<LatticeFermion> CG(1.0e-12, 5000);
ExactOneFlavourRatioPseudoFermionAction<WilsonImplR> Meofa(Lop, Rop, CG, Params, false);
Meofa.refresh(Umu, RNG5);
Meofa.refresh(Umu, sRNG,RNG5);
printf("<Phi|Meofa|Phi> = %1.15e\n", Meofa.S(Umu));
}
@@ -96,7 +98,7 @@ int main(int argc, char** argv)
ConjugateGradient<LatticeFermion> CG(1.0e-12, 5000);
ExactOneFlavourRatioPseudoFermionAction<WilsonImplR> Meofa(Lop, Rop, CG, Params, true);
Meofa.refresh(Umu, RNG5);
Meofa.refresh(Umu, sRNG,RNG5);
printf("<Phi|Meofa|Phi> = %1.15e\n", Meofa.S(Umu));
}

8
tests/debug/Test_heatbath_mobius_eofa_gparity.cc
View File

@@ -73,13 +73,15 @@ int main(int argc, char** argv)
std::vector<int> seeds4({1, 2, 3, 4});
std::vector<int> seeds5({5, 6, 7, 8});
GridParallelRNG RNG5(FGrid);
GridSerialRNG sRNG;
RNG5.SeedFixedIntegers(seeds5);
sRNG.SeedFixedIntegers(seeds5);
GridParallelRNG RNG4(UGrid);
RNG4.SeedFixedIntegers(seeds4);
// Random gauge field
LatticeGaugeField Umu(UGrid);
SU3::HotConfiguration(RNG4, Umu);
SU<Nc>::HotConfiguration(RNG4, Umu);
FermionAction::ImplParams params;
FermionAction Lop(Umu, *FGrid, *FrbGrid, *UGrid, *UrbGrid, mf, mf, mpv, 0.0, -1, M5, b, c, params);
@@ -91,7 +93,7 @@ int main(int argc, char** argv)
ConjugateGradient<FermionField> CG(1.0e-12, 5000);
ExactOneFlavourRatioPseudoFermionAction<FermionImplPolicy> Meofa(Lop, Rop, CG, Params, false);
Meofa.refresh(Umu, RNG5);
Meofa.refresh(Umu, sRNG, RNG5);
printf("<Phi|Meofa|Phi> = %1.15e\n", Meofa.S(Umu));
}
@@ -101,7 +103,7 @@ int main(int argc, char** argv)
ConjugateGradient<FermionField> CG(1.0e-12, 5000);
ExactOneFlavourRatioPseudoFermionAction<FermionImplPolicy> Meofa(Lop, Rop, CG, Params, true);
Meofa.refresh(Umu, RNG5);
Meofa.refresh(Umu, sRNG, RNG5);
printf("<Phi|Meofa|Phi> = %1.15e\n", Meofa.S(Umu));
}

2
tests/debug/Test_reweight_dwf_eofa.cc
View File

@@ -102,7 +102,7 @@ int main(int argc, char **argv)
// Random gauge field
LatticeGaugeField Umu(UGrid);
SU3::HotConfiguration(RNG4, Umu);
SU<Nc>::HotConfiguration(RNG4, Umu);
// Initialize RHMC fermion operators
DomainWallFermionR Ddwf_f(Umu, *FGrid, *FrbGrid, *UGrid, *UrbGrid, mf, M5);

2
tests/debug/Test_reweight_dwf_eofa_gparity.cc
View File

@@ -104,7 +104,7 @@ int main(int argc, char **argv)
// Random gauge field
LatticeGaugeField Umu(UGrid);
SU3::HotConfiguration(RNG4, Umu);
SU<Nc>::HotConfiguration(RNG4, Umu);
// Initialize RHMC fermion operators
GparityDomainWallFermionR::ImplParams params;

2
tests/debug/Test_reweight_mobius_eofa.cc
View File

@@ -104,7 +104,7 @@ int main(int argc, char **argv)
// Random gauge field
LatticeGaugeField Umu(UGrid);
SU3::HotConfiguration(RNG4, Umu);
SU<Nc>::HotConfiguration(RNG4, Umu);
// Initialize RHMC fermion operators
MobiusFermionR Ddwf_f(Umu, *FGrid, *FrbGrid, *UGrid, *UrbGrid, mf, M5, b, c);

2
tests/debug/Test_reweight_mobius_eofa_gparity.cc
View File

@@ -106,7 +106,7 @@ int main(int argc, char **argv)
// Random gauge field
LatticeGaugeField Umu(UGrid);
SU3::HotConfiguration(RNG4, Umu);
SU<Nc>::HotConfiguration(RNG4, Umu);
// Initialize RHMC fermion operators
GparityDomainWallFermionR::ImplParams params;

4
tests/forces/Test_contfrac_force.cc
View File

@@ -59,7 +59,7 @@ int main (int argc, char ** argv)
LatticeGaugeField U(UGrid);
SU3::HotConfiguration(RNG4,U);
SU<Nc>::HotConfiguration(RNG4,U);
////////////////////////////////////
// Unmodified matrix element
@@ -93,7 +93,7 @@ int main (int argc, char ** argv)
for(int mu=0;mu<Nd;mu++){
SU3::GaussianFundamentalLieAlgebraMatrix(RNG4, mommu); // Traceless antihermitian momentum; gaussian in lie alg
SU<Nc>::GaussianFundamentalLieAlgebraMatrix(RNG4, mommu); // Traceless antihermitian momentum; gaussian in lie alg
PokeIndex<LorentzIndex>(mom,mommu,mu);

4
tests/forces/Test_dwf_force.cc
View File

@@ -60,7 +60,7 @@ int main (int argc, char ** argv)
LatticeGaugeField U(UGrid);
SU3::HotConfiguration(RNG4,U);
SU<Nc>::HotConfiguration(RNG4,U);
////////////////////////////////////
// Unmodified matrix element
@@ -94,7 +94,7 @@ int main (int argc, char ** argv)
for(int mu=0;mu<Nd;mu++){
SU3::GaussianFundamentalLieAlgebraMatrix(RNG4, mommu); // Traceless antihermitian momentum; gaussian in lie alg
SU<Nc>::GaussianFundamentalLieAlgebraMatrix(RNG4, mommu); // Traceless antihermitian momentum; gaussian in lie alg
PokeIndex<LorentzIndex>(mom,mommu,mu);

8
tests/forces/Test_dwf_force_eofa.cc
View File

@@ -72,7 +72,7 @@ int main (int argc, char** argv)
LatticeFermion MphiPrime (FGrid);
LatticeGaugeField U(UGrid);
SU3::HotConfiguration(RNG4,U);
SU<Nc>::HotConfiguration(RNG4,U);
////////////////////////////////////
// Unmodified matrix element
@@ -86,7 +86,9 @@ int main (int argc, char** argv)
ConjugateGradient<LatticeFermion> CG(1.0e-12, 5000);
ExactOneFlavourRatioPseudoFermionAction<WilsonImplR> Meofa(Lop, Rop, CG, CG, CG, CG, CG, Params, true);
Meofa.refresh(U, RNG5);
GridSerialRNG sRNG; sRNG.SeedFixedIntegers(seeds4);
Meofa.refresh(U, sRNG, RNG5 );
RealD S = Meofa.S(U); // pdag M p
// get the deriv of phidag M phi with respect to "U"
@@ -105,7 +107,7 @@ int main (int argc, char** argv)
for(int mu=0; mu<Nd; mu++){
SU3::GaussianFundamentalLieAlgebraMatrix(RNG4, mommu); // Traceless antihermitian momentum; gaussian in lie alg
SU<Nc>::GaussianFundamentalLieAlgebraMatrix(RNG4, mommu); // Traceless antihermitian momentum; gaussian in lie alg
PokeIndex<LorentzIndex>(mom, mommu, mu);

129
tests/forces/Test_dwf_gpforce.cc
View File

@@ -63,8 +63,8 @@ int main (int argc, char ** argv)
LatticeGaugeField U(UGrid);
SU3::HotConfiguration(RNG4,U);
// SU3::ColdConfiguration(pRNG,U);
SU<Nc>::HotConfiguration(RNG4,U);
// SU<Nc>::ColdConfiguration(pRNG,U);
////////////////////////////////////
// Unmodified matrix element
@@ -84,6 +84,13 @@ int main (int argc, char ** argv)
GparityDomainWallFermionR::ImplParams params;
params.twists = twists;
/*
params.boundary_phases[0] = 1.0;
params.boundary_phases[1] = 1.0;
params.boundary_phases[2] = 1.0;
params.boundary_phases[3] =- 1.0;
*/
GparityDomainWallFermionR Dw(U,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,params);
Dw.M (phi,Mphi);
@@ -96,6 +103,16 @@ int main (int argc, char ** argv)
Dw.MDeriv(tmp , Mphi, phi,DaggerNo ); UdSdU=tmp;
Dw.MDeriv(tmp , phi, Mphi,DaggerYes ); UdSdU=(UdSdU+tmp);
// *****************************************************************************************
// *** There is a funny negative sign in all derivatives. This is - UdSdU. ***
// *** ***
// *** Deriv in both Wilson gauge action and the TwoFlavour.h seems to miss a minus sign ***
// *** UdSdU is negated relative to what I think - call what is returned mUdSdU, ***
// *** and insert minus sign ***
// *****************************************************************************************
UdSdU = - UdSdU ; // Follow sign convention of actions in Grid. Seems crazy.
FermionField Ftmp (FGrid);
@@ -106,18 +123,28 @@ int main (int argc, char ** argv)
RealD Hmom = 0.0;
RealD Hmomprime = 0.0;
LatticeColourMatrix mommu(UGrid);
LatticeColourMatrix forcemu(UGrid);
LatticeColourMatrix mUdSdUmu(UGrid);
LatticeGaugeField mom(UGrid);
LatticeGaugeField Uprime(UGrid);
for(int mu=0;mu<Nd;mu++){
SU3::GaussianFundamentalLieAlgebraMatrix(RNG4, mommu); // Traceless antihermitian momentum; gaussian in lie alg
SU<Nc>::GaussianFundamentalLieAlgebraMatrix(RNG4, mommu); // Traceless antihermitian momentum; gaussian in lie alg
Hmom -= real(sum(trace(mommu*mommu)));
// Momentum Hamiltonian is - trace(p^2)/HMC_MOM_DENOMINATOR
//
// Integrator.h: RealD H = - FieldImplementation::FieldSquareNorm(P)/HMC_MOMENTUM_DENOMINATOR; // - trace (P*P)/denom // GaugeImplTypes.h: Hloc += trace(Pmu * Pmu);
// Sign comes from a sneaky multiply by "i" in GaussianFundemantalLie algebra
// P is i P^a_\mu T^a, not Pa Ta
//
// Integrator.h: H = Hmom + sum S(action)
Hmom -= real(sum(trace(mommu*mommu)))/ HMC_MOMENTUM_DENOMINATOR;
PokeIndex<LorentzIndex>(mom,mommu,mu);
// -- Drops factor of "i" in the U update: U' = e^{P dt} U [ _not_ e^{iPdt}U ]. P is anti hermitian already
// -- Udot = p U
// fourth order exponential approx
autoView( mom_v, mom, CpuRead);
autoView( U_v , U, CpuRead);
@@ -134,8 +161,8 @@ int main (int argc, char ** argv)
;
});
}
std::cout << GridLogMessage <<"Initial mom hamiltonian is "<< Hmom <<std::endl;
Dw.ImportGauge(Uprime);
Dw.M (phi,MphiPrime);
@@ -145,53 +172,60 @@ int main (int argc, char ** argv)
// Use derivative to estimate dS
//////////////////////////////////////////////
for(int mu=0;mu<Nd;mu++){
std::cout << "" <<std::endl;
mommu = PeekIndex<LorentzIndex>(mom,mu);
std::cout << GridLogMessage<< " Mommu " << norm2(mommu)<<std::endl;
mommu = mommu+adj(mommu);
std::cout << GridLogMessage<< " Mommu + Mommudag " << norm2(mommu)<<std::endl;
mommu = PeekIndex<LorentzIndex>(UdSdU,mu);
std::cout << GridLogMessage<< " dsdumu " << norm2(mommu)<<std::endl;
mommu = mommu+adj(mommu);
std::cout << GridLogMessage<< " dsdumu + dag " << norm2(mommu)<<std::endl;
}
//
// Ta has 1/2([ F - adj(F) ])_traceless and want the UdSdU _and_ UdagdSdUdag terms so 2x.
//
LatticeComplex dS(UGrid); dS = Zero();
LatticeComplex dSmom(UGrid); dSmom = Zero();
LatticeComplex dSmom2(UGrid); dSmom2 = Zero();
for(int mu=0;mu<Nd;mu++){
mommu = PeekIndex<LorentzIndex>(UdSdU,mu);
mommu=Ta(mommu)*2.0;
mommu=Ta(mommu); // projectForce , GaugeImplTypes.h
PokeIndex<LorentzIndex>(UdSdU,mommu,mu);
}
for(int mu=0;mu<Nd;mu++){
mommu = PeekIndex<LorentzIndex>(mom,mu);
std::cout << GridLogMessage<< " Mommu " << norm2(mommu)<<std::endl;
mommu = mommu+adj(mommu);
std::cout << GridLogMessage<< " Mommu + Mommudag " << norm2(mommu)<<std::endl;
mommu = PeekIndex<LorentzIndex>(UdSdU,mu);
std::cout << GridLogMessage<< " dsdumu " << norm2(mommu)<<std::endl;
mommu = mommu+adj(mommu);
std::cout << GridLogMessage<< " dsdumu + dag " << norm2(mommu)<<std::endl;
}
for(int mu=0;mu<Nd;mu++){
forcemu = PeekIndex<LorentzIndex>(UdSdU,mu);
mUdSdUmu= PeekIndex<LorentzIndex>(UdSdU,mu);
mommu = PeekIndex<LorentzIndex>(mom,mu);
// Update PF action density
dS = dS+trace(mommu*forcemu)*dt;
//
// Derive HMC eom:
//
// Sdot = - 2 trace( p p^dot ) / D - trace( p [ mUdSdU - h.c. ] ) = 0
//
//
// Sdot = 0 = - 2 trace( p p^dot ) / D - 2 trace( p Ta( mUdSdU ) = 0
//
// EOM:
//
// pdot = - D Ta( mUdSdU ) -- source of sign is the "funny sign" above
//
// dSqcd_dt = - 2.0*trace(mommu* Ta(mUdSdU) )*dt -- i.e. mUdSdU with adjoint term -> force has a 2x implicit
//
// dH_mom/dt = - 2 trace (p pdot)/Denom
//
// dH_tot / dt = 0 <= pdot = - Denom * mUdSdU
//
// dH_mom/dt = 2 trace (p mUdSdU )
//
// True Momentum delta H has a dt^2 piece
//
// dSmom = [ trace mom*mom - trace ( (mom-Denom*f*dt)(mom-Denom*f*dt) ) ] / Denom
// = 2*trace(mom*f) dt - Denom*dt*dt * trace(f*f).
// = dSmom + dSmom2
//
dSmom = dSmom - trace(mommu*forcemu) * dt;
dSmom2 = dSmom2 - trace(forcemu*forcemu) *(0.25* dt*dt);
dS = dS - 2.0*trace(mommu*mUdSdUmu)*dt; // U and Udagger derivs hence 2x.
// Update mom action density
mommu = mommu + forcemu*(dt*0.5);
dSmom = dSmom + 2.0*trace(mommu*mUdSdUmu) * dt; // this 2.0 coms from derivative of p^2
dSmom2 = dSmom2 - trace(mUdSdUmu*mUdSdUmu) * dt*dt* HMC_MOMENTUM_DENOMINATOR; // Remnant
Hmomprime -= real(sum(trace(mommu*mommu)));
// Update mom action density . Verbatim update_P in Integrator.h
mommu = mommu - mUdSdUmu * dt* HMC_MOMENTUM_DENOMINATOR;;
Hmomprime -= real(sum(trace(mommu*mommu))) / HMC_MOMENTUM_DENOMINATOR;
}
@@ -199,20 +233,25 @@ int main (int argc, char ** argv)
ComplexD dSm = sum(dSmom);
ComplexD dSm2 = sum(dSmom2);
std::cout << GridLogMessage <<"dSm "<< dSm<<std::endl;
std::cout << GridLogMessage <<"dSm2 "<< dSm2<<std::endl;
std::cout << GridLogMessage <<"Initial mom hamiltonian is "<< Hmom <<std::endl;
std::cout << GridLogMessage <<"Final mom hamiltonian is "<< Hmomprime <<std::endl;
std::cout << GridLogMessage <<"Delta mom hamiltonian is "<< Hmomprime-Hmom <<std::endl;
std::cout << GridLogMessage << " S "<<S<<std::endl;
std::cout << GridLogMessage << " Sprime "<<Sprime<<std::endl;
std::cout << GridLogMessage << "dS "<<Sprime-S<<std::endl;
std::cout << GridLogMessage << "predict dS "<< dSpred <<std::endl;
std::cout << GridLogMessage <<"dSm "<< dSm<<std::endl;
std::cout << GridLogMessage <<"dSm2"<< dSm2<<std::endl;
std::cout << GridLogMessage <<"Delta mom hamiltonian is "<< Hmomprime-Hmom <<std::endl;
std::cout << GridLogMessage <<"predict Delta mom hamiltonian is "<< dSm+dSm2 <<std::endl;
std::cout << GridLogMessage << "Initial S "<<S<<std::endl;
std::cout << GridLogMessage << "Final S "<<Sprime<<std::endl;
std::cout << GridLogMessage << "Delta S "<<Sprime-S<<std::endl;
std::cout << GridLogMessage << "predict delta S"<< dSpred <<std::endl;
std::cout << GridLogMessage << "defect "<< Sprime-S-dSpred <<std::endl;
std::cout << GridLogMessage << "Total dS "<< Hmomprime - Hmom + Sprime - S <<std::endl;
std::cout << GridLogMessage << "dS - dt^2 term "<< Hmomprime - Hmom + Sprime - S - dSm2 <<std::endl;
assert( fabs(real(Sprime-S-dSpred)) < 5.0 ) ;
std::cout<< GridLogMessage << "Done" <<std::endl;

7
tests/forces/Test_dwf_gpforce_eofa.cc
View File

@@ -75,7 +75,7 @@ int main (int argc, char** argv)
FermionField MphiPrime (FGrid);
LatticeGaugeField U(UGrid);
SU3::HotConfiguration(RNG4,U);
SU<Nc>::HotConfiguration(RNG4,U);
////////////////////////////////////
// Unmodified matrix element
@@ -90,7 +90,8 @@ int main (int argc, char** argv)
ConjugateGradient<FermionField> CG(1.0e-12, 5000);
ExactOneFlavourRatioPseudoFermionAction<FermionImplPolicy> Meofa(Lop, Rop, CG, CG, CG, CG, CG, Params, true);
Meofa.refresh(U, RNG5);
GridSerialRNG sRNG; sRNG.SeedFixedIntegers(seeds4);
Meofa.refresh(U, sRNG, RNG5);
RealD S = Meofa.S(U); // pdag M p
// get the deriv of phidag M phi with respect to "U"
@@ -109,7 +110,7 @@ int main (int argc, char** argv)
for(int mu=0; mu<Nd; mu++){
SU3::GaussianFundamentalLieAlgebraMatrix(RNG4, mommu); // Traceless antihermitian momentum; gaussian in lie alg
SU<Nc>::GaussianFundamentalLieAlgebraMatrix(RNG4, mommu); // Traceless antihermitian momentum; gaussian in lie alg
PokeIndex<LorentzIndex>(mom, mommu, mu);

4
tests/forces/Test_gp_plaq_force.cc
View File

@@ -51,7 +51,7 @@ int main (int argc, char ** argv)
LatticeGaugeField U(&Grid);
SU3::HotConfiguration(pRNG,U);
SU<Nc>::HotConfiguration(pRNG,U);
double beta = 1.0;
ConjugateWilsonGaugeActionR Action(beta);
@@ -80,7 +80,7 @@ int main (int argc, char ** argv)
for(int mu=0;mu<Nd;mu++){
SU3::GaussianFundamentalLieAlgebraMatrix(pRNG, mommu); // Traceless antihermitian momentum; gaussian in lie alg
SU<Nc>::GaussianFundamentalLieAlgebraMatrix(pRNG, mommu); // Traceless antihermitian momentum; gaussian in lie alg
PokeIndex<LorentzIndex>(mom,mommu,mu);

9
tests/forces/Test_gp_rect_force.cc
View File

@@ -29,7 +29,6 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
using namespace std;
using namespace Grid;
;
@@ -54,11 +53,15 @@ int main (int argc, char ** argv)
LatticeGaugeField U(&Grid);
SU3::HotConfiguration(pRNG,U);
SU<Nc>::HotConfiguration(pRNG,U);
double beta = 1.0;
double c1 = 0.331;
const int nu = 1;
std::vector<int> twists(Nd,0);
twists[nu] = 1;
ConjugateGimplD::setDirections(twists);
ConjugatePlaqPlusRectangleActionR Action(beta,c1);
//ConjugateWilsonGaugeActionR Action(beta);
//WilsonGaugeActionR Action(beta);
@@ -82,7 +85,7 @@ int main (int argc, char ** argv)
for(int mu=0;mu<Nd;mu++){
SU3::GaussianFundamentalLieAlgebraMatrix(pRNG, mommu); // Traceless antihermitian momentum; gaussian in lie alg
SU<Nc>::GaussianFundamentalLieAlgebraMatrix(pRNG, mommu); // Traceless antihermitian momentum; gaussian in lie alg
PokeIndex<LorentzIndex>(mom,mommu,mu);

4
tests/forces/Test_gpdwf_force.cc
View File

@@ -63,7 +63,7 @@ int main (int argc, char ** argv)
LatticeGaugeField U(UGrid);
SU3::HotConfiguration(RNG4,U);
SU<Nc>::HotConfiguration(RNG4,U);
////////////////////////////////////
// Unmodified matrix element
@@ -100,7 +100,7 @@ int main (int argc, char ** argv)
for(int mu=0;mu<Nd;mu++){
SU3::GaussianFundamentalLieAlgebraMatrix(RNG4, mommu); // Traceless antihermitian momentum; gaussian in lie alg
SU<Nc>::GaussianFundamentalLieAlgebraMatrix(RNG4, mommu); // Traceless antihermitian momentum; gaussian in lie alg
PokeIndex<LorentzIndex>(mom,mommu,mu);

4
tests/forces/Test_gpwilson_force.cc
View File

@@ -57,7 +57,7 @@ int main (int argc, char ** argv)
LatticeGaugeField U(UGrid);
SU3::HotConfiguration(RNG4,U);
SU<Nc>::HotConfiguration(RNG4,U);
////////////////////////////////////
// Unmodified matrix element
@@ -94,7 +94,7 @@ int main (int argc, char ** argv)
for(int mu=0;mu<Nd;mu++){
// Traceless antihermitian momentum; gaussian in lie alg
SU3::GaussianFundamentalLieAlgebraMatrix(RNG4, mommu);
SU<Nc>::GaussianFundamentalLieAlgebraMatrix(RNG4, mommu);
PokeIndex<LorentzIndex>(mom,mommu,mu);

7
tests/forces/Test_laplacian_force.cc
View File

@@ -46,6 +46,7 @@ int main (int argc, char ** argv)
std::vector<int> seeds({1,2,3,4});
GridSerialRNG sRNG; sRNG.SeedFixedIntegers({4,5,6,7});
GridParallelRNG pRNG(&Grid);
pRNG.SeedFixedIntegers(std::vector<int>({15,91,21,3}));
@@ -58,7 +59,7 @@ int main (int argc, char ** argv)
PokeIndex<LorentzIndex>(P, P_mu, mu);
}
SU3::HotConfiguration(pRNG,U);
SU<Nc>::HotConfiguration(pRNG,U);
ConjugateGradient<LatticeGaugeField> CG(1.0e-8, 10000);
@@ -67,7 +68,7 @@ int main (int argc, char ** argv)
LaplacianAdjointField<PeriodicGimplR> Laplacian(&Grid, CG, LapPar, Kappa);
GeneralisedMomenta<PeriodicGimplR> LaplacianMomenta(&Grid, Laplacian);
LaplacianMomenta.M.ImportGauge(U);
LaplacianMomenta.MomentaDistribution(pRNG);// fills the Momenta with the correct distr
LaplacianMomenta.MomentaDistribution(sRNG,pRNG);// fills the Momenta with the correct distr
std::cout << std::setprecision(15);
@@ -95,7 +96,7 @@ int main (int argc, char ** argv)
std::cout << GridLogMessage << "Update the U " << std::endl;
for(int mu=0;mu<Nd;mu++){
// Traceless antihermitian momentum; gaussian in lie algebra
SU3::GaussianFundamentalLieAlgebraMatrix(pRNG, mommu);
SU<Nc>::GaussianFundamentalLieAlgebraMatrix(pRNG, mommu);
auto Umu = PeekIndex<LorentzIndex>(U, mu);
PokeIndex<LorentzIndex>(mom,mommu,mu);
Umu = expMat(mommu, dt, 12) * Umu;

113
tests/forces/Test_mobius_force.cc
View File

@@ -60,7 +60,7 @@ int main (int argc, char ** argv)
LatticeGaugeField U(UGrid);
SU3::HotConfiguration(RNG4,U);
SU<Nc>::HotConfiguration(RNG4,U);
////////////////////////////////////
// Unmodified matrix element
@@ -69,7 +69,14 @@ int main (int argc, char ** argv)
RealD M5=1.8;
RealD b=0.5;
RealD c=0.5;
MobiusFermionR Ddwf(U,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,b,c);
WilsonImplParams p;
p.boundary_phases[0] = 1.0;
p.boundary_phases[1] = 1.0;
p.boundary_phases[2] = 1.0;
p.boundary_phases[3] =- 1.0;
MobiusFermionR Ddwf(U,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,b,c,p);
Ddwf.M (phi,Mphi);
ComplexD S = innerProduct(Mphi,Mphi); // pdag MdagM p
@@ -82,24 +89,44 @@ int main (int argc, char ** argv)
Ddwf.MDeriv(tmp , Mphi, phi,DaggerNo ); UdSdU=tmp;
Ddwf.MDeriv(tmp , phi, Mphi,DaggerYes ); UdSdU=(UdSdU+tmp);
// *****************************************************************************************
// *** There is a funny negative sign in all derivatives. This is - UdSdU. ***
// *** ***
// *** Deriv in both Wilson gauge action and the TwoFlavour.h seems to miss a minus sign ***
// *** UdSdU is negated relative to what I think - call what is returned mUdSdU, ***
// *** and insert minus sign ***
// *****************************************************************************************
UdSdU = - UdSdU ; // Follow sign convention of actions in Grid. Seems crazy.
LatticeFermion Ftmp (FGrid);
////////////////////////////////////
// Modify the gauge field a little
////////////////////////////////////
RealD dt = 0.0001;
RealD dt = 0.001;
RealD Hmom = 0.0;
RealD Hmomprime = 0.0;
LatticeColourMatrix mommu(UGrid);
LatticeColourMatrix forcemu(UGrid);
LatticeColourMatrix mUdSdUmu(UGrid);
LatticeGaugeField mom(UGrid);
LatticeGaugeField Uprime(UGrid);
for(int mu=0;mu<Nd;mu++){
SU3::GaussianFundamentalLieAlgebraMatrix(RNG4, mommu); // Traceless antihermitian momentum; gaussian in lie alg
SU<Nc>::GaussianFundamentalLieAlgebraMatrix(RNG4, mommu); // Traceless antihermitian momentum; gaussian in lie alg
PokeIndex<LorentzIndex>(mom,mommu,mu);
// Momentum Hamiltonian is - trace(p^2)/HMC_MOM_DENOMINATOR
//
// Integrator.h: RealD H = - FieldImplementation::FieldSquareNorm(P)/HMC_MOMENTUM_DENOMINATOR; // - trace (P*P)/denom // GaugeImplTypes.h: Hloc += trace(Pmu * Pmu);
// Sign comes from a sneaky multiply by "i" in GaussianFundemantalLie algebra
// P is i P^a_\mu T^a, not Pa Ta
//
// Integrator.h: H = Hmom + sum S(action)
Hmom -= real(sum(trace(mommu*mommu)))/ HMC_MOMENTUM_DENOMINATOR;
// fourth order exponential approx
autoView( U_v , U, CpuRead);
autoView( mom_v, mom, CpuRead);
@@ -115,6 +142,7 @@ int main (int argc, char ** argv)
;
});
}
std::cout << GridLogMessage <<"Initial mom hamiltonian is "<< Hmom <<std::endl;
Ddwf.ImportGauge(Uprime);
Ddwf.M (phi,MphiPrime);
@@ -125,32 +153,87 @@ int main (int argc, char ** argv)
// Use derivative to estimate dS
//////////////////////////////////////////////
LatticeComplex dS(UGrid); dS = Zero();
LatticeComplex dSmom(UGrid); dSmom = Zero();
LatticeComplex dSmom2(UGrid); dSmom2 = Zero();
for(int mu=0;mu<Nd;mu++){
mommu = PeekIndex<LorentzIndex>(UdSdU,mu);
mommu=Ta(mommu)*2.0;
mommu=Ta(mommu);
PokeIndex<LorentzIndex>(UdSdU,mommu,mu);
}
for(int mu=0;mu<Nd;mu++){
forcemu = PeekIndex<LorentzIndex>(UdSdU,mu);
mUdSdUmu= PeekIndex<LorentzIndex>(UdSdU,mu);
mommu = PeekIndex<LorentzIndex>(mom,mu);
// Update PF action density
dS = dS+trace(mommu*forcemu)*dt;
//
// Derive HMC eom:
//
// Sdot = - 2 trace( p p^dot ) / D - trace( p [ mUdSdU - h.c. ] ) = 0
//
//
// Sdot = 0 = - 2 trace( p p^dot ) / D - 2 trace( p Ta( mUdSdU ) = 0
//
// EOM:
//
// pdot = - D Ta( mUdSdU ) -- source of sign is the "funny sign" above
//
// dSqcd_dt = - 2.0*trace(mommu* Ta(mUdSdU) )*dt -- i.e. mUdSdU with adjoint term -> force has a 2x implicit
//
// dH_mom/dt = - 2 trace (p pdot)/Denom
//
// dH_tot / dt = 0 <= pdot = - Denom * mUdSdU
//
// dH_mom/dt = 2 trace (p mUdSdU )
//
// True Momentum delta H has a dt^2 piece
//
// dSmom = [ trace mom*mom - trace ( (mom-Denom*f*dt)(mom-Denom*f*dt) ) ] / Denom
// = 2*trace(mom*f) dt - Denom*dt*dt * trace(f*f).
// = dSmom + dSmom2
//
dS = dS - 2.0*trace(mommu*mUdSdUmu)*dt; // U and Udagger derivs hence 2x.
dSmom = dSmom + 2.0*trace(mommu*mUdSdUmu) * dt; // this 2.0 coms from derivative of p^2
dSmom2 = dSmom2 - trace(mUdSdUmu*mUdSdUmu) * dt*dt* HMC_MOMENTUM_DENOMINATOR; // Remnant
mommu = mommu - mUdSdUmu * dt* HMC_MOMENTUM_DENOMINATOR;;
Hmomprime -= real(sum(trace(mommu*mommu))) / HMC_MOMENTUM_DENOMINATOR;
}
ComplexD dSpred = sum(dS);
ComplexD dSm = sum(dSmom);
ComplexD dSm2 = sum(dSmom2);
std::cout << GridLogMessage << " -- S "<<S<<std::endl;
std::cout << GridLogMessage << " -- Sprime "<<Sprime<<std::endl;
std::cout << GridLogMessage << "dS "<<Sprime-S<<std::endl;
std::cout << GridLogMessage << "predict dS "<< dSpred <<std::endl;
std::cout << GridLogMessage <<"dSm "<< dSm<<std::endl;
std::cout << GridLogMessage <<"dSm2 "<< dSm2<<std::endl;
std::cout << GridLogMessage <<"Initial mom hamiltonian is "<< Hmom <<std::endl;
std::cout << GridLogMessage <<"Final mom hamiltonian is "<< Hmomprime <<std::endl;
std::cout << GridLogMessage <<"Delta mom hamiltonian is "<< Hmomprime-Hmom <<std::endl;
std::cout << GridLogMessage <<"predict Delta mom hamiltonian is "<< dSm+dSm2 <<std::endl;
std::cout << GridLogMessage << "Initial S "<<S<<std::endl;
std::cout << GridLogMessage << "Final S "<<Sprime<<std::endl;
std::cout << GridLogMessage << "Delta S "<<Sprime-S<<std::endl;
std::cout << GridLogMessage << "predict delta S"<< dSpred <<std::endl;
std::cout << GridLogMessage << "defect "<< Sprime-S-dSpred <<std::endl;
std::cout << GridLogMessage << "Total dS "<< Hmomprime - Hmom + Sprime - S <<std::endl;
std::cout << GridLogMessage << "dS - dt^2 term "<< Hmomprime - Hmom + Sprime - S - dSm2 <<std::endl;
assert( fabs(real(Sprime-S-dSpred)) < 1.0 ) ;
std::cout<< GridLogMessage << "Done" <<std::endl;
Grid_finalize();
}

7
tests/forces/Test_mobius_force_eofa.cc
View File

@@ -72,7 +72,7 @@ int main (int argc, char** argv)
LatticeFermion MphiPrime (FGrid);
LatticeGaugeField U(UGrid);
SU3::HotConfiguration(RNG4,U);
SU<Nc>::HotConfiguration(RNG4,U);
////////////////////////////////////
// Unmodified matrix element
@@ -88,7 +88,8 @@ int main (int argc, char** argv)
ConjugateGradient<LatticeFermion> CG(1.0e-12, 5000);
ExactOneFlavourRatioPseudoFermionAction<WilsonImplR> Meofa(Lop, Rop, CG, CG, CG, CG, CG, Params, false);
Meofa.refresh(U, RNG5);
GridSerialRNG sRNG; sRNG.SeedFixedIntegers(seeds4);
Meofa.refresh(U, sRNG, RNG5 );
RealD S = Meofa.S(U); // pdag M p
// get the deriv of phidag M phi with respect to "U"
@@ -107,7 +108,7 @@ int main (int argc, char** argv)
for(int mu=0; mu<Nd; mu++){
SU3::GaussianFundamentalLieAlgebraMatrix(RNG4, mommu); // Traceless antihermitian momentum; gaussian in lie alg
SU<Nc>::GaussianFundamentalLieAlgebraMatrix(RNG4, mommu); // Traceless antihermitian momentum; gaussian in lie alg
PokeIndex<LorentzIndex>(mom, mommu, mu);

7
tests/forces/Test_mobius_gpforce_eofa.cc
View File

@@ -76,7 +76,7 @@ int main (int argc, char** argv)
FermionField MphiPrime (FGrid);
LatticeGaugeField U(UGrid);
SU3::HotConfiguration(RNG4,U);
SU<Nc>::HotConfiguration(RNG4,U);
////////////////////////////////////
// Unmodified matrix element
@@ -93,7 +93,8 @@ int main (int argc, char** argv)
ConjugateGradient<FermionField> CG(1.0e-12, 5000);
ExactOneFlavourRatioPseudoFermionAction<FermionImplPolicy> Meofa(Lop, Rop, CG, CG, CG, CG, CG, Params, false);
Meofa.refresh(U, RNG5);
GridSerialRNG sRNG; sRNG.SeedFixedIntegers(seeds4);
Meofa.refresh(U, sRNG, RNG5 );
RealD S = Meofa.S(U); // pdag M p
// get the deriv of phidag M phi with respect to "U"
@@ -112,7 +113,7 @@ int main (int argc, char** argv)
for(int mu=0; mu<Nd; mu++){
SU3::GaussianFundamentalLieAlgebraMatrix(RNG4, mommu); // Traceless antihermitian momentum; gaussian in lie alg
SU<Nc>::GaussianFundamentalLieAlgebraMatrix(RNG4, mommu); // Traceless antihermitian momentum; gaussian in lie alg
PokeIndex<LorentzIndex>(mom, mommu, mu);
autoView( U_v , U, CpuRead);

156
tests/forces/Test_momentum_filter.cc Normal file
View File

@@ -0,0 +1,156 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: ./tests/Test_wilson_force.cc
Copyright (C) 2015
Author: Christopher Kelly <ckelly@bnl.gov>
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License along
with this program; if not, write to the Free Software Foundation, Inc.,
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
See the full license in the file "LICENSE" in the top level distribution directory
*************************************************************************************/
/* END LEGAL */
#include <Grid/Grid.h>
using namespace std;
using namespace Grid;
//Get the mu-directected links on the upper boundary and the bulk remainder
template<typename Field>
void getLinksBoundaryBulk(Field &bound, Field &bulk, Field &from, const Coordinate &latt_size){
bound = Zero(); bulk = Zero();
for(int mu=0;mu<Nd;mu++){
LatticeInteger mucoor(bound.Grid());
LatticeCoordinate(mucoor, mu);
bound = where( mucoor == (Integer)(latt_size[mu] - 1), from, bound );
bulk = where( mucoor != (Integer)(latt_size[mu] - 1), from, bulk );
}
}
int main (int argc, char ** argv)
{
Grid_init(&argc,&argv);
Coordinate latt_size = GridDefaultLatt();
Coordinate simd_layout = GridDefaultSimd(Nd,vComplex::Nsimd());
Coordinate mpi_layout = GridDefaultMpi();
GridCartesian Grid(latt_size,simd_layout,mpi_layout);
GridRedBlackCartesian RBGrid(&Grid);
int threads = GridThread::GetThreads();
std::cout<<GridLogMessage << "Grid is setup to use "<<threads<<" threads"<<std::endl;
std::vector<int> seeds({1,2,3,4});
GridParallelRNG pRNG(&Grid);
GridSerialRNG sRNG;
pRNG.SeedFixedIntegers(seeds);
sRNG.SeedFixedIntegers(seeds);
typedef PeriodicGimplR Gimpl;
typedef WilsonGaugeAction<Gimpl> GaugeAction;
typedef NoHirep Representation; //fundamental
typedef NoSmearing<Gimpl> Smearing;
typedef MinimumNorm2<Gimpl, Smearing> Omelyan;
typedef Gimpl::Field Field;
typedef MomentumFilterApplyPhase<Field> Filter;
Filter filter(&Grid);
//Setup a filter that disables link update on links passing through the global lattice boundary
typedef Filter::LatticeLorentzScalarType MaskType;
typedef Filter::LorentzScalarType MaskSiteType;
MaskSiteType zero, one;
for(int mu=0;mu<Nd;mu++){
zero(mu)()() = 0.;
one(mu)()() = 1.;
}
MaskType zeroField(&Grid), oneField(&Grid);
zeroField = zero;
oneField = one;
filter.phase = oneField; //make every site 1.0
//Zero mu-directed links at upper boundary
for(int mu=0;mu<Nd;mu++){
LatticeInteger mucoor(&Grid);
LatticeCoordinate(mucoor, mu);
filter.phase = where( mucoor == (Integer)(latt_size[mu] - 1) , zeroField, filter.phase );
}
//Start with a random gauge field
Field U(&Grid);
SU<Nc>::HotConfiguration(pRNG,U);
//Get the original links on the bulk and boundary for later use
Field Ubnd_orig(&Grid), Ubulk_orig(&Grid);
getLinksBoundaryBulk(Ubnd_orig, Ubulk_orig, U, latt_size);
ActionSet<Field,Representation> actions(1);
double beta=6;
GaugeAction gauge_action(beta);
actions[0].push_back(&gauge_action);
Smearing smear;
IntegratorParameters params(1,1.); //1 MD step
Omelyan integrator(&Grid, params, actions, smear);
integrator.setMomentumFilter(filter);
integrator.refresh(U, sRNG, pRNG); //doesn't actually change the gauge field
//Check the momentum is zero on the boundary
const auto &P = integrator.getMomentum();
Field Pbnd(&Grid), Pbulk(&Grid);
getLinksBoundaryBulk(Pbnd, Pbulk, const_cast<Field&>(P), latt_size);
RealD Pbnd_nrm = norm2(Pbnd); //expect zero
std::cout << GridLogMessage << "After refresh, norm2 of mu-directed conjugate momentum on boundary is: " << Pbnd_nrm << " (expect 0)" << std::endl;
RealD Pbulk_nrm = norm2(Pbulk); //expect non-zero
std::cout << GridLogMessage << "After refresh, norm2 of bulk conjugate momentum is: " << Pbulk_nrm << " (expect non-zero)" << std::endl;
//Evolve the gauge field
integrator.integrate(U);
//Check momentum is still zero on boundary
getLinksBoundaryBulk(Pbnd, Pbulk, const_cast<Field&>(P), latt_size);
Pbnd_nrm = norm2(Pbnd); //expect zero
std::cout << GridLogMessage << "After integrate, norm2 of mu-directed conjugate momentum on boundary is: " << Pbnd_nrm << " (expect 0)" << std::endl;
Pbulk_nrm = norm2(Pbulk); //expect non-zero
std::cout << GridLogMessage << "After integrate, norm2 of bulk conjugate momentum is: " << Pbulk_nrm << " (expect non-zero)" << std::endl;
//Get the new bulk and bound links
Field Ubnd_new(&Grid), Ubulk_new(&Grid);
getLinksBoundaryBulk(Ubnd_new, Ubulk_new, U, latt_size);
Field Ubnd_diff = Ubnd_new - Ubnd_orig;
Field Ubulk_diff = Ubulk_new - Ubulk_orig;
RealD Ubnd_change = norm2( Ubnd_diff );
RealD Ubulk_change = norm2( Ubulk_diff );
std::cout << GridLogMessage << "After integrate, norm2 of change in mu-directed boundary links is : " << Ubnd_change << " (expect 0)" << std::endl;
std::cout << GridLogMessage << "After integrate, norm2 of change in bulk links is : " << Ubulk_change << " (expect non-zero)" << std::endl;
Grid_finalize();
}

4
tests/forces/Test_partfrac_force.cc
View File

@@ -62,7 +62,7 @@ int main (int argc, char ** argv)
LatticeGaugeField U(UGrid);
SU3::HotConfiguration(RNG4,U);
SU<Nc>::HotConfiguration(RNG4,U);
////////////////////////////////////
// Unmodified matrix element
@@ -96,7 +96,7 @@ int main (int argc, char ** argv)
for(int mu=0;mu<Nd;mu++){
SU3::GaussianFundamentalLieAlgebraMatrix(RNG4, mommu); // Traceless antihermitian momentum; gaussian in lie alg
SU<Nc>::GaussianFundamentalLieAlgebraMatrix(RNG4, mommu); // Traceless antihermitian momentum; gaussian in lie alg
PokeIndex<LorentzIndex>(mom,mommu,mu);

4
tests/forces/Test_rect_force.cc
View File

@@ -54,7 +54,7 @@ int main (int argc, char ** argv)
LatticeGaugeField U(&Grid);
SU3::HotConfiguration(pRNG,U);
SU<Nc>::HotConfiguration(pRNG,U);
double beta = 1.0;
double c1 = -0.331;
@@ -82,7 +82,7 @@ int main (int argc, char ** argv)
for(int mu=0;mu<Nd;mu++){
SU3::GaussianFundamentalLieAlgebraMatrix(pRNG, mommu); // Traceless antihermitian momentum; gaussian in lie alg
SU<Nc>::GaussianFundamentalLieAlgebraMatrix(pRNG, mommu); // Traceless antihermitian momentum; gaussian in lie alg
PokeIndex<LorentzIndex>(mom,mommu,mu);

4
tests/forces/Test_wilson_force.cc
View File

@@ -61,7 +61,7 @@ int main (int argc, char ** argv)
LatticeGaugeField U(&Grid);
//SU2::HotConfiguration(pRNG,U);
SU3::ColdConfiguration(pRNG,U);
SU<Nc>::ColdConfiguration(pRNG,U);
////////////////////////////////////
// Unmodified matrix element
@@ -98,7 +98,7 @@ int main (int argc, char ** argv)
for(int mu=0;mu<Nd;mu++){
// Traceless antihermitian momentum; gaussian in lie alg
SU3::GaussianFundamentalLieAlgebraMatrix(pRNG, mommu);
SU<Nc>::GaussianFundamentalLieAlgebraMatrix(pRNG, mommu);
Hmom -= real(sum(trace(mommu*mommu)));

6
tests/forces/Test_wilsonclover_force.cc
View File

@@ -62,8 +62,8 @@ int main(int argc, char **argv)
LatticeGaugeField U(&Grid);
SU3::HotConfiguration(pRNG, U);
//SU3::ColdConfiguration(pRNG, U);// Clover term Zero()
SU<Nc>::HotConfiguration(pRNG, U);
//SU<Nc>::ColdConfiguration(pRNG, U);// Clover term Zero()
////////////////////////////////////
// Unmodified matrix element
@@ -101,7 +101,7 @@ int main(int argc, char **argv)
for (int mu = 0; mu < Nd; mu++)
{
// Traceless antihermitian momentum; gaussian in lie alg
SU3::GaussianFundamentalLieAlgebraMatrix(pRNG, mommu);
SU<Nc>::GaussianFundamentalLieAlgebraMatrix(pRNG, mommu);
Hmom -= real(sum(trace(mommu * mommu)));
PokeIndex<LorentzIndex>(mom, mommu, mu);

4
tests/forces/Test_zmobius_force.cc
View File

@@ -59,7 +59,7 @@ int main (int argc, char ** argv)
LatticeGaugeField U(UGrid);
SU3::HotConfiguration(RNG4,U);
SU<Nc>::HotConfiguration(RNG4,U);
////////////////////////////////////
// Unmodified matrix element
@@ -109,7 +109,7 @@ int main (int argc, char ** argv)
for(int mu=0;mu<Nd;mu++){
SU3::GaussianFundamentalLieAlgebraMatrix(RNG4, mommu); // Traceless antihermitian momentum; gaussian in lie alg
SU<Nc>::GaussianFundamentalLieAlgebraMatrix(RNG4, mommu); // Traceless antihermitian momentum; gaussian in lie alg
PokeIndex<LorentzIndex>(mom,mommu,mu);

7
tests/hmc/Test_hmc_EODWFRatio_Gparity.cc
View File

@@ -81,6 +81,10 @@ int main(int argc, char **argv) {
// that have a complex construction
// standard
RealD beta = 5.6 ;
const int nu = 3;
std::vector<int> twists(Nd,0);
twists[nu] = 1;
ConjugateGimplD::setDirections(twists);
ConjugateWilsonGaugeActionR Waction(beta);
const int Ls = 8;
@@ -93,9 +97,6 @@ int main(int argc, char **argv) {
// temporarily need a gauge field
LatticeGaugeField U(GridPtr);
const int nu = 3;
std::vector<int> twists(Nd,0);
twists[nu] = 1;
FermionAction::ImplParams params;
params.twists = twists;
Real mass=0.04;

4
tests/hmc/Test_hmc_GparityIwasakiGauge.cc
View File

@@ -79,6 +79,10 @@ int main(int argc, char **argv) {
// that have a complex construction
// standard
RealD beta = 2.6 ;
const int nu = 3;
std::vector<int> twists(Nd,0);
twists[nu] = 1;
ConjugateGimplD::setDirections(twists);
ConjugateIwasakiGaugeActionR Waction(beta);

3
tests/hmc/Test_hmc_GparityWilsonGauge.cc
View File

@@ -80,6 +80,9 @@ int main(int argc, char **argv) {
// that have a complex construction
// standard
RealD beta = 5.6 ;
std::vector<int> twists(Nd,0);
twists[3] = 1;
ConjugateGimplD::setDirections(twists);
ConjugateWilsonGaugeActionR Waction(beta);

2
tests/lanczos/Test_dwf_compressed_lanczos_reorg_synthetic.cc
View File

@@ -293,7 +293,7 @@ int main (int argc, char ** argv) {
{
std::vector<int> seeds4({1,2,3,4});
GridParallelRNG RNG4(UGrid); RNG4.SeedFixedIntegers(seeds4);
SU3::HotConfiguration(RNG4, Umu);
SU<Nc>::HotConfiguration(RNG4, Umu);
}
std::cout << GridLogMessage << "Lattice dimensions: " << GridDefaultLatt() << " Ls: " << Ls << std::endl;

2
tests/lanczos/Test_dwf_lanczos.cc
View File

@@ -54,7 +54,7 @@ int main (int argc, char ** argv)
GridParallelRNG RNG5rb(FrbGrid); RNG5.SeedFixedIntegers(seeds5);
LatticeGaugeField Umu(UGrid);
SU3::HotConfiguration(RNG4, Umu);
SU<Nc>::HotConfiguration(RNG4, Umu);
std::vector<LatticeColourMatrix> U(4,UGrid);
for(int mu=0;mu<Nd;mu++){

2
tests/lanczos/Test_wilson_lanczos.cc
View File

@@ -61,7 +61,7 @@ int main(int argc, char** argv) {
RNG5.SeedFixedIntegers(seeds5);
LatticeGaugeField Umu(UGrid);
SU3::HotConfiguration(RNG4, Umu);
SU<Nc>::HotConfiguration(RNG4, Umu);
/*
std::vector<LatticeColourMatrix> U(4, UGrid);

2
tests/qdpxx/Test_qdpxx_baryon.cc
View File

@@ -280,7 +280,7 @@ void make_gauge(GaugeField &Umu, Grid::LatticePropagator &q1,Grid::LatticePropag
Grid::GridCartesian *UGrid = (Grid::GridCartesian *)Umu.Grid();
Grid::GridParallelRNG RNG4(UGrid);
RNG4.SeedFixedIntegers(seeds4);
Grid::SU3::HotConfiguration(RNG4, Umu);
Grid::SU<Nc>::HotConfiguration(RNG4, Umu);
// Propagator
Grid::gaussian(RNG4, q1);

2
tests/qdpxx/Test_qdpxx_loops_staples.cc
View File

@@ -277,7 +277,7 @@ double calc_grid_p(Grid::LatticeGaugeField & Umu)
Grid::GridCartesian * UGrid = (Grid::GridCartesian *) Umu.Grid();
Grid::GridParallelRNG RNG4(UGrid); RNG4.SeedFixedIntegers(seeds4);
Grid::SU3::HotConfiguration(RNG4,Umu);
Grid::SU<Nc>::HotConfiguration(RNG4,Umu);
Grid::LatticeColourMatrix tmp(UGrid);
tmp = Grid::zero;

2
tests/qdpxx/Test_qdpxx_munprec.cc
View File

@@ -502,7 +502,7 @@ void calc_grid(ChromaAction action,Grid::LatticeGaugeField & Umu, Grid::LatticeF
Grid::gaussian(RNG5,src);
Grid::gaussian(RNG5,res);
Grid::SU3::HotConfiguration(RNG4,Umu);
Grid::SU<Nc>::HotConfiguration(RNG4,Umu);
/*
Grid::LatticeColourMatrix U(UGrid);

2
tests/qdpxx/Test_qdpxx_stag.cc
View File

@@ -333,7 +333,7 @@ void make_gauge(GaugeField & Umu,FermionField &src)
Grid::GridCartesian * UGrid = (Grid::GridCartesian *) Umu.Grid();
Grid::GridParallelRNG RNG4(UGrid); RNG4.SeedFixedIntegers(seeds4);
Grid::SU3::HotConfiguration(RNG4,Umu);
Grid::SU<Nc>::HotConfiguration(RNG4,Umu);
Grid::gaussian(RNG4,src);
}

2
tests/qdpxx/Test_qdpxx_wilson.cc
View File

@@ -348,7 +348,7 @@ void make_gauge(GaugeField &Umu, FermionField &src)
Grid::GridCartesian *UGrid = (Grid::GridCartesian *)Umu._grid;
Grid::GridParallelRNG RNG4(UGrid);
RNG4.SeedFixedIntegers(seeds4);
Grid::SU3::HotConfiguration(RNG4, Umu);
Grid::SU<Nc>::HotConfiguration(RNG4, Umu);
// Fermion field
Grid::gaussian(RNG4, src);

4
tests/smearing/Test_smearing.cc
View File

@@ -47,8 +47,8 @@ int main (int argc, char ** argv)
RealD nrm = norm2(src);
LatticeFermion result(&Grid); result=Zero();
LatticeGaugeField Umu(&Grid);
// SU3::HotConfiguration(pRNG,Umu);
SU3::ColdConfiguration(Umu);
// SU<Nc>::HotConfiguration(pRNG,Umu);
SU<Nc>::ColdConfiguration(Umu);
std::vector<LatticeColourMatrix> U(4,&Grid);
for(int mu=0;mu<Nd;mu++){

2
tests/solver/Test_cf_cr_unprec.cc
View File

@@ -61,7 +61,7 @@ int main (int argc, char ** argv)
LatticeFermion src(FGrid); random(RNG5,src);
LatticeFermion result(FGrid); result=Zero();
LatticeGaugeField Umu(UGrid); SU3::HotConfiguration(RNG4,Umu);
LatticeGaugeField Umu(UGrid); SU<Nc>::HotConfiguration(RNG4,Umu);
std::vector<LatticeColourMatrix> U(4,UGrid);
for(int mu=0;mu<Nd;mu++){

475
tests/solver/Test_coarse_even_odd.cc Normal file
View File

@@ -0,0 +1,475 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: ./tests/solver/Test_coarse_even_odd.cc
Copyright (C) 2015-2020
Author: Daniel Richtmann <daniel.richtmann@ur.de>
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License along
with this program; if not, write to the Free Software Foundation, Inc.,
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
See the full license in the file "LICENSE" in the top level distribution directory
*************************************************************************************/
/* END LEGAL */
#include <Grid/Grid.h>
using namespace Grid;
#ifndef NBASIS
#define NBASIS 40
#endif
// NOTE: The tests in this file are written in analogy to
// - tests/core/Test_wilson_even_odd.cc
// - tests/core/Test_wilson_clover.cc
std::vector<int> readFromCommandlineIvec(int* argc,
char*** argv,
std::string&& option,
const std::vector<int>& defaultValue) {
std::string arg;
std::vector<int> ret(defaultValue);
if(GridCmdOptionExists(*argv, *argv + *argc, option)) {
arg = GridCmdOptionPayload(*argv, *argv + *argc, option);
GridCmdOptionIntVector(arg, ret);
}
return ret;
}
int main(int argc, char** argv) {
Grid_init(&argc, &argv);
/////////////////////////////////////////////////////////////////////////////
// Read from command line //
/////////////////////////////////////////////////////////////////////////////
const int nbasis = NBASIS; static_assert((nbasis & 0x1) == 0, "");
const int nb = nbasis/2;
Coordinate blockSize = readFromCommandlineIvec(&argc, &argv, "--blocksize", {2, 2, 2, 2});
std::cout << GridLogMessage << "Compiled with nbasis = " << nbasis << " -> nb = " << nb << std::endl;
/////////////////////////////////////////////////////////////////////////////
// General setup //
/////////////////////////////////////////////////////////////////////////////
Coordinate clatt = GridDefaultLatt();
for(int d=0; d<clatt.size(); d++) clatt[d] = clatt[d] / blockSize[d];
GridCartesian* Grid_f = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd, vComplex::Nsimd()), GridDefaultMpi());
GridCartesian* Grid_c = SpaceTimeGrid::makeFourDimGrid(clatt, GridDefaultSimd(Nd, vComplex::Nsimd()), GridDefaultMpi());
GridRedBlackCartesian* RBGrid_f = SpaceTimeGrid::makeFourDimRedBlackGrid(Grid_f);
GridRedBlackCartesian* RBGrid_c = SpaceTimeGrid::makeFourDimRedBlackGrid(Grid_c);
std::cout << GridLogMessage << "Grid_f:" << std::endl; Grid_f->show_decomposition();
std::cout << GridLogMessage << "Grid_c:" << std::endl; Grid_c->show_decomposition();
std::cout << GridLogMessage << "RBGrid_f:" << std::endl; RBGrid_f->show_decomposition();
std::cout << GridLogMessage << "RBGrid_c:" << std::endl; RBGrid_c->show_decomposition();
GridParallelRNG pRNG_f(Grid_f);
GridParallelRNG pRNG_c(Grid_c);
std::vector<int> seeds({1, 2, 3, 4});
pRNG_f.SeedFixedIntegers(seeds);
pRNG_c.SeedFixedIntegers(seeds);
/////////////////////////////////////////////////////////////////////////////
// Setup of Dirac Matrix and Operator //
/////////////////////////////////////////////////////////////////////////////
LatticeGaugeField Umu(Grid_f); SU3::HotConfiguration(pRNG_f, Umu);
RealD checkTolerance = (getPrecision<LatticeFermion>::value == 1) ? 1e-7 : 1e-15;
RealD mass = -0.30;
RealD csw = 1.9192;
WilsonCloverFermionR Dwc(Umu, *Grid_f, *RBGrid_f, mass, csw, csw);
MdagMLinearOperator<WilsonCloverFermionR, LatticeFermion> MdagMOp_Dwc(Dwc);
/////////////////////////////////////////////////////////////////////////////
// Type definitions //
/////////////////////////////////////////////////////////////////////////////
typedef Aggregation<vSpinColourVector, vTComplex, nbasis> Aggregates;
typedef CoarsenedMatrix<vSpinColourVector, vTComplex, nbasis> CoarseDiracMatrix;
typedef CoarseDiracMatrix::CoarseVector CoarseVector;
/////////////////////////////////////////////////////////////////////////////
// Setup of Aggregation //
/////////////////////////////////////////////////////////////////////////////
Aggregates Aggs(Grid_c, Grid_f, 0);
{
LatticeFermion tmp(Aggs.subspace[0].Grid());
for(int n = 0; n < nb; n++) {
gaussian(pRNG_f, Aggs.subspace[n]);
G5C(tmp, Aggs.subspace[n]);
axpby(Aggs.subspace[n + nb], 0.5, -0.5, Aggs.subspace[n], tmp);
axpby(Aggs.subspace[n], 0.5, 0.5, Aggs.subspace[n], tmp);
}
}
/////////////////////////////////////////////////////////////////////////////
// Setup of CoarsenedMatrix and Operator //
/////////////////////////////////////////////////////////////////////////////
const int hermitian = 0;
CoarseDiracMatrix Dc(*Grid_c, *RBGrid_c, hermitian);
Dc.CoarsenOperator(Grid_f, MdagMOp_Dwc, Aggs);
MdagMLinearOperator<CoarseDiracMatrix, CoarseVector> MdagMOp_Dc(Dc);
/////////////////////////////////////////////////////////////////////////////
// Setup vectors used in all tests //
/////////////////////////////////////////////////////////////////////////////
CoarseVector src(Grid_c); random(pRNG_c, src);
CoarseVector diff(Grid_c); diff = Zero();
/////////////////////////////////////////////////////////////////////////////
// Start of tests //
/////////////////////////////////////////////////////////////////////////////
{
std::cout << GridLogMessage << "===========================================================================" << std::endl;
std::cout << GridLogMessage << "= Test Dhop + Mdiag = Munprec" << std::endl;
std::cout << GridLogMessage << "===========================================================================" << std::endl;
CoarseVector phi(Grid_c); phi = Zero();
CoarseVector chi(Grid_c); chi = Zero();
CoarseVector res(Grid_c); res = Zero();
CoarseVector ref(Grid_c); ref = Zero();
Dc.Mdiag(src, phi); std::cout << GridLogMessage << "Applied Mdiag" << std::endl;
Dc.Dhop(src, chi, DaggerNo); std::cout << GridLogMessage << "Applied Dhop" << std::endl;
Dc.M(src, ref); std::cout << GridLogMessage << "Applied M" << std::endl;
res = phi + chi;
diff = ref - res;
auto absDev = norm2(diff);
auto relDev = absDev / norm2(ref);
std::cout << GridLogMessage << "norm2(Munprec), norm2(Dhop + Mdiag), abs. deviation, rel. deviation: "
<< norm2(ref) << " " << norm2(res) << " " << absDev << " " << relDev << " -> check "
<< ((relDev < checkTolerance) ? "passed" : "failed") << std::endl;
assert(relDev <= checkTolerance);
}
{
std::cout << GridLogMessage << "===========================================================================" << std::endl;
std::cout << GridLogMessage << "= Test Meo + Moe = Dhop" << std::endl;
std::cout << GridLogMessage << "===========================================================================" << std::endl;
CoarseVector src_e(RBGrid_c); src_e = Zero();
CoarseVector src_o(RBGrid_c); src_o = Zero();
CoarseVector res_e(RBGrid_c); res_e = Zero();
CoarseVector res_o(RBGrid_c); res_o = Zero();
CoarseVector res(Grid_c); res = Zero();
CoarseVector ref(Grid_c); ref = Zero();
pickCheckerboard(Even, src_e, src);
pickCheckerboard(Odd, src_o, src);
Dc.Meooe(src_e, res_o); std::cout << GridLogMessage << "Applied Meo" << std::endl;
Dc.Meooe(src_o, res_e); std::cout << GridLogMessage << "Applied Moe" << std::endl;
Dc.Dhop(src, ref, DaggerNo); std::cout << GridLogMessage << "Applied Dhop" << std::endl;
setCheckerboard(res, res_o);
setCheckerboard(res, res_e);
diff = ref - res;
auto absDev = norm2(diff);
auto relDev = absDev / norm2(ref);
std::cout << GridLogMessage << "norm2(Dhop), norm2(Meo + Moe), abs. deviation, rel. deviation: "
<< norm2(ref) << " " << norm2(res) << " " << absDev << " " << relDev
<< " -> check " << ((relDev < checkTolerance) ? "passed" : "failed") << std::endl;
assert(relDev <= checkTolerance);
}
{
std::cout << GridLogMessage << "===========================================================================" << std::endl;
std::cout << GridLogMessage << "= Test |(Im(v^dag M^dag M v)| = 0" << std::endl;
std::cout << GridLogMessage << "===========================================================================" << std::endl;
CoarseVector tmp(Grid_c); tmp = Zero();
CoarseVector phi(Grid_c); phi = Zero();
Dc.M(src, tmp); std::cout << GridLogMessage << "Applied M" << std::endl;
Dc.Mdag(tmp, phi); std::cout << GridLogMessage << "Applied Mdag" << std::endl;
std::cout << GridLogMessage << "src = " << norm2(src) << " tmp = " << norm2(tmp) << " phi = " << norm2(phi) << std::endl;
ComplexD dot = innerProduct(src, phi);
auto relDev = abs(imag(dot)) / abs(real(dot));
std::cout << GridLogMessage << "Re(v^dag M^dag M v), Im(v^dag M^dag M v), rel.deviation: "
<< real(dot) << " " << imag(dot) << " " << relDev
<< " -> check " << ((relDev < checkTolerance) ? "passed" : "failed") << std::endl;
assert(relDev <= checkTolerance);
}
{
std::cout << GridLogMessage << "===========================================================================" << std::endl;
std::cout << GridLogMessage << "= Test |(Im(v^dag Mooee^dag Mooee v)| = 0 (full grid)" << std::endl;
std::cout << GridLogMessage << "===========================================================================" << std::endl;
CoarseVector tmp(Grid_c); tmp = Zero();
CoarseVector phi(Grid_c); phi = Zero();
Dc.Mooee(src, tmp); std::cout << GridLogMessage << "Applied Mooee" << std::endl;
Dc.MooeeDag(tmp, phi); std::cout << GridLogMessage << "Applied MooeeDag" << std::endl;
ComplexD dot = innerProduct(src, phi);
auto relDev = abs(imag(dot)) / abs(real(dot));
std::cout << GridLogMessage << "Re(v^dag Mooee^dag Mooee v), Im(v^dag Mooee^dag Mooee v), rel.deviation: "
<< real(dot) << " " << imag(dot) << " " << relDev
<< " -> check " << ((relDev < checkTolerance) ? "passed" : "failed") << std::endl;
assert(relDev <= checkTolerance);
}
{
std::cout << GridLogMessage << "===========================================================================" << std::endl;
std::cout << GridLogMessage << "= Test MooeeInv Mooee = 1 (full grid)" << std::endl;
std::cout << GridLogMessage << "===========================================================================" << std::endl;
CoarseVector tmp(Grid_c); tmp = Zero();
CoarseVector phi(Grid_c); phi = Zero();
Dc.Mooee(src, tmp); std::cout << GridLogMessage << "Applied Mooee" << std::endl;
Dc.MooeeInv(tmp, phi); std::cout << GridLogMessage << "Applied MooeeInv" << std::endl;
diff = src - phi;
auto absDev = norm2(diff);
auto relDev = absDev / norm2(src);
std::cout << GridLogMessage << "norm2(src), norm2(MooeeInv Mooee src), abs. deviation, rel. deviation: "
<< norm2(src) << " " << norm2(phi) << " " << absDev << " " << relDev
<< " -> check " << ((relDev < checkTolerance) ? "passed" : "failed") << std::endl;
assert(relDev <= checkTolerance);
}
{
std::cout << GridLogMessage << "===========================================================================" << std::endl;
std::cout << GridLogMessage << "= Test MeooeDagger is the dagger of Meooe by requiring" << std::endl;
std::cout << GridLogMessage << "= < phi | Meooe | chi > * = < chi | Meooe^dag| phi>" << std::endl;
std::cout << GridLogMessage << "===========================================================================" << std::endl;
// clang-format off
CoarseVector phi(Grid_c); random(pRNG_c, phi);
CoarseVector chi(Grid_c); random(pRNG_c, chi);
CoarseVector chi_e(RBGrid_c); chi_e = Zero();
CoarseVector chi_o(RBGrid_c); chi_o = Zero();
CoarseVector dchi_e(RBGrid_c); dchi_e = Zero();
CoarseVector dchi_o(RBGrid_c); dchi_o = Zero();
CoarseVector phi_e(RBGrid_c); phi_e = Zero();
CoarseVector phi_o(RBGrid_c); phi_o = Zero();
CoarseVector dphi_e(RBGrid_c); dphi_e = Zero();
CoarseVector dphi_o(RBGrid_c); dphi_o = Zero();
// clang-format on
pickCheckerboard(Even, chi_e, chi);
pickCheckerboard(Odd, chi_o, chi);
pickCheckerboard(Even, phi_e, phi);
pickCheckerboard(Odd, phi_o, phi);
Dc.Meooe(chi_e, dchi_o); std::cout << GridLogMessage << "Applied Meo" << std::endl;
Dc.Meooe(chi_o, dchi_e); std::cout << GridLogMessage << "Applied Moe" << std::endl;
Dc.MeooeDag(phi_e, dphi_o); std::cout << GridLogMessage << "Applied MeoDag" << std::endl;
Dc.MeooeDag(phi_o, dphi_e); std::cout << GridLogMessage << "Applied MoeDag" << std::endl;
ComplexD phiDchi_e = innerProduct(phi_e, dchi_e);
ComplexD phiDchi_o = innerProduct(phi_o, dchi_o);
ComplexD chiDphi_e = innerProduct(chi_e, dphi_e);
ComplexD chiDphi_o = innerProduct(chi_o, dphi_o);
std::cout << GridLogDebug << "norm dchi_e = " << norm2(dchi_e) << " norm dchi_o = " << norm2(dchi_o) << " norm dphi_e = " << norm2(dphi_e)
<< " norm dphi_o = " << norm2(dphi_e) << std::endl;
std::cout << GridLogMessage << "e " << phiDchi_e << " " << chiDphi_e << std::endl;
std::cout << GridLogMessage << "o " << phiDchi_o << " " << chiDphi_o << std::endl;
std::cout << GridLogMessage << "phiDchi_e - conj(chiDphi_o) " << phiDchi_e - conj(chiDphi_o) << std::endl;
std::cout << GridLogMessage << "phiDchi_o - conj(chiDphi_e) " << phiDchi_o - conj(chiDphi_e) << std::endl;
}
{
std::cout << GridLogMessage << "===========================================================================" << std::endl;
std::cout << GridLogMessage << "= Test MooeeInv Mooee = 1 (checkerboards separately)" << std::endl;
std::cout << GridLogMessage << "===========================================================================" << std::endl;
CoarseVector chi(Grid_c); random(pRNG_c, chi);
CoarseVector tmp(Grid_c); tmp = Zero();
CoarseVector phi(Grid_c); phi = Zero();
CoarseVector chi_e(RBGrid_c); chi_e = Zero();
CoarseVector chi_o(RBGrid_c); chi_o = Zero();
CoarseVector phi_e(RBGrid_c); phi_e = Zero();
CoarseVector phi_o(RBGrid_c); phi_o = Zero();
CoarseVector tmp_e(RBGrid_c); tmp_e = Zero();
CoarseVector tmp_o(RBGrid_c); tmp_o = Zero();
pickCheckerboard(Even, chi_e, chi);
pickCheckerboard(Odd, chi_o, chi);
pickCheckerboard(Even, tmp_e, tmp);
pickCheckerboard(Odd, tmp_o, tmp);
Dc.Mooee(chi_e, tmp_e); std::cout << GridLogMessage << "Applied Mee" << std::endl;
Dc.MooeeInv(tmp_e, phi_e); std::cout << GridLogMessage << "Applied MeeInv" << std::endl;
Dc.Mooee(chi_o, tmp_o); std::cout << GridLogMessage << "Applied Moo" << std::endl;
Dc.MooeeInv(tmp_o, phi_o); std::cout << GridLogMessage << "Applied MooInv" << std::endl;
setCheckerboard(phi, phi_e);
setCheckerboard(phi, phi_o);
diff = chi - phi;
auto absDev = norm2(diff);
auto relDev = absDev / norm2(chi);
std::cout << GridLogMessage << "norm2(chi), norm2(MeeInv Mee chi), abs. deviation, rel. deviation: "
<< norm2(chi) << " " << norm2(phi) << " " << absDev << " " << relDev
<< " -> check " << ((relDev < checkTolerance) ? "passed" : "failed") << std::endl;
assert(relDev <= checkTolerance);
}
{
std::cout << GridLogMessage << "===========================================================================" << std::endl;
std::cout << GridLogMessage << "= Test MooeeDag MooeeInvDag = 1 (checkerboards separately)" << std::endl;
std::cout << GridLogMessage << "===========================================================================" << std::endl;
CoarseVector chi(Grid_c); random(pRNG_c, chi);
CoarseVector tmp(Grid_c); tmp = Zero();
CoarseVector phi(Grid_c); phi = Zero();
CoarseVector chi_e(RBGrid_c); chi_e = Zero();
CoarseVector chi_o(RBGrid_c); chi_o = Zero();
CoarseVector phi_e(RBGrid_c); phi_e = Zero();
CoarseVector phi_o(RBGrid_c); phi_o = Zero();
CoarseVector tmp_e(RBGrid_c); tmp_e = Zero();
CoarseVector tmp_o(RBGrid_c); tmp_o = Zero();
pickCheckerboard(Even, chi_e, chi);
pickCheckerboard(Odd, chi_o, chi);
pickCheckerboard(Even, tmp_e, tmp);
pickCheckerboard(Odd, tmp_o, tmp);
Dc.MooeeDag(chi_e, tmp_e); std::cout << GridLogMessage << "Applied MeeDag" << std::endl;
Dc.MooeeInvDag(tmp_e, phi_e); std::cout << GridLogMessage << "Applied MeeInvDag" << std::endl;
Dc.MooeeDag(chi_o, tmp_o); std::cout << GridLogMessage << "Applied MooDag" << std::endl;
Dc.MooeeInvDag(tmp_o, phi_o); std::cout << GridLogMessage << "Applied MooInvDag" << std::endl;
setCheckerboard(phi, phi_e);
setCheckerboard(phi, phi_o);
diff = chi - phi;
auto absDev = norm2(diff);
auto relDev = absDev / norm2(chi);
std::cout << GridLogMessage << "norm2(chi), norm2(MeeDag MeeInvDag chi), abs. deviation, rel. deviation: "
<< norm2(chi) << " " << norm2(phi) << " " << absDev << " " << relDev
<< " -> check " << ((relDev < checkTolerance) ? "passed" : "failed") << std::endl;
assert(relDev <= checkTolerance);
}
{
std::cout << GridLogMessage << "===========================================================================" << std::endl;
std::cout << GridLogMessage << "= Test Meo + Moe + Moo + Mee = Munprec" << std::endl;
std::cout << GridLogMessage << "===========================================================================" << std::endl;
CoarseVector chi(Grid_c); chi = Zero();
CoarseVector phi(Grid_c); phi = Zero();
CoarseVector ref(Grid_c); ref = Zero();
CoarseVector src_e(RBGrid_c); src_e = Zero();
CoarseVector src_o(RBGrid_c); src_o = Zero();
CoarseVector phi_e(RBGrid_c); phi_e = Zero();
CoarseVector phi_o(RBGrid_c); phi_o = Zero();
CoarseVector chi_e(RBGrid_c); chi_e = Zero();
CoarseVector chi_o(RBGrid_c); chi_o = Zero();
pickCheckerboard(Even, src_e, src);
pickCheckerboard(Odd, src_o, src);
pickCheckerboard(Even, phi_e, phi);
pickCheckerboard(Odd, phi_o, phi);
pickCheckerboard(Even, chi_e, chi);
pickCheckerboard(Odd, chi_o, chi);
// M phi = (Mooee src_e + Meooe src_o , Mooee src_o + Meooe src_e)
Dc.M(src, ref); // Reference result from the unpreconditioned operator
// EO matrix
Dc.Mooee(src_e, chi_e); std::cout << GridLogMessage << "Applied Mee" << std::endl;
Dc.Mooee(src_o, chi_o); std::cout << GridLogMessage << "Applied Moo" << std::endl;
Dc.Meooe(src_o, phi_e); std::cout << GridLogMessage << "Applied Moe" << std::endl;
Dc.Meooe(src_e, phi_o); std::cout << GridLogMessage << "Applied Meo" << std::endl;
phi_o += chi_o;
phi_e += chi_e;
setCheckerboard(phi, phi_e);
setCheckerboard(phi, phi_o);
std::cout << GridLogDebug << "norm phi_e = " << norm2(phi_e) << " norm phi_o = " << norm2(phi_o) << " norm phi = " << norm2(phi) << std::endl;
diff = ref - phi;
auto absDev = norm2(diff);
auto relDev = absDev / norm2(ref);
std::cout << GridLogMessage << "norm2(Dunprec), norm2(Deoprec), abs. deviation, rel. deviation: "
<< norm2(ref) << " " << norm2(phi) << " " << absDev << " " << relDev
<< " -> check " << ((relDev < checkTolerance) ? "passed" : "failed") << std::endl;
assert(relDev <= checkTolerance);
}
{
std::cout << GridLogMessage << "===========================================================================" << std::endl;
std::cout << GridLogMessage << "= Test MpcDagMpc is hermitian" << std::endl;
std::cout << GridLogMessage << "===========================================================================" << std::endl;
CoarseVector phi(Grid_c); random(pRNG_c, phi);
CoarseVector chi(Grid_c); random(pRNG_c, chi);
CoarseVector chi_e(RBGrid_c); chi_e = Zero();
CoarseVector chi_o(RBGrid_c); chi_o = Zero();
CoarseVector dchi_e(RBGrid_c); dchi_e = Zero();
CoarseVector dchi_o(RBGrid_c); dchi_o = Zero();
CoarseVector phi_e(RBGrid_c); phi_e = Zero();
CoarseVector phi_o(RBGrid_c); phi_o = Zero();
CoarseVector dphi_e(RBGrid_c); dphi_e = Zero();
CoarseVector dphi_o(RBGrid_c); dphi_o = Zero();
pickCheckerboard(Even, chi_e, chi);
pickCheckerboard(Odd, chi_o, chi);
pickCheckerboard(Even, phi_e, phi);
pickCheckerboard(Odd, phi_o, phi);
SchurDiagMooeeOperator<CoarseDiracMatrix,CoarseVector> HermOpEO(Dc);
HermOpEO.MpcDagMpc(chi_e, dchi_e); std::cout << GridLogMessage << "Applied MpcDagMpc to chi_e" << std::endl;
HermOpEO.MpcDagMpc(chi_o, dchi_o); std::cout << GridLogMessage << "Applied MpcDagMpc to chi_o" << std::endl;
HermOpEO.MpcDagMpc(phi_e, dphi_e); std::cout << GridLogMessage << "Applied MpcDagMpc to phi_e" << std::endl;
HermOpEO.MpcDagMpc(phi_o, dphi_o); std::cout << GridLogMessage << "Applied MpcDagMpc to phi_o" << std::endl;
ComplexD phiDchi_e = innerProduct(phi_e, dchi_e);
ComplexD phiDchi_o = innerProduct(phi_o, dchi_o);
ComplexD chiDphi_e = innerProduct(chi_e, dphi_e);
ComplexD chiDphi_o = innerProduct(chi_o, dphi_o);
std::cout << GridLogMessage << "e " << phiDchi_e << " " << chiDphi_e << std::endl;
std::cout << GridLogMessage << "o " << phiDchi_o << " " << chiDphi_o << std::endl;
std::cout << GridLogMessage << "phiDchi_e - conj(chiDphi_e) " << phiDchi_e - conj(chiDphi_e) << std::endl;
std::cout << GridLogMessage << "phiDchi_o - conj(chiDphi_o) " << phiDchi_o - conj(chiDphi_o) << std::endl;
}
Grid_finalize();
}

2
tests/solver/Test_contfrac_cg.cc
View File

@@ -94,7 +94,7 @@ int main (int argc, char ** argv)
GridParallelRNG RNG5(FGrid); RNG5.SeedFixedIntegers(seeds5);
GridParallelRNG RNG4(UGrid); RNG4.SeedFixedIntegers(seeds4);
LatticeGaugeField Umu(UGrid); SU3::HotConfiguration(RNG4,Umu);
LatticeGaugeField Umu(UGrid); SU<Nc>::HotConfiguration(RNG4,Umu);
std::vector<LatticeColourMatrix> U(4,UGrid);
RealD mass=0.1;

2
tests/solver/Test_dwf_cg_prec.cc
View File

@@ -67,7 +67,7 @@ int main(int argc, char** argv) {
result = Zero();
LatticeGaugeField Umu(UGrid);
SU3::HotConfiguration(RNG4, Umu);
SU<Nc>::HotConfiguration(RNG4, Umu);
std::cout << GridLogMessage << "Lattice dimensions: " << GridDefaultLatt()
<< " Ls: " << Ls << std::endl;

2
tests/solver/Test_dwf_cg_schur.cc
View File

@@ -61,7 +61,7 @@ int main (int argc, char ** argv)
LatticeFermion src(FGrid); random(RNG5,src);
LatticeFermion result(FGrid); result=Zero();
LatticeGaugeField Umu(UGrid); SU3::HotConfiguration(RNG4,Umu);
LatticeGaugeField Umu(UGrid); SU<Nc>::HotConfiguration(RNG4,Umu);
std::vector<LatticeColourMatrix> U(4,UGrid);
for(int mu=0;mu<Nd;mu++){

2
tests/solver/Test_dwf_cg_unprec.cc
View File

@@ -61,7 +61,7 @@ int main (int argc, char ** argv)
LatticeFermion src(FGrid); random(RNG5,src);
LatticeFermion result(FGrid); result=Zero();
LatticeGaugeField Umu(UGrid); SU3::HotConfiguration(RNG4,Umu);
LatticeGaugeField Umu(UGrid); SU<Nc>::HotConfiguration(RNG4,Umu);
std::vector<LatticeColourMatrix> U(4,UGrid);
for(int mu=0;mu<Nd;mu++){

2
tests/solver/Test_dwf_cr_unprec.cc
View File

@@ -65,7 +65,7 @@ int main (int argc, char ** argv)
LatticeFermion src(FGrid); random(RNG5,src);
LatticeFermion result(FGrid); result=Zero();
LatticeGaugeField Umu(UGrid); SU3::HotConfiguration(RNG4,Umu);
LatticeGaugeField Umu(UGrid); SU<Nc>::HotConfiguration(RNG4,Umu);
std::vector<LatticeColourMatrix> U(4,UGrid);

2
tests/solver/Test_dwf_fpgcr.cc
View File

@@ -68,7 +68,7 @@ int main (int argc, char ** argv)
LatticeFermion result(FGrid); result=Zero();
LatticeGaugeField Umu(UGrid);
SU3::HotConfiguration(RNG4,Umu);
SU<Nc>::HotConfiguration(RNG4,Umu);
ConjugateResidual<LatticeFermion> CR(1.0e-6,10000);

15
tests/solver/Test_dwf_hdcr.cc
View File

@@ -222,9 +222,16 @@ int main (int argc, char ** argv)
GridCartesian *Coarse4d = SpaceTimeGrid::makeFourDimGrid(clatt, GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi());;
GridCartesian *Coarse5d = SpaceTimeGrid::makeFiveDimGrid(1,Coarse4d);
GridCartesian *CoarseCoarse4d = SpaceTimeGrid::makeFourDimGrid(cclatt, GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi());;
GridCartesian *CoarseCoarse4d = SpaceTimeGrid::makeFourDimGrid(cclatt, GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi());
GridCartesian *CoarseCoarse5d = SpaceTimeGrid::makeFiveDimGrid(1,CoarseCoarse4d);
GridRedBlackCartesian * Coarse4dRB = SpaceTimeGrid::makeFourDimRedBlackGrid(Coarse4d);
GridRedBlackCartesian * Coarse5dRB = SpaceTimeGrid::makeFiveDimRedBlackGrid(1,Coarse4d);
GridRedBlackCartesian *CoarseCoarse4dRB = SpaceTimeGrid::makeFourDimRedBlackGrid(CoarseCoarse4d);
GridRedBlackCartesian *CoarseCoarse5dRB = SpaceTimeGrid::makeFiveDimRedBlackGrid(1,CoarseCoarse4d);
std::vector<int> seeds4({1,2,3,4});
std::vector<int> seeds5({5,6,7,8});
std::vector<int> cseeds({5,6,7,8});
@@ -282,8 +289,7 @@ int main (int argc, char ** argv)
Gamma5R5HermitianLinearOperator<DomainWallFermionR,LatticeFermion> HermIndefOp(Ddwf);
Level1Op LDOp(*Coarse5d,1); LDOp.CoarsenOperator(FGrid,HermIndefOp,Aggregates);
Level1Op LDOp(*Coarse5d,*Coarse5dRB,1); LDOp.CoarsenOperator(FGrid,HermIndefOp,Aggregates);
//////////////////////////////////////////////////
// Deflate the course space. Recursive multigrid?
@@ -311,12 +317,11 @@ int main (int argc, char ** argv)
}
}
Level2Op L2Op(*CoarseCoarse5d,1); // Hermitian matrix
Level2Op L2Op(*CoarseCoarse5d,*CoarseCoarse5dRB,1); // Hermitian matrix
typedef Level2Op::CoarseVector CoarseCoarseVector;
HermitianLinearOperator<Level1Op,CoarseVector> L1LinOp(LDOp);
L2Op.CoarsenOperator(Coarse5d,L1LinOp,CoarseAggregates);
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
std::cout<<GridLogMessage << " Running CoarseCoarse grid Lanczos "<< std::endl;
std::cout<<GridLogMessage << "**************************************************"<< std::endl;

397
tests/solver/Test_dwf_hdcr_16_rb.cc Normal file
View File

@@ -0,0 +1,397 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: ./tests/Test_dwf_hdcr.cc
Copyright (C) 2015
Author: Antonin Portelli <antonin.portelli@me.com>
Author: Peter Boyle <paboyle@ph.ed.ac.uk>
Author: paboyle <paboyle@ph.ed.ac.uk>
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License along
with this program; if not, write to the Free Software Foundation, Inc.,
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
See the full license in the file "LICENSE" in the top level distribution directory
*************************************************************************************/
/* END LEGAL */
#include <Grid/Grid.h>
#include <Grid/algorithms/iterative/PrecGeneralisedConjugateResidual.h>
using namespace std;
using namespace Grid;
/* Params
* Grid:
* block1(4)
* block2(4)
*
* Subspace
* * Fine : Subspace(nbasis,hi,lo,order,first,step) -- 32, 60,0.02,500,100,100
* * Coarse: Subspace(nbasis,hi,lo,order,first,step) -- 32, 18,0.02,500,100,100
* Smoother:
* * Fine: Cheby(hi, lo, order) -- 60,0.5,10
* * Coarse: Cheby(hi, lo, order) -- 12,0.1,4
* Lanczos:
* CoarseCoarse IRL( Nk, Nm, Nstop, poly(lo,hi,order)) 24,36,24,0.002,4.0,61
*/
RealD InverseApproximation(RealD x){
return 1.0/x;
}
template<class Field> class SolverWrapper : public LinearFunction<Field> {
private:
CheckerBoardedSparseMatrixBase<Field> & _Matrix;
SchurRedBlackBase<Field> & _Solver;
public:
/////////////////////////////////////////////////////
// Wrap the usual normal equations trick
/////////////////////////////////////////////////////
SolverWrapper(CheckerBoardedSparseMatrixBase<Field> &Matrix,
SchurRedBlackBase<Field> &Solver)
: _Matrix(Matrix), _Solver(Solver) {};
void operator() (const Field &in, Field &out){
_Solver(_Matrix,in,out); // Mdag M out = Mdag in
}
};
template<class Field,class Matrix> class ChebyshevSmoother : public LinearFunction<Field>
{
public:
typedef LinearOperatorBase<Field> FineOperator;
Matrix & _SmootherMatrix;
FineOperator & _SmootherOperator;
Chebyshev<Field> Cheby;
ChebyshevSmoother(RealD _lo,RealD _hi,int _ord, FineOperator &SmootherOperator,Matrix &SmootherMatrix) :
_SmootherOperator(SmootherOperator),
_SmootherMatrix(SmootherMatrix),
Cheby(_lo,_hi,_ord,InverseApproximation)
{};
void operator() (const Field &in, Field &out)
{
Field tmp(in.Grid());
MdagMLinearOperator<Matrix,Field> MdagMOp(_SmootherMatrix);
_SmootherOperator.AdjOp(in,tmp);
Cheby(MdagMOp,tmp,out);
}
};
template<class Field,class Matrix> class MirsSmoother : public LinearFunction<Field>
{
public:
typedef LinearOperatorBase<Field> FineOperator;
Matrix & SmootherMatrix;
FineOperator & SmootherOperator;
RealD tol;
RealD shift;
int maxit;
MirsSmoother(RealD _shift,RealD _tol,int _maxit,FineOperator &_SmootherOperator,Matrix &_SmootherMatrix) :
shift(_shift),tol(_tol),maxit(_maxit),
SmootherOperator(_SmootherOperator),
SmootherMatrix(_SmootherMatrix)
{};
void operator() (const Field &in, Field &out)
{
ZeroGuesser<Field> Guess;
ConjugateGradient<Field> CG(tol,maxit,false);
Field src(in.Grid());
ShiftedMdagMLinearOperator<SparseMatrixBase<Field>,Field> MdagMOp(SmootherMatrix,shift);
SmootherOperator.AdjOp(in,src);
Guess(src,out);
CG(MdagMOp,src,out);
}
};
template<class Fobj,class CComplex,int nbasis, class Matrix, class Guesser, class CoarseSolver>
class MultiGridPreconditioner : public LinearFunction< Lattice<Fobj> > {
public:
typedef Aggregation<Fobj,CComplex,nbasis> Aggregates;
typedef CoarsenedMatrix<Fobj,CComplex,nbasis> CoarseOperator;
typedef typename Aggregation<Fobj,CComplex,nbasis>::CoarseVector CoarseVector;
typedef typename Aggregation<Fobj,CComplex,nbasis>::CoarseMatrix CoarseMatrix;
typedef typename Aggregation<Fobj,CComplex,nbasis>::FineField FineField;
typedef LinearOperatorBase<FineField> FineOperator;
typedef LinearFunction <FineField> FineSmoother;
Aggregates & _Aggregates;
CoarseOperator & _CoarseOperator;
Matrix & _FineMatrix;
FineOperator & _FineOperator;
Guesser & _Guess;
FineSmoother & _Smoother;
CoarseSolver & _CoarseSolve;
int level; void Level(int lv) {level = lv; };
#define GridLogLevel std::cout << GridLogMessage <<std::string(level,'\t')<< " Level "<<level <<" "
MultiGridPreconditioner(Aggregates &Agg, CoarseOperator &Coarse,
FineOperator &Fine,Matrix &FineMatrix,
FineSmoother &Smoother,
Guesser &Guess_,
CoarseSolver &CoarseSolve_)
: _Aggregates(Agg),
_CoarseOperator(Coarse),
_FineOperator(Fine),
_FineMatrix(FineMatrix),
_Smoother(Smoother),
_Guess(Guess_),
_CoarseSolve(CoarseSolve_),
level(1) { }
virtual void operator()(const FineField &in, FineField & out)
{
CoarseVector Csrc(_CoarseOperator.Grid());
CoarseVector Csol(_CoarseOperator.Grid());
FineField vec1(in.Grid());
FineField vec2(in.Grid());
double t;
// Fine Smoother
t=-usecond();
_Smoother(in,out);
t+=usecond();
GridLogLevel << "Smoother took "<< t/1000.0<< "ms" <<std::endl;
// Update the residual
_FineOperator.Op(out,vec1); sub(vec1, in ,vec1);
// Fine to Coarse
t=-usecond();
_Aggregates.ProjectToSubspace (Csrc,vec1);
t+=usecond();
GridLogLevel << "Project to coarse took "<< t/1000.0<< "ms" <<std::endl;
// Coarse correction
t=-usecond();
_CoarseSolve(Csrc,Csol);
t+=usecond();
GridLogLevel << "Coarse solve took "<< t/1000.0<< "ms" <<std::endl;
// Coarse to Fine
t=-usecond();
_Aggregates.PromoteFromSubspace(Csol,vec1);
add(out,out,vec1);
t+=usecond();
GridLogLevel << "Promote to this level took "<< t/1000.0<< "ms" <<std::endl;
// Residual
_FineOperator.Op(out,vec1); sub(vec1 ,in , vec1);
// Fine Smoother
t=-usecond();
_Smoother(vec1,vec2);
t+=usecond();
GridLogLevel << "Smoother took "<< t/1000.0<< "ms" <<std::endl;
add( out,out,vec2);
}
};
int main (int argc, char ** argv)
{
Grid_init(&argc,&argv);
const int Ls=16;
GridCartesian * UGrid = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi());
GridRedBlackCartesian * UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
GridCartesian * FGrid = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid);
GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid);
///////////////////////////////////////////////////
// Construct a coarsened grid; utility for this?
///////////////////////////////////////////////////
std::vector<int> block ({2,2,2,2});
std::vector<int> blockc ({2,2,2,2});
const int nbasis= 32;
const int nbasisc= 32;
auto clatt = GridDefaultLatt();
for(int d=0;d<clatt.size();d++){
clatt[d] = clatt[d]/block[d];
}
auto cclatt = clatt;
for(int d=0;d<clatt.size();d++){
cclatt[d] = clatt[d]/blockc[d];
}
GridCartesian *Coarse4d = SpaceTimeGrid::makeFourDimGrid(clatt, GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi());;
GridCartesian *Coarse5d = SpaceTimeGrid::makeFiveDimGrid(1,Coarse4d);
// GridCartesian *CoarseCoarse4d = SpaceTimeGrid::makeFourDimGrid(cclatt, GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi());;
// GridCartesian *CoarseCoarse5d = SpaceTimeGrid::makeFiveDimGrid(1,CoarseCoarse4d);
std::vector<int> seeds4({1,2,3,4});
std::vector<int> seeds5({5,6,7,8});
std::vector<int> cseeds({5,6,7,8});
GridParallelRNG RNG5(FGrid); RNG5.SeedFixedIntegers(seeds5);
GridParallelRNG RNG4(UGrid); RNG4.SeedFixedIntegers(seeds4);
GridParallelRNG CRNG(Coarse5d);CRNG.SeedFixedIntegers(cseeds);
LatticeFermion src(FGrid); gaussian(RNG5,src);// src=src+g5*src;
LatticeFermion result(FGrid);
LatticeGaugeField Umu(UGrid);
FieldMetaData header;
std::string file("./ckpoint_lat.4000");
//std::string file("./ckpoint_lat.1000");
NerscIO::readConfiguration(Umu,header,file);
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
std::cout<<GridLogMessage << "Building g5R5 hermitian DWF operator" <<std::endl;
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
RealD mass=0.001;
RealD M5=1.8;
DomainWallFermionR Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5);
typedef Aggregation<vSpinColourVector,vTComplex,nbasis> Subspace;
typedef CoarsenedMatrix<vSpinColourVector,vTComplex,nbasis> CoarseOperator;
typedef CoarseOperator::CoarseVector CoarseVector;
typedef CoarseOperator::siteVector siteVector;
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
std::cout<<GridLogMessage << "Calling Aggregation class to build subspace" <<std::endl;
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
MdagMLinearOperator<DomainWallFermionR,LatticeFermion> HermDefOp(Ddwf);
Subspace Aggregates(Coarse5d,FGrid,0);
assert ( (nbasis & 0x1)==0);
{
int nb=nbasis/2;
LatticeFermion A(FGrid);
LatticeFermion B(FGrid);
// Aggregates.CreateSubspaceChebyshev(RNG5,HermDefOp,nb,60.0,0.002,1000,800,100,0.0);
// Aggregates.CreateSubspaceChebyshev(RNG5,HermDefOp,nb,60.0,0.02,1000,800,100,0.0);
Aggregates.CreateSubspaceChebyshev(RNG5,HermDefOp,nb,60.0,0.05,500,200,150,0.0);//
// Aggregates.CreateSubspaceChebyshev(RNG5,HermDefOp,nb,60.0,0.01,1000,100,100,0.0); // Slightly faster
for(int n=0;n<nb;n++){
std::cout << GridLogMessage << " G5R5 "<<n<<std::endl;
G5R5(Aggregates.subspace[n+nb],Aggregates.subspace[n]);
std::cout << GridLogMessage << " Projection "<<n<<std::endl;
A = Aggregates.subspace[n];
B = Aggregates.subspace[n+nb];
std::cout << GridLogMessage << " Copy "<<n<<std::endl;
Aggregates.subspace[n] = A+B; // 1+G5 // eigen value of G5R5 is +1
std::cout << GridLogMessage << " P+ "<<n<<std::endl;
Aggregates.subspace[n+nb]= A-B; // 1-G5 // eigen value of G5R5 is -1
std::cout << GridLogMessage << " P- "<<n<<std::endl;
}
}
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
std::cout<<GridLogMessage << "Building coarse representation of Indef operator" <<std::endl;
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
typedef CoarsenedMatrix<vSpinColourVector,vTComplex,nbasis> Level1Op;
typedef CoarsenedMatrix<siteVector,iScalar<vTComplex>,nbasisc> Level2Op;
Gamma5R5HermitianLinearOperator<DomainWallFermionR,LatticeFermion> HermIndefOp(Ddwf);
GridRedBlackCartesian * Coarse4dRB = SpaceTimeGrid::makeFourDimRedBlackGrid(Coarse4d);
std::cout << " Making 5D coarse RB grid " <<std::endl;
GridRedBlackCartesian * Coarse5dRB = SpaceTimeGrid::makeFiveDimRedBlackGrid(1,Coarse4d);
std::cout << " Made 5D coarse RB grid " <<std::endl;
Level1Op LDOp(*Coarse5d,*Coarse5dRB,1); LDOp.CoarsenOperator(FGrid,HermIndefOp,Aggregates);
//////////////////////////////////////////////////
// Deflate the course space. Recursive multigrid?
//////////////////////////////////////////////////
typedef Aggregation<siteVector,iScalar<vTComplex>,nbasisc> CoarseSubspace;
// CoarseSubspace CoarseAggregates(CoarseCoarse5d,Coarse5d,0);
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
std::cout<<GridLogMessage << "Build deflation space in coarse operator "<< std::endl;
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
MdagMLinearOperator<CoarseOperator,CoarseVector> PosdefLdop(LDOp);
typedef Level2Op::CoarseVector CoarseCoarseVector;
CoarseVector c_src(Coarse5d); c_src=1.0;
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
std::cout<<GridLogMessage << "Building 3 level Multigrid "<< std::endl;
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
typedef MultiGridPreconditioner<vSpinColourVector, vTComplex,nbasis, DomainWallFermionR,ZeroGuesser<CoarseVector> , SolverWrapper<CoarseVector> > TwoLevelMG;
typedef MultiGridPreconditioner<siteVector,iScalar<vTComplex>,nbasisc,Level1Op, DeflatedGuesser<CoarseCoarseVector>, NormalEquations<CoarseCoarseVector> > CoarseMG;
typedef MultiGridPreconditioner<vSpinColourVector, vTComplex,nbasis, DomainWallFermionR,ZeroGuesser<CoarseVector>, LinearFunction<CoarseVector> > ThreeLevelMG;
ChebyshevSmoother<LatticeFermion,DomainWallFermionR> FineSmoother(0.5,60.0,12,HermIndefOp,Ddwf);
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
std::cout<<GridLogMessage << "Calling 2 level Multigrid "<< std::endl;
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
result=Zero();
ZeroGuesser<CoarseVector> CoarseZeroGuesser;
ConjugateGradient<CoarseVector> CoarseCG(0.005,1000);
// SchurDiagMooeeOperator<CoarseOperator,CoarseVector> CoarseMpcDagMpc(LDOp);
SchurRedBlackDiagMooeeSolve<CoarseVector> CoarseRBCG(CoarseCG);
SolverWrapper<CoarseVector> CoarseSolver(LDOp,CoarseRBCG);
// NormalEquations<CoarseVector> CoarseCGNE(LDOp,CoarseCG,CoarseZeroGuesser);
TwoLevelMG TwoLevelPrecon(Aggregates, LDOp,
HermIndefOp,Ddwf,
FineSmoother,
CoarseZeroGuesser,
CoarseSolver);
TwoLevelPrecon.Level(1);
PrecGeneralisedConjugateResidual<LatticeFermion> l1PGCR(1.0e-8,20,HermIndefOp,TwoLevelPrecon,16,16);
l1PGCR.Level(1);
l1PGCR(src,result);
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
std::cout<<GridLogMessage << "Calling CG "<< std::endl;
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
ConjugateGradient<LatticeFermion> pCG(1.0e-8,60000);
result=Zero();
// pCG(HermDefOp,src,result);
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
std::cout<<GridLogMessage << "Calling red black CG "<< std::endl;
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
result=Zero();
LatticeFermion src_o(FrbGrid);
LatticeFermion result_o(FrbGrid);
pickCheckerboard(Odd,src_o,src);
result_o=Zero();
SchurDiagMooeeOperator<DomainWallFermionR,LatticeFermion> HermOpEO(Ddwf);
// pCG(HermOpEO,src_o,result_o);
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
std::cout<<GridLogMessage << " Fine PowerMethod "<< std::endl;
PowerMethod<LatticeFermion> PM; PM(HermDefOp,src);
std::cout<<GridLogMessage << " Coarse PowerMethod "<< std::endl;
PowerMethod<CoarseVector> cPM; cPM(PosdefLdop,c_src);
// std::cout<<GridLogMessage << " CoarseCoarse PowerMethod "<< std::endl;
// PowerMethod<CoarseCoarseVector> ccPM; ccPM(IRLHermOpL2,cc_src);
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
std::cout<<GridLogMessage << "Done "<< std::endl;
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
Grid_finalize();
}

477
tests/solver/Test_dwf_hdcr_24_regression.cc Normal file
View File

@@ -0,0 +1,477 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: ./tests/Test_dwf_hdcr.cc
Copyright (C) 2015
Author: Antonin Portelli <antonin.portelli@me.com>
Author: Peter Boyle <paboyle@ph.ed.ac.uk>
Author: paboyle <paboyle@ph.ed.ac.uk>
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License along
with this program; if not, write to the Free Software Foundation, Inc.,
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
See the full license in the file "LICENSE" in the top level distribution directory
*************************************************************************************/
/* END LEGAL */
#include <Grid/Grid.h>
#include <Grid/algorithms/iterative/PrecGeneralisedConjugateResidual.h>
using namespace std;
using namespace Grid;
/* Params
* Grid:
* block1(4)
* block2(4)
*
* Subspace
* * Fine : Subspace(nbasis,hi,lo,order,first,step) -- 32, 60,0.02,500,100,100
* * Coarse: Subspace(nbasis,hi,lo,order,first,step) -- 32, 18,0.02,500,100,100
* Smoother:
* * Fine: Cheby(hi, lo, order) -- 60,0.5,10
* * Coarse: Cheby(hi, lo, order) -- 12,0.1,4
* Lanczos:
* CoarseCoarse IRL( Nk, Nm, Nstop, poly(lo,hi,order)) 24,36,24,0.002,4.0,61
*/
RealD InverseApproximation(RealD x){
return 1.0/x;
}
template<class Field,class Matrix> class ChebyshevSmoother : public LinearFunction<Field>
{
public:
typedef LinearOperatorBase<Field> FineOperator;
Matrix & _SmootherMatrix;
FineOperator & _SmootherOperator;
Chebyshev<Field> Cheby;
ChebyshevSmoother(RealD _lo,RealD _hi,int _ord, FineOperator &SmootherOperator,Matrix &SmootherMatrix) :
_SmootherOperator(SmootherOperator),
_SmootherMatrix(SmootherMatrix),
Cheby(_lo,_hi,_ord,InverseApproximation)
{};
void operator() (const Field &in, Field &out)
{
Field tmp(in.Grid());
MdagMLinearOperator<Matrix,Field> MdagMOp(_SmootherMatrix);
_SmootherOperator.AdjOp(in,tmp);
Cheby(MdagMOp,tmp,out);
}
};
template<class Field,class Matrix> class MirsSmoother : public LinearFunction<Field>
{
public:
typedef LinearOperatorBase<Field> FineOperator;
Matrix & SmootherMatrix;
FineOperator & SmootherOperator;
RealD tol;
RealD shift;
int maxit;
MirsSmoother(RealD _shift,RealD _tol,int _maxit,FineOperator &_SmootherOperator,Matrix &_SmootherMatrix) :
shift(_shift),tol(_tol),maxit(_maxit),
SmootherOperator(_SmootherOperator),
SmootherMatrix(_SmootherMatrix)
{};
void operator() (const Field &in, Field &out)
{
ZeroGuesser<Field> Guess;
ConjugateGradient<Field> CG(tol,maxit,false);
Field src(in.Grid());
ShiftedMdagMLinearOperator<SparseMatrixBase<Field>,Field> MdagMOp(SmootherMatrix,shift);
SmootherOperator.AdjOp(in,src);
Guess(src,out);
CG(MdagMOp,src,out);
}
};
template<class Fobj,class CComplex,int nbasis, class Matrix, class Guesser, class CoarseSolver>
class MultiGridPreconditioner : public LinearFunction< Lattice<Fobj> > {
public:
typedef Aggregation<Fobj,CComplex,nbasis> Aggregates;
typedef CoarsenedMatrix<Fobj,CComplex,nbasis> CoarseOperator;
typedef typename Aggregation<Fobj,CComplex,nbasis>::CoarseVector CoarseVector;
typedef typename Aggregation<Fobj,CComplex,nbasis>::CoarseMatrix CoarseMatrix;
typedef typename Aggregation<Fobj,CComplex,nbasis>::FineField FineField;
typedef LinearOperatorBase<FineField> FineOperator;
typedef LinearFunction <FineField> FineSmoother;
Aggregates & _Aggregates;
CoarseOperator & _CoarseOperator;
Matrix & _FineMatrix;
FineOperator & _FineOperator;
Guesser & _Guess;
FineSmoother & _Smoother;
CoarseSolver & _CoarseSolve;
int level; void Level(int lv) {level = lv; };
#define GridLogLevel std::cout << GridLogMessage <<std::string(level,'\t')<< " Level "<<level <<" "
MultiGridPreconditioner(Aggregates &Agg, CoarseOperator &Coarse,
FineOperator &Fine,Matrix &FineMatrix,
FineSmoother &Smoother,
Guesser &Guess_,
CoarseSolver &CoarseSolve_)
: _Aggregates(Agg),
_CoarseOperator(Coarse),
_FineOperator(Fine),
_FineMatrix(FineMatrix),
_Smoother(Smoother),
_Guess(Guess_),
_CoarseSolve(CoarseSolve_),
level(1) { }
virtual void operator()(const FineField &in, FineField & out)
{
CoarseVector Csrc(_CoarseOperator.Grid());
CoarseVector Csol(_CoarseOperator.Grid());
FineField vec1(in.Grid());
FineField vec2(in.Grid());
double t;
// Fine Smoother
t=-usecond();
_Smoother(in,out);
t+=usecond();
GridLogLevel << "Smoother took "<< t/1000.0<< "ms" <<std::endl;
// Update the residual
_FineOperator.Op(out,vec1); sub(vec1, in ,vec1);
// Fine to Coarse
t=-usecond();
_Aggregates.ProjectToSubspace (Csrc,vec1);
t+=usecond();
GridLogLevel << "Project to coarse took "<< t/1000.0<< "ms" <<std::endl;
// Coarse correction
t=-usecond();
_CoarseSolve(Csrc,Csol);
t+=usecond();
GridLogLevel << "Coarse solve took "<< t/1000.0<< "ms" <<std::endl;
// Coarse to Fine
t=-usecond();
_Aggregates.PromoteFromSubspace(Csol,vec1);
add(out,out,vec1);
t+=usecond();
GridLogLevel << "Promote to this level took "<< t/1000.0<< "ms" <<std::endl;
// Residual
_FineOperator.Op(out,vec1); sub(vec1 ,in , vec1);
// Fine Smoother
t=-usecond();
_Smoother(vec1,vec2);
t+=usecond();
GridLogLevel << "Smoother took "<< t/1000.0<< "ms" <<std::endl;
add( out,out,vec2);
}
};
int main (int argc, char ** argv)
{
Grid_init(&argc,&argv);
const int Ls=24;
GridCartesian * UGrid = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi());
GridRedBlackCartesian * UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
GridCartesian * FGrid = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid);
GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid);
///////////////////////////////////////////////////
// Construct a coarsened grid; utility for this?
///////////////////////////////////////////////////
std::vector<int> block ({2,2,2,2});
std::vector<int> blockc ({2,2,2,2});
const int nbasis= 40;
const int nbasisc= 40;
auto clatt = GridDefaultLatt();
for(int d=0;d<clatt.size();d++){
clatt[d] = clatt[d]/block[d];
}
auto cclatt = clatt;
for(int d=0;d<clatt.size();d++){
cclatt[d] = clatt[d]/blockc[d];
}
GridCartesian *Coarse4d = SpaceTimeGrid::makeFourDimGrid(clatt, GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi());;
GridCartesian *Coarse5d = SpaceTimeGrid::makeFiveDimGrid(1,Coarse4d);
// GridCartesian *CoarseCoarse4d = SpaceTimeGrid::makeFourDimGrid(cclatt, GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi());;
// GridCartesian *CoarseCoarse5d = SpaceTimeGrid::makeFiveDimGrid(1,CoarseCoarse4d);
std::vector<int> seeds4({1,2,3,4});
std::vector<int> seeds5({5,6,7,8});
std::vector<int> cseeds({5,6,7,8});
GridParallelRNG RNG5(FGrid); RNG5.SeedFixedIntegers(seeds5);
GridParallelRNG RNG4(UGrid); RNG4.SeedFixedIntegers(seeds4);
GridParallelRNG CRNG(Coarse5d);CRNG.SeedFixedIntegers(cseeds);
LatticeFermion src(FGrid); gaussian(RNG5,src);// src=src+g5*src;
LatticeFermion result(FGrid);
LatticeGaugeField Umu(UGrid);
FieldMetaData header;
// std::string file("./ckpoint_lat.4000");
// std::string file("./ckpoint_lat.1000");
// NerscIO::readConfiguration(Umu,header,file);
SU<Nc>::HotConfiguration(RNG4,Umu);
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
std::cout<<GridLogMessage << "Building g5R5 hermitian DWF operator" <<std::endl;
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
RealD mass=0.00078;
RealD M5=1.8;
DomainWallFermionR Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5);
typedef Aggregation<vSpinColourVector,vTComplex,nbasis> Subspace;
typedef CoarsenedMatrix<vSpinColourVector,vTComplex,nbasis> CoarseOperator;
typedef CoarseOperator::CoarseVector CoarseVector;
typedef CoarseOperator::siteVector siteVector;
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
std::cout<<GridLogMessage << "Calling Aggregation class to build subspace" <<std::endl;
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
MdagMLinearOperator<DomainWallFermionR,LatticeFermion> HermDefOp(Ddwf);
Subspace Aggregates(Coarse5d,FGrid,0);
assert ( (nbasis & 0x1)==0);
{
int nb=nbasis/2;
LatticeFermion A(FGrid);
LatticeFermion B(FGrid);
// Aggregates.CreateSubspaceChebyshev(RNG5,HermDefOp,nb,60.0,0.002,1000,800,100,0.0);
// Aggregates.CreateSubspaceChebyshev(RNG5,HermDefOp,nb,60.0,0.02,1000,800,100,0.0);
Aggregates.CreateSubspaceChebyshev(RNG5,HermDefOp,nb,60.0,0.01,400,50,50,0.0); // Slightly faster
for(int n=0;n<nb;n++){
std::cout << GridLogMessage << " G5R5 "<<n<<std::endl;
G5R5(Aggregates.subspace[n+nb],Aggregates.subspace[n]);
std::cout << GridLogMessage << " Projection "<<n<<std::endl;
A = Aggregates.subspace[n];
B = Aggregates.subspace[n+nb];
std::cout << GridLogMessage << " Copy "<<n<<std::endl;
Aggregates.subspace[n] = A+B; // 1+G5 // eigen value of G5R5 is +1
std::cout << GridLogMessage << " P+ "<<n<<std::endl;
Aggregates.subspace[n+nb]= A-B; // 1-G5 // eigen value of G5R5 is -1
std::cout << GridLogMessage << " P- "<<n<<std::endl;
}
}
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
std::cout<<GridLogMessage << "Building coarse representation of Indef operator" <<std::endl;
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
typedef CoarsenedMatrix<vSpinColourVector,vTComplex,nbasis> Level1Op;
typedef CoarsenedMatrix<siteVector,iScalar<vTComplex>,nbasisc> Level2Op;
Gamma5R5HermitianLinearOperator<DomainWallFermionR,LatticeFermion> HermIndefOp(Ddwf);
GridRedBlackCartesian * Coarse4dRB = SpaceTimeGrid::makeFourDimRedBlackGrid(Coarse4d);
std::cout << " Making 5D coarse RB grid " <<std::endl;
GridRedBlackCartesian * Coarse5dRB = SpaceTimeGrid::makeFiveDimRedBlackGrid(1,Coarse4d);
std::cout << " Made 5D coarse RB grid " <<std::endl;
Level1Op LDOp(*Coarse5d,*Coarse5dRB,1);
std::cout << " LDOp.CoarsenOperator " <<std::endl;
LDOp.CoarsenOperator(FGrid,HermIndefOp,Aggregates);
std::cout << " Coarsened Operator " <<std::endl;
//////////////////////////////////////////////////
// Deflate the course space. Recursive multigrid?
//////////////////////////////////////////////////
typedef Aggregation<siteVector,iScalar<vTComplex>,nbasisc> CoarseSubspace;
// CoarseSubspace CoarseAggregates(CoarseCoarse5d,Coarse5d,0);
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
std::cout<<GridLogMessage << "Build deflation space in coarse operator "<< std::endl;
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
MdagMLinearOperator<CoarseOperator,CoarseVector> PosdefLdop(LDOp);
/*
{
int nb=nbasisc/2;
CoarseAggregates.CreateSubspaceChebyshev(CRNG,PosdefLdop,nb,15.0,0.02,1000,800,100,0.0);
for(int n=0;n<nb;n++){
autoView( subspace , CoarseAggregates.subspace[n],CpuWrite);
autoView( subspace_g5, CoarseAggregates.subspace[n+nb],CpuWrite);
for(int nn=0;nn<nb;nn++){
for(int site=0;site<Coarse5d->oSites();site++){
subspace_g5[site](nn) = subspace[site](nn);
subspace_g5[site](nn+nb)=-subspace[site](nn+nb);
}
}
}
}
*/
typedef Level2Op::CoarseVector CoarseCoarseVector;
/*
Level2Op L2Op(*CoarseCoarse5d,1); // Hermitian matrix
HermitianLinearOperator<Level1Op,CoarseVector> L1LinOp(LDOp);
L2Op.CoarsenOperator(Coarse5d,L1LinOp,CoarseAggregates);
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
std::cout<<GridLogMessage << " Running CoarseCoarse grid Lanczos "<< std::endl;
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
MdagMLinearOperator<Level2Op,CoarseCoarseVector> IRLHermOpL2(L2Op);
CoarseCoarseVector cc_src(CoarseCoarse5d); cc_src=1.0;
*/
/*
Chebyshev<CoarseCoarseVector> IRLChebyL2(0.001,15.0,301);
FunctionHermOp<CoarseCoarseVector> IRLOpChebyL2(IRLChebyL2,IRLHermOpL2);
PlainHermOp<CoarseCoarseVector> IRLOpL2 (IRLHermOpL2);
int cNk=24;
int cNm=36;
int cNstop=24;
ImplicitlyRestartedLanczos<CoarseCoarseVector> IRLL2(IRLOpChebyL2,IRLOpL2,cNstop,cNk,cNm,1.0e-3,20);
int cNconv;
std::vector<RealD> eval2(cNm);
std::vector<CoarseCoarseVector> evec2(cNm,CoarseCoarse5d);
IRLL2.calc(eval2,evec2,cc_src,cNconv);
ConjugateGradient<CoarseCoarseVector> CoarseCoarseCG(0.1,1000);
DeflatedGuesser<CoarseCoarseVector> DeflCoarseCoarseGuesser(evec2,eval2);
NormalEquations<CoarseCoarseVector> DeflCoarseCoarseCGNE(L2Op,CoarseCoarseCG,DeflCoarseCoarseGuesser);
*/
/*
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
std::cout<<GridLogMessage << " Running Coarse grid Lanczos "<< std::endl;
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
MdagMLinearOperator<Level1Op,CoarseVector> IRLHermOp(LDOp);
// Chebyshev<CoarseVector> IRLCheby(0.001,15.0,301);
Chebyshev<CoarseVector> IRLCheby(0.03,12.0,101);
FunctionHermOp<CoarseVector> IRLOpCheby(IRLCheby,IRLHermOp);
PlainHermOp<CoarseVector> IRLOp (IRLHermOp);
int Nk=64;
int Nm=128;
int Nstop=Nk;
ImplicitlyRestartedLanczos<CoarseVector> IRL(IRLOpCheby,IRLOp,Nstop,Nk,Nm,1.0e-3,20);
int Nconv;
std::vector<RealD> eval(Nm);
std::vector<CoarseVector> evec(Nm,Coarse5d);
IRL.calc(eval,evec,c_src,Nconv);
*/
CoarseVector c_src(Coarse5d); c_src=1.0;
// DeflatedGuesser<CoarseVector> DeflCoarseGuesser(evec,eval);
// NormalEquations<CoarseVector> DeflCoarseCGNE(LDOp,CoarseCG,DeflCoarseGuesser);
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
std::cout<<GridLogMessage << "Building 3 level Multigrid "<< std::endl;
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
// typedef MultiGridPreconditioner<vSpinColourVector, vTComplex,nbasis, DomainWallFermionR,DeflatedGuesser<CoarseVector> , NormalEquations<CoarseVector> > TwoLevelMG;
typedef MultiGridPreconditioner<vSpinColourVector, vTComplex,nbasis, DomainWallFermionR,ZeroGuesser<CoarseVector> , NormalEquations<CoarseVector> > TwoLevelMG;
typedef MultiGridPreconditioner<siteVector,iScalar<vTComplex>,nbasisc,Level1Op, DeflatedGuesser<CoarseCoarseVector>, NormalEquations<CoarseCoarseVector> > CoarseMG;
typedef MultiGridPreconditioner<vSpinColourVector, vTComplex,nbasis, DomainWallFermionR,ZeroGuesser<CoarseVector>, LinearFunction<CoarseVector> > ThreeLevelMG;
ChebyshevSmoother<LatticeFermion,DomainWallFermionR> FineSmoother(0.25,60.0,12,HermIndefOp,Ddwf);
/*
// MultiGrid preconditioner acting on the coarse space <-> coarsecoarse space
ChebyshevSmoother<CoarseVector, Level1Op > CoarseSmoother(0.1,15.0,3,L1LinOp,LDOp);
// MirsSmoother<CoarseVector, Level1Op > CoarseCGSmoother(0.1,0.1,4,L1LinOp,LDOp);
// MirsSmoother<LatticeFermion,DomainWallFermionR> FineCGSmoother(0.0,0.01,8,HermIndefOp,Ddwf);
CoarseMG Level2Precon (CoarseAggregates, L2Op,
L1LinOp,LDOp,
CoarseSmoother,
DeflCoarseCoarseGuesser,
DeflCoarseCoarseCGNE);
Level2Precon.Level(2);
// PGCR Applying this solver to solve the coarse space problem
PrecGeneralisedConjugateResidual<CoarseVector> l2PGCR(0.1, 100, L1LinOp,Level2Precon,16,16);
l2PGCR.Level(2);
// Wrap the 2nd level solver in a MultiGrid preconditioner acting on the fine space
ZeroGuesser<CoarseVector> CoarseZeroGuesser;
ThreeLevelMG ThreeLevelPrecon(Aggregates, LDOp,
HermIndefOp,Ddwf,
FineSmoother,
CoarseZeroGuesser,
l2PGCR);
ThreeLevelPrecon.Level(1);
// Apply the fine-coarse-coarsecoarse 2 deep MG preconditioner in an outer PGCR on the fine fgrid
PrecGeneralisedConjugateResidual<LatticeFermion> l1PGCR(1.0e-8,1000,HermIndefOp,ThreeLevelPrecon,16,16);
l1PGCR.Level(1);
*/
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
std::cout<<GridLogMessage << "Calling 2 level Multigrid "<< std::endl;
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
result=Zero();
ZeroGuesser<CoarseVector> CoarseZeroGuesser;
ConjugateGradient<CoarseVector> CoarseCG(0.01,1000);
NormalEquations<CoarseVector> CoarseCGNE(LDOp,CoarseCG,CoarseZeroGuesser);
TwoLevelMG TwoLevelPrecon(Aggregates, LDOp,
HermIndefOp,Ddwf,
FineSmoother,
CoarseZeroGuesser,
CoarseCGNE);
TwoLevelPrecon.Level(1);
PrecGeneralisedConjugateResidual<LatticeFermion> l1PGCR(1.0e-8,20,HermIndefOp,TwoLevelPrecon,16,16);
l1PGCR.Level(1);
l1PGCR(src,result);
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
std::cout<<GridLogMessage << "Calling CG "<< std::endl;
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
ConjugateGradient<LatticeFermion> pCG(1.0e-8,60000);
result=Zero();
// pCG(HermDefOp,src,result);
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
std::cout<<GridLogMessage << "Calling red black CG "<< std::endl;
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
result=Zero();
LatticeFermion src_o(FrbGrid);
LatticeFermion result_o(FrbGrid);
pickCheckerboard(Odd,src_o,src);
result_o=Zero();
SchurDiagMooeeOperator<DomainWallFermionR,LatticeFermion> HermOpEO(Ddwf);
// pCG(HermOpEO,src_o,result_o);
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
std::cout<<GridLogMessage << " Fine PowerMethod "<< std::endl;
PowerMethod<LatticeFermion> PM; PM(HermDefOp,src);
std::cout<<GridLogMessage << " Coarse PowerMethod "<< std::endl;
PowerMethod<CoarseVector> cPM; cPM(PosdefLdop,c_src);
// std::cout<<GridLogMessage << " CoarseCoarse PowerMethod "<< std::endl;
// PowerMethod<CoarseCoarseVector> ccPM; ccPM(IRLHermOpL2,cc_src);
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
std::cout<<GridLogMessage << "Done "<< std::endl;
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
Grid_finalize();
}

8
tests/solver/Test_dwf_hdcr_2level.cc
View File

@@ -262,6 +262,8 @@ int main (int argc, char ** argv)
GridCartesian *Coarse4d = SpaceTimeGrid::makeFourDimGrid(clatt, GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi());;
GridCartesian *Coarse5d = SpaceTimeGrid::makeFiveDimGrid(1,Coarse4d);
GridRedBlackCartesian * Coarse4dRB = SpaceTimeGrid::makeFourDimRedBlackGrid(Coarse4d);
GridRedBlackCartesian * Coarse5dRB = SpaceTimeGrid::makeFiveDimRedBlackGrid(1,Coarse4d);
std::vector<int> seeds4({1,2,3,4});
std::vector<int> seeds5({5,6,7,8});
@@ -328,7 +330,7 @@ int main (int argc, char ** argv)
Gamma5R5HermitianLinearOperator<DomainWallFermionR,LatticeFermion> HermIndefOp(Ddwf);
Level1Op LDOp(*Coarse5d,1); LDOp.CoarsenOperator(FGrid,HermIndefOp,Aggregates);
Level1Op LDOp(*Coarse5d,*Coarse5dRB,1); LDOp.CoarsenOperator(FGrid,HermIndefOp,Aggregates);
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
std::cout<<GridLogMessage << " Running Coarse grid Lanczos "<< std::endl;
@@ -352,7 +354,9 @@ int main (int argc, char ** argv)
// ConjugateGradient<CoarseVector> CoarseCG(0.01,1000);
ConjugateGradient<CoarseVector> CoarseCG(0.02,1000);// 14.7s
ConjugateGradient<CoarseVector> CoarseCG(0.01,2000);// 14.7s
eval.resize(0);
evec.resize(0,Coarse5d);
DeflatedGuesser<CoarseVector> DeflCoarseGuesser(evec,eval);
NormalEquations<CoarseVector> DeflCoarseCGNE(LDOp,CoarseCG,DeflCoarseGuesser);

397
tests/solver/Test_dwf_hdcr_48_rb.cc Normal file
View File

@@ -0,0 +1,397 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: ./tests/Test_dwf_hdcr.cc
Copyright (C) 2015
Author: Antonin Portelli <antonin.portelli@me.com>
Author: Peter Boyle <paboyle@ph.ed.ac.uk>
Author: paboyle <paboyle@ph.ed.ac.uk>
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License along
with this program; if not, write to the Free Software Foundation, Inc.,
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
See the full license in the file "LICENSE" in the top level distribution directory
*************************************************************************************/
/* END LEGAL */
#include <Grid/Grid.h>
#include <Grid/algorithms/iterative/PrecGeneralisedConjugateResidual.h>
using namespace std;
using namespace Grid;
/* Params
* Grid:
* block1(4)
* block2(4)
*
* Subspace
* * Fine : Subspace(nbasis,hi,lo,order,first,step) -- 32, 60,0.02,500,100,100
* * Coarse: Subspace(nbasis,hi,lo,order,first,step) -- 32, 18,0.02,500,100,100
* Smoother:
* * Fine: Cheby(hi, lo, order) -- 60,0.5,10
* * Coarse: Cheby(hi, lo, order) -- 12,0.1,4
* Lanczos:
* CoarseCoarse IRL( Nk, Nm, Nstop, poly(lo,hi,order)) 24,36,24,0.002,4.0,61
*/
RealD InverseApproximation(RealD x){
return 1.0/x;
}
template<class Field> class SolverWrapper : public LinearFunction<Field> {
private:
CheckerBoardedSparseMatrixBase<Field> & _Matrix;
SchurRedBlackBase<Field> & _Solver;
public:
/////////////////////////////////////////////////////
// Wrap the usual normal equations trick
/////////////////////////////////////////////////////
SolverWrapper(CheckerBoardedSparseMatrixBase<Field> &Matrix,
SchurRedBlackBase<Field> &Solver)
: _Matrix(Matrix), _Solver(Solver) {};
void operator() (const Field &in, Field &out){
_Solver(_Matrix,in,out); // Mdag M out = Mdag in
}
};
template<class Field,class Matrix> class ChebyshevSmoother : public LinearFunction<Field>
{
public:
typedef LinearOperatorBase<Field> FineOperator;
Matrix & _SmootherMatrix;
FineOperator & _SmootherOperator;
Chebyshev<Field> Cheby;
ChebyshevSmoother(RealD _lo,RealD _hi,int _ord, FineOperator &SmootherOperator,Matrix &SmootherMatrix) :
_SmootherOperator(SmootherOperator),
_SmootherMatrix(SmootherMatrix),
Cheby(_lo,_hi,_ord,InverseApproximation)
{};
void operator() (const Field &in, Field &out)
{
Field tmp(in.Grid());
MdagMLinearOperator<Matrix,Field> MdagMOp(_SmootherMatrix);
_SmootherOperator.AdjOp(in,tmp);
Cheby(MdagMOp,tmp,out);
}
};
template<class Field,class Matrix> class MirsSmoother : public LinearFunction<Field>
{
public:
typedef LinearOperatorBase<Field> FineOperator;
Matrix & SmootherMatrix;
FineOperator & SmootherOperator;
RealD tol;
RealD shift;
int maxit;
MirsSmoother(RealD _shift,RealD _tol,int _maxit,FineOperator &_SmootherOperator,Matrix &_SmootherMatrix) :
shift(_shift),tol(_tol),maxit(_maxit),
SmootherOperator(_SmootherOperator),
SmootherMatrix(_SmootherMatrix)
{};
void operator() (const Field &in, Field &out)
{
ZeroGuesser<Field> Guess;
ConjugateGradient<Field> CG(tol,maxit,false);
Field src(in.Grid());
ShiftedMdagMLinearOperator<SparseMatrixBase<Field>,Field> MdagMOp(SmootherMatrix,shift);
SmootherOperator.AdjOp(in,src);
Guess(src,out);
CG(MdagMOp,src,out);
}
};
template<class Fobj,class CComplex,int nbasis, class Matrix, class Guesser, class CoarseSolver>
class MultiGridPreconditioner : public LinearFunction< Lattice<Fobj> > {
public:
typedef Aggregation<Fobj,CComplex,nbasis> Aggregates;
typedef CoarsenedMatrix<Fobj,CComplex,nbasis> CoarseOperator;
typedef typename Aggregation<Fobj,CComplex,nbasis>::CoarseVector CoarseVector;
typedef typename Aggregation<Fobj,CComplex,nbasis>::CoarseMatrix CoarseMatrix;
typedef typename Aggregation<Fobj,CComplex,nbasis>::FineField FineField;
typedef LinearOperatorBase<FineField> FineOperator;
typedef LinearFunction <FineField> FineSmoother;
Aggregates & _Aggregates;
CoarseOperator & _CoarseOperator;
Matrix & _FineMatrix;
FineOperator & _FineOperator;
Guesser & _Guess;
FineSmoother & _Smoother;
CoarseSolver & _CoarseSolve;
int level; void Level(int lv) {level = lv; };
#define GridLogLevel std::cout << GridLogMessage <<std::string(level,'\t')<< " Level "<<level <<" "
MultiGridPreconditioner(Aggregates &Agg, CoarseOperator &Coarse,
FineOperator &Fine,Matrix &FineMatrix,
FineSmoother &Smoother,
Guesser &Guess_,
CoarseSolver &CoarseSolve_)
: _Aggregates(Agg),
_CoarseOperator(Coarse),
_FineOperator(Fine),
_FineMatrix(FineMatrix),
_Smoother(Smoother),
_Guess(Guess_),
_CoarseSolve(CoarseSolve_),
level(1) { }
virtual void operator()(const FineField &in, FineField & out)
{
CoarseVector Csrc(_CoarseOperator.Grid());
CoarseVector Csol(_CoarseOperator.Grid());
FineField vec1(in.Grid());
FineField vec2(in.Grid());
double t;
// Fine Smoother
t=-usecond();
_Smoother(in,out);
t+=usecond();
GridLogLevel << "Smoother took "<< t/1000.0<< "ms" <<std::endl;
// Update the residual
_FineOperator.Op(out,vec1); sub(vec1, in ,vec1);
// Fine to Coarse
t=-usecond();
_Aggregates.ProjectToSubspace (Csrc,vec1);
t+=usecond();
GridLogLevel << "Project to coarse took "<< t/1000.0<< "ms" <<std::endl;
// Coarse correction
t=-usecond();
_CoarseSolve(Csrc,Csol);
t+=usecond();
GridLogLevel << "Coarse solve took "<< t/1000.0<< "ms" <<std::endl;
// Coarse to Fine
t=-usecond();
_Aggregates.PromoteFromSubspace(Csol,vec1);
add(out,out,vec1);
t+=usecond();
GridLogLevel << "Promote to this level took "<< t/1000.0<< "ms" <<std::endl;
// Residual
_FineOperator.Op(out,vec1); sub(vec1 ,in , vec1);
// Fine Smoother
t=-usecond();
_Smoother(vec1,vec2);
t+=usecond();
GridLogLevel << "Smoother took "<< t/1000.0<< "ms" <<std::endl;
add( out,out,vec2);
}
};
int main (int argc, char ** argv)
{
Grid_init(&argc,&argv);
const int Ls=24;
GridCartesian * UGrid = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi());
GridRedBlackCartesian * UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
GridCartesian * FGrid = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid);
GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid);
///////////////////////////////////////////////////
// Construct a coarsened grid; utility for this?
///////////////////////////////////////////////////
std::vector<int> block ({2,2,2,2});
//std::vector<int> block ({2,2,2,2});
const int nbasis= 40;
const int nbasisc= 40;
auto clatt = GridDefaultLatt();
for(int d=0;d<clatt.size();d++){
clatt[d] = clatt[d]/block[d];
}
GridCartesian *Coarse4d = SpaceTimeGrid::makeFourDimGrid(clatt, GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi());;
GridCartesian *Coarse5d = SpaceTimeGrid::makeFiveDimGrid(1,Coarse4d);
// GridCartesian *CoarseCoarse4d = SpaceTimeGrid::makeFourDimGrid(cclatt, GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi());;
// GridCartesian *CoarseCoarse5d = SpaceTimeGrid::makeFiveDimGrid(1,CoarseCoarse4d);
std::vector<int> seeds4({1,2,3,4});
std::vector<int> seeds5({5,6,7,8});
std::vector<int> cseeds({5,6,7,8});
GridParallelRNG RNG5(FGrid); RNG5.SeedFixedIntegers(seeds5);
GridParallelRNG RNG4(UGrid); RNG4.SeedFixedIntegers(seeds4);
GridParallelRNG CRNG(Coarse5d);CRNG.SeedFixedIntegers(cseeds);
LatticeFermion src(FGrid); gaussian(RNG5,src);// src=src+g5*src;
LatticeFermion result(FGrid);
LatticeGaugeField Umu(UGrid);
FieldMetaData header;
//std::string file("./ckpoint_lat.4000");
std::string file("./ckpoint_lat.1000");
NerscIO::readConfiguration(Umu,header,file);
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
std::cout<<GridLogMessage << "Building g5R5 hermitian DWF operator" <<std::endl;
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
RealD mass=0.00078;
RealD M5=1.8;
DomainWallFermionR Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5);
typedef Aggregation<vSpinColourVector,vTComplex,nbasis> Subspace;
typedef CoarsenedMatrix<vSpinColourVector,vTComplex,nbasis> CoarseOperator;
typedef CoarseOperator::CoarseVector CoarseVector;
typedef CoarseOperator::siteVector siteVector;
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
std::cout<<GridLogMessage << "Calling Aggregation class to build subspace" <<std::endl;
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
MdagMLinearOperator<DomainWallFermionR,LatticeFermion> HermDefOp(Ddwf);
Subspace Aggregates(Coarse5d,FGrid,0);
assert ( (nbasis & 0x1)==0);
{
int nb=nbasis/2;
LatticeFermion A(FGrid);
LatticeFermion B(FGrid);
// Aggregates.CreateSubspaceChebyshev(RNG5,HermDefOp,nb,60.0,0.002,1000,800,100,0.0);
// Aggregates.CreateSubspaceChebyshev(RNG5,HermDefOp,nb,60.0,0.02,1000,800,100,0.0);
Aggregates.CreateSubspaceChebyshev(RNG5,HermDefOp,nb,60.0,0.01,1000,100,100,0.0); // Slightly faster
for(int n=0;n<nb;n++){
std::cout << GridLogMessage << " G5R5 "<<n<<std::endl;
G5R5(Aggregates.subspace[n+nb],Aggregates.subspace[n]);
std::cout << GridLogMessage << " Projection "<<n<<std::endl;
A = Aggregates.subspace[n];
B = Aggregates.subspace[n+nb];
std::cout << GridLogMessage << " Copy "<<n<<std::endl;
Aggregates.subspace[n] = A+B; // 1+G5 // eigen value of G5R5 is +1
std::cout << GridLogMessage << " P+ "<<n<<std::endl;
Aggregates.subspace[n+nb]= A-B; // 1-G5 // eigen value of G5R5 is -1
std::cout << GridLogMessage << " P- "<<n<<std::endl;
}
}
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
std::cout<<GridLogMessage << "Building coarse representation of Indef operator" <<std::endl;
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
typedef CoarsenedMatrix<vSpinColourVector,vTComplex,nbasis> Level1Op;
typedef CoarsenedMatrix<siteVector,iScalar<vTComplex>,nbasisc> Level2Op;
Gamma5R5HermitianLinearOperator<DomainWallFermionR,LatticeFermion> HermIndefOp(Ddwf);
GridRedBlackCartesian * Coarse4dRB = SpaceTimeGrid::makeFourDimRedBlackGrid(Coarse4d);
GridRedBlackCartesian * Coarse5dRB = SpaceTimeGrid::makeFiveDimRedBlackGrid(1,Coarse4d);
Level1Op LDOp(*Coarse5d,*Coarse5dRB,1); LDOp.CoarsenOperator(FGrid,HermIndefOp,Aggregates);
//////////////////////////////////////////////////
// Deflate the course space. Recursive multigrid?
//////////////////////////////////////////////////
typedef Aggregation<siteVector,iScalar<vTComplex>,nbasisc> CoarseSubspace;
// CoarseSubspace CoarseAggregates(CoarseCoarse5d,Coarse5d,0);
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
std::cout<<GridLogMessage << "Build deflation space in coarse operator "<< std::endl;
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
MdagMLinearOperator<CoarseOperator,CoarseVector> PosdefLdop(LDOp);
typedef Level2Op::CoarseVector CoarseCoarseVector;
CoarseVector c_src(Coarse5d); c_src=1.0;
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
std::cout<<GridLogMessage << "Building 3 level Multigrid "<< std::endl;
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
typedef MultiGridPreconditioner<vSpinColourVector, vTComplex,nbasis, DomainWallFermionR,ZeroGuesser<CoarseVector> , SolverWrapper<CoarseVector> > TwoLevelMG;
typedef MultiGridPreconditioner<siteVector,iScalar<vTComplex>,nbasisc,Level1Op, DeflatedGuesser<CoarseCoarseVector>, NormalEquations<CoarseCoarseVector> > CoarseMG;
typedef MultiGridPreconditioner<vSpinColourVector, vTComplex,nbasis, DomainWallFermionR,ZeroGuesser<CoarseVector>, LinearFunction<CoarseVector> > ThreeLevelMG;
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
std::cout<<GridLogMessage << "Calling 2 level Multigrid "<< std::endl;
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
std::vector<RealD> tols({0.015});
std::vector<int> ords({12});
std::vector<RealD> los({0.8});
for(int l=0;l<los.size();l++){
for(int o=0;o<ords.size();o++){
for(int t=0;t<tols.size();t++){
result=Zero();
std::cout << GridLogMessage <<" tol " << tols[t] << " cheby order " <<ords[o]<< " lo "<<los[l] <<std::endl;
ChebyshevSmoother<LatticeFermion,DomainWallFermionR> FineSmoother(los[l],60.0,ords[o],HermIndefOp,Ddwf);
ZeroGuesser<CoarseVector> CoarseZeroGuesser;
ConjugateGradient<CoarseVector> CoarseCG(tols[t],10000);
SchurRedBlackDiagMooeeSolve<CoarseVector> CoarseRBCG(CoarseCG);
SolverWrapper<CoarseVector> CoarseSolver(LDOp,CoarseRBCG);
TwoLevelMG TwoLevelPrecon(Aggregates, LDOp,
HermIndefOp,Ddwf,
FineSmoother,
CoarseZeroGuesser,
CoarseSolver);
TwoLevelPrecon.Level(1);
PrecGeneralisedConjugateResidual<LatticeFermion> l1PGCR(1.0e-8,20,HermIndefOp,TwoLevelPrecon,16,16);
l1PGCR.Level(1);
l1PGCR(src,result);
}}}
ConjugateGradient<LatticeFermion> pCG(1.0e-8,60000);
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
std::cout<<GridLogMessage << "Calling red black CG "<< std::endl;
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
result=Zero();
LatticeFermion src_o(FrbGrid);
LatticeFermion result_o(FrbGrid);
pickCheckerboard(Odd,src_o,src);
result_o=Zero();
SchurDiagMooeeOperator<DomainWallFermionR,LatticeFermion> HermOpEO(Ddwf);
pCG(HermOpEO,src_o,result_o);
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
std::cout<<GridLogMessage << "Calling CG "<< std::endl;
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
result=Zero();
pCG(HermDefOp,src,result);
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
std::cout<<GridLogMessage << " Fine PowerMethod "<< std::endl;
PowerMethod<LatticeFermion> PM; PM(HermDefOp,src);
std::cout<<GridLogMessage << " Coarse PowerMethod "<< std::endl;
PowerMethod<CoarseVector> cPM; cPM(PosdefLdop,c_src);
// std::cout<<GridLogMessage << " CoarseCoarse PowerMethod "<< std::endl;
// PowerMethod<CoarseCoarseVector> ccPM; ccPM(IRLHermOpL2,cc_src);
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
std::cout<<GridLogMessage << "Done "<< std::endl;
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
Grid_finalize();
}

473
tests/solver/Test_dwf_hdcr_48_regression.cc Normal file
View File

@@ -0,0 +1,473 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: ./tests/Test_dwf_hdcr.cc
Copyright (C) 2015
Author: Antonin Portelli <antonin.portelli@me.com>
Author: Peter Boyle <paboyle@ph.ed.ac.uk>
Author: paboyle <paboyle@ph.ed.ac.uk>
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License along
with this program; if not, write to the Free Software Foundation, Inc.,
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
See the full license in the file "LICENSE" in the top level distribution directory
*************************************************************************************/
/* END LEGAL */
#include <Grid/Grid.h>
#include <Grid/algorithms/iterative/PrecGeneralisedConjugateResidual.h>
using namespace std;
using namespace Grid;
/* Params
* Grid:
* block1(4)
* block2(4)
*
* Subspace
* * Fine : Subspace(nbasis,hi,lo,order,first,step) -- 32, 60,0.02,500,100,100
* * Coarse: Subspace(nbasis,hi,lo,order,first,step) -- 32, 18,0.02,500,100,100
* Smoother:
* * Fine: Cheby(hi, lo, order) -- 60,0.5,10
* * Coarse: Cheby(hi, lo, order) -- 12,0.1,4
* Lanczos:
* CoarseCoarse IRL( Nk, Nm, Nstop, poly(lo,hi,order)) 24,36,24,0.002,4.0,61
*/
RealD InverseApproximation(RealD x){
return 1.0/x;
}
template<class Field,class Matrix> class ChebyshevSmoother : public LinearFunction<Field>
{
public:
typedef LinearOperatorBase<Field> FineOperator;
Matrix & _SmootherMatrix;
FineOperator & _SmootherOperator;
Chebyshev<Field> Cheby;
ChebyshevSmoother(RealD _lo,RealD _hi,int _ord, FineOperator &SmootherOperator,Matrix &SmootherMatrix) :
_SmootherOperator(SmootherOperator),
_SmootherMatrix(SmootherMatrix),
Cheby(_lo,_hi,_ord,InverseApproximation)
{};
void operator() (const Field &in, Field &out)
{
Field tmp(in.Grid());
MdagMLinearOperator<Matrix,Field> MdagMOp(_SmootherMatrix);
_SmootherOperator.AdjOp(in,tmp);
Cheby(MdagMOp,tmp,out);
}
};
template<class Field,class Matrix> class MirsSmoother : public LinearFunction<Field>
{
public:
typedef LinearOperatorBase<Field> FineOperator;
Matrix & SmootherMatrix;
FineOperator & SmootherOperator;
RealD tol;
RealD shift;
int maxit;
MirsSmoother(RealD _shift,RealD _tol,int _maxit,FineOperator &_SmootherOperator,Matrix &_SmootherMatrix) :
shift(_shift),tol(_tol),maxit(_maxit),
SmootherOperator(_SmootherOperator),
SmootherMatrix(_SmootherMatrix)
{};
void operator() (const Field &in, Field &out)
{
ZeroGuesser<Field> Guess;
ConjugateGradient<Field> CG(tol,maxit,false);
Field src(in.Grid());
ShiftedMdagMLinearOperator<SparseMatrixBase<Field>,Field> MdagMOp(SmootherMatrix,shift);
SmootherOperator.AdjOp(in,src);
Guess(src,out);
CG(MdagMOp,src,out);
}
};
template<class Fobj,class CComplex,int nbasis, class Matrix, class Guesser, class CoarseSolver>
class MultiGridPreconditioner : public LinearFunction< Lattice<Fobj> > {
public:
typedef Aggregation<Fobj,CComplex,nbasis> Aggregates;
typedef CoarsenedMatrix<Fobj,CComplex,nbasis> CoarseOperator;
typedef typename Aggregation<Fobj,CComplex,nbasis>::CoarseVector CoarseVector;
typedef typename Aggregation<Fobj,CComplex,nbasis>::CoarseMatrix CoarseMatrix;
typedef typename Aggregation<Fobj,CComplex,nbasis>::FineField FineField;
typedef LinearOperatorBase<FineField> FineOperator;
typedef LinearFunction <FineField> FineSmoother;
Aggregates & _Aggregates;
CoarseOperator & _CoarseOperator;
Matrix & _FineMatrix;
FineOperator & _FineOperator;
Guesser & _Guess;
FineSmoother & _Smoother;
CoarseSolver & _CoarseSolve;
int level; void Level(int lv) {level = lv; };
#define GridLogLevel std::cout << GridLogMessage <<std::string(level,'\t')<< " Level "<<level <<" "
MultiGridPreconditioner(Aggregates &Agg, CoarseOperator &Coarse,
FineOperator &Fine,Matrix &FineMatrix,
FineSmoother &Smoother,
Guesser &Guess_,
CoarseSolver &CoarseSolve_)
: _Aggregates(Agg),
_CoarseOperator(Coarse),
_FineOperator(Fine),
_FineMatrix(FineMatrix),
_Smoother(Smoother),
_Guess(Guess_),
_CoarseSolve(CoarseSolve_),
level(1) { }
virtual void operator()(const FineField &in, FineField & out)
{
CoarseVector Csrc(_CoarseOperator.Grid());
CoarseVector Csol(_CoarseOperator.Grid());
FineField vec1(in.Grid());
FineField vec2(in.Grid());
double t;
// Fine Smoother
t=-usecond();
_Smoother(in,out);
t+=usecond();
GridLogLevel << "Smoother took "<< t/1000.0<< "ms" <<std::endl;
// Update the residual
_FineOperator.Op(out,vec1); sub(vec1, in ,vec1);
// Fine to Coarse
t=-usecond();
_Aggregates.ProjectToSubspace (Csrc,vec1);
t+=usecond();
GridLogLevel << "Project to coarse took "<< t/1000.0<< "ms" <<std::endl;
// Coarse correction
t=-usecond();
_CoarseSolve(Csrc,Csol);
t+=usecond();
GridLogLevel << "Coarse solve took "<< t/1000.0<< "ms" <<std::endl;
// Coarse to Fine
t=-usecond();
_Aggregates.PromoteFromSubspace(Csol,vec1);
add(out,out,vec1);
t+=usecond();
GridLogLevel << "Promote to this level took "<< t/1000.0<< "ms" <<std::endl;
// Residual
_FineOperator.Op(out,vec1); sub(vec1 ,in , vec1);
// Fine Smoother
t=-usecond();
_Smoother(vec1,vec2);
t+=usecond();
GridLogLevel << "Smoother took "<< t/1000.0<< "ms" <<std::endl;
add( out,out,vec2);
}
};
int main (int argc, char ** argv)
{
Grid_init(&argc,&argv);
const int Ls=24;
GridCartesian * UGrid = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi());
GridRedBlackCartesian * UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
GridCartesian * FGrid = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid);
GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid);
///////////////////////////////////////////////////
// Construct a coarsened grid; utility for this?
///////////////////////////////////////////////////
std::vector<int> block ({2,2,2,2});
std::vector<int> blockc ({2,2,2,2});
const int nbasis= 40;
const int nbasisc= 40;
auto clatt = GridDefaultLatt();
for(int d=0;d<clatt.size();d++){
clatt[d] = clatt[d]/block[d];
}
auto cclatt = clatt;
for(int d=0;d<clatt.size();d++){
cclatt[d] = clatt[d]/blockc[d];
}
GridCartesian *Coarse4d = SpaceTimeGrid::makeFourDimGrid(clatt, GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi());;
GridCartesian *Coarse5d = SpaceTimeGrid::makeFiveDimGrid(1,Coarse4d);
// GridCartesian *CoarseCoarse4d = SpaceTimeGrid::makeFourDimGrid(cclatt, GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi());;
// GridCartesian *CoarseCoarse5d = SpaceTimeGrid::makeFiveDimGrid(1,CoarseCoarse4d);
std::vector<int> seeds4({1,2,3,4});
std::vector<int> seeds5({5,6,7,8});
std::vector<int> cseeds({5,6,7,8});
GridParallelRNG RNG5(FGrid); RNG5.SeedFixedIntegers(seeds5);
GridParallelRNG RNG4(UGrid); RNG4.SeedFixedIntegers(seeds4);
GridParallelRNG CRNG(Coarse5d);CRNG.SeedFixedIntegers(cseeds);
LatticeFermion src(FGrid); gaussian(RNG5,src);// src=src+g5*src;
LatticeFermion result(FGrid);
LatticeGaugeField Umu(UGrid);
FieldMetaData header;
// std::string file("./ckpoint_lat.4000");
std::string file("./ckpoint_lat.1000");
NerscIO::readConfiguration(Umu,header,file);
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
std::cout<<GridLogMessage << "Building g5R5 hermitian DWF operator" <<std::endl;
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
RealD mass=0.00078;
RealD M5=1.8;
DomainWallFermionR Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5);
typedef Aggregation<vSpinColourVector,vTComplex,nbasis> Subspace;
typedef CoarsenedMatrix<vSpinColourVector,vTComplex,nbasis> CoarseOperator;
typedef CoarseOperator::CoarseVector CoarseVector;
typedef CoarseOperator::siteVector siteVector;
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
std::cout<<GridLogMessage << "Calling Aggregation class to build subspace" <<std::endl;
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
MdagMLinearOperator<DomainWallFermionR,LatticeFermion> HermDefOp(Ddwf);
Subspace Aggregates(Coarse5d,FGrid,0);
assert ( (nbasis & 0x1)==0);
{
int nb=nbasis/2;
LatticeFermion A(FGrid);
LatticeFermion B(FGrid);
// Aggregates.CreateSubspaceChebyshev(RNG5,HermDefOp,nb,60.0,0.002,1000,800,100,0.0);
// Aggregates.CreateSubspaceChebyshev(RNG5,HermDefOp,nb,60.0,0.02,1000,800,100,0.0);
Aggregates.CreateSubspaceChebyshev(RNG5,HermDefOp,nb,60.0,0.01,1000,100,100,0.0); // Slightly faster
for(int n=0;n<nb;n++){
std::cout << GridLogMessage << " G5R5 "<<n<<std::endl;
G5R5(Aggregates.subspace[n+nb],Aggregates.subspace[n]);
std::cout << GridLogMessage << " Projection "<<n<<std::endl;
A = Aggregates.subspace[n];
B = Aggregates.subspace[n+nb];
std::cout << GridLogMessage << " Copy "<<n<<std::endl;
Aggregates.subspace[n] = A+B; // 1+G5 // eigen value of G5R5 is +1
std::cout << GridLogMessage << " P+ "<<n<<std::endl;
Aggregates.subspace[n+nb]= A-B; // 1-G5 // eigen value of G5R5 is -1
std::cout << GridLogMessage << " P- "<<n<<std::endl;
}
}
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
std::cout<<GridLogMessage << "Building coarse representation of Indef operator" <<std::endl;
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
typedef CoarsenedMatrix<vSpinColourVector,vTComplex,nbasis> Level1Op;
typedef CoarsenedMatrix<siteVector,iScalar<vTComplex>,nbasisc> Level2Op;
Gamma5R5HermitianLinearOperator<DomainWallFermionR,LatticeFermion> HermIndefOp(Ddwf);
GridRedBlackCartesian * Coarse4dRB = SpaceTimeGrid::makeFourDimRedBlackGrid(Coarse4d);
std::cout << " Making 5D coarse RB grid " <<std::endl;
GridRedBlackCartesian * Coarse5dRB = SpaceTimeGrid::makeFiveDimRedBlackGrid(1,Coarse4d);
std::cout << " Made 5D coarse RB grid " <<std::endl;
Level1Op LDOp(*Coarse5d,*Coarse5dRB,1); LDOp.CoarsenOperator(FGrid,HermIndefOp,Aggregates);
//////////////////////////////////////////////////
// Deflate the course space. Recursive multigrid?
//////////////////////////////////////////////////
typedef Aggregation<siteVector,iScalar<vTComplex>,nbasisc> CoarseSubspace;
// CoarseSubspace CoarseAggregates(CoarseCoarse5d,Coarse5d,0);
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
std::cout<<GridLogMessage << "Build deflation space in coarse operator "<< std::endl;
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
MdagMLinearOperator<CoarseOperator,CoarseVector> PosdefLdop(LDOp);
/*
{
int nb=nbasisc/2;
CoarseAggregates.CreateSubspaceChebyshev(CRNG,PosdefLdop,nb,15.0,0.02,1000,800,100,0.0);
for(int n=0;n<nb;n++){
autoView( subspace , CoarseAggregates.subspace[n],CpuWrite);
autoView( subspace_g5, CoarseAggregates.subspace[n+nb],CpuWrite);
for(int nn=0;nn<nb;nn++){
for(int site=0;site<Coarse5d->oSites();site++){
subspace_g5[site](nn) = subspace[site](nn);
subspace_g5[site](nn+nb)=-subspace[site](nn+nb);
}
}
}
}
*/
typedef Level2Op::CoarseVector CoarseCoarseVector;
/*
Level2Op L2Op(*CoarseCoarse5d,1); // Hermitian matrix
HermitianLinearOperator<Level1Op,CoarseVector> L1LinOp(LDOp);
L2Op.CoarsenOperator(Coarse5d,L1LinOp,CoarseAggregates);
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
std::cout<<GridLogMessage << " Running CoarseCoarse grid Lanczos "<< std::endl;
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
MdagMLinearOperator<Level2Op,CoarseCoarseVector> IRLHermOpL2(L2Op);
CoarseCoarseVector cc_src(CoarseCoarse5d); cc_src=1.0;
*/
/*
Chebyshev<CoarseCoarseVector> IRLChebyL2(0.001,15.0,301);
FunctionHermOp<CoarseCoarseVector> IRLOpChebyL2(IRLChebyL2,IRLHermOpL2);
PlainHermOp<CoarseCoarseVector> IRLOpL2 (IRLHermOpL2);
int cNk=24;
int cNm=36;
int cNstop=24;
ImplicitlyRestartedLanczos<CoarseCoarseVector> IRLL2(IRLOpChebyL2,IRLOpL2,cNstop,cNk,cNm,1.0e-3,20);
int cNconv;
std::vector<RealD> eval2(cNm);
std::vector<CoarseCoarseVector> evec2(cNm,CoarseCoarse5d);
IRLL2.calc(eval2,evec2,cc_src,cNconv);
ConjugateGradient<CoarseCoarseVector> CoarseCoarseCG(0.1,1000);
DeflatedGuesser<CoarseCoarseVector> DeflCoarseCoarseGuesser(evec2,eval2);
NormalEquations<CoarseCoarseVector> DeflCoarseCoarseCGNE(L2Op,CoarseCoarseCG,DeflCoarseCoarseGuesser);
*/
/*
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
std::cout<<GridLogMessage << " Running Coarse grid Lanczos "<< std::endl;
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
MdagMLinearOperator<Level1Op,CoarseVector> IRLHermOp(LDOp);
// Chebyshev<CoarseVector> IRLCheby(0.001,15.0,301);
Chebyshev<CoarseVector> IRLCheby(0.03,12.0,101);
FunctionHermOp<CoarseVector> IRLOpCheby(IRLCheby,IRLHermOp);
PlainHermOp<CoarseVector> IRLOp (IRLHermOp);
int Nk=64;
int Nm=128;
int Nstop=Nk;
ImplicitlyRestartedLanczos<CoarseVector> IRL(IRLOpCheby,IRLOp,Nstop,Nk,Nm,1.0e-3,20);
int Nconv;
std::vector<RealD> eval(Nm);
std::vector<CoarseVector> evec(Nm,Coarse5d);
IRL.calc(eval,evec,c_src,Nconv);
*/
CoarseVector c_src(Coarse5d); c_src=1.0;
// DeflatedGuesser<CoarseVector> DeflCoarseGuesser(evec,eval);
// NormalEquations<CoarseVector> DeflCoarseCGNE(LDOp,CoarseCG,DeflCoarseGuesser);
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
std::cout<<GridLogMessage << "Building 3 level Multigrid "<< std::endl;
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
// typedef MultiGridPreconditioner<vSpinColourVector, vTComplex,nbasis, DomainWallFermionR,DeflatedGuesser<CoarseVector> , NormalEquations<CoarseVector> > TwoLevelMG;
typedef MultiGridPreconditioner<vSpinColourVector, vTComplex,nbasis, DomainWallFermionR,ZeroGuesser<CoarseVector> , NormalEquations<CoarseVector> > TwoLevelMG;
typedef MultiGridPreconditioner<siteVector,iScalar<vTComplex>,nbasisc,Level1Op, DeflatedGuesser<CoarseCoarseVector>, NormalEquations<CoarseCoarseVector> > CoarseMG;
typedef MultiGridPreconditioner<vSpinColourVector, vTComplex,nbasis, DomainWallFermionR,ZeroGuesser<CoarseVector>, LinearFunction<CoarseVector> > ThreeLevelMG;
ChebyshevSmoother<LatticeFermion,DomainWallFermionR> FineSmoother(0.25,60.0,12,HermIndefOp,Ddwf);
/*
// MultiGrid preconditioner acting on the coarse space <-> coarsecoarse space
ChebyshevSmoother<CoarseVector, Level1Op > CoarseSmoother(0.1,15.0,3,L1LinOp,LDOp);
// MirsSmoother<CoarseVector, Level1Op > CoarseCGSmoother(0.1,0.1,4,L1LinOp,LDOp);
// MirsSmoother<LatticeFermion,DomainWallFermionR> FineCGSmoother(0.0,0.01,8,HermIndefOp,Ddwf);
CoarseMG Level2Precon (CoarseAggregates, L2Op,
L1LinOp,LDOp,
CoarseSmoother,
DeflCoarseCoarseGuesser,
DeflCoarseCoarseCGNE);
Level2Precon.Level(2);
// PGCR Applying this solver to solve the coarse space problem
PrecGeneralisedConjugateResidual<CoarseVector> l2PGCR(0.1, 100, L1LinOp,Level2Precon,16,16);
l2PGCR.Level(2);
// Wrap the 2nd level solver in a MultiGrid preconditioner acting on the fine space
ZeroGuesser<CoarseVector> CoarseZeroGuesser;
ThreeLevelMG ThreeLevelPrecon(Aggregates, LDOp,
HermIndefOp,Ddwf,
FineSmoother,
CoarseZeroGuesser,
l2PGCR);
ThreeLevelPrecon.Level(1);
// Apply the fine-coarse-coarsecoarse 2 deep MG preconditioner in an outer PGCR on the fine fgrid
PrecGeneralisedConjugateResidual<LatticeFermion> l1PGCR(1.0e-8,1000,HermIndefOp,ThreeLevelPrecon,16,16);
l1PGCR.Level(1);
*/
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
std::cout<<GridLogMessage << "Calling 2 level Multigrid "<< std::endl;
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
result=Zero();
ZeroGuesser<CoarseVector> CoarseZeroGuesser;
ConjugateGradient<CoarseVector> CoarseCG(0.01,1000);
NormalEquations<CoarseVector> CoarseCGNE(LDOp,CoarseCG,CoarseZeroGuesser);
TwoLevelMG TwoLevelPrecon(Aggregates, LDOp,
HermIndefOp,Ddwf,
FineSmoother,
CoarseZeroGuesser,
CoarseCGNE);
TwoLevelPrecon.Level(1);
PrecGeneralisedConjugateResidual<LatticeFermion> l1PGCR(1.0e-8,20,HermIndefOp,TwoLevelPrecon,16,16);
l1PGCR.Level(1);
l1PGCR(src,result);
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
std::cout<<GridLogMessage << "Calling CG "<< std::endl;
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
ConjugateGradient<LatticeFermion> pCG(1.0e-8,60000);
result=Zero();
// pCG(HermDefOp,src,result);
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
std::cout<<GridLogMessage << "Calling red black CG "<< std::endl;
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
result=Zero();
LatticeFermion src_o(FrbGrid);
LatticeFermion result_o(FrbGrid);
pickCheckerboard(Odd,src_o,src);
result_o=Zero();
SchurDiagMooeeOperator<DomainWallFermionR,LatticeFermion> HermOpEO(Ddwf);
pCG(HermOpEO,src_o,result_o);
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
std::cout<<GridLogMessage << " Fine PowerMethod "<< std::endl;
PowerMethod<LatticeFermion> PM; PM(HermDefOp,src);
std::cout<<GridLogMessage << " Coarse PowerMethod "<< std::endl;
PowerMethod<CoarseVector> cPM; cPM(PosdefLdop,c_src);
// std::cout<<GridLogMessage << " CoarseCoarse PowerMethod "<< std::endl;
// PowerMethod<CoarseCoarseVector> ccPM; ccPM(IRLHermOpL2,cc_src);
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
std::cout<<GridLogMessage << "Done "<< std::endl;
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
Grid_finalize();
}

2
tests/solver/Test_dwf_mrhs_cg.cc
View File

@@ -93,7 +93,7 @@ int main (int argc, char ** argv)
for(int s=0;s<nrhs;s++) random(pRNG5,src[s]);
for(int s=0;s<nrhs;s++) result[s]=Zero();
LatticeGaugeField Umu(UGrid); SU3::HotConfiguration(pRNG,Umu);
LatticeGaugeField Umu(UGrid); SU<Nc>::HotConfiguration(pRNG,Umu);
///////////////////////////////////////////////////////////////
// Bounce these fields to disk

4
tests/solver/Test_dwf_mrhs_cg_mpi.cc
View File

@@ -136,11 +136,11 @@ int main (int argc, char ** argv)
std::cout << GridLogMessage << "Intialising 4D RNG "<<std::endl;
pRNG.SeedFixedIntegers(seeds);
std::cout << GridLogMessage << "Intialised 4D RNG "<<std::endl;
SU3::HotConfiguration(pRNG,Umu);
SU<Nc>::HotConfiguration(pRNG,Umu);
std::cout << GridLogMessage << "Intialised the HOT Gauge Field"<<std::endl;
// std::cout << " Site zero "<< Umu[0] <<std::endl;
} else {
SU3::ColdConfiguration(Umu);
SU<Nc>::ColdConfiguration(Umu);
std::cout << GridLogMessage << "Intialised the COLD Gauge Field"<<std::endl;
}
/////////////////

2
tests/solver/Test_dwf_mrhs_cg_mpieo.cc
View File

@@ -87,7 +87,7 @@ int main (int argc, char ** argv)
for(int s=0;s<nrhs;s++) random(pRNG5,src[s]);
for(int s=0;s<nrhs;s++) result[s]=Zero();
LatticeGaugeField Umu(UGrid); SU3::HotConfiguration(pRNG,Umu);
LatticeGaugeField Umu(UGrid); SU<Nc>::HotConfiguration(pRNG,Umu);
/////////////////
// MPI only sends

9
tests/solver/Test_dwf_multigrid.cc
View File

@@ -370,6 +370,11 @@ int main (int argc, char ** argv)
GridCartesian *CoarseCoarse4d = SpaceTimeGrid::makeFourDimGrid(cclatt, GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi());;
GridCartesian *CoarseCoarse5d = SpaceTimeGrid::makeFiveDimGrid(1,CoarseCoarse4d);
GridRedBlackCartesian * Coarse4dRB = SpaceTimeGrid::makeFourDimRedBlackGrid(Coarse4d);
GridRedBlackCartesian * Coarse5dRB = SpaceTimeGrid::makeFiveDimRedBlackGrid(1,Coarse4d);
GridRedBlackCartesian *CoarseCoarse4dRB = SpaceTimeGrid::makeFourDimRedBlackGrid(CoarseCoarse4d);
GridRedBlackCartesian *CoarseCoarse5dRB = SpaceTimeGrid::makeFiveDimRedBlackGrid(1,CoarseCoarse4d);
std::vector<int> seeds4({1,2,3,4});
std::vector<int> seeds5({5,6,7,8});
std::vector<int> cseeds({5,6,7,8});
@@ -434,8 +439,8 @@ int main (int argc, char ** argv)
std::cout<<GridLogMessage << "Building coarse representation of Indef operator" <<std::endl;
std::cout<<GridLogMessage << "**************************************************"<< std::endl;
Level1Op LDOp(*Coarse5d,1); LDOp.CoarsenOperator(FGrid,HermIndefOp,Aggregates);
Level1Op LDOpPV(*Coarse5d,1); LDOpPV.CoarsenOperator(FGrid,HermIndefOpPV,Aggregates);
Level1Op LDOp(*Coarse5d,*Coarse5dRB,1); LDOp.CoarsenOperator(FGrid,HermIndefOp,Aggregates);
Level1Op LDOpPV(*Coarse5d,*Coarse5dRB,1); LDOpPV.CoarsenOperator(FGrid,HermIndefOpPV,Aggregates);
std::cout<<GridLogMessage << "**************************************************"<< std::endl;

2
tests/solver/Test_dwf_qmr_unprec.cc
View File

@@ -51,7 +51,7 @@ int main (int argc, char ** argv)
LatticeFermion src(FGrid); random(RNG5,src);
LatticeFermion result(FGrid); result=Zero();
LatticeGaugeField Umu(UGrid); SU3::HotConfiguration(RNG4,Umu);
LatticeGaugeField Umu(UGrid); SU<Nc>::HotConfiguration(RNG4,Umu);
std::vector<LatticeColourMatrix> U(4,UGrid);

4
tests/solver/Test_hw_multigrid.cc
View File

@@ -274,6 +274,8 @@ int main (int argc, char ** argv)
GridCartesian *Coarse4d = SpaceTimeGrid::makeFourDimGrid(clatt, GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi());;
GridCartesian *Coarse5d = SpaceTimeGrid::makeFiveDimGrid(Ls,Coarse4d);
GridRedBlackCartesian * Coarse4dRB = SpaceTimeGrid::makeFourDimRedBlackGrid(Coarse4d);
GridRedBlackCartesian * Coarse5dRB = SpaceTimeGrid::makeFiveDimRedBlackGrid(1,Coarse4d);
std::vector<int> seeds({1,2,3,4});
GridParallelRNG RNG5(FGrid); RNG5.SeedFixedIntegers(seeds);
@@ -335,7 +337,7 @@ int main (int argc, char ** argv)
NonHermitianLinearOperator<DomainWallFermionR,LatticeFermion> LinOpDwf(Ddwf);
Level1Op LDOp (*Coarse5d,0);
Level1Op LDOp (*Coarse5d,*Coarse5dRB,0);
std::cout<<GridLogMessage << " Callinig Coarsen the operator " <<std::endl;
LDOp.CoarsenOperator(FGrid,LinOpDwf,Aggregates5D);

1287
tests/solver/Test_hw_multigrid_mixed_48.cc Normal file
View File

File diff suppressed because it is too large Load Diff

1326
tests/solver/Test_hw_multigrid_mixed_48_rb.cc Normal file
View File

File diff suppressed because it is too large Load Diff

2
tests/solver/Test_mobius_bcg.cc
View File

@@ -128,7 +128,7 @@ int main (int argc, char ** argv)
std::cout << GridLogMessage << "Intialising 4D RNG "<<std::endl;
pRNG.SeedFixedIntegers(seeds);
std::cout << GridLogMessage << "Intialised 4D RNG "<<std::endl;
SU3::HotConfiguration(pRNG,Umu);
SU<Nc>::HotConfiguration(pRNG,Umu);
std::cout << GridLogMessage << "Intialised the HOT Gauge Field"<<std::endl;
}

4
tests/solver/Test_mobius_bcg_nosplit.cc
View File

@@ -93,10 +93,10 @@ int main (int argc, char ** argv)
GridParallelRNG pRNG(UGrid );
pRNG.SeedFixedIntegers(seeds);
SU3::HotConfiguration(pRNG,Umu);
SU<Nc>::HotConfiguration(pRNG,Umu);
std::cout << GridLogMessage << "Intialised the HOT Gauge Field"<<std::endl;
} else {
SU3::ColdConfiguration(Umu);
SU<Nc>::ColdConfiguration(Umu);
std::cout << GridLogMessage << "Intialised the COLD Gauge Field"<<std::endl;
}

4
tests/solver/Test_mobius_bcg_phys_nosplit.cc
View File

@@ -94,10 +94,10 @@ int main (int argc, char ** argv)
GridParallelRNG pRNG(UGrid );
pRNG.SeedFixedIntegers(seeds);
SU3::HotConfiguration(pRNG,Umu);
SU<Nc>::HotConfiguration(pRNG,Umu);
std::cout << GridLogMessage << "Intialised the HOT Gauge Field"<<std::endl;
} else {
SU3::ColdConfiguration(Umu);
SU<Nc>::ColdConfiguration(Umu);
std::cout << GridLogMessage << "Intialised the COLD Gauge Field"<<std::endl;
}

Some files were not shown because too many files have changed in this diff Show More