Command line args and a general clean up

TODO

@@ -5,13 +5,17 @@
 *** Hacks and bug fixes to clean up
 * Had to hack assignment to 1.0 in the tests/Grid_gamma test
 * norm2l is a hack. figure out syntax error and make this norm2 c.f. tests/Grid_gamma.cc
-* Reduce implemention is poor
+
+* Reduce implemention is poor ; need threaded reductions; OMP isn't able to do it for generic objects.
+  
 * Bug in SeedFixedIntegers gives same output on each site.
 * Bug in RNG with complex numbers ; only filling real values; need helper function -- DONE
 * Conformable test in Cshift routines.
 
 *** Functionality
 
+* Implement where to take template scheme.
+
 * Command line args for geometry, simd, etc. layout. Is it necessary to have
   user pass these? Is this a QCD specific?
 

benchmarks/Grid_comms.cc

@@ -8,9 +8,11 @@ int main (int argc, char ** argv)
 {
   Grid_init(&argc,&argv);
 
-  std::vector<int> simd_layout({1,1,2,2});
+  std::vector<int> latt_size;
-  std::vector<int> mpi_layout ({1,2,2,1});
+  std::vector<int> simd_layout;
+  std::vector<int> mpi_layout;
 
+  GridParseLayout(argv,argc,mpi_layout,simd_layout,latt_size);
 
   int Nloop=10;
   int nmu=0;

benchmarks/Grid_memory_bandwidth.cc

@@ -8,8 +8,11 @@ int main (int argc, char ** argv)
 {
   Grid_init(&argc,&argv);
 
-  std::vector<int> simd_layout({1,2,2,2});
+  std::vector<int> tmp_latt_size;
-  std::vector<int> mpi_layout ({1,1,1,1});
+  std::vector<int> simd_layout;
+  std::vector<int> mpi_layout;
+
+  GridParseLayout(argv,argc,mpi_layout,simd_layout,tmp_latt_size);
 
   const int Nvec=8;
   typedef Lattice< iVector< vReal,Nvec> > LatticeVec;

benchmarks/Grid_wilson.cc

@@ -20,9 +20,11 @@ int main (int argc, char ** argv)
 {
   Grid_init(&argc,&argv);
 
-  std::vector<int> simd_layout({1,1,2,2});
+  std::vector<int> latt_size;
-  std::vector<int> mpi_layout ({1,1,1,1});
+  std::vector<int> simd_layout;
-  std::vector<int> latt_size  ({8,8,8,8});
+  std::vector<int> mpi_layout;
+
+  GridParseLayout(argv,argc,mpi_layout,simd_layout,latt_size);
     
   GridCartesian     Grid(latt_size,simd_layout,mpi_layout);
   std::vector<int> seeds({1,2,3,4});

lib/Grid.h

@@ -57,6 +57,7 @@
 
 namespace Grid {
 
+
   void Grid_init(int *argc,char ***argv);
   void Grid_finalize(void);
   void Grid_sa_signal_handler(int sig,siginfo_t *si,void * ptr);
@@ -68,6 +69,13 @@ namespace Grid {
   double usecond(void);
 
 
+  // Common parsing chores
+  std::string GridCmdOptionPayload(char ** begin, char ** end, const std::string & option);
+  bool        GridCmdOptionExists(char** begin, char** end, const std::string& option);
+  void        GridParseIntVector(std::string &str,std::vector<int> & vec);
+  void        GridParseLayout(char **argv,int argc,std::vector<int> &mpi,std::vector<int> &simd,std::vector<int> &latt);
+
+
 };
 
 #endif

lib/Grid_aligned_allocator.h

@@ -20,31 +20,36 @@ public:
   typedef _Tp        value_type;
 
   template<typename _Tp1>  struct rebind { typedef alignedAllocator<_Tp1> other; };
+
   alignedAllocator() throw() { }
+
   alignedAllocator(const alignedAllocator&) throw() { }
+
   template<typename _Tp1> alignedAllocator(const alignedAllocator<_Tp1>&) throw() { }
+
   ~alignedAllocator() throw() { }
+
   pointer       address(reference __x)       const { return &__x; }
   const_pointer address(const_reference __x) const { return &__x; }
+
   size_type  max_size() const throw() { return size_t(-1) / sizeof(_Tp); }
-  // Should override allocate and deallocate
+
   pointer allocate(size_type __n, const void* = 0)
   { 
-    //_Tp * ptr = (_Tp *) memalign(sizeof(_Tp),__n*sizeof(_Tp));
-    // _Tp * ptr = (_Tp *) memalign(128,__n*sizeof(_Tp));
 #ifdef AVX512
     _Tp * ptr = (_Tp *) memalign(128,__n*sizeof(_Tp));
 #else
     _Tp * ptr = (_Tp *) _mm_malloc(__n*sizeof(_Tp),128);
 #endif
-
     return ptr;
   }
+
   void deallocate(pointer __p, size_type) { 
     free(__p); 
   }
   void construct(pointer __p, const _Tp& __val) { };
   void construct(pointer __p) { };
+
   void destroy(pointer __p) { };
 };
 

lib/Grid_comparison.h

@@ -3,6 +3,12 @@
 
 namespace Grid {
 
+  /////////////////////////////////////////
+  // This implementation is a bit poor.
+  // Only support logical operations (== etc)
+  // on scalar objects. Strip any tensor structures.
+  // Should guard this with isGridTensor<> enable if?
+  /////////////////////////////////////////
     // Generic list of functors
     template<class lobj,class robj> class veq {
     public:

lib/Grid_cshift.h

@@ -7,10 +7,6 @@
 #include <cshift/Grid_cshift_none.h>
 #endif
 
-#ifdef GRID_COMMS_FAKE
-#include <cshift/Grid_cshift_fake.h>
-#endif
-
 #ifdef GRID_COMMS_MPI
 #include <cshift/Grid_cshift_mpi.h>
 #endif 

lib/Grid_init.cc

@@ -12,6 +12,7 @@
 #include <signal.h>
 #include <iostream>
 #include <Grid.h>
+#include <algorithm>
 
 #undef __X86_64
 #define MAC
@@ -22,21 +23,6 @@
 
 namespace Grid {
 
-  std::streambuf *Grid_saved_stream_buf;
-#if 0
-  void Grid_quiesce_nodes(void)
-  {
-#ifdef GRID_COMMS_MPI
-    int me;
-    MPI_Comm_rank(MPI_COMM_WORLD,&me);
-    std::streambuf* Grid_saved_stream_buf = std::cout.rdbuf();
-    if ( me ) { 
-      std::ofstream file("log.node");
-      std::cout.rdbuf(file.rdbuf());
-    }
-#endif
-  }
-#endif
   void Grid_quiesce_nodes(void)
   {
 #ifdef GRID_COMMS_MPI
@@ -54,18 +40,95 @@ namespace Grid {
 #endif
   }
 
+std::string GridCmdOptionPayload(char ** begin, char ** end, const std::string & option)
+{
+  char ** itr = std::find(begin, end, option);
+  if (itr != end && ++itr != end) {
+    std::string payload(*itr);
+    return payload;
+  }
+  return std::string("");
+}
+bool GridCmdOptionExists(char** begin, char** end, const std::string& option)
+{
+  return std::find(begin, end, option) != end;
+}
 void Grid_init(int *argc,char ***argv)
 {
 #ifdef GRID_COMMS_MPI
   MPI_Init(argc,argv);
 #endif
-  Grid_debug_handler_init();
+  // Parse command line args.
   Grid_quiesce_nodes();
+
 }
+
+void GridCmdOptionIntVector(std::string &str,std::vector<int> & vec)
+{
+  vec.resize(0);
+  std::stringstream ss(str);
+  int i;
+  while (ss >> i){
+    vec.push_back(i);
+    if (ss.peek() == ',')
+      ss.ignore();
+  }    
+  return;
+}
+
+void GridParseLayout(char **argv,int argc,std::vector<int> &mpi,std::vector<int> &simd,std::vector<int> &latt)
+{
+  mpi =std::vector<int>({1,1,1,1});
+#if defined(AVX) || defined (AVX2)
+  simd=std::vector<int>({1,1,2,2});
+#endif
+#if defined(SSE4)
+  simd=std::vector<int>({1,1,1,2});
+#endif
+#if defined(AVX512)
+  simd=std::vector<int>({1,2,2,2});
+#endif
+  latt=std::vector<int>({8,8,8,8});
+  
+  std::string arg;
+  if( GridCmdOptionExists(argv,argv+argc,"--mpi-layout") ){
+    arg = GridCmdOptionPayload(argv,argv+argc,"--mpi-layout");
+    GridCmdOptionIntVector(arg,mpi);
+  }
+  if( GridCmdOptionExists(argv,argv+argc,"--simd-layout") ){
+    arg= GridCmdOptionPayload(argv,argv+argc,"--simd-layout");
+    GridCmdOptionIntVector(arg,simd);
+  }
+  if( GridCmdOptionExists(argv,argv+argc,"--lattice") ){
+    arg= GridCmdOptionPayload(argv,argv+argc,"--lattice");
+    GridCmdOptionIntVector(arg,latt);
+  }
+  std::cout<<"MPI layout";
+  for(int i=0;i<mpi.size();i++){
+    std::cout<<mpi[i]<<" ";
+  }
+  std::cout<<std::endl;
+  
+  std::cout<<"SIMD layout";
+  for(int i=0;i<simd.size();i++){
+    std::cout<<simd[i]<<" ";
+  }
+  std::cout<<std::endl;
+  
+  std::cout<<"Grid ";
+  for(int i=0;i<latt.size();i++){
+    std::cout<<latt[i]<<" ";
+  }
+  std::cout<<std::endl;
+  
+}
+
+
 void Grid_finalize(void)
 {
 #ifdef GRID_COMMS_MPI
   MPI_Finalize();
+  Grid_unquiesce_nodes();
 #endif
 }
 double usecond(void) {
@@ -74,7 +137,6 @@ double usecond(void) {
   return 1.0*tv.tv_usec + 1.0e6*tv.tv_sec;
 }
 
-
 #define _NBACKTRACE (256)
 void * Grid_backtrace_buffer[_NBACKTRACE];
 

lib/Makefile.am

@@ -7,7 +7,7 @@ if BUILD_COMMS_MPI
 endif
 
 if BUILD_COMMS_NONE
-  extra_sources+=communicator/Grid_communicator_fake.cc
+  extra_sources+=communicator/Grid_communicator_none.cc
 endif
 
 #

lib/communicator/Grid_communicator_none.cc

@@ -29,7 +29,7 @@ void CartesianCommunicator::SendToRecvFrom(void *xmit,
 					   int from,
 					   int bytes)
 {
-  exit(-1);
+  assert(0);
 }
 void CartesianCommunicator::SendToRecvFromBegin(std::vector<CommsRequest_t> &list,
 						void *xmit,
@@ -38,11 +38,11 @@ void CartesianCommunicator::SendToRecvFromBegin(std::vector<CommsRequest_t> &lis
 						int from,
 						int bytes)
 {
-  exit(-1);
+  assert(0);
 }
 void CartesianCommunicator::SendToRecvFromComplete(std::vector<CommsRequest_t> &list)
 {
-  exit(-1);
+  assert(0);
 }
 
 void CartesianCommunicator::Barrier(void)

lib/math/Grid_math_arith_sub.h

@@ -68,24 +68,6 @@ template<class vtype,class ltype,class rtype, int N> inline void sub(iMatrix<vty
     return;
 }
 
-template<class v> void vprefetch(const iScalar<v> &vv)
-{
-  vprefetch(vv._internal);
-}
-template<class v,int N> void vprefetch(const iVector<v,N> &vv)
-{
-  for(int i=0;i<N;i++){
-    vprefetch(vv._internal[i]);
-  }
-}
-template<class v,int N> void vprefetch(const iMatrix<v,N> &vv)
-{
-  for(int i=0;i<N;i++){
-  for(int j=0;j<N;j++){
-    vprefetch(vv._internal[i][j]);
-  }}
-}
-
     // - operator for scalar, vector, matrix
 template<class ltype,class rtype> inline auto
 operator - (const iScalar<ltype>& lhs, const iScalar<rtype>& rhs) -> iScalar<decltype(lhs._internal - rhs._internal)>

lib/math/Grid_math_tensors.h

@@ -16,9 +16,9 @@ namespace Grid {
 // However note that doing this eliminates some syntactical sugar such as 
 // calling the constructor explicitly or implicitly
 //
-#undef TENSOR_IS_POD
+class GridTensorBase {};
 
-template<class vtype> class iScalar
+template<class vtype> class iScalar :public GridTensorBase
 {
 public:
   vtype _internal;
@@ -34,12 +34,9 @@ public:
 
   // Scalar no action
   //  template<int Level> using tensor_reduce_level = typename iScalar<GridTypeMapper<vtype>::tensor_reduce_level<Level> >;
-
-#ifndef TENSOR_IS_POD
   iScalar()=default;
   iScalar(scalar_type s) : _internal(s) {};// recurse down and hit the constructor for vector_type
   iScalar(const Zero &z){ *this = zero; };
-#endif
 
     iScalar<vtype> & operator= (const Zero &hero){
       zeroit(*this);
@@ -87,15 +84,12 @@ public:
     inline const vtype & operator ()(void) const {
       return _internal;
     }
-    //    inline vtype && operator ()(void) {
-    //      return _internal;
-    //    }
 
     operator ComplexD () const { return(TensorRemove(_internal)); };
     operator RealD () const { return(real(TensorRemove(_internal))); }
 
     // convert from a something to a scalar
-    template<class T,typename std::enable_if<isGridTensor<T>::notvalue, T>::type* = nullptr > inline auto operator = (T arg) -> iScalar<vtype>
+    template<class T,typename std::enable_if<!isGridTensor<T>::value, T>::type* = nullptr > inline auto operator = (T arg) -> iScalar<vtype>
     { 
       _internal = vtype(arg);
       return *this;
@@ -105,13 +99,13 @@ public:
 ///////////////////////////////////////////////////////////
 // Allows to turn scalar<scalar<scalar<double>>>> back to double.
 ///////////////////////////////////////////////////////////
-template<class T>     inline typename std::enable_if<isGridTensor<T>::notvalue, T>::type TensorRemove(T arg) { return arg;}
+template<class T>     inline typename std::enable_if<!isGridTensor<T>::value, T>::type TensorRemove(T arg) { return arg;}
 template<class vtype> inline auto TensorRemove(iScalar<vtype> arg) -> decltype(TensorRemove(arg._internal))
 {
   return TensorRemove(arg._internal);
 }
     
-template<class vtype,int N> class iVector
+template<class vtype,int N> class iVector :public GridTensorBase
 {
 public:
   vtype _internal[N];
@@ -125,11 +119,8 @@ public:
 
 
   enum { TensorLevel = GridTypeMapper<vtype>::TensorLevel + 1};
-
-#ifndef TENSOR_IS_POD
   iVector(const Zero &z){ *this = zero; };
   iVector() =default;
-#endif
 
     iVector<vtype,N> & operator= (const Zero &hero){
         zeroit(*this);
@@ -184,7 +175,7 @@ public:
     //    }
 };
     
-template<class vtype,int N> class iMatrix
+template<class vtype,int N> class iMatrix :public GridTensorBase
 {
 public:
   vtype _internal[N][N];
@@ -198,16 +189,16 @@ public:
 
   enum { TensorLevel = GridTypeMapper<vtype>::TensorLevel + 1};
 
-#ifndef TENSOR_IS_POD
+
   iMatrix(const Zero &z){ *this = zero; };
   iMatrix() =default;
-#endif
+
 
   iMatrix<vtype,N> & operator= (const Zero &hero){
     zeroit(*this);
     return *this;
   }
-  template<class T,typename std::enable_if<isGridTensor<T>::notvalue, T>::type* = nullptr > inline auto operator = (T arg) -> iMatrix<vtype,N>
+  template<class T,typename std::enable_if<!isGridTensor<T>::value, T>::type* = nullptr > inline auto operator = (T arg) -> iMatrix<vtype,N>
     { 
       zeroit(*this);
       for(int i=0;i<N;i++)
@@ -278,6 +269,23 @@ public:
 
 };
 
+template<class v> void vprefetch(const iScalar<v> &vv)
+{
+  vprefetch(vv._internal);
+}
+template<class v,int N> void vprefetch(const iVector<v,N> &vv)
+{
+  for(int i=0;i<N;i++){
+    vprefetch(vv._internal[i]);
+  }
+}
+template<class v,int N> void vprefetch(const iMatrix<v,N> &vv)
+{
+  for(int i=0;i<N;i++){
+  for(int j=0;j<N;j++){
+    vprefetch(vv._internal[i][j]);
+  }}
+}
 
 
 }

lib/qcd/Grid_qcd_wilson_dop.cc

@@ -103,6 +103,7 @@ void WilsonMatrix::Dhop(const LatticeFermion &in, LatticeFermion &out)
   vHalfSpinColourVector *chi_p;
   int offset,local,perm, ptype;
 
+#pragma omp parallel for
   for(int sss=0;sss<grid->oSites();sss++){
 
     int ss = sss;

tests/Grid_gamma.cc

@@ -14,9 +14,11 @@ int main (int argc, char ** argv)
 {
   Grid_init(&argc,&argv);
 
-  std::vector<int> simd_layout({1,1,2,2});
+  std::vector<int> latt_size;
-  std::vector<int> mpi_layout ({1,1,1,1});
+  std::vector<int> simd_layout;
-  std::vector<int> latt_size  ({8,8,8,8});
+  std::vector<int> mpi_layout;
+
+  GridParseLayout(argv,argc,mpi_layout,simd_layout,latt_size);
     
   GridCartesian     Grid(latt_size,simd_layout,mpi_layout);
 

tests/Grid_main.cc

@@ -25,15 +25,12 @@ int main (int argc, char ** argv)
 {
   Grid_init(&argc,&argv);
 
-  std::vector<int> latt_size(4);
+  std::vector<int> latt_size;
+  std::vector<int> simd_layout;
+  std::vector<int> mpi_layout;
 
-  std::vector<int> simd_layout(4);
+  GridParseLayout(argv,argc,mpi_layout,simd_layout,latt_size);
-  
+  latt_size.resize(4);
-  std::vector<int> mpi_layout(4);
-  mpi_layout[0]=1;
-  mpi_layout[1]=1;
-  mpi_layout[2]=1;
-  mpi_layout[3]=1;
 
 #ifdef AVX512
  for(int omp=128;omp<236;omp+=16){
@@ -52,25 +49,6 @@ int main (int argc, char ** argv)
     latt_size[3] = lat;
     double volume = latt_size[0]*latt_size[1]*latt_size[2]*latt_size[3];
     
-#ifdef AVX512
-    simd_layout[0] = 1;
-    simd_layout[1] = 2;
-    simd_layout[2] = 2;
-    simd_layout[3] = 2;
-#endif
-#if defined (AVX1)|| defined (AVX2)
-    simd_layout[0] = 1;
-    simd_layout[1] = 1;
-    simd_layout[2] = 2;
-    simd_layout[3] = 2;
-#endif
-#if defined (SSE4)
-    simd_layout[0] = 1;
-    simd_layout[1] = 1;
-    simd_layout[2] = 1;
-    simd_layout[3] = 2;
-#endif
-    
     GridCartesian Fine(latt_size,simd_layout,mpi_layout);
     GridRedBlackCartesian rbFine(latt_size,simd_layout,mpi_layout);
     GridParallelRNG       FineRNG(&Fine);

tests/Grid_nersc_io.cc

@@ -10,8 +10,12 @@ int main (int argc, char ** argv)
 {
   Grid_init(&argc,&argv);
 
-  std::vector<int> simd_layout({1,1,2,2});
+  std::vector<int> tmp_latt_size;
-  std::vector<int> mpi_layout ({1,1,1,1});
+  std::vector<int> simd_layout;
+  std::vector<int> mpi_layout;
+
+  GridParseLayout(argv,argc,mpi_layout,simd_layout,tmp_latt_size);
+
   std::vector<int> latt_size  ({16,16,16,32});
   std::vector<int> clatt_size  ({4,4,4,8});
   int orthodir=3;

tests/Grid_simd.cc

@@ -106,9 +106,11 @@ int main (int argc, char ** argv)
 {
   Grid_init(&argc,&argv);
 
-  std::vector<int> simd_layout({1,1,2,2});
+  std::vector<int> latt_size;
-  std::vector<int> mpi_layout ({1,1,1,1});
+  std::vector<int> simd_layout;
-  std::vector<int> latt_size  ({8,8,8,8});
+  std::vector<int> mpi_layout;
+
+  GridParseLayout(argv,argc,mpi_layout,simd_layout,latt_size);
     
   GridCartesian     Grid(latt_size,simd_layout,mpi_layout);
   std::vector<int> seeds({1,2,3,4});

tests/Grid_stencil.cc

@@ -8,9 +8,11 @@ int main (int argc, char ** argv)
 {
   Grid_init(&argc,&argv);
 
-  std::vector<int> simd_layout({1,1,2,2});
+  std::vector<int> latt_size;
-  std::vector<int> mpi_layout ({2,2,1,2});
+  std::vector<int> simd_layout;
-  std::vector<int> latt_size  ({8,8,8,8});
+  std::vector<int> mpi_layout;
+
+  GridParseLayout(argv,argc,mpi_layout,simd_layout,latt_size);
 
   double volume = latt_size[0]*latt_size[1]*latt_size[2]*latt_size[3];