aarch64 libunwind compatibility fix

Grid/Makefile.am

@@ -54,24 +54,22 @@ Version.h: version-cache
 include Make.inc
 include Eigen.inc
 
-if BUILD_FERMION_INSTANTIATIONS
+extra_sources+=$(WILS_FERMION_FILES)
-  extra_sources+=$(WILS_FERMION_FILES)
+extra_sources+=$(STAG_FERMION_FILES)
-  extra_sources+=$(STAG_FERMION_FILES)
 if BUILD_ZMOBIUS
-    extra_sources+=$(ZWILS_FERMION_FILES)
+  extra_sources+=$(ZWILS_FERMION_FILES)
 endif
 if BUILD_GPARITY
-    extra_sources+=$(GP_FERMION_FILES)
+  extra_sources+=$(GP_FERMION_FILES)
 endif
 if BUILD_FERMION_REPS
-    extra_sources+=$(ADJ_FERMION_FILES)
+  extra_sources+=$(ADJ_FERMION_FILES)
-    extra_sources+=$(TWOIND_FERMION_FILES)
+  extra_sources+=$(TWOIND_FERMION_FILES)
 endif
 if BUILD_SP
     extra_sources+=$(SP_FERMION_FILES)
 if BUILD_FERMION_REPS
-      extra_sources+=$(SP_TWOIND_FERMION_FILES)
+    extra_sources+=$(SP_TWOIND_FERMION_FILES)
-endif
 endif
 endif
 

Grid/algorithms/blas/BatchedBlas.h

@@ -28,7 +28,6 @@ Author: Peter Boyle <pboyle@bnl.gov>
 #pragma once
 
 #ifdef GRID_HIP
-#include <hip/hip_version.h>
 #include <hipblas/hipblas.h>
 #endif
 #ifdef GRID_CUDA
@@ -256,29 +255,16 @@ public:
     if ( OpB == GridBLAS_OP_N ) hOpB = HIPBLAS_OP_N;
     if ( OpB == GridBLAS_OP_T ) hOpB = HIPBLAS_OP_T;
     if ( OpB == GridBLAS_OP_C ) hOpB = HIPBLAS_OP_C;
-#if defined(HIP_VERSION_MAJOR) && (HIP_VERSION_MAJOR >=7)
     auto err = hipblasZgemmBatched(gridblasHandle,
 				   hOpA,
 				   hOpB,
 				   m,n,k,
-				   (hipDoubleComplex *) &alpha_p[0],
+				   (hipblasDoubleComplex *) &alpha_p[0],
-				   (hipDoubleComplex **)&Amk[0], lda,
+				   (hipblasDoubleComplex **)&Amk[0], lda,
-				   (hipDoubleComplex **)&Bkn[0], ldb,
+				   (hipblasDoubleComplex **)&Bkn[0], ldb,
-				   (hipDoubleComplex *) &beta_p[0],
+				   (hipblasDoubleComplex *) &beta_p[0],
-				   (hipDoubleComplex **)&Cmn[0], ldc,
+				   (hipblasDoubleComplex **)&Cmn[0], ldc,
 				   batchCount);
-#else
-    auto err = hipblasZgemmBatched(gridblasHandle,
-                                   hOpA,
-                                   hOpB,
-                                   m,n,k,
-                                   (hipblasDoubleComplex *) &alpha_p[0],
-                                   (hipblasDoubleComplex **)&Amk[0], lda,
-                                   (hipblasDoubleComplex **)&Bkn[0], ldb,
-                                   (hipblasDoubleComplex *) &beta_p[0],
-                                   (hipblasDoubleComplex **)&Cmn[0], ldc,
-                                   batchCount);
-#endif
     //	 std::cout << " hipblas return code " <<(int)err<<std::endl;
     GRID_ASSERT(err==HIPBLAS_STATUS_SUCCESS);
 #endif
@@ -517,30 +503,17 @@ public:
     if ( OpB == GridBLAS_OP_N ) hOpB = HIPBLAS_OP_N;
     if ( OpB == GridBLAS_OP_T ) hOpB = HIPBLAS_OP_T;
     if ( OpB == GridBLAS_OP_C ) hOpB = HIPBLAS_OP_C;
-#if defined(HIP_VERSION_MAJOR) && (HIP_VERSION_MAJOR >=7)
     auto err = hipblasCgemmBatched(gridblasHandle,
 				   hOpA,
 				   hOpB,
 				   m,n,k,
-				   (hipComplex *) &alpha_p[0],
+				   (hipblasComplex *) &alpha_p[0],
-				   (hipComplex **)&Amk[0], lda,
+				   (hipblasComplex **)&Amk[0], lda,
-				   (hipComplex **)&Bkn[0], ldb,
+				   (hipblasComplex **)&Bkn[0], ldb,
-				   (hipComplex *) &beta_p[0],
+				   (hipblasComplex *) &beta_p[0],
-				   (hipComplex **)&Cmn[0], ldc,
+				   (hipblasComplex **)&Cmn[0], ldc,
 				   batchCount);
-#else
-    auto err = hipblasCgemmBatched(gridblasHandle,
-                                   hOpA,
-                                   hOpB,
-                                   m,n,k,
-                                   (hipblasComplex *) &alpha_p[0],
-                                   (hipblasComplex **)&Amk[0], lda,
-                                   (hipblasComplex **)&Bkn[0], ldb,
-                                   (hipblasComplex *) &beta_p[0],
-                                   (hipblasComplex **)&Cmn[0], ldc,
-                                   batchCount);
 
-#endif
     GRID_ASSERT(err==HIPBLAS_STATUS_SUCCESS);
 #endif
 #ifdef GRID_CUDA
@@ -1121,19 +1094,11 @@ public:
     GRID_ASSERT(info.size()==batchCount);
 
 #ifdef GRID_HIP
-#if defined(HIP_VERSION_MAJOR) && (HIP_VERSION_MAJOR >=7)
     auto err = hipblasZgetrfBatched(gridblasHandle,(int)n,
-				    (hipDoubleComplex **)&Ann[0], (int)n,
+				    (hipblasDoubleComplex **)&Ann[0], (int)n,
 				    (int*) &ipiv[0],
 				    (int*) &info[0],
 				    (int)batchCount);
-#else
-    auto err = hipblasZgetrfBatched(gridblasHandle,(int)n,
-                                    (hipblasDoubleComplex **)&Ann[0], (int)n,
-                                    (int*) &ipiv[0],
-                                    (int*) &info[0],
-                                    (int)batchCount);
-#endif
     GRID_ASSERT(err==HIPBLAS_STATUS_SUCCESS);
 #endif
 #ifdef GRID_CUDA
@@ -1159,20 +1124,11 @@ public:
     GRID_ASSERT(info.size()==batchCount);
 
 #ifdef GRID_HIP
-#if defined(HIP_VERSION_MAJOR) && (HIP_VERSION_MAJOR >=7)
     auto err = hipblasCgetrfBatched(gridblasHandle,(int)n,
-				    (hipComplex **)&Ann[0], (int)n,
+				    (hipblasComplex **)&Ann[0], (int)n,
 				    (int*) &ipiv[0],
 				    (int*) &info[0],
 				    (int)batchCount);
-#else
-    auto err = hipblasCgetrfBatched(gridblasHandle,(int)n,
-                                    (hipblasComplex **)&Ann[0], (int)n,
-                                    (int*) &ipiv[0],
-                                    (int*) &info[0],
-                                    (int)batchCount);
-#endif
-
     GRID_ASSERT(err==HIPBLAS_STATUS_SUCCESS);
 #endif
 #ifdef GRID_CUDA
@@ -1245,22 +1201,12 @@ public:
     GRID_ASSERT(Cnn.size()==batchCount);
 
 #ifdef GRID_HIP
-#if defined(HIP_VERSION_MAJOR) && (HIP_VERSION_MAJOR >=7)
     auto err = hipblasZgetriBatched(gridblasHandle,(int)n,
-				    (hipDoubleComplex **)&Ann[0], (int)n,
+				    (hipblasDoubleComplex **)&Ann[0], (int)n,
 				    (int*) &ipiv[0],
-				    (hipDoubleComplex **)&Cnn[0], (int)n,
+				    (hipblasDoubleComplex **)&Cnn[0], (int)n,
 				    (int*) &info[0],
 				    (int)batchCount);
-#else
-    auto err = hipblasZgetriBatched(gridblasHandle,(int)n,
-                                    (hipblasDoubleComplex **)&Ann[0], (int)n,
-                                    (int*) &ipiv[0],
-                                    (hipblasDoubleComplex **)&Cnn[0], (int)n,
-                                    (int*) &info[0],
-                                    (int)batchCount);
-
-#endif
     GRID_ASSERT(err==HIPBLAS_STATUS_SUCCESS);
 #endif
 #ifdef GRID_CUDA
@@ -1289,21 +1235,12 @@ public:
     GRID_ASSERT(Cnn.size()==batchCount);
 
 #ifdef GRID_HIP
-#if defined(HIP_VERSION_MAJOR) && (HIP_VERSION_MAJOR >=7)
     auto err = hipblasCgetriBatched(gridblasHandle,(int)n,
-				    (hipComplex **)&Ann[0], (int)n,
+				    (hipblasComplex **)&Ann[0], (int)n,
 				    (int*) &ipiv[0],
-				    (hipComplex **)&Cnn[0], (int)n,
+				    (hipblasComplex **)&Cnn[0], (int)n,
 				    (int*) &info[0],
 				    (int)batchCount);
-#else
-    auto err = hipblasCgetriBatched(gridblasHandle,(int)n,
-                                    (hipblasComplex **)&Ann[0], (int)n,
-                                    (int*) &ipiv[0],
-                                    (hipblasComplex **)&Cnn[0], (int)n,
-                                    (int*) &info[0],
-                                    (int)batchCount);
-#endif
     GRID_ASSERT(err==HIPBLAS_STATUS_SUCCESS);
 #endif
 #ifdef GRID_CUDA

Grid/algorithms/iterative/AdefMrhs.h

@@ -92,8 +92,8 @@ class TwoLevelCGmrhs
   // Vector case
   virtual void operator() (std::vector<Field> &src, std::vector<Field> &x)
   {
-    SolveSingleSystem(src,x);
+    //    SolveSingleSystem(src,x);
-	// SolvePrecBlockCG(src,x);
+    SolvePrecBlockCG(src,x);
   }
 
 ////////////////////////////////////////////////////////////////////////////////////////////////////

Grid/algorithms/multigrid/Aggregates.h

@@ -97,7 +97,7 @@ public:
 
     RealD scale;
 
-    ConjugateGradient<FineField> CG(1.0e-4,2000,false);
+    ConjugateGradient<FineField> CG(1.0e-3,400,false);
     FineField noise(FineGrid);
     FineField Mn(FineGrid);
 
@@ -131,10 +131,7 @@ public:
     RealD scale;
 
     TrivialPrecon<FineField> simple_fine;
-    //    PrecGeneralisedConjugateResidualNonHermitian<FineField> GCR(0.001,10,DiracOp,simple_fine,30,30);
+    PrecGeneralisedConjugateResidualNonHermitian<FineField> GCR(0.001,30,DiracOp,simple_fine,12,12);
-    //    PrecGeneralisedConjugateResidualNonHermitian<FineField> GCR(0.001,10,DiracOp,simple_fine,12,12);
-    //    PrecGeneralisedConjugateResidualNonHermitian<FineField> GCR(0.001,30,DiracOp,simple_fine,12,12);
-    PrecGeneralisedConjugateResidualNonHermitian<FineField> GCR(0.001,30,DiracOp,simple_fine,10,10);
     FineField noise(FineGrid);
     FineField src(FineGrid);
     FineField guess(FineGrid);
@@ -149,16 +146,16 @@ public:
       
       DiracOp.Op(noise,Mn); std::cout<<GridLogMessage << "noise   ["<<b<<"] <n|Op|n> "<<innerProduct(noise,Mn)<<std::endl;
 
-      for(int i=0;i<3;i++){
+      for(int i=0;i<2;i++){
 	//  void operator() (const Field &src, Field &psi){
 #if 1
-	if (i==0)std::cout << GridLogMessage << " inverting on noise "<<std::endl;
+	std::cout << GridLogMessage << " inverting on noise "<<std::endl;
 	src = noise;
 	guess=Zero();
 	GCR(src,guess);
 	subspace[b] = guess;
 #else
-	if (i==0)std::cout << GridLogMessage << " inverting on zero "<<std::endl;
+	std::cout << GridLogMessage << " inverting on zero "<<std::endl;
 	src=Zero();
 	guess = noise;
 	GCR(src,guess);
@@ -170,7 +167,7 @@ public:
 
       }
 
-      DiracOp.Op(noise,Mn); std::cout<<GridLogMessage << "filtered["<<b<<"] <f|Op|f> "<<innerProduct(noise,Mn)<<" <f|OpDagOp|f>"<<norm2(Mn)<<std::endl;
+      DiracOp.Op(noise,Mn); std::cout<<GridLogMessage << "filtered["<<b<<"] <f|Op|f> "<<innerProduct(noise,Mn)<<std::endl;
       subspace[b]   = noise;
 
     }

Grid/communicator/Communicator_none.cc

@@ -27,8 +27,6 @@ Author: Peter Boyle <paboyle@ph.ed.ac.uk>
 /*  END LEGAL */
 #include <Grid/GridCore.h>
 
-void GridAbort(void) { abort(); }
-
 NAMESPACE_BEGIN(Grid);
 
 ///////////////////////////////////////////////////////////////////////////////////////////////////
@@ -36,6 +34,7 @@ NAMESPACE_BEGIN(Grid);
 ///////////////////////////////////////////////////////////////////////////////////////////////////
 Grid_MPI_Comm       CartesianCommunicator::communicator_world;
 
+void GridAbort(void) { abort(); }
 
 void CartesianCommunicator::Init(int *argc, char *** arv)
 {

Grid/cshift/Cshift_mpi.h

@@ -124,7 +124,7 @@ template<class vobj> void Cshift_simple(Lattice<vobj>& ret,const Lattice<vobj> &
     void *hsend_buf = (void *)&hrhs[0];
     void *hrecv_buf = (void *)&hret[0];
 
-    acceleratorCopyFromDevice(send_buf,hsend_buf,bytes);
+    acceleratorCopyFromDevice(&send_buf[0],&hsend_buf[0],bytes);
 
     grid->SendToRecvFrom(hsend_buf,
 			 xmit_to_rank,
@@ -132,7 +132,7 @@ template<class vobj> void Cshift_simple(Lattice<vobj>& ret,const Lattice<vobj> &
 			 recv_from_rank,
 			 bytes);
 
-    acceleratorCopyToDevice(hrecv_buf,recv_buf,bytes);
+    acceleratorCopyToDevice(&hrecv_buf[0],&recv_buf[0],bytes);
 #endif
   }
 }

Grid/qcd/QCD.h

@@ -596,32 +596,16 @@ template<int Index,class vobj> inline vobj transposeColour(const vobj &lhs){
 //////////////////////////////////////////
 // Trace lattice and non-lattice
 //////////////////////////////////////////
-#define GRID_UNOP(name)   name
-#define GRID_DEF_UNOP(op, name)						\
-  template <typename T1, typename std::enable_if<is_lattice<T1>::value||is_lattice_expr<T1>::value,T1>::type * = nullptr> \
-  inline auto op(const T1 &arg) ->decltype(LatticeUnaryExpression<GRID_UNOP(name),T1>(GRID_UNOP(name)(), arg)) \
-  {									\
-    return     LatticeUnaryExpression<GRID_UNOP(name),T1>(GRID_UNOP(name)(), arg); \
-  }
-
 template<int Index,class vobj>
 inline auto traceSpin(const Lattice<vobj> &lhs) -> Lattice<decltype(traceIndex<SpinIndex>(vobj()))>
 {
   return traceIndex<SpinIndex>(lhs);
 }
-
-GridUnopClass(UnaryTraceSpin, traceIndex<SpinIndex>(a));
-GRID_DEF_UNOP(traceSpin, UnaryTraceSpin);
-
 template<int Index,class vobj>
 inline auto traceColour(const Lattice<vobj> &lhs) -> Lattice<decltype(traceIndex<ColourIndex>(vobj()))>
 {
   return traceIndex<ColourIndex>(lhs);
 }
-
-GridUnopClass(UnaryTraceColour, traceIndex<ColourIndex>(a));
-GRID_DEF_UNOP(traceColour, UnaryTraceColour);
-
 template<int Index,class vobj>
 inline auto traceSpin(const vobj &lhs) -> Lattice<decltype(traceIndex<SpinIndex>(lhs))>
 {
@@ -633,8 +617,6 @@ inline auto traceColour(const vobj &lhs) -> Lattice<decltype(traceIndex<ColourIn
   return traceIndex<ColourIndex>(lhs);
 }
 
-#undef GRID_UNOP
-#undef GRID_DEF_UNOP
 //////////////////////////////////////////
 // Current types
 //////////////////////////////////////////

Grid/qcd/action/fermion/implementation/WilsonKernelsImplementation.h

@@ -411,7 +411,7 @@ void WilsonKernels<Impl>::DhopDirKernel( StencilImpl &st, DoubledGaugeField &U,S
 #undef LoopBody
 }
 
-#if 0
+#ifdef GRID_SYCL
 extern "C" {
     ulong SYCL_EXTERNAL __attribute__((overloadable)) intel_get_cycle_counter( void );
     uint  SYCL_EXTERNAL __attribute__((overloadable)) intel_get_active_channel_mask( void );

Grid/qcd/action/gauge/GaugeImplTypes.h

@@ -138,13 +138,10 @@ public:
     //auto start = std::chrono::high_resolution_clock::now();
     autoView(U_v,U,AcceleratorWrite);
     autoView(P_v,P,AcceleratorRead);
-    typedef typename Field::vector_object vobj;
+    accelerator_for(ss, P.Grid()->oSites(),1,{
-    const int Nsimd = vobj::Nsimd();
-    accelerator_for(ss, P.Grid()->oSites(),Nsimd,{
       for (int mu = 0; mu < Nd; mu++) {
-	auto tmp = Exponentiate(P_v(ss)(mu), ep, Nexp) * U_v(ss)(mu);
+          U_v[ss](mu) = Exponentiate(P_v[ss](mu), ep, Nexp) * U_v[ss](mu);
-	tmp      = Group::ProjectOnGeneralGroup(tmp);
+          U_v[ss](mu) = Group::ProjectOnGeneralGroup(U_v[ss](mu));
-	coalescedWrite(U_v[ss](mu),tmp);
       }
     });
    //auto end = std::chrono::high_resolution_clock::now();

Grid/qcd/modules/ObservableModules.h

@@ -103,18 +103,6 @@ class PolyakovMod: public ObservableModule<PolyakovLogger<Impl>, NoParameters>{
   PolyakovMod(): ObsBase(NoParameters()){}
 };
 
-template < class Impl >
-class SpatialPolyakovMod: public ObservableModule<SpatialPolyakovLogger<Impl>, NoParameters>{
-  typedef ObservableModule<SpatialPolyakovLogger<Impl>, NoParameters> ObsBase;
-  using ObsBase::ObsBase; // for constructors
-
-  // acquire resource
-  virtual void initialize(){
-    this->ObservablePtr.reset(new SpatialPolyakovLogger<Impl>());
-  }
-  public:
-  SpatialPolyakovMod(): ObsBase(NoParameters()){}
-};
 
 template < class Impl >
 class TopologicalChargeMod: public ObservableModule<TopologicalCharge<Impl>, TopologyObsParameters>{

Grid/qcd/observables/polyakov_loop.h

@@ -2,12 +2,11 @@
 
 Grid physics library, www.github.com/paboyle/Grid
 
-Source file: ./Grid/qcd/observables/polyakov_loop.h
+Source file: ./lib/qcd/modules/polyakov_line.h
 
-Copyright (C) 2025
+Copyright (C) 2017
 
 Author: David Preti <david.preti@csic.es>
-Author: Alexis Verney-Provatas <2414441@swansea.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
@@ -61,43 +60,4 @@ class PolyakovLogger : public HmcObservable<typename Impl::Field> {
   }
 };
 
-template <class Impl>
-class SpatialPolyakovLogger : public HmcObservable<typename Impl::Field> {
- public:
-    // here forces the Impl to be of gauge fields
-    // if not the compiler will complain
-    INHERIT_GIMPL_TYPES(Impl);
-
-     // necessary for HmcObservable compatibility
-    typedef typename Impl::Field Field;
-
-    void TrajectoryComplete(int traj,
-                            Field &U,
-                            GridSerialRNG &sRNG,
-                            GridParallelRNG &pRNG) {
-
-    // Save current numerical output precision
-    int def_prec = std::cout.precision();
-
-    // Assume that the dimensions are D=3+1
-    int Ndim = 3;
-    ComplexD polyakov;
-   
-    // Iterate over the spatial directions and print the average spatial polyakov loop
-    // over them 
-    for (int idx=0; idx<Ndim; idx++) {
-        polyakov = WilsonLoops<Impl>::avgPolyakovLoop(U, idx);
-    
-        std::cout << GridLogMessage
-            << std::setprecision(std::numeric_limits<Real>::digits10 + 1)
-            << "Polyakov Loop in the " << idx << " spatial direction : [ " << traj << " ] "<< polyakov << std::endl;
-
-    }
-
-    // Return to original output precision
-    std::cout.precision(def_prec);
-
-  }
-};
-
 NAMESPACE_END(Grid);

Grid/qcd/smearing/GaugeConfigurationMasked.h

@@ -291,8 +291,8 @@ public:
     int idx=0;
     for(int mu=0;mu<4;mu++){
     for(int nu=0;nu<4;nu++){
-      if ( mu!=nu) assert(this->StoutSmearing->SmearRho[idx]==rho);
+      if ( mu!=nu) GRID_ASSERT(this->StoutSmearing->SmearRho[idx]==rho);
-      else         assert(this->StoutSmearing->SmearRho[idx]==0.0);
+      else         GRID_ASSERT(this->StoutSmearing->SmearRho[idx]==0.0);
       idx++;
     }}
     //////////////////////////////////////////////////////////////////
@@ -825,7 +825,6 @@ public:
   virtual void fill_smearedSet(GaugeField &U)
   {
     this->ThinLinks = &U;  // attach the smearing routine to the field U
-    std::cout << GridLogMessage << " fill_smearedSet " << WilsonLoops<PeriodicGimplR>::avgPlaquette(U) << std::endl;
 
     // check the pointer is not null
     if (this->ThinLinks == NULL)
@@ -847,8 +846,6 @@ public:
 	ApplyMask(smeared_A,smearLvl);
 	smeared_B = previous_u;
 	ApplyMask(smeared_B,smearLvl);
-	std::cout << GridLogMessage << " smeared_A " << norm2(smeared_A) << std::endl;
-	std::cout << GridLogMessage << " smeared_B " << norm2(smeared_B) << std::endl;
 	// Replace only the masked portion
 	this->SmearedSet[smearLvl] = previous_u-smeared_B + smeared_A;
         previous_u = this->SmearedSet[smearLvl];
@@ -937,10 +934,10 @@ public:
   SmearedConfigurationMasked(GridCartesian* _UGrid, unsigned int Nsmear, Smear_Stout<Gimpl>& Stout)
     : SmearedConfiguration<Gimpl>(_UGrid, Nsmear,Stout)
   {
-    assert(Nsmear%(2*Nd)==0); // Or multiply by 8??
+    GRID_ASSERT(Nsmear%(2*Nd)==0); // Or multiply by 8??
 
     // was resized in base class
-    assert(this->SmearedSet.size()==Nsmear);
+    GRID_ASSERT(this->SmearedSet.size()==Nsmear);
     
     GridRedBlackCartesian * UrbGrid;
     UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(_UGrid);

Grid/qcd/smearing/JacobianAction.h

@@ -54,7 +54,7 @@ public:
   // Usual cases are not used
   //////////////////////////////////
   virtual void refresh(const GaugeField &U, GridSerialRNG &sRNG, GridParallelRNG &pRNG){ GRID_ASSERT(0);};
-  virtual RealD S(const GaugeField &U) { GRID_ASSERT(0); return 0; }
+  virtual RealD S(const GaugeField &U) { GRID_ASSERT(0); }
   virtual void deriv(const GaugeField &U, GaugeField &dSdU) { GRID_ASSERT(0);  }
 
   //////////////////////////////////

Grid/qcd/utils/Sp2n.impl.h

@@ -254,9 +254,9 @@ static void testGenerators(GroupName::Sp) {
   }
 }
 
-template <class vtype, int N>
+template <int N>
-static Lattice<iScalar<iScalar<iMatrix<vtype, N> > > >
+static Lattice<iScalar<iScalar<iMatrix<vComplexD, N> > > >
-ProjectOnGeneralGroup(const Lattice<iScalar<iScalar<iMatrix<vtype, N> > > > &Umu, GroupName::Sp) {
+ProjectOnGeneralGroup(const Lattice<iScalar<iScalar<iMatrix<vComplexD, N> > > > &Umu, GroupName::Sp) {
   return ProjectOnSpGroup(Umu);
 }
 

Grid/qcd/utils/WilsonLoops.h

@@ -177,43 +177,25 @@ public:
   }
 
   //////////////////////////////////////////////////
-  // average Polyakov loop in mu direction over all directions != mu
+  // average over all x,y,z the temporal loop
   //////////////////////////////////////////////////
-  static ComplexD avgPolyakovLoop(const GaugeField &Umu, const int mu) {  //assume Nd=4
+  static ComplexD avgPolyakovLoop(const GaugeField &Umu) {  //assume Nd=4
-    
+    GaugeMat Ut(Umu.Grid()), P(Umu.Grid());
-    // Protect against bad value of mu [0, 3]
-    if ((mu < 0 ) || (mu > 3)) {
-      std::cout << GridLogError << "Index is not an integer inclusively between 0 and 3." << std::endl;
-      exit(1);
-    }
-
-    // U_loop is U_{mu}
-    GaugeMat U_loop(Umu.Grid()), P(Umu.Grid());
     ComplexD out;
     int T = Umu.Grid()->GlobalDimensions()[3];
     int X = Umu.Grid()->GlobalDimensions()[0];
     int Y = Umu.Grid()->GlobalDimensions()[1];
     int Z = Umu.Grid()->GlobalDimensions()[2];
 
-    // Number of sites in mu direction
+    Ut = peekLorentz(Umu,3); //Select temporal direction
-    int N_mu = Umu.Grid()->GlobalDimensions()[mu];
+    P = Ut;
-
+    for (int t=1;t<T;t++){ 
-    U_loop = peekLorentz(Umu, mu); //Select direction
+      P = Gimpl::CovShiftForward(Ut,3,P);
-    P = U_loop;
-    for (int t=1;t<N_mu;t++){ 
-      P = Gimpl::CovShiftForward(U_loop,mu,P);
     }
    RealD norm = 1.0/(Nc*X*Y*Z*T);
    out = sum(trace(P))*norm;
    return out;   
-  }  
+}
-
-  /////////////////////////////////////////////////
-  // overload for temporal Polyakov loop
-  /////////////////////////////////////////////////
-  static ComplexD avgPolyakovLoop(const GaugeField &Umu) { 
-    return avgPolyakovLoop(Umu, 3);
-  }
 
   //////////////////////////////////////////////////
   // average over traced single links

Grid/stencil/Stencil.h

@@ -751,7 +751,7 @@ public:
     obj.xbytes      = xbytes;
     obj.rbytes      = rbytes;
     obj.cb          = cb;
-    
+
     for(int i=0;i<CachedTransfers.size();i++){
       if (   (CachedTransfers[i].direction  ==direction)
 	   &&(CachedTransfers[i].OrthogPlane==OrthogPlane)
@@ -763,13 +763,11 @@ public:
 	     ){
 	// FIXME worry about duplicate with partial compression
 	// Wont happen as DWF has no duplicates, but...
-	//	AddCopy(CachedTransfers[i].recv_buf,recv_buf,rbytes);
+	AddCopy(CachedTransfers[i].recv_buf,recv_buf,rbytes);
-	//	std::cout << "Duplicate dir " <<direction<<" "<<" OrthogPlane "<<OrthogPlane<<" Dest"<<DestProc <<" xbytes " <<xbytes<<" lane "<< lane<<" cb "<<cb<<std::endl;
+	return 1;
-	return 0;
-	
-	//	return 1;
       }
     }
+
     CachedTransfers.push_back(obj);
     return 0;
   }

Grid/util/Init.cc

@@ -755,7 +755,7 @@ void Grid_generic_handler(int sig,siginfo_t *si,void * ptr)
   sig_print_uint(si->si_code);
   SIGLOG("\n");
 
-  ucontext_t *uc= (ucontext_t *)ptr;
+  unw_context_t *uc= (unw_context_t *)ptr;
   
   SIGLOG("Backtrace:\n");
 #ifdef HAVE_UNWIND

HMC/FTHMC2p1f.cc

@@ -24,11 +24,7 @@ See the full license in the file "LICENSE" in the top level distribution
 directory
 *************************************************************************************/
 /*  END LEGAL */
-
+#include <Grid/Grid.h>
-#include "disable_examples_without_instantiations.h"
-#ifdef ENABLE_FERMION_INSTANTIATIONS
-
-#include<Grid/Grid.h>
 
 #if Nc == 3
 #include <Grid/qcd/smearing/GaugeConfigurationMasked.h>
@@ -234,4 +230,3 @@ int main(int argc, char **argv)
 #endif
 } // main
 
-#endif

HMC/FTHMC2p1f_3GeV.cc

@@ -25,11 +25,7 @@ directory
 *************************************************************************************/
 /*  END LEGAL */
 
-
+#include <Grid/Grid.h>
-#include "disable_examples_without_instantiations.h"
-#ifdef ENABLE_FERMION_INSTANTIATIONS
-
-#include<Grid/Grid.h>
 
 #if Nc == 3
 #include <Grid/qcd/smearing/GaugeConfigurationMasked.h>
@@ -235,4 +231,5 @@ int main(int argc, char **argv)
 #endif
 } // main
 
-#endif
+
+

HMC/HMC2p1f_3GeV.cc

@@ -24,11 +24,7 @@ See the full license in the file "LICENSE" in the top level distribution
 directory
 *************************************************************************************/
 /*  END LEGAL */
-
+#include <Grid/Grid.h>
-#include "disable_examples_without_instantiations.h"
-#ifdef ENABLE_FERMION_INSTANTIATIONS
-
-#include<Grid/Grid.h>
 
 #if Nc == 3
 #include <Grid/qcd/smearing/GaugeConfigurationMasked.h>
@@ -234,4 +230,5 @@ int main(int argc, char **argv)
 #endif
 } // main
 
-#endif
+
+

HMC/Mobius2p1f.cc

@@ -27,11 +27,7 @@ See the full license in the file "LICENSE" in the top level distribution
 directory
 *************************************************************************************/
 /*  END LEGAL */
-
+#include <Grid/Grid.h>
-#include "disable_examples_without_instantiations.h"
-#ifdef ENABLE_FERMION_INSTANTIATIONS
-
-#include<Grid/Grid.h>
 
 int main(int argc, char **argv) {
   using namespace Grid;
@@ -199,4 +195,5 @@ int main(int argc, char **argv) {
   Grid_finalize();
 } // main
 
-#endif
+
+

HMC/Mobius2p1fEOFA.cc

@@ -28,11 +28,7 @@ See the full license in the file "LICENSE" in the top level distribution
 directory
 *************************************************************************************/
 /*  END LEGAL */
-
+#include <Grid/Grid.h>
-#include "disable_examples_without_instantiations.h"
-#ifdef ENABLE_FERMION_INSTANTIATIONS
-
-#include<Grid/Grid.h>
 
 #ifdef GRID_DEFAULT_PRECISION_DOUBLE
 #define MIXED_PRECISION
@@ -453,4 +449,5 @@ int main(int argc, char **argv) {
   Grid_finalize();
 } // main
 
-#endif
+
+

HMC/Mobius2p1fEOFA_F1.cc

@@ -28,11 +28,7 @@ See the full license in the file "LICENSE" in the top level distribution
 directory
 *************************************************************************************/
 /*  END LEGAL */
-
+#include <Grid/Grid.h>
-#include "disable_examples_without_instantiations.h"
-#ifdef ENABLE_FERMION_INSTANTIATIONS
-
-#include<Grid/Grid.h>
 
 #ifdef GRID_DEFAULT_PRECISION_DOUBLE
 #define MIXED_PRECISION
@@ -446,4 +442,5 @@ int main(int argc, char **argv) {
   Grid_finalize();
 } // main
 
-#endif
+
+

HMC/Mobius2p1fIDSDRGparityEOFA_40ID.cc

@@ -28,11 +28,7 @@ See the full license in the file "LICENSE" in the top level distribution
 directory
 *************************************************************************************/
 /*  END LEGAL */
-
+#include <Grid/Grid.h>
-#include "disable_examples_without_instantiations.h"
-#ifdef ENABLE_FERMION_INSTANTIATIONS
-
-#include<Grid/Grid.h>
 
 using namespace Grid;
 
@@ -922,5 +918,3 @@ int main(int argc, char **argv) {
   return 0;
 #endif
 } // main
-
-#endif

HMC/Mobius2p1fIDSDRGparityEOFA_48ID.cc

@@ -28,11 +28,7 @@ See the full license in the file "LICENSE" in the top level distribution
 directory
 *************************************************************************************/
 /*  END LEGAL */
-
+#include <Grid/Grid.h>
-#include "disable_examples_without_instantiations.h"
-#ifdef ENABLE_FERMION_INSTANTIATIONS
-
-#include<Grid/Grid.h>
 
 using namespace Grid;
 
@@ -877,5 +873,3 @@ int main(int argc, char **argv) {
   return 0;
 #endif
 } // main
-
-#endif

HMC/Mobius2p1fRHMC.cc

@@ -27,11 +27,7 @@ See the full license in the file "LICENSE" in the top level distribution
 directory
 *************************************************************************************/
 /*  END LEGAL */
-
+#include <Grid/Grid.h>
-#include "disable_examples_without_instantiations.h"
-#ifdef ENABLE_FERMION_INSTANTIATIONS
-
-#include<Grid/Grid.h>
 
 int main(int argc, char **argv) {
   using namespace Grid;
@@ -197,4 +193,5 @@ int main(int argc, char **argv) {
   Grid_finalize();
 } // main
 
-#endif
+
+

HMC/Mobius2p1f_DD_EOFA_96I_3level.cc

@@ -27,11 +27,7 @@ See the full license in the file "LICENSE" in the top level distribution
 directory
 *************************************************************************************/
 /*  END LEGAL */
-
+#include <Grid/Grid.h>
-#include "disable_examples_without_instantiations.h"
-#ifdef ENABLE_FERMION_INSTANTIATIONS
-
-#include<Grid/Grid.h>
 
 NAMESPACE_BEGIN(Grid);
 
@@ -516,4 +512,5 @@ int main(int argc, char **argv) {
   Grid_finalize();
 } // main
 
-#endif
+
+

HMC/Mobius2p1f_DD_EOFA_96I_double.cc

@@ -27,11 +27,7 @@ See the full license in the file "LICENSE" in the top level distribution
 directory
 *************************************************************************************/
 /*  END LEGAL */
-
+#include <Grid/Grid.h>
-#include "disable_examples_without_instantiations.h"
-#ifdef ENABLE_FERMION_INSTANTIATIONS
-
-#include<Grid/Grid.h>
 
 int main(int argc, char **argv) {
   using namespace Grid;
@@ -349,4 +345,5 @@ int main(int argc, char **argv) {
   Grid_finalize();
 } // main
 
-#endif
+
+

HMC/Mobius2p1f_DD_EOFA_96I_mixed.cc

@@ -27,11 +27,7 @@ See the full license in the file "LICENSE" in the top level distribution
 directory
 *************************************************************************************/
 /*  END LEGAL */
-
+#include <Grid/Grid.h>
-#include "disable_examples_without_instantiations.h"
-#ifdef ENABLE_FERMION_INSTANTIATIONS
-
-#include<Grid/Grid.h>
 
 NAMESPACE_BEGIN(Grid);
 
@@ -520,4 +516,5 @@ int main(int argc, char **argv) {
   Grid_finalize();
 } // main
 
-#endif
+
+

HMC/Mobius2p1f_DD_EOFA_96I_mshift.cc

@@ -27,11 +27,7 @@ See the full license in the file "LICENSE" in the top level distribution
 directory
 *************************************************************************************/
 /*  END LEGAL */
-
+#include <Grid/Grid.h>
-#include "disable_examples_without_instantiations.h"
-#ifdef ENABLE_FERMION_INSTANTIATIONS
-
-#include<Grid/Grid.h>
 
 NAMESPACE_BEGIN(Grid);
 
@@ -571,4 +567,5 @@ int main(int argc, char **argv) {
   Grid_finalize();
 } // main
 
-#endif
+
+

HMC/Mobius2p1f_DD_RHMC.cc

@@ -27,11 +27,7 @@ See the full license in the file "LICENSE" in the top level distribution
 directory
 *************************************************************************************/
 /*  END LEGAL */
-
+#include <Grid/Grid.h>
-#include "disable_examples_without_instantiations.h"
-#ifdef ENABLE_FERMION_INSTANTIATIONS
-
-#include<Grid/Grid.h>
 
 int main(int argc, char **argv) {
   using namespace Grid;
@@ -267,4 +263,5 @@ int main(int argc, char **argv) {
   Grid_finalize();
 } // main
 
-#endif
+
+

HMC/Mobius2p1f_DD_RHMC_96I.cc

@@ -27,11 +27,7 @@ See the full license in the file "LICENSE" in the top level distribution
 directory
 *************************************************************************************/
 /*  END LEGAL */
-
+#include <Grid/Grid.h>
-#include "disable_examples_without_instantiations.h"
-#ifdef ENABLE_FERMION_INSTANTIATIONS
-
-#include<Grid/Grid.h>
 
 int main(int argc, char **argv) {
   using namespace Grid;
@@ -421,4 +417,5 @@ int main(int argc, char **argv) {
   Grid_finalize();
 } // main
 
-#endif
+
+

HMC/Mobius2p1f_DD_RHMC_96I_mixed.cc

@@ -27,11 +27,7 @@ See the full license in the file "LICENSE" in the top level distribution
 directory
 *************************************************************************************/
 /*  END LEGAL */
-
+#include <Grid/Grid.h>
-#include "disable_examples_without_instantiations.h"
-#ifdef ENABLE_FERMION_INSTANTIATIONS
-
-#include<Grid/Grid.h>
 
 NAMESPACE_BEGIN(Grid);
 
@@ -456,4 +452,5 @@ int main(int argc, char **argv) {
   Grid_finalize();
 } // main
 
-#endif
+
+

HMC/Mobius2p1f_EOFA_96I_hmc.cc

@@ -27,11 +27,7 @@ See the full license in the file "LICENSE" in the top level distribution
 directory
 *************************************************************************************/
 /*  END LEGAL */
-
+#include <Grid/Grid.h>
-#include "disable_examples_without_instantiations.h"
-#ifdef ENABLE_FERMION_INSTANTIATIONS
-
-#include<Grid/Grid.h>
 
 NAMESPACE_BEGIN(Grid);
 
@@ -466,4 +462,5 @@ int main(int argc, char **argv) {
   Grid_finalize();
 } // main
 
-#endif
+
+

HMC/Mobius2p1f_EOFA_96I_hmc_double.cc

@@ -27,11 +27,7 @@ See the full license in the file "LICENSE" in the top level distribution
 directory
 *************************************************************************************/
 /*  END LEGAL */
-
+#include <Grid/Grid.h>
-#include "disable_examples_without_instantiations.h"
-#ifdef ENABLE_FERMION_INSTANTIATIONS
-
-#include<Grid/Grid.h>
 
 
 
@@ -268,4 +264,5 @@ int main(int argc, char **argv) {
   Grid_finalize();
 } // main
 
-#endif
+
+

HMC/disable_examples_without_instantiations.h

@@ -1,16 +0,0 @@
-#include <Grid/Grid.h>
-#pragma once
-
-
-#ifndef ENABLE_FERMION_INSTANTIATIONS
-#include <iostream>
-
-int main(void) {
-  std::cout << "This build of Grid was configured to exclude fermion instantiations, "
-	    << "which this example relies on. "
-	    << "Please reconfigure and rebuild Grid with --enable-fermion-instantiations"
-	    << "to run this example."
-	    << std::endl;
-  return 1;
-}
-#endif

benchmarks/Benchmark_ITT.cc

@@ -26,9 +26,6 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
     See the full license in the file "LICENSE" in the top level distribution directory
     *************************************************************************************/
     /*  END LEGAL */
-#include "disable_benchmarks_without_instantiations.h"
-#ifdef ENABLE_FERMION_INSTANTIATIONS
-
 #include <Grid/Grid.h>
 
 using namespace Grid;
@@ -734,5 +731,3 @@ int main (int argc, char ** argv)
 
   Grid_finalize();
 }
-
-#endif

benchmarks/Benchmark_dwf.cc

@@ -20,9 +20,6 @@
     See the full license in the file "LICENSE" in the top level distribution directory
     *************************************************************************************/
     /*  END LEGAL */
-#include "disable_benchmarks_without_instantiations.h"
-#ifdef ENABLE_FERMION_INSTANTIATIONS
-
 #include <Grid/Grid.h>
 #ifdef GRID_CUDA
 #define CUDA_PROFILE
@@ -442,4 +439,3 @@ void Benchmark(int Ls, Coordinate Dirichlet,bool sloppy)
   GRID_ASSERT(norm2(src_e)<1.0e-4);
   GRID_ASSERT(norm2(src_o)<1.0e-4);
 }
-#endif

benchmarks/Benchmark_dwf_fp32.cc

@@ -20,10 +20,6 @@
     See the full license in the file "LICENSE" in the top level distribution directory
     *************************************************************************************/
     /*  END LEGAL */
-#include "disable_benchmarks_without_instantiations.h"
-
-#ifdef ENABLE_FERMION_INSTANTIATIONS
-
 #include <Grid/Grid.h>
 #ifdef GRID_CUDA
 #define CUDA_PROFILE
@@ -443,5 +439,3 @@ void Benchmark(int Ls, Coordinate Dirichlet,bool sloppy)
   GRID_ASSERT(norm2(src_e)<1.0e-4);
   GRID_ASSERT(norm2(src_o)<1.0e-4);
 }
-
-#endif

benchmarks/Benchmark_dwf_fp32_paranoid.cc

@@ -20,9 +20,6 @@
     See the full license in the file "LICENSE" in the top level distribution directory
     *************************************************************************************/
     /*  END LEGAL */
-#include "disable_benchmarks_without_instantiations.h"
-#ifdef ENABLE_FERMION_INSTANTIATIONS
-
 #include <Grid/Grid.h>
 #ifdef GRID_CUDA
 #define CUDA_PROFILE
@@ -388,5 +385,3 @@ int main (int argc, char ** argv)
   Grid_finalize();
   exit(0);
 }
-
-#endif

benchmarks/Benchmark_dwf_sweep.cc

@@ -26,9 +26,6 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
     See the full license in the file "LICENSE" in the top level distribution directory
     *************************************************************************************/
     /*  END LEGAL */
-#include "disable_benchmarks_without_instantiations.h"
-#ifdef ENABLE_FERMION_INSTANTIATIONS
-
 #include <Grid/Grid.h>
 
 using namespace std;
@@ -241,4 +238,5 @@ void benchDw(std::vector<int> & latt4, int Ls, int threads,int report )
   }
 }
 
-#endif
+
+

benchmarks/Benchmark_gparity.cc

@@ -1,7 +1,3 @@
-#include "disable_benchmarks_without_instantiations.h"
-#ifdef ENABLE_FERMION_INSTANTIATIONS
-
-
 #include <Grid/Grid.h>
 #include <sstream>
 using namespace std;
@@ -159,4 +155,3 @@ int main (int argc, char ** argv)
   Grid_finalize();
 }
 
-#endif

benchmarks/Benchmark_halo.cc

@@ -20,9 +20,6 @@
     See the full license in the file "LICENSE" in the top level distribution directory
     *************************************************************************************/
     /*  END LEGAL */
-#include "disable_benchmarks_without_instantiations.h"
-#ifdef ENABLE_FERMION_INSTANTIATIONS
-
 #include <Grid/Grid.h>
 #ifdef GRID_CUDA
 #define CUDA_PROFILE
@@ -132,5 +129,3 @@ int main (int argc, char ** argv)
   Grid_finalize();
   exit(0);
 }
-
-#endif

benchmarks/Benchmark_mooee.cc

@@ -26,9 +26,6 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
     See the full license in the file "LICENSE" in the top level distribution directory
     *************************************************************************************/
     /*  END LEGAL */
-#include "disable_benchmarks_without_instantiations.h"
-#ifdef ENABLE_FERMION_INSTANTIATIONS
-
 #include <Grid/Grid.h>
 
 using namespace std;
@@ -152,5 +149,3 @@ int main (int argc, char ** argv)
 
   Grid_finalize();
 }
-
-#endif

benchmarks/Benchmark_schur.cc

@@ -26,9 +26,6 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
     See the full license in the file "LICENSE" in the top level distribution directory
     *************************************************************************************/
     /*  END LEGAL */
-#include "disable_benchmarks_without_instantiations.h"
-#ifdef ENABLE_FERMION_INSTANTIATIONS
-
 #include <Grid/Grid.h>
 
 using namespace std;
@@ -175,4 +172,5 @@ void benchDw(std::vector<int> & latt4, int Ls)
   //  Dw.Report();
 }
 
-#endif
+
+

benchmarks/Benchmark_staggered.cc

@@ -26,9 +26,6 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
     See the full license in the file "LICENSE" in the top level distribution directory
     *************************************************************************************/
     /*  END LEGAL */
-#include "disable_benchmarks_without_instantiations.h"
-#ifdef ENABLE_FERMION_INSTANTIATIONS
-
 #include <Grid/Grid.h>
 
 using namespace std;
@@ -113,5 +110,3 @@ int main (int argc, char ** argv)
 
   Grid_finalize();
 }
-
-#endif

benchmarks/Benchmark_staggeredF.cc

@@ -26,9 +26,6 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
     See the full license in the file "LICENSE" in the top level distribution directory
     *************************************************************************************/
     /*  END LEGAL */
-#include "disable_benchmarks_without_instantiations.h"
-#ifdef ENABLE_FERMION_INSTANTIATIONS
-
 #include <Grid/Grid.h>
 
 using namespace std;
@@ -115,5 +112,3 @@ int main (int argc, char ** argv)
 
   Grid_finalize();
 }
-
-#endif

benchmarks/Benchmark_usqcd.cc

@@ -26,10 +26,6 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
     See the full license in the file "LICENSE" in the top level distribution directory
     *************************************************************************************/
     /*  END LEGAL */
-#include "disable_benchmarks_without_instantiations.h"
-#ifdef ENABLE_FERMION_INSTANTIATIONS
-
-
 #include <Grid/Grid.h>
 #include <Grid/algorithms/blas/BatchedBlas.h>
 
@@ -982,5 +978,3 @@ int main (int argc, char ** argv)
   Grid_finalize();
   fclose(FP);
 }
-
-#endif

benchmarks/Benchmark_wilson.cc

@@ -26,9 +26,6 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
     See the full license in the file "LICENSE" in the top level distribution directory
     *************************************************************************************/
     /*  END LEGAL */
-#include "disable_benchmarks_without_instantiations.h"
-#ifdef ENABLE_FERMION_INSTANTIATIONS
-
 #include <Grid/Grid.h>
 
 using namespace std;
@@ -261,5 +258,3 @@ int main (int argc, char ** argv)
 
   Grid_finalize();
 }
-
-#endif

benchmarks/Benchmark_wilson_sweep.cc

@@ -19,9 +19,6 @@ Author: Richard Rollins <rprollins@users.noreply.github.com>
     See the full license in the file "LICENSE" in the top level distribution directory
 *************************************************************************************/
 /*  END LEGAL */
-#include "disable_benchmarks_without_instantiations.h"
-#ifdef ENABLE_FERMION_INSTANTIATIONS
-
 #include <Grid/Grid.h>
 
 using namespace std;
@@ -164,5 +161,3 @@ void bench_wilson_eo (
   double flops = (single_site_flops * volume * ncall)/2.0;
   std::cout << flops/(t1-t0) << "\t\t";
 }
-
-#endif

benchmarks/disable_benchmarks_without_instantiations.h

@@ -1,16 +0,0 @@
-#include <Grid/Grid.h>
-
-#pragma once
-
-#ifndef ENABLE_FERMION_INSTANTIATIONS
-#include <iostream>
-
-int main(void) {
-  std::cout << "This build of Grid was configured to exclude fermion instantiations, "
-	    << "which this benchmark relies on. "
-	    << "Please reconfigure and rebuild Grid with --enable-fermion-instantiations"
-	    << "to run this benchmark."
-	    << std::endl;
-  return 1;
-}
-#endif

configure.ac

@@ -172,12 +172,6 @@ case ${ac_TRACING} in
 esac
 
 ############### fermions
-AC_ARG_ENABLE([fermion-instantiations],
-     [AS_HELP_STRING([--enable-fermion-instantiations=yes|no],[enable fermion instantiations])],
-     [ac_FERMION_REPS=${enable_fermion_instantiations}], [ac_FERMION_INSTANTIATIONS=yes])
-
-AM_CONDITIONAL(BUILD_FERMION_INSTANTIATIONS, [ test "${ac_FERMION_INSTANTIATIONS}X" == "yesX" ])
-
 AC_ARG_ENABLE([fermion-reps],
      [AS_HELP_STRING([--enable-fermion-reps=yes|no],[enable extra fermion representation support])],
      [ac_FERMION_REPS=${enable_fermion_reps}], [ac_FERMION_REPS=yes])
@@ -200,9 +194,6 @@ AM_CONDITIONAL(BUILD_ZMOBIUS, [ test "${ac_ZMOBIUS}X" == "yesX" ])
 case ${ac_FERMION_REPS} in
    yes) AC_DEFINE([ENABLE_FERMION_REPS],[1],[non QCD fermion reps]);;
 esac
-case ${ac_FERMION_INSTANTIATIONS} in
-   yes) AC_DEFINE([ENABLE_FERMION_INSTANTIATIONS],[1],[enable fermions]);;
-esac
 case ${ac_GPARITY} in
    yes) AC_DEFINE([ENABLE_GPARITY],[1],[fermion actions with GPARITY BCs]);;
 esac
@@ -418,10 +409,16 @@ AC_SEARCH_LIBS([unw_backtrace], [unwind],
                [have_unwind=true],
 	             [AC_MSG_WARN(libunwind library was not found in your system.)])
 
-AC_SEARCH_LIBS([_Ux86_64_step], [unwind-x86_64],
+AS_CASE([$host_cpu], [x86_64],
-               [AC_DEFINE([HAVE_UNWIND_X86_64], [1], [Define to 1 if you have the `libunwind-x86_64' library])]
+  [AC_SEARCH_LIBS([_Ux86_64_step], [unwind-x86_64],
-               [have_unwind_x86_64=true],
+                [AC_DEFINE([HAVE_UNWIND_X86_64], [1], [Define to 1 if you have the `libunwind-x86_64' library])]
-	             [AC_MSG_WARN(libunwind library was not found in your system.)])
+                [have_unwind_x86_64=true],
+                [AC_MSG_WARN(libunwind library was not found in your system.)])],
+  [aarch64],
+  [AC_SEARCH_LIBS([_Uaarch64_step], [unwind-aarch64],
+                [AC_DEFINE([HAVE_UNWIND_AARCH64], [1], [Define to 1 if you have the `libunwind-aarch64' library])]
+                [have_unwind_aarch64=true],
+                [AC_MSG_WARN(libunwind library was not found in your system.)])])
 
 AC_SEARCH_LIBS([SHA256_Init], [crypto],
                [AC_DEFINE([HAVE_CRYPTO], [1], [Define to 1 if you have the `OpenSSL' library])]

examples/Example_Mobius_spectrum.cc

@@ -3,9 +3,6 @@
  * without regression / tests being applied
  */
 
-#include "disable_examples_without_instantiations.h"
-#ifdef ENABLE_FERMION_INSTANTIATIONS
-
 #include <Grid/Grid.h>
 
 using namespace std;
@@ -313,4 +310,5 @@ int main (int argc, char ** argv)
   Grid_finalize();
 }
 
-#endif
+
+

examples/Example_christoph.cc

@@ -3,9 +3,6 @@
  * without regression / tests being applied
  */
 
-#include "disable_examples_without_instantiations.h"
-#ifdef ENABLE_FERMION_INSTANTIATIONS
-
 #include <Grid/Grid.h>
 
 using namespace std;
@@ -435,4 +432,5 @@ int main (int argc, char ** argv)
   Grid_finalize();
 }
 
-#endif
+
+

examples/Example_wall_wall_3pt.cc

@@ -3,9 +3,6 @@
  * without regression / tests being applied
  */
 
-#include "disable_examples_without_instantiations.h"
-#ifdef ENABLE_FERMION_INSTANTIATIONS
-
 #include <Grid/Grid.h>
 
 using namespace std;
@@ -538,4 +535,5 @@ int main (int argc, char ** argv)
   Grid_finalize();
 }
 
-#endif
+
+

examples/Example_wall_wall_spectrum.cc

@@ -3,9 +3,6 @@
  * without regression / tests being applied
  */
 
-#include "disable_examples_without_instantiations.h"
-#ifdef ENABLE_FERMION_INSTANTIATIONS
-
 #include <Grid/Grid.h>
 
 using namespace std;
@@ -432,4 +429,5 @@ int main (int argc, char ** argv)
   Grid_finalize();
 }
 
-#endif
+
+

examples/disable_examples_without_instantiations.h

@@ -1,15 +0,0 @@
-#include <Grid/Grid.h>
-#pragma once
-
-#ifndef ENABLE_FERMION_INSTANTIATIONS
-#include <iostream>
-
-int main(void) {
-  std::cout << "This build of Grid was configured to exclude fermion instantiations, "
-	    << "which this example relies on. "
-	    << "Please reconfigure and rebuild Grid with --enable-fermion-instantiations"
-	    << "to run this example."
-	    << std::endl;
-  return 1;
-}
-#endif

systems/Aurora/config-command

@@ -10,7 +10,7 @@ export HDF5=/opt/cray/pe/hdf5/1.12.2.3/gnu/9.1
 	--disable-gparity \
 	--disable-fermion-reps \
 	--enable-shm=nvlink \
-	--enable-checksum-comms=no \
+	--enable-checksum-comms=yes \
 	--enable-log-views=yes \
 	--enable-accelerator=sycl \
 	--enable-accelerator-aware-mpi=no \

tests/Test_cayley_even_odd_vec.cc

@@ -25,9 +25,6 @@ Author: Peter Boyle <paboyle@ph.ed.ac.uk>
     See the full license in the file "LICENSE" in the top level distribution directory
     *************************************************************************************/
     /*  END LEGAL */
-#include "disable_tests_without_instantiations.h"
-#ifdef ENABLE_FERMION_INSTANTIATIONS
-
 #include <Grid/Grid.h>
 
 using namespace std;
@@ -276,6 +273,8 @@ void  TestWhat(What & Ddwf,
 
   err = phi-chi;
   std::cout<<GridLogMessage << "norm diff   "<< norm2(err)<< std::endl;
+
+  
 }
 
-#endif
+

tests/Test_compressed_lanczos_hot_start.cc

@@ -30,9 +30,6 @@ Author: Peter Boyle <paboyle@ph.ed.ac.uk>
  *  Reimplement the badly named "multigrid" lanczos as compressed Lanczos using the features 
  *  in Grid that were intended to be used to support blocked Aggregates, from
  */
-#include "disable_tests_without_instantiations.h"
-#ifdef ENABLE_FERMION_INSTANTIATIONS
-
 #include <Grid/Grid.h>
 #include <Grid/algorithms/iterative/ImplicitlyRestartedLanczos.h>
 #include <Grid/algorithms/iterative/LocalCoherenceLanczos.h>
@@ -259,4 +256,3 @@ int main (int argc, char ** argv) {
   Grid_finalize();
 }
 
-#endif

tests/Test_dwf_dslash_repro.cc

@@ -25,9 +25,6 @@ Author: Peter Boyle <paboyle@ph.ed.ac.uk>
     See the full license in the file "LICENSE" in the top level distribution directory
     *************************************************************************************/
     /*  END LEGAL */
-#include "disable_tests_without_instantiations.h"
-#ifdef ENABLE_FERMION_INSTANTIATIONS
-
 #include <Grid/Grid.h>
 
 using namespace std;
@@ -240,5 +237,3 @@ int main (int argc, char ** argv)
   
   Grid_finalize();
 }
-
-#endif

tests/Test_dwf_mixedcg_prec.cc

@@ -25,9 +25,6 @@ Author: Peter Boyle <paboyle@ph.ed.ac.uk>
     See the full license in the file "LICENSE" in the top level distribution directory
     *************************************************************************************/
     /*  END LEGAL */
-#include "disable_tests_without_instantiations.h"
-#ifdef ENABLE_FERMION_INSTANTIATIONS
-
 #include <Grid/Grid.h>
 
 using namespace std;
@@ -225,5 +222,3 @@ int main (int argc, char ** argv)
   
   Grid_finalize();
 }
-
-#endif

tests/Test_dwf_mixedcg_prec_halfcomms.cc

@@ -25,9 +25,6 @@ Author: Peter Boyle <paboyle@ph.ed.ac.uk>
     See the full license in the file "LICENSE" in the top level distribution directory
     *************************************************************************************/
     /*  END LEGAL */
-#include "disable_tests_without_instantiations.h"
-#ifdef ENABLE_FERMION_INSTANTIATIONS
-
 #include <Grid/Grid.h>
 
 using namespace std;
@@ -121,4 +118,3 @@ int main (int argc, char ** argv)
   Grid_finalize();
 }
 #endif
-#endif

tests/Test_meson_field.cc

@@ -24,8 +24,6 @@ with this program; if not, write to the Free Software Foundation, Inc.,
 
 See the full license in the file "LICENSE" in the top level distribution directory
 *************************************************************************************/
-#include "disable_tests_without_instantiations.h"
-#ifdef ENABLE_FERMION_INSTANTIATIONS
 
 #include <Grid/Grid.h>
 #include <Grid/qcd/utils/A2Autils.h>
@@ -159,5 +157,3 @@ int main(int argc, char *argv[])
   
   return EXIT_SUCCESS;
 }
-
-#endif

tests/debug/Test_general_coarse_hdcg.cc

@@ -128,10 +128,6 @@ int main (int argc, char ** argv)
   typedef HermOpAdaptor<LatticeFermionD> HermFineMatrix;
   HermFineMatrix FineHermOp(HermOpEO);
 
-  LatticeFermionD src(FrbGrid); 
-  src = ComplexD(1.0);
-  PowerMethod<LatticeFermionD>       PM;   PM(HermOpEO,src);
-
   ////////////////////////////////////////////////////////////
   ///////////// Coarse basis and Little Dirac Operator ///////
   ////////////////////////////////////////////////////////////
@@ -154,7 +150,7 @@ int main (int argc, char ** argv)
   std::cout << "**************************************"<<std::endl;
   std::cout << "Create Subspace"<<std::endl;
   std::cout << "**************************************"<<std::endl;
-  Aggregates.CreateSubspaceChebyshev(RNG5,HermOpEO,nbasis,35.,0.01,500);// <== last run
+  Aggregates.CreateSubspaceChebyshevNew(RNG5,HermOpEO,95.); 
 
   std::cout << "**************************************"<<std::endl;
   std::cout << "Refine Subspace"<<std::endl;
@@ -189,7 +185,7 @@ int main (int argc, char ** argv)
   std::cout << "**************************************"<<std::endl;
 
   typedef HermitianLinearOperator<MultiGeneralCoarsenedMatrix_t,CoarseVector> MrhsHermMatrix;
-  Chebyshev<CoarseVector>      IRLCheby(0.01,16.0,201);  // 1 iter
+  Chebyshev<CoarseVector>      IRLCheby(0.05,40.0,101);  // 1 iter
   MrhsHermMatrix MrhsCoarseOp     (mrhs);
 
   CoarseVector pm_src(CoarseMrhs);
@@ -197,10 +193,10 @@ int main (int argc, char ** argv)
   PowerMethod<CoarseVector>       cPM;
   cPM(MrhsCoarseOp,pm_src);
 
-  //  int Nk=16;
+  int Nk=nrhs;
-  //  int Nm=Nk*3;
+  int Nm=Nk*3;
-  int Nk=32;
+  //  int Nk=36;
-  int Nm=128;
+  //  int Nm=144;
   int Nstop=Nk;
   int Nconv_test_interval=1;
   
@@ -214,7 +210,7 @@ int main (int argc, char ** argv)
 							  nrhs,
 							  Nk,
 							  Nm,
-							  1e-4,100);
+							  1e-4,10);
 
   int Nconv;
   std::vector<RealD>            eval(Nm);
@@ -235,6 +231,8 @@ int main (int argc, char ** argv)
   std::cout << "**************************************"<<std::endl;
   std::cout << " Recompute coarse evecs  "<<std::endl;
   std::cout << "**************************************"<<std::endl;
+  evec.resize(Nm,Coarse5d);
+  eval.resize(Nm);
   for(int r=0;r<nrhs;r++){
     random(CRNG,c_src[r]);
   }
@@ -245,7 +243,7 @@ int main (int argc, char ** argv)
   // Deflation guesser object
   ///////////////////////
   std::cout << "**************************************"<<std::endl;
-  std::cout << " Reimport coarse evecs "<<evec.size()<<" "<<eval.size()<<std::endl;
+  std::cout << " Reimport coarse evecs  "<<std::endl;
   std::cout << "**************************************"<<std::endl;
   MultiRHSDeflation<CoarseVector> MrhsGuesser;
   MrhsGuesser.ImportEigenBasis(evec,eval);
@@ -254,11 +252,9 @@ int main (int argc, char ** argv)
   // Extra HDCG parameters
   //////////////////////////
   int maxit=3000;
-  //  ConjugateGradient<CoarseVector>  CG(2.0e-1,maxit,false);
+  ConjugateGradient<CoarseVector>  CG(2.0e-1,maxit,false);
-  //  ConjugateGradient<CoarseVector>  CG(1.0e-2,maxit,false);
+  RealD lo=2.0;
-  ConjugateGradient<CoarseVector>  CG(5.0e-2,maxit,false);
+  int ord = 9;
-  RealD lo=0.2;
-  int ord = 7;
 
   DoNothingGuesser<CoarseVector> DoNothing;
   HPDSolver<CoarseVector> HPDSolveMrhs(MrhsCoarseOp,CG,DoNothing);
@@ -304,19 +300,6 @@ int main (int argc, char ** argv)
     ConjugateGradient<LatticeFermionD>  CGfine(1.0e-8,30000,false);
     CGfine(HermOpEO, src, result);
   }
-  {
-    std::cout << "**************************************"<<std::endl;
-    std::cout << "Calling MdagM CG"<<std::endl;
-    std::cout << "**************************************"<<std::endl;
-      
-    LatticeFermion result(FGrid); result=Zero();
-    LatticeFermion    src(FGrid); random(RNG5,src);
-    result=Zero();
-
-    MdagMLinearOperator<MobiusFermionD, LatticeFermionD> HermOp(Ddwf);
-    ConjugateGradient<LatticeFermionD>  CGfine(1.0e-8,30000,false);
-    CGfine(HermOp, src, result);
-  }
 #endif  
   Grid_finalize();
   return 0;

tests/debug/Test_general_coarse_pvdagm.cc

@@ -368,10 +368,7 @@ int main (int argc, char ** argv)
   TrivialPrecon<CoarseVector> simple;
   NonHermitianLinearOperator<LittleDiracOperator,CoarseVector> LinOpCoarse(LittleDiracOpPV);
   //  PrecGeneralisedConjugateResidualNonHermitian<CoarseVector>  L2PGCR(1.0e-4, 100, LinOpCoarse,simple,10,10); 
-  //  PrecGeneralisedConjugateResidualNonHermitian<CoarseVector>  L2PGCR(3.0e-2, 100, LinOpCoarse,simple,12,12);  // 35 outer
+  PrecGeneralisedConjugateResidualNonHermitian<CoarseVector>  L2PGCR(3.0e-2, 100, LinOpCoarse,simple,10,10); 
-  //  PrecGeneralisedConjugateResidualNonHermitian<CoarseVector>  L2PGCR(5.0e-2, 100, LinOpCoarse,simple,12,12);  // 36 outer, 12s
-  //  PrecGeneralisedConjugateResidualNonHermitian<CoarseVector>  L2PGCR(1.0e-1, 100, LinOpCoarse,simple,12,12);  // 36 ; 11s   
-  PrecGeneralisedConjugateResidualNonHermitian<CoarseVector>  L2PGCR(3.0e-1, 100, LinOpCoarse,simple,12,12);     
   L2PGCR.Level(3);
   c_res=Zero();
   L2PGCR(c_src,c_res);
@@ -403,7 +400,7 @@ int main (int argc, char ** argv)
 			    LinOpCoarse,
 			    L2PGCR);
   
-  PrecGeneralisedConjugateResidualNonHermitian<LatticeFermion> L1PGCR(1.0e-8,100,PVdagM,TwoLevelPrecon,10,10);
+  PrecGeneralisedConjugateResidualNonHermitian<LatticeFermion> L1PGCR(1.0e-8,1000,PVdagM,TwoLevelPrecon,16,16);
   L1PGCR.Level(1);
 
   f_res=Zero();

tests/debug/Test_general_coarse_pvdagm_svd.cc

@@ -1,493 +0,0 @@
-/*************************************************************************************
-
-    Grid physics library, www.github.com/paboyle/Grid 
-
-    Source file: ./tests/Test_padded_cell.cc
-
-    Copyright (C) 2023
-
-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>
-#include <Grid/lattice/PaddedCell.h>
-#include <Grid/stencil/GeneralLocalStencil.h>
-
-#include <Grid/algorithms/iterative/PrecGeneralisedConjugateResidual.h>
-#include <Grid/algorithms/iterative/PrecGeneralisedConjugateResidualNonHermitian.h>
-#include <Grid/algorithms/iterative/BiCGSTAB.h>
-
-using namespace std;
-using namespace Grid;
-
-template<class Matrix,class Field>
-class PVdagMLinearOperator : public LinearOperatorBase<Field> {
-  Matrix &_Mat;
-  Matrix &_PV;
-public:
-  PVdagMLinearOperator(Matrix &Mat,Matrix &PV): _Mat(Mat),_PV(PV){};
-
-  void OpDiag (const Field &in, Field &out) {    assert(0);  }
-  void OpDir  (const Field &in, Field &out,int dir,int disp) {    assert(0);  }
-  void OpDirAll  (const Field &in, std::vector<Field> &out){    assert(0);  };
-  void Op     (const Field &in, Field &out){
-    //    std::cout << GridLogMessage<< "Op: PVdag M "<<std::endl;
-    Field tmp(in.Grid());
-    _Mat.M(in,tmp);
-    _PV.Mdag(tmp,out);
-  }
-  void AdjOp     (const Field &in, Field &out){
-    //    std::cout << GridLogMessage<<"AdjOp: Mdag PV "<<std::endl;
-    Field tmp(in.Grid());
-    _PV.M(in,tmp);
-    _Mat.Mdag(tmp,out);
-  }
-  void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){
-    assert(0);
-  }
-  void HermOp(const Field &in, Field &out){
-    //    std::cout <<GridLogMessage<< "HermOp: Mdag PV PVdag M"<<std::endl;
-    Field tmp(in.Grid());
-    Op(in,tmp);
-    AdjOp(tmp,out);
-    //    std::cout << "HermOp done "<<norm2(out)<<std::endl;
-  }
-};
-template<class Matrix,class Field>
-class MdagPVLinearOperator : public LinearOperatorBase<Field> {
-  Matrix &_Mat;
-  Matrix &_PV;
-public:
-  MdagPVLinearOperator(Matrix &Mat,Matrix &PV): _Mat(Mat),_PV(PV){};
-
-  void OpDiag (const Field &in, Field &out) {    assert(0);  }
-  void OpDir  (const Field &in, Field &out,int dir,int disp) {    assert(0);  }
-  void OpDirAll  (const Field &in, std::vector<Field> &out){    assert(0);  };
-  void Op     (const Field &in, Field &out){
-    Field tmp(in.Grid());
-    //    std::cout <<GridLogMessage<< "Op: PVdag M "<<std::endl;
-    _PV.M(in,tmp);
-    _Mat.Mdag(tmp,out);
-  }
-  void AdjOp     (const Field &in, Field &out){
-    //    std::cout <<GridLogMessage<< "AdjOp: Mdag PV "<<std::endl;
-    Field tmp(in.Grid());
-    _Mat.M(in,tmp);
-    _PV.Mdag(tmp,out);
-  }
-  void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){
-    assert(0);
-  }
-  void HermOp(const Field &in, Field &out){
-    //    std::cout << GridLogMessage<<"HermOp: PVdag M Mdag PV "<<std::endl;
-    Field tmp(in.Grid());
-    Op(in,tmp);
-    AdjOp(tmp,out);
-    //    std::cout << "HermOp done "<<norm2(out)<<std::endl;
-  }
-};
-template<class Matrix,class Field>
-class ShiftedPVdagMLinearOperator : public LinearOperatorBase<Field> {
-  Matrix &_Mat;
-  Matrix &_PV;
-  RealD shift;
-public:
-  ShiftedPVdagMLinearOperator(RealD _shift,Matrix &Mat,Matrix &PV): shift(_shift),_Mat(Mat),_PV(PV){};
-
-  void OpDiag (const Field &in, Field &out) {    assert(0);  }
-  void OpDir  (const Field &in, Field &out,int dir,int disp) {    assert(0);  }
-  void OpDirAll  (const Field &in, std::vector<Field> &out){    assert(0);  };
-  void Op     (const Field &in, Field &out){
-    //    std::cout << "Op: PVdag M "<<std::endl;
-    Field tmp(in.Grid());
-    _Mat.M(in,tmp);
-    _PV.Mdag(tmp,out);
-    out = out + shift * in;
-  }
-  void AdjOp     (const Field &in, Field &out){
-    //    std::cout << "AdjOp: Mdag PV "<<std::endl;
-    Field tmp(in.Grid());
-    _PV.M(tmp,out);
-    _Mat.Mdag(in,tmp);
-    out = out + shift * in;
-  }
-  void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){    assert(0);  }
-  void HermOp(const Field &in, Field &out){
-    //    std::cout << "HermOp: Mdag PV PVdag M"<<std::endl;
-    Field tmp(in.Grid());
-    Op(in,tmp);
-    AdjOp(tmp,out);
-  }
-};
-template<class Fobj,class CComplex,int nbasis>
-class MGPreconditionerSVD : public LinearFunction< Lattice<Fobj> > {
-public:
-  using LinearFunction<Lattice<Fobj> >::operator();
-
-  typedef Aggregation<Fobj,CComplex,nbasis> Aggregates;
-  typedef typename Aggregation<Fobj,CComplex,nbasis>::FineField    FineField;
-  typedef typename Aggregation<Fobj,CComplex,nbasis>::CoarseVector CoarseVector;
-  typedef typename Aggregation<Fobj,CComplex,nbasis>::CoarseMatrix CoarseMatrix;
-  typedef LinearOperatorBase<FineField>                            FineOperator;
-  typedef LinearFunction    <FineField>                            FineSmoother;
-  typedef LinearOperatorBase<CoarseVector>                         CoarseOperator;
-  typedef LinearFunction    <CoarseVector>                         CoarseSolver;
-  Aggregates     & _FineToCoarse;
-  Aggregates     & _CoarseToFine;
-  FineOperator   & _FineOperator;
-  FineSmoother   & _PreSmoother;
-  FineSmoother   & _PostSmoother;
-  CoarseOperator & _CoarseOperator;
-  CoarseSolver   & _CoarseSolve;
-
-  int    level;  void Level(int lv) {level = lv; };
-
-  MGPreconditionerSVD(Aggregates &FtoC,
-		      Aggregates &CtoF,
-		      FineOperator &Fine,
-		      FineSmoother &PreSmoother,
-		      FineSmoother &PostSmoother,
-		      CoarseOperator &CoarseOperator_,
-		      CoarseSolver &CoarseSolve_)
-    : _FineToCoarse(FtoC),
-      _CoarseToFine(CtoF),
-      _FineOperator(Fine),
-      _PreSmoother(PreSmoother),
-      _PostSmoother(PostSmoother),
-      _CoarseOperator(CoarseOperator_),
-      _CoarseSolve(CoarseSolve_),
-      level(1)  {  }
-
-  virtual void operator()(const FineField &in, FineField & out) 
-  {
-    GridBase *CoarseGrid = _FineToCoarse.CoarseGrid;
-    //    auto CoarseGrid = _CoarseOperator.Grid();
-    CoarseVector Csrc(CoarseGrid);
-    CoarseVector Csol(CoarseGrid);
-    FineField vec1(in.Grid());
-    FineField vec2(in.Grid());
-
-    std::cout<<GridLogMessage << "Calling PreSmoother " <<std::endl;
-
-    //    std::cout<<GridLogMessage << "Calling PreSmoother input residual "<<norm2(in) <<std::endl;
-    double t;
-    // Fine Smoother
-    //    out = in;
-    out = Zero();
-    t=-usecond();
-    _PreSmoother(in,out);
-    t+=usecond();
-
-    std::cout<<GridLogMessage << "PreSmoother took "<< t/1000.0<< "ms" <<std::endl;
-
-    // Update the residual
-    _FineOperator.Op(out,vec1);  sub(vec1, in ,vec1);   
-    //    std::cout<<GridLogMessage <<"Residual-1 now " <<norm2(vec1)<<std::endl;
-
-    // Fine to Coarse 
-    t=-usecond();
-    _FineToCoarse.ProjectToSubspace  (Csrc,vec1);
-    t+=usecond();
-    std::cout<<GridLogMessage << "Project to coarse took "<< t/1000.0<< "ms" <<std::endl;
-
-    // Coarse correction
-    t=-usecond();
-    Csol = Zero();
-    _CoarseSolve(Csrc,Csol);
-    //Csol=Zero();
-    t+=usecond();
-    std::cout<<GridLogMessage << "Coarse solve took "<< t/1000.0<< "ms" <<std::endl;
-
-    // Coarse to Fine
-    t=-usecond();  
-    //    _CoarseOperator.PromoteFromSubspace(_Aggregates,Csol,vec1);
-    _CoarseToFine.PromoteFromSubspace(Csol,vec1); 
-    add(out,out,vec1);
-    t+=usecond();
-    std::cout<<GridLogMessage << "Promote to this level took "<< t/1000.0<< "ms" <<std::endl;
-
-    // Residual
-    _FineOperator.Op(out,vec1);  sub(vec1 ,in , vec1);  
-    //    std::cout<<GridLogMessage <<"Residual-2 now " <<norm2(vec1)<<std::endl;
-
-    // Fine Smoother
-    t=-usecond();
-    //    vec2=vec1;
-    vec2=Zero();
-    _PostSmoother(vec1,vec2);
-    t+=usecond();
-    std::cout<<GridLogMessage << "PostSmoother took "<< t/1000.0<< "ms" <<std::endl;
-
-    add( out,out,vec2);
-    std::cout<<GridLogMessage << "Done " <<std::endl;
-  }
-};
-
-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
-  Coordinate clatt = GridDefaultLatt();
-  for(int d=0;d<clatt.size();d++){
-    clatt[d] = clatt[d]/2;
-    //    clatt[d] = clatt[d]/4;
-  }
-  GridCartesian *Coarse4d =  SpaceTimeGrid::makeFourDimGrid(clatt, GridDefaultSimd(Nd,vComplex::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); random(RNG5,src);
-  LatticeFermion result(FGrid); result=Zero();
-  LatticeFermion    ref(FGrid); ref=Zero();
-  LatticeFermion    tmp(FGrid);
-  LatticeFermion    err(FGrid);
-  LatticeGaugeField Umu(UGrid);
-
-  FieldMetaData header;
-  std::string file("ckpoint_lat.4000");
-  NerscIO::readConfiguration(Umu,header,file);
-  
-  RealD mass=0.01;
-  RealD M5=1.8;
-
-  DomainWallFermionD Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5);
-  DomainWallFermionD Dpv(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,1.0,M5);
-
-  const int nbasis = 30;
-  const int cb = 0 ;
-
-
-  NextToNearestStencilGeometry5D geom(Coarse5d);
-
-  std::cout<<GridLogMessage<<std::endl;
-  std::cout<<GridLogMessage<<"*******************************************"<<std::endl;
-  std::cout<<GridLogMessage<<std::endl;
-
-  typedef PVdagMLinearOperator<DomainWallFermionD,LatticeFermionD> PVdagM_t;
-  typedef MdagPVLinearOperator<DomainWallFermionD,LatticeFermionD> MdagPV_t;
-  typedef ShiftedPVdagMLinearOperator<DomainWallFermionD,LatticeFermionD> ShiftedPVdagM_t;
-  PVdagM_t PVdagM(Ddwf,Dpv);
-  MdagPV_t MdagPV(Ddwf,Dpv);
-  //  ShiftedPVdagM_t ShiftedPVdagM(2.0,Ddwf,Dpv); // 355
-  //  ShiftedPVdagM_t ShiftedPVdagM(1.0,Ddwf,Dpv); // 246
-  //  ShiftedPVdagM_t ShiftedPVdagM(0.5,Ddwf,Dpv); // 183
-  //  ShiftedPVdagM_t ShiftedPVdagM(0.25,Ddwf,Dpv); // 145
-  //  ShiftedPVdagM_t ShiftedPVdagM(0.1,Ddwf,Dpv); // 134
-  //  ShiftedPVdagM_t ShiftedPVdagM(0.1,Ddwf,Dpv); // 127 -- NULL space via inverse iteration
-  //  ShiftedPVdagM_t ShiftedPVdagM(0.1,Ddwf,Dpv); // 57 -- NULL space via inverse iteration; 3 iterations
-  //  ShiftedPVdagM_t ShiftedPVdagM(0.25,Ddwf,Dpv); // 57 , tighter inversion
-  //  ShiftedPVdagM_t ShiftedPVdagM(0.25,Ddwf,Dpv); // nbasis 20 -- 49 iters
-  //  ShiftedPVdagM_t ShiftedPVdagM(0.25,Ddwf,Dpv); // nbasis 20 -- 70 iters; asymmetric 
-  //  ShiftedPVdagM_t ShiftedPVdagM(0.25,Ddwf,Dpv); // 58; Loosen coarse, tighten fine
-  //  ShiftedPVdagM_t ShiftedPVdagM(0.1,Ddwf,Dpv); // 56 ... 
-  //  ShiftedPVdagM_t ShiftedPVdagM(0.1,Ddwf,Dpv); // 51 ...  with 24 vecs
-  //  ShiftedPVdagM_t ShiftedPVdagM(0.1,Ddwf,Dpv); // 31 ...  with 24 vecs and 2^4 blocking
-  //  ShiftedPVdagM_t ShiftedPVdagM(0.1,Ddwf,Dpv); // 43 ...  with 16 vecs and 2^4 blocking, sloppier
-  //  ShiftedPVdagM_t ShiftedPVdagM(0.1,Ddwf,Dpv); // 35  ...  with 20 vecs and 2^4 blocking
-  //  ShiftedPVdagM_t ShiftedPVdagM(0.1,Ddwf,Dpv); // 35  ...  with 20 vecs and 2^4 blocking, looser coarse
-  //  ShiftedPVdagM_t ShiftedPVdagM(0.1,Ddwf,Dpv); // 64  ...  with 20 vecs, Christoph setup, and 2^4 blocking, looser coarse
-  ShiftedPVdagM_t ShiftedPVdagM(0.01,Ddwf,Dpv); // 
-
-
-  // Run power method on HOA??
-  PowerMethod<LatticeFermion>       PM;
-  //  PM(PVdagM,src);
-  //  PM(MdagPV,src);
- 
-  // Warning: This routine calls PVdagM.Op, not PVdagM.HermOp
-  typedef Aggregation<vSpinColourVector,vTComplex,nbasis> Subspace;
-  Subspace V(Coarse5d,FGrid,cb);
-  Subspace U(Coarse5d,FGrid,cb);
-
-  // Breeds right singular vectors with call to HermOp (V)
-  V.CreateSubspaceChebyshev(RNG5,PVdagM,
-			    nbasis,
-			    4000.0,0.003,
-			    500);
-
-  // Breeds left singular vectors with call to HermOp (U)
-  //  U.CreateSubspaceChebyshev(RNG5,PVdagM,
-  U.CreateSubspaceChebyshev(RNG5,MdagPV,
-			    nbasis,
-			    4000.0,0.003,
-			    500);
-  
-  typedef Aggregation<vSpinColourVector,vTComplex,2*nbasis> CombinedSubspace;
-  CombinedSubspace CombinedUV(Coarse5d,FGrid,cb);
-  for(int b=0;b<nbasis;b++){
-    CombinedUV.subspace[b]        = V.subspace[b];
-    CombinedUV.subspace[b+nbasis] = U.subspace[b];
-  }
-
-  int bl, br;
-  std::cout <<" <V| PVdagM| V> " <<std::endl;
-  for(bl=0;bl<nbasis;bl++){
-  for(br=0;br<nbasis;br++){
-    PVdagM.Op(V.subspace[br],src);
-    std::cout <<bl<<" "<<br<<"\t"<<innerProduct(V.subspace[bl],src)<<std::endl;
-  }}
-  std::cout <<" <V| PVdagM| U> " <<std::endl;
-  for(bl=0;bl<nbasis;bl++){
-  for(br=0;br<nbasis;br++){
-    PVdagM.Op(U.subspace[br],src);
-    std::cout <<bl<<" "<<br<<"\t"<<innerProduct(V.subspace[bl],src)<<std::endl;
-  }}
-  std::cout <<" <U| PVdagM| V> " <<std::endl;
-  for(bl=0;bl<nbasis;bl++){
-  for(br=0;br<nbasis;br++){
-    PVdagM.Op(V.subspace[br],src);
-    std::cout <<bl<<" "<<br<<"\t"<<innerProduct(U.subspace[bl],src)<<std::endl;
-  }}
-  std::cout <<" <U| PVdagM| U> " <<std::endl;
-  for(bl=0;bl<nbasis;bl++){
-  for(br=0;br<nbasis;br++){
-    PVdagM.Op(U.subspace[br],src);
-    std::cout <<bl<<" "<<br<<"\t"<<innerProduct(U.subspace[bl],src)<<std::endl;
-  }}
-
-  typedef GeneralCoarsenedMatrix<vSpinColourVector,vTComplex,nbasis> LittleDiracOperatorV;
-  typedef LittleDiracOperatorV::CoarseVector CoarseVectorV;
-  
-  typedef GeneralCoarsenedMatrix<vSpinColourVector,vTComplex,2*nbasis> LittleDiracOperator;
-  typedef LittleDiracOperator::CoarseVector CoarseVector;
-
-  V.Orthogonalise();
-  for(int b =0 ; b<nbasis;b++){
-    CoarseVectorV c_src (Coarse5d);
-    V.ProjectToSubspace  (c_src,U.subspace[b]);
-    V.PromoteFromSubspace(c_src,src);
-    std::cout << " Completeness of U in V ["<< b<<"] "<< std::sqrt(norm2(src)/norm2(U.subspace[b]))<<std::endl;
-  }
-  
-  CoarseVector c_src (Coarse5d);
-  CoarseVector c_res (Coarse5d);
-  CoarseVector c_proj(Coarse5d);
-  LittleDiracOperator LittleDiracOpPV(geom,FGrid,Coarse5d);
-  LittleDiracOpPV.CoarsenOperator(PVdagM,CombinedUV,CombinedUV);
-
-  std::cout<<GridLogMessage<<std::endl;
-  std::cout<<GridLogMessage<<"*******************************************"<<std::endl;
-  std::cout<<GridLogMessage<<std::endl;
-  std::cout<<GridLogMessage<<"Testing coarsened operator "<<std::endl;
-  
-  Complex one(1.0);
-  c_src = one;  // 1 in every element for vector 1.
-
-  blockPromote(c_src,err,CombinedUV.subspace);
-
-  LatticeFermion prom(FGrid);
-  prom=Zero();
-  for(int b=0;b<nbasis*2;b++){
-    prom=prom+CombinedUV.subspace[b];
-  }
-
-  std::cout<<GridLogMessage<<"c_src "<<norm2(c_src)<<std::endl;
-  std::cout<<GridLogMessage<<"prom  "<<norm2(prom)<<std::endl;
-
-  PVdagM.Op(prom,tmp);
-  blockProject(c_proj,tmp,CombinedUV.subspace);
-  std::cout<<GridLogMessage<<" Called Big Dirac Op "<<norm2(tmp)<<std::endl;
-
-  LittleDiracOpPV.M(c_src,c_res);
-  std::cout<<GridLogMessage<<" Called Little Dirac Op c_src "<< norm2(c_src) << "  c_res "<< norm2(c_res) <<std::endl;
-
-  std::cout<<GridLogMessage<<"Little dop : "<<norm2(c_res)<<std::endl;
-
-  std::cout<<GridLogMessage<<"Big dop in subspace : "<<norm2(c_proj)<<std::endl;
-
-  c_proj = c_proj - c_res;
-  std::cout<<GridLogMessage<<" ldop error: "<<norm2(c_proj)<<std::endl;
-
-
-  /**********
-   * Some solvers
-   **********
-   */
-
-  ///////////////////////////////////////
-  // Coarse grid solver test
-  ///////////////////////////////////////
-
-  std::cout<<GridLogMessage<<"******************* "<<std::endl;
-  std::cout<<GridLogMessage<<" Coarse Grid Solve -- Level 3 "<<std::endl;
-  std::cout<<GridLogMessage<<"******************* "<<std::endl;
-  TrivialPrecon<CoarseVector> simple;
-  NonHermitianLinearOperator<LittleDiracOperator,CoarseVector> LinOpCoarse(LittleDiracOpPV);
-  //  PrecGeneralisedConjugateResidualNonHermitian<CoarseVector>  L2PGCR(1.0e-4, 100, LinOpCoarse,simple,10,10); 
-  PrecGeneralisedConjugateResidualNonHermitian<CoarseVector>  L2PGCR(1.0e-2, 10, LinOpCoarse,simple,20,20); 
-  L2PGCR.Level(3);
-  c_res=Zero();
-  L2PGCR(c_src,c_res);
-
-  ////////////////////////////////////////
-  // Fine grid smoother
-  ////////////////////////////////////////
-  std::cout<<GridLogMessage<<"******************* "<<std::endl;
-  std::cout<<GridLogMessage<<" Fine Grid Smoother -- Level 2 "<<std::endl;
-  std::cout<<GridLogMessage<<"******************* "<<std::endl;
-  TrivialPrecon<LatticeFermionD> simple_fine;
-  //  NonHermitianLinearOperator<PVdagM_t,LatticeFermionD> LinOpSmooth(PVdagM);
-  PrecGeneralisedConjugateResidualNonHermitian<LatticeFermionD> SmootherGCR(0.01,1,ShiftedPVdagM,simple_fine,16,16);
-  SmootherGCR.Level(2);
-  
-  LatticeFermionD f_src(FGrid);
-  LatticeFermionD f_res(FGrid);
-
-  f_src = one;  // 1 in every element for vector 1.
-  f_res=Zero();
-  SmootherGCR(f_src,f_res);
-
-  typedef MGPreconditionerSVD<vSpinColourVector,  vTComplex,nbasis*2> TwoLevelMG;
-
-  TwoLevelMG TwoLevelPrecon(CombinedUV,CombinedUV,
-			    PVdagM,
-			    simple_fine,
-			    SmootherGCR,
-			    LinOpCoarse,
-			    L2PGCR);
-  
-  PrecGeneralisedConjugateResidualNonHermitian<LatticeFermion> L1PGCR(1.0e-8,1000,PVdagM,TwoLevelPrecon,20,20);
-  L1PGCR.Level(1);
-
-  f_res=Zero();
-  L1PGCR(f_src,f_res);
-
-  std::cout<<GridLogMessage<<std::endl;
-  std::cout<<GridLogMessage<<"*******************************************"<<std::endl;
-  std::cout<<GridLogMessage<<std::endl;
-  std::cout<<GridLogMessage << "Done "<< std::endl;
-
-  Grid_finalize();
-  return 0;
-}

tests/debug/Test_general_coarse_pvdagm_svd_cg.cc

@@ -1,492 +0,0 @@
-/*************************************************************************************
-
-    Grid physics library, www.github.com/paboyle/Grid 
-
-    Source file: ./tests/Test_padded_cell.cc
-
-    Copyright (C) 2023
-
-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>
-#include <Grid/lattice/PaddedCell.h>
-#include <Grid/stencil/GeneralLocalStencil.h>
-
-#include <Grid/algorithms/iterative/PrecGeneralisedConjugateResidual.h>
-#include <Grid/algorithms/iterative/PrecGeneralisedConjugateResidualNonHermitian.h>
-#include <Grid/algorithms/iterative/BiCGSTAB.h>
-
-using namespace std;
-using namespace Grid;
-
-template<class Matrix,class Field>
-class PVdagMLinearOperator : public LinearOperatorBase<Field> {
-  Matrix &_Mat;
-  Matrix &_PV;
-public:
-  PVdagMLinearOperator(Matrix &Mat,Matrix &PV): _Mat(Mat),_PV(PV){};
-
-  void OpDiag (const Field &in, Field &out) {    assert(0);  }
-  void OpDir  (const Field &in, Field &out,int dir,int disp) {    assert(0);  }
-  void OpDirAll  (const Field &in, std::vector<Field> &out){    assert(0);  };
-  void Op     (const Field &in, Field &out){
-    //    std::cout << GridLogMessage<< "Op: PVdag M "<<std::endl;
-    Field tmp(in.Grid());
-    _Mat.M(in,tmp);
-    _PV.Mdag(tmp,out);
-  }
-  void AdjOp     (const Field &in, Field &out){
-    //    std::cout << GridLogMessage<<"AdjOp: Mdag PV "<<std::endl;
-    Field tmp(in.Grid());
-    _PV.M(in,tmp);
-    _Mat.Mdag(tmp,out);
-  }
-  void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){
-    HermOp(in,out);
-    ComplexD dot = innerProduct(in,out);
-    n1=real(dot);
-    n2=norm2(out);
-  }
-  void HermOp(const Field &in, Field &out){
-    //    std::cout <<GridLogMessage<< "HermOp: Mdag PV PVdag M"<<std::endl;
-    Field tmp(in.Grid());
-    Op(in,tmp);
-    AdjOp(tmp,out);
-    //    std::cout << "HermOp done "<<norm2(out)<<std::endl;
-  }
-};
-template<class Matrix,class Field>
-class MdagPVLinearOperator : public LinearOperatorBase<Field> {
-  Matrix &_Mat;
-  Matrix &_PV;
-public:
-  MdagPVLinearOperator(Matrix &Mat,Matrix &PV): _Mat(Mat),_PV(PV){};
-
-  void OpDiag (const Field &in, Field &out) {    assert(0);  }
-  void OpDir  (const Field &in, Field &out,int dir,int disp) {    assert(0);  }
-  void OpDirAll  (const Field &in, std::vector<Field> &out){    assert(0);  };
-  void Op     (const Field &in, Field &out){
-    Field tmp(in.Grid());
-    //    std::cout <<GridLogMessage<< "Op: PVdag M "<<std::endl;
-    _PV.M(in,tmp);
-    _Mat.Mdag(tmp,out);
-  }
-  void AdjOp     (const Field &in, Field &out){
-    //    std::cout <<GridLogMessage<< "AdjOp: Mdag PV "<<std::endl;
-    Field tmp(in.Grid());
-    _Mat.M(in,tmp);
-    _PV.Mdag(tmp,out);
-  }
-  void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){
-    ComplexD dot = innerProduct(in,out);
-    n1=real(dot);
-    n2=norm2(out);
-  }
-  void HermOp(const Field &in, Field &out){
-    //    std::cout << GridLogMessage<<"HermOp: PVdag M Mdag PV "<<std::endl;
-    Field tmp(in.Grid());
-    Op(in,tmp);
-    AdjOp(tmp,out);
-    //    std::cout << "HermOp done "<<norm2(out)<<std::endl;
-  }
-};
-template<class Matrix,class Field>
-class ShiftedPVdagMLinearOperator : public LinearOperatorBase<Field> {
-  Matrix &_Mat;
-  Matrix &_PV;
-  RealD shift;
-public:
-  ShiftedPVdagMLinearOperator(RealD _shift,Matrix &Mat,Matrix &PV): shift(_shift),_Mat(Mat),_PV(PV){};
-
-  void OpDiag (const Field &in, Field &out) {    assert(0);  }
-  void OpDir  (const Field &in, Field &out,int dir,int disp) {    assert(0);  }
-  void OpDirAll  (const Field &in, std::vector<Field> &out){    assert(0);  };
-  void Op     (const Field &in, Field &out){
-    //    std::cout << "Op: PVdag M "<<std::endl;
-    Field tmp(in.Grid());
-    _Mat.M(in,tmp);
-    _PV.Mdag(tmp,out);
-    out = out + shift * in;
-  }
-  void AdjOp     (const Field &in, Field &out){
-    //    std::cout << "AdjOp: Mdag PV "<<std::endl;
-    Field tmp(in.Grid());
-    _PV.M(tmp,out);
-    _Mat.Mdag(in,tmp);
-    out = out + shift * in;
-  }
-  void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){    assert(0);  }
-  void HermOp(const Field &in, Field &out){
-    //    std::cout << "HermOp: Mdag PV PVdag M"<<std::endl;
-    Field tmp(in.Grid());
-    Op(in,tmp);
-    AdjOp(tmp,out);
-  }
-};
-template<class Fobj,class CComplex,int nbasis>
-class MGPreconditionerSVD : public LinearFunction< Lattice<Fobj> > {
-public:
-  using LinearFunction<Lattice<Fobj> >::operator();
-
-  typedef Aggregation<Fobj,CComplex,nbasis> Aggregates;
-  typedef typename Aggregation<Fobj,CComplex,nbasis>::FineField    FineField;
-  typedef typename Aggregation<Fobj,CComplex,nbasis>::CoarseVector CoarseVector;
-  typedef typename Aggregation<Fobj,CComplex,nbasis>::CoarseMatrix CoarseMatrix;
-  typedef LinearOperatorBase<FineField>                            FineOperator;
-  typedef LinearFunction    <FineField>                            FineSmoother;
-  typedef LinearOperatorBase<CoarseVector>                         CoarseOperator;
-  typedef LinearFunction    <CoarseVector>                         CoarseSolver;
-  Aggregates     & _FineToCoarse;
-  Aggregates     & _CoarseToFine;
-  FineOperator   & _FineOperator;
-  FineSmoother   & _PreSmoother;
-  FineSmoother   & _PostSmoother;
-  CoarseOperator & _CoarseOperator;
-  CoarseSolver   & _CoarseSolve;
-
-  int    level;  void Level(int lv) {level = lv; };
-
-  MGPreconditionerSVD(Aggregates &FtoC,
-		      Aggregates &CtoF,
-		      FineOperator &Fine,
-		      FineSmoother &PreSmoother,
-		      FineSmoother &PostSmoother,
-		      CoarseOperator &CoarseOperator_,
-		      CoarseSolver &CoarseSolve_)
-    : _FineToCoarse(FtoC),
-      _CoarseToFine(CtoF),
-      _FineOperator(Fine),
-      _PreSmoother(PreSmoother),
-      _PostSmoother(PostSmoother),
-      _CoarseOperator(CoarseOperator_),
-      _CoarseSolve(CoarseSolve_),
-      level(1)  {  }
-
-  virtual void operator()(const FineField &in, FineField & out) 
-  {
-    GridBase *CoarseGrid = _FineToCoarse.CoarseGrid;
-    //    auto CoarseGrid = _CoarseOperator.Grid();
-    CoarseVector Csrc(CoarseGrid);
-    CoarseVector Csol(CoarseGrid);
-    FineField vec1(in.Grid());
-    FineField vec2(in.Grid());
-
-    std::cout<<GridLogMessage << "Calling PreSmoother " <<std::endl;
-
-    //    std::cout<<GridLogMessage << "Calling PreSmoother input residual "<<norm2(in) <<std::endl;
-    double t;
-    // Fine Smoother
-    //    out = in;
-    out = Zero();
-    t=-usecond();
-    _PreSmoother(in,out);
-    t+=usecond();
-
-    std::cout<<GridLogMessage << "PreSmoother took "<< t/1000.0<< "ms" <<std::endl;
-
-    // Update the residual
-    _FineOperator.Op(out,vec1);  sub(vec1, in ,vec1);   
-    //    std::cout<<GridLogMessage <<"Residual-1 now " <<norm2(vec1)<<std::endl;
-
-    // Fine to Coarse 
-    t=-usecond();
-    _FineToCoarse.ProjectToSubspace  (Csrc,vec1);
-    t+=usecond();
-    std::cout<<GridLogMessage << "Project to coarse took "<< t/1000.0<< "ms" <<std::endl;
-
-    // Coarse correction
-    t=-usecond();
-    Csol = Zero();
-    _CoarseSolve(Csrc,Csol);
-    //Csol=Zero();
-    t+=usecond();
-    std::cout<<GridLogMessage << "Coarse solve took "<< t/1000.0<< "ms" <<std::endl;
-
-    // Coarse to Fine
-    t=-usecond();  
-    //    _CoarseOperator.PromoteFromSubspace(_Aggregates,Csol,vec1);
-    _CoarseToFine.PromoteFromSubspace(Csol,vec1); 
-    add(out,out,vec1);
-    t+=usecond();
-    std::cout<<GridLogMessage << "Promote to this level took "<< t/1000.0<< "ms" <<std::endl;
-
-    // Residual
-    _FineOperator.Op(out,vec1);  sub(vec1 ,in , vec1);  
-    //    std::cout<<GridLogMessage <<"Residual-2 now " <<norm2(vec1)<<std::endl;
-
-    // Fine Smoother
-    t=-usecond();
-    //    vec2=vec1;
-    vec2=Zero();
-    _PostSmoother(vec1,vec2);
-    t+=usecond();
-    std::cout<<GridLogMessage << "PostSmoother took "<< t/1000.0<< "ms" <<std::endl;
-
-    add( out,out,vec2);
-    std::cout<<GridLogMessage << "Done " <<std::endl;
-  }
-};
-
-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
-  Coordinate clatt = GridDefaultLatt();
-  for(int d=0;d<clatt.size();d++){
-    clatt[d] = clatt[d]/2;
-    //    clatt[d] = clatt[d]/4;
-  }
-  GridCartesian *Coarse4d =  SpaceTimeGrid::makeFourDimGrid(clatt, GridDefaultSimd(Nd,vComplex::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); random(RNG5,src);
-  LatticeFermion result(FGrid); result=Zero();
-  LatticeFermion    ref(FGrid); ref=Zero();
-  LatticeFermion    tmp(FGrid);
-  LatticeFermion    err(FGrid);
-  LatticeGaugeField Umu(UGrid);
-
-  FieldMetaData header;
-  std::string file("ckpoint_lat.4000");
-  NerscIO::readConfiguration(Umu,header,file);
-  
-  RealD mass=0.01;
-  RealD M5=1.8;
-
-  DomainWallFermionD Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5);
-  DomainWallFermionD Dpv(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,1.0,M5);
-
-  const int nbasis = 20;
-  const int cb = 0 ;
-
-
-  NextToNearestStencilGeometry5D geom(Coarse5d);
-
-  std::cout<<GridLogMessage<<std::endl;
-  std::cout<<GridLogMessage<<"*******************************************"<<std::endl;
-  std::cout<<GridLogMessage<<std::endl;
-
-  typedef PVdagMLinearOperator<DomainWallFermionD,LatticeFermionD> PVdagM_t;
-  typedef MdagPVLinearOperator<DomainWallFermionD,LatticeFermionD> MdagPV_t;
-  typedef ShiftedPVdagMLinearOperator<DomainWallFermionD,LatticeFermionD> ShiftedPVdagM_t;
-  PVdagM_t PVdagM(Ddwf,Dpv);
-  MdagPV_t MdagPV(Ddwf,Dpv);
-  //  ShiftedPVdagM_t ShiftedPVdagM(2.0,Ddwf,Dpv); // 355
-  //  ShiftedPVdagM_t ShiftedPVdagM(1.0,Ddwf,Dpv); // 246
-  //  ShiftedPVdagM_t ShiftedPVdagM(0.5,Ddwf,Dpv); // 183
-  //  ShiftedPVdagM_t ShiftedPVdagM(0.25,Ddwf,Dpv); // 145
-  //  ShiftedPVdagM_t ShiftedPVdagM(0.1,Ddwf,Dpv); // 134
-  //  ShiftedPVdagM_t ShiftedPVdagM(0.1,Ddwf,Dpv); // 127 -- NULL space via inverse iteration
-  //  ShiftedPVdagM_t ShiftedPVdagM(0.1,Ddwf,Dpv); // 57 -- NULL space via inverse iteration; 3 iterations
-  //  ShiftedPVdagM_t ShiftedPVdagM(0.25,Ddwf,Dpv); // 57 , tighter inversion
-  //  ShiftedPVdagM_t ShiftedPVdagM(0.25,Ddwf,Dpv); // nbasis 20 -- 49 iters
-  //  ShiftedPVdagM_t ShiftedPVdagM(0.25,Ddwf,Dpv); // nbasis 20 -- 70 iters; asymmetric 
-  //  ShiftedPVdagM_t ShiftedPVdagM(0.25,Ddwf,Dpv); // 58; Loosen coarse, tighten fine
-  //  ShiftedPVdagM_t ShiftedPVdagM(0.1,Ddwf,Dpv); // 56 ... 
-  //  ShiftedPVdagM_t ShiftedPVdagM(0.1,Ddwf,Dpv); // 51 ...  with 24 vecs
-  //  ShiftedPVdagM_t ShiftedPVdagM(0.1,Ddwf,Dpv); // 31 ...  with 24 vecs and 2^4 blocking
-  //  ShiftedPVdagM_t ShiftedPVdagM(0.1,Ddwf,Dpv); // 43 ...  with 16 vecs and 2^4 blocking, sloppier
-  //  ShiftedPVdagM_t ShiftedPVdagM(0.1,Ddwf,Dpv); // 35  ...  with 20 vecs and 2^4 blocking
-  //  ShiftedPVdagM_t ShiftedPVdagM(0.1,Ddwf,Dpv); // 35  ...  with 20 vecs and 2^4 blocking, looser coarse
-  //  ShiftedPVdagM_t ShiftedPVdagM(0.1,Ddwf,Dpv); // 64  ...  with 20 vecs, Christoph setup, and 2^4 blocking, looser coarse
-  ShiftedPVdagM_t ShiftedPVdagM(0.01,Ddwf,Dpv); // 
-
-
-  // Run power method on HOA??
-  PowerMethod<LatticeFermion>       PM;
-  //  PM(PVdagM,src);
-  //  PM(MdagPV,src);
- 
-  // Warning: This routine calls PVdagM.Op, not PVdagM.HermOp
-  typedef Aggregation<vSpinColourVector,vTComplex,nbasis> Subspace;
-  Subspace V(Coarse5d,FGrid,cb);
-  Subspace U(Coarse5d,FGrid,cb);
-
-  // Breeds right singular vectors with call to HermOp (V)
-  V.CreateSubspace(RNG5,PVdagM,nbasis);
-
-  // Breeds left singular vectors with call to HermOp (U)
-  //  U.CreateSubspaceChebyshev(RNG5,MdagPV,
-  U.CreateSubspace(RNG5,PVdagM,nbasis);
-  
-  typedef Aggregation<vSpinColourVector,vTComplex,2*nbasis> CombinedSubspace;
-  CombinedSubspace CombinedUV(Coarse5d,FGrid,cb);
-  for(int b=0;b<nbasis;b++){
-    CombinedUV.subspace[b]        = V.subspace[b];
-    CombinedUV.subspace[b+nbasis] = U.subspace[b];
-  }
-
-  int bl, br;
-  std::cout <<" <V| PVdagM| V> " <<std::endl;
-  for(bl=0;bl<nbasis;bl++){
-  for(br=0;br<nbasis;br++){
-    PVdagM.Op(V.subspace[br],src);
-    std::cout <<bl<<" "<<br<<"\t"<<innerProduct(V.subspace[bl],src)<<std::endl;
-  }}
-  std::cout <<" <V| PVdagM| U> " <<std::endl;
-  for(bl=0;bl<nbasis;bl++){
-  for(br=0;br<nbasis;br++){
-    PVdagM.Op(U.subspace[br],src);
-    std::cout <<bl<<" "<<br<<"\t"<<innerProduct(V.subspace[bl],src)<<std::endl;
-  }}
-  std::cout <<" <U| PVdagM| V> " <<std::endl;
-  for(bl=0;bl<nbasis;bl++){
-  for(br=0;br<nbasis;br++){
-    PVdagM.Op(V.subspace[br],src);
-    std::cout <<bl<<" "<<br<<"\t"<<innerProduct(U.subspace[bl],src)<<std::endl;
-  }}
-  std::cout <<" <U| PVdagM| U> " <<std::endl;
-  for(bl=0;bl<nbasis;bl++){
-  for(br=0;br<nbasis;br++){
-    PVdagM.Op(U.subspace[br],src);
-    std::cout <<bl<<" "<<br<<"\t"<<innerProduct(U.subspace[bl],src)<<std::endl;
-  }}
-
-  typedef GeneralCoarsenedMatrix<vSpinColourVector,vTComplex,nbasis> LittleDiracOperatorV;
-  typedef LittleDiracOperatorV::CoarseVector CoarseVectorV;
-  
-  typedef GeneralCoarsenedMatrix<vSpinColourVector,vTComplex,2*nbasis> LittleDiracOperator;
-  typedef LittleDiracOperator::CoarseVector CoarseVector;
-
-  V.Orthogonalise();
-  for(int b =0 ; b<nbasis;b++){
-    CoarseVectorV c_src (Coarse5d);
-    V.ProjectToSubspace  (c_src,U.subspace[b]);
-    V.PromoteFromSubspace(c_src,src);
-    std::cout << " Completeness of U in V ["<< b<<"] "<< std::sqrt(norm2(src)/norm2(U.subspace[b]))<<std::endl;
-  }
-  
-  CoarseVector c_src (Coarse5d);
-  CoarseVector c_res (Coarse5d);
-  CoarseVector c_proj(Coarse5d);
-  LittleDiracOperator LittleDiracOpPV(geom,FGrid,Coarse5d);
-  LittleDiracOpPV.CoarsenOperator(PVdagM,CombinedUV,CombinedUV);
-
-  std::cout<<GridLogMessage<<std::endl;
-  std::cout<<GridLogMessage<<"*******************************************"<<std::endl;
-  std::cout<<GridLogMessage<<std::endl;
-  std::cout<<GridLogMessage<<"Testing coarsened operator "<<std::endl;
-  
-  Complex one(1.0);
-  c_src = one;  // 1 in every element for vector 1.
-
-  blockPromote(c_src,err,CombinedUV.subspace);
-
-  LatticeFermion prom(FGrid);
-  prom=Zero();
-  for(int b=0;b<nbasis*2;b++){
-    prom=prom+CombinedUV.subspace[b];
-  }
-
-  std::cout<<GridLogMessage<<"c_src "<<norm2(c_src)<<std::endl;
-  std::cout<<GridLogMessage<<"prom  "<<norm2(prom)<<std::endl;
-
-  PVdagM.Op(prom,tmp);
-  blockProject(c_proj,tmp,CombinedUV.subspace);
-  std::cout<<GridLogMessage<<" Called Big Dirac Op "<<norm2(tmp)<<std::endl;
-
-  LittleDiracOpPV.M(c_src,c_res);
-  std::cout<<GridLogMessage<<" Called Little Dirac Op c_src "<< norm2(c_src) << "  c_res "<< norm2(c_res) <<std::endl;
-
-  std::cout<<GridLogMessage<<"Little dop : "<<norm2(c_res)<<std::endl;
-
-  std::cout<<GridLogMessage<<"Big dop in subspace : "<<norm2(c_proj)<<std::endl;
-
-  c_proj = c_proj - c_res;
-  std::cout<<GridLogMessage<<" ldop error: "<<norm2(c_proj)<<std::endl;
-
-
-  /**********
-   * Some solvers
-   **********
-   */
-
-  ///////////////////////////////////////
-  // Coarse grid solver test
-  ///////////////////////////////////////
-
-  std::cout<<GridLogMessage<<"******************* "<<std::endl;
-  std::cout<<GridLogMessage<<" Coarse Grid Solve -- Level 3 "<<std::endl;
-  std::cout<<GridLogMessage<<"******************* "<<std::endl;
-  TrivialPrecon<CoarseVector> simple;
-  NonHermitianLinearOperator<LittleDiracOperator,CoarseVector> LinOpCoarse(LittleDiracOpPV);
-  //  PrecGeneralisedConjugateResidualNonHermitian<CoarseVector>  L2PGCR(1.0e-4, 100, LinOpCoarse,simple,10,10); 
-  PrecGeneralisedConjugateResidualNonHermitian<CoarseVector>  L2PGCR(1.0e-2, 10, LinOpCoarse,simple,20,20); 
-  L2PGCR.Level(3);
-  c_res=Zero();
-  L2PGCR(c_src,c_res);
-
-  ////////////////////////////////////////
-  // Fine grid smoother
-  ////////////////////////////////////////
-  std::cout<<GridLogMessage<<"******************* "<<std::endl;
-  std::cout<<GridLogMessage<<" Fine Grid Smoother -- Level 2 "<<std::endl;
-  std::cout<<GridLogMessage<<"******************* "<<std::endl;
-  TrivialPrecon<LatticeFermionD> simple_fine;
-  //  NonHermitianLinearOperator<PVdagM_t,LatticeFermionD> LinOpSmooth(PVdagM);
-  PrecGeneralisedConjugateResidualNonHermitian<LatticeFermionD> SmootherGCR(0.01,1,ShiftedPVdagM,simple_fine,16,16);
-  SmootherGCR.Level(2);
-  
-  LatticeFermionD f_src(FGrid);
-  LatticeFermionD f_res(FGrid);
-
-  f_src = one;  // 1 in every element for vector 1.
-  f_res=Zero();
-  SmootherGCR(f_src,f_res);
-
-  typedef MGPreconditionerSVD<vSpinColourVector,  vTComplex,nbasis*2> TwoLevelMG;
-
-  TwoLevelMG TwoLevelPrecon(CombinedUV,CombinedUV,
-			    PVdagM,
-			    simple_fine,
-			    SmootherGCR,
-			    LinOpCoarse,
-			    L2PGCR);
-  
-  PrecGeneralisedConjugateResidualNonHermitian<LatticeFermion> L1PGCR(1.0e-8,1000,PVdagM,TwoLevelPrecon,20,20);
-  L1PGCR.Level(1);
-
-  f_res=Zero();
-  L1PGCR(f_src,f_res);
-
-  std::cout<<GridLogMessage<<std::endl;
-  std::cout<<GridLogMessage<<"*******************************************"<<std::endl;
-  std::cout<<GridLogMessage<<std::endl;
-  std::cout<<GridLogMessage << "Done "<< std::endl;
-
-  Grid_finalize();
-  return 0;
-}

tests/debug/Test_general_coarse_pvdagm_svd_uv.cc

@@ -1,479 +0,0 @@
-/*************************************************************************************
-
-    Grid physics library, www.github.com/paboyle/Grid 
-
-    Source file: ./tests/Test_padded_cell.cc
-
-    Copyright (C) 2023
-
-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>
-#include <Grid/lattice/PaddedCell.h>
-#include <Grid/stencil/GeneralLocalStencil.h>
-
-#include <Grid/algorithms/iterative/PrecGeneralisedConjugateResidual.h>
-#include <Grid/algorithms/iterative/PrecGeneralisedConjugateResidualNonHermitian.h>
-#include <Grid/algorithms/iterative/BiCGSTAB.h>
-
-using namespace std;
-using namespace Grid;
-
-template<class Matrix,class Field>
-class PVdagMLinearOperator : public LinearOperatorBase<Field> {
-  Matrix &_Mat;
-  Matrix &_PV;
-public:
-  PVdagMLinearOperator(Matrix &Mat,Matrix &PV): _Mat(Mat),_PV(PV){};
-
-  void OpDiag (const Field &in, Field &out) {    assert(0);  }
-  void OpDir  (const Field &in, Field &out,int dir,int disp) {    assert(0);  }
-  void OpDirAll  (const Field &in, std::vector<Field> &out){    assert(0);  };
-  void Op     (const Field &in, Field &out){
-    //    std::cout << GridLogMessage<< "Op: PVdag M "<<std::endl;
-    Field tmp(in.Grid());
-    _Mat.M(in,tmp);
-    _PV.Mdag(tmp,out);
-  }
-  void AdjOp     (const Field &in, Field &out){
-    //    std::cout << GridLogMessage<<"AdjOp: Mdag PV "<<std::endl;
-    Field tmp(in.Grid());
-    _PV.M(in,tmp);
-    _Mat.Mdag(tmp,out);
-  }
-  void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){
-    assert(0);
-  }
-  void HermOp(const Field &in, Field &out){
-    //    std::cout <<GridLogMessage<< "HermOp: Mdag PV PVdag M"<<std::endl;
-    Field tmp(in.Grid());
-    Op(in,tmp);
-    AdjOp(tmp,out);
-    //    std::cout << "HermOp done "<<norm2(out)<<std::endl;
-  }
-};
-template<class Matrix,class Field>
-class MdagPVLinearOperator : public LinearOperatorBase<Field> {
-  Matrix &_Mat;
-  Matrix &_PV;
-public:
-  MdagPVLinearOperator(Matrix &Mat,Matrix &PV): _Mat(Mat),_PV(PV){};
-
-  void OpDiag (const Field &in, Field &out) {    assert(0);  }
-  void OpDir  (const Field &in, Field &out,int dir,int disp) {    assert(0);  }
-  void OpDirAll  (const Field &in, std::vector<Field> &out){    assert(0);  };
-  void Op     (const Field &in, Field &out){
-    Field tmp(in.Grid());
-    //    std::cout <<GridLogMessage<< "Op: PVdag M "<<std::endl;
-    _PV.M(in,tmp);
-    _Mat.Mdag(tmp,out);
-  }
-  void AdjOp     (const Field &in, Field &out){
-    //    std::cout <<GridLogMessage<< "AdjOp: Mdag PV "<<std::endl;
-    Field tmp(in.Grid());
-    _Mat.M(in,tmp);
-    _PV.Mdag(tmp,out);
-  }
-  void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){
-    assert(0);
-  }
-  void HermOp(const Field &in, Field &out){
-    //    std::cout << GridLogMessage<<"HermOp: PVdag M Mdag PV "<<std::endl;
-    Field tmp(in.Grid());
-    Op(in,tmp);
-    AdjOp(tmp,out);
-    //    std::cout << "HermOp done "<<norm2(out)<<std::endl;
-  }
-};
-template<class Matrix,class Field>
-class ShiftedPVdagMLinearOperator : public LinearOperatorBase<Field> {
-  Matrix &_Mat;
-  Matrix &_PV;
-  RealD shift;
-public:
-  ShiftedPVdagMLinearOperator(RealD _shift,Matrix &Mat,Matrix &PV): shift(_shift),_Mat(Mat),_PV(PV){};
-
-  void OpDiag (const Field &in, Field &out) {    assert(0);  }
-  void OpDir  (const Field &in, Field &out,int dir,int disp) {    assert(0);  }
-  void OpDirAll  (const Field &in, std::vector<Field> &out){    assert(0);  };
-  void Op     (const Field &in, Field &out){
-    //    std::cout << "Op: PVdag M "<<std::endl;
-    Field tmp(in.Grid());
-    _Mat.M(in,tmp);
-    _PV.Mdag(tmp,out);
-    out = out + shift * in;
-  }
-  void AdjOp     (const Field &in, Field &out){
-    //    std::cout << "AdjOp: Mdag PV "<<std::endl;
-    Field tmp(in.Grid());
-    _PV.M(tmp,out);
-    _Mat.Mdag(in,tmp);
-    out = out + shift * in;
-  }
-  void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){    assert(0);  }
-  void HermOp(const Field &in, Field &out){
-    //    std::cout << "HermOp: Mdag PV PVdag M"<<std::endl;
-    Field tmp(in.Grid());
-    Op(in,tmp);
-    AdjOp(tmp,out);
-  }
-};
-template<class Fobj,class CComplex,int nbasis>
-class MGPreconditionerSVD : public LinearFunction< Lattice<Fobj> > {
-public:
-  using LinearFunction<Lattice<Fobj> >::operator();
-
-  typedef Aggregation<Fobj,CComplex,nbasis> Aggregates;
-  typedef typename Aggregation<Fobj,CComplex,nbasis>::FineField    FineField;
-  typedef typename Aggregation<Fobj,CComplex,nbasis>::CoarseVector CoarseVector;
-  typedef typename Aggregation<Fobj,CComplex,nbasis>::CoarseMatrix CoarseMatrix;
-  typedef LinearOperatorBase<FineField>                            FineOperator;
-  typedef LinearFunction    <FineField>                            FineSmoother;
-  typedef LinearOperatorBase<CoarseVector>                         CoarseOperator;
-  typedef LinearFunction    <CoarseVector>                         CoarseSolver;
-  ///////////////////////////////
-  // SVD is M = U S Vdag
-  //
-  // Define a subset of Vc and Uc in Complex_f,c  matrix
-  // - these are the coarsening, non-square matrices
-  //
-  // Solve a coarse approx to
-  //
-  //      M psi = eta
-  //
-  //  via
-  //
-  //  Uc^dag U S Vdag Vc Vc^dag psi = Uc^dag eta
-  //
-  //  M_coarse Vc^dag psi = M_coarse psi_c = eta_c
-  //  
-  ///////////////////////////////
-  Aggregates     & _U;
-  Aggregates     & _V;
-  FineOperator   & _FineOperator;
-  FineSmoother   & _PreSmoother;
-  FineSmoother   & _PostSmoother;
-  CoarseOperator & _CoarseOperator;
-  CoarseSolver   & _CoarseSolve;
-
-  int    level;  void Level(int lv) {level = lv; };
-
-  MGPreconditionerSVD(Aggregates &U,
-		      Aggregates &V,
-		      FineOperator &Fine,
-		      FineSmoother &PreSmoother,
-		      FineSmoother &PostSmoother,
-		      CoarseOperator &CoarseOperator_,
-		      CoarseSolver &CoarseSolve_)
-    : _U(U),
-      _V(V),
-      _FineOperator(Fine),
-      _PreSmoother(PreSmoother),
-      _PostSmoother(PostSmoother),
-      _CoarseOperator(CoarseOperator_),
-      _CoarseSolve(CoarseSolve_),
-      level(1)  {  }
-
-  virtual void operator()(const FineField &in, FineField & out) 
-  {
-    GridBase *CoarseGrid = _U.CoarseGrid;
-    //    auto CoarseGrid = _CoarseOperator.Grid();
-    CoarseVector Csrc(CoarseGrid);
-    CoarseVector Csol(CoarseGrid);
-    FineField vec1(in.Grid());
-    FineField vec2(in.Grid());
-
-    std::cout<<GridLogMessage << "Calling PreSmoother " <<std::endl;
-
-    //    std::cout<<GridLogMessage << "Calling PreSmoother input residual "<<norm2(in) <<std::endl;
-    double t;
-    // Fine Smoother
-    //    out = in;
-    out = Zero();
-    t=-usecond();
-    _PreSmoother(in,out);
-    t+=usecond();
-
-    std::cout<<GridLogMessage << "PreSmoother took "<< t/1000.0<< "ms" <<std::endl;
-
-    // Update the residual
-    _FineOperator.Op(out,vec1);  sub(vec1, in ,vec1);   
-    //    std::cout<<GridLogMessage <<"Residual-1 now " <<norm2(vec1)<<std::endl;
-
-    //  Uc^dag U S Vdag Vc Vc^dag psi = Uc^dag eta
-    // Fine to Coarse 
-    t=-usecond();
-    _U.ProjectToSubspace  (Csrc,vec1);
-    t+=usecond();
-    std::cout<<GridLogMessage << "Project to coarse took "<< t/1000.0<< "ms" <<std::endl;
-
-    // Coarse correction
-    t=-usecond();
-    Csol = Zero();
-    _CoarseSolve(Csrc,Csol);
-    //Csol=Zero();
-    t+=usecond();
-    std::cout<<GridLogMessage << "Coarse solve took "<< t/1000.0<< "ms" <<std::endl;
-
-    // Coarse to Fine
-    t=-usecond();  
-    //    _CoarseOperator.PromoteFromSubspace(_Aggregates,Csol,vec1);
-    _V.PromoteFromSubspace(Csol,vec1); 
-    add(out,out,vec1);
-    t+=usecond();
-    std::cout<<GridLogMessage << "Promote to this level took "<< t/1000.0<< "ms" <<std::endl;
-
-    // Residual
-    _FineOperator.Op(out,vec1);  sub(vec1 ,in , vec1);  
-    //    std::cout<<GridLogMessage <<"Residual-2 now " <<norm2(vec1)<<std::endl;
-
-    // Fine Smoother
-    t=-usecond();
-    //    vec2=vec1;
-    vec2=Zero();
-    _PostSmoother(vec1,vec2);
-    t+=usecond();
-    std::cout<<GridLogMessage << "PostSmoother took "<< t/1000.0<< "ms" <<std::endl;
-
-    add( out,out,vec2);
-    std::cout<<GridLogMessage << "Done " <<std::endl;
-  }
-};
-
-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
-  Coordinate clatt = GridDefaultLatt();
-  for(int d=0;d<clatt.size();d++){
-    clatt[d] = clatt[d]/2;
-    //    clatt[d] = clatt[d]/4;
-  }
-  GridCartesian *Coarse4d =  SpaceTimeGrid::makeFourDimGrid(clatt, GridDefaultSimd(Nd,vComplex::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); random(RNG5,src);
-  LatticeFermion result(FGrid); result=Zero();
-  LatticeFermion    ref(FGrid); ref=Zero();
-  LatticeFermion    tmp(FGrid);
-  LatticeFermion    err(FGrid);
-  LatticeGaugeField Umu(UGrid);
-
-  FieldMetaData header;
-  std::string file("ckpoint_lat.4000");
-  NerscIO::readConfiguration(Umu,header,file);
-  
-  RealD mass=0.01;
-  RealD M5=1.8;
-
-  DomainWallFermionD Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5);
-  DomainWallFermionD Dpv(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,1.0,M5);
-
-  const int nbasis = 60;
-  const int cb = 0 ;
-
-  NextToNearestStencilGeometry5D geom(Coarse5d);
-
-  std::cout<<GridLogMessage<<std::endl;
-  std::cout<<GridLogMessage<<"*******************************************"<<std::endl;
-  std::cout<<GridLogMessage<<std::endl;
-
-  typedef PVdagMLinearOperator<DomainWallFermionD,LatticeFermionD> PVdagM_t;
-  typedef MdagPVLinearOperator<DomainWallFermionD,LatticeFermionD> MdagPV_t;
-  typedef ShiftedPVdagMLinearOperator<DomainWallFermionD,LatticeFermionD> ShiftedPVdagM_t;
-  PVdagM_t PVdagM(Ddwf,Dpv);
-  MdagPV_t MdagPV(Ddwf,Dpv);
-  //  ShiftedPVdagM_t ShiftedPVdagM(2.0,Ddwf,Dpv); // 355
-  //  ShiftedPVdagM_t ShiftedPVdagM(1.0,Ddwf,Dpv); // 246
-  //  ShiftedPVdagM_t ShiftedPVdagM(0.5,Ddwf,Dpv); // 183
-  //  ShiftedPVdagM_t ShiftedPVdagM(0.25,Ddwf,Dpv); // 145
-  //  ShiftedPVdagM_t ShiftedPVdagM(0.1,Ddwf,Dpv); // 134
-  //  ShiftedPVdagM_t ShiftedPVdagM(0.1,Ddwf,Dpv); // 127 -- NULL space via inverse iteration
-  //  ShiftedPVdagM_t ShiftedPVdagM(0.1,Ddwf,Dpv); // 57 -- NULL space via inverse iteration; 3 iterations
-  //  ShiftedPVdagM_t ShiftedPVdagM(0.25,Ddwf,Dpv); // 57 , tighter inversion
-  //  ShiftedPVdagM_t ShiftedPVdagM(0.25,Ddwf,Dpv); // nbasis 20 -- 49 iters
-  //  ShiftedPVdagM_t ShiftedPVdagM(0.25,Ddwf,Dpv); // nbasis 20 -- 70 iters; asymmetric 
-  //  ShiftedPVdagM_t ShiftedPVdagM(0.25,Ddwf,Dpv); // 58; Loosen coarse, tighten fine
-  //  ShiftedPVdagM_t ShiftedPVdagM(0.1,Ddwf,Dpv); // 56 ... 
-  //  ShiftedPVdagM_t ShiftedPVdagM(0.1,Ddwf,Dpv); // 51 ...  with 24 vecs
-  //  ShiftedPVdagM_t ShiftedPVdagM(0.1,Ddwf,Dpv); // 31 ...  with 24 vecs and 2^4 blocking
-  //  ShiftedPVdagM_t ShiftedPVdagM(0.1,Ddwf,Dpv); // 43 ...  with 16 vecs and 2^4 blocking, sloppier
-  //  ShiftedPVdagM_t ShiftedPVdagM(0.1,Ddwf,Dpv); // 35  ...  with 20 vecs and 2^4 blocking
-  //  ShiftedPVdagM_t ShiftedPVdagM(0.1,Ddwf,Dpv); // 35  ...  with 20 vecs and 2^4 blocking, looser coarse
-  //  ShiftedPVdagM_t ShiftedPVdagM(0.1,Ddwf,Dpv); // 64  ...  with 20 vecs, Christoph setup, and 2^4 blocking, looser coarse
-  ShiftedPVdagM_t ShiftedPVdagM(0.01,Ddwf,Dpv); // 
-
-
-  // Run power method on HOA??
-  PowerMethod<LatticeFermion>       PM;
-  PM(PVdagM,src);
-  PM(MdagPV,src);
-
- 
-  // Warning: This routine calls PVdagM.Op, not PVdagM.HermOp
-  typedef Aggregation<vSpinColourVector,vTComplex,nbasis> Subspace;
-  Subspace V(Coarse5d,FGrid,cb);
-  //  Subspace U(Coarse5d,FGrid,cb);
-
-  // Breeds right singular vectors with call to HermOp
-  V.CreateSubspaceChebyshev(RNG5,PVdagM,
-			    nbasis,
-			    4000.0,0.003,
-			    300);
-
-  // Breeds left singular vectors with call to HermOp
-  //  U.CreateSubspaceChebyshev(RNG5,MdagPV,
-  //			    nbasis,
-  //			    4000.0,0.003,
-  //			    300);
-  //  U.subspace=V.subspace;
-  
-  //  typedef Aggregation<vSpinColourVector,vTComplex,2*nbasis> CombinedSubspace;
-  //  CombinedSubspace CombinedUV(Coarse5d,FGrid,cb);
-  //  for(int b=0;b<nbasis;b++){
-  //    CombinedUV.subspace[b]        = V.subspace[b];
-  //    CombinedUV.subspace[b+nbasis] = U.subspace[b];
-  //  }
-  
-  
-  //  typedef GeneralCoarsenedMatrix<vSpinColourVector,vTComplex,2*nbasis> LittleDiracOperator;
-  typedef GeneralCoarsenedMatrix<vSpinColourVector,vTComplex,nbasis> LittleDiracOperator;
-  typedef LittleDiracOperator::CoarseVector CoarseVector;
-  LittleDiracOperator LittleDiracOpPV(geom,FGrid,Coarse5d);
-  LittleDiracOpPV.CoarsenOperator(PVdagM,V,V);
-
-  std::cout<<GridLogMessage<<std::endl;
-  std::cout<<GridLogMessage<<"*******************************************"<<std::endl;
-  std::cout<<GridLogMessage<<std::endl;
-  std::cout<<GridLogMessage<<"Testing coarsened operator "<<std::endl;
-  
-  CoarseVector c_src (Coarse5d);
-  CoarseVector c_res (Coarse5d);
-  CoarseVector c_proj(Coarse5d);
-
-
-  Complex one(1.0);
-  c_src = one;  // 1 in every element for vector 1.
-
-  //  blockPromote(c_src,err,CoarseToFine.subspace);
-
-  LatticeFermion prom(FGrid);
-  prom=Zero();
-  for(int b=0;b<nbasis;b++){
-    prom=prom+V.subspace[b];
-  }
-
-  std::cout<<GridLogMessage<<"c_src "<<norm2(c_src)<<std::endl;
-  std::cout<<GridLogMessage<<"prom  "<<norm2(prom)<<std::endl;
-
-  PVdagM.Op(prom,tmp);
-  blockProject(c_proj,tmp,V.subspace);
-  std::cout<<GridLogMessage<<" Called Big Dirac Op "<<norm2(tmp)<<std::endl;
-
-  LittleDiracOpPV.M(c_src,c_res);
-  std::cout<<GridLogMessage<<" Called Little Dirac Op c_src "<< norm2(c_src) << "  c_res "<< norm2(c_res) <<std::endl;
-
-  std::cout<<GridLogMessage<<"Little dop : "<<norm2(c_res)<<std::endl;
-
-  std::cout<<GridLogMessage<<"Big dop in subspace : "<<norm2(c_proj)<<std::endl;
-
-  c_proj = c_proj - c_res;
-  std::cout<<GridLogMessage<<" ldop error: "<<norm2(c_proj)<<std::endl;
-
-
-  /**********
-   * Some solvers
-   **********
-   */
-
-  ///////////////////////////////////////
-  // Coarse grid solver test
-  ///////////////////////////////////////
-
-  std::cout<<GridLogMessage<<"******************* "<<std::endl;
-  std::cout<<GridLogMessage<<" Coarse Grid Solve -- Level 3 "<<std::endl;
-  std::cout<<GridLogMessage<<"******************* "<<std::endl;
-  TrivialPrecon<CoarseVector> simple;
-  NonHermitianLinearOperator<LittleDiracOperator,CoarseVector> LinOpCoarse(LittleDiracOpPV);
-  //  PrecGeneralisedConjugateResidualNonHermitian<CoarseVector>  L2PGCR(1.0e-4, 100, LinOpCoarse,simple,10,10); 
-  PrecGeneralisedConjugateResidualNonHermitian<CoarseVector>  L3PGCR(1.0e-4, 10, LinOpCoarse,simple,20,20); 
-  L3PGCR.Level(3);
-  c_res=Zero();
-  L3PGCR(c_src,c_res);
-
-  ////////////////////////////////////////
-  // Fine grid smoother
-  ////////////////////////////////////////
-  std::cout<<GridLogMessage<<"******************* "<<std::endl;
-  std::cout<<GridLogMessage<<" Fine Grid Smoother -- Level 2 "<<std::endl;
-  std::cout<<GridLogMessage<<"******************* "<<std::endl;
-  TrivialPrecon<LatticeFermionD> simple_fine;
-  //  NonHermitianLinearOperator<PVdagM_t,LatticeFermionD> LinOpSmooth(PVdagM);
-  PrecGeneralisedConjugateResidualNonHermitian<LatticeFermionD> SmootherGCR(0.01,1,ShiftedPVdagM,simple_fine,16,16);
-  SmootherGCR.Level(2);
-  
-  LatticeFermionD f_src(FGrid);
-  LatticeFermionD f_res(FGrid);
-
-  f_src = one;  // 1 in every element for vector 1.
-  f_res=Zero();
-  SmootherGCR(f_src,f_res);
-
-  //  typedef MGPreconditionerSVD<vSpinColourVector,  vTComplex,nbasis*2> TwoLevelMG;
-  typedef MGPreconditionerSVD<vSpinColourVector,  vTComplex,nbasis> TwoLevelMG;
-
-  //  TwoLevelMG TwoLevelPrecon(CombinedUV,CombinedUV,
-  TwoLevelMG TwoLevelPrecon(V,V,
-			    PVdagM,
-			    simple_fine,
-			    SmootherGCR,
-			    LinOpCoarse,
-			    L3PGCR);
-  
-  PrecGeneralisedConjugateResidualNonHermitian<LatticeFermion> L1PGCR(1.0e-8,1000,PVdagM,TwoLevelPrecon,16,16);
-  L1PGCR.Level(1);
-
-  f_res=Zero();
-  L1PGCR(f_src,f_res);
-
-  std::cout<<GridLogMessage<<std::endl;
-  std::cout<<GridLogMessage<<"*******************************************"<<std::endl;
-  std::cout<<GridLogMessage<<std::endl;
-  std::cout<<GridLogMessage << "Done "<< std::endl;
-
-  Grid_finalize();
-  return 0;
-}

tests/debug/Test_general_coarse_wilson.cc

@@ -1,333 +0,0 @@
-/*************************************************************************************
-
-    Grid physics library, www.github.com/paboyle/Grid 
-
-    Source file: ./tests/Test_padded_cell.cc
-
-    Copyright (C) 2023
-
-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>
-#include <Grid/lattice/PaddedCell.h>
-#include <Grid/stencil/GeneralLocalStencil.h>
-
-#include <Grid/algorithms/iterative/PrecGeneralisedConjugateResidual.h>
-#include <Grid/algorithms/iterative/PrecGeneralisedConjugateResidualNonHermitian.h>
-#include <Grid/algorithms/iterative/BiCGSTAB.h>
-
-using namespace std;
-using namespace Grid;
-
-template<class Fobj,class CComplex,int nbasis>
-class MGPreconditioner : public LinearFunction< Lattice<Fobj> > {
-public:
-  using LinearFunction<Lattice<Fobj> >::operator();
-
-  typedef Aggregation<Fobj,CComplex,nbasis> Aggregates;
-  typedef typename Aggregation<Fobj,CComplex,nbasis>::FineField    FineField;
-  typedef typename Aggregation<Fobj,CComplex,nbasis>::CoarseVector CoarseVector;
-  typedef typename Aggregation<Fobj,CComplex,nbasis>::CoarseMatrix CoarseMatrix;
-  typedef LinearOperatorBase<FineField>                            FineOperator;
-  typedef LinearFunction    <FineField>                            FineSmoother;
-  typedef LinearOperatorBase<CoarseVector>                         CoarseOperator;
-  typedef LinearFunction    <CoarseVector>                         CoarseSolver;
-  Aggregates     & _Aggregates;
-  FineOperator   & _FineOperator;
-  FineSmoother   & _PreSmoother;
-  FineSmoother   & _PostSmoother;
-  CoarseOperator & _CoarseOperator;
-  CoarseSolver   & _CoarseSolve;
-
-  int    level;  void Level(int lv) {level = lv; };
-
-  MGPreconditioner(Aggregates &Agg,
-		   FineOperator &Fine,
-		   FineSmoother &PreSmoother,
-		   FineSmoother &PostSmoother,
-		   CoarseOperator &CoarseOperator_,
-		   CoarseSolver &CoarseSolve_)
-    : _Aggregates(Agg),
-      _FineOperator(Fine),
-      _PreSmoother(PreSmoother),
-      _PostSmoother(PostSmoother),
-      _CoarseOperator(CoarseOperator_),
-      _CoarseSolve(CoarseSolve_),
-      level(1)  {  }
-
-  virtual void operator()(const FineField &in, FineField & out) 
-  {
-    GridBase *CoarseGrid = _Aggregates.CoarseGrid;
-    //    auto CoarseGrid = _CoarseOperator.Grid();
-    CoarseVector Csrc(CoarseGrid);
-    CoarseVector Csol(CoarseGrid);
-    FineField vec1(in.Grid());
-    FineField vec2(in.Grid());
-
-    std::cout<<GridLogMessage << "Calling PreSmoother " <<std::endl;
-
-    //    std::cout<<GridLogMessage << "Calling PreSmoother input residual "<<norm2(in) <<std::endl;
-    double t;
-    // Fine Smoother
-    //    out = in;
-    out = Zero();
-    t=-usecond();
-    _PreSmoother(in,out);
-    t+=usecond();
-
-    std::cout<<GridLogMessage << "PreSmoother took "<< t/1000.0<< "ms" <<std::endl;
-
-    // Update the residual
-    _FineOperator.Op(out,vec1);  sub(vec1, in ,vec1);   
-    //    std::cout<<GridLogMessage <<"Residual-1 now " <<norm2(vec1)<<std::endl;
-
-    // Fine to Coarse 
-    t=-usecond();
-    _Aggregates.ProjectToSubspace  (Csrc,vec1);
-    t+=usecond();
-    std::cout<<GridLogMessage << "Project to coarse took "<< t/1000.0<< "ms" <<std::endl;
-
-    // Coarse correction
-    t=-usecond();
-    Csol = Zero();
-    _CoarseSolve(Csrc,Csol);
-    //Csol=Zero();
-    t+=usecond();
-    std::cout<<GridLogMessage << "Coarse solve took "<< t/1000.0<< "ms" <<std::endl;
-
-    // Coarse to Fine
-    t=-usecond();  
-    //    _CoarseOperator.PromoteFromSubspace(_Aggregates,Csol,vec1);
-    _Aggregates.PromoteFromSubspace(Csol,vec1); 
-    add(out,out,vec1);
-    t+=usecond();
-    std::cout<<GridLogMessage << "Promote to this level took "<< t/1000.0<< "ms" <<std::endl;
-
-    // Residual
-    _FineOperator.Op(out,vec1);  sub(vec1 ,in , vec1);  
-    //    std::cout<<GridLogMessage <<"Residual-2 now " <<norm2(vec1)<<std::endl;
-
-    // Fine Smoother
-    t=-usecond();
-    //    vec2=vec1;
-    vec2=Zero();
-    _PostSmoother(vec1,vec2);
-    t+=usecond();
-    std::cout<<GridLogMessage << "PostSmoother took "<< t/1000.0<< "ms" <<std::endl;
-
-    add( out,out,vec2);
-    std::cout<<GridLogMessage << "Done " <<std::endl;
-  }
-};
-
-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   = UGrid;
-  GridRedBlackCartesian * FrbGrid = UrbGrid;
-
-  // Construct a coarsened grid
-  Coordinate clatt = GridDefaultLatt();
-  for(int d=0;d<clatt.size();d++){
-    clatt[d] = clatt[d]/2;
-    //clatt[d] = clatt[d]/4;
-  }
-  GridCartesian *Coarse4d =  SpaceTimeGrid::makeFourDimGrid(clatt, GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi());;
-
-  std::vector<int> seeds4({1,2,3,4});
-  std::vector<int> cseeds({5,6,7,8});
-  GridParallelRNG          RNG4(UGrid);   RNG4.SeedFixedIntegers(seeds4);
-  GridParallelRNG          CRNG(Coarse4d);CRNG.SeedFixedIntegers(cseeds);
-
-  Complex one(1.0);
-
-  LatticeFermion    src(FGrid); src=one;
-  LatticeFermion result(FGrid); result=Zero();
-  LatticeFermion    ref(FGrid); ref=Zero();
-  LatticeFermion    tmp(FGrid);
-  LatticeFermion    err(FGrid);
-  LatticeFermion    precsrc(FGrid);
-  LatticeGaugeField Umu(UGrid);
-
-  FieldMetaData header;
-  std::string file("ckpoint_lat");
-  NerscIO::readConfiguration(Umu,header,file);
-
-  RealD csw =0.0;
-  RealD mass=-0.92;
-
-  WilsonCloverFermionD Dw(Umu,*UGrid,*UrbGrid,mass,csw,csw);
-
-  const int nbasis = 20;
-  const int cb = 0 ;
-  LatticeFermion prom(FGrid);
-
-  typedef GeneralCoarsenedMatrix<vSpinColourVector,vTComplex,2*nbasis> LittleDiracOperator;
-  typedef LittleDiracOperator::CoarseVector CoarseVector;
-
-  NearestStencilGeometry4D geom(Coarse4d);
-
-  std::cout<<GridLogMessage<<std::endl;
-  std::cout<<GridLogMessage<<"*******************************************"<<std::endl;
-  std::cout<<GridLogMessage<<std::endl;
- 
-  // Warning: This routine calls Linop.Op, not LinOpo.HermOp
-  typedef Aggregation<vSpinColourVector,vTComplex,nbasis> Subspace;
-  Subspace Aggregates(Coarse4d,FGrid,cb);
-
-  NonHermitianLinearOperator<WilsonCloverFermionD,LatticeFermion> LinOpDw(Dw);
-  ShiftedNonHermitianLinearOperator<WilsonCloverFermionD,LatticeFermion> ShiftedLinOpDw(Dw,0.01);
-
-  Aggregates.CreateSubspaceGCR(RNG4,
-			       LinOpDw,
-			       nbasis);
-  
-  typedef Aggregation<vSpinColourVector,vTComplex,2*nbasis> CombinedSubspace;
-  CombinedSubspace CombinedUV(Coarse4d,UGrid,cb);
-  for(int b=0;b<nbasis;b++){
-    Gamma G5(Gamma::Algebra::Gamma5);
-    CombinedUV.subspace[b]        = Aggregates.subspace[b];
-    CombinedUV.subspace[b+nbasis] = G5*Aggregates.subspace[b];
-  }
-
-  LittleDiracOperator LittleDiracOp(geom,FGrid,Coarse4d);
-  LittleDiracOp.CoarsenOperator(LinOpDw,CombinedUV);
-
-  std::cout<<GridLogMessage<<std::endl;
-  std::cout<<GridLogMessage<<"*******************************************"<<std::endl;
-  std::cout<<GridLogMessage<<std::endl;
-  std::cout<<GridLogMessage<<"Testing coarsened operator "<<std::endl;
-  
-  CoarseVector c_src (Coarse4d);
-  CoarseVector c_res (Coarse4d);
-  CoarseVector c_proj(Coarse4d);
-
-  std::vector<LatticeFermion> subspace(2*nbasis,FGrid);
-  subspace=CombinedUV.subspace;
-
-  c_src = one;  // 1 in every element for vector 1.
-  blockPromote(c_src,err,subspace);
-
-  prom=Zero();
-  for(int b=0;b<2*nbasis;b++){
-    prom=prom+subspace[b];
-  }
-  err=err-prom; 
-  std::cout<<GridLogMessage<<"Promoted back from subspace: err "<<norm2(err)<<std::endl;
-  std::cout<<GridLogMessage<<"c_src "<<norm2(c_src)<<std::endl;
-  std::cout<<GridLogMessage<<"prom  "<<norm2(prom)<<std::endl;
-
-  LinOpDw.Op(prom,tmp);
-  blockProject(c_proj,tmp,subspace);
-  std::cout<<GridLogMessage<<" Called Big Dirac Op "<<norm2(tmp)<<std::endl;
-
-  LittleDiracOp.M(c_src,c_res);
-  std::cout<<GridLogMessage<<" Called Little Dirac Op c_src "<< norm2(c_src) << "  c_res "<< norm2(c_res) <<std::endl;
-
-  std::cout<<GridLogMessage<<"Little dop : "<<norm2(c_res)<<std::endl;
-  //  std::cout<<GridLogMessage<<" Little "<< c_res<<std::endl;
-
-  std::cout<<GridLogMessage<<"Big dop in subspace : "<<norm2(c_proj)<<std::endl;
-  //  std::cout<<GridLogMessage<<" Big "<< c_proj<<std::endl;
-  c_proj = c_proj - c_res;
-  std::cout<<GridLogMessage<<" ldop error: "<<norm2(c_proj)<<std::endl;
-  //  std::cout<<GridLogMessage<<" error "<< c_proj<<std::endl;
-
-
-  /**********
-   * Some solvers
-   **********
-   */
-
-  // CG
-  {
-    MdagMLinearOperator<WilsonFermionD,LatticeFermion> HermOp(Dw);
-    ConjugateGradient<LatticeFermion> CG(1.0e-8,10000);
-    Dw.Mdag(src,precsrc);
-    CG(HermOp,precsrc,result);
-    result=Zero();
-  }
- 
-  ///////////////////////////////////////
-  // Coarse grid solver test
-  ///////////////////////////////////////
-
-  std::cout<<GridLogMessage<<"******************* "<<std::endl;
-  std::cout<<GridLogMessage<<" Coarse Grid Solve -- Level 3 "<<std::endl;
-  std::cout<<GridLogMessage<<"******************* "<<std::endl;
-  TrivialPrecon<CoarseVector> simple;
-  NonHermitianLinearOperator<LittleDiracOperator,CoarseVector> LinOpCoarse(LittleDiracOp);
-  ShiftedNonHermitianLinearOperator<LittleDiracOperator,CoarseVector> ShiftedLinOpCoarse(LittleDiracOp,0.001);
-  //  ShiftedNonHermitianLinearOperator<LittleDiracOperator,CoarseVector> ShiftedLinOpCoarse(LittleDiracOp,0.01);
-  //  ShiftedNonHermitianLinearOperator<LittleDiracOperator,CoarseVector> ShiftedLinOpCoarse(LinOpCoarse,0.001);
-  //  PrecGeneralisedConjugateResidualNonHermitian<CoarseVector>  L2PGCR(1.0e-4, 100, LinOpCoarse,simple,10,10); 
-  //  PrecGeneralisedConjugateResidualNonHermitian<CoarseVector>  L2PGCR(1.0e-1, 100, LinOpCoarse,simple,30,30); 
-  PrecGeneralisedConjugateResidualNonHermitian<CoarseVector>  L2PGCR(2.0e-1, 50, ShiftedLinOpCoarse,simple,50,50); 
-  L2PGCR.Level(3);
-  c_res=Zero();
-  L2PGCR(c_src,c_res);
-
-  ////////////////////////////////////////
-  // Fine grid smoother
-  ////////////////////////////////////////
-  std::cout<<GridLogMessage<<"******************* "<<std::endl;
-  std::cout<<GridLogMessage<<" Fine Grid Smoother -- Level 2 "<<std::endl;
-  std::cout<<GridLogMessage<<"******************* "<<std::endl;
-  TrivialPrecon<LatticeFermionD> simple_fine;
-
-  PrecGeneralisedConjugateResidualNonHermitian<LatticeFermionD> SmootherGCR(0.1,1,ShiftedLinOpDw,simple_fine,4,4);
-  SmootherGCR.Level(2);
-  
-  LatticeFermionD f_src(FGrid);
-  LatticeFermionD f_res(FGrid);
-
-  f_src = one;  // 1 in every element for vector 1.
-  f_res=Zero();
-  SmootherGCR(f_src,f_res);
-
-  typedef MGPreconditioner<vSpinColourVector,  vTComplex,2*nbasis> TwoLevelMG;
-
-  TwoLevelMG TwoLevelPrecon(CombinedUV,
-			    LinOpDw,
-			    simple_fine,
-			    SmootherGCR,
-			    LinOpCoarse,
-			    L2PGCR);
-  
-  PrecGeneralisedConjugateResidualNonHermitian<LatticeFermion> L1PGCR(1.0e-8,1000,LinOpDw,TwoLevelPrecon,16,16);
-  L1PGCR.Level(1);
-
-  f_res=Zero();
-  L1PGCR(f_src,f_res);
-
-  std::cout<<GridLogMessage<<std::endl;
-  std::cout<<GridLogMessage<<"*******************************************"<<std::endl;
-  std::cout<<GridLogMessage<<std::endl;
-  std::cout<<GridLogMessage << "Done "<< std::endl;
-
-  Grid_finalize();
-  return 0;
-}

tests/debug/Test_general_coarse_wilson_nog5.cc

@@ -1,326 +0,0 @@
-/*************************************************************************************
-
-    Grid physics library, www.github.com/paboyle/Grid 
-
-    Source file: ./tests/Test_padded_cell.cc
-
-    Copyright (C) 2023
-
-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>
-#include <Grid/lattice/PaddedCell.h>
-#include <Grid/stencil/GeneralLocalStencil.h>
-
-#include <Grid/algorithms/iterative/PrecGeneralisedConjugateResidual.h>
-#include <Grid/algorithms/iterative/PrecGeneralisedConjugateResidualNonHermitian.h>
-#include <Grid/algorithms/iterative/BiCGSTAB.h>
-
-using namespace std;
-using namespace Grid;
-
-template<class Fobj,class CComplex,int nbasis>
-class MGPreconditioner : public LinearFunction< Lattice<Fobj> > {
-public:
-  using LinearFunction<Lattice<Fobj> >::operator();
-
-  typedef Aggregation<Fobj,CComplex,nbasis> Aggregates;
-  typedef typename Aggregation<Fobj,CComplex,nbasis>::FineField    FineField;
-  typedef typename Aggregation<Fobj,CComplex,nbasis>::CoarseVector CoarseVector;
-  typedef typename Aggregation<Fobj,CComplex,nbasis>::CoarseMatrix CoarseMatrix;
-  typedef LinearOperatorBase<FineField>                            FineOperator;
-  typedef LinearFunction    <FineField>                            FineSmoother;
-  typedef LinearOperatorBase<CoarseVector>                         CoarseOperator;
-  typedef LinearFunction    <CoarseVector>                         CoarseSolver;
-  Aggregates     & _Aggregates;
-  FineOperator   & _FineOperator;
-  FineSmoother   & _PreSmoother;
-  FineSmoother   & _PostSmoother;
-  CoarseOperator & _CoarseOperator;
-  CoarseSolver   & _CoarseSolve;
-
-  int    level;  void Level(int lv) {level = lv; };
-
-  MGPreconditioner(Aggregates &Agg,
-		   FineOperator &Fine,
-		   FineSmoother &PreSmoother,
-		   FineSmoother &PostSmoother,
-		   CoarseOperator &CoarseOperator_,
-		   CoarseSolver &CoarseSolve_)
-    : _Aggregates(Agg),
-      _FineOperator(Fine),
-      _PreSmoother(PreSmoother),
-      _PostSmoother(PostSmoother),
-      _CoarseOperator(CoarseOperator_),
-      _CoarseSolve(CoarseSolve_),
-      level(1)  {  }
-
-  virtual void operator()(const FineField &in, FineField & out) 
-  {
-    GridBase *CoarseGrid = _Aggregates.CoarseGrid;
-    //    auto CoarseGrid = _CoarseOperator.Grid();
-    CoarseVector Csrc(CoarseGrid);
-    CoarseVector Csol(CoarseGrid);
-    FineField vec1(in.Grid());
-    FineField vec2(in.Grid());
-
-    std::cout<<GridLogMessage << "Calling PreSmoother " <<std::endl;
-
-    //    std::cout<<GridLogMessage << "Calling PreSmoother input residual "<<norm2(in) <<std::endl;
-    double t;
-    // Fine Smoother
-    //    out = in;
-    out = Zero();
-    t=-usecond();
-    _PreSmoother(in,out);
-    t+=usecond();
-
-    std::cout<<GridLogMessage << "PreSmoother took "<< t/1000.0<< "ms" <<std::endl;
-
-    // Update the residual
-    _FineOperator.Op(out,vec1);  sub(vec1, in ,vec1);   
-    //    std::cout<<GridLogMessage <<"Residual-1 now " <<norm2(vec1)<<std::endl;
-
-    // Fine to Coarse 
-    t=-usecond();
-    _Aggregates.ProjectToSubspace  (Csrc,vec1);
-    t+=usecond();
-    std::cout<<GridLogMessage << "Project to coarse took "<< t/1000.0<< "ms" <<std::endl;
-
-    // Coarse correction
-    t=-usecond();
-    Csol = Zero();
-    _CoarseSolve(Csrc,Csol);
-    //Csol=Zero();
-    t+=usecond();
-    std::cout<<GridLogMessage << "Coarse solve took "<< t/1000.0<< "ms" <<std::endl;
-
-    // Coarse to Fine
-    t=-usecond();  
-    //    _CoarseOperator.PromoteFromSubspace(_Aggregates,Csol,vec1);
-    _Aggregates.PromoteFromSubspace(Csol,vec1); 
-    add(out,out,vec1);
-    t+=usecond();
-    std::cout<<GridLogMessage << "Promote to this level took "<< t/1000.0<< "ms" <<std::endl;
-
-    // Residual
-    _FineOperator.Op(out,vec1);  sub(vec1 ,in , vec1);  
-    //    std::cout<<GridLogMessage <<"Residual-2 now " <<norm2(vec1)<<std::endl;
-
-    // Fine Smoother
-    t=-usecond();
-    //    vec2=vec1;
-    vec2=Zero();
-    _PostSmoother(vec1,vec2);
-    t+=usecond();
-    std::cout<<GridLogMessage << "PostSmoother took "<< t/1000.0<< "ms" <<std::endl;
-
-    add( out,out,vec2);
-    std::cout<<GridLogMessage << "Done " <<std::endl;
-  }
-};
-
-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   = UGrid;
-  GridRedBlackCartesian * FrbGrid = UrbGrid;
-
-  // Construct a coarsened grid
-  Coordinate clatt = GridDefaultLatt();
-  for(int d=0;d<clatt.size();d++){
-    clatt[d] = clatt[d]/2;
-    //    clatt[d] = clatt[d]/4;
-  }
-  GridCartesian *Coarse4d =  SpaceTimeGrid::makeFourDimGrid(clatt, GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi());;
-
-  std::vector<int> seeds4({1,2,3,4});
-  std::vector<int> cseeds({5,6,7,8});
-  GridParallelRNG          RNG4(UGrid);   RNG4.SeedFixedIntegers(seeds4);
-  GridParallelRNG          CRNG(Coarse4d);CRNG.SeedFixedIntegers(cseeds);
-
-  Complex one(1.0);
-
-  LatticeFermion    src(FGrid); src=one;
-  LatticeFermion result(FGrid); result=Zero();
-  LatticeFermion    ref(FGrid); ref=Zero();
-  LatticeFermion    tmp(FGrid);
-  LatticeFermion    err(FGrid);
-  LatticeFermion    precsrc(FGrid);
-  LatticeGaugeField Umu(UGrid);
-
-  FieldMetaData header;
-  std::string file("ckpoint_lat");
-  NerscIO::readConfiguration(Umu,header,file);
-
-  RealD csw =0.0;
-  RealD mass=-0.92;
-
-  WilsonCloverFermionD Dw(Umu,*UGrid,*UrbGrid,mass,csw,csw);
-
-  const int nbasis = 40;
-  const int cb = 0 ;
-  LatticeFermion prom(FGrid);
-
-  typedef GeneralCoarsenedMatrix<vSpinColourVector,vTComplex,nbasis> LittleDiracOperator;
-  typedef LittleDiracOperator::CoarseVector CoarseVector;
-
-  NearestStencilGeometry4D geom(Coarse4d);
-
-  std::cout<<GridLogMessage<<std::endl;
-  std::cout<<GridLogMessage<<"*******************************************"<<std::endl;
-  std::cout<<GridLogMessage<<std::endl;
- 
-  // Warning: This routine calls Linop.Op, not LinOpo.HermOp
-  typedef Aggregation<vSpinColourVector,vTComplex,nbasis> Subspace;
-  Subspace Aggregates(Coarse4d,FGrid,cb);
-
-  NonHermitianLinearOperator<WilsonCloverFermionD,LatticeFermion> LinOpDw(Dw);
-  ShiftedNonHermitianLinearOperator<WilsonCloverFermionD,LatticeFermion> ShiftedLinOpDw(Dw,0.01);
-
-  Aggregates.CreateSubspaceGCR(RNG4,
-			       LinOpDw,
-			       nbasis);
-  
-
-  LittleDiracOperator LittleDiracOp(geom,FGrid,Coarse4d);
-  LittleDiracOp.CoarsenOperator(LinOpDw,Aggregates);
-
-  std::cout<<GridLogMessage<<std::endl;
-  std::cout<<GridLogMessage<<"*******************************************"<<std::endl;
-  std::cout<<GridLogMessage<<std::endl;
-  std::cout<<GridLogMessage<<"Testing coarsened operator "<<std::endl;
-  
-  CoarseVector c_src (Coarse4d);
-  CoarseVector c_res (Coarse4d);
-  CoarseVector c_proj(Coarse4d);
-
-  std::vector<LatticeFermion> subspace(nbasis,FGrid);
-  subspace=Aggregates.subspace;
-
-  c_src = one;  // 1 in every element for vector 1.
-  blockPromote(c_src,err,subspace);
-
-  prom=Zero();
-  for(int b=0;b<nbasis;b++){
-    prom=prom+subspace[b];
-  }
-  err=err-prom; 
-  std::cout<<GridLogMessage<<"Promoted back from subspace: err "<<norm2(err)<<std::endl;
-  std::cout<<GridLogMessage<<"c_src "<<norm2(c_src)<<std::endl;
-  std::cout<<GridLogMessage<<"prom  "<<norm2(prom)<<std::endl;
-
-  LinOpDw.Op(prom,tmp);
-  blockProject(c_proj,tmp,subspace);
-  std::cout<<GridLogMessage<<" Called Big Dirac Op "<<norm2(tmp)<<std::endl;
-
-  LittleDiracOp.M(c_src,c_res);
-  std::cout<<GridLogMessage<<" Called Little Dirac Op c_src "<< norm2(c_src) << "  c_res "<< norm2(c_res) <<std::endl;
-
-  std::cout<<GridLogMessage<<"Little dop : "<<norm2(c_res)<<std::endl;
-  //  std::cout<<GridLogMessage<<" Little "<< c_res<<std::endl;
-
-  std::cout<<GridLogMessage<<"Big dop in subspace : "<<norm2(c_proj)<<std::endl;
-  //  std::cout<<GridLogMessage<<" Big "<< c_proj<<std::endl;
-  c_proj = c_proj - c_res;
-  std::cout<<GridLogMessage<<" ldop error: "<<norm2(c_proj)<<std::endl;
-  //  std::cout<<GridLogMessage<<" error "<< c_proj<<std::endl;
-
-
-  /**********
-   * Some solvers
-   **********
-   */
-
-  // CG
-  {
-    MdagMLinearOperator<WilsonFermionD,LatticeFermion> HermOp(Dw);
-    ConjugateGradient<LatticeFermion> CG(1.0e-8,10000);
-    Dw.Mdag(src,precsrc);
-    CG(HermOp,precsrc,result);
-    result=Zero();
-  }
- 
-  ///////////////////////////////////////
-  // Coarse grid solver test
-  ///////////////////////////////////////
-
-  std::cout<<GridLogMessage<<"******************* "<<std::endl;
-  std::cout<<GridLogMessage<<" Coarse Grid Solve -- Level 3 "<<std::endl;
-  std::cout<<GridLogMessage<<"******************* "<<std::endl;
-  TrivialPrecon<CoarseVector> simple;
-  NonHermitianLinearOperator<LittleDiracOperator,CoarseVector> LinOpCoarse(LittleDiracOp);
-  ShiftedNonHermitianLinearOperator<LittleDiracOperator,CoarseVector> ShiftedLinOpCoarse(LittleDiracOp,0.001);
-  //  ShiftedNonHermitianLinearOperator<LittleDiracOperator,CoarseVector> ShiftedLinOpCoarse(LittleDiracOp,0.01);
-  //  ShiftedNonHermitianLinearOperator<LittleDiracOperator,CoarseVector> ShiftedLinOpCoarse(LinOpCoarse,0.001);
-  //  PrecGeneralisedConjugateResidualNonHermitian<CoarseVector>  L2PGCR(1.0e-4, 100, LinOpCoarse,simple,10,10); 
-  PrecGeneralisedConjugateResidualNonHermitian<CoarseVector>  L2PGCR(1.0e-1, 100, LinOpCoarse,simple,30,30); 
-  //  PrecGeneralisedConjugateResidualNonHermitian<CoarseVector>  L2PGCR(2.0e-1, 50, ShiftedLinOpCoarse,simple,50,50); 
-  L2PGCR.Level(3);
-  c_res=Zero();
-  L2PGCR(c_src,c_res);
-
-  ////////////////////////////////////////
-  // Fine grid smoother
-  ////////////////////////////////////////
-  std::cout<<GridLogMessage<<"******************* "<<std::endl;
-  std::cout<<GridLogMessage<<" Fine Grid Smoother -- Level 2 "<<std::endl;
-  std::cout<<GridLogMessage<<"******************* "<<std::endl;
-  TrivialPrecon<LatticeFermionD> simple_fine;
-
-  PrecGeneralisedConjugateResidualNonHermitian<LatticeFermionD> SmootherGCR(0.1,1,ShiftedLinOpDw,simple_fine,6,6);
-  SmootherGCR.Level(2);
-  
-  LatticeFermionD f_src(FGrid);
-  LatticeFermionD f_res(FGrid);
-
-  f_src = one;  // 1 in every element for vector 1.
-  f_res=Zero();
-  SmootherGCR(f_src,f_res);
-
-  typedef MGPreconditioner<vSpinColourVector,  vTComplex,nbasis> TwoLevelMG;
-
-  TwoLevelMG TwoLevelPrecon(Aggregates,
-			    LinOpDw,
-			    simple_fine,
-			    SmootherGCR,
-			    LinOpCoarse,
-			    L2PGCR);
-  
-  PrecGeneralisedConjugateResidualNonHermitian<LatticeFermion> L1PGCR(1.0e-8,1000,LinOpDw,TwoLevelPrecon,16,16);
-  L1PGCR.Level(1);
-
-  f_res=Zero();
-  L1PGCR(f_src,f_res);
-
-  std::cout<<GridLogMessage<<std::endl;
-  std::cout<<GridLogMessage<<"*******************************************"<<std::endl;
-  std::cout<<GridLogMessage<<std::endl;
-  std::cout<<GridLogMessage << "Done "<< std::endl;
-
-  Grid_finalize();
-  return 0;
-}

tests/debug/Test_general_coarse_wilson_svd.cc

@@ -1,320 +0,0 @@
-/*************************************************************************************
-
-    Grid physics library, www.github.com/paboyle/Grid 
-
-    Source file: ./tests/Test_padded_cell.cc
-
-    Copyright (C) 2023
-
-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>
-#include <Grid/lattice/PaddedCell.h>
-#include <Grid/stencil/GeneralLocalStencil.h>
-
-#include <Grid/algorithms/iterative/PrecGeneralisedConjugateResidual.h>
-#include <Grid/algorithms/iterative/PrecGeneralisedConjugateResidualNonHermitian.h>
-#include <Grid/algorithms/iterative/BiCGSTAB.h>
-
-using namespace std;
-using namespace Grid;
-
-template<class Fobj,class CComplex,int nbasis>
-class MGPreconditioner : public LinearFunction< Lattice<Fobj> > {
-public:
-  using LinearFunction<Lattice<Fobj> >::operator();
-
-  typedef Aggregation<Fobj,CComplex,nbasis> Aggregates;
-  typedef typename Aggregation<Fobj,CComplex,nbasis>::FineField    FineField;
-  typedef typename Aggregation<Fobj,CComplex,nbasis>::CoarseVector CoarseVector;
-  typedef typename Aggregation<Fobj,CComplex,nbasis>::CoarseMatrix CoarseMatrix;
-  typedef LinearOperatorBase<FineField>                            FineOperator;
-  typedef LinearFunction    <FineField>                            FineSmoother;
-  typedef LinearOperatorBase<CoarseVector>                         CoarseOperator;
-  typedef LinearFunction    <CoarseVector>                         CoarseSolver;
-  Aggregates     & _Aggregates;
-  FineOperator   & _FineOperator;
-  FineSmoother   & _PreSmoother;
-  FineSmoother   & _PostSmoother;
-  CoarseOperator & _CoarseOperator;
-  CoarseSolver   & _CoarseSolve;
-
-  int    level;  void Level(int lv) {level = lv; };
-
-  MGPreconditioner(Aggregates &Agg,
-		   FineOperator &Fine,
-		   FineSmoother &PreSmoother,
-		   FineSmoother &PostSmoother,
-		   CoarseOperator &CoarseOperator_,
-		   CoarseSolver &CoarseSolve_)
-    : _Aggregates(Agg),
-      _FineOperator(Fine),
-      _PreSmoother(PreSmoother),
-      _PostSmoother(PostSmoother),
-      _CoarseOperator(CoarseOperator_),
-      _CoarseSolve(CoarseSolve_),
-      level(1)  {  }
-
-  virtual void operator()(const FineField &in, FineField & out) 
-  {
-    GridBase *CoarseGrid = _Aggregates.CoarseGrid;
-    //    auto CoarseGrid = _CoarseOperator.Grid();
-    CoarseVector Csrc(CoarseGrid);
-    CoarseVector Csol(CoarseGrid);
-    FineField vec1(in.Grid());
-    FineField vec2(in.Grid());
-
-    std::cout<<GridLogMessage << "Calling PreSmoother " <<std::endl;
-
-    //    std::cout<<GridLogMessage << "Calling PreSmoother input residual "<<norm2(in) <<std::endl;
-    double t;
-    // Fine Smoother
-    //    out = in;
-    out = Zero();
-    t=-usecond();
-    _PreSmoother(in,out);
-    t+=usecond();
-
-    std::cout<<GridLogMessage << "PreSmoother took "<< t/1000.0<< "ms" <<std::endl;
-
-    // Update the residual
-    _FineOperator.Op(out,vec1);  sub(vec1, in ,vec1);   
-    //    std::cout<<GridLogMessage <<"Residual-1 now " <<norm2(vec1)<<std::endl;
-
-    // Fine to Coarse 
-    t=-usecond();
-    _Aggregates.ProjectToSubspace  (Csrc,vec1);
-    t+=usecond();
-    std::cout<<GridLogMessage << "Project to coarse took "<< t/1000.0<< "ms" <<std::endl;
-
-    // Coarse correction
-    t=-usecond();
-    Csol = Zero();
-    _CoarseSolve(Csrc,Csol);
-    //Csol=Zero();
-    t+=usecond();
-    std::cout<<GridLogMessage << "Coarse solve took "<< t/1000.0<< "ms" <<std::endl;
-
-    // Coarse to Fine
-    t=-usecond();  
-    //    _CoarseOperator.PromoteFromSubspace(_Aggregates,Csol,vec1);
-    _Aggregates.PromoteFromSubspace(Csol,vec1); 
-    add(out,out,vec1);
-    t+=usecond();
-    std::cout<<GridLogMessage << "Promote to this level took "<< t/1000.0<< "ms" <<std::endl;
-
-    // Residual
-    _FineOperator.Op(out,vec1);  sub(vec1 ,in , vec1);  
-    //    std::cout<<GridLogMessage <<"Residual-2 now " <<norm2(vec1)<<std::endl;
-
-    // Fine Smoother
-    t=-usecond();
-    //    vec2=vec1;
-    vec2=Zero();
-    _PostSmoother(vec1,vec2);
-    t+=usecond();
-    std::cout<<GridLogMessage << "PostSmoother took "<< t/1000.0<< "ms" <<std::endl;
-
-    add( out,out,vec2);
-    std::cout<<GridLogMessage << "Done " <<std::endl;
-  }
-};
-
-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   = UGrid;
-  GridRedBlackCartesian * FrbGrid = UrbGrid;
-
-  // Construct a coarsened grid
-  Coordinate clatt = GridDefaultLatt();
-  for(int d=0;d<clatt.size();d++){
-    clatt[d] = clatt[d]/2;
-    //    clatt[d] = clatt[d]/4;
-  }
-  GridCartesian *Coarse4d =  SpaceTimeGrid::makeFourDimGrid(clatt, GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi());;
-
-  std::vector<int> seeds4({1,2,3,4});
-  std::vector<int> cseeds({5,6,7,8});
-  GridParallelRNG          RNG4(UGrid);   RNG4.SeedFixedIntegers(seeds4);
-  GridParallelRNG          CRNG(Coarse4d);CRNG.SeedFixedIntegers(cseeds);
-
-  LatticeFermion    src(FGrid); random(RNG4,src);
-  LatticeFermion result(FGrid); result=Zero();
-  LatticeFermion    ref(FGrid); ref=Zero();
-  LatticeFermion    tmp(FGrid);
-  LatticeFermion    err(FGrid);
-  LatticeGaugeField Umu(UGrid);
-
-  FieldMetaData header;
-  std::string file("ckpoint_lat");
-  NerscIO::readConfiguration(Umu,header,file);
-  
-  RealD csw =0.0;
-  RealD mass=-0.92;
-
-  WilsonCloverFermionD Dw(Umu,*UGrid,*UrbGrid,mass,csw,csw);
-
-  const int nbasis = 20;
-  const int cb = 0 ;
-  LatticeFermion prom(FGrid);
-
-  typedef GeneralCoarsenedMatrix<vSpinColourVector,vTComplex,2*nbasis> LittleDiracOperator;
-  typedef LittleDiracOperator::CoarseVector CoarseVector;
-
-  NearestStencilGeometry4D geom(Coarse4d);
-
-  std::cout<<GridLogMessage<<std::endl;
-  std::cout<<GridLogMessage<<"*******************************************"<<std::endl;
-  std::cout<<GridLogMessage<<std::endl;
- 
-  // Warning: This routine calls Linop.Op, not LinOpo.HermOp
-  typedef Aggregation<vSpinColourVector,vTComplex,nbasis> Subspace;
-  Subspace Aggregates(Coarse4d,FGrid,cb);
-
-  MdagMLinearOperator<WilsonCloverFermionD,LatticeFermion> MdagMOpDw(Dw);
-  NonHermitianLinearOperator<WilsonCloverFermionD,LatticeFermion> LinOpDw(Dw);
-  ShiftedNonHermitianLinearOperator<WilsonCloverFermionD,LatticeFermion> ShiftedLinOpDw(Dw,0.5);
-
-  //  Aggregates.CreateSubspaceGCR(RNG4,
-  //			       LinOpDw,
-  //			       nbasis);
-  Aggregates.CreateSubspace(RNG4,MdagMOpDw,nbasis);
-  
-  typedef Aggregation<vSpinColourVector,vTComplex,2*nbasis> CombinedSubspace;
-  CombinedSubspace CombinedUV(Coarse4d,UGrid,cb);
-  for(int b=0;b<nbasis;b++){
-    Gamma G5(Gamma::Algebra::Gamma5);
-    CombinedUV.subspace[b]        = Aggregates.subspace[b];
-    CombinedUV.subspace[b+nbasis] = G5*Aggregates.subspace[b];
-  }
-
-  LittleDiracOperator LittleDiracOp(geom,FGrid,Coarse4d);
-  LittleDiracOp.CoarsenOperator(LinOpDw,CombinedUV);
-
-  std::cout<<GridLogMessage<<std::endl;
-  std::cout<<GridLogMessage<<"*******************************************"<<std::endl;
-  std::cout<<GridLogMessage<<std::endl;
-  std::cout<<GridLogMessage<<"Testing coarsened operator "<<std::endl;
-  
-  CoarseVector c_src (Coarse4d);
-  CoarseVector c_res (Coarse4d);
-  CoarseVector c_proj(Coarse4d);
-
-  std::vector<LatticeFermion> subspace(2*nbasis,FGrid);
-  subspace=CombinedUV.subspace;
-
-  Complex one(1.0);
-  c_src = one;  // 1 in every element for vector 1.
-  blockPromote(c_src,err,subspace);
-
-  prom=Zero();
-  for(int b=0;b<2*nbasis;b++){
-    prom=prom+subspace[b];
-  }
-  err=err-prom; 
-  std::cout<<GridLogMessage<<"Promoted back from subspace: err "<<norm2(err)<<std::endl;
-  std::cout<<GridLogMessage<<"c_src "<<norm2(c_src)<<std::endl;
-  std::cout<<GridLogMessage<<"prom  "<<norm2(prom)<<std::endl;
-
-  LinOpDw.Op(prom,tmp);
-  blockProject(c_proj,tmp,subspace);
-  std::cout<<GridLogMessage<<" Called Big Dirac Op "<<norm2(tmp)<<std::endl;
-
-  LittleDiracOp.M(c_src,c_res);
-  std::cout<<GridLogMessage<<" Called Little Dirac Op c_src "<< norm2(c_src) << "  c_res "<< norm2(c_res) <<std::endl;
-
-  std::cout<<GridLogMessage<<"Little dop : "<<norm2(c_res)<<std::endl;
-  //  std::cout<<GridLogMessage<<" Little "<< c_res<<std::endl;
-
-  std::cout<<GridLogMessage<<"Big dop in subspace : "<<norm2(c_proj)<<std::endl;
-  //  std::cout<<GridLogMessage<<" Big "<< c_proj<<std::endl;
-  c_proj = c_proj - c_res;
-  std::cout<<GridLogMessage<<" ldop error: "<<norm2(c_proj)<<std::endl;
-  //  std::cout<<GridLogMessage<<" error "<< c_proj<<std::endl;
-
-
-  /**********
-   * Some solvers
-   **********
-   */
-
-  ///////////////////////////////////////
-  // Coarse grid solver test
-  ///////////////////////////////////////
-
-  std::cout<<GridLogMessage<<"******************* "<<std::endl;
-  std::cout<<GridLogMessage<<" Coarse Grid Solve -- Level 3 "<<std::endl;
-  std::cout<<GridLogMessage<<"******************* "<<std::endl;
-  TrivialPrecon<CoarseVector> simple;
-  NonHermitianLinearOperator<LittleDiracOperator,CoarseVector> LinOpCoarse(LittleDiracOp);
-  //  PrecGeneralisedConjugateResidualNonHermitian<CoarseVector>  L2PGCR(1.0e-4, 100, LinOpCoarse,simple,10,10); 
-  PrecGeneralisedConjugateResidualNonHermitian<CoarseVector>  L2PGCR(1.0e-2, 100, LinOpCoarse,simple,30,30); 
-  L2PGCR.Level(3);
-  c_res=Zero();
-  L2PGCR(c_src,c_res);
-
-  ////////////////////////////////////////
-  // Fine grid smoother
-  ////////////////////////////////////////
-  std::cout<<GridLogMessage<<"******************* "<<std::endl;
-  std::cout<<GridLogMessage<<" Fine Grid Smoother -- Level 2 "<<std::endl;
-  std::cout<<GridLogMessage<<"******************* "<<std::endl;
-  TrivialPrecon<LatticeFermionD> simple_fine;
-
-  PrecGeneralisedConjugateResidualNonHermitian<LatticeFermionD> SmootherGCR(0.01,1,ShiftedLinOpDw,simple_fine,4,4);
-  SmootherGCR.Level(2);
-  
-  LatticeFermionD f_src(FGrid);
-  LatticeFermionD f_res(FGrid);
-
-  f_src = one;  // 1 in every element for vector 1.
-  f_res=Zero();
-  SmootherGCR(f_src,f_res);
-
-  typedef MGPreconditioner<vSpinColourVector,  vTComplex,2*nbasis> TwoLevelMG;
-
-  TwoLevelMG TwoLevelPrecon(CombinedUV,
-			    LinOpDw,
-			    simple_fine,
-			    SmootherGCR,
-			    LinOpCoarse,
-			    L2PGCR);
-  
-  PrecGeneralisedConjugateResidualNonHermitian<LatticeFermion> L1PGCR(1.0e-8,1000,LinOpDw,TwoLevelPrecon,32,32);
-  L1PGCR.Level(1);
-
-  f_res=Zero();
-  L1PGCR(f_src,f_res);
-
-  std::cout<<GridLogMessage<<std::endl;
-  std::cout<<GridLogMessage<<"*******************************************"<<std::endl;
-  std::cout<<GridLogMessage<<std::endl;
-  std::cout<<GridLogMessage << "Done "<< std::endl;
-
-  Grid_finalize();
-  return 0;
-}

tests/debug/Test_general_coarse_wilson_svd_no5g.cc

@@ -1,312 +0,0 @@
-/*************************************************************************************
-
-    Grid physics library, www.github.com/paboyle/Grid 
-
-    Source file: ./tests/Test_padded_cell.cc
-
-    Copyright (C) 2023
-
-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>
-#include <Grid/lattice/PaddedCell.h>
-#include <Grid/stencil/GeneralLocalStencil.h>
-
-#include <Grid/algorithms/iterative/PrecGeneralisedConjugateResidual.h>
-#include <Grid/algorithms/iterative/PrecGeneralisedConjugateResidualNonHermitian.h>
-#include <Grid/algorithms/iterative/BiCGSTAB.h>
-
-using namespace std;
-using namespace Grid;
-
-template<class Fobj,class CComplex,int nbasis>
-class MGPreconditioner : public LinearFunction< Lattice<Fobj> > {
-public:
-  using LinearFunction<Lattice<Fobj> >::operator();
-
-  typedef Aggregation<Fobj,CComplex,nbasis> Aggregates;
-  typedef typename Aggregation<Fobj,CComplex,nbasis>::FineField    FineField;
-  typedef typename Aggregation<Fobj,CComplex,nbasis>::CoarseVector CoarseVector;
-  typedef typename Aggregation<Fobj,CComplex,nbasis>::CoarseMatrix CoarseMatrix;
-  typedef LinearOperatorBase<FineField>                            FineOperator;
-  typedef LinearFunction    <FineField>                            FineSmoother;
-  typedef LinearOperatorBase<CoarseVector>                         CoarseOperator;
-  typedef LinearFunction    <CoarseVector>                         CoarseSolver;
-  Aggregates     & _Aggregates;
-  FineOperator   & _FineOperator;
-  FineSmoother   & _PreSmoother;
-  FineSmoother   & _PostSmoother;
-  CoarseOperator & _CoarseOperator;
-  CoarseSolver   & _CoarseSolve;
-
-  int    level;  void Level(int lv) {level = lv; };
-
-  MGPreconditioner(Aggregates &Agg,
-		   FineOperator &Fine,
-		   FineSmoother &PreSmoother,
-		   FineSmoother &PostSmoother,
-		   CoarseOperator &CoarseOperator_,
-		   CoarseSolver &CoarseSolve_)
-    : _Aggregates(Agg),
-      _FineOperator(Fine),
-      _PreSmoother(PreSmoother),
-      _PostSmoother(PostSmoother),
-      _CoarseOperator(CoarseOperator_),
-      _CoarseSolve(CoarseSolve_),
-      level(1)  {  }
-
-  virtual void operator()(const FineField &in, FineField & out) 
-  {
-    GridBase *CoarseGrid = _Aggregates.CoarseGrid;
-    //    auto CoarseGrid = _CoarseOperator.Grid();
-    CoarseVector Csrc(CoarseGrid);
-    CoarseVector Csol(CoarseGrid);
-    FineField vec1(in.Grid());
-    FineField vec2(in.Grid());
-
-    std::cout<<GridLogMessage << "Calling PreSmoother " <<std::endl;
-
-    //    std::cout<<GridLogMessage << "Calling PreSmoother input residual "<<norm2(in) <<std::endl;
-    double t;
-    // Fine Smoother
-    //    out = in;
-    out = Zero();
-    t=-usecond();
-    _PreSmoother(in,out);
-    t+=usecond();
-
-    std::cout<<GridLogMessage << "PreSmoother took "<< t/1000.0<< "ms" <<std::endl;
-
-    // Update the residual
-    _FineOperator.Op(out,vec1);  sub(vec1, in ,vec1);   
-    //    std::cout<<GridLogMessage <<"Residual-1 now " <<norm2(vec1)<<std::endl;
-
-    // Fine to Coarse 
-    t=-usecond();
-    _Aggregates.ProjectToSubspace  (Csrc,vec1);
-    t+=usecond();
-    std::cout<<GridLogMessage << "Project to coarse took "<< t/1000.0<< "ms" <<std::endl;
-
-    // Coarse correction
-    t=-usecond();
-    Csol = Zero();
-    _CoarseSolve(Csrc,Csol);
-    //Csol=Zero();
-    t+=usecond();
-    std::cout<<GridLogMessage << "Coarse solve took "<< t/1000.0<< "ms" <<std::endl;
-
-    // Coarse to Fine
-    t=-usecond();  
-    //    _CoarseOperator.PromoteFromSubspace(_Aggregates,Csol,vec1);
-    _Aggregates.PromoteFromSubspace(Csol,vec1); 
-    add(out,out,vec1);
-    t+=usecond();
-    std::cout<<GridLogMessage << "Promote to this level took "<< t/1000.0<< "ms" <<std::endl;
-
-    // Residual
-    _FineOperator.Op(out,vec1);  sub(vec1 ,in , vec1);  
-    //    std::cout<<GridLogMessage <<"Residual-2 now " <<norm2(vec1)<<std::endl;
-
-    // Fine Smoother
-    t=-usecond();
-    //    vec2=vec1;
-    vec2=Zero();
-    _PostSmoother(vec1,vec2);
-    t+=usecond();
-    std::cout<<GridLogMessage << "PostSmoother took "<< t/1000.0<< "ms" <<std::endl;
-
-    add( out,out,vec2);
-    std::cout<<GridLogMessage << "Done " <<std::endl;
-  }
-};
-
-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   = UGrid;
-  GridRedBlackCartesian * FrbGrid = UrbGrid;
-
-  // Construct a coarsened grid
-  Coordinate clatt = GridDefaultLatt();
-  for(int d=0;d<clatt.size();d++){
-    clatt[d] = clatt[d]/2;
-    //    clatt[d] = clatt[d]/4;
-  }
-  GridCartesian *Coarse4d =  SpaceTimeGrid::makeFourDimGrid(clatt, GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi());;
-
-  std::vector<int> seeds4({1,2,3,4});
-  std::vector<int> cseeds({5,6,7,8});
-  GridParallelRNG          RNG4(UGrid);   RNG4.SeedFixedIntegers(seeds4);
-  GridParallelRNG          CRNG(Coarse4d);CRNG.SeedFixedIntegers(cseeds);
-
-  LatticeFermion    src(FGrid); random(RNG4,src);
-  LatticeFermion result(FGrid); result=Zero();
-  LatticeFermion    ref(FGrid); ref=Zero();
-  LatticeFermion    tmp(FGrid);
-  LatticeFermion    err(FGrid);
-  LatticeGaugeField Umu(UGrid);
-
-  FieldMetaData header;
-  std::string file("ckpoint_lat");
-  NerscIO::readConfiguration(Umu,header,file);
-  
-  RealD csw =0.0;
-  RealD mass=-0.92;
-
-  WilsonCloverFermionD Dw(Umu,*UGrid,*UrbGrid,mass,csw,csw);
-
-  const int nbasis = 40;
-  const int cb = 0 ;
-  LatticeFermion prom(FGrid);
-
-  typedef GeneralCoarsenedMatrix<vSpinColourVector,vTComplex,nbasis> LittleDiracOperator;
-  typedef LittleDiracOperator::CoarseVector CoarseVector;
-
-  NearestStencilGeometry4D geom(Coarse4d);
-
-  std::cout<<GridLogMessage<<std::endl;
-  std::cout<<GridLogMessage<<"*******************************************"<<std::endl;
-  std::cout<<GridLogMessage<<std::endl;
- 
-  // Warning: This routine calls Linop.Op, not LinOpo.HermOp
-  typedef Aggregation<vSpinColourVector,vTComplex,nbasis> Subspace;
-  Subspace Aggregates(Coarse4d,FGrid,cb);
-
-  MdagMLinearOperator<WilsonCloverFermionD,LatticeFermion> MdagMOpDw(Dw);
-  NonHermitianLinearOperator<WilsonCloverFermionD,LatticeFermion> LinOpDw(Dw);
-  ShiftedNonHermitianLinearOperator<WilsonCloverFermionD,LatticeFermion> ShiftedLinOpDw(Dw,0.5);
-
-  //  Aggregates.CreateSubspaceGCR(RNG4,
-  //			       LinOpDw,
-  //			       nbasis);
-  Aggregates.CreateSubspace(RNG4,MdagMOpDw,nbasis);
-
-  LittleDiracOperator LittleDiracOp(geom,FGrid,Coarse4d);
-  LittleDiracOp.CoarsenOperator(LinOpDw,Aggregates);
-
-  std::cout<<GridLogMessage<<std::endl;
-  std::cout<<GridLogMessage<<"*******************************************"<<std::endl;
-  std::cout<<GridLogMessage<<std::endl;
-  std::cout<<GridLogMessage<<"Testing coarsened operator "<<std::endl;
-  
-  CoarseVector c_src (Coarse4d);
-  CoarseVector c_res (Coarse4d);
-  CoarseVector c_proj(Coarse4d);
-
-  std::vector<LatticeFermion> subspace(nbasis,FGrid);
-  subspace=Aggregates.subspace;
-
-  Complex one(1.0);
-  c_src = one;  // 1 in every element for vector 1.
-  blockPromote(c_src,err,subspace);
-
-  prom=Zero();
-  for(int b=0;b<nbasis;b++){
-    prom=prom+subspace[b];
-  }
-  err=err-prom; 
-  std::cout<<GridLogMessage<<"Promoted back from subspace: err "<<norm2(err)<<std::endl;
-  std::cout<<GridLogMessage<<"c_src "<<norm2(c_src)<<std::endl;
-  std::cout<<GridLogMessage<<"prom  "<<norm2(prom)<<std::endl;
-
-  LinOpDw.Op(prom,tmp);
-  blockProject(c_proj,tmp,subspace);
-  std::cout<<GridLogMessage<<" Called Big Dirac Op "<<norm2(tmp)<<std::endl;
-
-  LittleDiracOp.M(c_src,c_res);
-  std::cout<<GridLogMessage<<" Called Little Dirac Op c_src "<< norm2(c_src) << "  c_res "<< norm2(c_res) <<std::endl;
-
-  std::cout<<GridLogMessage<<"Little dop : "<<norm2(c_res)<<std::endl;
-  //  std::cout<<GridLogMessage<<" Little "<< c_res<<std::endl;
-
-  std::cout<<GridLogMessage<<"Big dop in subspace : "<<norm2(c_proj)<<std::endl;
-  //  std::cout<<GridLogMessage<<" Big "<< c_proj<<std::endl;
-  c_proj = c_proj - c_res;
-  std::cout<<GridLogMessage<<" ldop error: "<<norm2(c_proj)<<std::endl;
-  //  std::cout<<GridLogMessage<<" error "<< c_proj<<std::endl;
-
-
-  /**********
-   * Some solvers
-   **********
-   */
-
-  ///////////////////////////////////////
-  // Coarse grid solver test
-  ///////////////////////////////////////
-
-  std::cout<<GridLogMessage<<"******************* "<<std::endl;
-  std::cout<<GridLogMessage<<" Coarse Grid Solve -- Level 3 "<<std::endl;
-  std::cout<<GridLogMessage<<"******************* "<<std::endl;
-  TrivialPrecon<CoarseVector> simple;
-  NonHermitianLinearOperator<LittleDiracOperator,CoarseVector> LinOpCoarse(LittleDiracOp);
-  //  PrecGeneralisedConjugateResidualNonHermitian<CoarseVector>  L2PGCR(1.0e-4, 100, LinOpCoarse,simple,10,10); 
-  PrecGeneralisedConjugateResidualNonHermitian<CoarseVector>  L2PGCR(1.0e-2, 100, LinOpCoarse,simple,30,30); 
-  L2PGCR.Level(3);
-  c_res=Zero();
-  L2PGCR(c_src,c_res);
-
-  ////////////////////////////////////////
-  // Fine grid smoother
-  ////////////////////////////////////////
-  std::cout<<GridLogMessage<<"******************* "<<std::endl;
-  std::cout<<GridLogMessage<<" Fine Grid Smoother -- Level 2 "<<std::endl;
-  std::cout<<GridLogMessage<<"******************* "<<std::endl;
-  TrivialPrecon<LatticeFermionD> simple_fine;
-
-  PrecGeneralisedConjugateResidualNonHermitian<LatticeFermionD> SmootherGCR(0.01,1,ShiftedLinOpDw,simple_fine,6,6);
-  SmootherGCR.Level(2);
-  
-  LatticeFermionD f_src(FGrid);
-  LatticeFermionD f_res(FGrid);
-
-  f_src = one;  // 1 in every element for vector 1.
-  f_res=Zero();
-  SmootherGCR(f_src,f_res);
-
-  typedef MGPreconditioner<vSpinColourVector,  vTComplex,nbasis> TwoLevelMG;
-
-  TwoLevelMG TwoLevelPrecon(Aggregates,
-			    LinOpDw,
-			    simple_fine,
-			    SmootherGCR,
-			    LinOpCoarse,
-			    L2PGCR);
-  
-  PrecGeneralisedConjugateResidualNonHermitian<LatticeFermion> L1PGCR(1.0e-8,1000,LinOpDw,TwoLevelPrecon,32,32);
-  L1PGCR.Level(1);
-
-  f_res=Zero();
-  L1PGCR(f_src,f_res);
-
-  std::cout<<GridLogMessage<<std::endl;
-  std::cout<<GridLogMessage<<"*******************************************"<<std::endl;
-  std::cout<<GridLogMessage<<std::endl;
-  std::cout<<GridLogMessage << "Done "<< std::endl;
-
-  Grid_finalize();
-  return 0;
-}

tests/debug/Test_padded_cell_staple.cc

@@ -490,7 +490,7 @@ public:
 
 	    }
 	  }
-	  assert(s==nshift);
+	  GRID_ASSERT(s==nshift);
 	  coalescedWrite(gStaple_v[ss],stencil_ss);
 	}
 	);

tests/disable_tests_without_instantiations.h

@@ -1,16 +0,0 @@
-#include <Grid/Grid.h>
-
-#pragma once
-
-#ifndef ENABLE_FERMION_INSTANTIATIONS
-#include <iostream>
-
-int main(void) {
-  std::cout << "This build of Grid was configured to exclude fermion instantiations, "
-	    << "which this test relies on. "
-	    << "Please reconfigure and rebuild Grid with --enable-fermion-instantiations"
-	    << "to run this test."
-	    << std::endl;
-  return 1;
-}
-#endif

tests/forces/Test_fthmc.cc

@@ -62,7 +62,6 @@ void ForceTest(Action<LatticeGaugeField> &action,ConfigurationBase<LatticeGaugeF
   
   Gimpl::generate_momenta(P,sRNG,RNG4);
   //  Filter.applyFilter(P);
-  std::cout << GridLogMessage << "Initial momenta " << norm2(P) << std::endl;
 
   action.refresh(smU,sRNG,RNG4);
 
@@ -71,8 +70,6 @@ void ForceTest(Action<LatticeGaugeField> &action,ConfigurationBase<LatticeGaugeF
   std::cout << GridLogMessage << "+++++++++++++++++++++++++++++++++++++++++++++++++++++++++"<<std::endl;
 
   RealD S1 = action.S(smU);
-  std::cout << GridLogMessage << "Initial action " << S1 << std::endl;
-
 
   Gimpl::update_field(P,U,eps);
   smU.set_Field(U);
@@ -83,7 +80,6 @@ void ForceTest(Action<LatticeGaugeField> &action,ConfigurationBase<LatticeGaugeF
   action.deriv(smU,UdSdU);
   UdSdU = Ta(UdSdU);
   //  Filter.applyFilter(UdSdU);
-  std::cout << GridLogMessage << "Derivative " << norm2(UdSdU) << std::endl;
 
   DumpSliceNorm("Force",UdSdU,Nd-1);
   
@@ -95,7 +91,6 @@ void ForceTest(Action<LatticeGaugeField> &action,ConfigurationBase<LatticeGaugeF
   std::cout << GridLogMessage << "+++++++++++++++++++++++++++++++++++++++++++++++++++++++++"<<std::endl;
   
   RealD S2 = action.S(smU);
-  std::cout << GridLogMessage << "Final action " << S1 << std::endl;
 
   // Use the derivative
   LatticeComplex dS(UGrid); dS = Zero();
@@ -150,8 +145,6 @@ int main (int argc, char ** argv)
   GridParallelRNG          RNG4(UGrid);  RNG4.SeedFixedIntegers(seeds);
   SU<Nc>::HotConfiguration(RNG4,U);
 #endif
-  std::cout << GridLogMessage << "Initial plaquette: " << WilsonLoops<PeriodicGimplR>::avgPlaquette(U) << std::endl;
-
 
   
   WilsonGaugeActionR  PlaqAction(6.0);

visualisation/Visualise5D.cxx

@@ -1,729 +0,0 @@
-// Derived from VTK/Examples/Cxx/Medical2.cxx
-// The example reads a volume dataset, extracts two isosurfaces that
-// represent the skin and bone, and then displays them.
-//
-// Modified heavily by Peter Boyle to display lattice field theory data as movies and compare multiple files
-
-#include <vtkActor.h>
-#include <vtkCamera.h>
-#include <vtkMetaImageReader.h>
-#include <vtkNamedColors.h>
-#include <vtkNew.h>
-#include <vtkOutlineFilter.h>
-#include <vtkPolyDataMapper.h>
-#include <vtkProperty.h>
-#include <vtkRenderWindow.h>
-#include <vtkRenderWindowInteractor.h>
-#include <vtkRenderer.h>
-#include <vtkStripper.h>
-#include <vtkImageData.h>
-#include <vtkVersion.h>
-#include <vtkCallbackCommand.h>
-#include <vtkTextActor.h>
-#include <vtkTextProperty.h>
-
-#define MPEG
-#ifdef MPEG
-#include <vtkFFMPEGWriter.h>
-#endif
-
-#include <vtkProperty2D.h>
-#include <vtkSliderWidget.h>
-#include <vtkSliderRepresentation2D.h>
-#include <vtkWindowToImageFilter.h>
-
-#include <array>
-#include <string>
-
-#include <Grid/Grid.h>
-
-#define USE_FLYING_EDGES
-#ifdef USE_FLYING_EDGES
-#include <vtkFlyingEdges3D.h>
-typedef vtkFlyingEdges3D isosurface;
-#else
-#include <vtkMarchingCubes.h>
-typedef vtkMarchingCubes isosurface;
-#endif
-
-int mpeg = 0 ;
-int framerate = 10;
-
-template <class T> void readFile(T& out, std::string const fname){
-  Grid::emptyUserRecord record;
-  Grid::ScidacReader RD;
-  RD.open(fname);
-  RD.readScidacFieldRecord(out,record);
-  RD.close();
-}
-using namespace Grid;
-
-class FrameUpdater : public vtkCallbackCommand
-{
-public:
-
-  FrameUpdater() {
-    ffile=0;
-    TimerCount = 0;
-    xoff       = 0;
-    t          = 0;
-    imageData = nullptr;
-    timerId = 0;
-    maxCount = -1;
-    old_file=-1;
-  }
-  
-  static FrameUpdater* New()
-  {
-    FrameUpdater* cb = new FrameUpdater;
-    cb->TimerCount = 0;
-    return cb;
-  }
-
-  //
-  // Must map a x,y,z + frame index into
-  // i)  a d-dimensional site Coordinate
-  // ii) a file name
-  // Need a:
-  //     loop_ranges
-  //     sum_ranges
-  //     loop_vol  -- map loop_idx -> loop_coor
-  //     sum_vol   -- map sum_idx -> sum_coor with Lexicographic
-  //
-  /*
-   * Just set this up
-   */
-  int old_file ; // Cache, avoid reread
-  
-  Coordinate          latt;
-  Coordinate          xyz_dims    ; // List lattice dimensions corresponding to xyz_dims displayed 
-  Coordinate          xyz_ranges  ; // 3-vector
-  Coordinate          g_xyz_ranges; // Nd-vector
-  uint64_t            xyz_vol     ;
-  
-  Coordinate          loop_dims;    // List lattice dimensions put into movie time
-  Coordinate          loop_ranges;  // movie time ranges
-  uint64_t            loop_vol;
-
-  Coordinate          sum_dims;      // List lattice dimensions summed
-  Coordinate          sum_ranges;    // summation ranges
-  uint64_t            sum_vol;
-
-  Coordinate          slice_dims;      // List slice dimensions 
-
-  Coordinate          Slice;
-  
-  std::vector<std::string> files;        // file list that is looped over
-
-  int Nd;
-  GridBase *grid;
-  Grid::LatticeComplexD *grid_data;
-
-  void SetGrid(GridBase *_grid)
-  {
-    grid = _grid;
-    Nd=grid->Nd();
-    latt = grid->GlobalDimensions();
-    grid_data = new Grid::LatticeComplexD(grid);
-  }
-
-  void SetFiles(std::vector<std::string> list)      { files = list; old_file = -1; }
-  void SetSlice(Coordinate _Slice)             { Slice = _Slice;} // Offset / skew for lattice coords
-  void SetSumDimensions (Coordinate _SumDims ) {
-    sum_ranges=Coordinate(Nd);
-    sum_dims = _SumDims; // 1 hot for dimensions summed
-    sum_vol = 1;
-    for(int d=0;d<sum_dims.size();d++){
-      if ( sum_dims[d] == 1 ) sum_ranges[d] = latt[d];
-      else                    sum_ranges[d] = 1;
-      sum_vol*=sum_ranges[d];
-    }
-  }
-  void SetLoopDimensions(Coordinate _LoopDims) {
-    loop_ranges=Coordinate(Nd);
-    loop_dims= _LoopDims;
-    loop_vol = 1;
-    for(int d=0;d<loop_dims.size();d++){
-      if ( loop_dims[d] == 1 ) loop_ranges[d] = latt[d];
-      else                     loop_ranges[d] = 1;
-      loop_vol*=loop_ranges[d];
-    }
-  } // 
-  void SetDisplayDimensions(Coordinate _xyz_dims   ) {
-    g_xyz_ranges=Coordinate(Nd);
-    xyz_ranges=Coordinate(3);
-    xyz_dims = _xyz_dims;
-    xyz_vol  = 1;
-    for(int d=0;d<3;d++){
-      xyz_ranges[d] = latt[xyz_dims[d]];
-      xyz_vol *= xyz_ranges[d];
-    }
-    // Find dim extents for grid
-    int dd=0;
-    for(int d=0;d<Nd;d++){
-      g_xyz_ranges[d] = 1;
-      for(int dd=0;dd<3;dd++) {
-	if ( xyz_dims[dd]==d ) {
-	  g_xyz_ranges[d] = latt[d];
-	}
-      }
-    }
-  }
-  void SetSliceDimensions(void) {
-    Coordinate _slice_dims;
-    for ( int d=0;d<Nd;d++){
-      int is_slice = 1;
-      if(g_xyz_ranges[d]>1) is_slice = 0;
-      if(loop_dims[d]) is_slice = 0;
-      if(sum_dims[d] ) is_slice = 0;
-      if(is_slice) _slice_dims.push_back(d);
-    }
-    slice_dims = _slice_dims;
-    std::cout << " Setting Slice Dimensions to "<<slice_dims<<std::endl;
-  }
-  virtual void Execute(vtkObject* caller, unsigned long eventId,void* vtkNotUsed(callData))
-  {
-    const int max=256;
-    char text_string[max];
-
-    auto latt_size = grid->GlobalDimensions();
-    
-    if ( vtkCommand::KeyPressEvent == eventId ) {
-
-      vtkRenderWindowInteractor* iren = static_cast<vtkRenderWindowInteractor*>(caller);
-      std::string key = iren->GetKeySym();
-      std::cout << "Pressed: " << key << std::endl;
-      if (slice_dims.size()>0) {
-	int vert = slice_dims[slice_dims.size()-1];
-	int horz = slice_dims[0];
-	if ( key == "Up" ) {
-	  Slice[vert] = (Slice[vert]+1)%latt[vert];
-	}
-	if ( key == "Down" ) {
-	  Slice[vert] = (Slice[vert]+latt[vert]-1)%latt[vert];
-	}
-	if ( key == "Right" ) {
-	  Slice[horz] = (Slice[horz]+1)%latt[horz];
-	}
-	if ( key == "Left" ) {
-	  Slice[horz] = (Slice[horz]+latt[horz]-1)%latt[horz];
-	}
-      }
-      if ( key == "greater" ) {
-	ffile = (ffile + 1) % files.size();
-      }
-      if ( key == "less" ) {
-	ffile = (ffile - 1 + files.size()) % files.size();
-      }
-      std::cout <<"Slice " <<Slice <<std::endl;
-      std::cout <<"File  " <<ffile <<std::endl;
-    }
-      
-    // Make a new frame for frame index TimerCount
-    if ( vtkCommand::TimerEvent == eventId || vtkCommand::KeyPressEvent == eventId)
-      {
-	int file     = ((this->TimerCount / loop_vol) + ffile )%files.size();
-
-	if ( file != old_file ) {
-	  readFile(*grid_data,files[file]);
-	  old_file = file;
-	}
-
-	RealD max, min, max_abs,min_abs;
-	Coordinate max_site;
-	Coordinate min_site;
-	Coordinate max_abs_site;
-	Coordinate min_abs_site;
-	for(int idx=0;idx<grid->gSites();idx++){
-	  Coordinate site;
-	  Lexicographic::CoorFromIndex (site,idx,latt);
-	  RealD val=real(peekSite(*grid_data,site));
-	  if (idx==0){
-	    max = min = val;
-	    max_abs = min_abs = fabs(val);
-	    max_site = site;
-	    min_site = site;
-	    min_abs_site = site;
-	    max_abs_site = site;
-	  } else {
-	    if ( val > max ) {
-	      max=val;
-	      max_site = site;
-	    }
-	    if ( fabs(val) > max_abs ) {
-	      max_abs=fabs(val);
-	      max_abs_site = site;
-	    }
-	    if ( val < min ) {
-	      min=val;
-	      min_site = site;
-	    }	    
-	    if ( fabs(val) < min_abs ) {
-	      min_abs=fabs(val);
-	      min_abs_site = site;
-	    }	    
-	  }
-	}
-	std::cout << " abs_max "<<max_abs<<" at " << max_abs_site<<std::endl;
-	std::cout << " abs_min "<<min_abs<<" at " << min_abs_site<<std::endl;
-	std::cout << " max "<<max<<" at " << max_site<<std::endl;
-	std::cout << " min "<<min<<" at " << min_site<<std::endl;
-	// Looped dimensions, map index to coordinate
-	int loop_idx = this->TimerCount % loop_vol;
-	Coordinate loop_coor;
-	Lexicographic::CoorFromIndex (loop_coor,loop_idx,loop_ranges);
-	
-	// Loop over xyz sites
-	Coordinate xyz_coor(3);
-	Coordinate g_xyz_coor(Nd);
-	Coordinate sum_coor(Nd);
-	for(uint64_t xyz = 0 ; xyz< xyz_vol; xyz++){
-	  Lexicographic::CoorFromIndex (xyz_coor,xyz,xyz_ranges);
-	  Lexicographic::CoorFromIndex (g_xyz_coor,xyz,g_xyz_ranges);
-	  
-	  RealD sum_value = 0.0;
-	  for(uint64_t sum_idx = 0 ; sum_idx< sum_vol; sum_idx++){
-	    Lexicographic::CoorFromIndex (sum_coor,sum_idx,sum_ranges);
-	    Coordinate site(Nd);
-	    for(int d=0;d<Nd;d++){
-	      site[d] = (sum_coor[d] + loop_coor[d] + g_xyz_coor[d] + Slice[d])%latt[d];
-	    }
-	    sum_value+= real(peekSite(*grid_data,site));
-	    if(xyz==0) std::cout << "sum "<<sum_idx<<" "<<sum_value<<std::endl;
-	  }
-	  imageData->SetScalarComponentFromDouble(xyz_coor[0],xyz_coor[1],xyz_coor[2],0,sum_value);
-	}
-	imageData->Modified();
-
-	std::stringstream ss;
-	ss<< files[file] <<"\nSlice "<<Slice << "\nLoop  " <<loop_coor<<"\nSummed "<<sum_dims;
-
-	text->SetInput(ss.str().c_str());
-	
-
-	vtkRenderWindowInteractor* iren = dynamic_cast<vtkRenderWindowInteractor*>(caller);
-	iren->GetRenderWindow()->Render();
-	  
-      }
-
-    if ( vtkCommand::TimerEvent == eventId ) {
-      ++this->TimerCount;
-      std::cout << " This was a timer event count "<<this->TimerCount << std::endl;
-    }
-  
-    
-    if (this->TimerCount >= this->maxCount) {
-      vtkRenderWindowInteractor* iren = dynamic_cast<vtkRenderWindowInteractor*>(caller);
-      if (this->timerId > -1)
-	{
-        iren->DestroyTimer(this->timerId);
-	}
-    }
-  }
-  
-
-private:
-  int TimerCount;
-  int ffile;
-  int xoff;
-  int t;
-public:
-  vtkImageData* imageData = nullptr;
-  vtkTextActor* text = nullptr;
-  vtkFFMPEGWriter *writer = nullptr;
-  int timerId ;
-  int maxCount ;
-  double rms;
-  
-  isosurface * posExtractor;
-  isosurface * negExtractor;
-};
-class SliderCallback : public vtkCommand
-{
-public:
-    static SliderCallback* New()
-    {
-        return new SliderCallback;
-    }
-    virtual void Execute(vtkObject* caller, unsigned long eventId, void* callData)
-    {
-        vtkSliderWidget *sliderWidget = vtkSliderWidget::SafeDownCast(caller);
-        if (sliderWidget)
-        {
-	  contour = ((vtkSliderRepresentation *)sliderWidget->GetRepresentation())->GetValue();
-        }
-	fu->posExtractor->SetValue(0,  SliderCallback::contour*fu->rms);
-	fu->negExtractor->SetValue(0, -SliderCallback::contour*fu->rms);
-	fu->posExtractor->Modified();
-	fu->negExtractor->Modified();
-    }    
-public:
-  static double contour;
-  FrameUpdater * fu;
-};
-
-FrameUpdater * KBfu;
-void KeypressCallbackFunction(vtkObject* caller, long unsigned int eventId,
-                              void* clientData, void* callData)
-{
-  std::cout << "Keypress callback" << std::endl;
-  vtkRenderWindowInteractor* iren = static_cast<vtkRenderWindowInteractor*>(caller);
-  std::cout << "Pressed: " << iren->GetKeySym() << std::endl;
-  //  imageData->Modified();
-}
-
-double SliderCallback::contour;
-
-int main(int argc, char* argv[])
-{
-  using namespace Grid;
-
-  Grid_init(&argc, &argv);
-  GridLogLayout();
-
-  auto latt_size   = GridDefaultLatt();
-  auto simd_layout = GridDefaultSimd(latt_size.size(), vComplex::Nsimd());
-  auto mpi_layout  = GridDefaultMpi();
-  GridCartesian    Grid(latt_size, simd_layout, mpi_layout);
-
-  double default_contour = 1.0;
-  std::string arg;
- 
-  std::cout << argc << " command Line arguments "<<std::endl;
-  for(int c=0;c<argc;c++) {
-    std::cout << " - "<<argv[c]<<std::endl;
-  }
-
-  std::vector<std::string> file_list({
-      "file1",
-      "file2",
-      "file3",
-      "file4",
-      "file5",
-      "file6",
-      "file7",
-      "file8"
-    });
-
-  if( GridCmdOptionExists(argv,argv+argc,"--files") ){
-    arg=GridCmdOptionPayload(argv,argv+argc,"--files");
-    GridCmdOptionCSL(arg, file_list);
-  }
-
-#ifdef MPEG
-  if( GridCmdOptionExists(argv,argv+argc,"--mpeg") ){
-    mpeg = 1;
-  }
-#endif
-
-  if( GridCmdOptionExists(argv,argv+argc,"--fps") ){
-    arg=GridCmdOptionPayload(argv,argv+argc,"--fps");
-    GridCmdOptionInt(arg,framerate);
-  }
-  if( GridCmdOptionExists(argv,argv+argc,"--isosurface") ){
-    arg=GridCmdOptionPayload(argv,argv+argc,"--isosurface");
-    GridCmdOptionFloat(arg,default_contour);
-  }
-  for(int c=0;c<file_list.size();c++) {
-    std::cout << " file: "<<file_list[c]<<std::endl;
-  }
-
-  int NoTime = 0;
-  int Nd; Nd = Grid.Nd();
-  Coordinate    Slice(Nd,0);
-  Coordinate  SumDims(Nd,0);
-  Coordinate LoopDims(Nd,0);
-  Coordinate  XYZDims({0,1,2});
-  if( GridCmdOptionExists(argv,argv+argc,"--slice") ){
-    arg=GridCmdOptionPayload(argv,argv+argc,"--slice");
-    GridCmdOptionIntVector(arg,Slice);
-  }
-  if( GridCmdOptionExists(argv,argv+argc,"--sum") ){
-    arg=GridCmdOptionPayload(argv,argv+argc,"--sum");
-    GridCmdOptionIntVector(arg,SumDims);
-  }
-  if( GridCmdOptionExists(argv,argv+argc,"--loop") ){
-    arg=GridCmdOptionPayload(argv,argv+argc,"--loop");
-    GridCmdOptionIntVector(arg,LoopDims);
-  }
-  if( GridCmdOptionExists(argv,argv+argc,"--xyz") ){
-    arg=GridCmdOptionPayload(argv,argv+argc,"--xyz");
-    GridCmdOptionIntVector(arg,XYZDims);
-    std::cout << "xyz : "<<XYZDims<<std::endl;
-  }
-  if( GridCmdOptionExists(argv,argv+argc,"--notime") ){
-    NoTime = 1;
-    std::cout << "Suppressing time loop"<<std::endl;
-  }
-
-  // Common things:
-  vtkNew<vtkNamedColors> colors;
-  std::array<unsigned char, 4> posColor{{240, 184, 160, 255}};  colors->SetColor("posColor", posColor.data());
-  std::array<unsigned char, 4> bkg{{51, 77, 102, 255}};         colors->SetColor("BkgColor", bkg.data());
-
-  // Create the renderer, the render window, and the interactor. The renderer
-  // draws into the render window, the interactor enables mouse- and
-  // keyboard-based interaction with the data within the render window.
-  //
-  vtkNew<vtkRenderWindow> renWin;
-  vtkNew<vtkRenderWindowInteractor> iren;
-  iren->SetRenderWindow(renWin);
-
-  
-  //  std::vector<LatticeComplexD> data(file_list.size(),&Grid);
-  //  FieldMetaData header;
-
-  
-
-  int frameCount = file_list.size();
-  for(int d=0;d<Grid.Nd();d++) {
-    if ( LoopDims[d] ) frameCount*= latt_size[d];
-
-  }
-  // It is convenient to create an initial view of the data. The FocalPoint
-  // and Position form a vector direction. Later on (ResetCamera() method)
-  // this vector is used to position the camera to look at the data in
-  // this direction.
-  vtkNew<vtkCamera> aCamera;
-  aCamera->SetViewUp(0, 0, -1);
-  aCamera->SetPosition(0, -1000, 0);
-  aCamera->SetFocalPoint(0, 0, 0);
-  aCamera->ComputeViewPlaneNormal();
-  aCamera->Azimuth(30.0);
-  aCamera->Elevation(30.0);
-
-    
-  vtkNew<vtkRenderer> aRenderer;
-  renWin->AddRenderer(aRenderer);
-    
-  double vol = Grid.gSites();
-  std::cout << "Reading "<<file_list[0]<<std::endl;
-  double nrm, nrmbar,rms, contour;
-  {
-    LatticeComplexD data(&Grid);
-    readFile(data,file_list[0]);
-    nrm    = norm2(data);
-  }
-  nrmbar = nrm/vol;
-  rms    = sqrt(nrmbar);
-  
-  contour = default_contour * rms; // default to 1 x RMS
-
-  // The following reader is used to read a series of 2D slices (images)
-  // that compose the volume. The slice dimensions are set, and the
-  // pixel spacing. The data Endianness must also be specified. The reader
-  // uses the FilePrefix in combination with the slice number to construct
-  // filenames using the format FilePrefix.%d. (In this case the FilePrefix
-  // is the root name of the file: quarter.)
-  vtkNew<vtkImageData> imageData;
-  imageData->SetDimensions(latt_size[0],latt_size[1],latt_size[2]);
-  imageData->AllocateScalars(VTK_DOUBLE, 1);
-  for(int xx=0;xx<latt_size[0];xx++){
-    for(int yy=0;yy<latt_size[1];yy++){
-      for(int zz=0;zz<latt_size[2];zz++){
-	Coordinate site({xx,yy,zz,0});
-	RealD value = 0;
-	imageData->SetScalarComponentFromDouble(xx,yy,zz,0,value);
-   }}}
-
-  
-  vtkNew<isosurface> posExtractor;
-  posExtractor->SetInputData(imageData);
-  posExtractor->SetValue(0, contour);
-  
-  vtkNew<vtkStripper> posStripper;
-  posStripper->SetInputConnection(posExtractor->GetOutputPort());
-
-  vtkNew<vtkPolyDataMapper> posMapper;
-  posMapper->SetInputConnection(posStripper->GetOutputPort());
-  posMapper->ScalarVisibilityOff();
-
-  vtkNew<vtkActor> pos;
-  pos->SetMapper(posMapper);
-  pos->GetProperty()->SetDiffuseColor(colors->GetColor3d("posColor").GetData());
-  pos->GetProperty()->SetSpecular(0.3);
-  pos->GetProperty()->SetSpecularPower(20);
-  pos->GetProperty()->SetOpacity(0.5);
-
-  // An isosurface, or contour value is set
-  // The triangle stripper is used to create triangle strips from the
-  // isosurface; these render much faster on may systems.
-  vtkNew<isosurface> negExtractor;
-  negExtractor->SetInputData(imageData);
-  negExtractor->SetValue(0, -contour);
-
-  vtkNew<vtkStripper> negStripper;
-  negStripper->SetInputConnection(negExtractor->GetOutputPort());
-
-  vtkNew<vtkPolyDataMapper> negMapper;
-  negMapper->SetInputConnection(negStripper->GetOutputPort());
-  negMapper->ScalarVisibilityOff();
-    
-  vtkNew<vtkActor> neg;
-  neg->SetMapper(negMapper);
-  neg->GetProperty()->SetDiffuseColor(colors->GetColor3d("Ivory").GetData());
-
-  // An outline provides context around the data.
-  vtkNew<vtkOutlineFilter> outlineData;
-  outlineData->SetInputData(imageData);
-
-  vtkNew<vtkPolyDataMapper> mapOutline;
-  mapOutline->SetInputConnection(outlineData->GetOutputPort());
-  
-  vtkNew<vtkActor> outline;
-  outline->SetMapper(mapOutline);
-  outline->GetProperty()->SetColor(colors->GetColor3d("Black").GetData());
-
-  vtkNew<vtkTextActor> Text;
-  //  Text->SetInput(file_list[f].c_str());
-  Text->SetPosition2(0,0);
-  Text->GetTextProperty()->SetFontSize(24);
-  Text->GetTextProperty()->SetColor(colors->GetColor3d("Gold").GetData());
-
-  vtkNew<vtkTextActor> TextT;
-  TextT->SetInput("T=0");
-  TextT->SetPosition(0,.7*1025);
-  TextT->GetTextProperty()->SetFontSize(24);
-  TextT->GetTextProperty()->SetColor(colors->GetColor3d("Gold").GetData());
-    
-  
-  // Actors are added to the renderer. An initial camera view is created.
-  // The Dolly() method moves the camera towards the FocalPoint,
-  // thereby enlarging the image.
-  //    aRenderer->AddActor(Text);
-  aRenderer->AddActor(TextT);
-  aRenderer->AddActor(outline);
-  aRenderer->AddActor(pos);
-  aRenderer->AddActor(neg);
-
-  // Sign up to receive TimerEvent
-  vtkNew<FrameUpdater> fu;
-
-  fu->SetGrid(&Grid);
-  fu->SetFiles(file_list);
-  fu->SetSlice(Slice);
-  fu->SetSumDimensions (SumDims);
-  fu->SetLoopDimensions(LoopDims);
-  fu->SetDisplayDimensions(XYZDims);
-  fu->SetSliceDimensions();
-  
-
-  fu->imageData = imageData;
-  //    fu->grid_data = &data[f];
-  fu->text      = TextT;
-  fu->maxCount = frameCount;
-  fu->posExtractor = posExtractor;
-  fu->negExtractor = negExtractor;
-  fu->rms = rms;
-      
-  iren->AddObserver(vtkCommand::TimerEvent, fu);
-  iren->AddObserver(vtkCommand::KeyPressEvent, fu);
-
-  aRenderer->SetActiveCamera(aCamera);
-  aRenderer->ResetCamera();
-  aRenderer->SetBackground(colors->GetColor3d("BkgColor").GetData());
-  aCamera->Dolly(1.0);
-
-  //    double nf = file_list.size();
-  //    std::cout << " Adding renderer " <<f<<" of "<<nf<<std::endl;
-  aRenderer->SetViewport(0.0, 0.0,1.0 , 1.0);
-
-  // Note that when camera movement occurs (as it does in the Dolly()
-  // method), the clipping planes often need adjusting. Clipping planes
-  // consist of two planes: near and far along the view direction. The
-  // near plane clips out objects in front of the plane; the far plane
-  // clips out objects behind the plane. This way only what is drawn
-  // between the planes is actually rendered.
-  aRenderer->ResetCameraClippingRange();
-    
-
-  // Set a background color for the renderer and set the size of the
-  // render window (expressed in pixels).
-  // Initialize the event loop and then start it.
-  renWin->SetSize(1024, 1024);
-  renWin->SetWindowName("FieldDensity");
-  renWin->Render();
-
-  // Take a pointer to the FrameUpdater for keypress mgt.
-  //  KBfu = fu;
-  //  vtkNew<vtkCallbackCommand> keypressCallback;
-  //  keypressCallback->SetCallback(KeypressCallbackFunction);
-  //  iren->AddObserver(vtkCommand::KeyPressEvent,keypressCallback);
-
-  iren->Initialize();
-
-  if ( mpeg ) {
-#ifdef MPEG
-    vtkWindowToImageFilter *imageFilter = vtkWindowToImageFilter::New();
-    imageFilter->SetInput( renWin );
-    imageFilter->SetInputBufferTypeToRGB();
-    
-    vtkFFMPEGWriter *writer = vtkFFMPEGWriter::New();
-    writer->SetFileName("movie.avi");
-    writer->SetRate(framerate);
-    writer->SetInputConnection(imageFilter->GetOutputPort());
-    writer->Start();
-
-    for(int i=0;i<fu->maxCount;i++){
-	fu->Execute(iren,vtkCommand::TimerEvent,nullptr);
-	imageFilter->Modified();
-	writer->Write();
-    }
-    writer->End();
-    writer->Delete();
-#else
-    assert(-1 && "MPEG support not compiled");
-#endif
-  } else { 
-  
-    // Add control of contour threshold
-    // Create a slider widget
-    vtkSmartPointer<vtkSliderRepresentation2D> sliderRep = vtkSmartPointer<vtkSliderRepresentation2D>::New();
-    sliderRep->SetMinimumValue(0.0);
-    sliderRep->SetMaximumValue(10.0);
-    sliderRep->SetValue(1.0);
-    sliderRep->SetTitleText("Fraction RMS");
-    // Set color properties:
-
-    // Change the color of the knob that slides
-    //  sliderRep->GetSliderProperty()->SetColor(colors->GetColor3d("Green").GetData());
-    sliderRep->GetTitleProperty()->SetColor(colors->GetColor3d("AliceBlue").GetData());
-    sliderRep->GetLabelProperty()->SetColor(colors->GetColor3d("AliceBlue").GetData());
-    sliderRep->GetSelectedProperty()->SetColor(colors->GetColor3d("DeepPink").GetData());
-
-    // Change the color of the bar
-    sliderRep->GetTubeProperty()->SetColor(colors->GetColor3d("MistyRose").GetData());
-    sliderRep->GetCapProperty()->SetColor(colors->GetColor3d("Yellow").GetData());
-    sliderRep->SetSliderLength(0.05);
-    sliderRep->SetSliderWidth(0.025);
-    sliderRep->SetEndCapLength(0.02);
-
-    sliderRep->GetPoint1Coordinate()->SetCoordinateSystemToNormalizedDisplay();
-    sliderRep->GetPoint1Coordinate()->SetValue(0.1, 0.1);
-    sliderRep->GetPoint2Coordinate()->SetCoordinateSystemToNormalizedDisplay();
-    sliderRep->GetPoint2Coordinate()->SetValue(0.9, 0.1);
-  
-    vtkSmartPointer<vtkSliderWidget> sliderWidget = vtkSmartPointer<vtkSliderWidget>::New();
-    sliderWidget->SetInteractor(iren);
-    sliderWidget->SetRepresentation(sliderRep);
-    sliderWidget->SetAnimationModeToAnimate();
-    sliderWidget->EnabledOn();
-    
-    // Create the slider callback
-    vtkSmartPointer<SliderCallback> slidercallback = vtkSmartPointer<SliderCallback>::New();
-    slidercallback->fu = fu;
-    sliderWidget->AddObserver(vtkCommand::InteractionEvent, slidercallback);
-
-    if ( NoTime==0 ) {
-      int timerId = iren->CreateRepeatingTimer(10000/framerate);
-      std::cout << "timerId "<<timerId<<std::endl;
-    }
-
-    // Start the interaction and timer
-    iren->Start();
-  }
-
-  Grid_finalize();
-
-  return EXIT_SUCCESS;
-}