Feature/fthmc optimise

Grid/Grid.h

@@ -40,11 +40,11 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
 
 #include <Grid/GridCore.h>
 #include <Grid/GridQCDcore.h>
-#include <Grid/qcd/action/Action.h>
+#include <Grid/qcd/action/Action.h>// now contains #include <Grid/parallelIO/MetaData.h>
 #include <Grid/qcd/utils/GaugeFix.h>
 #include <Grid/qcd/utils/CovariantSmearing.h>
 #include <Grid/qcd/smearing/Smearing.h>
-#include <Grid/parallelIO/MetaData.h>
+//#include <Grid/parallelIO/MetaData.h>//moved to ActionCore.h 
 #include <Grid/qcd/hmc/HMC_aggregate.h>
 
 #endif

Grid/GridCore.h

@@ -42,12 +42,12 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
 #include <Grid/DisableWarnings.h>
 #include <Grid/Namespace.h>
 #include <Grid/GridStd.h>
-#include <Grid/threads/Pragmas.h>
 #include <Grid/perfmon/Timer.h>
-//#include <Grid/perfmon/PerfCount.h>
-#include <Grid/util/Util.h>
 #include <Grid/log/Log.h>
 #include <Grid/perfmon/Tracing.h>
+#include <Grid/threads/Pragmas.h>
+//#include <Grid/perfmon/PerfCount.h>
+#include <Grid/util/Util.h>
 #include <Grid/allocator/Allocator.h>
 #include <Grid/simd/Simd.h>
 #include <Grid/threads/ThreadReduction.h>

Grid/Grid_Eigen_Dense.h

@@ -39,6 +39,7 @@
 #endif
 
 /* HIP save and restore compile environment*/
+
 #ifdef GRID_HIP
 #pragma push
 #pragma push_macro("__HIP_DEVICE_COMPILE__")
@@ -63,11 +64,12 @@
 #endif
 
 /*HIP restore*/
+/*
 #ifdef __HIP__REDEFINE__
 #pragma pop_macro("__HIP_DEVICE_COMPILE__")
 #pragma pop
 #endif
-
+*/
 #if defined __GNUC__
 #pragma GCC diagnostic pop
 #endif

Grid/cartesian/Cartesian_base.h

@@ -92,6 +92,7 @@ class GridBase : public CartesianCommunicator , public GridThread {
   ////////////////////////////////////////////////////////////////
   virtual int CheckerBoarded(int dim) =0;
   virtual int CheckerBoard(const Coordinate &site)=0;
+  virtual int CheckerDim(void){ return 0; };
   virtual int CheckerBoardDestination(int source_cb,int shift,int dim)=0;
   virtual int CheckerBoardShift(int source_cb,int dim,int shift,int osite)=0;
   virtual int CheckerBoardShiftForCB(int source_cb,int dim,int shift,int cb)=0;

Grid/cartesian/Cartesian_red_black.h

@@ -60,6 +60,7 @@ class GridRedBlackCartesian : public GridBase
   //  int              _checker_dim;
   std::vector<int> _checker_board;
 
+  virtual int CheckerDim(void){ return _checker_dim; };
   virtual int isCheckerBoarded(void) const { return 1; };
   virtual int CheckerBoarded(int dim){
     if( dim==_checker_dim) return 1;

Grid/cshift/Cshift_common.h

@@ -52,6 +52,7 @@ inline std::pair<int,int> *MapCshiftTable(void)
 template<class vobj> void 
 Gather_plane_simple (const Lattice<vobj> &rhs,deviceVector<vobj> &buffer,int dimension,int plane,int cbmask, int off=0)
 {
+  GRID_TRACE("Gather_plane_simple");
   int rd = rhs.Grid()->_rdimensions[dimension];
 
   if ( !rhs.Grid()->CheckerBoarded(dimension) ) {
@@ -105,6 +106,7 @@ Gather_plane_extract(const Lattice<vobj> &rhs,
 		     ExtractPointerArray<typename vobj::scalar_object> pointers,
 		     int dimension,int plane,int cbmask)
 {
+  GRID_TRACE("Gather_plane_extract");
   int rd = rhs.Grid()->_rdimensions[dimension];
 
   if ( !rhs.Grid()->CheckerBoarded(dimension) ) {
@@ -160,6 +162,7 @@ Gather_plane_extract(const Lattice<vobj> &rhs,
 //////////////////////////////////////////////////////
 template<class vobj> void Scatter_plane_simple (Lattice<vobj> &rhs,deviceVector<vobj> &buffer, int dimension,int plane,int cbmask)
 {
+  GRID_TRACE("Scatter_plane_simple");
   int rd = rhs.Grid()->_rdimensions[dimension];
 
   if ( !rhs.Grid()->CheckerBoarded(dimension) ) {
@@ -214,6 +217,7 @@ template<class vobj> void Scatter_plane_simple (Lattice<vobj> &rhs,deviceVector<
 //////////////////////////////////////////////////////
 template<class vobj> void Scatter_plane_merge(Lattice<vobj> &rhs,ExtractPointerArray<typename vobj::scalar_object> pointers,int dimension,int plane,int cbmask)
 {
+  GRID_TRACE("Scatter_plane_merge");
   int rd = rhs.Grid()->_rdimensions[dimension];
 
   if ( !rhs.Grid()->CheckerBoarded(dimension) ) {
@@ -263,6 +267,7 @@ template<class vobj> void Scatter_plane_merge(Lattice<vobj> &rhs,ExtractPointerA
 
 template<class vobj> void Copy_plane(Lattice<vobj>& lhs,const Lattice<vobj> &rhs, int dimension,int lplane,int rplane,int cbmask)
 {
+  GRID_TRACE("Copy_plane");
   int rd = rhs.Grid()->_rdimensions[dimension];
 
   if ( !rhs.Grid()->CheckerBoarded(dimension) ) {
@@ -311,6 +316,7 @@ template<class vobj> void Copy_plane(Lattice<vobj>& lhs,const Lattice<vobj> &rhs
 
 template<class vobj> void Copy_plane_permute(Lattice<vobj>& lhs,const Lattice<vobj> &rhs, int dimension,int lplane,int rplane,int cbmask,int permute_type)
 {
+  GRID_TRACE("Copy_plane_permute");
   int rd = rhs.Grid()->_rdimensions[dimension];
 
   if ( !rhs.Grid()->CheckerBoarded(dimension) ) {
@@ -358,6 +364,7 @@ template<class vobj> void Copy_plane_permute(Lattice<vobj>& lhs,const Lattice<vo
 //////////////////////////////////////////////////////
 template<class vobj> void Cshift_local(Lattice<vobj>& ret,const Lattice<vobj> &rhs,int dimension,int shift)
 {
+  GRID_TRACE("Cshift_local");
   int sshift[2];
 
   sshift[0] = rhs.Grid()->CheckerBoardShiftForCB(rhs.Checkerboard(),dimension,shift,Even);

Grid/cshift/Cshift_mpi.h

@@ -38,6 +38,7 @@ extern uint64_t checksum_index;
 const int Cshift_verbose=0;
 template<class vobj> Lattice<vobj> Cshift(const Lattice<vobj> &rhs,int dimension,int shift)
 {
+  GRID_TRACE("Cshift");
   typedef typename vobj::vector_type vector_type;
   typedef typename vobj::scalar_type scalar_type;
 
@@ -138,6 +139,7 @@ template<class vobj> void Cshift_simple(Lattice<vobj>& ret,const Lattice<vobj> &
 }
 template<class vobj> void Cshift_comms(Lattice<vobj>& ret,const Lattice<vobj> &rhs,int dimension,int shift)
 {
+  GRID_TRACE("Cshift_comms");
   int sshift[2];
 
   sshift[0] = rhs.Grid()->CheckerBoardShiftForCB(rhs.Checkerboard(),dimension,shift,Even);
@@ -156,6 +158,7 @@ template<class vobj> void Cshift_comms(Lattice<vobj>& ret,const Lattice<vobj> &r
 
 template<class vobj> void Cshift_comms_simd(Lattice<vobj>& ret,const Lattice<vobj> &rhs,int dimension,int shift)
 {
+  GRID_TRACE("Cshift_comms_simd");
   int sshift[2];
 
   sshift[0] = rhs.Grid()->CheckerBoardShiftForCB(rhs.Checkerboard(),dimension,shift,Even);
@@ -173,6 +176,7 @@ template<class vobj> void Cshift_comms_simd(Lattice<vobj>& ret,const Lattice<vob
 }
 template<class vobj> void Cshift_comms(Lattice<vobj> &ret,const Lattice<vobj> &rhs,int dimension,int shift,int cbmask)
 {
+  GRID_TRACE("Cshift_comms_cb");
   typedef typename vobj::vector_type vector_type;
   typedef typename vobj::scalar_type scalar_type;
 
@@ -305,6 +309,7 @@ template<class vobj> void Cshift_comms(Lattice<vobj> &ret,const Lattice<vobj> &r
 
 template<class vobj> void  Cshift_comms_simd(Lattice<vobj> &ret,const Lattice<vobj> &rhs,int dimension,int shift,int cbmask)
 {
+  GRID_TRACE("Cshift_comms_simd_cb");
   GridBase *grid=rhs.Grid();
   const int Nsimd = grid->Nsimd();
   typedef typename vobj::vector_type vector_type;

Grid/lattice/Lattice_trace.h

@@ -97,6 +97,7 @@ Lattice<iScalar<iScalar<iScalar<Vec> > > > Determinant(const Lattice<iScalar<iSc
 template<int N>
 Lattice<iScalar<iScalar<iMatrix<vComplexD, N> > > > Inverse(const Lattice<iScalar<iScalar<iMatrix<vComplexD, N> > > > &Umu)
 {
+#if 0
   GridBase *grid=Umu.Grid();
   auto lvol = grid->lSites();
   Lattice<iScalar<iScalar<iMatrix<vComplexD, N> > > > ret(grid);
@@ -121,9 +122,163 @@ Lattice<iScalar<iScalar<iMatrix<vComplexD, N> > > > Inverse(const Lattice<iScala
       }}
     pokeLocalSite(Ui,ret_v,lcoor);
   });
+#else
+  GridBase *grid=Umu.Grid();
+  auto osites = grid->oSites();
+  const int Nsimd=grid->Nsimd();
+  Lattice<iScalar<iScalar<iMatrix<vComplexD, N> > > > ret(grid);
+  autoView(Umu_v,Umu,CpuRead);
+  autoView(ret_v,ret,CpuWrite);
+  thread_for(site,osites,{
+    Eigen::MatrixXcd EigenU = Eigen::MatrixXcd::Zero(N,N);
+
+    iScalar<iScalar<iMatrix<ComplexD, N> > > Us;
+    iScalar<iScalar<iMatrix<ComplexD, N> > > Ui;
+
+    for(int lane=0;lane<Nsimd;lane++){
+      Us = extractLane(lane,Umu_v[site]);
+      for(int i=0;i<N;i++){
+	for(int j=0;j<N;j++){
+	  EigenU(i,j) = Us()()(i,j);
+	}}
+      Eigen::MatrixXcd EigenUinv = EigenU.inverse();
+      for(int i=0;i<N;i++){
+	for(int j=0;j<N;j++){
+	  Ui()()(i,j) = EigenUinv(i,j);
+	}}
+      insertLane(lane,ret_v[site],Ui);
+    }
+  });
+#endif  
   return ret;
 }
 
+/* Helper functions for inversion of real matrix on GPU based on Nobu's code*/
+template<class type1,class type2,int N>
+accelerator_inline void LUdcmp( iMatrix<type1,N> &LU, iVector<type2,N> &P)
+{
+
+  const RealD TINY=1.0e-40;
+
+  type1 vv[N], tmp, _max;
+  for(int i=0; i<N; i++) vv[i]=0.0;
+
+  for(int i=0; i<N; i++){
+    _max=0.0;
+    for(int j=0; j<N; j++) if( (tmp=abs(LU(i,j))) > _max ) _max = tmp;
+    assert( abs(_max) > TINY );
+    vv[i] = abs(1.0/_max);
+  }
+
+  int imax;
+  for(int k=0; k<N; k++){
+    _max=0.0;
+    for(int i=k; i<N; i++){
+      tmp = vv[i] * abs( LU(i,k) );
+      if(tmp>_max) {
+	_max = tmp;
+	imax = i;
+      }
+    }
+    if(k!=imax){
+      for(int j=0; j<N; j++){
+	tmp = LU(imax,j);
+        LU(imax,j) = LU(k,j);
+	LU(k,j) = tmp;
+      }
+      vv[imax] = vv[k];
+    }
+    P(k)=imax;
+
+    for(int i=k+1; i<N; i++){
+      LU(i,k) = LU(i,k) / LU(k,k);
+      tmp = LU(i,k);
+      for(int j=k+1; j<N; j++) LU(i,j) = LU(i,j) - tmp * LU(k,j);
+    } // end i
+  } // end k
+};
+
+template<class type1,class type2,int N>
+accelerator_inline void solve( iVector<type1,N> &x, const iMatrix<type1,N> LU, const iVector<type2,N> P){
+
+  type1 sum = 0.0;
+
+  int ii=0;
+  for(int i=0; i<N; i++){
+    int ip = P(i);
+    sum = x(ip);
+    x(ip) = x(i);
+    if(ii!=0)for(int j=ii-1;j<i;j++) sum = sum - LU(i,j)*x(j);
+    else if (abs(sum)>0.0) ii=i+1;
+    x(i) = sum;
+  }
+  for(int i=N-1; i>=0; i--){
+    sum = x(i);
+    for(int j=i+1; j<N; j++) sum = sum - LU(i,j)*x(j);
+    x(i) = sum/LU(i,i);
+  }
+};
+
+template<int N>
+Lattice<iScalar<iScalar<iMatrix<vComplexD, N> > > > Inverse_RealPart(const Lattice<iScalar<iScalar<iMatrix<vComplexD, N> > > > &Umu)
+{
+  GridBase *grid=Umu.Grid();
+  auto osites = grid->oSites();
+  const int Nsimd=grid->Nsimd();
+  Lattice<iScalar<iScalar<iMatrix<vComplexD, N> > > > ret(grid);
+#if 0 // CPU version
+  autoView(Umu_v,Umu,CpuRead);
+  autoView(ret_v,ret,CpuWrite);
+  thread_for(site,osites,{
+    Eigen::MatrixXd EigenU = Eigen::MatrixXd::Zero(N,N);
+
+    iScalar<iScalar<iMatrix<ComplexD, N> > > Us;
+    iScalar<iScalar<iMatrix<ComplexD, N> > > Ui;
+
+    for(int lane=0;lane<Nsimd;lane++){
+      Us = extractLane(lane,Umu_v[site]);
+      for(int i=0;i<N;i++){
+        for(int j=0;j<N;j++){
+          EigenU(i,j) = real(Us()()(i,j));
+        }}
+      Eigen::MatrixXd EigenUinv = EigenU.inverse();
+      for(int i=0;i<N;i++){
+        for(int j=0;j<N;j++){
+          Ui()()(i,j) = EigenUinv(i,j);
+        }}
+      insertLane(lane,ret_v[site],Ui);
+    }
+  });
+#else //GPU version
+  autoView(Umu_v,Umu,AcceleratorRead);
+  autoView(ret_v,ret,AcceleratorWrite);
+  accelerator_for(ss,grid->oSites(),vComplex::Nsimd(),{
+      iMatrix<RealD, N>  LU;
+      iVector<Integer, N> P;
+      iVector<RealD, N> e;
+      // scalar layout won't coalesce
+#ifdef GRID_SIMT
+      {
+	int blane=acceleratorSIMTlane(Nsimd); // buffer lane
+#else
+      for(int blane=0;blane<Nsimd;blane++) {
+#endif
+
+	for(int i=0;i<N;i++){
+	  for(int j=0;j<N;j++){
+	    LU(i,j) = getlane(toReal(TensorRemove(Umu_v(ss)()()(i,j))),blane);
+	  }}
+	LUdcmp(LU,P);
+	for(int j=0; j<N; j++){
+	  for(int i=0; i<N; i++) e(i) = (i==j);
+	  solve(e,LU,P);
+	  for(int i=0; i<N; i++) putlane(ret_v[ss]()()(i,j),(ComplexD) e(i),blane);
+	}
+      }
+    });
+#endif
+ return ret;
+}
 
 NAMESPACE_END(Grid);
 #endif