Enable interafaces with both coarse and fine being either single higher dim multiRHS lattice OR vector of lower dim fields

This commit is contained in:
Peter Boyle
2026-08-19 19:33:55 -04:00
parent 9c4a61850e
commit c10d451f4d
+559 -47
View File
@@ -61,6 +61,17 @@ public:
uint64_t coarse_vol;
uint64_t words;
////////////////////////////////////////////////////////////////////////////
// Blocking geometry in full local coordinates. Addressing in lSites rather
// than (lane,oSite) lets the fine and coarse spaces carry different SIMD
// layouts, so an unvectorised coarse space may block a vectorised fine one.
////////////////////////////////////////////////////////////////////////////
Coordinate fine_ldimensions;
Coordinate coarse_ldimensions;
Coordinate block_ldimensions;
Coordinate fine_simd;
Coordinate coarse_simd;
// Row major layout "C" order:
// BLAS_V[coarse_vol][nbasis][block_vol][words]
// BLAS_F[coarse_vol][nrhs][block_vol][words]
@@ -120,9 +131,26 @@ public:
fine_vol = fine_grid->lSites();
coarse_vol = coarse_grid->lSites();
block_vol = fine_vol/coarse_vol;
words = sizeof(scalar_object)/sizeof(scalar);
int nd = coarse_grid->_ndimension;
GRID_ASSERT(fine_grid->_ndimension == nd);
fine_ldimensions.resize(nd);
coarse_ldimensions.resize(nd);
block_ldimensions.resize(nd);
fine_simd = fine_grid->_simd_layout;
coarse_simd = coarse_grid->_simd_layout;
for(int d=0;d<nd;d++){
GRID_ASSERT(fine_grid->_processors[d] == coarse_grid->_processors[d]);
fine_ldimensions [d] = fine_grid->_rdimensions [d]*fine_grid->_simd_layout [d];
coarse_ldimensions[d] = coarse_grid->_rdimensions[d]*coarse_grid->_simd_layout[d];
block_ldimensions [d] = fine_ldimensions[d]/coarse_ldimensions[d];
GRID_ASSERT(block_ldimensions[d]*coarse_ldimensions[d] == fine_ldimensions[d]);
}
GRID_ASSERT(block_vol == fine_vol/coarse_vol);
BLAS_V.resize (fine_vol * words * nbasis );
}
void ImportFineGridVectors(std::vector <Field > &vecs, deviceVector<scalar> &blas)
@@ -133,22 +161,16 @@ public:
GRID_ASSERT(vecs[0].Grid()==fine_grid);
subdivides(coarse_grid,fine_grid); // require they map
int _ndimension = coarse_grid->_ndimension;
GRID_ASSERT(block_vol == fine_grid->oSites() / coarse_grid->oSites());
Coordinate block_r (_ndimension);
for(int d=0 ; d<_ndimension;d++){
block_r[d] = fine_grid->_rdimensions[d] / coarse_grid->_rdimensions[d];
}
uint64_t sz = blas.size();
acceleratorMemSet(&blas[0],0,blas.size()*sizeof(scalar));
Coordinate fine_rdimensions = fine_grid->_rdimensions;
Coordinate coarse_rdimensions = coarse_grid->_rdimensions;
Coordinate coarse_l = coarse_ldimensions;
Coordinate block_l = block_ldimensions;
Coordinate fsimd = fine_simd;
int64_t bv= block_vol;
for(int v=0;v<vecs.size();v++){
@@ -162,9 +184,7 @@ public:
// loop over fine sites
const int Nsimd = vobj::Nsimd();
// std::cout << "sz "<<sz<<std::endl;
// std::cout << "prod "<<Nsimd * coarse_grid->oSites() * block_vol * nvec * words<<std::endl;
GRID_ASSERT(sz == Nsimd * coarse_grid->oSites() * block_vol * nvec * words);
GRID_ASSERT(sz == coarse_vol * block_vol * nvec * words);
uint64_t lwords= words; // local variable for copy in to GPU
accelerator_for(sf,osites,Nsimd,{
#ifdef GRID_SIMT
@@ -175,27 +195,27 @@ public:
#endif
// One thread per fine site
Coordinate coor_f(_ndimension);
Coordinate coor_l(_ndimension);
Coordinate coor_b(_ndimension);
Coordinate coor_c(_ndimension);
// Fine site to fine coor
// Fine (oSite,lane) to full local coor
Lexicographic::CoorFromIndex(coor_f,sf,fine_rdimensions);
Lexicographic::CoorFromIndex(coor_l,lane,fsimd);
for(int d=0;d<_ndimension;d++) coor_f[d] += fine_rdimensions[d]*coor_l[d];
for(int d=0;d<_ndimension;d++) coor_b[d] = coor_f[d]%block_l[d];
for(int d=0;d<_ndimension;d++) coor_c[d] = coor_f[d]/block_l[d];
for(int d=0;d<_ndimension;d++) coor_b[d] = coor_f[d]%block_r[d];
for(int d=0;d<_ndimension;d++) coor_c[d] = coor_f[d]/block_r[d];
int sc;// coarse site
int sb;// block site
Lexicographic::IndexFromCoor(coor_c,sc,coarse_rdimensions);
Lexicographic::IndexFromCoor(coor_b,sb,block_r);
Lexicographic::IndexFromCoor(coor_c,sc,coarse_l);
Lexicographic::IndexFromCoor(coor_b,sb,block_l);
scalar_object data = extractLane(lane,fineData[sf]);
// BLAS layout address calculation
// words * block_vol * nbasis x coarse_vol
// coarse oSite x block vole x lanes
int64_t site = (lane*osites + sc*bv)*nvec
+ v*bv
// BLAS_F[coarse_vol][nvec][block_vol][words]
int64_t site = (sc*nvec + v)*bv
+ sb;
// GRID_ASSERT(site*lwords<sz);
@@ -213,6 +233,318 @@ public:
// std::cout << " BlockProjector imported vector"<<v<<std::endl;
}
}
////////////////////////////////////////////////////////////////////////////
// Import direct from multiRHS fine order, avoiding the unpack to a vector
// of single RHS fields.
//
// fine mrhs grid : rhs is dimension 0, undistributed, unvectorised
// Grid order : F[fine_vol][nrhs][words]
// BLAS order : BLAS_F[lane][coarse_vol][nrhs][block_vol][words]
//
// The gather of block_vol from fine_vol, and the transpose of nrhs against
// block_vol, are the irreducible part: this is not the identity.
////////////////////////////////////////////////////////////////////////////
void ImportFineGridMrhsVectors(Field &vec_mrhs, deviceVector<scalar> &blas)
{
typedef typename Field::vector_object vobj;
GridBase *fine_mrhs_grid = vec_mrhs.Grid();
int _ndimension = coarse_grid->_ndimension;
GRID_ASSERT(fine_mrhs_grid->_ndimension == _ndimension+1);
GRID_ASSERT(fine_mrhs_grid->_simd_layout[0] == 1);
GRID_ASSERT(fine_mrhs_grid->_processors[0] == 1);
for(int d=0;d<_ndimension;d++){
GRID_ASSERT(fine_mrhs_grid->_rdimensions[d+1] == fine_grid->_rdimensions[d]);
GRID_ASSERT(fine_mrhs_grid->_simd_layout[d+1] == fine_grid->_simd_layout[d]);
}
int nvec = fine_mrhs_grid->_rdimensions[0]; // nrhs
uint64_t sz = blas.size();
acceleratorMemSet(&blas[0],0,blas.size()*sizeof(scalar));
Coordinate fine_mrhs_rdimensions = fine_mrhs_grid->_rdimensions;
Coordinate fine_rdimensions = fine_grid->_rdimensions;
Coordinate coarse_l = coarse_ldimensions;
Coordinate block_l = block_ldimensions;
Coordinate fsimd = fine_simd;
int64_t bv= block_vol;
autoView( fineData , vec_mrhs, AcceleratorRead);
auto blasData_p = &blas[0];
auto fineData_p = &fineData[0];
int64_t osites = fine_grid->oSites(); // D dimensional
int64_t osites_hi = fine_mrhs_grid->oSites(); // nvec * osites
const int Nsimd = vobj::Nsimd();
GRID_ASSERT(sz == coarse_vol * block_vol * nvec * words);
uint64_t lwords= words;
int64_t lnvec = nvec;
accelerator_for(sfr,osites_hi,Nsimd,{
#ifdef GRID_SIMT
{
int lane=acceleratorSIMTlane(Nsimd); // buffer lane
#else
for(int lane=0;lane<Nsimd;lane++) {
#endif
Coordinate coor_hi(_ndimension+1);
Coordinate coor_f(_ndimension);
Coordinate coor_l(_ndimension);
Coordinate coor_b(_ndimension);
Coordinate coor_c(_ndimension);
// rhs is dimension 0 of the D+1 grid
Lexicographic::CoorFromIndex(coor_hi,sfr,fine_mrhs_rdimensions);
int v = coor_hi[0];
for(int d=0;d<_ndimension;d++) coor_f[d] = coor_hi[d+1];
Lexicographic::CoorFromIndex(coor_l,lane,fsimd);
for(int d=0;d<_ndimension;d++) coor_f[d] += fine_rdimensions[d]*coor_l[d];
for(int d=0;d<_ndimension;d++) coor_b[d] = coor_f[d]%block_l[d];
for(int d=0;d<_ndimension;d++) coor_c[d] = coor_f[d]/block_l[d];
int sc;// coarse site
int sb;// block site
Lexicographic::IndexFromCoor(coor_c,sc,coarse_l);
Lexicographic::IndexFromCoor(coor_b,sb,block_l);
scalar_object data = extractLane(lane,fineData[sfr]);
int64_t site = (sc*lnvec + v)*bv
+ sb;
scalar_object * ptr = (scalar_object *)&blasData_p[site*lwords];
*ptr = data;
#ifdef GRID_SIMT
}
#else
}
#endif
});
}
////////////////////////////////////////////////////////////////////////////
// Export direct to multiRHS coarse order.
//
// Grid order : C[coarse_vol][nrhs][nbasis]
// BLAS order : BLAS_C[lane][coarse_vol][nrhs][nbasis]
//
// At Nsimd()==1 these are the same sequence of addresses.
////////////////////////////////////////////////////////////////////////////
template<class vobj>
void ExportCoarseGridMrhsVectors(Lattice<vobj> &vec_mrhs, deviceVector<scalar> &blas)
{
typedef typename vobj::scalar_object coarse_scalar_object;
GridBase *coarse_mrhs_grid = vec_mrhs.Grid();
int _ndimension = coarse_grid->_ndimension;
GRID_ASSERT(coarse_mrhs_grid->_ndimension == _ndimension+1);
GRID_ASSERT(coarse_mrhs_grid->_simd_layout[0] == 1);
GRID_ASSERT(coarse_mrhs_grid->_processors[0] == 1);
for(int d=0;d<_ndimension;d++){
GRID_ASSERT(coarse_mrhs_grid->_rdimensions[d+1] == coarse_grid->_rdimensions[d]);
GRID_ASSERT(coarse_mrhs_grid->_simd_layout[d+1] == coarse_grid->_simd_layout[d]);
}
int nvec = coarse_mrhs_grid->_rdimensions[0]; // nrhs
Coordinate coarse_mrhs_rdimensions = coarse_mrhs_grid->_rdimensions;
Coordinate coarse_rdimensions = coarse_grid->_rdimensions;
Coordinate coarse_l = coarse_ldimensions;
Coordinate csimd = coarse_simd;
autoView( coarseData , vec_mrhs, AcceleratorWrite);
auto blasData_p = &blas[0];
auto coarseData_p = &coarseData[0];
int64_t osites = coarse_grid->oSites(); // D dimensional
int64_t osites_hi = coarse_mrhs_grid->oSites(); // nvec * osites
const int Nsimd = vobj::Nsimd();
uint64_t cwords=sizeof(typename vobj::scalar_object)/sizeof(scalar);
GRID_ASSERT(cwords==nbasis);
int64_t lnvec = nvec;
accelerator_for(scr,osites_hi,Nsimd,{
#ifdef GRID_SIMT
{
int lane=acceleratorSIMTlane(Nsimd); // buffer lane
#else
for(int lane=0;lane<Nsimd;lane++) {
#endif
Coordinate coor_hi(_ndimension+1);
Coordinate coor_l(_ndimension);
Coordinate coor_c(_ndimension);
Lexicographic::CoorFromIndex(coor_hi,scr,coarse_mrhs_rdimensions);
int v = coor_hi[0];
for(int d=0;d<_ndimension;d++) coor_c[d] = coor_hi[d+1];
Lexicographic::CoorFromIndex(coor_l,lane,csimd);
for(int d=0;d<_ndimension;d++) coor_c[d] += coarse_rdimensions[d]*coor_l[d];
int sc;
Lexicographic::IndexFromCoor(coor_c,sc,coarse_l);
int64_t blas_site = (sc*lnvec + v)*cwords;
coarse_scalar_object * ptr = (coarse_scalar_object *)&blasData_p[blas_site];
coarse_scalar_object data = *ptr;
insertLane(lane,coarseData[scr],data);
#ifdef GRID_SIMT
}
#else
}
#endif
});
}
////////////////////////////////////////////////////////////////////////////
// Reverse directions: BLAS_F -> fine mrhs field, coarse mrhs field -> BLAS_C
////////////////////////////////////////////////////////////////////////////
void ExportFineGridMrhsVectors(Field &vec_mrhs, deviceVector<scalar> &blas)
{
typedef typename Field::vector_object vobj;
GridBase *fine_mrhs_grid = vec_mrhs.Grid();
int _ndimension = coarse_grid->_ndimension;
GRID_ASSERT(fine_mrhs_grid->_ndimension == _ndimension+1);
GRID_ASSERT(fine_mrhs_grid->_simd_layout[0] == 1);
for(int d=0;d<_ndimension;d++){
GRID_ASSERT(fine_mrhs_grid->_rdimensions[d+1] == fine_grid->_rdimensions[d]);
GRID_ASSERT(fine_mrhs_grid->_simd_layout[d+1] == fine_grid->_simd_layout[d]);
}
int nvec = fine_mrhs_grid->_rdimensions[0];
Coordinate fine_mrhs_rdimensions = fine_mrhs_grid->_rdimensions;
Coordinate fine_rdimensions = fine_grid->_rdimensions;
Coordinate coarse_l = coarse_ldimensions;
Coordinate block_l = block_ldimensions;
Coordinate fsimd = fine_simd;
int64_t bv= block_vol;
autoView( fineData , vec_mrhs, AcceleratorWrite);
auto blasData_p = &blas[0];
auto fineData_p = &fineData[0];
int64_t osites = fine_grid->oSites();
int64_t osites_hi = fine_mrhs_grid->oSites();
const int Nsimd = vobj::Nsimd();
uint64_t lwords= words;
int64_t lnvec = nvec;
accelerator_for(sfr,osites_hi,Nsimd,{
#ifdef GRID_SIMT
{
int lane=acceleratorSIMTlane(Nsimd);
#else
for(int lane=0;lane<Nsimd;lane++) {
#endif
Coordinate coor_hi(_ndimension+1);
Coordinate coor_f(_ndimension);
Coordinate coor_l(_ndimension);
Coordinate coor_b(_ndimension);
Coordinate coor_c(_ndimension);
Lexicographic::CoorFromIndex(coor_hi,sfr,fine_mrhs_rdimensions);
int v = coor_hi[0];
for(int d=0;d<_ndimension;d++) coor_f[d] = coor_hi[d+1];
Lexicographic::CoorFromIndex(coor_l,lane,fsimd);
for(int d=0;d<_ndimension;d++) coor_f[d] += fine_rdimensions[d]*coor_l[d];
for(int d=0;d<_ndimension;d++) coor_b[d] = coor_f[d]%block_l[d];
for(int d=0;d<_ndimension;d++) coor_c[d] = coor_f[d]/block_l[d];
int sc,sb;
Lexicographic::IndexFromCoor(coor_c,sc,coarse_l);
Lexicographic::IndexFromCoor(coor_b,sb,block_l);
int64_t site = (sc*lnvec + v)*bv
+ sb;
scalar_object * ptr = (scalar_object *)&blasData_p[site*lwords];
scalar_object data = *ptr;
insertLane(lane,fineData[sfr],data);
#ifdef GRID_SIMT
}
#else
}
#endif
});
}
template<class vobj>
void ImportCoarseGridMrhsVectors(Lattice<vobj> &vec_mrhs, deviceVector<scalar> &blas)
{
typedef typename vobj::scalar_object coarse_scalar_object;
GridBase *coarse_mrhs_grid = vec_mrhs.Grid();
int _ndimension = coarse_grid->_ndimension;
GRID_ASSERT(coarse_mrhs_grid->_ndimension == _ndimension+1);
GRID_ASSERT(coarse_mrhs_grid->_simd_layout[0] == 1);
for(int d=0;d<_ndimension;d++){
GRID_ASSERT(coarse_mrhs_grid->_rdimensions[d+1] == coarse_grid->_rdimensions[d]);
GRID_ASSERT(coarse_mrhs_grid->_simd_layout[d+1] == coarse_grid->_simd_layout[d]);
}
int nvec = coarse_mrhs_grid->_rdimensions[0];
Coordinate coarse_mrhs_rdimensions = coarse_mrhs_grid->_rdimensions;
Coordinate coarse_rdimensions = coarse_grid->_rdimensions;
Coordinate coarse_l = coarse_ldimensions;
Coordinate csimd = coarse_simd;
autoView( coarseData , vec_mrhs, AcceleratorRead);
auto blasData_p = &blas[0];
auto coarseData_p = &coarseData[0];
int64_t osites = coarse_grid->oSites();
int64_t osites_hi = coarse_mrhs_grid->oSites();
const int Nsimd = vobj::Nsimd();
uint64_t cwords=sizeof(typename vobj::scalar_object)/sizeof(scalar);
GRID_ASSERT(cwords==nbasis);
int64_t lnvec = nvec;
accelerator_for(scr,osites_hi,Nsimd,{
#ifdef GRID_SIMT
{
int lane=acceleratorSIMTlane(Nsimd);
#else
for(int lane=0;lane<Nsimd;lane++) {
#endif
Coordinate coor_hi(_ndimension+1);
Coordinate coor_l(_ndimension);
Coordinate coor_c(_ndimension);
Lexicographic::CoorFromIndex(coor_hi,scr,coarse_mrhs_rdimensions);
int v = coor_hi[0];
for(int d=0;d<_ndimension;d++) coor_c[d] = coor_hi[d+1];
Lexicographic::CoorFromIndex(coor_l,lane,csimd);
for(int d=0;d<_ndimension;d++) coor_c[d] += coarse_rdimensions[d]*coor_l[d];
int sc;
Lexicographic::IndexFromCoor(coor_c,sc,coarse_l);
coarse_scalar_object data = extractLane(lane,coarseData[scr]);
int64_t blas_site = (sc*lnvec + v)*cwords;
coarse_scalar_object * ptr = (coarse_scalar_object *)&blasData_p[blas_site];
*ptr = data;
#ifdef GRID_SIMT
}
#else
}
#endif
});
}
void ExportFineGridVectors(std::vector <Field> &vecs, deviceVector<scalar> &blas)
{
typedef typename Field::vector_object vobj;
@@ -221,17 +553,12 @@ public:
GRID_ASSERT(vecs[0].Grid()==fine_grid);
subdivides(coarse_grid,fine_grid); // require they map
int _ndimension = coarse_grid->_ndimension;
GRID_ASSERT(block_vol == fine_grid->oSites() / coarse_grid->oSites());
Coordinate block_r (_ndimension);
for(int d=0 ; d<_ndimension;d++){
block_r[d] = fine_grid->_rdimensions[d] / coarse_grid->_rdimensions[d];
}
Coordinate fine_rdimensions = fine_grid->_rdimensions;
Coordinate coarse_rdimensions = coarse_grid->_rdimensions;
Coordinate coarse_l = coarse_ldimensions;
Coordinate block_l = block_ldimensions;
Coordinate fsimd = fine_simd;
// std::cout << " export fine Blas norm "<<blasNorm2(blas)<<std::endl;
@@ -259,23 +586,24 @@ public:
#endif
// One thread per fine site
Coordinate coor_f(_ndimension);
Coordinate coor_l(_ndimension);
Coordinate coor_b(_ndimension);
Coordinate coor_c(_ndimension);
Lexicographic::CoorFromIndex(coor_f,sf,fine_rdimensions);
Lexicographic::CoorFromIndex(coor_l,lane,fsimd);
for(int d=0;d<_ndimension;d++) coor_f[d] += fine_rdimensions[d]*coor_l[d];
for(int d=0;d<_ndimension;d++) coor_b[d] = coor_f[d]%block_l[d];
for(int d=0;d<_ndimension;d++) coor_c[d] = coor_f[d]/block_l[d];
for(int d=0;d<_ndimension;d++) coor_b[d] = coor_f[d]%block_r[d];
for(int d=0;d<_ndimension;d++) coor_c[d] = coor_f[d]/block_r[d];
int sc;
int sb;
Lexicographic::IndexFromCoor(coor_c,sc,coarse_rdimensions);
Lexicographic::IndexFromCoor(coor_b,sb,block_r);
Lexicographic::IndexFromCoor(coor_c,sc,coarse_l);
Lexicographic::IndexFromCoor(coor_b,sb,block_l);
// BLAS layout address calculation
// words * block_vol * nbasis x coarse_vol
int64_t site = (lane*osites + sc*bv)*nvec
+ v*bv
// BLAS_F[coarse_vol][nvec][block_vol][words]
int64_t site = (sc*nvec + v)*bv
+ sb;
scalar_object * ptr = (scalar_object *)&blasData_p[site*lwords];
@@ -306,7 +634,9 @@ public:
uint64_t sz = blas.size();
Coordinate coarse_rdimensions = coarse_grid->_rdimensions;
Coordinate coarse_l = coarse_ldimensions;
Coordinate csimd = coarse_simd;
for(int v=0;v<vecs.size();v++){
// std::cout << " BlockProjector importing coarse vector"<<v<<" "<<norm2(vecs[v])<<std::endl;
@@ -330,8 +660,17 @@ public:
for(int lane=0;lane<Nsimd;lane++) {
#endif
// C_br per site
int64_t blas_site = (lane*osites + sc)*nvec*cwords + v*cwords;
Coordinate coor_c(_ndimension);
Coordinate coor_l(_ndimension);
Lexicographic::CoorFromIndex(coor_c,sc,coarse_rdimensions);
Lexicographic::CoorFromIndex(coor_l,lane,csimd);
for(int d=0;d<_ndimension;d++) coor_c[d] += coarse_rdimensions[d]*coor_l[d];
int scl;
Lexicographic::IndexFromCoor(coor_c,scl,coarse_l);
int64_t blas_site = (scl*nvec + v)*cwords;
coarse_scalar_object data = extractLane(lane,coarseData[sc]);
coarse_scalar_object * ptr = (coarse_scalar_object *)&blasData_p[blas_site];
@@ -360,7 +699,9 @@ public:
uint64_t sz = blas.size();
Coordinate coarse_rdimensions = coarse_grid->_rdimensions;
Coordinate coarse_l = coarse_ldimensions;
Coordinate csimd = coarse_simd;
// std::cout << " export coarsee Blas norm "<<blasNorm2(blas)<<std::endl;
for(int v=0;v<vecs.size();v++){
@@ -385,7 +726,16 @@ public:
#else
for(int lane=0;lane<Nsimd;lane++) {
#endif
int64_t blas_site = (lane*osites + sc)*nvec*cwords + v*cwords;
Coordinate coor_c(_ndimension);
Coordinate coor_l(_ndimension);
Lexicographic::CoorFromIndex(coor_c,sc,coarse_rdimensions);
Lexicographic::CoorFromIndex(coor_l,lane,csimd);
for(int d=0;d<_ndimension;d++) coor_c[d] += coarse_rdimensions[d]*coor_l[d];
int scl;
Lexicographic::IndexFromCoor(coor_c,scl,coarse_l);
int64_t blas_site = (scl*nvec + v)*cwords;
coarse_scalar_object * ptr = (coarse_scalar_object *)&blasData_p[blas_site];
coarse_scalar_object data = *ptr;
insertLane(lane,coarseData[sc],data);
@@ -508,6 +858,168 @@ public:
ExportFineGridVectors(fine, BLAS_F);
// std::cout << " exported "<<std::endl;
}
////////////////////////////////////////////////////////////////////////////
// multiRHS ordered interfaces. The GEMM is identical; only the import and
// export differ, so those are the whole of the layout question.
////////////////////////////////////////////////////////////////////////////
template<class cobj>
void blockProject(Field &fine_mrhs,Lattice<cobj> &coarse_mrhs)
{
int nrhs = fine_mrhs.Grid()->_rdimensions[0];
int _nbasis = sizeof(typename cobj::scalar_object)/sizeof(scalar);
GRID_ASSERT(nbasis==_nbasis);
GRID_ASSERT(coarse_mrhs.Grid()->_rdimensions[0]==nrhs);
BLAS_F.resize (fine_vol * words * nrhs );
BLAS_C.resize (coarse_vol * nbasis * nrhs );
ImportFineGridMrhsVectors(fine_mrhs,BLAS_F);
ProjectBLAS(nrhs);
ExportCoarseGridMrhsVectors(coarse_mrhs,BLAS_C);
}
template<class cobj>
void blockPromote(Field &fine_mrhs,Lattice<cobj> &coarse_mrhs)
{
int nrhs = fine_mrhs.Grid()->_rdimensions[0];
int _nbasis = sizeof(typename cobj::scalar_object)/sizeof(scalar);
GRID_ASSERT(nbasis==_nbasis);
GRID_ASSERT(coarse_mrhs.Grid()->_rdimensions[0]==nrhs);
BLAS_F.resize (fine_vol * words * nrhs );
BLAS_C.resize (coarse_vol * nbasis * nrhs );
ImportCoarseGridMrhsVectors(coarse_mrhs,BLAS_C);
PromoteBLAS(nrhs);
ExportFineGridMrhsVectors(fine_mrhs,BLAS_F);
}
////////////////////////////////////////////////////////////////////////////
// Mixed orderings. A single RHS fine operator produces a vector of fine
// fields with no packing; the coarse side is still wanted in mrhs order.
////////////////////////////////////////////////////////////////////////////
template<class cobj>
void blockProject(std::vector<Field> &fine,Lattice<cobj> &coarse_mrhs)
{
int nrhs = fine.size();
int _nbasis = sizeof(typename cobj::scalar_object)/sizeof(scalar);
GRID_ASSERT(nbasis==_nbasis);
GRID_ASSERT(coarse_mrhs.Grid()->_rdimensions[0]==nrhs);
BLAS_F.resize (fine_vol * words * nrhs );
BLAS_C.resize (coarse_vol * nbasis * nrhs );
ImportFineGridVectors(fine,BLAS_F);
ProjectBLAS(nrhs);
ExportCoarseGridMrhsVectors(coarse_mrhs,BLAS_C);
}
template<class cobj>
void blockProject(Field &fine_mrhs,std::vector< Lattice<cobj> > &coarse)
{
int nrhs = fine_mrhs.Grid()->_rdimensions[0];
int _nbasis = sizeof(typename cobj::scalar_object)/sizeof(scalar);
GRID_ASSERT(nbasis==_nbasis);
GRID_ASSERT(coarse.size()==nrhs);
BLAS_F.resize (fine_vol * words * nrhs );
BLAS_C.resize (coarse_vol * nbasis * nrhs );
ImportFineGridMrhsVectors(fine_mrhs,BLAS_F);
ProjectBLAS(nrhs);
ExportCoarseGridVectors(coarse,BLAS_C);
}
template<class cobj>
void blockPromote(std::vector<Field> &fine,Lattice<cobj> &coarse_mrhs)
{
int nrhs = fine.size();
int _nbasis = sizeof(typename cobj::scalar_object)/sizeof(scalar);
GRID_ASSERT(nbasis==_nbasis);
GRID_ASSERT(coarse_mrhs.Grid()->_rdimensions[0]==nrhs);
BLAS_F.resize (fine_vol * words * nrhs );
BLAS_C.resize (coarse_vol * nbasis * nrhs );
ImportCoarseGridMrhsVectors(coarse_mrhs,BLAS_C);
PromoteBLAS(nrhs);
ExportFineGridVectors(fine,BLAS_F);
}
template<class cobj>
void blockPromote(Field &fine_mrhs,std::vector< Lattice<cobj> > &coarse)
{
int nrhs = fine_mrhs.Grid()->_rdimensions[0];
int _nbasis = sizeof(typename cobj::scalar_object)/sizeof(scalar);
GRID_ASSERT(nbasis==_nbasis);
GRID_ASSERT(coarse.size()==nrhs);
BLAS_F.resize (fine_vol * words * nrhs );
BLAS_C.resize (coarse_vol * nbasis * nrhs );
ImportCoarseGridVectors(coarse,BLAS_C);
PromoteBLAS(nrhs);
ExportFineGridMrhsVectors(fine_mrhs,BLAS_F);
}
////////////////////////////////////////////////////////////////////////////
// Pointer tables and the batched GEMM, shared by both orderings
////////////////////////////////////////////////////////////////////////////
void BLASPointers(int nrhs,
deviceVector<scalar *> &Vd,
deviceVector<scalar *> &Fd,
deviceVector<scalar *> &Cd)
{
for(int c=0;c<coarse_vol;c++){
scalar * Vh = & BLAS_V[c*nbasis*block_vol*words];
scalar * Fh = & BLAS_F[c*nrhs*block_vol*words];
scalar * Ch = & BLAS_C[c*nrhs*nbasis];
acceleratorPut(Vd[c],Vh);
acceleratorPut(Fd[c],Fh);
acceleratorPut(Cd[c],Ch);
}
}
// C_br = V^dag F
void ProjectBLAS(int nrhs)
{
deviceVector<scalar *> Vd(coarse_vol);
deviceVector<scalar *> Fd(coarse_vol);
deviceVector<scalar *> Cd(coarse_vol);
BLASPointers(nrhs,Vd,Fd,Cd);
GridBLAS BLAS;
int64_t vw = block_vol * words;
BLAS.gemmBatched(GridBLAS_OP_C,GridBLAS_OP_N,
nbasis,nrhs,vw,
scalar(1.0),
Vd,
Fd,
scalar(0.0), // wipe out C
Cd);
BLAS.synchronise();
}
// F_xr = Vxb Cbr
void PromoteBLAS(int nrhs)
{
deviceVector<scalar *> Vd(coarse_vol);
deviceVector<scalar *> Fd(coarse_vol);
deviceVector<scalar *> Cd(coarse_vol);
BLASPointers(nrhs,Vd,Fd,Cd);
GridBLAS BLAS;
int64_t vw = block_vol * words;
BLAS.gemmBatched(GridBLAS_OP_N,GridBLAS_OP_N,
vw,nrhs,nbasis,
scalar(1.0),
Vd,
Cd,
scalar(0.0), // wipe out F
Fd);
BLAS.synchronise();
}
};
NAMESPACE_END(Grid);