Boosted fermion attempt

Grid/algorithms/CoarsenedMatrix.h

@@ -158,18 +158,6 @@ public:
     blockPromote(CoarseVec,FineVec,subspace);
   }
 
-  virtual void CreateSubspaceRandom(GridParallelRNG  &RNG) {
-    int nn=nbasis;
-    RealD scale;
-    FineField noise(FineGrid);
-    for(int b=0;b<nn;b++){
-      subspace[b] = Zero();
-      gaussian(RNG,noise);
-      scale = std::pow(norm2(noise),-0.5); 
-      noise=noise*scale;
-      subspace[b] = noise;
-    }
-  }
   virtual void CreateSubspace(GridParallelRNG  &RNG,LinearOperatorBase<FineField> &hermop,int nn=nbasis) {
 
     RealD scale;

Grid/algorithms/GeneralCoarsenedMatrix.h

@@ -1,431 +0,0 @@
-/*************************************************************************************
-
-    Grid physics library, www.github.com/paboyle/Grid 
-
-    Source file: ./lib/algorithms/GeneralCoarsenedMatrix.h
-
-    Copyright (C) 2015
-
-Author: Peter Boyle <pboyle@bnl.gov>
-
-    This program is free software; you can redistribute it and/or modify
-    it under the terms of the GNU General Public License as published by
-    the Free Software Foundation; either version 2 of the License, or
-    (at your option) any later version.
-
-    This program is distributed in the hope that it will be useful,
-    but WITHOUT ANY WARRANTY; without even the implied warranty of
-    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-    GNU General Public License for more details.
-
-    You should have received a copy of the GNU General Public License along
-    with this program; if not, write to the Free Software Foundation, Inc.,
-    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
-
-    See the full license in the file "LICENSE" in the top level distribution directory
-*************************************************************************************/
-/*  END LEGAL */
-#pragma once
-
-#include <Grid/qcd/QCD.h> // needed for Dagger(Yes|No), Inverse(Yes|No)
-
-#include <Grid/lattice/PaddedCell.h>
-#include <Grid/stencil/GeneralLocalStencil.h>
-
-NAMESPACE_BEGIN(Grid);
-
-template<class vobj> void gpermute(vobj & inout,int perm){
-  vobj tmp=inout;
-  if (perm & 0x1 ) { permute(inout,tmp,0); tmp=inout;}
-  if (perm & 0x2 ) { permute(inout,tmp,1); tmp=inout;}
-  if (perm & 0x4 ) { permute(inout,tmp,2); tmp=inout;}
-  if (perm & 0x8 ) { permute(inout,tmp,3); tmp=inout;}
-}
-
-/////////////////////////////////////////////////////////////////
-// Reuse Aggregation class from CoarsenedMatrix for now
-// Might think about *smoothed* Aggregation
-// Equivalent of Geometry class in cartesian case
-/////////////////////////////////////////////////////////////////
-class NonLocalStencilGeometry {
-public:
-  int depth;
-  int npoint;
-  std::vector<Coordinate> shifts;
-  virtual void BuildShifts(void) { assert(0); } ;
-  int Depth(void){return depth;};
-  NonLocalStencilGeometry(int _depth) : depth(_depth)
-  {
-  };
-  virtual ~NonLocalStencilGeometry() {};
-};
-// Need to worry about red-black now
-class NextToNearestStencilGeometry4D : public NonLocalStencilGeometry {
-public:
-  NextToNearestStencilGeometry4D(void) : NonLocalStencilGeometry(2)
-  {
-      this->BuildShifts();
-  };
-  virtual ~NextToNearestStencilGeometry4D() {};
-  virtual void BuildShifts(void)
-  {
-    this->shifts.resize(0);
-    // Like HDCG:  81 point stencil including self connection
-    this->shifts.push_back(Coordinate({0,0,0,0}));
-    // +-x, +-y, +-z, +-t  : 8
-    for(int s=-1;s<=1;s+=2){
-      this->shifts.push_back(Coordinate({s,0,0,0}));
-      this->shifts.push_back(Coordinate({0,s,0,0}));
-      this->shifts.push_back(Coordinate({0,0,s,0}));
-      this->shifts.push_back(Coordinate({0,0,0,s}));
-    }
-    // +-x+-y, +-x+-z, +-x+-t, +-y+-z, +-y+-t, +-z+-t : 24
-    for(int s1=-1;s1<=1;s1+=2){
-    for(int s2=-1;s2<=1;s2+=2){
-      this->shifts.push_back(Coordinate({s1,s2,0,0}));
-      this->shifts.push_back(Coordinate({s1,0,s2,0}));
-      this->shifts.push_back(Coordinate({s1,0,0,s2}));
-      this->shifts.push_back(Coordinate({0,s1,s2,0}));
-      this->shifts.push_back(Coordinate({0,s1,0,s2}));
-      this->shifts.push_back(Coordinate({0,0,s1,s2}));
-    }}
-    this->npoint = this->shifts.size();
-  }
-};
-// Need to worry about red-black now
-class NextToNextToNextToNearestStencilGeometry4D : public NonLocalStencilGeometry {
-public:
-  NextToNextToNextToNearestStencilGeometry4D(void) : NonLocalStencilGeometry(4)
-  {
-    this->BuildShifts();
-  };
-  virtual ~NextToNextToNextToNearestStencilGeometry4D() {}
-  virtual void BuildShifts(void)
-  {
-    this->shifts.resize(0);
-    // Like HDCG:  81 point stencil including self connection
-    this->shifts.push_back(Coordinate({0,0,0,0}));
-    // +-x, +-y, +-z, +-t  : 8
-    for(int s=-1;s<=1;s+=2){
-      this->shifts.push_back(Coordinate({s,0,0,0}));
-      this->shifts.push_back(Coordinate({0,s,0,0}));
-      this->shifts.push_back(Coordinate({0,0,s,0}));
-      this->shifts.push_back(Coordinate({0,0,0,s}));
-    }
-    // +-x+-y, +-x+-z, +-x+-t, +-y+-z, +-y+-t, +-z+-t : 24
-    for(int s1=-1;s1<=1;s1+=2){
-    for(int s2=-1;s2<=1;s2+=2){
-      this->shifts.push_back(Coordinate({s1,s2,0,0}));
-      this->shifts.push_back(Coordinate({s1,0,s2,0}));
-      this->shifts.push_back(Coordinate({s1,0,0,s2}));
-      this->shifts.push_back(Coordinate({0,s1,s2,0}));
-      this->shifts.push_back(Coordinate({0,s1,0,s2}));
-      this->shifts.push_back(Coordinate({0,0,s1,s2}));
-    }}
-    // +-x+-y+-z, +-x+-y+-z, +-x+-y+-z,
-    for(int s1=-1;s1<=1;s1+=2){
-    for(int s2=-1;s2<=1;s2+=2){
-    for(int s3=-1;s3<=1;s3+=2){
-      this->shifts.push_back(Coordinate({s1,s2,s3,0})); // 8x4 = 32
-      this->shifts.push_back(Coordinate({s1,s2,0,s3}));
-      this->shifts.push_back(Coordinate({s1,0,s2,s3}));
-      this->shifts.push_back(Coordinate({0,s1,s2,s3}));
-    }}}
-    for(int s1=-1;s1<=1;s1+=2){
-    for(int s2=-1;s2<=1;s2+=2){
-    for(int s3=-1;s3<=1;s3+=2){
-    for(int s4=-1;s4<=1;s4+=2){
-      this->shifts.push_back(Coordinate({s1,s2,s3,s4})); // 16 
-    }}}}
-    this->npoint = this->shifts.size();
-  }
-};
-class NextToNearestStencilGeometry5D : public NonLocalStencilGeometry {
-public:
-  NextToNearestStencilGeometry5D(void) : NonLocalStencilGeometry(2)
-  {
-      this->BuildShifts();
-  };
-  virtual ~NextToNearestStencilGeometry5D() {};
-  virtual void BuildShifts(void)
-  {
-    this->shifts.resize(0);
-    // Like HDCG:  81 point stencil including self connection
-    this->shifts.push_back(Coordinate({0,0,0,0,0}));
-    // +-x, +-y, +-z, +-t  : 8
-    for(int s=-1;s<=1;s+=2){
-      this->shifts.push_back(Coordinate({0,s,0,0,0}));
-      this->shifts.push_back(Coordinate({0,0,s,0,0}));
-      this->shifts.push_back(Coordinate({0,0,0,s,0}));
-      this->shifts.push_back(Coordinate({0,0,0,0,s}));
-    }
-    // +-x+-y, +-x+-z, +-x+-t, +-y+-z, +-y+-t, +-z+-t : 24
-    for(int s1=-1;s1<=1;s1+=2){
-    for(int s2=-1;s2<=1;s2+=2){
-      this->shifts.push_back(Coordinate({0,s1,s2,0,0}));
-      this->shifts.push_back(Coordinate({0,s1,0,s2,0}));
-      this->shifts.push_back(Coordinate({0,s1,0,0,s2}));
-      this->shifts.push_back(Coordinate({0,0,s1,s2,0}));
-      this->shifts.push_back(Coordinate({0,0,s1,0,s2}));
-      this->shifts.push_back(Coordinate({0,0,0,s1,s2}));
-    }}
-    this->npoint = this->shifts.size();
-  }
-};
-// Need to worry about red-black now
-class NextToNextToNextToNearestStencilGeometry5D : public NonLocalStencilGeometry {
-public:
-  NextToNextToNextToNearestStencilGeometry5D(void) : NonLocalStencilGeometry(4)
-  {
-    this->BuildShifts();
-  };
-  virtual ~NextToNextToNextToNearestStencilGeometry5D() {}
-  virtual void BuildShifts(void)
-  {
-    this->shifts.resize(0);
-    // Like HDCG:  81 point stencil including self connection
-    this->shifts.push_back(Coordinate({0,0,0,0,0}));
-    // +-x, +-y, +-z, +-t  : 8
-    for(int s=-1;s<=1;s+=2){
-      this->shifts.push_back(Coordinate({0,s,0,0,0}));
-      this->shifts.push_back(Coordinate({0,0,s,0,0}));
-      this->shifts.push_back(Coordinate({0,0,0,s,0}));
-      this->shifts.push_back(Coordinate({0,0,0,0,s}));
-    }
-    // +-x+-y, +-x+-z, +-x+-t, +-y+-z, +-y+-t, +-z+-t : 24
-    for(int s1=-1;s1<=1;s1+=2){
-    for(int s2=-1;s2<=1;s2+=2){
-      this->shifts.push_back(Coordinate({0,s1,s2,0,0}));
-      this->shifts.push_back(Coordinate({0,s1,0,s2,0}));
-      this->shifts.push_back(Coordinate({0,s1,0,0,s2}));
-      this->shifts.push_back(Coordinate({0,0,s1,s2,0}));
-      this->shifts.push_back(Coordinate({0,0,s1,0,s2}));
-      this->shifts.push_back(Coordinate({0,0,0,s1,s2}));
-    }}
-    // +-x+-y+-z, +-x+-y+-z, +-x+-y+-z,
-    for(int s1=-1;s1<=1;s1+=2){
-    for(int s2=-1;s2<=1;s2+=2){
-    for(int s3=-1;s3<=1;s3+=2){
-      this->shifts.push_back(Coordinate({0,s1,s2,s3,0})); // 8x4 = 32
-      this->shifts.push_back(Coordinate({0,s1,s2,0,s3}));
-      this->shifts.push_back(Coordinate({0,s1,0,s2,s3}));
-      this->shifts.push_back(Coordinate({0,0,s1,s2,s3}));
-    }}}
-    for(int s1=-1;s1<=1;s1+=2){
-    for(int s2=-1;s2<=1;s2+=2){
-    for(int s3=-1;s3<=1;s3+=2){
-    for(int s4=-1;s4<=1;s4+=2){
-      this->shifts.push_back(Coordinate({0,s1,s2,s3,s4})); // 16 
-    }}}}
-    this->npoint = this->shifts.size();
-  }
-};
-
-// Fine Object == (per site) type of fine field
-// nbasis      == number of deflation vectors
-template<class Fobj,class CComplex,int nbasis>
-class GeneralCoarsenedMatrix : public SparseMatrixBase<Lattice<iVector<CComplex,nbasis > > >  {
-public:
-    
-  typedef iVector<CComplex,nbasis >           siteVector;
-  typedef Lattice<CComplex >                  CoarseComplexField;
-  typedef Lattice<siteVector>                 CoarseVector;
-  typedef Lattice<iMatrix<CComplex,nbasis > > CoarseMatrix;
-  typedef iMatrix<CComplex,nbasis >  Cobj;
-  typedef Lattice< CComplex >   CoarseScalar; // used for inner products on fine field
-  typedef Lattice<Fobj >        FineField;
-  typedef CoarseVector Field;
-  ////////////////////
-  // Data members
-  ////////////////////
-  int hermitian;
-  GridCartesian *       _FineGrid; 
-  GridCartesian *       _CoarseGrid; 
-  NonLocalStencilGeometry &geom;
-  PaddedCell Cell;
-  GeneralLocalStencil Stencil;
-  
-  std::vector<CoarseMatrix> _A;
-  std::vector<CoarseMatrix> _Adag;
-
-  ///////////////////////
-  // Interface
-  ///////////////////////
-  GridCartesian * Grid(void)           { return _FineGrid; };   // this is all the linalg routines need to know
-  GridCartesian * FineGrid(void)       { return _FineGrid; };   // this is all the linalg routines need to know
-  GridCartesian * CoarseGrid(void)     { return _CoarseGrid; };   // this is all the linalg routines need to know
-
-  GeneralCoarsenedMatrix(NonLocalStencilGeometry &_geom,GridCartesian *FineGrid, GridCartesian * CoarseGrid)
-    : geom(_geom),
-      _FineGrid(FineGrid),
-      _CoarseGrid(CoarseGrid),
-      hermitian(1),
-      Cell(_geom.Depth(),_CoarseGrid),
-      Stencil(Cell.grids.back(),geom.shifts)
-  {
-    _A.resize(geom.npoint,CoarseGrid);
-    _Adag.resize(geom.npoint,CoarseGrid);
-  }
-  void M (const CoarseVector &in, CoarseVector &out)
-  {
-    Mult(_A,in,out);
-  }
-  void Mdag (const CoarseVector &in, CoarseVector &out)
-  {
-    Mult(_Adag,in,out);
-  }
-  void Mult (std::vector<CoarseMatrix> &A,const CoarseVector &in, CoarseVector &out)
-  {
-    conformable(CoarseGrid(),in.Grid());
-    conformable(in.Grid(),out.Grid());
-    out.Checkerboard() = in.Checkerboard();
-    CoarseVector tin=in;
-    std::cout << "Calling Exchange"<<std::endl;
-    CoarseVector pin  = Cell.Exchange(tin);
-    //    std::cout << "Called Exchange"<<std::endl;
-    CoarseVector pout(pin.Grid());
-    
-    autoView( in_v , pin, AcceleratorRead);
-    autoView( out_v , pout, AcceleratorWrite);
-    autoView( Stencil_v  , Stencil, AcceleratorRead);
-    int npoint = geom.npoint;
-    typedef LatticeView<Cobj> Aview;
-      
-    Vector<Aview> AcceleratorViewContainer;
-  
-    for(int p=0;p<npoint;p++) AcceleratorViewContainer.push_back(A[p].View(AcceleratorRead));
-    
-    Aview *Aview_p = & AcceleratorViewContainer[0];
-
-    const int Nsimd = CComplex::Nsimd();
-    typedef siteVector calcVector;
-    typedef CComplex   calcComplex;
-
-    int osites=pin.Grid()->oSites();
-
-    for(int point=0;point<npoint;point++){
-      conformable(_A[point],pin);
-    }
-
-    // Should also exchange "A" and "Adag"
-    accelerator_for(sss, osites*nbasis, 1, {
-      int ss = sss/nbasis;
-      int b  = sss%nbasis;
-      assert(ss<osites);
-      calcComplex res;
-      res = Zero();
-      calcVector nbr;
-      int ptype;
-      StencilEntry *SE;
-
-      // FIXME -- exchange the A and the A dag
-      
-      for(int point=0;point<npoint;point++){
-
-	auto SE = Stencil_v.GetEntry(point,ss);
-	
-	int o = SE->_offset;
-
-	// gpermute etc..
-	nbr = in_v[o];
-	assert( o< osites);
-	gpermute(nbr,SE->_permute);
-
-	for(int bb=0;bb<nbasis;bb++) {
-	  res = res + Aview_p[point][ss](b,bb)*nbr(bb);
-	}
-      }
-      out_v[ss](b)=res;
-    });
-
-    for(int p=0;p<geom.npoint;p++) AcceleratorViewContainer[p].ViewClose();
-
-    out = Cell.Extract(pout);
-  };
-
-  void Test(LinearOperatorBase<Lattice<Fobj> > &linop,
-	    Aggregation<Fobj,CComplex,nbasis> & Subspace)
-  {
-    // Create a random
-    GridCartesian *grid = FineGrid();
-    FineField MbV(grid);
-    FineField tmp(grid);
-    FineField f_src(grid);
-    FineField f_res(grid);
-    FineField f_ref(grid);
-    CoarseVector c_src(CoarseGrid());
-    CoarseVector c_res(CoarseGrid());
-    CoarseVector coarseInner(CoarseGrid());
-    GridParallelRNG RNG(CoarseGrid());  RNG.SeedUniqueString(std::string("Coarse RNG"));
-    random(RNG,c_src);
-    blockPromote(c_src,f_src,Subspace.subspace);
-    linop.op(f_src,f_ref);
-    this->Mult (_A,c_src,c_res);
-    blockPromote(c_res,f_res,Subspace.subspace);
-    std::cout << " GeneralCoarsenedMatrix comparison res  "<<norm2(f_res)<<std::endl;
-    std::cout << " GeneralCoarsenedMatrix comparison ref  "<<norm2(f_ref)<<std::endl;
-    f_res = f_res - f_ref;
-    std::cout << " GeneralCoarsenedMatrix comparison diff "<<norm2(f_res)<<std::endl;
-  }
-  void CoarsenOperator(LinearOperatorBase<Lattice<Fobj> > &linop,
-		       Aggregation<Fobj,CComplex,nbasis> & Subspace)
-  {
-    std::cout << GridLogMessage<< "CoarsenMatrix "<< std::endl;
-    GridCartesian *grid = FineGrid();
-    // Orthogonalise the subblocks over the basis
-    CoarseScalar InnerProd(CoarseGrid()); 
-    for(int b=0;b<nbasis;b++){
-      std::cout << "subspace["<<b<<"] " <<norm2(Subspace.subspace[b])<<std::endl;
-    }
-    blockOrthogonalise(InnerProd,Subspace.subspace);
-
-    // Now compute the matrix elements of linop between this orthonormal
-    // set of vectors.
-    FineField bV(grid);
-    FineField MbV(grid);
-    FineField tmp(grid);
-    CoarseVector coarseInner(CoarseGrid());
-
-    // Very inefficient loop of order coarse volume.
-    // First pass hack
-    // Could replace with a coloring scheme our phase scheme
-    // as in BFM
-    for(int bidx=0;bidx<CoarseGrid()->gSites() ;bidx++){
-      Coordinate bcoor;
-      CoarseGrid()->GlobalIndexToGlobalCoor(bidx,bcoor);
-      std::cout << GridLogMessage<< "CoarsenMatrix block "<< bcoor << std::endl;
-      for(int b=0;b<nbasis;b++){
-	blockPick(CoarseGrid(),Subspace.subspace[b],bV,bcoor);
-	linop.HermOp(bV,MbV);
-	blockProject(coarseInner,MbV,Subspace.subspace);
-	for(int p=0;p<geom.npoint;p++){
-	  Coordinate scoor = bcoor;
-	  for(int mu=0;mu<bcoor.size();mu++){
-	    int L = CoarseGrid()->GlobalDimensions()[mu];
-	    scoor[mu] = (bcoor[mu] - geom.shifts[p][mu] + L) % L; // Modulo arithmetic
-	  }
-	  auto ip    = peekSite(coarseInner,scoor);
-	  std::cout << "A["<<b<<"]["<<p<<"]"<<scoor<<" "<<" shift "<<geom.shifts[p]<<" "<< ip <<std::endl;
-	  auto Ab    = peekSite(_A[p],scoor);
-	  auto Adagb = peekSite(_Adag[p],bcoor);
-	  for(int bb=0;bb<nbasis;bb++){
-	    Ab(bb,b)    = ip(bb);
-	    Adagb(b,bb) = conjugate(ip(bb));
-	  }
-	  pokeSite(Ab,_A[p],scoor);
-	  pokeSite(Adagb,_Adag[p],bcoor);
-	}
-      }
-    }
-    std::cout << " Exchanging _A " <<std::endl;
-    for(int p=0;p<geom.npoint;p++){
-      _A[p] = Cell.Exchange(_A[p]);
-      _Adag[p] = Cell.Exchange(_Adag[p]);
-    }
-  }
-  virtual  void Mdiag    (const Field &in, Field &out){ assert(0);};
-  virtual  void Mdir     (const Field &in, Field &out,int dir, int disp){assert(0);};
-  virtual  void MdirAll  (const Field &in, std::vector<Field> &out){assert(0);};
-};
-
-NAMESPACE_END(Grid);

Grid/algorithms/approx/Chebyshev.h

@@ -90,8 +90,9 @@ public:
     order=_order;
       
     if(order < 2) exit(-1);
-    Coeffs.resize(order,0.0);
+    Coeffs.resize(order);
-    Coeffs[order-1] = 1.0;
+    Coeffs.assign(0.,order);
+    Coeffs[order-1] = 1.;
   };
   
   // PB - more efficient low pass drops high modes above the low as 1/x uses all Chebyshev's.

Grid/lattice/PaddedCell.h

@@ -45,9 +45,8 @@ public:
     dims=_grid->Nd();
     AllocateGrids();
     Coordinate local     =unpadded_grid->LocalDimensions();
-    Coordinate procs     =unpadded_grid->ProcessorGrid();
     for(int d=0;d<dims;d++){
-      if ( procs[d] > 1 ) assert(local[d]>=depth);
+      assert(local[d]>=depth);
     }
   }
   void DeleteGrids(void)
@@ -112,7 +111,7 @@ public:
     if(dim==0) conformable(old_grid,unpadded_grid);
     else       conformable(old_grid,grids[dim-1]);
 
-    //    std::cout << " dim "<<dim<<" local "<<local << " padding to "<<plocal<<std::endl;
+    std::cout << " dim "<<dim<<" local "<<local << " padding to "<<plocal<<std::endl;
     // Middle bit
     for(int x=0;x<local[dim];x++){
       InsertSliceLocal(in,padded,x,depth+x,dim);

Grid/qcd/action/fermion/CayleyFermion5D.h

@@ -124,6 +124,11 @@ public:
   RealD                _b;
   RealD                _c;
 
+  // possible boost
+  std::vector<ComplexD> qmu;
+  void set_qmu(std::vector<ComplexD> _qmu) { qmu=_qmu; assert(qmu.size()==Nd);};
+  void addQmu(const FermionField &in, FermionField &out, int dag);
+  
   // Cayley form Moebius (tanh and zolotarev)
   Vector<Coeff_t> omega;
   Vector<Coeff_t> bs;    // S dependent coeffs

Grid/qcd/action/fermion/ContinuedFractionFermion5D.h

@@ -60,6 +60,50 @@ public:
   //      virtual void   Instantiatable(void)=0;
   virtual void   Instantiatable(void) =0;
 
+  void FreePropagator(const FermionField &in,FermionField &out,RealD mass,std::vector<Complex> boundary, std::vector<double> twist)
+  {
+    std::cout << "Free Propagator for PartialFraction"<<std::endl;
+    FermionField in_k(in.Grid());
+    FermionField prop_k(in.Grid());
+    
+    FFT theFFT((GridCartesian *) in.Grid());
+
+    //phase for boundary condition
+    ComplexField coor(in.Grid());
+    ComplexField ph(in.Grid());  ph = Zero();
+    FermionField in_buf(in.Grid()); in_buf = Zero();
+    typedef typename Simd::scalar_type Scalar;
+    Scalar ci(0.0,1.0);
+    assert(twist.size() == Nd);//check that twist is Nd
+    assert(boundary.size() == Nd);//check that boundary conditions is Nd
+    int shift = 0;
+    for(unsigned int nu = 0; nu < Nd; nu++)
+      {
+	// Shift coordinate lattice index by 1 to account for 5th dimension.
+	LatticeCoordinate(coor, nu + shift);
+	double boundary_phase = ::acos(real(boundary[nu]));
+	ph = ph + boundary_phase*coor*((1./(in.Grid()->_fdimensions[nu+shift])));
+	//momenta for propagator shifted by twist+boundary
+	twist[nu] = twist[nu] + boundary_phase/((2.0*M_PI));
+      }
+    in_buf = exp(ci*ph*(-1.0))*in;
+
+    theFFT.FFT_all_dim(in_k,in,FFT::forward);
+    this->MomentumSpacePropagatorHw(prop_k,in_k,mass,twist);
+    theFFT.FFT_all_dim(out,prop_k,FFT::backward);
+    
+    //phase for boundary condition
+    out = out * exp(ci*ph);
+  };
+
+  virtual void FreePropagator(const FermionField &in,FermionField &out,RealD mass) {
+    std::vector<double> twist(Nd,0.0); //default: periodic boundarys in all directions
+    std::vector<Complex> boundary;
+    for(int i=0;i<Nd;i++) boundary.push_back(1);//default: periodic boundary conditions
+    FreePropagator(in,out,mass,boundary,twist);
+  };
+
+  
   // Efficient support for multigrid coarsening
   virtual void  Mdir (const FermionField &in, FermionField &out,int dir,int disp);
   virtual void  MdirAll(const FermionField &in, std::vector<FermionField> &out);

Grid/qcd/action/fermion/PartialFractionFermion5D.h

@@ -39,7 +39,7 @@ class PartialFractionFermion5D : public WilsonFermion5D<Impl>
 public:
   INHERIT_IMPL_TYPES(Impl);
 
-  const int part_frac_chroma_convention=1;
+  const int part_frac_chroma_convention=0;
 
   void   Meooe_internal(const FermionField &in, FermionField &out,int dag);
   void   Mooee_internal(const FermionField &in, FermionField &out,int dag);
@@ -83,12 +83,63 @@ public:
 			   GridRedBlackCartesian &FourDimRedBlackGrid,
 			   RealD _mass,RealD M5,const ImplParams &p= ImplParams());
 
+  PartialFractionFermion5D(GaugeField &_Umu,
+			   GridCartesian         &FiveDimGrid,
+			   GridRedBlackCartesian &FiveDimRedBlackGrid,
+			   GridCartesian         &FourDimGrid,
+			   GridRedBlackCartesian &FourDimRedBlackGrid,
+			   RealD _mass,RealD M5,std::vector<RealD> &_qmu,const ImplParams &p= ImplParams());
+
+  void FreePropagator(const FermionField &in,FermionField &out,RealD mass,std::vector<Complex> boundary, std::vector<double> twist)
+  {
+    std::cout << "Free Propagator for PartialFraction"<<std::endl;
+    FermionField in_k(in.Grid());
+    FermionField prop_k(in.Grid());
+    
+    FFT theFFT((GridCartesian *) in.Grid());
+
+    //phase for boundary condition
+    ComplexField coor(in.Grid());
+    ComplexField ph(in.Grid());  ph = Zero();
+    FermionField in_buf(in.Grid()); in_buf = Zero();
+    typedef typename Simd::scalar_type Scalar;
+    Scalar ci(0.0,1.0);
+    assert(twist.size() == Nd);//check that twist is Nd
+    assert(boundary.size() == Nd);//check that boundary conditions is Nd
+    int shift = 0;
+    for(unsigned int nu = 0; nu < Nd; nu++)
+      {
+	// Shift coordinate lattice index by 1 to account for 5th dimension.
+	LatticeCoordinate(coor, nu + shift);
+	double boundary_phase = ::acos(real(boundary[nu]));
+	ph = ph + boundary_phase*coor*((1./(in.Grid()->_fdimensions[nu+shift])));
+	//momenta for propagator shifted by twist+boundary
+	twist[nu] = twist[nu] + boundary_phase/((2.0*M_PI));
+      }
+    in_buf = exp(ci*ph*(-1.0))*in;
+
+    theFFT.FFT_all_dim(in_k,in,FFT::forward);
+    this->MomentumSpacePropagatorHw(prop_k,in_k,mass,twist);
+    theFFT.FFT_all_dim(out,prop_k,FFT::backward);
+    
+    //phase for boundary condition
+    out = out * exp(ci*ph);
+  };
+
+  virtual void FreePropagator(const FermionField &in,FermionField &out,RealD mass) {
+    std::vector<double> twist(Nd,0.0); //default: periodic boundarys in all directions
+    std::vector<Complex> boundary;
+    for(int i=0;i<Nd;i++) boundary.push_back(1);//default: periodic boundary conditions
+    FreePropagator(in,out,mass,boundary,twist);
+  };
+  
 protected:
 
   virtual void SetCoefficientsTanh(Approx::zolotarev_data *zdata,RealD scale);
   virtual void SetCoefficientsZolotarev(RealD zolo_hi,Approx::zolotarev_data *zdata);
 
   // Part frac
+  std::vector<RealD> qmu;
   RealD mass;
   RealD dw_diag;
   RealD R;

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

@@ -48,7 +48,8 @@ CayleyFermion5D<Impl>::CayleyFermion5D(GaugeField &_Umu,
 			FourDimGrid,
 			FourDimRedBlackGrid,_M5,p),
   mass_plus(_mass), mass_minus(_mass)
-{ 
+{
+  // qmu defaults to zero size;
 }
 
 ///////////////////////////////////////////////////////////////
@@ -270,6 +271,34 @@ void CayleyFermion5D<Impl>::MeooeDag5D    (const FermionField &psi, FermionField
   M5Ddag(psi,psi,Din,lower,diag,upper);
 }
 
+template<class Impl>
+void CayleyFermion5D<Impl>::addQmu(const FermionField &psi,FermionField &chi, int dag)
+{
+  if ( qmu.size() ) {
+
+    Gamma::Algebra Gmu [] = {
+      Gamma::Algebra::GammaX,
+      Gamma::Algebra::GammaY,
+      Gamma::Algebra::GammaZ,
+      Gamma::Algebra::GammaT
+    };
+    std::vector<ComplexD> coeff(Nd);
+    ComplexD ci(0,1);
+
+    assert(qmu.size()==Nd);
+
+    for(int mu=0;mu<Nd;mu++){
+       coeff[mu] = ci*qmu[mu];
+       if ( dag ) coeff[mu] = conjugate(coeff[mu]);
+    }
+
+    chi = chi + Gamma(Gmu[0])*psi*coeff[0];
+    for(int mu=1;mu<Nd;mu++){
+      chi = chi + Gamma(Gmu[mu])*psi*coeff[mu];
+    }
+  }
+}
+
 template<class Impl>
 void CayleyFermion5D<Impl>::M    (const FermionField &psi, FermionField &chi)
 {
@@ -277,8 +306,12 @@ void CayleyFermion5D<Impl>::M    (const FermionField &psi, FermionField &chi)
   
   // Assemble Din
   Meooe5D(psi,Din);
-  
+
   this->DW(Din,chi,DaggerNo);
+
+  // add i q_mu gamma_mu here
+  addQmu(Din,chi,DaggerNo);
+  
   // ((b D_W + D_w hop terms +1) on s-diag
   axpby(chi,1.0,1.0,chi,psi); 
   
@@ -295,6 +328,9 @@ void CayleyFermion5D<Impl>::Mdag (const FermionField &psi, FermionField &chi)
   FermionField Din(psi.Grid());
   // Apply Dw
   this->DW(psi,Din,DaggerYes); 
+
+  // add -i conj(q_mu) gamma_mu here ... if qmu is real, gammm_5 hermitian, otherwise not.
+  addQmu(psi,Din,DaggerYes);
   
   MeooeDag5D(Din,chi);
   

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

@@ -42,13 +42,13 @@ template<class Impl>
 void ContinuedFractionFermion5D<Impl>::SetCoefficientsZolotarev(RealD zolo_hi,Approx::zolotarev_data *zdata)
 {
   // How to check Ls matches??
-  //      std::cout<<GridLogMessage << Ls << " Ls"<<std::endl;
+  std::cout<<GridLogMessage << zdata->n  << " - n"<<std::endl;
-  //      std::cout<<GridLogMessage << zdata->n  << " - n"<<std::endl;
+  std::cout<<GridLogMessage << zdata->da << " -da "<<std::endl;
-  //      std::cout<<GridLogMessage << zdata->da << " -da "<<std::endl;
+  std::cout<<GridLogMessage << zdata->db << " -db"<<std::endl;
-  //      std::cout<<GridLogMessage << zdata->db << " -db"<<std::endl;
+  std::cout<<GridLogMessage << zdata->dn << " -dn"<<std::endl;
-  //      std::cout<<GridLogMessage << zdata->dn << " -dn"<<std::endl;
+  std::cout<<GridLogMessage << zdata->dd << " -dd"<<std::endl;
-  //      std::cout<<GridLogMessage << zdata->dd << " -dd"<<std::endl;
   int Ls = this->Ls;
+  std::cout<<GridLogMessage << Ls << " Ls"<<std::endl;
   assert(zdata->db==Ls);// Beta has Ls coeffs
 
   R=(1+this->mass)/(1-this->mass);
@@ -320,7 +320,7 @@ ContinuedFractionFermion5D<Impl>::ContinuedFractionFermion5D(
       int Ls = this->Ls;
       conformable(solution5d.Grid(),this->FermionGrid());
       conformable(exported4d.Grid(),this->GaugeGrid());
-      ExtractSlice(exported4d, solution5d, Ls-1, Ls-1);
+      ExtractSlice(exported4d, solution5d, Ls-1, 0);
     }
     template<class Impl>
     void ContinuedFractionFermion5D<Impl>::ImportPhysicalFermionSource(const FermionField &input4d,FermionField &imported5d)
@@ -330,7 +330,7 @@ ContinuedFractionFermion5D<Impl>::ContinuedFractionFermion5D(
       conformable(input4d.Grid()   ,this->GaugeGrid());
       FermionField tmp(this->FermionGrid());
       tmp=Zero();
-      InsertSlice(input4d, tmp, Ls-1, Ls-1);
+      InsertSlice(input4d, tmp, Ls-1, 0);
       tmp=Gamma(Gamma::Algebra::Gamma5)*tmp;
       this->Dminus(tmp,imported5d);
     }

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

@@ -255,15 +255,76 @@ void   PartialFractionFermion5D<Impl>::M_internal(const FermionField &psi, Fermi
   }
 	
   {
+    // The 'conventional' Cayley overlap operator is
+    //
+    // Dov = (1+m)/2 + (1-m)/2 g5 sgn Hw
+    //
+    //
+    // With massless limit 1/2(1+g5 sgnHw)
+    //
+    // Luscher shows quite neatly that 1+g5 sgn Hw has tree level propagator i qslash +O(a^2)
+    //
+    // However, the conventional normalisation has both a leading order factor of 2 in Zq
+    // at tree level AND a mass dependent (1-m) that are convenient to absorb.
+    //
+    // In WilsonFermion5DImplementation.h, the tree level propagator for Hw is
+    //
+    // num = -i sin kmu gmu
+    //
+    // denom ( sqrt(sk^2 + (2shk^2 - 1)^2
+    //    b_k = sk2 - M5;
+    //     
+    //    w_k = sqrt(sk + b_k*b_k);
+    //
+    //    denom= ( w_k + b_k + mass*mass) ;
+    //
+    //    denom= one/denom;
+    //    out = num*denom;
+    //
+    // Chroma, and Grid define partial fraction via 4d operator
+    //
+    //   Dpf = 2/(1-m) x Dov = (1+m)/(1-m) + g5 sgn Hw
+    //
+    // Now since:
+    //
+    //      (1+m)/(1-m) = (1-m)/(1-m) + 2m/(1-m) = 1 + 2m/(1-m)
+    //
+    // This corresponds to a modified mass parameter
+    //
+    // It has an annoying 
+    //
+    // 
     double R=(1+this->mass)/(1-this->mass);
     //R g5 psi[Ls] + p[0] H
     ag5xpbg5y_ssp(chi,R*scale,psi,p[nblock]*scale/amax,D,Ls-1,Ls-1);
-	
+    
     for(int b=0;b<nblock;b++){
       int s = 2*b+1;
       double pp = p[nblock-1-b];
       axpby_ssp(chi,1.0,chi,-sqrt(amax*pp)*scale*sign,psi,Ls-1,s);
     }
+
+    if ( qmu.size() ) {
+
+      FermionField qslash_psi(psi.Grid());
+      
+      Gamma::Algebra Gmu [] = {
+			 Gamma::Algebra::GammaX,
+			 Gamma::Algebra::GammaY,
+			 Gamma::Algebra::GammaZ,
+			 Gamma::Algebra::GammaT
+      };
+      ComplexD ci(0,1);
+      assert(qmu.size()==Nd);
+      qslash_psi = Gamma(Gmu[0])*psi;
+      for(int mu=1;mu<Nd;mu++){
+	qslash_psi = Gamma(Gmu[mu])*psi;
+      }
+      //      RealD coeff = 1.0;
+      qslash_psi = Gamma(Gamma::Algebra::Gamma5)*qslash_psi*ci ; // i g5 qslash -- 1-m factor???
+      axpby_ssp(chi,1.0,chi,1.0, qslash_psi,Ls-1,Ls-1);
+    }
+    
   }
 
 }
@@ -411,7 +472,7 @@ void  PartialFractionFermion5D<Impl>::SetCoefficientsZolotarev(RealD zolo_hi,App
       int Ls = this->Ls;
       conformable(solution5d.Grid(),this->FermionGrid());
       conformable(exported4d.Grid(),this->GaugeGrid());
-      ExtractSlice(exported4d, solution5d, Ls-1, Ls-1);
+      ExtractSlice(exported4d, solution5d, Ls-1, 0);
     }
     template<class Impl>
     void PartialFractionFermion5D<Impl>::ImportPhysicalFermionSource(const FermionField &input4d,FermionField &imported5d)
@@ -421,7 +482,8 @@ void  PartialFractionFermion5D<Impl>::SetCoefficientsZolotarev(RealD zolo_hi,App
       conformable(input4d.Grid()   ,this->GaugeGrid());
       FermionField tmp(this->FermionGrid());
       tmp=Zero();
-      InsertSlice(input4d, tmp, Ls-1, Ls-1);
+      std::cout << " importing to slice " << Ls-1 <<std::endl;
+      InsertSlice(input4d, tmp, Ls-1, 0);
       tmp=Gamma(Gamma::Algebra::Gamma5)*tmp;
       this->Dminus(tmp,imported5d);
     }
@@ -442,7 +504,7 @@ PartialFractionFermion5D<Impl>::PartialFractionFermion5D(GaugeField &_Umu,
 
 {
   int Ls = this->Ls;
-
+  qmu.resize(0);
   assert((Ls&0x1)==1); // Odd Ls required
   int nrational=Ls-1;
 
@@ -460,6 +522,22 @@ PartialFractionFermion5D<Impl>::PartialFractionFermion5D(GaugeField &_Umu,
   Approx::zolotarev_free(zdata);
 
 }
+template<class Impl>
+PartialFractionFermion5D<Impl>::PartialFractionFermion5D(GaugeField &_Umu,
+							 GridCartesian         &FiveDimGrid,
+							 GridRedBlackCartesian &FiveDimRedBlackGrid,
+							 GridCartesian         &FourDimGrid,
+							 GridRedBlackCartesian &FourDimRedBlackGrid,
+							 RealD _mass,RealD M5,
+							 std::vector<RealD> &_qmu,
+							 const ImplParams &p)
+  : PartialFractionFermion5D<Impl>(_Umu,
+			     FiveDimGrid,FiveDimRedBlackGrid,
+			     FourDimGrid,FourDimRedBlackGrid,
+			     _mass,M5,p)
+{
+  qmu=_qmu;
+}
 
 NAMESPACE_END(Grid);
 

Grid/stencil/GeneralLocalStencil.h

@@ -46,7 +46,7 @@ class GeneralLocalStencilView {
   accelerator_inline GeneralStencilEntry * GetEntry(int point,int osite) { 
     return & this->_entries_p[point+this->_npoints*osite]; 
   }
-  void ViewClose(void){};
+
 };
 ////////////////////////////////////////
 // The Stencil Class itself
@@ -61,7 +61,7 @@ protected:
 public: 
   GridBase *Grid(void) const { return _grid; }
 
-  View_type View(int mode) const {
+  View_type View(void) const {
     View_type accessor(*( (View_type *) this));
     return accessor;
   }

systems/mac-arm/config-command-mpi

@@ -1,3 +1,4 @@
 BREW=/opt/local/
-CXXFLAGS=-fsanitize=address CXX=g++ ../../configure --enable-simd=NEONv8 --enable-comms=none --enable-unified=yes --prefix $HOME/QCD/GridInstall --with-lime=/Users/peterboyle/QCD/SciDAC/install/ --with-openssl=$BREW --disable-gparity --disable-fermion-reps 
+CXX=mpicxx-openmpi-mp ../../configure --enable-simd=GEN --enable-comms=mpi --enable-unified=yes --prefix $HOME/QCD/GridInstall --with-lime=/Users/peterboyle/QCD/SciDAC/install/ --with-openssl=$BREW --disable-fermion-reps --disable-gparity --disable-debug
+
 

tests/debug/Test_general_coarse.cc

@@ -1,329 +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/GeneralCoarsenedMatrix.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 Field>
-class HermOpAdaptor : public LinearOperatorBase<Field>
-{
-  LinearOperatorBase<Field> & wrapped;
-public:
-  HermOpAdaptor(LinearOperatorBase<Field> &wrapme) : wrapped(wrapme)  {};
-  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){
-    wrapped.HermOp(in,out);
-  }
-  void AdjOp     (const Field &in, Field &out){
-    wrapped.HermOp(in,out);
-  }
-  void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){    assert(0);  }
-  void HermOp(const Field &in, Field &out){
-    wrapped.HermOp(in,out);
-  }
-  
-};
-
-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){
-    Field tmp(in.Grid());
-    _Mat.M(in,tmp);
-    _PV.Mdag(tmp,out);
-  }
-  void AdjOp     (const Field &in, Field &out){
-    Field tmp(in.Grid());
-    _PV.M(tmp,out);
-    _Mat.Mdag(in,tmp);
-  }
-  void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){    assert(0);  }
-  void HermOp(const Field &in, Field &out){
-    std::cout << "HermOp"<<std::endl;
-    Field tmp(in.Grid());
-    _Mat.M(in,tmp);
-    _PV.Mdag(tmp,out);
-    _PV.M(out,tmp);
-    _Mat.Mdag(tmp,out);
-    std::cout << "HermOp done "<<norm2(out)<<std::endl;
-    
-  }
-};
-
-template<class Field> class DumbOperator  : public LinearOperatorBase<Field> {
-public:
-  LatticeComplex scale;
-  DumbOperator(GridBase *grid) : scale(grid)
-  {
-    scale = 0.0;
-    LatticeComplex scalesft(grid);
-    LatticeComplex scaletmp(grid);
-    for(int d=0;d<4;d++){
-      Lattice<iScalar<vInteger> > x(grid); LatticeCoordinate(x,d+1);
-      LatticeCoordinate(scaletmp,d+1);
-      scalesft = Cshift(scaletmp,d+1,1);
-      scale = 100.0*scale + where( mod(x    ,2)==(Integer)0, scalesft,scaletmp);
-    }
-    std::cout << " scale\n" << scale << std::endl;
-  }
-  // Support for coarsening to a multigrid
-  void OpDiag (const Field &in, Field &out) {};
-  void OpDir  (const Field &in, Field &out,int dir,int disp){};
-  void OpDirAll  (const Field &in, std::vector<Field> &out) {};
-
-  void Op     (const Field &in, Field &out){
-    out = scale * in;
-  }
-  void AdjOp  (const Field &in, Field &out){
-    out = scale * in;
-  }
-  void HermOp(const Field &in, Field &out){
-    double n1, n2;
-    HermOpAndNorm(in,out,n1,n2);
-  }
-  void HermOpAndNorm(const Field &in, Field &out,double &n1,double &n2){
-    ComplexD dot;
-
-    out = scale * in;
-
-    dot= innerProduct(in,out);
-    n1=real(dot);
-
-    dot = innerProduct(out,out);
-    n2=real(dot);
-  }
-};
-
-
-int main (int argc, char ** argv)
-{
-  Grid_init(&argc,&argv);
-
-  const int Ls=4;
-
-  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;
-  }
-  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);
-  //SU<Nc>::HotConfiguration(RNG4,Umu);
-  SU<Nc>::ColdConfiguration(Umu);
-  //  auto U = peekLorentz(Umu,0);
-  //  Umu=Zero(); // Make operator local for now
-  //  pokeLorentz(Umu,U,0);
-  
-  RealD mass=0.5;
-  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 ;
-  LatticeFermion prom(FGrid);
-
-  std::vector<LatticeFermion> subspace(nbasis,FGrid);
-
-  std::cout<<GridLogMessage<<"Calling Aggregation class" <<std::endl;
-
-  // Possible tactics -- with zero gauge field, verify block locality of dirac op
-  // Possible tactics -- with zero gauge field, take inner products
-  
-  // Squared operator
-  MdagMLinearOperator<DomainWallFermionD,LatticeFermion> HermDefOp(Ddwf);
-  DumbOperator<LatticeFermion> Diagonal(FGrid);
-  typedef Aggregation<vSpinColourVector,vTComplex,nbasis> Subspace;
-
-  Subspace Aggregates(Coarse5d,FGrid,cb);
-  Aggregates.CreateSubspaceRandom(RNG5);
-
-  std::cout<<GridLogMessage << "Called aggregation class"<< std::endl;
-
-  typedef GeneralCoarsenedMatrix<vSpinColourVector,vTComplex,nbasis> LittleDiracOperator;
-  typedef LittleDiracOperator::CoarseVector CoarseVector;
-
-  NextToNearestStencilGeometry5D geom;
-  LittleDiracOperator LittleDiracOp(geom,FGrid,Coarse5d);
-  LittleDiracOp.CoarsenOperator(HermDefOp,Aggregates);
-  //  LittleDiracOp.CoarsenOperator(Diagonal,Aggregates);
-
-  std::cout<<GridLogMessage<<"Coarsened operator "<<std::endl;
-  
-  CoarseVector c_src (Coarse5d);
-  CoarseVector c_res (Coarse5d);
-  CoarseVector c_proj(Coarse5d);
-
-  subspace=Aggregates.subspace;
-
-  Complex one(1.0);
-  c_src = one;  // 1 in every element for vector 1.
-  Coordinate coor(5,0);
-  
-  std::cout << "c_src"<< c_src<< std::endl;
-  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;
-
-  //  blockPick(Coarse5d,c_src,c_src,coor);
-  //  blockPromote(c_src,prom,subspace);
-
-  //  Diagonal.HermOp(prom,tmp);
-  HermDefOp.HermOp(prom,tmp);
-  //  HermDefOp.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;
-
-  std::cout<<GridLogMessage << "Testing Hermiticity stochastically "<< std::endl;
-  CoarseVector phi(Coarse5d);
-  CoarseVector chi(Coarse5d);
-  CoarseVector Aphi(Coarse5d);
-  CoarseVector Achi(Coarse5d);
-
-  random(CRNG,phi);
-  random(CRNG,chi);
-
-  std::cout<<GridLogMessage<<"Made randoms "<<norm2(phi)<<" " << norm2(chi)<<std::endl;
-
-  LittleDiracOp.M(phi,Aphi);
-  LittleDiracOp.Mdag(chi,Achi);
-  std::cout<<GridLogMessage<<"Aphi "<<norm2(Aphi)<<" Adag chi" << norm2(Achi)<<std::endl;
-
-  ComplexD pAc = innerProduct(chi,Aphi);
-  ComplexD cAp = innerProduct(phi,Achi);
-  ComplexD cAc = innerProduct(chi,Achi);
-  ComplexD pAp = innerProduct(phi,Aphi);
-
-  std::cout<<GridLogMessage<< "pAc "<<pAc<<" cAp "<< cAp<< " diff "<<pAc-adj(cAp)<<std::endl;
-
-  std::cout<<GridLogMessage<< "pAp "<<pAp<<" cAc "<< cAc<<"Should be real"<< std::endl;
-
-  std::cout<<GridLogMessage<<"Testing linearity"<<std::endl;
-  CoarseVector PhiPlusChi(Coarse5d);
-  CoarseVector APhiPlusChi(Coarse5d);
-  CoarseVector linerr(Coarse5d);
-  PhiPlusChi = phi+chi;
-  LittleDiracOp.M(PhiPlusChi,APhiPlusChi);
-
-  linerr= APhiPlusChi-Aphi;
-  linerr= linerr-Achi;
-  std::cout<<GridLogMessage<<"**Diff "<<norm2(linerr)<<std::endl;
-
-  std::cout<<GridLogMessage<<std::endl;
-  std::cout<<GridLogMessage<<std::endl;
-  std::cout<<GridLogMessage<<"*******************************************"<<std::endl;
-  std::cout<<GridLogMessage<<"*******************************************"<<std::endl;
-  std::cout<<GridLogMessage<<"*******************************************"<<std::endl;
-  std::cout<<GridLogMessage<<std::endl;
-  std::cout<<GridLogMessage<<std::endl;
-  
-  PVdagMLinearOperator<DomainWallFermionD,LatticeFermionD> PVdagM(Ddwf,Dpv);
-  HermOpAdaptor<LatticeFermionD> HOA(PVdagM);
-
-  // Run power method on HOA??
-  PowerMethod<LatticeFermion>       PM;   PM(HOA,src);
- 
-  // Warning: This routine calls PVdagM.Op, not PVdagM.HermOp
-  Subspace AggregatesPD(Coarse5d,FGrid,cb);
-  AggregatesPD.CreateSubspaceChebyshev(RNG5,
-				     HOA,
-				     nbasis,
-				     5000.0,
-				     0.02,
-				     100,
-				     50,
-				     50,
-				     0.0);
-  
-  LittleDiracOperator LittleDiracOpPV(geom,FGrid,Coarse5d);
-  LittleDiracOpPV.CoarsenOperator(PVdagM,AggregatesPD);
-
-  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.cc

@@ -1,226 +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/GeneralCoarsenedMatrix.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 Field>
-class HermOpAdaptor : public LinearOperatorBase<Field>
-{
-  LinearOperatorBase<Field> & wrapped;
-public:
-  HermOpAdaptor(LinearOperatorBase<Field> &wrapme) : wrapped(wrapme)  {};
-  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){
-    wrapped.HermOp(in,out);
-  }
-  void AdjOp     (const Field &in, Field &out){
-    wrapped.HermOp(in,out);
-  }
-  void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){    assert(0);  }
-  void HermOp(const Field &in, Field &out){
-    wrapped.HermOp(in,out);
-  }
-  
-};
-
-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){
-    Field tmp(in.Grid());
-    _Mat.M(in,tmp);
-    _PV.Mdag(tmp,out);
-  }
-  void AdjOp     (const Field &in, Field &out){
-    Field tmp(in.Grid());
-    _PV.M(tmp,out);
-    _Mat.Mdag(in,tmp);
-  }
-  void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){    assert(0);  }
-  void HermOp(const Field &in, Field &out){
-    std::cout << "HermOp"<<std::endl;
-    Field tmp(in.Grid());
-    _Mat.M(in,tmp);
-    _PV.Mdag(tmp,out);
-    _PV.M(out,tmp);
-    _Mat.Mdag(tmp,out);
-    std::cout << "HermOp done "<<norm2(out)<<std::endl;
-    
-  }
-};
-
-template<class Field> class DumbOperator  : public LinearOperatorBase<Field> {
-public:
-  LatticeComplex scale;
-  DumbOperator(GridBase *grid) : scale(grid)
-  {
-    scale = 0.0;
-    LatticeComplex scalesft(grid);
-    LatticeComplex scaletmp(grid);
-    for(int d=0;d<4;d++){
-      Lattice<iScalar<vInteger> > x(grid); LatticeCoordinate(x,d+1);
-      LatticeCoordinate(scaletmp,d+1);
-      scalesft = Cshift(scaletmp,d+1,1);
-      scale = 100.0*scale + where( mod(x    ,2)==(Integer)0, scalesft,scaletmp);
-    }
-    std::cout << " scale\n" << scale << std::endl;
-  }
-  // Support for coarsening to a multigrid
-  void OpDiag (const Field &in, Field &out) {};
-  void OpDir  (const Field &in, Field &out,int dir,int disp){};
-  void OpDirAll  (const Field &in, std::vector<Field> &out) {};
-
-  void Op     (const Field &in, Field &out){
-    out = scale * in;
-  }
-  void AdjOp  (const Field &in, Field &out){
-    out = scale * in;
-  }
-  void HermOp(const Field &in, Field &out){
-    double n1, n2;
-    HermOpAndNorm(in,out,n1,n2);
-  }
-  void HermOpAndNorm(const Field &in, Field &out,double &n1,double &n2){
-    ComplexD dot;
-
-    out = scale * in;
-
-    dot= innerProduct(in,out);
-    n1=real(dot);
-
-    dot = innerProduct(out,out);
-    n2=real(dot);
-  }
-};
-
-
-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;
-  }
-  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.5;
-  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 ;
-  LatticeFermion prom(FGrid);
-
-  typedef GeneralCoarsenedMatrix<vSpinColourVector,vTComplex,nbasis> LittleDiracOperator;
-  typedef LittleDiracOperator::CoarseVector CoarseVector;
-
-  NextToNearestStencilGeometry5D geom;
-
-  std::cout<<GridLogMessage<<std::endl;
-  std::cout<<GridLogMessage<<"*******************************************"<<std::endl;
-  std::cout<<GridLogMessage<<std::endl;
-  
-  PVdagMLinearOperator<DomainWallFermionD,LatticeFermionD> PVdagM(Ddwf,Dpv);
-  HermOpAdaptor<LatticeFermionD> HOA(PVdagM);
-
-  // Run power method on HOA??
-  PowerMethod<LatticeFermion>       PM;   PM(HOA,src);
- 
-  // Warning: This routine calls PVdagM.Op, not PVdagM.HermOp
-  typedef Aggregation<vSpinColourVector,vTComplex,nbasis> Subspace;
-  Subspace AggregatesPD(Coarse5d,FGrid,cb);
-  AggregatesPD.CreateSubspaceChebyshev(RNG5,
-				     HOA,
-				     nbasis,
-				     5000.0,
-				     0.02,
-				     100,
-				     50,
-				     50,
-				     0.0);
-  
-  LittleDiracOperator LittleDiracOpPV(geom,FGrid,Coarse5d);
-  LittleDiracOpPV.CoarsenOperator(PVdagM,AggregatesPD);
-
-  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/qdpxx/Test_qdpxx_munprec.cc

@@ -1,7 +1,6 @@
     /*************************************************************************************
 
     Grid physics library, www.github.com/paboyle/Grid 
-
     Source file: ./tests/qdpxx/Test_qdpxx_munprec.cc
 
     Copyright (C) 2015
@@ -26,13 +25,17 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
     See the full license in the file "LICENSE" in the top level distribution directory
     *************************************************************************************/
     /*  END LEGAL */
+#include <chroma.h>
+#include <actions/ferm/invert/syssolver_linop_cg_array.h>
+#include <actions/ferm/invert/syssolver_linop_aggregate.h>
+
 #include <Grid/Grid.h>
 
 int    Ls=8;
 double M5=1.6;
 double mq=0.01;
-double zolo_lo = 0.1;
+double zolo_lo = 0.01;
-double zolo_hi = 2.0;
+double zolo_hi = 7.0;
 double mobius_scale=2.0;
 
 enum ChromaAction {
@@ -55,11 +58,6 @@ enum ChromaAction {
 void calc_grid      (ChromaAction action,Grid::LatticeGaugeField & lat, Grid::LatticeFermion &src, Grid::LatticeFermion &res,int dag);
 void calc_chroma    (ChromaAction action,Grid::LatticeGaugeField & lat, Grid::LatticeFermion &src, Grid::LatticeFermion &res,int dag);
 
-#include <chroma.h>
-#include <actions/ferm/invert/syssolver_linop_cg_array.h>
-#include <actions/ferm/invert/syssolver_linop_aggregate.h>
-
-
 
 namespace Chroma { 
 
@@ -81,7 +79,7 @@ public:
 
     std::vector<int> x(4);
     QDP::multi1d<int> cx(4);
-    std::vector<int> gd= gr.Grid()->GlobalDimensions();
+    Grid::Coordinate gd = gr.Grid()->GlobalDimensions();
 
     for (x[0]=0;x[0]<gd[0];x[0]++){
     for (x[1]=0;x[1]<gd[1];x[1]++){
@@ -124,7 +122,7 @@ public:
 
     std::vector<int> x(5);
     QDP::multi1d<int> cx(4);
-    std::vector<int> gd= gr.Grid()->GlobalDimensions();
+    Grid::Coordinate gd= gr.Grid()->GlobalDimensions();
 
     for (x[0]=0;x[0]<gd[0];x[0]++){
     for (x[1]=0;x[1]<gd[1];x[1]++){
@@ -166,7 +164,7 @@ public:
 
     std::vector<int> x(5);
     QDP::multi1d<int> cx(4);
-    std::vector<int> gd= gr.Grid()->GlobalDimensions();
+    Grid::Coordinate gd= gr.Grid()->GlobalDimensions();
 
     for (x[0]=0;x[0]<gd[0];x[0]++){
     for (x[1]=0;x[1]<gd[1];x[1]++){
@@ -304,7 +302,30 @@ public:
      //     param.approximation_type=COEFF_TYPE_TANH_UNSCALED;
      //     param.approximation_type=COEFF_TYPE_TANH;
      param.tuning_strategy_xml=
-"<TuningStrategy><Name>OVEXT_CONSTANT_STRATEGY</Name></TuningStrategy>\n";
+"<TuningStrategy><Name>OVEXT_CONSTANT_STRATEGY</Name><TuningConstant>1.0</TuningConstant></TuningStrategy>\n";
+     UnprecOvExtFermActArray S_f(cfs,param);
+     Handle< FermState<T4,U,U> > fs( S_f.createState(u) );
+     Handle< LinearOperatorArray<T4> > M(S_f.linOp(fs));
+     return M;
+   }
+   if ( parms == HwPartFracTanh ) {
+     if ( Ls%2 == 0 ) { 
+       printf("Ls is not odd\n");
+       exit(-1);
+     }
+     UnprecOvExtFermActArrayParams param;
+     param.OverMass=M5; 
+     param.Mass=_mq;
+     param.RatPolyDeg = Ls;
+     param.ApproxMin =eps_lo;
+     param.ApproxMax =eps_hi;
+     param.b5 =1.0;
+     param.c5 =1.0;
+     //     param.approximation_type=COEFF_TYPE_ZOLOTAREV;
+     param.approximation_type=COEFF_TYPE_TANH_UNSCALED;
+     //param.approximation_type=COEFF_TYPE_TANH;
+     param.tuning_strategy_xml=
+       "<TuningStrategy><Name>OVEXT_CONSTANT_STRATEGY</Name><TuningConstant>1.0</TuningConstant></TuningStrategy>\n";
      UnprecOvExtFermActArray S_f(cfs,param);
      Handle< FermState<T4,U,U> > fs( S_f.createState(u) );
      Handle< LinearOperatorArray<T4> > M(S_f.linOp(fs));
@@ -316,7 +337,35 @@ public:
      param.ApproxMin=eps_lo;
      param.ApproxMax=eps_hi;
      param.approximation_type=COEFF_TYPE_ZOLOTAREV;
-     param.RatPolyDeg=Ls;
+     param.RatPolyDeg=Ls-1;
+     // The following is why I think Chroma made some directional errors:
+     param.AuxFermAct= std::string(
+"<AuxFermAct>\n"
+"  <FermAct>UNPRECONDITIONED_WILSON</FermAct>\n"
+"  <Mass>-1.8</Mass>\n"
+"  <b5>1</b5>\n"
+"  <c5>0</c5>\n"
+"  <MaxCG>1000</MaxCG>\n"
+"  <RsdCG>1.0e-9</RsdCG>\n"
+"  <FermionBC>\n"
+"      <FermBC>SIMPLE_FERMBC</FermBC>\n"
+"      <boundary>1 1 1 1</boundary>\n"
+"   </FermionBC> \n"
+"</AuxFermAct>"
+);
+     param.AuxFermActGrp= std::string("");
+     UnprecOvlapContFrac5DFermActArray S_f(fbc,param);
+     Handle< FermState<T4,U,U> > fs( S_f.createState(u) );
+     Handle< LinearOperatorArray<T4> > M(S_f.linOp(fs));
+     return  M;
+   }
+   if ( parms == HwContFracTanh ) {
+     UnprecOvlapContFrac5DFermActParams param;
+     param.Mass=_mq; // How is M5 set? Wilson mass In AuxFermAct
+     param.ApproxMin=eps_lo;
+     param.ApproxMax=eps_hi;
+     param.approximation_type=COEFF_TYPE_TANH_UNSCALED;
+     param.RatPolyDeg=Ls-1;
      // The following is why I think Chroma made some directional errors:
      param.AuxFermAct= std::string(
 "<AuxFermAct>\n"
@@ -378,7 +427,14 @@ int main (int argc,char **argv )
    * Setup QDP
    *********************************************************/
   Chroma::initialize(&argc,&argv);
-  Chroma::WilsonTypeFermActs4DEnv::registerAll(); 
+  //  Chroma::WilsonTypeFermActs4DEnv::registerAll(); 
+  Chroma::WilsonTypeFermActsEnv::registerAll(); 
+  //bool linkageHack(void)
+  //{
+  //  bool foo = true;
+  // Inline Measurements
+  //  InlineAggregateEnv::registerAll();
+  //  GaugeInitEnv::registerAll();
 
   /********************************************************
    * Setup Grid
@@ -388,26 +444,34 @@ int main (int argc,char **argv )
                                                                        Grid::GridDefaultSimd(Grid::Nd,Grid::vComplex::Nsimd()),
                                                                        Grid::GridDefaultMpi());
   
-  std::vector<int> gd = UGrid->GlobalDimensions();
+  Grid::Coordinate gd = UGrid->GlobalDimensions();
   QDP::multi1d<int> nrow(QDP::Nd);
   for(int mu=0;mu<4;mu++) nrow[mu] = gd[mu];
 
   QDP::Layout::setLattSize(nrow);
   QDP::Layout::create();
 
-  Grid::GridCartesian         * FGrid   = Grid::SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid);
-  Grid::LatticeGaugeField lat(UGrid);
-  Grid::LatticeFermion    src(FGrid);
-  Grid::LatticeFermion    res_chroma(FGrid);
-  Grid::LatticeFermion    res_grid  (FGrid);
-  
   std::vector<ChromaAction> ActionList({
 		 HtCayleyTanh, // Plain old DWF.
 		 HmCayleyTanh,
 		 HwCayleyTanh,
 		 HtCayleyZolo, // Plain old DWF.
 		 HmCayleyZolo,
-		 HwCayleyZolo
+		 HwCayleyZolo,
+		 HwPartFracZolo,
+		 HwContFracZolo,
+		 HwContFracTanh
+  });
+  std::vector<int> LsList({
+      8,//HtCayleyTanh, // Plain old DWF.
+      8,//HmCayleyTanh,
+      8,//HwCayleyTanh,
+      8,//HtCayleyZolo, // Plain old DWF.
+      8,//HmCayleyZolo,
+      8,//HwCayleyZolo,
+      9,//HwPartFracZolo
+      9, //HwContFracZolo
+      9 //HwContFracTanh
   });
   std::vector<std::string> ActionName({
         "HtCayleyTanh",
@@ -415,10 +479,19 @@ int main (int argc,char **argv )
 	"HwCayleyTanh",
 	"HtCayleyZolo",
 	"HmCayleyZolo",
-        "HwCayleyZolo"
+        "HwCayleyZolo",
+	"HwPartFracZolo",
+	"HwContFracZolo",
+	"HwContFracTanh"
   });
 
   for(int i=0;i<ActionList.size();i++) {
+    Ls = LsList[i];
+    Grid::GridCartesian      * FGrid   = Grid::SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid);
+    Grid::LatticeGaugeField lat(UGrid);
+    Grid::LatticeFermion    src(FGrid);
+    Grid::LatticeFermion    res_chroma(FGrid);
+    Grid::LatticeFermion    res_grid  (FGrid);
     std::cout << "*****************************"<<std::endl;
     std::cout << "Action "<<ActionName[i]<<std::endl;
     std::cout << "*****************************"<<std::endl;
@@ -439,6 +512,7 @@ int main (int argc,char **argv )
       
       std::cout << "Norm of difference "<<Grid::norm2(res_chroma)<<std::endl;
     }
+    delete FGrid;
   }
 
   std::cout << "Finished test "<<std::endl;
@@ -502,7 +576,7 @@ void calc_grid(ChromaAction action,Grid::LatticeGaugeField & Umu, Grid::LatticeF
   Grid::gaussian(RNG5,src);
   Grid::gaussian(RNG5,res);
 
-  Grid::SU<Nc>::HotConfiguration(RNG4,Umu);
+  Grid::SU<Grid::Nc>::HotConfiguration(RNG4,Umu);
 
   /*
   Grid::LatticeColourMatrix U(UGrid);
@@ -519,7 +593,7 @@ void calc_grid(ChromaAction action,Grid::LatticeGaugeField & Umu, Grid::LatticeF
 
   if ( action == HtCayleyTanh ) { 
 
-    Grid::DomainWallFermionR Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,_mass,_M5);
+    Grid::DomainWallFermionD Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,_mass,_M5);
 
     std::cout << Grid::GridLogMessage <<" Calling domain wall multiply "<<std::endl;
 
@@ -535,7 +609,7 @@ void calc_grid(ChromaAction action,Grid::LatticeGaugeField & Umu, Grid::LatticeF
 
     Grid::Real _b = 0.5*(mobius_scale +1.0);
     Grid::Real _c = 0.5*(mobius_scale -1.0);
-    Grid::MobiusZolotarevFermionR D(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,_mass,_M5,_b,_c,zolo_lo,zolo_hi);
+    Grid::MobiusZolotarevFermionD D(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,_mass,_M5,_b,_c,zolo_lo,zolo_hi);
 
     std::cout << Grid::GridLogMessage <<" Calling mobius zolo multiply "<<std::endl;
 
@@ -549,7 +623,7 @@ void calc_grid(ChromaAction action,Grid::LatticeGaugeField & Umu, Grid::LatticeF
 
   if ( action == HtCayleyZolo ) {
 
-    Grid::ShamirZolotarevFermionR D(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,_mass,_M5,zolo_lo,zolo_hi);
+    Grid::ShamirZolotarevFermionD D(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,_mass,_M5,zolo_lo,zolo_hi);
 
     std::cout << Grid::GridLogMessage <<" Calling shamir zolo multiply "<<std::endl;
 
@@ -561,6 +635,60 @@ void calc_grid(ChromaAction action,Grid::LatticeGaugeField & Umu, Grid::LatticeF
     return;
   }
 
+  if ( action == HwPartFracTanh ) {
+
+    Grid::OverlapWilsonPartialFractionTanhFermionD Dov(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,_mass,_M5,1.0);
+
+    std::cout << Grid::GridLogMessage <<" Calling part frac tanh multiply "<<std::endl;
+
+    if ( dag ) 
+      Dov.Mdag(src,res);  
+    else 
+      Dov.M(src,res);  
+
+    return;
+  }
+
+  if ( action == HwContFracTanh ) {
+
+    Grid::OverlapWilsonContFracTanhFermionD Dov(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,_mass,_M5,1.0);
+
+    std::cout << Grid::GridLogMessage <<" Calling cont frac tanh multiply "<<std::endl;
+
+    if ( dag ) 
+      Dov.Mdag(src,res);  
+    else 
+      Dov.M(src,res);  
+
+    return;
+  }
+  if ( action == HwContFracZolo ) {
+
+    Grid::OverlapWilsonContFracZolotarevFermionD Dov(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,_mass,_M5,zolo_lo,zolo_hi);
+
+    std::cout << Grid::GridLogMessage <<" Calling cont frac zolo multiply "<<std::endl;
+
+    if ( dag ) 
+      Dov.Mdag(src,res);  
+    else 
+      Dov.M(src,res);  
+
+    return;
+  }
+
+  if ( action == HwPartFracZolo ) {
+
+    Grid::OverlapWilsonPartialFractionZolotarevFermionD Dov(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,_mass,_M5,zolo_lo,zolo_hi);
+    std::cout << Grid::GridLogMessage <<" Calling part frac zolotarev multiply "<<std::endl;
+
+    if ( dag ) 
+      Dov.Mdag(src,res);  
+    else 
+      Dov.M(src,res);  
+
+    return;
+  }
+  
   /*
   if ( action == HmCayleyTanh ) {
     Grid::Real _b = 0.5*(mobius_scale +1.0);
@@ -581,7 +709,7 @@ void calc_grid(ChromaAction action,Grid::LatticeGaugeField & Umu, Grid::LatticeF
 
   if ( action == HmCayleyTanh ) {
 
-    Grid::ScaledShamirFermionR D(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,_mass,_M5,mobius_scale);
+    Grid::ScaledShamirFermionD D(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,_mass,_M5,mobius_scale);
 
     std::cout << Grid::GridLogMessage <<" Calling scaled shamir multiply "<<std::endl;
 
@@ -595,7 +723,7 @@ void calc_grid(ChromaAction action,Grid::LatticeGaugeField & Umu, Grid::LatticeF
 
   if ( action == HwCayleyTanh ) {
 
-    Grid::OverlapWilsonCayleyTanhFermionR D(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,_mass,_M5,1.0);
+    Grid::OverlapWilsonCayleyTanhFermionD D(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,_mass,_M5,1.0);
 
     if ( dag ) 
       D.Mdag(src,res);  
@@ -607,7 +735,7 @@ void calc_grid(ChromaAction action,Grid::LatticeGaugeField & Umu, Grid::LatticeF
 
   if ( action == HwCayleyZolo ) {
 
-    Grid::OverlapWilsonCayleyZolotarevFermionR D(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,_mass,_M5,zolo_lo,zolo_hi);
+    Grid::OverlapWilsonCayleyZolotarevFermionD D(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,_mass,_M5,zolo_lo,zolo_hi);
 
     if ( dag ) 
       D.Mdag(src,res);  

tests/solver/Test_dwf_cg_prec.cc

@@ -1,4 +1,4 @@
-/*************************************************************************************
+*************************************************************************************
 
 Grid physics library, www.github.com/paboyle/Grid
 
@@ -67,7 +67,13 @@ int main(int argc, char** argv) {
   result = Zero();
   LatticeGaugeField Umu(UGrid);
 
+#if 0
+  FieldMetaData header;
+  std::string file("ckpoint_lat.4000");
+  NerscIO::readConfiguration(Umu,header,file);
+#else  
   SU<Nc>::HotConfiguration(RNG4, Umu);
+#endif
 
   std::cout << GridLogMessage << "Lattice dimensions: " << GridDefaultLatt()
             << "   Ls: " << Ls << std::endl;

tests/solver/Test_dwf_cg_unprec.cc

@@ -54,15 +54,30 @@ int main (int argc, char ** argv)
   GridCartesian         * FGrid   = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid);
   GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid);
 
+  std::vector<ComplexD> qmu;
+  qmu.push_back(ComplexD(0.1,0.0));
+  qmu.push_back(ComplexD(0.0,0.0));
+  qmu.push_back(ComplexD(0.0,0.0));
+  qmu.push_back(ComplexD(0.0,0.01));
+  
+
   std::vector<int> seeds4({1,2,3,4});
   std::vector<int> seeds5({5,6,7,8});
   GridParallelRNG          RNG5(FGrid);  RNG5.SeedFixedIntegers(seeds5);
   GridParallelRNG          RNG4(UGrid);  RNG4.SeedFixedIntegers(seeds4);
 
+  LatticeFermion    tmp(FGrid);
   LatticeFermion    src(FGrid); random(RNG5,src);
   LatticeFermion result(FGrid); result=Zero();
-  LatticeGaugeField Umu(UGrid); SU<Nc>::HotConfiguration(RNG4,Umu);
+  LatticeGaugeField Umu(UGrid); 
-
+#if 0
+  FieldMetaData header;
+  std::string file("ckpoint_lat.4000");
+  NerscIO::readConfiguration(Umu,header,file);
+#else  
+  SU<Nc>::HotConfiguration(RNG4,Umu);
+#endif
+  
   std::vector<LatticeColourMatrix> U(4,UGrid);
   for(int mu=0;mu<Nd;mu++){
     U[mu] = PeekIndex<LorentzIndex>(Umu,mu);
@@ -71,8 +86,15 @@ int main (int argc, char ** argv)
   RealD mass=0.1;
   RealD M5=1.8;
   DomainWallFermionD Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5);
+  Ddwf.qmu = qmu;
 
+  Ddwf.M(src,tmp);
+  std::cout << " |M src|^2 "<<norm2(tmp)<<std::endl;
   MdagMLinearOperator<DomainWallFermionD,LatticeFermion> HermOp(Ddwf);
+  HermOp.HermOp(src,tmp);
+
+  std::cout << " <src|MdagM| src> "<<innerProduct(src,tmp)<<std::endl;
+  
   ConjugateGradient<LatticeFermion> CG(1.0e-6,10000);
   CG(HermOp,src,result);