Some small steps towards a multigrid

lib/Make.inc

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

lib/algorithms/CoarsenedMatrix.h

@@ -85,10 +85,6 @@ namespace Grid {
     void Orthogonalise(void){
       CoarseScalar InnerProd(CoarseGrid); 
       blockOrthogonalise(InnerProd,subspace);
-#if 1
-      //      CheckOrthogonal();
-#endif
-
     } 
     void CheckOrthogonal(void){
       CoarseVector iProj(CoarseGrid); 
@@ -125,7 +121,7 @@ namespace Grid {
 
       RealD scale;
 
-      ConjugateGradient<FineField> CG(1.0e-4,10000);
+      ConjugateGradient<FineField> CG(1.0e-3,10000);
       FineField noise(FineGrid);
       FineField Mn(FineGrid);
 

lib/algorithms/iterative/ConjugateGradient.h

@@ -15,7 +15,7 @@ public:
     Integer MaxIterations;
     int verbose;
     ConjugateGradient(RealD tol,Integer maxit) : Tolerance(tol), MaxIterations(maxit) { 
-      verbose=0;
+      verbose=1;
     };
 
 

lib/algorithms/iterative/PrecGeneralisedConjugateResidual.h

@@ -48,11 +48,11 @@ namespace Grid {
 	if(cp<rsq) {
 	  Linop.HermOp(psi,r);
 	  axpy(r,-1.0,src,r);
-	  RealD true_resid = norm2(r);
+	  RealD tr = norm2(r);
 	  std::cout<<"PrecGeneralisedConjugateResidual: Converged on iteration " <<steps
 		   << " computed residual "<<sqrt(cp/ssq)
-	           << " true residual "<<true_resid
+	           << " true residual "    <<sqrt(tr/ssq)
-	           << " target "       <<Tolerance <<std::endl;
+	           << " target "           <<Tolerance <<std::endl;
 	  return;
 	}
 

lib/qcd/action/fermion/WilsonFermion5D.cc

@@ -87,8 +87,8 @@ void WilsonFermion5D::DhopDir(const LatticeFermion &in, LatticeFermion &out,int
 {
   int dir = dir5-1; // Maps to the ordering above in "directions" that is passed to stencil
                     // we drop off the innermost fifth dimension
-  assert( (disp==1)||(disp==-1) );
+  //  assert( (disp==1)||(disp==-1) );
-  assert( (dir>=0)&&(dir<4) ); //must do x,y,z or t;
+  //  assert( (dir>=0)&&(dir<4) ); //must do x,y,z or t;
 
   WilsonCompressor compressor(DaggerNo);
   Stencil.HaloExchange<vSpinColourVector,vHalfSpinColourVector,WilsonCompressor>(in,comm_buf,compressor);
@@ -100,7 +100,7 @@ void WilsonFermion5D::DhopDir(const LatticeFermion &in, LatticeFermion &out,int
   assert(dirdisp<=7);
   assert(dirdisp>=0);
 
-PARALLEL_FOR_LOOP
+//PARALLEL_FOR_LOOP
   for(int ss=0;ss<Umu._grid->oSites();ss++){
     for(int s=0;s<Ls;s++){
       int sU=ss;
@@ -114,7 +114,7 @@ void WilsonFermion5D::DhopInternal(CartesianStencil & st, LebesgueOrder &lo,
 				   LatticeDoubledGaugeField & U,
 			   const LatticeFermion &in, LatticeFermion &out,int dag)
 {
-  assert((dag==DaggerNo) ||(dag==DaggerYes));
+  //  assert((dag==DaggerNo) ||(dag==DaggerYes));
 
   WilsonCompressor compressor(dag);
 
@@ -127,29 +127,32 @@ void WilsonFermion5D::DhopInternal(CartesianStencil & st, LebesgueOrder &lo,
   // - 8 linear access unit stride streams per thread for Fermion for hw prefetchable.
   if ( dag == DaggerYes ) {
     if( HandOptDslash ) {
-PARALLEL_FOR_LOOP
+#pragma parallel for
       for(int ss=0;ss<U._grid->oSites();ss++){
-	for(int s=0;s<Ls;s++){
+	{
-	  //int sU=lo.Reorder(ss);
+	  for(int s=0;s<Ls;s++){
-	  int sU=ss;
+	    int sU=ss;
-	  int sF = s+Ls*sU;
+	    int sF = s+Ls*sU;
-	  DiracOptHand::DhopSiteDag(st,U,comm_buf,sF,sU,in,out);
+	    DiracOptHand::DhopSiteDag(st,U,comm_buf,sF,sU,in,out);
+	  }
 	}
       }
     } else { 
-PARALLEL_FOR_LOOP
+#pragma parallel for
       for(int ss=0;ss<U._grid->oSites();ss++){
-	for(int s=0;s<Ls;s++){
+	{
-	  //	  int sU=lo.Reorder(ss);
+	  int sd;
-	  int sU=ss;
+	  for(sd=0;sd<Ls;sd++){
-	  int sF = s+Ls*sU;
+	    int sU=ss;
-	  DiracOpt::DhopSiteDag(st,U,comm_buf,sF,sU,in,out);
+	    int sF = sd+Ls*sU;
+	    DiracOpt::DhopSiteDag(st,U,comm_buf,sF,sU,in,out);
+	  }
 	}
       }
     }
   } else {
     if( HandOptDslash ) {
-PARALLEL_FOR_LOOP
+#pragma parallel for
       for(int ss=0;ss<U._grid->oSites();ss++){
 	for(int s=0;s<Ls;s++){
 	  //	  int sU=lo.Reorder(ss);
@@ -160,7 +163,7 @@ PARALLEL_FOR_LOOP
       }
 
     } else { 
-PARALLEL_FOR_LOOP
+#pragma parallel for
       for(int ss=0;ss<U._grid->oSites();ss++){
 	for(int s=0;s<Ls;s++){
 	  //	  int sU=lo.Reorder(ss);

lib/qcd/utils/SUn.h

@@ -166,7 +166,7 @@ public:
     su2SubGroupIndex(i0,i1,su2_index);
 
 PARALLEL_FOR_LOOP
-    for(int ss=0;ss!=grid->oSites();ss++){
+    for(int ss=0;ss<grid->oSites();ss++){
       subgroup._odata[ss]()()(0,0) = source._odata[ss]()()(i0,i0);
       subgroup._odata[ss]()()(0,1) = source._odata[ss]()()(i0,i1);
       subgroup._odata[ss]()()(1,0) = source._odata[ss]()()(i1,i0);
@@ -201,7 +201,7 @@ PARALLEL_FOR_LOOP
 
     dest = 1.0; // start out with identity
 PARALLEL_FOR_LOOP
-    for(int ss=0;ss!=grid->oSites();ss++){
+    for(int ss=0;ss<grid->oSites();ss++){
       dest._odata[ss]()()(i0,i0) = subgroup._odata[ss]()()(0,0);
       dest._odata[ss]()()(i0,i1) = subgroup._odata[ss]()()(0,1);
       dest._odata[ss]()()(i1,i0) = subgroup._odata[ss]()()(1,0);

tests/Make.inc

@@ -1,15 +1,11 @@
 
-bin_PROGRAMS = Test_GaugeAction Test_Metropolis Test_cayley_cg Test_cayley_coarsen_support Test_cayley_even_odd Test_cayley_ldop_cg Test_cayley_ldop_cr Test_cf_coarsen_support Test_cf_cr_unprec Test_contfrac_cg Test_contfrac_even_odd Test_cshift Test_cshift_red_black Test_dwf_cg_prec Test_dwf_cg_schur Test_dwf_cg_unprec Test_dwf_cr_unprec Test_dwf_even_odd Test_dwf_fpgcr Test_gamma Test_lie_generators Test_main Test_multishift_sqrt Test_nersc_io Test_quenched_update Test_remez Test_rng Test_rng_fixed Test_simd Test_stencil Test_wilson_cg_prec Test_wilson_cg_schur Test_wilson_cg_unprec Test_wilson_cr_unprec Test_wilson_even_odd
+bin_PROGRAMS = Test_GaugeAction Test_cayley_cg Test_cayley_coarsen_support Test_cayley_even_odd Test_cayley_ldop_cr Test_cf_coarsen_support Test_cf_cr_unprec Test_contfrac_cg Test_contfrac_even_odd Test_cshift Test_cshift_red_black Test_dwf_cg_prec Test_dwf_cg_schur Test_dwf_cg_unprec Test_dwf_cr_unprec Test_dwf_even_odd Test_dwf_fpgcr Test_dwf_hdcr Test_gamma Test_lie_generators Test_main Test_multishift_sqrt Test_nersc_io Test_quenched_update Test_remez Test_rng Test_rng_fixed Test_simd Test_stencil Test_wilson_cg_prec Test_wilson_cg_schur Test_wilson_cg_unprec Test_wilson_cr_unprec Test_wilson_even_odd
 
 
 Test_GaugeAction_SOURCES=Test_GaugeAction.cc
 Test_GaugeAction_LDADD=-lGrid
 
 
-Test_Metropolis_SOURCES=Test_Metropolis.cc
-Test_Metropolis_LDADD=-lGrid
-
-
 Test_cayley_cg_SOURCES=Test_cayley_cg.cc
 Test_cayley_cg_LDADD=-lGrid
 
@@ -22,10 +18,6 @@ Test_cayley_even_odd_SOURCES=Test_cayley_even_odd.cc
 Test_cayley_even_odd_LDADD=-lGrid
 
 
-Test_cayley_ldop_cg_SOURCES=Test_cayley_ldop_cg.cc
-Test_cayley_ldop_cg_LDADD=-lGrid
-
-
 Test_cayley_ldop_cr_SOURCES=Test_cayley_ldop_cr.cc
 Test_cayley_ldop_cr_LDADD=-lGrid
 
@@ -78,6 +70,10 @@ Test_dwf_fpgcr_SOURCES=Test_dwf_fpgcr.cc
 Test_dwf_fpgcr_LDADD=-lGrid
 
 
+Test_dwf_hdcr_SOURCES=Test_dwf_hdcr.cc
+Test_dwf_hdcr_LDADD=-lGrid
+
+
 Test_gamma_SOURCES=Test_gamma.cc
 Test_gamma_LDADD=-lGrid
 

tests/Sqrt.gnu

@@ -1,2 +0,0 @@
-f(x) = 6.81384+(-2.34645e-06/(x+0.000228091))+(-1.51593e-05/(x+0.00112084))+(-6.89254e-05/(x+0.003496))+(-0.000288983/(x+0.00954309))+(-0.00119277/(x+0.024928))+(-0.0050183/(x+0.0646627))+(-0.0226449/(x+0.171576))+(-0.123767/(x+0.491792))+(-1.1705/(x+1.78667))+(-102.992/(x+18.4866));
-f(x) = 0.14676+(0.00952992/(x+5.40933e-05))+(0.0115952/(x+0.000559699))+(0.0161824/(x+0.00203338))+(0.0243252/(x+0.00582831))+(0.0379533/(x+0.0154649))+(0.060699/(x+0.0401156))+(0.100345/(x+0.104788))+(0.178335/(x+0.286042))+(0.381586/(x+0.892189))+(1.42625/(x+4.38422));

tests/Test_cayley_ldop_cg.cc

@@ -1,155 +0,0 @@
-#include <Grid.h>
-#include <qcd/utils/WilsonLoops.h>
-#include <qcd/utils/SUn.h>
-
-using namespace std;
-using namespace Grid;
-using namespace Grid::QCD;
-
-
-template<class d>
-struct scal {
-  d internal;
-};
-
-  Gamma::GammaMatrix Gmu [] = {
-    Gamma::GammaX,
-    Gamma::GammaY,
-    Gamma::GammaZ,
-    Gamma::GammaT
-  };
-
-double lowpass(double x)
-{
-  return pow(x*x+1.0,-2);
-}
-
-int main (int argc, char ** argv)
-{
-  Grid_init(&argc,&argv);
-
-  Chebyshev<LatticeFermion> filter(-150.0, 150.0,16, lowpass);
-  ofstream csv(std::string("filter.dat"),std::ios::out|std::ios::trunc);
-  filter.csv(csv);
-  csv.close();
-
-  const int Ls=8;
-
-  GridCartesian         * UGrid   = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplexF::Nsimd()),GridDefaultMpi());
-  GridRedBlackCartesian * UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
-
-  GridCartesian         * FGrid   = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid);
-  GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid);
-
-  // Construct a coarsened grid
-  std::vector<int> clatt = GridDefaultLatt();
-  for(int d=0;d<clatt.size();d++){
-    clatt[d] = clatt[d]/2;
-  }
-  GridCartesian *Coarse4d =  SpaceTimeGrid::makeFourDimGrid(clatt, GridDefaultSimd(Nd,vComplexF::Nsimd()),GridDefaultMpi());;
-  GridCartesian *Coarse5d =  SpaceTimeGrid::makeFiveDimGrid(1,Coarse4d);
-
-  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);
-  LatticeFermion result(FGrid); result=zero;
-  LatticeFermion    ref(FGrid); ref=zero;
-  LatticeFermion    tmp(FGrid);
-  LatticeFermion    err(FGrid);
-  LatticeGaugeField Umu(UGrid); 
-
-
-
-  //gaussian(RNG4,Umu);
-  //random(RNG4,Umu);
-
-  NerscField header;
-  std::string file("./ckpoint_lat.400");
-  readNerscConfiguration(Umu,header,file);
-  //  SU3::ColdConfiguration(RNG4,Umu);
-  //  SU3::TepidConfiguration(RNG4,Umu);
-  //  SU3::HotConfiguration(RNG4,Umu);
-  //  Umu=zero;
-
-#if 0
-  LatticeColourMatrix U(UGrid);
-  for(int nn=0;nn<Nd;nn++){
-    U=peekIndex<LorentzIndex>(Umu,nn);
-    U=U*adj(U)-1.0;
-    std::cout<<"SU3 test "<<norm2(U)<<std::endl;
-  }
-#endif  
-  RealD mass=0.1;
-  RealD M5=1.5;
-
-  DomainWallFermion Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5);
-  Gamma5R5HermitianLinearOperator<DomainWallFermion,LatticeFermion> HermIndefOp(Ddwf);
-
-  const int nbasis = 8;
-
-#if 0
-  std::vector<LatticeFermion> subspace(nbasis,FGrid);
-  LatticeFermion noise(FGrid);
-  LatticeFermion ms(FGrid);
-  for(int b=0;b<nbasis;b++){
-
-    gaussian(RNG5,noise);
-    RealD scale = pow(norm2(noise),-0.5); 
-    noise=noise*scale;
-
-    HermIndefOp.Op(noise,ms); std::cout << "Noise    "<<b<<" Ms "<<norm2(ms)<< " "<< norm2(noise)<<std::endl;
-
-    //    filter(HermIndefOp,noise,subspace[b]);
-    // inverse iteration
-    MdagMLinearOperator<DomainWallFermion,LatticeFermion> HermDefOp(Ddwf);
-    ConjugateGradient<LatticeFermion> CG(1.0e-4,10000);
-
-    for(int i=0;i<1;i++){
-
-      CG(HermDefOp,noise,subspace[b]);
-      noise = subspace[b];
-
-      scale = pow(norm2(noise),-0.5); 
-      noise=noise*scale;
-      HermDefOp.Op(noise,ms); std::cout << "filt    "<<b<<" <u|H|u> "<<norm2(ms)<< " "<< norm2(noise)<<std::endl;
-    }
-
-    subspace[b] = noise;
-    HermIndefOp.Op(subspace[b],ms); std::cout << "Filtered "<<b<<" Ms "<<norm2(ms)<< " "<<norm2(subspace[b]) <<std::endl;
-
-  }
-  std::cout << "Computed randoms"<< std::endl;
-#else
-  std::cout<<"Calling Aggregation class" <<std::endl;
-  MdagMLinearOperator<DomainWallFermion,LatticeFermion> HermDefOp(Ddwf);
-  typedef Aggregation<vSpinColourVector,vTComplex,nbasis> Subspace;
-  Subspace Aggregates(Coarse5d,FGrid);
-  Aggregates.CreateSubspace(RNG5,HermDefOp);
-  std::cout << "Called aggregation class"<< std::endl;
-#endif
-
-
-  typedef CoarsenedMatrix<vSpinColourVector,vTComplex,nbasis> LittleDiracOperator;
-  typedef LittleDiracOperator::CoarseVector CoarseVector;
-
-  LittleDiracOperator LittleDiracOp(*Coarse5d);
-  LittleDiracOp.CoarsenOperator(FGrid,HermIndefOp,Aggregates);
-  
-  CoarseVector c_src (Coarse5d);
-  CoarseVector c_res (Coarse5d);
-  gaussian(CRNG,c_src);
-  c_res=zero;
-
-  std::cout << "Solving CG on coarse space "<< std::endl;
-  MdagMLinearOperator<LittleDiracOperator,CoarseVector> PosdefLdop(LittleDiracOp);
-  ConjugateGradient<CoarseVector> CG(1.0e-6,10000);
-  CG(PosdefLdop,c_src,c_res);
-
-  std::cout << "Done "<< std::endl;
-  Grid_finalize();
-}

tests/Test_cayley_ldop_cr.cc

@@ -1,38 +1,13 @@
 #include <Grid.h>
-#include <qcd/utils/WilsonLoops.h>
-#include <qcd/utils/SUn.h>
 
 using namespace std;
 using namespace Grid;
 using namespace Grid::QCD;
 
-
-template<class d>
-struct scal {
-  d internal;
-};
-
-  Gamma::GammaMatrix Gmu [] = {
-    Gamma::GammaX,
-    Gamma::GammaY,
-    Gamma::GammaZ,
-    Gamma::GammaT
-  };
-
-double lowpass(double x)
-{
-  return pow(x*x+1.0,-2);
-}
-
 int main (int argc, char ** argv)
 {
   Grid_init(&argc,&argv);
 
-  Chebyshev<LatticeFermion> filter(-150.0, 150.0,16, lowpass);
-  ofstream csv(std::string("filter.dat"),std::ios::out|std::ios::trunc);
-  filter.csv(csv);
-  csv.close();
-
   const int Ls=8;
 
   GridCartesian         * UGrid   = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplexF::Nsimd()),GridDefaultMpi());
@@ -41,7 +16,9 @@ int main (int argc, char ** argv)
   GridCartesian         * FGrid   = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid);
   GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid);
 
-  // Construct a coarsened grid
+  ///////////////////////////////////////////////////
+  // Construct a coarsened grid; utility for this?
+  ///////////////////////////////////////////////////
   std::vector<int> clatt = GridDefaultLatt();
   for(int d=0;d<clatt.size();d++){
     clatt[d] = clatt[d]/2;
@@ -63,80 +40,38 @@ int main (int argc, char ** argv)
   LatticeFermion    err(FGrid);
   LatticeGaugeField Umu(UGrid); 
 
-
-
-  //gaussian(RNG4,Umu);
-  //random(RNG4,Umu);
-
   NerscField header;
   std::string file("./ckpoint_lat.400");
   readNerscConfiguration(Umu,header,file);
+
   //  SU3::ColdConfiguration(RNG4,Umu);
   //  SU3::TepidConfiguration(RNG4,Umu);
   //  SU3::HotConfiguration(RNG4,Umu);
   //  Umu=zero;
 
-#if 0
-  LatticeColourMatrix U(UGrid);
-  for(int nn=0;nn<Nd;nn++){
-    U=peekIndex<LorentzIndex>(Umu,nn);
-    U=U*adj(U)-1.0;
-    std::cout<<"SU3 test "<<norm2(U)<<std::endl;
-  }
-#endif  
   RealD mass=0.1;
   RealD M5=1.5;
 
+  std::cout << "**************************************************"<< std::endl;
+  std::cout << "Building g5R5 hermitian DWF operator" <<std::endl;
+  std::cout << "**************************************************"<< std::endl;
   DomainWallFermion Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5);
   Gamma5R5HermitianLinearOperator<DomainWallFermion,LatticeFermion> HermIndefOp(Ddwf);
 
   const int nbasis = 8;
 
-#if 0
+  typedef Aggregation<vSpinColourVector,vTComplex,nbasis>     Subspace;
-  std::vector<LatticeFermion> subspace(nbasis,FGrid);
+  typedef CoarsenedMatrix<vSpinColourVector,vTComplex,nbasis> LittleDiracOperator;
-  LatticeFermion noise(FGrid);
+  typedef LittleDiracOperator::CoarseVector                   CoarseVector;
-  LatticeFermion ms(FGrid);
-  for(int b=0;b<nbasis;b++){
 
-    gaussian(RNG5,noise);
+  std::cout << "**************************************************"<< std::endl;
-    RealD scale = pow(norm2(noise),-0.5); 
+  std::cout << "Calling Aggregation class to build subspace" <<std::endl;
-    noise=noise*scale;
+  std::cout << "**************************************************"<< std::endl;
-
-    HermIndefOp.Op(noise,ms); std::cout << "Noise    "<<b<<" Ms "<<norm2(ms)<< " "<< norm2(noise)<<std::endl;
-
-    //    filter(HermIndefOp,noise,subspace[b]);
-    // inverse iteration
-    MdagMLinearOperator<DomainWallFermion,LatticeFermion> HermDefOp(Ddwf);
-    ConjugateGradient<LatticeFermion> CG(1.0e-4,10000);
-
-    for(int i=0;i<1;i++){
-
-      CG(HermDefOp,noise,subspace[b]);
-      noise = subspace[b];
-
-      scale = pow(norm2(noise),-0.5); 
-      noise=noise*scale;
-      HermDefOp.Op(noise,ms); std::cout << "filt    "<<b<<" <u|H|u> "<<norm2(ms)<< " "<< norm2(noise)<<std::endl;
-    }
-
-    subspace[b] = noise;
-    HermIndefOp.Op(subspace[b],ms); std::cout << "Filtered "<<b<<" Ms "<<norm2(ms)<< " "<<norm2(subspace[b]) <<std::endl;
-
-  }
-  std::cout << "Computed randoms"<< std::endl;
-#else
-  std::cout<<"Calling Aggregation class" <<std::endl;
   MdagMLinearOperator<DomainWallFermion,LatticeFermion> HermDefOp(Ddwf);
-  typedef Aggregation<vSpinColourVector,vTComplex,nbasis> Subspace;
   Subspace Aggregates(Coarse5d,FGrid);
   Aggregates.CreateSubspace(RNG5,HermDefOp);
-  std::cout << "Called aggregation class"<< std::endl;
-#endif
 
 
-  typedef CoarsenedMatrix<vSpinColourVector,vTComplex,nbasis> LittleDiracOperator;
-  typedef LittleDiracOperator::CoarseVector CoarseVector;
-
   LittleDiracOperator LittleDiracOp(*Coarse5d);
   LittleDiracOp.CoarsenOperator(FGrid,HermIndefOp,Aggregates);
   
@@ -145,18 +80,22 @@ int main (int argc, char ** argv)
   gaussian(CRNG,c_src);
   c_res=zero;
 
-
+  std::cout << "**************************************************"<< std::endl;
-  std::cout << "Solving CG on coarse space "<< std::endl;
+  std::cout << "Solving mdagm-CG on coarse space "<< std::endl;
+  std::cout << "**************************************************"<< std::endl;
   MdagMLinearOperator<LittleDiracOperator,CoarseVector> PosdefLdop(LittleDiracOp);
   ConjugateGradient<CoarseVector> CG(1.0e-6,10000);
   CG(PosdefLdop,c_src,c_res);
 
-
+  std::cout << "**************************************************"<< std::endl;
-  std::cout << "Solving MCR on coarse space "<< std::endl;
+  std::cout << "Solving indef-MCR on coarse space "<< std::endl;
+  std::cout << "**************************************************"<< std::endl;
   HermitianLinearOperator<LittleDiracOperator,CoarseVector> HermIndefLdop(LittleDiracOp);
   ConjugateResidual<CoarseVector> MCR(1.0e-6,10000);
   MCR(HermIndefLdop,c_src,c_res);
 
+  std::cout << "**************************************************"<< std::endl;
   std::cout << "Done "<< std::endl;
+  std::cout << "**************************************************"<< std::endl;
   Grid_finalize();
 }

tests/Test_dwf_hdcr.cc

@@ -0,0 +1,188 @@
+#include <Grid.h>
+#include <algorithms/iterative/PrecGeneralisedConjugateResidual.h>
+
+using namespace std;
+using namespace Grid;
+using namespace Grid::QCD;
+
+template<class Fobj,class CComplex,int nbasis>
+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>::siteVector     siteVector;
+  typedef typename Aggregation<Fobj,CComplex,nbasis>::CoarseScalar CoarseScalar;
+  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;
+
+  Aggregates     & _Aggregates;
+  CoarseOperator & _CoarseOperator;
+  FineOperator   & _FineOperator;
+
+  // Constructor
+  MultiGridPreconditioner(Aggregates &Agg, CoarseOperator &Coarse, FineOperator &Fine) 
+    : _Aggregates(Agg),
+      _CoarseOperator(Coarse),
+      _FineOperator(Fine)
+  {
+  }
+
+  void operator()(const FineField &in, FineField & out) {
+
+    FineField Min(in._grid);
+
+    CoarseVector Csrc(_CoarseOperator.Grid());
+    CoarseVector Ctmp(_CoarseOperator.Grid());
+    CoarseVector Csol(_CoarseOperator.Grid());
+
+    // Monitor completeness of low mode space
+    _Aggregates.ProjectToSubspace  (Csrc,in);
+    _Aggregates.PromoteFromSubspace(Csrc,out);
+    std::cout<<"Completeness: "<<std::sqrt(norm2(out)/norm2(in))<<std::endl;
+
+    ConjugateResidual<FineField>    MCR(1.0e-2,1000);
+    ConjugateGradient<CoarseVector>  CG(1.0e-2,10000);
+
+    // [PTM+Q] in = [1 - Q A] M in + Q in = Min + Q [ in -A Min]
+    // Smoothing step, followed by coarse grid correction
+
+    MCR(_FineOperator,in,Min);
+    _FineOperator.Op(Min,out);
+    out = in -out; // out = in - A Min
+
+    MdagMLinearOperator<CoarseOperator,CoarseVector>     MdagMOp(_CoarseOperator);
+    HermitianLinearOperator<CoarseOperator,CoarseVector> HermOp(_CoarseOperator);
+    Csol=zero;
+    _Aggregates.ProjectToSubspace  (Csrc,out);
+    HermOp.AdjOp(Csrc,Ctmp);// Normal equations
+    CG(MdagMOp  ,Ctmp,Csol);
+    _Aggregates.PromoteFromSubspace(Csol,out);
+
+    out = Min + out;;
+  }
+
+};
+
+int main (int argc, char ** argv)
+{
+  Grid_init(&argc,&argv);
+
+  const int Ls=8;
+
+  GridCartesian         * UGrid   = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplexF::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> clatt = GridDefaultLatt();
+  for(int d=0;d<clatt.size();d++){
+    clatt[d] = clatt[d]/4;
+  }
+  GridCartesian *Coarse4d =  SpaceTimeGrid::makeFourDimGrid(clatt, GridDefaultSimd(Nd,vComplexF::Nsimd()),GridDefaultMpi());;
+  GridCartesian *Coarse5d =  SpaceTimeGrid::makeFiveDimGrid(1,Coarse4d);
+
+  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);
+  LatticeFermion result(FGrid); result=zero;
+  LatticeFermion    ref(FGrid); ref=zero;
+  LatticeFermion    tmp(FGrid);
+  LatticeFermion    err(FGrid);
+  LatticeGaugeField Umu(UGrid); 
+
+  NerscField header;
+  std::string file("./ckpoint_lat.4000");
+  readNerscConfiguration(Umu,header,file);
+
+  //  SU3::ColdConfiguration(RNG4,Umu);
+  //  SU3::TepidConfiguration(RNG4,Umu);
+  //  SU3::HotConfiguration(RNG4,Umu);
+  //  Umu=zero;
+
+  RealD mass=0.04;
+  RealD M5=1.8;
+
+  std::cout << "**************************************************"<< std::endl;
+  std::cout << "Building g5R5 hermitian DWF operator" <<std::endl;
+  std::cout << "**************************************************"<< std::endl;
+  DomainWallFermion Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5);
+
+  const int nbasis = 4;
+
+  typedef Aggregation<vSpinColourVector,vTComplex,nbasis>              Subspace;
+  typedef CoarsenedMatrix<vSpinColourVector,vTComplex,nbasis>          CoarseOperator;
+  typedef CoarseOperator::CoarseVector                                 CoarseVector;
+
+  std::cout << "**************************************************"<< std::endl;
+  std::cout << "Calling Aggregation class to build subspace" <<std::endl;
+  std::cout << "**************************************************"<< std::endl;
+  MdagMLinearOperator<DomainWallFermion,LatticeFermion> HermDefOp(Ddwf);
+  Subspace Aggregates(Coarse5d,FGrid);
+  Aggregates.CreateSubspace(RNG5,HermDefOp);
+
+  std::cout << "**************************************************"<< std::endl;
+  std::cout << "Building coarse representation of Indef operator" <<std::endl;
+  std::cout << "**************************************************"<< std::endl;
+  Gamma5R5HermitianLinearOperator<DomainWallFermion,LatticeFermion> HermIndefOp(Ddwf);
+  CoarsenedMatrix<vSpinColourVector,vTComplex,nbasis> LDOp(*Coarse5d);
+  LDOp.CoarsenOperator(FGrid,HermIndefOp,Aggregates);
+
+  std::cout << "**************************************************"<< std::endl;
+  std::cout << "Testing some coarse space solvers  " <<std::endl;
+  std::cout << "**************************************************"<< std::endl;
+  CoarseVector c_src (Coarse5d);
+  CoarseVector c_res (Coarse5d);
+  gaussian(CRNG,c_src);
+  c_res=zero;
+
+  std::cout << "**************************************************"<< std::endl;
+  std::cout << "Solving posdef-CG on coarse space "<< std::endl;
+  std::cout << "**************************************************"<< std::endl;
+  MdagMLinearOperator<CoarseOperator,CoarseVector> PosdefLdop(LDOp);
+  ConjugateGradient<CoarseVector> CG(1.0e-6,10000);
+  CG(PosdefLdop,c_src,c_res);
+
+  std::cout << "**************************************************"<< std::endl;
+  std::cout << "Solving indef-MCR on coarse space "<< std::endl;
+  std::cout << "**************************************************"<< std::endl;
+  HermitianLinearOperator<CoarseOperator,CoarseVector> HermIndefLdop(LDOp);
+  ConjugateResidual<CoarseVector> MCR(1.0e-6,10000);
+  //MCR(HermIndefLdop,c_src,c_res);
+
+  std::cout << "**************************************************"<< std::endl;
+  std::cout << "Building deflation preconditioner "<< std::endl;
+  std::cout << "**************************************************"<< std::endl;
+
+  MultiGridPreconditioner <vSpinColourVector,vTComplex,nbasis> Precon(Aggregates, LDOp,HermIndefOp);
+
+  std::cout << "**************************************************"<< std::endl;
+  std::cout << "Building a one level PGCR "<< std::endl;
+  std::cout << "**************************************************"<< std::endl;
+  TrivialPrecon<LatticeFermion> simple;
+  PrecGeneralisedConjugateResidual<LatticeFermion> GCR(1.0e-6,10000,simple,8,64);
+  GCR(HermIndefOp,src,result);
+
+  std::cout << "**************************************************"<< std::endl;
+  std::cout << "Building a two level PGCR "<< std::endl;
+  std::cout << "**************************************************"<< std::endl;
+  PrecGeneralisedConjugateResidual<LatticeFermion> PGCR(1.0e-6,10000,Precon,8,64);
+  PGCR(HermIndefOp,src,result);
+
+  std::cout << "**************************************************"<< std::endl;
+  std::cout << "Done "<< std::endl;
+  std::cout << "**************************************************"<< std::endl;
+  Grid_finalize();
+}