portelli/Grid
1
0
Fork 0
mirror of https://github.com/paboyle/Grid.git synced 2025-06-10 19:36:56 +01:00
Code Releases Activity

Large scale change to support 5d fermion formulations.

Browse Source 
Have 5d replicated wilson with 4d gauge working and matrix regressing
to Ls copies of wilson.
This commit is contained in:
Peter Boyle
2015-05-31 15:09:02 +01:00
parent 59db857ad1
commit 5644ab1e19
45 changed files with 1549 additions and 692 deletions
Download Patch File Download Diff File Expand all files Collapse all files

103
lib/stencil/Grid_lebesgue.cc Normal file
View File

@ -0,0 +1,103 @@
#include <Grid.h>
namespace Grid {
int LebesgueOrder::UseLebesgueOrder;
LebesgueOrder::IndexInteger LebesgueOrder::alignup(IndexInteger n){
n--; // 1000 0011 --> 1000 0010
n |= n >> 1; // 1000 0010 | 0100 0001 = 1100 0011
n |= n >> 2; // 1100 0011 | 0011 0000 = 1111 0011
n |= n >> 4; // 1111 0011 | 0000 1111 = 1111 1111
n |= n >> 8; // ... (At this point all bits are 1, so further bitwise-or
n |= n >> 16; // operations produce no effect.)
n++; // 1111 1111 --> 1 0000 0000
return n;
};
LebesgueOrder::LebesgueOrder(GridBase *grid)
{
_LebesgueReorder.resize(0);
// Align up dimensions to power of two.
const IndexInteger one=1;
IndexInteger ND = grid->_ndimension;
std::vector<IndexInteger> dims(ND);
std::vector<IndexInteger> adims(ND);
std::vector<std::vector<IndexInteger> > bitlist(ND);
for(IndexInteger mu=0;mu<ND;mu++){
dims[mu] = grid->_rdimensions[mu];
assert ( dims[mu] != 0 );
adims[mu] = alignup(dims[mu]);
}
// List which bits of padded volume coordinate contribute; this strategy
// i) avoids recursion
// ii) has loop lengths at most the width of a 32 bit word.
int sitebit=0;
int split=2;
for(int mu=0;mu<ND;mu++){ // mu 0 takes bit 0; mu 1 takes bit 1 etc...
for(int bit=0;bit<split;bit++){
IndexInteger mask = one<<bit;
if ( mask&(adims[mu]-1) ){
bitlist[mu].push_back(sitebit);
sitebit++;
}
}
}
for(int bit=split;bit<32;bit++){
IndexInteger mask = one<<bit;
for(int mu=0;mu<ND;mu++){ // mu 0 takes bit 0; mu 1 takes bit 1 etc...
if ( mask&(adims[mu]-1) ){
bitlist[mu].push_back(sitebit);
sitebit++;
}
}
}
// Work out padded and unpadded volumes
IndexInteger avol = 1;
for(int mu=0;mu<ND;mu++) avol = avol * adims[mu];
IndexInteger vol = 1;
for(int mu=0;mu<ND;mu++) vol = vol * dims[mu];
// Loop over padded volume, following Lebesgue curve
// We interleave the bits from sequential "mu".
std::vector<IndexInteger> ax(ND);
for(IndexInteger asite=0;asite<avol;asite++){
// Start with zero and collect bits
for(int mu=0;mu<ND;mu++) ax[mu] = 0;
int contained = 1;
for(int mu=0;mu<ND;mu++){
// Build the coordinate on the aligned volume
for(int bit=0;bit<bitlist[mu].size();bit++){
int sbit=bitlist[mu][bit];
if(asite&(one<<sbit)){
ax[mu]|=one<<bit;
}
}
// Is it contained in original box
if ( ax[mu]>dims[mu]-1 ) contained = 0;
}
if ( contained ) {
int site = ax[0]
+ dims[0]*ax[1]
+dims[0]*dims[1]*ax[2]
+dims[0]*dims[1]*dims[2]*ax[3];
_LebesgueReorder.push_back(site);
}
}
assert( _LebesgueReorder.size() == vol );
}
}

29
lib/stencil/Grid_lebesgue.h Normal file
View File

@ -0,0 +1,29 @@
#ifndef GRID_LEBESGUE_H
#define GRID_LEBESGUE_H
#include<vector>
// Lebesgue, Morton, Z-graph ordering assistance
namespace Grid {
class LebesgueOrder {
public:
static int UseLebesgueOrder;
typedef uint32_t IndexInteger;
inline IndexInteger Reorder(IndexInteger ss) {
return UseLebesgueOrder ? _LebesgueReorder[ss] : ss;
};
IndexInteger alignup(IndexInteger n);
LebesgueOrder(GridBase *grid);
private:
std::vector<IndexInteger> _LebesgueReorder;
};
}
#endif

103
lib/stencil/Grid_stencil_common.cc
View File

@ -3,95 +3,6 @@
namespace Grid {
void CartesianStencil::LebesgueOrder(void)
{
_LebesgueReorder.resize(0);
// Align up dimensions to power of two.
const StencilInteger one=1;
StencilInteger ND = _grid->_ndimension;
std::vector<StencilInteger> dims(ND);
std::vector<StencilInteger> adims(ND);
std::vector<std::vector<StencilInteger> > bitlist(ND);
for(StencilInteger mu=0;mu<ND;mu++){
dims[mu] = _grid->_rdimensions[mu];
assert ( dims[mu] != 0 );
adims[mu] = alignup(dims[mu]);
}
// List which bits of padded volume coordinate contribute; this strategy
// i) avoids recursion
// ii) has loop lengths at most the width of a 32 bit word.
int sitebit=0;
int split=24;
for(int mu=0;mu<ND;mu++){ // mu 0 takes bit 0; mu 1 takes bit 1 etc...
for(int bit=0;bit<split;bit++){
StencilInteger mask = one<<bit;
if ( mask&(adims[mu]-1) ){
bitlist[mu].push_back(sitebit);
sitebit++;
}
}
}
for(int bit=split;bit<32;bit++){
StencilInteger mask = one<<bit;
for(int mu=0;mu<ND;mu++){ // mu 0 takes bit 0; mu 1 takes bit 1 etc...
if ( mask&(adims[mu]-1) ){
bitlist[mu].push_back(sitebit);
sitebit++;
}
}
}
// Work out padded and unpadded volumes
StencilInteger avol = 1;
for(int mu=0;mu<ND;mu++) avol = avol * adims[mu];
StencilInteger vol = 1;
for(int mu=0;mu<ND;mu++) vol = vol * dims[mu];
// Loop over padded volume, following Lebesgue curve
// We interleave the bits from sequential "mu".
std::vector<StencilInteger> ax(ND);
for(StencilInteger asite=0;asite<avol;asite++){
// Start with zero and collect bits
for(int mu=0;mu<ND;mu++) ax[mu] = 0;
int contained = 1;
for(int mu=0;mu<ND;mu++){
// Build the coordinate on the aligned volume
for(int bit=0;bit<bitlist[mu].size();bit++){
int sbit=bitlist[mu][bit];
if(asite&(one<<sbit)){
ax[mu]|=one<<bit;
}
}
// Is it contained in original box
if ( ax[mu]>dims[mu]-1 ) contained = 0;
}
if ( contained ) {
int site = ax[0]
+ dims[0]*ax[1]
+dims[0]*dims[1]*ax[2]
+dims[0]*dims[1]*dims[2]*ax[3];
_LebesgueReorder.push_back(site);
}
}
assert( _LebesgueReorder.size() == vol );
}
CartesianStencil::CartesianStencil(GridBase *grid,
int npoints,
int checkerboard,
@ -110,8 +21,6 @@ void CartesianStencil::LebesgueOrder(void)
_unified_buffer_size=0;
_request_count =0;
LebesgueOrder();
int osites = _grid->oSites();
for(int i=0;i<npoints;i++){
@ -143,8 +52,8 @@ void CartesianStencil::LebesgueOrder(void)
// up a table containing the npoint "neighbours" and whether they
// live in lattice or a comms buffer.
if ( !comm_dim ) {
sshift[0] = _grid->CheckerBoardShift(_checkerboard,dimension,shift,0);
sshift[1] = _grid->CheckerBoardShift(_checkerboard,dimension,shift,1);
sshift[0] = _grid->CheckerBoardShift(_checkerboard,dimension,shift,Even);
sshift[1] = _grid->CheckerBoardShift(_checkerboard,dimension,shift,Odd);
if ( sshift[0] == sshift[1] ) {
Local(point,dimension,shift,0x3);
@ -154,8 +63,8 @@ void CartesianStencil::LebesgueOrder(void)
}
} else { // All permute extract done in comms phase prior to Stencil application
// So tables are the same whether comm_dim or splice_dim
sshift[0] = _grid->CheckerBoardShift(_checkerboard,dimension,shift,0);
sshift[1] = _grid->CheckerBoardShift(_checkerboard,dimension,shift,1);
sshift[0] = _grid->CheckerBoardShift(_checkerboard,dimension,shift,Even);
sshift[1] = _grid->CheckerBoardShift(_checkerboard,dimension,shift,Odd);
if ( sshift[0] == sshift[1] ) {
Comms(point,dimension,shift,0x3);
} else {
@ -185,7 +94,7 @@ void CartesianStencil::LebesgueOrder(void)
int o = 0;
int bo = x * _grid->_ostride[dimension];
int cb= (cbmask==0x2)? 1 : 0;
int cb= (cbmask==0x2)? Odd : Even;
int sshift = _grid->CheckerBoardShift(_checkerboard,dimension,shift,cb);
int sx = (x+sshift)%rd;
@ -224,7 +133,7 @@ void CartesianStencil::LebesgueOrder(void)
_comm_buf_size[point] = buffer_size; // Size of _one_ plane. Multiple planes may be gathered and
// send to one or more remote nodes.
int cb= (cbmask==0x2)? 1 : 0;
int cb= (cbmask==0x2)? Odd : Even;
int sshift= _grid->CheckerBoardShift(_checkerboard,dimension,shift,cb);
for(int x=0;x<rd;x++){