Merge branch 'feature/distil' of https://github.com/mmphys/Grid into feature/distil

Grid/qcd/utils/A2Autils.h

@@ -101,6 +101,186 @@ public:
 #endif
 };
 
+/*
+template <class FImpl>
+template <typename TensorType>
+void A2Autils<FImpl>::BaryonField(TensorType &mat, 
+				 const FermionField *one,
+				 const FermionField *two,
+				 const FermionField *three,
+				 std::vector<Gamma::Algebra> gammaA,
+				 std::vector<Gamma::Algebra> gammaB,
+				 const std::vector<ComplexField > &mom,
+				 int orthogdim, double *t_kernel, double *t_gsum) 
+{
+  typedef typename FImpl::SiteSpinor vobj;
+
+  typedef typename vobj::scalar_object sobj;
+  typedef typename vobj::scalar_type scalar_type;
+  typedef typename vobj::vector_type vector_type;
+
+  typedef iSpinMatrix<vector_type> SpinMatrix_v;
+  typedef iSpinMatrix<scalar_type> SpinMatrix_s;
+  
+  int oneBlock = mat.dimension(3); 
+  int twoBlock = mat.dimension(4);
+  int threeBlock = mat.dimension(5);
+
+  GridBase *grid = one[0]._grid;
+  
+  const int    Nd = grid->_ndimension;
+  const int Nsimd = grid->Nsimd();
+
+  int Nt     = grid->GlobalDimensions()[orthogdim];
+  int Ngamma = gammaA.size();
+  assert (Ngamma = gammaB.size()); // Only combinatin of two gammas gives correct operator!
+  int Nmom   = mom.size();
+
+  int fd=grid->_fdimensions[orthogdim];
+  int ld=grid->_ldimensions[orthogdim];
+  int rd=grid->_rdimensions[orthogdim];
+
+  // will locally sum vectors first
+  // sum across these down to scalars
+  // splitting the SIMD
+  int MFrvol = rd*oneBlock*twoBlock*threeBlock*Nmom;
+  int MFlvol = ld*oneBlock*twoBlock*threeBlock*Nmom;
+
+  Vector<SpinMatrix_v > lvSum(MFrvol);
+  parallel_for (int r = 0; r < MFrvol; r++){
+    lvSum[r] = zero;
+  }
+
+  Vector<SpinMatrix_s > lsSum(MFlvol);             
+  parallel_for (int r = 0; r < MFlvol; r++){
+    lsSum[r]=scalar_type(0.0);
+  }
+
+  int e1=    grid->_slice_nblock[orthogdim];
+  int e2=    grid->_slice_block [orthogdim];
+  int stride=grid->_slice_stride[orthogdim];
+
+  // potentially wasting cores here if local time extent too small
+  if (t_kernel) *t_kernel = -usecond();
+  parallel_for(int r=0;r<rd;r++){
+
+    int so=r*grid->_ostride[orthogdim]; // base offset for start of plane 
+
+    // first, the diquark two*gammaB*three
+    for(int n=0;n<e1;n++){
+      for(int b=0;b<e2;b++){
+
+	int ss= so+n*stride+b;
+
+	  for(int j=0;j<twoBlock;j++){
+
+	    auto two_j = two[j]._odata[ss];
+
+	  for(int k=0;k<threeBlock;k++){
+
+	    auto three_k = three[j]._odata[ss];
+
+	    SpinMatrix_v vv;
+
+	    for(int s1=0;s1<Ns;s1++){
+	    for(int s2=0;s2<Ns;s2++){
+	      vv()(s1,s2)() = two_j()(s2)(0) * three_k()(s1)(0)   //make this a colorMatrix for the diquark???
+		+             two_j()(s2)(1) * three_k()(s1)(1)
+		+             two_j()(s2)(2) * three_k()(s1)(2);
+	    }}
+	    
+	    // After getting the sitewise product do the mom phase loop
+	    int base = Nmom*i+Nmom*Lblock*j+Nmom*Lblock*Rblock*r;
+	    for ( int m=0;m<Nmom;m++){
+	      int idx = m+base;
+	      auto phase = mom[m]._odata[ss];
+	      mac(&lvSum[idx],&vv,&phase);
+	    }
+	 
+	  }
+	}
+      }
+    }
+  }
+
+
+  // Sum across simd lanes in the plane, breaking out orthog dir.
+  parallel_for(int rt=0;rt<rd;rt++){
+
+    std::vector<int> icoor(Nd);
+    std::vector<SpinMatrix_s> extracted(Nsimd);               
+
+    for(int i=0;i<Lblock;i++){
+    for(int j=0;j<Rblock;j++){
+    for(int m=0;m<Nmom;m++){
+
+      int ij_rdx = m+Nmom*i+Nmom*Lblock*j+Nmom*Lblock*Rblock*rt;
+
+      extract(lvSum[ij_rdx],extracted);
+
+      for(int idx=0;idx<Nsimd;idx++){
+
+	grid->iCoorFromIindex(icoor,idx);
+
+	int ldx    = rt+icoor[orthogdim]*rd;
+
+	int ij_ldx = m+Nmom*i+Nmom*Lblock*j+Nmom*Lblock*Rblock*ldx;
+
+	lsSum[ij_ldx]=lsSum[ij_ldx]+extracted[idx];
+
+      }
+    }}}
+  }
+  if (t_kernel) *t_kernel += usecond();
+  assert(mat.dimension(0) == Nmom);
+  assert(mat.dimension(1) == Ngamma);
+  assert(mat.dimension(2) == Nt);
+
+  // ld loop and local only??
+  int pd = grid->_processors[orthogdim];
+  int pc = grid->_processor_coor[orthogdim];
+  parallel_for_nest2(int lt=0;lt<ld;lt++)
+  {
+    for(int pt=0;pt<pd;pt++){
+      int t = lt + pt*ld;
+      if (pt == pc){
+	for(int i=0;i<Lblock;i++){
+	  for(int j=0;j<Rblock;j++){
+	    for(int m=0;m<Nmom;m++){
+	      int ij_dx = m+Nmom*i + Nmom*Lblock * j + Nmom*Lblock * Rblock * lt;
+	      for(int mu=0;mu<Ngamma;mu++){
+		// this is a bit slow
+		mat(m,mu,t,i,j) = trace(lsSum[ij_dx]*Gamma(gammas[mu]));
+	      }
+	    }
+	  }
+	}
+      } else { 
+	const scalar_type zz(0.0);
+	for(int i=0;i<Lblock;i++){
+	  for(int j=0;j<Rblock;j++){
+	    for(int mu=0;mu<Ngamma;mu++){
+	      for(int m=0;m<Nmom;m++){
+		mat(m,mu,t,i,j) =zz;
+	      }
+	    }
+	  }
+	}
+      }
+    }
+  }
+
+  ////////////////////////////////////////////////////////////////////
+  // This global sum is taking as much as 50% of time on 16 nodes
+  // Vector size is 7 x 16 x 32 x 16 x 16 x sizeof(complex) = 2MB - 60MB depending on volume
+  // Healthy size that should suffice
+  ////////////////////////////////////////////////////////////////////
+  if (t_gsum) *t_gsum = -usecond();
+  grid->GlobalSumVector(&mat(0,0,0,0,0),Nmom*Ngamma*Nt*Lblock*Rblock);
+  if (t_gsum) *t_gsum += usecond();
+}
+*/
+
 template <class FImpl>
 template <typename TensorType>
 void A2Autils<FImpl>::MesonField(TensorType &mat, 

Grid/serialisation/BaseIO.h

@@ -172,7 +172,7 @@ namespace Grid {
   template <typename ETensor, typename Lambda>
   typename std::enable_if<EigenIO::is_tensor_of_scalar<ETensor>::value, void>::type
   for_all_do_lambda( Lambda lambda, typename ETensor::Scalar &scalar, typename ETensor::Index &Seq,
-                    std::array<std::size_t, ETensor::NumIndices + EigenIO::Traits<typename ETensor::Scalar>::rank_non_trivial> &MyIndex)
+                    std::array<std::size_t, ETensor::NumIndices + EigenIO::Traits<typename ETensor::Scalar>::rank> &MyIndex)
   {
     lambda( scalar, Seq++, MyIndex );
   }
@@ -180,8 +180,8 @@ namespace Grid {
   // for_all helper function to call the lambda
   template <typename ETensor, typename Lambda>
   typename std::enable_if<EigenIO::is_tensor_of_container<ETensor>::value, void>::type
-  for_all_do_lambda( Lambda lambda, typename ETensor::Scalar &scalar, typename ETensor::Index Seq,
+  for_all_do_lambda( Lambda lambda, typename ETensor::Scalar &scalar, typename ETensor::Index &Seq,
-                    std::array<std::size_t, ETensor::NumIndices + EigenIO::Traits<typename ETensor::Scalar>::rank_non_trivial> &MyIndex)
+                    std::array<std::size_t, ETensor::NumIndices + EigenIO::Traits<typename ETensor::Scalar>::rank> &MyIndex)
   {
     using Scalar = typename ETensor::Scalar; // This could be a Container - we'll check later
     const auto InnerRank = EigenIO::Traits<Scalar>::rank_non_trivial;
@@ -190,7 +190,7 @@ namespace Grid {
       lambda(Source, Seq++, MyIndex );
       // Now increment SubIndex
       for( auto i = InnerRank - 1; i != -1 && ++MyIndex[rank + i] == EigenIO::Traits<Scalar>::DimensionNT(i); i-- )
-        MyIndex[i] = 0;
+        MyIndex[rank + i] = 0;
     }
   }
   
@@ -206,7 +206,7 @@ namespace Grid {
     assert( NumScalars > 0 );
     using Index = typename ETensor::Index;
     Index ScalarElementCount{1};
-    const auto InnerRank = EigenIO::Traits<Scalar>::rank_non_trivial;
+    const auto InnerRank = EigenIO::Traits<Scalar>::rank;
     const auto rank{ETensor::NumIndices};
     std::array<std::size_t, rank + InnerRank> Dims;
     for(auto i = 0; i < rank; i++ ) {
@@ -218,11 +218,11 @@ namespace Grid {
     }
     // Check that the number of containers is correct ... and we didn't lose anything in conversions
     assert( NumScalars == ScalarElementCount );
-    // If the Scalar is actually a container, add the inner Scalar's non-trivial dimensions
+    // If the Scalar is actually a container, add the inner Scalar's dimensions
     size_t InnerScalarCount{1};
     for(auto i = 0; i < InnerRank; i++ ) {
-      auto dim = EigenIO::Traits<Scalar>::DimensionNT(i);
+      auto dim = EigenIO::Traits<Scalar>::Dimension(i);
-      assert( dim > 1 );
+      assert( dim > 0 );
       Dims[rank + i] = static_cast<std::size_t>(dim);
       assert( Dims[rank + i] == dim ); // check we didn't lose anything in the conversion
       InnerScalarCount *= dim;
@@ -242,11 +242,12 @@ namespace Grid {
       } else {
         for( auto i = 0; i < rank && ++MyIndex[i] == Dims[i]; i++ )
           MyIndex[i] = 0;
-        Seq = 0;
+        size_t NewSeq = 0;
         for( auto i = 0; i < rank + InnerRank ; i++ ) {
-          Seq *= Dims[i];
+          NewSeq *= Dims[i];
-          Seq += MyIndex[i];
+          NewSeq += MyIndex[i];
         }
+        Seq = static_cast<Index>( NewSeq );
       }
       pScalar++;
     }
@@ -271,7 +272,7 @@ namespace Grid {
   {
     using Traits = EigenIO::Traits<typename ETensor::Scalar>;
     using scalar_type = typename Traits::scalar_type;
-    for_all( ET, [&](scalar_type &c, typename ETensor::Index n, const std::array<size_t, ETensor::NumIndices + Traits::rank_non_trivial> &Dims ) {
+    for_all( ET, [&](scalar_type &c, typename ETensor::Index n, const std::array<size_t, ETensor::NumIndices + Traits::rank> &Dims ) {
       c = Inc * static_cast<typename RealType<scalar_type>::type>(n);
     } );
   }
@@ -291,7 +292,7 @@ namespace Grid {
       std::cout << pName;
     for( auto i = 0 ; i < rank; i++ ) std::cout << "[" << dims[i] << "]";
     std::cout << " in memory order:" << std::endl;
-    for_all( t, [&](typename Traits::scalar_type &c, typename T::Index index, const std::array<size_t, T::NumIndices + Traits::rank_non_trivial> &Dims ){
+    for_all( t, [&](typename Traits::scalar_type &c, typename T::Index index, const std::array<size_t, T::NumIndices + Traits::rank> &Dims ){
       std::cout << "  ";
       for( auto dim : Dims )
         std::cout << "[" << dim << "]";

Grid/serialisation/TextIO.h

@@ -51,6 +51,8 @@ namespace Grid
     void writeDefault(const std::string &s, const U &x);
     template <typename U>
     void writeDefault(const std::string &s, const std::vector<U> &x);
+    template <typename U>
+    void writeMultiDim(const std::string &s, const std::vector<size_t> & Dimensions, const U * pDataRowMajor, size_t NumElements);
   private:
     void indent(void);
   private:
@@ -69,6 +71,8 @@ namespace Grid
     void readDefault(const std::string &s, U &output);
     template <typename U>
     void readDefault(const std::string &s, std::vector<U> &output);
+    template <typename U>
+    void readMultiDim(const std::string &s, std::vector<U> &buf, std::vector<size_t> &dim);
   private:
     void checkIndent(void);
   private:
@@ -95,7 +99,18 @@ namespace Grid
       write(s, x[i]);
     }
   }
-  
+
+  template <typename U>
+  void TextWriter::writeMultiDim(const std::string &s, const std::vector<size_t> & Dimensions, const U * pDataRowMajor, size_t NumElements)
+  {
+    uint64_t Rank = Dimensions.size();
+    write(s, Rank);
+    for( uint64_t d : Dimensions )
+      write(s, d);
+    while( NumElements-- )
+      write(s, *pDataRowMajor++);
+  }
+
   // Reader template implementation ////////////////////////////////////////////
   template <typename U>
   void TextReader::readDefault(const std::string &s, U &output)
@@ -121,6 +136,23 @@ namespace Grid
       read("", output[i]);
     }
   }
+
+  template <typename U>
+  void TextReader::readMultiDim(const std::string &s, std::vector<U> &buf, std::vector<size_t> &dim)
+  {
+    const char sz[] = "";
+    uint64_t Rank;
+    read(sz, Rank);
+    dim.resize( Rank );
+    size_t NumElements = 1;
+    for( auto &d : dim ) {
+      read(sz, d);
+      NumElements *= d;
+    }
+    buf.resize( NumElements );
+    for( auto &x : buf )
+      read(s, x);
+  }
 }
 
 #endif

Grid/serialisation/XmlIO.h

@@ -57,6 +57,8 @@ namespace Grid
     void writeDefault(const std::string &s, const U &x);
     template <typename U>
     void writeDefault(const std::string &s, const std::vector<U> &x);
+    template <typename U>
+    void writeMultiDim(const std::string &s, const std::vector<size_t> & Dimensions, const U * pDataRowMajor, size_t NumElements);
     std::string docString(void);
     std::string string(void);
   private:
@@ -79,6 +81,8 @@ namespace Grid
     void readDefault(const std::string &s, U &output);
     template <typename U>
     void readDefault(const std::string &s, std::vector<U> &output);
+    template <typename U>
+    void readMultiDim(const std::string &s, std::vector<U> &buf, std::vector<size_t> &dim);
     void readCurrentSubtree(std::string &s);
   private:
     void checkParse(const pugi::xml_parse_result &result, const std::string name);
@@ -122,13 +126,30 @@ namespace Grid
   void XmlWriter::writeDefault(const std::string &s, const std::vector<U> &x)
   {
     push(s);
-    for (auto &x_i: x)
+    for( auto &u : x )
     {
-      write("elem", x_i);
+      write("elem", u);
     }
     pop();
   }
-  
+
+  template <typename U>
+  void XmlWriter::writeMultiDim(const std::string &s, const std::vector<size_t> & Dimensions, const U * pDataRowMajor, size_t NumElements)
+  {
+    push(s);
+    size_t count = 1;
+    const size_t Rank = Dimensions.size();
+    write("rank", Rank );
+    for( auto d : Dimensions ) {
+      write("dim", d);
+      count *= d;
+    }
+    assert( count == NumElements && "XmlIO : element count doesn't match dimensions" );
+    while (NumElements--)
+      write("elem", *pDataRowMajor++);
+    pop();
+  }
+
   // Reader template implementation ////////////////////////////////////////////
   template <typename U>
   void XmlReader::readDefault(const std::string &s, U &output)
@@ -145,25 +166,47 @@ namespace Grid
   template <typename U>
   void XmlReader::readDefault(const std::string &s, std::vector<U> &output)
   {
-    std::string    buf;
-    unsigned int   i = 0;
-    
     if (!push(s))
     {
       std::cout << GridLogWarning << "XML: cannot open node '" << s << "'";
       std::cout << std::endl;
-
+    } else {
-      return; 
+      for(unsigned int i = 0; node_.child("elem"); )
+      {
+        output.resize(i + 1);
+        read("elem", output[i++]);
+        node_.child("elem").set_name("elem-done");
+      }
+      pop();
+    }
+  }
+
+  template <typename U>
+  void XmlReader::readMultiDim(const std::string &s, std::vector<U> &buf, std::vector<size_t> &dim)
+  {
+    if (!push(s))
+    {
+      std::cout << GridLogWarning << "XML: cannot open node '" << s << "'";
+      std::cout << std::endl;
+    } else {
+      size_t Rank;
+      read("rank", Rank);
+      dim.resize( Rank );
+      size_t NumElements = 1;
+      for( auto &d : dim )
+      {
+        read("dim", d);
+        node_.child("dim").set_name("dim-done");
+        NumElements *= d;
+      }
+      buf.resize( NumElements );
+      for( auto &x : buf )
+      {
+        read("elem", x);
+        node_.child("elem").set_name("elem-done");
+      }
+      pop();
     }
-    while (node_.child("elem"))
-    {
-      output.resize(i + 1);
-      read("elem", output[i]);
-      node_.child("elem").set_name("elem-done");
-      i++;
-    }
-    pop();
   }
-  
 }
 #endif

tests/IO/Test_serialisation.cc

@@ -86,7 +86,7 @@ void ioTest(const std::string &filename, const O &object, const std::string &nam
   // writer needs to be destroyed so that writing physically happens
   {
     W writer(filename);
-
+    writer.setPrecision(std::numeric_limits<double>::digits10 + 1);
     write(writer, tag , object);
   }
 
@@ -113,7 +113,7 @@ typedef Eigen::Tensor<TestScalar, 3, Eigen::StorageOptions::RowMajor> TensorRank
 typedef Eigen::TensorFixedSize<TestScalar, Eigen::Sizes<9,4,2>, Eigen::StorageOptions::RowMajor> Tensor_9_4_2;
 typedef std::vector<Tensor_9_4_2> aTensor_9_4_2;
 typedef Eigen::TensorFixedSize<SpinColourVector, Eigen::Sizes<6,5>> LSCTensor;
-#ifdef DEBUG
+#ifndef DEBUG
 typedef Eigen::TensorFixedSize<iMatrix<iVector<iMatrix<iVector<LorentzColourMatrix,5>,2>,7>,3>, Eigen::Sizes<2,4,11,10,9>, Eigen::StorageOptions::RowMajor> LCMTensor;
 #endif
 
@@ -141,14 +141,14 @@ public:
   , tensorRank3(7,3,2)
   , atensor_9_4_2(3)
   {
-    Grid_complex<double> Flag{1,-3.1415927};
+    //Grid_complex<double> Flag{1,-3.1415927}; // Gives errors on readback for text types
+    Grid_complex<double> Flag{1,-1};
     SequentialInit(Perambulator,  Flag);
     SequentialInit(Perambulator2, Flag);
     SequentialInit(tensorRank5UShort);
     SequentialInit(tensorRank3, Flag);
     SequentialInit(tensor_9_4_2, Flag);
-    for( auto &t : atensor_9_4_2 )
+    for( auto &t : atensor_9_4_2 ) SequentialInit(t, Flag);
-      SequentialInit(t, Flag);
     SequentialInit( MyLSCTensor );
   }
 };
@@ -193,7 +193,7 @@ void EigenHdf5IOTest(const char * pszExtension)
   }
   TensorWriteRead ( LSCTensor )
   TensorWriteReadLarge( PerambIOTestClass )
-#ifdef DEBUG
+#ifndef DEBUG
   std::cout << "sizeof( LCMTensor ) = " << sizeof( LCMTensor ) / 1024 / 1024 << " MB" << std::endl;
   TensorWriteReadLarge ( LCMTensor )
   // Also write > 4GB of complex numbers (I suspect this will fail inside Hdf5)
@@ -297,6 +297,10 @@ int main(int argc,char **argv)
 #endif
   std::cout << "\n==== detailed binary tensor tests (Grid::EigenIO)" << std::endl;
   EigenHdf5IOTest<BinaryWriter, BinaryReader>(".bin");
+  std::cout << "\n==== detailed text tensor tests (Grid::EigenIO)" << std::endl;
+  EigenHdf5IOTest<TextWriter, TextReader>(".dat");
+  std::cout << "\n==== detailed xml tensor tests (Grid::EigenIO)" << std::endl;
+  EigenHdf5IOTest<XmlWriter, XmlReader>(".xml");
 
   std::cout << "\n==== vector flattening/reconstruction" << std::endl;
   typedef std::vector<std::vector<std::vector<double>>> vec3d;