From 3b2d805398c0691dcb61e9aa4cbbadb7881aa4b6 Mon Sep 17 00:00:00 2001
From: Daniel Richtmann <daniel.richtmann@gmail.com>
Date: Sat, 20 Jan 2018 18:48:53 +0100
Subject: [PATCH] WilsonMG: Some first steps towards coarse spin dofs; not
 compiling yet

A failing conversion from the innermost type (Grid::Simd<...>) to a coarse
scalar (triple iScalar) in blockPromote prohibits this commit from working.
---
 lib/algorithms/CoarsenedMatrix.h | 35 +++++++++++++++++++-------------
 lib/lattice/Lattice_transfer.h   | 26 ++++++++++++++----------
 2 files changed, 36 insertions(+), 25 deletions(-)

diff --git a/lib/algorithms/CoarsenedMatrix.h b/lib/algorithms/CoarsenedMatrix.h
index 8af8d7ac..04ac18ae 100644
--- a/lib/algorithms/CoarsenedMatrix.h
+++ b/lib/algorithms/CoarsenedMatrix.h
@@ -93,12 +93,16 @@ namespace Grid {
   template<class Fobj,class CComplex,int nbasis>
   class Aggregation   {
   public:
-    typedef iVector<CComplex,nbasis >             siteVector;
-    typedef Lattice<siteVector>                 CoarseVector;
-    typedef Lattice<iMatrix<CComplex,nbasis > > CoarseMatrix;
 
-    typedef Lattice< CComplex >   CoarseScalar; // used for inner products on fine field
-    typedef Lattice<Fobj >        FineField;
+    typedef typename GridTypeMapper<CComplex>::vector_type     innerType;
+    typedef iScalar<iScalar<iScalar<innerType         >    > > siteScalar; // used for inner products on fine field
+    typedef iScalar<iVector<iVector<innerType, nbasis >, 1 > > siteVector;
+    typedef iScalar<iMatrix<iMatrix<innerType, nbasis >, 1 > > siteMatrix;
+    typedef Lattice<siteScalar> CoarseScalar; // used for inner products on fine field
+    typedef Lattice<siteVector> CoarseVector;
+    typedef Lattice<siteMatrix> CoarseMatrix;
+
+    typedef Lattice<Fobj>       FineField;
 
     GridBase *CoarseGrid;
     GridBase *FineGrid;
@@ -129,7 +133,7 @@ namespace Grid {
 	blockProject(iProj,subspace[i],subspace);
 	eProj=zero; 
 	parallel_for(int ss=0;ss<CoarseGrid->oSites();ss++){
-	  eProj._odata[ss](i)=CComplex(1.0);
+	  eProj._odata[ss]()(0)(i)=innerType(1.0);
 	}
 	eProj=eProj - iProj;
 	std::cout<<GridLogMessage<<"Orthog check error "<<i<<" " << norm2(eProj)<<std::endl;
@@ -239,15 +243,18 @@ namespace Grid {
   // Fine Object == (per site) type of fine field
   // nbasis      == number of deflation vectors
   template<class Fobj,class CComplex,int nbasis>
-  class CoarsenedMatrix : public SparseMatrixBase<Lattice<iVector<CComplex,nbasis > > >  {
+  class CoarsenedMatrix : public SparseMatrixBase<Lattice<iScalar<iVector<iVector<typename GridTypeMapper<CComplex>::vector_type, nbasis >, 1 > > > >  {
   public:
     
-    typedef iVector<CComplex,nbasis >             siteVector;
-    typedef Lattice<siteVector>                 CoarseVector;
-    typedef Lattice<iMatrix<CComplex,nbasis > > CoarseMatrix;
+    typedef typename GridTypeMapper<CComplex>::vector_type     innerType;
+    typedef iScalar<iScalar<iScalar<innerType         >    > > siteScalar;
+    typedef iScalar<iVector<iVector<innerType, nbasis >, 1 > > siteVector;
+    typedef iScalar<iMatrix<iMatrix<innerType, nbasis >, 1 > > siteMatrix;
+    typedef Lattice<siteScalar> CoarseScalar; // used for inner products on fine field
+    typedef Lattice<siteVector> CoarseVector;
+    typedef Lattice<siteMatrix> CoarseMatrix;
 
-    typedef Lattice< CComplex >   CoarseScalar; // used for inner products on fine field
-    typedef Lattice<Fobj >        FineField;
+    typedef Lattice<Fobj>       FineField;
 
     ////////////////////
     // Data members
@@ -387,9 +394,9 @@ namespace Grid {
 	  parallel_for(int ss=0;ss<Grid()->oSites();ss++){
 	    for(int j=0;j<nbasis;j++){
 	      if( disp!= 0 ) {
-		A[p]._odata[ss](j,i) = oProj._odata[ss](j);
+		A[p]._odata[ss]()(0,0)(j,i) = oProj._odata[ss]()(0)(j);
 	      }
-	      A[self_stencil]._odata[ss](j,i) =	A[self_stencil]._odata[ss](j,i) + iProj._odata[ss](j);
+	      A[self_stencil]._odata[ss]()(0,0)(j,i) =	A[self_stencil]._odata[ss]()(0,0)(j,i) + iProj._odata[ss]()(0)(j);
 	    }
 	  }
 	}
diff --git a/lib/lattice/Lattice_transfer.h b/lib/lattice/Lattice_transfer.h
index 32c15d22..5a04d551 100644
--- a/lib/lattice/Lattice_transfer.h
+++ b/lib/lattice/Lattice_transfer.h
@@ -80,8 +80,8 @@ inline void subdivides(GridBase *coarse,GridBase *fine)
   }
   
 
-template<class vobj,class CComplex,int nbasis>
-inline void blockProject(Lattice<iVector<CComplex,nbasis > > &coarseData,
+template<class vobj,class vobjC>
+inline void blockProject(Lattice<vobjC>                    &coarseData,
 			 const             Lattice<vobj>   &fineData,
 			 const std::vector<Lattice<vobj> > &Basis)
 {
@@ -90,7 +90,8 @@ inline void blockProject(Lattice<iVector<CComplex,nbasis > > &coarseData,
   int  _ndimension = coarse->_ndimension;
 
   // checks
-  assert( nbasis == Basis.size() );
+  assert((Basis.size() != 0) && ((Basis.size() & 0x1) == 0));
+  auto nbasis = Basis.size();
   subdivides(coarse,fine); 
   for(int i=0;i<nbasis;i++){
     conformable(Basis[i],fineData);
@@ -118,8 +119,8 @@ inline void blockProject(Lattice<iVector<CComplex,nbasis > > &coarseData,
 PARALLEL_CRITICAL
     for(int i=0;i<nbasis;i++) {
 
-      coarseData._odata[sc](i)=coarseData._odata[sc](i)
-	+ innerProduct(Basis[i]._odata[sf],fineData._odata[sf]);
+      coarseData._odata[sc]()(0)(i)=coarseData._odata[sc]()(0)(i)
+        + TensorRemove(innerProduct(Basis[i]._odata[sf],fineData._odata[sf]));
 
     }
   }
@@ -285,9 +286,9 @@ inline void blockOrthogonalise(Lattice<CComplex> &ip,std::vector<Lattice<vobj> >
   }
 }
 
-template<class vobj,class CComplex,int nbasis>
-inline void blockPromote(const Lattice<iVector<CComplex,nbasis > > &coarseData,
-			 Lattice<vobj>   &fineData,
+template<class vobj,class vobjC>
+inline void blockPromote(const Lattice<vobjC>              &coarseData,
+			 Lattice<vobj>                     &fineData,
 			 const std::vector<Lattice<vobj> > &Basis)
 {
   GridBase * fine  = fineData._grid;
@@ -295,7 +296,9 @@ inline void blockPromote(const Lattice<iVector<CComplex,nbasis > > &coarseData,
   int  _ndimension = coarse->_ndimension;
 
   // checks
-  assert( nbasis == Basis.size() );
+  assert((Basis.size() != 0) && ((Basis.size() & 0x1) == 0));
+  auto nbasis = Basis.size();
+
   subdivides(coarse,fine); 
   for(int i=0;i<nbasis;i++){
     conformable(Basis[i]._grid,fine);
@@ -319,9 +322,10 @@ inline void blockPromote(const Lattice<iVector<CComplex,nbasis > > &coarseData,
       for(int d=0;d<_ndimension;d++) coor_c[d]=coor_f[d]/block_r[d];
       Lexicographic::IndexFromCoor(coor_c,sc,coarse->_rdimensions);
       
+      // TODO: These lines here prevent this commit from working
       for(int i=0;i<nbasis;i++) {
-	if(i==0) fineData._odata[sf]=coarseData._odata[sc](i) * Basis[i]._odata[sf];
-	else     fineData._odata[sf]=fineData._odata[sf]+coarseData._odata[sc](i)*Basis[i]._odata[sf];
+	if(i==0) fineData._odata[sf]=coarseData._odata[sc]()(0)(i) * Basis[i]._odata[sf];
+	else     fineData._odata[sf]=fineData._odata[sf]+coarseData._odata[sc]()(0)(i)*Basis[i]._odata[sf];
       }
     }
   }