Dirichlet filters running on AMD and now integrated in Fermion op

node boundary or at the single rank boundary

Grid/qcd/action/ActionCore.h

@@ -37,6 +37,9 @@ NAMESPACE_CHECK(ActionSet);
 #include <Grid/qcd/action/ActionParams.h>
 NAMESPACE_CHECK(ActionParams);
 
+#include <Grid/qcd/action/filters/MomentumFilter.h>
+#include <Grid/qcd/action/filters/DirichletFilter.h>
+
 ////////////////////////////////////////////
 // Gauge Actions
 ////////////////////////////////////////////

Grid/qcd/action/fermion/FermionOperator.h

@@ -49,6 +49,8 @@ public:
 
   virtual FermionField &tmp(void) = 0;
 
+  virtual void DirichletBlock(Coordinate & _Block) { assert(0); };
+  
   GridBase * Grid(void)   { return FermionGrid(); };   // this is all the linalg routines need to know
   GridBase * RedBlackGrid(void) { return FermionRedBlackGrid(); };
 

Grid/qcd/action/fermion/WilsonFermion5D.h

@@ -75,6 +75,10 @@ public:
   FermionField _tmp;
   FermionField &tmp(void) { return _tmp; }
 
+  int Dirichlet;
+  Coordinate Block; 
+
+  /********** Deprecate timers **********/
   void Report(void);
   void ZeroCounters(void);
   double DhopCalls;
@@ -174,10 +178,16 @@ public:
 		  GridRedBlackCartesian &FourDimRedBlackGrid,
 		  double _M5,const ImplParams &p= ImplParams());
 
-  void DirichletBlock(std::vector<int> & block){
-    Stencil.DirichletBlock(block); 
-    StencilEven.DirichletBlock(block); 
-    StencilOdd.DirichletBlock(block); 
+  virtual void DirichletBlock(Coordinate & block)
+  {
+    assert(block.size()==Nd+1);
+    if ( block[0] || block[1] || block[2] || block[3] || block[4] ){
+      Dirichlet = 1;
+      Block = block;
+      Stencil.DirichletBlock(block); 
+      StencilEven.DirichletBlock(block); 
+      StencilOdd.DirichletBlock(block);
+    }
   }
   // Constructors
   /*

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

@@ -60,7 +60,8 @@ WilsonFermion5D<Impl>::WilsonFermion5D(GaugeField &_Umu,
   UmuOdd (_FourDimRedBlackGrid),
   Lebesgue(_FourDimGrid),
   LebesgueEvenOdd(_FourDimRedBlackGrid),
-  _tmp(&FiveDimRedBlackGrid)
+  _tmp(&FiveDimRedBlackGrid),
+  Dirichlet(0)
 {
   // some assertions
   assert(FiveDimGrid._ndimension==5);
@@ -218,6 +219,14 @@ void WilsonFermion5D<Impl>::ImportGauge(const GaugeField &_Umu)
 {
   GaugeField HUmu(_Umu.Grid());
   HUmu = _Umu*(-0.5);
+  if ( Dirichlet ) {
+    std::cout << GridLogMessage << " Dirichlet BCs 5d " <<Block<<std::endl;
+    Coordinate GaugeBlock(Nd);
+    for(int d=0;d<Nd;d++) GaugeBlock[d] = Block[d+1];
+    std::cout << GridLogMessage << " Dirichlet BCs 4d " <<GaugeBlock<<std::endl;
+    DirichletFilter<GaugeField> Filter(GaugeBlock);
+    Filter.applyFilter(HUmu);
+  }
   Impl::DoubleStore(GaugeGrid(),Umu,HUmu);
   pickCheckerboard(Even,UmuEven,Umu);
   pickCheckerboard(Odd ,UmuOdd,Umu);

Grid/qcd/action/filters/DirichletFilter.h

@@ -0,0 +1,71 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid
+
+Source file: ./lib/qcd/hmc/integrators/DirichletFilter.h
+
+Copyright (C) 2015
+
+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 */
+//--------------------------------------------------------------------
+#pragma once
+
+NAMESPACE_BEGIN(Grid);
+
+template<typename MomentaField>
+struct DirichletFilter: public MomentumFilterBase<MomentaField>
+{
+  typedef typename MomentaField::vector_type vector_type; //SIMD-vectorized complex type
+  typedef typename MomentaField::scalar_type scalar_type; //scalar complex type
+
+  typedef iScalar<iScalar<iScalar<vector_type> > >            ScalarType; //complex phase for each site
+  
+  Coordinate Block;
+  
+  DirichletFilter(const Coordinate &_Block): Block(_Block){}
+
+  void applyFilter(MomentaField &P) const override
+  {
+    GridBase *grid = P.Grid();
+    typedef decltype(PeekIndex<LorentzIndex>(P, 0)) LatCM;
+    ////////////////////////////////////////////////////
+    // Zero strictly links crossing between domains
+    ////////////////////////////////////////////////////
+    LatticeInteger coor(grid); 
+    LatCM zz(grid); zz = Zero();
+    for(int mu=0;mu<Nd;mu++) {
+      if ( (Block[mu]) && (Block[mu] < grid->GlobalDimensions()[mu] ) ) {
+	// If costly could provide Grid earlier and precompute masks
+	std::cout << " Dirichlet in mu="<<mu<<std::endl;
+	LatticeCoordinate(coor,mu);
+	auto P_mu = PeekIndex<LorentzIndex>(P, mu);
+	P_mu = where(mod(coor,Block[mu])==Integer(Block[mu]-1),zz,P_mu);
+	PokeIndex<LorentzIndex>(P, P_mu, mu);
+      }
+    }
+  }
+};
+
+
+
+NAMESPACE_END(Grid);
+

Grid/qcd/hmc/integrators/Integrator.h

@@ -33,7 +33,6 @@ directory
 #define INTEGRATOR_INCLUDED
 
 #include <memory>
-#include "MomentumFilter.h"
 
 NAMESPACE_BEGIN(Grid);
 

Grid/stencil/Stencil.h

@@ -648,7 +648,7 @@ public:
     }
   }
   /// Introduce a block structure and switch off comms on boundaries
-  void DirichletBlock(const std::vector<int> &dirichlet_block)
+  void DirichletBlock(const Coordinate &dirichlet_block)
   {
     this->_dirichlet = 1;
     for(int ii=0;ii<this->_npoints;ii++){

benchmarks/Benchmark_dwf_fp32.cc

@@ -36,100 +36,6 @@ using namespace Grid;
 /// Move to domains ////
 ////////////////////////
 
-struct DomainDecomposition
-{
-  Coordinate Block;
-
-  DomainDecomposition(const Coordinate &_Block): Block(_Block){ assert(Block.size()==Nd);};
-  
-  template<class Field>
-  void ProjectDomain(Field &f,Integer domain)
-  {
-    GridBase *grid = f.Grid();
-    int dims = grid->Nd();
-    int isDWF= (dims==Nd+1);
-    assert((dims==Nd)||(dims==Nd+1));
-
-    Field   zz(grid);  zz = Zero();
-    LatticeInteger coor(grid);
-    LatticeInteger domaincoor(grid);
-    LatticeInteger mask(grid); mask = Integer(1);
-    LatticeInteger zi(grid);     zi = Integer(0);
-    for(int d=0;d<Nd;d++){
-      Integer B= Block[d];
-      if ( B ) {
-	LatticeCoordinate(coor,d+isDWF);
-	domaincoor = mod(coor,B);
-	mask = where(domaincoor==Integer(0),zi,mask);
-	mask = where(domaincoor==Integer(B-1),zi,mask);
-      }
-    }
-    if ( !domain )
-      f = where(mask==Integer(1),f,zz);
-    else 
-      f = where(mask==Integer(0),f,zz);
-  };
-};
-
-template<typename MomentaField>
-struct DirichletFilter: public MomentumFilterBase<MomentaField>
-{
-  Coordinate Block;
-  
-  DirichletFilter(const Coordinate &_Block): Block(_Block) {}
-
-  // Edge detect using domain projectors
-  void applyFilter (MomentaField &U) const override
-  {
-    DomainDecomposition Domains(Block);
-    GridBase *grid = U.Grid();
-    LatticeInteger  coor(grid);
-    LatticeInteger  face(grid);
-    LatticeInteger  one(grid);   one = 1;
-    LatticeInteger  zero(grid); zero = 0;
-    LatticeInteger  omega(grid);
-    LatticeInteger  omegabar(grid);
-    LatticeInteger  tmp(grid);
-
-    omega=one;    Domains.ProjectDomain(omega,0);
-    omegabar=one; Domains.ProjectDomain(omegabar,1);
-    
-    LatticeInteger nface(grid); nface=Zero();
-    
-    MomentaField projected(grid); projected=Zero();
-    typedef decltype(PeekIndex<LorentzIndex>(U,0)) MomentaLinkField;
-    MomentaLinkField  Umu(grid);
-    MomentaLinkField   zz(grid); zz=Zero();
-
-    int dims = grid->Nd();
-    Coordinate Global=grid->GlobalDimensions();
-    assert(dims==Nd);
-
-    for(int mu=0;mu<Nd;mu++){
-
-      if ( Block[mu]!=0 ) {
-
-	Umu = PeekIndex<LorentzIndex>(U,mu);
-
-	// Upper face 
- 	tmp = Cshift(omegabar,mu,1);
-	tmp = tmp + omega;
-	face = where(tmp == Integer(2),one,zero );
-
- 	tmp = Cshift(omega,mu,1);
-	tmp = tmp + omegabar;
-	face = where(tmp == Integer(2),one,face );
-
-	Umu = where(face,zz,Umu);
-
-	PokeIndex<LorentzIndex>(U, Umu, mu);
-      }
-    }
-  }
-};
-
-
-
 Gamma::Algebra Gmu [] = {
 			 Gamma::Algebra::GammaX,
 			 Gamma::Algebra::GammaY,
@@ -137,7 +43,7 @@ Gamma::Algebra Gmu [] = {
 			 Gamma::Algebra::GammaT
 };
 
-void Benchmark(int Ls, std::vector<int> Dirichlet);
+void Benchmark(int Ls, Coordinate Dirichlet);
 
 int main (int argc, char ** argv)
 {
@@ -152,31 +58,65 @@ int main (int argc, char ** argv)
       std::stringstream ss(argv[i+1]); ss >> Ls;
     }
   }
-  std::vector<int> Dirichlet(5,0);
+
+  //////////////////
+  // With comms
+  //////////////////
+  Coordinate Dirichlet(Nd+1,0);
+
+  std::cout << "\n\n\n\n\n\n" <<std::endl;
+  std::cout << GridLogMessage<< "++++++++++++++++++++++++++++++++++++++++++++++++" <<std::endl;
+  std::cout << GridLogMessage<< " Testing with full communication " <<std::endl;
+  std::cout << GridLogMessage<< "++++++++++++++++++++++++++++++++++++++++++++++++" <<std::endl;
+  
   Benchmark(Ls,Dirichlet);
+
+  //////////////////
+  // Domain decomposed
+  //////////////////
   Coordinate latt4  = GridDefaultLatt();
   Coordinate mpi    = GridDefaultMpi();
+  Coordinate CommDim(Nd);
   Coordinate shm;
   GlobalSharedMemory::GetShmDims(mpi,shm);
-  /*
-  Dirichlet = std::vector<int>({0,
-				latt4[0]/mpi[0] * shm[0],
-				latt4[1]/mpi[1] * shm[1],
-				latt4[2]/mpi[2] * shm[2],
-				latt4[3]/mpi[3] * shm[3]});
-  */
-  Dirichlet = std::vector<int>({0,
-				latt4[0]/mpi[0] ,
-				latt4[1]/mpi[1] ,
-				latt4[2]/mpi[2] ,
-				latt4[3]/mpi[3] });
-  
-  std::cout << " Dirichlet block "<< Dirichlet<< std::endl;
+
+
+  //////////////////////
+  // Node level
+  //////////////////////
+  std::cout << "\n\n\n\n\n\n" <<std::endl;
+  std::cout << GridLogMessage<< "++++++++++++++++++++++++++++++++++++++++++++++++" <<std::endl;
+  std::cout << GridLogMessage<< " Testing without internode communication " <<std::endl;
+  std::cout << GridLogMessage<< "++++++++++++++++++++++++++++++++++++++++++++++++" <<std::endl;
+
+  for(int d=0;d<Nd;d++) CommDim[d]= (mpi[d]/shm[d])>1 ? 1 : 0;
+  Dirichlet[0] = 0;
+  Dirichlet[1] = CommDim[0]*latt4[0]/mpi[0] * shm[0];
+  Dirichlet[2] = CommDim[1]*latt4[1]/mpi[1] * shm[1];
+  Dirichlet[3] = CommDim[2]*latt4[2]/mpi[2] * shm[2];
+  Dirichlet[4] = CommDim[3]*latt4[3]/mpi[3] * shm[3];
+
   Benchmark(Ls,Dirichlet);
+
+  std::cout << "\n\n\n\n\n\n" <<std::endl;
+
+  std::cout << GridLogMessage<< "++++++++++++++++++++++++++++++++++++++++++++++++" <<std::endl;
+  std::cout << GridLogMessage<< " Testing without intranode communication " <<std::endl;
+  std::cout << GridLogMessage<< "++++++++++++++++++++++++++++++++++++++++++++++++" <<std::endl;
+
+  for(int d=0;d<Nd;d++) CommDim[d]= mpi[d]>1 ? 1 : 0;
+  Dirichlet[0] = 0;
+  Dirichlet[1] = CommDim[0]*latt4[0]/mpi[0];
+  Dirichlet[2] = CommDim[1]*latt4[1]/mpi[1];
+  Dirichlet[3] = CommDim[2]*latt4[2]/mpi[2];
+  Dirichlet[4] = CommDim[3]*latt4[3]/mpi[3];
+  
+  Benchmark(Ls,Dirichlet);
+
   Grid_finalize();
   exit(0);
 }
-void Benchmark(int Ls, std::vector<int> Dirichlet)
+void Benchmark(int Ls, Coordinate Dirichlet)
 {
   Coordinate latt4 = GridDefaultLatt();
   GridLogLayout();
@@ -203,8 +143,20 @@ void Benchmark(int Ls, std::vector<int> Dirichlet)
   GridParallelRNG          RNG5(FGrid);  RNG5.SeedUniqueString(std::string("The 5D RNG"));
 
   LatticeFermionF src   (FGrid); random(RNG5,src);
+#if 1
+  src = Zero();
+  {
+    Coordinate origin({0,0,0,latt4[2]-1,0});
+    SpinColourVectorF tmp;
+    tmp=Zero();
+    tmp()(0)(0)=Complex(-2.0,0.0);
+    std::cout << " source site 0 " << tmp<<std::endl;
+    pokeSite(tmp,src,origin);
+  }
+#else
   RealD N2 = 1.0/::sqrt(norm2(src));
   src = src*N2;
+#endif
 
   LatticeFermionF result(FGrid); result=Zero();
   LatticeFermionF    ref(FGrid);    ref=Zero();
@@ -213,29 +165,31 @@ void Benchmark(int Ls, std::vector<int> Dirichlet)
 
   std::cout << GridLogMessage << "Drawing gauge field" << std::endl;
   LatticeGaugeFieldF Umu(UGrid);
+  LatticeGaugeFieldF UmuCopy(UGrid);
   SU<Nc>::HotConfiguration(RNG4,Umu);
+  UmuCopy=Umu;
   std::cout << GridLogMessage << "Random gauge initialised " << std::endl;
 
   ////////////////////////////////////
   // Apply BCs
   ////////////////////////////////////
-  std::cout << GridLogMessage << "Applying BCs " << std::endl;
   Coordinate Block(4);
   for(int d=0;d<4;d++)  Block[d]= Dirichlet[d+1];
 
-  std::cout << GridLogMessage << "Dirichlet Block " << Block<< std::endl;
+  std::cout << GridLogMessage << "Applying BCs for Dirichlet Block5 " << Dirichlet << std::endl;
+  std::cout << GridLogMessage << "Applying BCs for Dirichlet Block4 " << Block << std::endl;
+
   DirichletFilter<LatticeGaugeFieldF> Filter(Block);
   Filter.applyFilter(Umu);
   
   ////////////////////////////////////
   // Naive wilson implementation
   ////////////////////////////////////
-  // replicate across fifth dimension
-  //  LatticeGaugeFieldF Umu5d(FGrid);
   std::vector<LatticeColourMatrixF> U(4,UGrid);
   for(int mu=0;mu<Nd;mu++){
     U[mu] = PeekIndex<LorentzIndex>(Umu,mu);
   }
+
   std::cout << GridLogMessage << "Setting up Cshift based reference " << std::endl;
 
   if (1)
@@ -297,6 +251,8 @@ void Benchmark(int Ls, std::vector<int> Dirichlet)
 
   DomainWallFermionF Dw(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5);
   Dw.DirichletBlock(Dirichlet);
+  Dw.ImportGauge(Umu);
+  
   int ncall =300;
 
   if (1) {
@@ -328,8 +284,8 @@ void Benchmark(int Ls, std::vector<int> Dirichlet)
     std::cout<<GridLogMessage << "mflop/s =   "<< flops/(t1-t0)<<std::endl;
     std::cout<<GridLogMessage << "mflop/s per rank =  "<< flops/(t1-t0)/NP<<std::endl;
     std::cout<<GridLogMessage << "mflop/s per node =  "<< flops/(t1-t0)/NN<<std::endl;
-    std::cout<<GridLogMessage << "RF  GiB/s (base 2) =   "<< 1000000. * data_rf/((t1-t0))<<std::endl;
-    std::cout<<GridLogMessage << "mem GiB/s (base 2) =   "<< 1000000. * data_mem/((t1-t0))<<std::endl;
+    //    std::cout<<GridLogMessage << "RF  GiB/s (base 2) =   "<< 1000000. * data_rf/((t1-t0))<<std::endl;
+    //    std::cout<<GridLogMessage << "mem GiB/s (base 2) =   "<< 1000000. * data_mem/((t1-t0))<<std::endl;
     err = ref-result;
     std::cout<<GridLogMessage << "norm diff   "<< norm2(err)<<std::endl;
 
@@ -382,20 +338,18 @@ void Benchmark(int Ls, std::vector<int> Dirichlet)
     }
     ref = -0.5*ref;
   }
-  //  dump=1;
-  Dw.Dhop(src,result,1);
+
+  Dw.Dhop(src,result,DaggerYes);
+
+  std::cout << GridLogMessage << "----------------------------------------------------------------" << std::endl;
   std::cout << GridLogMessage << "Compare to naive wilson implementation Dag to verify correctness" << std::endl;
+  std::cout << GridLogMessage << "----------------------------------------------------------------" << std::endl;
+
   std::cout<<GridLogMessage << "Called DwDag"<<std::endl;
   std::cout<<GridLogMessage << "norm dag result "<< norm2(result)<<std::endl;
   std::cout<<GridLogMessage << "norm dag ref    "<< norm2(ref)<<std::endl;
   err = ref-result;
   std::cout<<GridLogMessage << "norm dag diff   "<< norm2(err)<<std::endl;
-
-  if (  norm2(err) > 1.0e-4 ) {
-    std::cout << "Error vector is\n" <<err << std::endl;
-    std::cout << "Ref   vector is\n" <<ref << std::endl;
-    std::cout << "Result  vector is\n" <<result << std::endl;
-  }
   assert((norm2(err)<1.0e-4));
   
   LatticeFermionF src_e (FrbGrid);

systems/Spock/sourceme.sh

@@ -1,5 +1,5 @@
 module load PrgEnv-gnu
-module load rocm/4.3.0
+module load rocm/4.5.0
 module load gmp
 module load cray-fftw
 module load craype-accel-amd-gfx908