Possible trash commit

extras/Hadrons/AllToAllVectors.hpp

@@ -102,7 +102,6 @@ class A2AModesSchurDiagTwo
 
     void low_mode_v(Matrix &action, const Field &evec, const RealD &eval, Field &vout_5d, Field &vout_4d)
     {
-
         GridBase *grid = action.RedBlackGrid();
         Field src_o(grid);
         Field sol_e(grid);

extras/Hadrons/Modules/MContraction/MesonFieldGamma.hpp

@@ -39,7 +39,8 @@ class TMesonFieldGamma : public Module<MesonFieldPar>
       public:
         GRID_SERIALIZABLE_CLASS_MEMBERS(Result,
                                         Gamma::Algebra, gamma,
-                                        std::vector<std::vector<std::vector<ComplexD>>>, MesonField);
+                                        std::vector<std::vector<std::vector<ComplexD>>>, MesonField,
+                                        ComplexD, last);
     };
 
   public:
@@ -178,7 +179,7 @@ void TMesonFieldGamma<FImpl>::execute(void)
             LOG(Message) << "MF for i = " << i << " of " << N << std::endl;
         }
     }
-
+    result[0].last = MesonField_ij[7];
     saveResult(par().output, "meson", result);
 }
 

lib/lattice/Lattice_reduction.h

@@ -39,8 +39,9 @@ template<class vobj> inline RealD norm2(const Lattice<vobj> &arg){
 
 // Double inner product
 template<class vobj>
-inline ComplexD innerProduct(const Lattice<vobj> &left,const Lattice<vobj> &right) 
+inline ComplexD innerProduct(const Lattice<vobj> &left,const Lattice<vobj> &right)
 {
+  std::cout << GridLogMessage << "Start alloc innerProduct" << std::endl;
   typedef typename vobj::scalar_type scalar_type;
   typedef typename vobj::vector_typeD vector_type;
   GridBase *grid = left._grid;
@@ -49,6 +50,8 @@ inline ComplexD innerProduct(const Lattice<vobj> &left,const Lattice<vobj> &righ
   ComplexD  inner;
   Vector<ComplexD> sumarray(grid->SumArraySize()*pad);
 
+  std::cout << GridLogMessage << "End alloc innerProduct" << std::endl;
+  std::cout << GridLogMessage << "Start parallel for innerProduct" << std::endl;
   parallel_for(int thr=0;thr<grid->SumArraySize();thr++){
     int nwork, mywork, myoff;
     GridThread::GetWork(left._grid->oSites(),thr,mywork,myoff);
@@ -62,13 +65,18 @@ inline ComplexD innerProduct(const Lattice<vobj> &left,const Lattice<vobj> &righ
     ComplexD tmp = Reduce(TensorRemove(vinner)) ;
     vstream(sumarray[thr*pad],tmp);
   }
-  
+  std::cout << GridLogMessage << "End parallel for innerProduct" << std::endl;
+
+  std::cout << GridLogMessage << "Start inner sum innerProduct" << std::endl;
   inner=0.0;
   for(int i=0;i<grid->SumArraySize();i++){
     inner = inner+sumarray[i*pad];
   } 
   right._grid->GlobalSum(inner);
   return inner;
+  std::cout << GridLogMessage << "End inner sum innerProduct" << std::endl;
+
+  std::cout << GridLogMessage << "End innerProduct" << std::endl;
 }
 
 /////////////////////////
@@ -274,6 +282,115 @@ template<class vobj> inline void sliceSum(const Lattice<vobj> &Data,std::vector<
   }
 }
 
+template<class vobj>
+static void mySliceInnerProductVector( std::vector<ComplexD> & result, const Lattice<vobj> &lhs,const Lattice<vobj> &rhs,int orthogdim)
+{
+  std::cout << GridLogMessage << "Start mySsliceInnerProductVector" << std::endl;
+
+  typedef typename vobj::scalar_type scalar_type;
+  std::vector<scalar_type> lsSum;
+  localSliceInnerProductVector(result, lhs, rhs, lsSum, orthogdim);
+  globalSliceInnerProductVector(result, lhs, lsSum, orthogdim);
+  std::cout << GridLogMessage << "End mySliceInnerProductVector" << std::endl;
+}
+
+template <class vobj>
+static void localSliceInnerProductVector(std::vector<ComplexD> &result, const Lattice<vobj> &lhs, const Lattice<vobj> &rhs, std::vector<typename vobj::scalar_type> &lsSum, int orthogdim)
+{
+  std::cout << GridLogMessage << "Start prep" << std::endl;
+  typedef typename vobj::vector_type   vector_type;
+  typedef typename vobj::scalar_type   scalar_type;
+  GridBase  *grid = lhs._grid;
+  assert(grid!=NULL);
+  conformable(grid,rhs._grid);
+
+  const int    Nd = grid->_ndimension;
+  const int Nsimd = grid->Nsimd();
+
+  assert(orthogdim >= 0);
+  assert(orthogdim < Nd);
+
+  int fd=grid->_fdimensions[orthogdim];
+  int ld=grid->_ldimensions[orthogdim];
+  int rd=grid->_rdimensions[orthogdim];
+  std::cout << GridLogMessage << "Start alloc" << std::endl;
+
+  std::vector<vector_type,alignedAllocator<vector_type> > lvSum(rd); // will locally sum vectors first
+  lsSum.resize(ld,scalar_type(0.0));                    // sum across these down to scalars
+  std::vector<iScalar<scalar_type>> extracted(Nsimd);   // splitting the SIMD  
+  std::cout << GridLogMessage << "End alloc" << std::endl;
+
+  result.resize(fd); // And then global sum to return the same vector to every node for IO to file
+  for(int r=0;r<rd;r++){
+    lvSum[r]=zero;
+  }
+
+  int e1=    grid->_slice_nblock[orthogdim];
+  int e2=    grid->_slice_block [orthogdim];
+  int stride=grid->_slice_stride[orthogdim];
+  std::cout << GridLogMessage << "End prep" << std::endl;
+  std::cout << GridLogMessage << "Start parallel inner product, _rd = " << rd << std::endl;
+  vector_type vv;
+  parallel_for(int r=0;r<rd;r++)
+  {
+
+    int so=r*grid->_ostride[orthogdim]; // base offset for start of plane 
+
+    for(int n=0;n<e1;n++){
+      for(int b=0;b<e2;b++){
+        int ss = so + n * stride + b;
+        vv = TensorRemove(innerProduct(lhs._odata[ss]._internal, rhs._odata[ss]._internal));
+        lvSum[r] = lvSum[r] + vv;
+      }
+    }
+  }
+  std::cout << GridLogMessage << "End parallel inner product" << std::endl;
+
+  // Sum across simd lanes in the plane, breaking out orthog dir.
+  std::vector<int> icoor(Nd);
+  for(int rt=0;rt<rd;rt++){
+
+    iScalar<vector_type> temp; 
+    temp._internal = lvSum[rt];
+    extract(temp,extracted);
+
+    for(int idx=0;idx<Nsimd;idx++){
+
+      grid->iCoorFromIindex(icoor,idx);
+
+      int ldx =rt+icoor[orthogdim]*rd;
+
+      lsSum[ldx]=lsSum[ldx]+extracted[idx]._internal;
+
+    }
+  }
+  std::cout << GridLogMessage << "End sum over simd lanes" << std::endl;
+}
+template <class vobj>
+static void globalSliceInnerProductVector(std::vector<ComplexD> &result, const Lattice<vobj> &lhs, std::vector<typename vobj::scalar_type> &lsSum, int orthogdim)
+{
+  typedef typename vobj::scalar_type scalar_type;
+  GridBase *grid = lhs._grid;
+  int fd = result.size();
+  int ld = lsSum.size();
+  // sum over nodes.
+  std::vector<scalar_type> gsum;
+  gsum.resize(fd, scalar_type(0.0));
+  std::cout << GridLogMessage << "Start of gsum[t] creation:" << std::endl;
+  for(int t=0;t<fd;t++){
+    int pt = t/ld; // processor plane
+    int lt = t%ld;
+    if ( pt == grid->_processor_coor[orthogdim] ) {
+      gsum[t]=lsSum[lt];
+    }
+  }
+  std::cout << GridLogMessage << "End of gsum[t] creation:" << std::endl;
+  std::cout << GridLogMessage << "Start of GlobalSumVector:" << std::endl;
+  grid->GlobalSumVector(&gsum[0], fd);
+  std::cout << GridLogMessage << "End of GlobalSumVector:" << std::endl;
+
+  result = gsum;
+}
 template<class vobj>
 static void sliceInnerProductVector( std::vector<ComplexD> & result, const Lattice<vobj> &lhs,const Lattice<vobj> &rhs,int orthogdim) 
 {

lib/tensors/Tensor_inner.h

@@ -106,6 +106,7 @@ inline vRealD innerProductD(const vRealF &l,const vRealF &r){
     typedef decltype(innerProduct(lhs._internal[0],rhs._internal[0])) ret_t;
     iScalar<ret_t> ret;
     ret=zero;
+    // std::cout << GridLogMessage << "innerProduct iVector" << std::endl;
     for(int c1=0;c1<N;c1++){
       ret._internal += innerProduct(lhs._internal[c1],rhs._internal[c1]);
     }
@@ -129,9 +130,34 @@ inline vRealD innerProductD(const vRealF &l,const vRealF &r){
   {
     typedef decltype(innerProduct(lhs._internal,rhs._internal)) ret_t;
     iScalar<ret_t> ret;
+    // std::cout << GridLogMessage << "innerProduct iScalar" << std::endl;
+
     ret._internal = innerProduct(lhs._internal,rhs._internal);
     return ret;
   }
+  template<class l,class r,int N> inline
+  auto myInnerProduct (const iVector<l,N>& lhs,const iVector<r,N>& rhs) -> iScalar<decltype(innerProduct(lhs._internal[0],rhs._internal[0]))>
+  {
+    typedef decltype(innerProduct(lhs._internal[0],rhs._internal[0])) ret_t;
+    iScalar<ret_t> ret;
+    ret=zero;
+    std::cout << GridLogMessage << "myInnerProduct iVector, N = " << N << std::endl;
+    for(int c1=0;c1<N;c1++){
+      ret._internal += innerProduct(lhs._internal[c1],rhs._internal[c1]);
+    }
+    return ret;
+  }
+
+  template<class l,class r> inline
+  auto myInnerProduct (const iScalar<l>& lhs,const iScalar<r>& rhs) -> iScalar<decltype(innerProduct(lhs._internal,rhs._internal))>
+  {
+    typedef decltype(innerProduct(lhs._internal,rhs._internal)) ret_t;
+    iScalar<ret_t> ret;
+    std::cout << GridLogMessage << "myInnerProduct iScalar" << std::endl;
+
+    ret._internal = myInnerProduct(lhs._internal,rhs._internal);
+    return ret;
+  }
 
 }
 #endif