From df8c208f5c1aa9283a98e04935f6ca6f765baa80 Mon Sep 17 00:00:00 2001 From: Daniel Richtmann Date: Wed, 16 May 2018 16:02:17 +0200 Subject: [PATCH] WilsonMG: Revert CoarsenedMatrix.h and Lattice_transfer.h back to state of develop branch --- lib/algorithms/CoarsenedMatrix.h | 102 +++++++++---------------------- lib/lattice/Lattice_transfer.h | 29 ++++----- 2 files changed, 40 insertions(+), 91 deletions(-) diff --git a/lib/algorithms/CoarsenedMatrix.h b/lib/algorithms/CoarsenedMatrix.h index 15db315a..8af8d7ac 100644 --- a/lib/algorithms/CoarsenedMatrix.h +++ b/lib/algorithms/CoarsenedMatrix.h @@ -93,16 +93,12 @@ namespace Grid { template class Aggregation { public: + typedef iVector siteVector; + typedef Lattice CoarseVector; + typedef Lattice > CoarseMatrix; - typedef typename CComplex::vector_type innerType; - typedef iScalar > > siteScalar; // used for inner products on fine field - typedef iScalar, 1 > > siteVector; - typedef iScalar, 1 > > siteMatrix; - typedef Lattice CoarseScalar; // used for inner products on fine field - typedef Lattice CoarseVector; - typedef Lattice CoarseMatrix; - - typedef Lattice FineField; + typedef Lattice< CComplex > CoarseScalar; // used for inner products on fine field + typedef Lattice FineField; GridBase *CoarseGrid; GridBase *FineGrid; @@ -119,12 +115,12 @@ namespace Grid { void Orthogonalise(void){ CoarseScalar InnerProd(CoarseGrid); - std::cout << GridLogMessage <<"Gram-Schmidt pass 1"<oSites();ss++){ - eProj._odata[ss]()(0)(i)=innerType(1.0); + eProj._odata[ss](i)=CComplex(1.0); } eProj=eProj - iProj; std::cout< TrivialPrec; - FlexibleGeneralisedMinimalResidual FGMRES(1.0e-14,1,TrivialPrec,1,false); // TODO: need to use GMRES as long as Mdag doesn't work on coarser levels (i.e., MdagM isn't hermitian) + ConjugateGradient CG(1.0e-2,10000); FineField noise(FineGrid); FineField Mn(FineGrid); @@ -222,9 +217,9 @@ namespace Grid { hermop.Op(noise,Mn); std::cout< "< - class CoarsenedMatrix : public SparseMatrixBase, 1 > > > > { + class CoarsenedMatrix : public SparseMatrixBase > > { public: + + typedef iVector siteVector; + typedef Lattice CoarseVector; + typedef Lattice > CoarseMatrix; - typedef typename CComplex::vector_type innerType; - typedef iScalar > > siteScalar; - typedef iScalar, 1 > > siteVector; - typedef iScalar, 1 > > siteMatrix; - typedef Lattice CoarseScalar; // used for inner products on fine field - typedef Lattice CoarseVector; - typedef Lattice CoarseMatrix; - - typedef Lattice FineField; + typedef Lattice< CComplex > CoarseScalar; // used for inner products on fine field + typedef Lattice FineField; //////////////////// // Data members @@ -303,50 +295,13 @@ namespace Grid { return norm2(out); }; - RealD Mdag (const CoarseVector &in, CoarseVector &out){ // TODO: get this correct + RealD Mdag (const CoarseVector &in, CoarseVector &out){ return M(in,out); }; - void Mdir(const CoarseVector &in, CoarseVector &out, int dir, int disp) { - - conformable(_grid,in._grid); - conformable(in._grid,out._grid); - - SimpleCompressor compressor; - Stencil.HaloExchange(in,compressor); - - auto point = [dir, disp](){ - if(dir == 0 and disp == 0) - return 8; - else - return (4 * dir + 1 - disp) / 2; - }(); - - parallel_for(int ss=0;ssoSites();ss++){ - siteVector res = zero; - siteVector nbr; - int ptype; - StencilEntry *SE; - - SE=Stencil.GetEntry(ptype,point,ss); - - if(SE->_is_local&&SE->_permute) { - permute(nbr,in._odata[SE->_offset],ptype); - } else if(SE->_is_local) { - nbr = in._odata[SE->_offset]; - } else { - nbr = Stencil.CommBuf()[SE->_offset]; - } - - res = res + A[point]._odata[ss]*nbr; - - vstream(out._odata[ss],res); - } - }; - - void Mdiag(const CoarseVector &in, CoarseVector &out) { - Mdir(in, out, 0, 0); // use the self coupling (= last) point of the stencil - }; + // Defer support for further coarsening for now + void Mdiag (const CoarseVector &in, CoarseVector &out){}; + void Mdir (const CoarseVector &in, CoarseVector &out,int dir, int disp){}; CoarsenedMatrix(GridCartesian &CoarseGrid) : @@ -372,9 +327,10 @@ namespace Grid { CoarseVector iProj(Grid()); CoarseVector oProj(Grid()); + CoarseScalar InnerProd(Grid()); // Orthogonalise the subblocks over the basis - Subspace.Orthogonalise(); + blockOrthogonalise(InnerProd,Subspace.subspace); // Compute the matrix elements of linop between this orthonormal // set of vectors. @@ -431,9 +387,9 @@ namespace Grid { parallel_for(int ss=0;ssoSites();ss++){ for(int j=0;j -inline void blockProject(Lattice &coarseData, +template +inline void blockProject(Lattice > &coarseData, const Lattice &fineData, const std::vector > &Basis) { @@ -90,8 +90,7 @@ inline void blockProject(Lattice &coarseData, int _ndimension = coarse->_ndimension; // checks - assert((Basis.size() != 0) && ((Basis.size() & 0x1) == 0)); - auto nbasis = Basis.size(); + assert( nbasis == Basis.size() ); subdivides(coarse,fine); for(int i=0;i &coarseData, PARALLEL_CRITICAL for(int i=0;i &ip,std::vector > } } -template -inline void blockPromote(const Lattice &coarseData, - Lattice &fineData, +template +inline void blockPromote(const Lattice > &coarseData, + Lattice &fineData, const std::vector > &Basis) { GridBase * fine = fineData._grid; @@ -296,9 +295,7 @@ inline void blockPromote(const Lattice &coarseData, int _ndimension = coarse->_ndimension; // checks - assert((Basis.size() != 0) && ((Basis.size() & 0x1) == 0)); - auto nbasis = Basis.size(); - + assert( nbasis == Basis.size() ); subdivides(coarse,fine); for(int i=0;i &coarseData, for(int d=0;d<_ndimension;d++) coor_c[d]=coor_f[d]/block_r[d]; Lexicographic::IndexFromCoor(coor_c,sc,coarse->_rdimensions); - // The temporary is necessary, since a pure instance of Grid::simd<...> is - // not a valid argument to operator+ with an iVector, we need an an iScalar - typename vobjC::tensor_reduced tmp; // iScalar>> -> iScalar>> for(int i=0;i