From e82ddcff5df342c4941a205f296af179387cd51a Mon Sep 17 00:00:00 2001 From: Peter Boyle Date: Thu, 7 Sep 2023 10:53:51 -0400 Subject: [PATCH] Working getting closer to HDCG but some low level engineering work still needed + MUCH work on optimisation --- Grid/algorithms/GeneralCoarsenedMatrix.h | 132 +++++++++++++++++++---- 1 file changed, 110 insertions(+), 22 deletions(-) diff --git a/Grid/algorithms/GeneralCoarsenedMatrix.h b/Grid/algorithms/GeneralCoarsenedMatrix.h index b8527cfb..4762bf73 100644 --- a/Grid/algorithms/GeneralCoarsenedMatrix.h +++ b/Grid/algorithms/GeneralCoarsenedMatrix.h @@ -62,7 +62,7 @@ public: // Need to worry about red-black now class NextToNearestStencilGeometry4D : public NonLocalStencilGeometry { public: - NextToNearestStencilGeometry4D(void) : NonLocalStencilGeometry(2) + NextToNearestStencilGeometry4D(void) : NonLocalStencilGeometry(1) { this->BuildShifts(); }; @@ -95,7 +95,7 @@ public: // Need to worry about red-black now class NextToNextToNextToNearestStencilGeometry4D : public NonLocalStencilGeometry { public: - NextToNextToNextToNearestStencilGeometry4D(void) : NonLocalStencilGeometry(4) + NextToNextToNextToNearestStencilGeometry4D(void) : NonLocalStencilGeometry(1) { this->BuildShifts(); }; @@ -142,7 +142,7 @@ public: }; class NextToNearestStencilGeometry5D : public NonLocalStencilGeometry { public: - NextToNearestStencilGeometry5D(void) : NonLocalStencilGeometry(2) + NextToNearestStencilGeometry5D(void) : NonLocalStencilGeometry(1) { this->BuildShifts(); }; @@ -175,7 +175,7 @@ public: // Need to worry about red-black now class NextToNextToNextToNearestStencilGeometry5D : public NonLocalStencilGeometry { public: - NextToNextToNextToNearestStencilGeometry5D(void) : NonLocalStencilGeometry(4) + NextToNextToNextToNearestStencilGeometry5D(void) : NonLocalStencilGeometry(1) { this->BuildShifts(); }; @@ -239,7 +239,7 @@ public: // Data members //////////////////// int hermitian; - GridCartesian * _FineGrid; + GridBase * _FineGrid; GridCartesian * _CoarseGrid; NonLocalStencilGeometry &geom; PaddedCell Cell; @@ -251,11 +251,11 @@ public: /////////////////////// // Interface /////////////////////// - GridCartesian * Grid(void) { return _FineGrid; }; // this is all the linalg routines need to know - GridCartesian * FineGrid(void) { return _FineGrid; }; // this is all the linalg routines need to know + GridBase * Grid(void) { return _FineGrid; }; // this is all the linalg routines need to know + GridBase * FineGrid(void) { return _FineGrid; }; // this is all the linalg routines need to know GridCartesian * CoarseGrid(void) { return _CoarseGrid; }; // this is all the linalg routines need to know - GeneralCoarsenedMatrix(NonLocalStencilGeometry &_geom,GridCartesian *FineGrid, GridCartesian * CoarseGrid) + GeneralCoarsenedMatrix(NonLocalStencilGeometry &_geom,GridBase *FineGrid, GridCartesian * CoarseGrid) : geom(_geom), _FineGrid(FineGrid), _CoarseGrid(CoarseGrid), @@ -263,6 +263,20 @@ public: Cell(_geom.Depth(),_CoarseGrid), Stencil(Cell.grids.back(),geom.shifts) { + { + int npoint = _geom.npoint; + StencilEntry *SE; + autoView( Stencil_v , Stencil, AcceleratorRead); + int osites=Stencil.Grid()->oSites(); + for(int ss=0;ss_offset; + assert( o< osites); + } + } + } + _A.resize(geom.npoint,CoarseGrid); _Adag.resize(geom.npoint,CoarseGrid); } @@ -280,9 +294,9 @@ public: conformable(in.Grid(),out.Grid()); out.Checkerboard() = in.Checkerboard(); CoarseVector tin=in; - std::cout << "Calling Exchange"<oSites(); for(int point=0;point_offset; + assert( o < osites); // gpermute etc.. nbr = in_v[o]; - assert( o< osites); gpermute(nbr,SE->_permute); for(int bb=0;bb & Subspace) { // Create a random - GridCartesian *grid = FineGrid(); + GridBase *grid = FineGrid(); FineField MbV(grid); FineField tmp(grid); FineField f_src(grid); @@ -370,8 +382,12 @@ public: void CoarsenOperator(LinearOperatorBase > &linop, Aggregation & Subspace) { + RealD tproj=0.0; + RealD tpick=0.0; + RealD tmat=0.0; + RealD tpeek=0.0; std::cout << GridLogMessage<< "CoarsenMatrix "<< std::endl; - GridCartesian *grid = FineGrid(); + GridBase *grid = FineGrid(); // Orthogonalise the subblocks over the basis CoarseScalar InnerProd(CoarseGrid()); for(int b=0;bgSites() ;bidx++){ Coordinate bcoor; CoarseGrid()->GlobalIndexToGlobalCoor(bidx,bcoor); - std::cout << GridLogMessage<< "CoarsenMatrix block "<< bcoor << std::endl; + for(int b=0;bGlobalDimensions()[mu]; scoor[mu] = (bcoor[mu] - geom.shifts[p][mu] + L) % L; // Modulo arithmetic } + // Flip to peekLocalSite + // Flip to pokeLocalSite auto ip = peekSite(coarseInner,scoor); - std::cout << "A["< > &linop, + Aggregation & Subspace) + { + std::cout << GridLogMessage<< "CoarsenMatrixColoured "<< std::endl; + GridBase *grid = FineGrid(); + + ///////////////////////////////////////////////////////////// + // Orthogonalise the subblocks over the basis + ///////////////////////////////////////////////////////////// + CoarseScalar InnerProd(CoarseGrid()); + blockOrthogonalise(InnerProd,Subspace.subspace); + + const int npoint = geom.npoint + + // Now compute the matrix elements of linop between this orthonormal + // set of vectors. + FineField bV(grid); + FineField MbV(grid); + FineField tmp(grid); + CoarseVector coarseInner(CoarseGrid()); + Coordinate clatt = CoarseGrid()->GlobalDimensions(); + Coordinate steclatt = CoarseGrid()->GlobalDimensions(); + for(int v=0;v dmom = + for(int mu=0;mu<5;mu++){ + dmom[mu] = imom[mu]*2*M_PI/global_nb[mu]; + } + /* + * Here, k,l index which possible momentum/shift within the N-points connected by MdagM. + * Matrix index i is mapped to this shift via + * GetDelta(imat_ball_idx[i]). + * + * conj(phases[block]) proj[k][ block*Nvec+j ] = \sum_{l in ball} e^{i q_k . delta_l} < phi_{block,j} | MdagM | phi_{(block+delta_l),i} > + * = \sum_{l in ball} e^{iqk.delta_l} A_ji^{b.b+l} + * = M_{kl} A_ji^{b.b+l} + * + * Must assemble and invert matrix M_k,l = e^[i q_k . delta_l] + * + * Where q_k = delta_k . (2*M_PI/global_nb[mu]) + * + * Then A{ji}^{b,b+l} = M^{-1}_{lm} ComputeProj_{m,b,i,j} + */ + } + } + for(int p=0;p &out){assert(0);};