2017-04-09 15:41:04 +01:00
|
|
|
/*************************************************************************************
|
2016-01-02 14:51:32 +00:00
|
|
|
|
|
|
|
Grid physics library, www.github.com/paboyle/Grid
|
|
|
|
|
|
|
|
Source file: ./lib/lattice/Lattice_transfer.h
|
|
|
|
|
|
|
|
Copyright (C) 2015
|
|
|
|
|
|
|
|
Author: Peter Boyle <paboyle@ph.ed.ac.uk>
|
|
|
|
|
|
|
|
This program is free software; you can redistribute it and/or modify
|
|
|
|
it under the terms of the GNU General Public License as published by
|
|
|
|
the Free Software Foundation; either version 2 of the License, or
|
|
|
|
(at your option) any later version.
|
|
|
|
|
|
|
|
This program is distributed in the hope that it will be useful,
|
|
|
|
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
|
|
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
|
|
|
GNU General Public License for more details.
|
|
|
|
|
|
|
|
You should have received a copy of the GNU General Public License along
|
|
|
|
with this program; if not, write to the Free Software Foundation, Inc.,
|
|
|
|
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
|
|
|
|
|
|
|
|
See the full license in the file "LICENSE" in the top level distribution directory
|
|
|
|
*************************************************************************************/
|
|
|
|
/* END LEGAL */
|
2015-04-18 20:44:19 +01:00
|
|
|
#ifndef GRID_LATTICE_TRANSFER_H
|
|
|
|
#define GRID_LATTICE_TRANSFER_H
|
|
|
|
|
|
|
|
namespace Grid {
|
|
|
|
|
2015-04-24 18:40:44 +01:00
|
|
|
inline void subdivides(GridBase *coarse,GridBase *fine)
|
|
|
|
{
|
|
|
|
assert(coarse->_ndimension == fine->_ndimension);
|
|
|
|
|
|
|
|
int _ndimension = coarse->_ndimension;
|
|
|
|
|
|
|
|
// local and global volumes subdivide cleanly after SIMDization
|
|
|
|
for(int d=0;d<_ndimension;d++){
|
|
|
|
assert(coarse->_processors[d] == fine->_processors[d]);
|
|
|
|
assert(coarse->_simd_layout[d] == fine->_simd_layout[d]);
|
|
|
|
assert((fine->_rdimensions[d] / coarse->_rdimensions[d])* coarse->_rdimensions[d]==fine->_rdimensions[d]);
|
|
|
|
}
|
|
|
|
}
|
|
|
|
|
2016-05-11 14:12:02 +01:00
|
|
|
|
2015-04-18 20:44:19 +01:00
|
|
|
////////////////////////////////////////////////////////////////////////////////////////////
|
|
|
|
// remove and insert a half checkerboard
|
|
|
|
////////////////////////////////////////////////////////////////////////////////////////////
|
|
|
|
template<class vobj> inline void pickCheckerboard(int cb,Lattice<vobj> &half,const Lattice<vobj> &full){
|
|
|
|
half.checkerboard = cb;
|
|
|
|
int ssh=0;
|
2017-02-21 10:24:27 +00:00
|
|
|
//parallel_for
|
2015-04-18 20:44:19 +01:00
|
|
|
for(int ss=0;ss<full._grid->oSites();ss++){
|
|
|
|
std::vector<int> coor;
|
|
|
|
int cbos;
|
|
|
|
|
|
|
|
full._grid->oCoorFromOindex(coor,ss);
|
|
|
|
cbos=half._grid->CheckerBoard(coor);
|
|
|
|
|
|
|
|
if (cbos==cb) {
|
|
|
|
half._odata[ssh] = full._odata[ss];
|
|
|
|
ssh++;
|
|
|
|
}
|
|
|
|
}
|
|
|
|
}
|
|
|
|
template<class vobj> inline void setCheckerboard(Lattice<vobj> &full,const Lattice<vobj> &half){
|
|
|
|
int cb = half.checkerboard;
|
|
|
|
int ssh=0;
|
2017-02-21 10:24:27 +00:00
|
|
|
//parallel_for
|
2015-04-18 20:44:19 +01:00
|
|
|
for(int ss=0;ss<full._grid->oSites();ss++){
|
|
|
|
std::vector<int> coor;
|
|
|
|
int cbos;
|
2015-05-31 15:09:02 +01:00
|
|
|
|
2015-04-18 20:44:19 +01:00
|
|
|
full._grid->oCoorFromOindex(coor,ss);
|
|
|
|
cbos=half._grid->CheckerBoard(coor);
|
|
|
|
|
|
|
|
if (cbos==cb) {
|
|
|
|
full._odata[ss]=half._odata[ssh];
|
|
|
|
ssh++;
|
|
|
|
}
|
|
|
|
}
|
|
|
|
}
|
|
|
|
|
2015-04-24 18:40:44 +01:00
|
|
|
|
2015-06-08 12:04:59 +01:00
|
|
|
template<class vobj,class CComplex,int nbasis>
|
|
|
|
inline void blockProject(Lattice<iVector<CComplex,nbasis > > &coarseData,
|
|
|
|
const Lattice<vobj> &fineData,
|
|
|
|
const std::vector<Lattice<vobj> > &Basis)
|
2015-04-24 18:40:44 +01:00
|
|
|
{
|
|
|
|
GridBase * fine = fineData._grid;
|
|
|
|
GridBase * coarse= coarseData._grid;
|
|
|
|
int _ndimension = coarse->_ndimension;
|
|
|
|
|
|
|
|
// checks
|
|
|
|
assert( nbasis == Basis.size() );
|
|
|
|
subdivides(coarse,fine);
|
|
|
|
for(int i=0;i<nbasis;i++){
|
2015-06-08 12:04:59 +01:00
|
|
|
conformable(Basis[i],fineData);
|
2015-04-24 18:40:44 +01:00
|
|
|
}
|
|
|
|
|
|
|
|
std::vector<int> block_r (_ndimension);
|
|
|
|
|
|
|
|
for(int d=0 ; d<_ndimension;d++){
|
|
|
|
block_r[d] = fine->_rdimensions[d] / coarse->_rdimensions[d];
|
2015-06-08 12:04:59 +01:00
|
|
|
assert(block_r[d]*coarse->_rdimensions[d] == fine->_rdimensions[d]);
|
2015-04-24 18:40:44 +01:00
|
|
|
}
|
|
|
|
|
|
|
|
coarseData=zero;
|
|
|
|
|
|
|
|
// Loop with a cache friendly loop ordering
|
|
|
|
for(int sf=0;sf<fine->oSites();sf++){
|
|
|
|
|
|
|
|
int sc;
|
|
|
|
std::vector<int> coor_c(_ndimension);
|
|
|
|
std::vector<int> coor_f(_ndimension);
|
2016-02-11 13:37:39 +00:00
|
|
|
Lexicographic::CoorFromIndex(coor_f,sf,fine->_rdimensions);
|
2015-04-24 18:40:44 +01:00
|
|
|
for(int d=0;d<_ndimension;d++) coor_c[d]=coor_f[d]/block_r[d];
|
2016-02-11 13:37:39 +00:00
|
|
|
Lexicographic::IndexFromCoor(coor_c,sc,coarse->_rdimensions);
|
2015-04-24 18:40:44 +01:00
|
|
|
|
|
|
|
for(int i=0;i<nbasis;i++) {
|
|
|
|
|
2015-06-08 12:04:59 +01:00
|
|
|
coarseData._odata[sc](i)=coarseData._odata[sc](i)
|
2015-06-09 22:39:13 +01:00
|
|
|
+ innerProduct(Basis[i]._odata[sf],fineData._odata[sf]);
|
2015-04-24 18:40:44 +01:00
|
|
|
|
|
|
|
}
|
|
|
|
}
|
|
|
|
return;
|
|
|
|
}
|
|
|
|
|
2015-06-08 12:04:59 +01:00
|
|
|
template<class vobj,class CComplex>
|
|
|
|
inline void blockZAXPY(Lattice<vobj> &fineZ,
|
|
|
|
const Lattice<CComplex> &coarseA,
|
|
|
|
const Lattice<vobj> &fineX,
|
|
|
|
const Lattice<vobj> &fineY)
|
2015-04-24 18:40:44 +01:00
|
|
|
{
|
2015-06-08 12:04:59 +01:00
|
|
|
GridBase * fine = fineZ._grid;
|
|
|
|
GridBase * coarse= coarseA._grid;
|
2015-04-24 18:40:44 +01:00
|
|
|
|
2015-06-08 12:04:59 +01:00
|
|
|
fineZ.checkerboard=fineX.checkerboard;
|
|
|
|
subdivides(coarse,fine); // require they map
|
|
|
|
conformable(fineX,fineY);
|
|
|
|
conformable(fineX,fineZ);
|
2015-04-24 18:40:44 +01:00
|
|
|
|
2015-06-08 12:04:59 +01:00
|
|
|
int _ndimension = coarse->_ndimension;
|
2015-04-24 18:40:44 +01:00
|
|
|
|
2015-06-08 12:04:59 +01:00
|
|
|
std::vector<int> block_r (_ndimension);
|
|
|
|
|
|
|
|
// FIXME merge with subdivide checking routine as this is redundant
|
2015-04-24 18:40:44 +01:00
|
|
|
for(int d=0 ; d<_ndimension;d++){
|
|
|
|
block_r[d] = fine->_rdimensions[d] / coarse->_rdimensions[d];
|
2015-06-08 12:04:59 +01:00
|
|
|
assert(block_r[d]*coarse->_rdimensions[d]==fine->_rdimensions[d]);
|
2015-04-24 18:40:44 +01:00
|
|
|
}
|
|
|
|
|
2017-02-21 10:24:27 +00:00
|
|
|
parallel_for(int sf=0;sf<fine->oSites();sf++){
|
2015-06-08 12:04:59 +01:00
|
|
|
|
2015-04-24 18:40:44 +01:00
|
|
|
int sc;
|
|
|
|
std::vector<int> coor_c(_ndimension);
|
|
|
|
std::vector<int> coor_f(_ndimension);
|
|
|
|
|
2016-02-11 13:37:39 +00:00
|
|
|
Lexicographic::CoorFromIndex(coor_f,sf,fine->_rdimensions);
|
2015-04-24 18:40:44 +01:00
|
|
|
for(int d=0;d<_ndimension;d++) coor_c[d]=coor_f[d]/block_r[d];
|
2016-02-11 13:37:39 +00:00
|
|
|
Lexicographic::IndexFromCoor(coor_c,sc,coarse->_rdimensions);
|
2015-04-24 18:40:44 +01:00
|
|
|
|
2015-06-08 12:04:59 +01:00
|
|
|
// z = A x + y
|
|
|
|
fineZ._odata[sf]=coarseA._odata[sc]*fineX._odata[sf]+fineY._odata[sf];
|
2015-04-24 18:40:44 +01:00
|
|
|
|
|
|
|
}
|
2015-06-08 12:04:59 +01:00
|
|
|
|
2015-04-24 18:40:44 +01:00
|
|
|
return;
|
|
|
|
}
|
2015-06-08 12:04:59 +01:00
|
|
|
template<class vobj,class CComplex>
|
|
|
|
inline void blockInnerProduct(Lattice<CComplex> &CoarseInner,
|
|
|
|
const Lattice<vobj> &fineX,
|
|
|
|
const Lattice<vobj> &fineY)
|
|
|
|
{
|
|
|
|
typedef decltype(innerProduct(fineX._odata[0],fineY._odata[0])) dotp;
|
2015-04-24 18:40:44 +01:00
|
|
|
|
2015-06-08 12:04:59 +01:00
|
|
|
GridBase *coarse(CoarseInner._grid);
|
|
|
|
GridBase *fine (fineX._grid);
|
|
|
|
|
|
|
|
Lattice<dotp> fine_inner(fine);
|
|
|
|
Lattice<dotp> coarse_inner(coarse);
|
|
|
|
|
|
|
|
fine_inner = localInnerProduct(fineX,fineY);
|
|
|
|
blockSum(coarse_inner,fine_inner);
|
2017-02-21 10:24:27 +00:00
|
|
|
parallel_for(int ss=0;ss<coarse->oSites();ss++){
|
2015-06-08 12:04:59 +01:00
|
|
|
CoarseInner._odata[ss] = coarse_inner._odata[ss];
|
|
|
|
}
|
|
|
|
}
|
|
|
|
template<class vobj,class CComplex>
|
|
|
|
inline void blockNormalise(Lattice<CComplex> &ip,Lattice<vobj> &fineX)
|
|
|
|
{
|
|
|
|
GridBase *coarse = ip._grid;
|
2015-06-09 10:26:19 +01:00
|
|
|
Lattice<vobj> zz(fineX._grid); zz=zero;
|
2015-06-08 12:04:59 +01:00
|
|
|
blockInnerProduct(ip,fineX,fineX);
|
2015-06-20 22:22:56 +01:00
|
|
|
ip = pow(ip,-0.5);
|
2015-06-09 10:26:19 +01:00
|
|
|
blockZAXPY(fineX,ip,fineX,zz);
|
2015-06-08 12:04:59 +01:00
|
|
|
}
|
2015-04-24 18:40:44 +01:00
|
|
|
// useful in multigrid project;
|
|
|
|
// Generic name : Coarsen?
|
|
|
|
template<class vobj>
|
2015-06-08 12:04:59 +01:00
|
|
|
inline void blockSum(Lattice<vobj> &coarseData,const Lattice<vobj> &fineData)
|
2015-04-24 18:40:44 +01:00
|
|
|
{
|
|
|
|
GridBase * fine = fineData._grid;
|
|
|
|
GridBase * coarse= coarseData._grid;
|
|
|
|
|
|
|
|
subdivides(coarse,fine); // require they map
|
|
|
|
|
|
|
|
int _ndimension = coarse->_ndimension;
|
|
|
|
|
|
|
|
std::vector<int> block_r (_ndimension);
|
|
|
|
|
|
|
|
for(int d=0 ; d<_ndimension;d++){
|
|
|
|
block_r[d] = fine->_rdimensions[d] / coarse->_rdimensions[d];
|
|
|
|
}
|
|
|
|
|
|
|
|
coarseData=zero;
|
|
|
|
for(int sf=0;sf<fine->oSites();sf++){
|
|
|
|
|
|
|
|
int sc;
|
|
|
|
std::vector<int> coor_c(_ndimension);
|
|
|
|
std::vector<int> coor_f(_ndimension);
|
|
|
|
|
2016-02-11 13:37:39 +00:00
|
|
|
Lexicographic::CoorFromIndex(coor_f,sf,fine->_rdimensions);
|
2015-04-24 18:40:44 +01:00
|
|
|
for(int d=0;d<_ndimension;d++) coor_c[d]=coor_f[d]/block_r[d];
|
2016-02-11 13:37:39 +00:00
|
|
|
Lexicographic::IndexFromCoor(coor_c,sc,coarse->_rdimensions);
|
2015-04-24 18:40:44 +01:00
|
|
|
|
|
|
|
coarseData._odata[sc]=coarseData._odata[sc]+fineData._odata[sf];
|
|
|
|
|
|
|
|
}
|
|
|
|
return;
|
|
|
|
}
|
|
|
|
|
2015-06-09 10:26:19 +01:00
|
|
|
template<class vobj>
|
|
|
|
inline void blockPick(GridBase *coarse,const Lattice<vobj> &unpicked,Lattice<vobj> &picked,std::vector<int> coor)
|
|
|
|
{
|
|
|
|
GridBase * fine = unpicked._grid;
|
|
|
|
|
|
|
|
Lattice<vobj> zz(fine);
|
|
|
|
Lattice<iScalar<vInteger> > fcoor(fine);
|
|
|
|
|
|
|
|
zz = zero;
|
|
|
|
|
|
|
|
picked = unpicked;
|
|
|
|
for(int d=0;d<fine->_ndimension;d++){
|
|
|
|
LatticeCoordinate(fcoor,d);
|
|
|
|
int block= fine->_rdimensions[d] / coarse->_rdimensions[d];
|
|
|
|
int lo = (coor[d])*block;
|
|
|
|
int hi = (coor[d]+1)*block;
|
|
|
|
picked = where( (fcoor<hi) , picked, zz);
|
|
|
|
picked = where( (fcoor>=lo), picked, zz);
|
|
|
|
}
|
|
|
|
}
|
2015-06-08 12:04:59 +01:00
|
|
|
|
|
|
|
template<class vobj,class CComplex>
|
|
|
|
inline void blockOrthogonalise(Lattice<CComplex> &ip,std::vector<Lattice<vobj> > &Basis)
|
|
|
|
{
|
|
|
|
GridBase *coarse = ip._grid;
|
|
|
|
GridBase *fine = Basis[0]._grid;
|
|
|
|
|
|
|
|
int nbasis = Basis.size() ;
|
|
|
|
int _ndimension = coarse->_ndimension;
|
|
|
|
|
|
|
|
// checks
|
|
|
|
subdivides(coarse,fine);
|
|
|
|
for(int i=0;i<nbasis;i++){
|
|
|
|
conformable(Basis[i]._grid,fine);
|
|
|
|
}
|
|
|
|
|
|
|
|
for(int v=0;v<nbasis;v++) {
|
|
|
|
for(int u=0;u<v;u++) {
|
|
|
|
//Inner product & remove component
|
|
|
|
blockInnerProduct(ip,Basis[u],Basis[v]);
|
|
|
|
ip = -ip;
|
|
|
|
blockZAXPY<vobj,CComplex> (Basis[v],ip,Basis[u],Basis[v]);
|
|
|
|
}
|
|
|
|
blockNormalise(ip,Basis[v]);
|
|
|
|
}
|
|
|
|
}
|
|
|
|
|
|
|
|
template<class vobj,class CComplex,int nbasis>
|
|
|
|
inline void blockPromote(const Lattice<iVector<CComplex,nbasis > > &coarseData,
|
|
|
|
Lattice<vobj> &fineData,
|
|
|
|
const std::vector<Lattice<vobj> > &Basis)
|
|
|
|
{
|
|
|
|
GridBase * fine = fineData._grid;
|
|
|
|
GridBase * coarse= coarseData._grid;
|
|
|
|
int _ndimension = coarse->_ndimension;
|
|
|
|
|
|
|
|
// checks
|
|
|
|
assert( nbasis == Basis.size() );
|
|
|
|
subdivides(coarse,fine);
|
|
|
|
for(int i=0;i<nbasis;i++){
|
2015-06-09 22:39:13 +01:00
|
|
|
conformable(Basis[i]._grid,fine);
|
2015-06-08 12:04:59 +01:00
|
|
|
}
|
|
|
|
|
|
|
|
std::vector<int> block_r (_ndimension);
|
|
|
|
|
|
|
|
for(int d=0 ; d<_ndimension;d++){
|
|
|
|
block_r[d] = fine->_rdimensions[d] / coarse->_rdimensions[d];
|
|
|
|
}
|
|
|
|
|
|
|
|
// Loop with a cache friendly loop ordering
|
|
|
|
for(int sf=0;sf<fine->oSites();sf++){
|
|
|
|
|
|
|
|
int sc;
|
|
|
|
std::vector<int> coor_c(_ndimension);
|
|
|
|
std::vector<int> coor_f(_ndimension);
|
|
|
|
|
2016-02-11 13:37:39 +00:00
|
|
|
Lexicographic::CoorFromIndex(coor_f,sf,fine->_rdimensions);
|
2015-06-08 12:04:59 +01:00
|
|
|
for(int d=0;d<_ndimension;d++) coor_c[d]=coor_f[d]/block_r[d];
|
2016-02-11 13:37:39 +00:00
|
|
|
Lexicographic::IndexFromCoor(coor_c,sc,coarse->_rdimensions);
|
2015-06-08 12:04:59 +01:00
|
|
|
|
|
|
|
for(int i=0;i<nbasis;i++) {
|
2015-06-09 22:39:13 +01:00
|
|
|
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];
|
2015-06-08 12:04:59 +01:00
|
|
|
|
|
|
|
}
|
|
|
|
}
|
|
|
|
return;
|
|
|
|
|
|
|
|
}
|
|
|
|
|
2016-05-11 15:18:47 +01:00
|
|
|
// Useful for precision conversion, or indeed anything where an operator= does a conversion on scalars.
|
|
|
|
// Simd layouts need not match since we use peek/poke Local
|
|
|
|
template<class vobj,class vvobj>
|
|
|
|
void localConvert(const Lattice<vobj> &in,Lattice<vvobj> &out)
|
|
|
|
{
|
|
|
|
typedef typename vobj::scalar_object sobj;
|
|
|
|
typedef typename vvobj::scalar_object ssobj;
|
|
|
|
|
|
|
|
GridBase *ig = in._grid;
|
|
|
|
GridBase *og = out._grid;
|
|
|
|
|
|
|
|
int ni = ig->_ndimension;
|
|
|
|
int no = og->_ndimension;
|
|
|
|
|
|
|
|
assert(ni == no);
|
|
|
|
|
|
|
|
for(int d=0;d<no;d++){
|
|
|
|
assert(ig->_processors[d] == og->_processors[d]);
|
|
|
|
assert(ig->_ldimensions[d] == og->_ldimensions[d]);
|
2017-02-07 06:16:39 +00:00
|
|
|
assert(ig->lSites() == og->lSites());
|
2016-05-11 15:18:47 +01:00
|
|
|
}
|
|
|
|
|
2017-02-21 10:24:27 +00:00
|
|
|
parallel_for(int idx=0;idx<ig->lSites();idx++){
|
2017-02-07 06:16:39 +00:00
|
|
|
sobj s;
|
|
|
|
ssobj ss;
|
|
|
|
|
2016-05-11 15:18:47 +01:00
|
|
|
std::vector<int> lcoor(ni);
|
|
|
|
ig->LocalIndexToLocalCoor(idx,lcoor);
|
|
|
|
peekLocalSite(s,in,lcoor);
|
|
|
|
ss=s;
|
|
|
|
pokeLocalSite(ss,out,lcoor);
|
|
|
|
}
|
|
|
|
}
|
|
|
|
|
|
|
|
|
2016-05-11 14:12:02 +01:00
|
|
|
template<class vobj>
|
2017-04-09 15:41:04 +01:00
|
|
|
void InsertSlice(const Lattice<vobj> &lowDim,Lattice<vobj> & higherDim,int slice, int orthog)
|
2016-05-11 14:12:02 +01:00
|
|
|
{
|
|
|
|
typedef typename vobj::scalar_object sobj;
|
|
|
|
|
2016-05-11 15:06:54 +01:00
|
|
|
GridBase *lg = lowDim._grid;
|
|
|
|
GridBase *hg = higherDim._grid;
|
2016-05-11 14:12:02 +01:00
|
|
|
int nl = lg->_ndimension;
|
|
|
|
int nh = hg->_ndimension;
|
|
|
|
|
|
|
|
assert(nl+1 == nh);
|
|
|
|
assert(orthog<nh);
|
|
|
|
assert(orthog>=0);
|
2016-05-12 18:35:38 +01:00
|
|
|
assert(hg->_processors[orthog]==1);
|
2016-05-11 14:12:02 +01:00
|
|
|
|
|
|
|
int dl; dl = 0;
|
|
|
|
for(int d=0;d<nh;d++){
|
|
|
|
if ( d != orthog) {
|
|
|
|
assert(lg->_processors[dl] == hg->_processors[d]);
|
|
|
|
assert(lg->_ldimensions[dl] == hg->_ldimensions[d]);
|
|
|
|
dl++;
|
|
|
|
}
|
|
|
|
}
|
|
|
|
|
|
|
|
// the above should guarantee that the operations are local
|
2017-02-21 10:24:27 +00:00
|
|
|
parallel_for(int idx=0;idx<lg->lSites();idx++){
|
2017-02-07 06:16:39 +00:00
|
|
|
sobj s;
|
2016-05-11 15:18:47 +01:00
|
|
|
std::vector<int> lcoor(nl);
|
|
|
|
std::vector<int> hcoor(nh);
|
2016-05-11 14:12:02 +01:00
|
|
|
lg->LocalIndexToLocalCoor(idx,lcoor);
|
2017-02-22 17:19:09 +00:00
|
|
|
int ddl=0;
|
2016-05-11 15:18:47 +01:00
|
|
|
hcoor[orthog] = slice;
|
2016-05-11 14:12:02 +01:00
|
|
|
for(int d=0;d<nh;d++){
|
2016-05-11 15:06:54 +01:00
|
|
|
if ( d!=orthog ) {
|
2017-02-22 17:19:09 +00:00
|
|
|
hcoor[d]=lcoor[ddl++];
|
2016-05-11 14:12:02 +01:00
|
|
|
}
|
|
|
|
}
|
2016-05-12 18:35:38 +01:00
|
|
|
peekLocalSite(s,lowDim,lcoor);
|
|
|
|
pokeLocalSite(s,higherDim,hcoor);
|
2016-05-11 14:12:02 +01:00
|
|
|
}
|
|
|
|
}
|
|
|
|
|
|
|
|
template<class vobj>
|
2017-04-09 15:41:04 +01:00
|
|
|
void ExtractSlice(Lattice<vobj> &lowDim,const Lattice<vobj> & higherDim,int slice, int orthog)
|
2016-05-11 14:12:02 +01:00
|
|
|
{
|
|
|
|
typedef typename vobj::scalar_object sobj;
|
|
|
|
|
2016-05-11 15:06:54 +01:00
|
|
|
GridBase *lg = lowDim._grid;
|
|
|
|
GridBase *hg = higherDim._grid;
|
2016-05-11 14:12:02 +01:00
|
|
|
int nl = lg->_ndimension;
|
|
|
|
int nh = hg->_ndimension;
|
|
|
|
|
|
|
|
assert(nl+1 == nh);
|
|
|
|
assert(orthog<nh);
|
|
|
|
assert(orthog>=0);
|
2016-05-12 18:35:38 +01:00
|
|
|
assert(hg->_processors[orthog]==1);
|
2016-05-11 14:12:02 +01:00
|
|
|
|
|
|
|
int dl; dl = 0;
|
2016-07-13 15:49:18 +01:00
|
|
|
for(int d=0;d<nh;d++){
|
|
|
|
if ( d != orthog) {
|
|
|
|
assert(lg->_processors[dl] == hg->_processors[d]);
|
|
|
|
assert(lg->_ldimensions[dl] == hg->_ldimensions[d]);
|
|
|
|
dl++;
|
2016-05-11 14:12:02 +01:00
|
|
|
}
|
|
|
|
}
|
|
|
|
// the above should guarantee that the operations are local
|
2017-02-21 10:24:27 +00:00
|
|
|
parallel_for(int idx=0;idx<lg->lSites();idx++){
|
2017-02-07 06:16:39 +00:00
|
|
|
sobj s;
|
2016-05-11 15:18:47 +01:00
|
|
|
std::vector<int> lcoor(nl);
|
|
|
|
std::vector<int> hcoor(nh);
|
2016-05-11 14:12:02 +01:00
|
|
|
lg->LocalIndexToLocalCoor(idx,lcoor);
|
2017-02-23 00:25:29 +00:00
|
|
|
int ddl=0;
|
2016-05-11 15:18:47 +01:00
|
|
|
hcoor[orthog] = slice;
|
2016-05-11 14:12:02 +01:00
|
|
|
for(int d=0;d<nh;d++){
|
2016-05-11 15:06:54 +01:00
|
|
|
if ( d!=orthog ) {
|
2017-02-22 17:19:09 +00:00
|
|
|
hcoor[d]=lcoor[ddl++];
|
2016-05-11 14:12:02 +01:00
|
|
|
}
|
|
|
|
}
|
2016-05-12 18:35:38 +01:00
|
|
|
peekLocalSite(s,higherDim,hcoor);
|
|
|
|
pokeLocalSite(s,lowDim,lcoor);
|
2016-05-11 14:12:02 +01:00
|
|
|
}
|
|
|
|
|
|
|
|
}
|
2015-06-08 12:04:59 +01:00
|
|
|
|
2016-08-17 01:33:55 +01:00
|
|
|
|
|
|
|
template<class vobj>
|
2017-04-09 15:41:04 +01:00
|
|
|
void InsertSliceLocal(const Lattice<vobj> &lowDim, Lattice<vobj> & higherDim,int slice_lo,int slice_hi, int orthog)
|
2016-08-17 01:33:55 +01:00
|
|
|
{
|
|
|
|
typedef typename vobj::scalar_object sobj;
|
|
|
|
|
|
|
|
GridBase *lg = lowDim._grid;
|
|
|
|
GridBase *hg = higherDim._grid;
|
|
|
|
int nl = lg->_ndimension;
|
|
|
|
int nh = hg->_ndimension;
|
|
|
|
|
|
|
|
assert(nl == nh);
|
|
|
|
assert(orthog<nh);
|
|
|
|
assert(orthog>=0);
|
|
|
|
|
|
|
|
for(int d=0;d<nh;d++){
|
|
|
|
assert(lg->_processors[d] == hg->_processors[d]);
|
|
|
|
assert(lg->_ldimensions[d] == hg->_ldimensions[d]);
|
|
|
|
}
|
|
|
|
|
|
|
|
// the above should guarantee that the operations are local
|
2017-02-21 10:24:27 +00:00
|
|
|
parallel_for(int idx=0;idx<lg->lSites();idx++){
|
2017-02-07 06:16:39 +00:00
|
|
|
sobj s;
|
2016-08-17 01:33:55 +01:00
|
|
|
std::vector<int> lcoor(nl);
|
|
|
|
std::vector<int> hcoor(nh);
|
|
|
|
lg->LocalIndexToLocalCoor(idx,lcoor);
|
|
|
|
if( lcoor[orthog] == slice_lo ) {
|
|
|
|
hcoor=lcoor;
|
|
|
|
hcoor[orthog] = slice_hi;
|
|
|
|
peekLocalSite(s,lowDim,lcoor);
|
|
|
|
pokeLocalSite(s,higherDim,hcoor);
|
|
|
|
}
|
|
|
|
}
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
|
|
template<class vobj>
|
|
|
|
void ExtractSliceLocal(Lattice<vobj> &lowDim, Lattice<vobj> & higherDim,int slice_lo,int slice_hi, int orthog)
|
|
|
|
{
|
|
|
|
typedef typename vobj::scalar_object sobj;
|
|
|
|
|
|
|
|
GridBase *lg = lowDim._grid;
|
|
|
|
GridBase *hg = higherDim._grid;
|
|
|
|
int nl = lg->_ndimension;
|
|
|
|
int nh = hg->_ndimension;
|
|
|
|
|
|
|
|
assert(nl == nh);
|
|
|
|
assert(orthog<nh);
|
|
|
|
assert(orthog>=0);
|
|
|
|
|
|
|
|
for(int d=0;d<nh;d++){
|
|
|
|
assert(lg->_processors[d] == hg->_processors[d]);
|
|
|
|
assert(lg->_ldimensions[d] == hg->_ldimensions[d]);
|
|
|
|
}
|
|
|
|
|
|
|
|
// the above should guarantee that the operations are local
|
2017-02-21 10:24:27 +00:00
|
|
|
parallel_for(int idx=0;idx<lg->lSites();idx++){
|
2017-02-07 06:16:39 +00:00
|
|
|
sobj s;
|
2016-08-17 01:33:55 +01:00
|
|
|
std::vector<int> lcoor(nl);
|
|
|
|
std::vector<int> hcoor(nh);
|
|
|
|
lg->LocalIndexToLocalCoor(idx,lcoor);
|
|
|
|
if( lcoor[orthog] == slice_lo ) {
|
|
|
|
hcoor=lcoor;
|
|
|
|
hcoor[orthog] = slice_hi;
|
|
|
|
peekLocalSite(s,higherDim,hcoor);
|
|
|
|
pokeLocalSite(s,lowDim,lcoor);
|
|
|
|
}
|
|
|
|
}
|
|
|
|
}
|
|
|
|
|
|
|
|
|
2015-08-19 10:26:07 +01:00
|
|
|
template<class vobj>
|
|
|
|
void Replicate(Lattice<vobj> &coarse,Lattice<vobj> & fine)
|
|
|
|
{
|
|
|
|
typedef typename vobj::scalar_object sobj;
|
|
|
|
|
|
|
|
GridBase *cg = coarse._grid;
|
|
|
|
GridBase *fg = fine._grid;
|
|
|
|
|
|
|
|
int nd = cg->_ndimension;
|
|
|
|
|
|
|
|
subdivides(cg,fg);
|
|
|
|
|
|
|
|
assert(cg->_ndimension==fg->_ndimension);
|
|
|
|
|
|
|
|
std::vector<int> ratio(cg->_ndimension);
|
|
|
|
|
|
|
|
for(int d=0;d<cg->_ndimension;d++){
|
|
|
|
ratio[d] = fg->_fdimensions[d]/cg->_fdimensions[d];
|
|
|
|
}
|
|
|
|
|
|
|
|
std::vector<int> fcoor(nd);
|
|
|
|
std::vector<int> ccoor(nd);
|
|
|
|
for(int g=0;g<fg->gSites();g++){
|
|
|
|
|
|
|
|
fg->GlobalIndexToGlobalCoor(g,fcoor);
|
|
|
|
for(int d=0;d<nd;d++){
|
|
|
|
ccoor[d] = fcoor[d]%cg->_gdimensions[d];
|
|
|
|
}
|
|
|
|
|
|
|
|
sobj tmp;
|
|
|
|
peekSite(tmp,coarse,ccoor);
|
|
|
|
pokeSite(tmp,fine,fcoor);
|
|
|
|
}
|
|
|
|
|
|
|
|
}
|
|
|
|
|
2016-07-06 20:57:04 +01:00
|
|
|
//Copy SIMD-vectorized lattice to array of scalar objects in lexicographic order
|
|
|
|
template<typename vobj, typename sobj>
|
|
|
|
typename std::enable_if<isSIMDvectorized<vobj>::value && !isSIMDvectorized<sobj>::value, void>::type unvectorizeToLexOrdArray(std::vector<sobj> &out, const Lattice<vobj> &in){
|
|
|
|
typedef typename vobj::vector_type vtype;
|
|
|
|
|
|
|
|
GridBase* in_grid = in._grid;
|
|
|
|
out.resize(in_grid->lSites());
|
|
|
|
|
|
|
|
int ndim = in_grid->Nd();
|
|
|
|
int in_nsimd = vtype::Nsimd();
|
|
|
|
|
2016-07-06 23:15:15 +01:00
|
|
|
std::vector<std::vector<int> > in_icoor(in_nsimd);
|
2016-07-06 20:57:04 +01:00
|
|
|
|
|
|
|
for(int lane=0; lane < in_nsimd; lane++){
|
|
|
|
in_icoor[lane].resize(ndim);
|
|
|
|
in_grid->iCoorFromIindex(in_icoor[lane], lane);
|
|
|
|
}
|
|
|
|
|
2017-02-21 10:24:27 +00:00
|
|
|
parallel_for(int in_oidx = 0; in_oidx < in_grid->oSites(); in_oidx++){ //loop over outer index
|
2016-07-06 20:57:04 +01:00
|
|
|
//Assemble vector of pointers to output elements
|
|
|
|
std::vector<sobj*> out_ptrs(in_nsimd);
|
|
|
|
|
|
|
|
std::vector<int> in_ocoor(ndim);
|
|
|
|
in_grid->oCoorFromOindex(in_ocoor, in_oidx);
|
|
|
|
|
|
|
|
std::vector<int> lcoor(in_grid->Nd());
|
|
|
|
|
|
|
|
for(int lane=0; lane < in_nsimd; lane++){
|
|
|
|
for(int mu=0;mu<ndim;mu++)
|
|
|
|
lcoor[mu] = in_ocoor[mu] + in_grid->_rdimensions[mu]*in_icoor[lane][mu];
|
|
|
|
|
|
|
|
int lex;
|
|
|
|
Lexicographic::IndexFromCoor(lcoor, lex, in_grid->_ldimensions);
|
|
|
|
out_ptrs[lane] = &out[lex];
|
|
|
|
}
|
|
|
|
|
|
|
|
//Unpack into those ptrs
|
|
|
|
const vobj & in_vobj = in._odata[in_oidx];
|
|
|
|
extract1(in_vobj, out_ptrs, 0);
|
|
|
|
}
|
|
|
|
}
|
|
|
|
|
|
|
|
//Convert a Lattice from one precision to another
|
|
|
|
template<class VobjOut, class VobjIn>
|
|
|
|
void precisionChange(Lattice<VobjOut> &out, const Lattice<VobjIn> &in){
|
|
|
|
assert(out._grid->Nd() == in._grid->Nd());
|
|
|
|
out.checkerboard = in.checkerboard;
|
|
|
|
GridBase *in_grid=in._grid;
|
|
|
|
GridBase *out_grid = out._grid;
|
2015-08-19 10:26:07 +01:00
|
|
|
|
2016-07-06 20:57:04 +01:00
|
|
|
typedef typename VobjOut::scalar_object SobjOut;
|
|
|
|
typedef typename VobjIn::scalar_object SobjIn;
|
|
|
|
|
|
|
|
int ndim = out._grid->Nd();
|
|
|
|
int out_nsimd = out_grid->Nsimd();
|
|
|
|
|
2016-07-06 23:22:15 +01:00
|
|
|
std::vector<std::vector<int> > out_icoor(out_nsimd);
|
2016-07-06 20:57:04 +01:00
|
|
|
|
|
|
|
for(int lane=0; lane < out_nsimd; lane++){
|
|
|
|
out_icoor[lane].resize(ndim);
|
|
|
|
out_grid->iCoorFromIindex(out_icoor[lane], lane);
|
|
|
|
}
|
|
|
|
|
|
|
|
std::vector<SobjOut> in_slex_conv(in_grid->lSites());
|
|
|
|
unvectorizeToLexOrdArray(in_slex_conv, in);
|
|
|
|
|
2017-02-21 10:24:27 +00:00
|
|
|
parallel_for(int out_oidx=0;out_oidx<out_grid->oSites();out_oidx++){
|
2016-07-06 20:57:04 +01:00
|
|
|
std::vector<int> out_ocoor(ndim);
|
|
|
|
out_grid->oCoorFromOindex(out_ocoor, out_oidx);
|
|
|
|
|
|
|
|
std::vector<SobjOut*> ptrs(out_nsimd);
|
|
|
|
|
|
|
|
std::vector<int> lcoor(out_grid->Nd());
|
|
|
|
|
|
|
|
for(int lane=0; lane < out_nsimd; lane++){
|
|
|
|
for(int mu=0;mu<ndim;mu++)
|
|
|
|
lcoor[mu] = out_ocoor[mu] + out_grid->_rdimensions[mu]*out_icoor[lane][mu];
|
|
|
|
|
|
|
|
int llex; Lexicographic::IndexFromCoor(lcoor, llex, out_grid->_ldimensions);
|
|
|
|
ptrs[lane] = &in_slex_conv[llex];
|
|
|
|
}
|
|
|
|
merge(out._odata[out_oidx], ptrs, 0);
|
|
|
|
}
|
|
|
|
}
|
|
|
|
|
2015-04-18 20:44:19 +01:00
|
|
|
}
|
|
|
|
#endif
|