Below is what I did with the ordering. The issue is that the array of 4d spinors is
a full parity spinor, with the even parity part in the first half of the 4d spinor and
the odd parity part in the second half. If the 5d quda spinor is a full parity spinor
then we have to reorder the parity parts so all Ls evens go into the first half of
the 5d spinor and all Ls odds go in the second half of the 5d array. That gives one
type of mapping for the variable x below. On the other hand, if the 5d quda spinor
is a parity spinor, then we have to know which parity so that we know which
part of the 4d spinor to use. Also the mapping for x is different. The result
of this is that the function needs to know whether it is a full parity spinor or
a half, and if half, it needs to know which parity. This would require additions
to the cpuColorSpinorField so that the information gets stored in field.
Joel
QOPDomainWallOrder(cpuColorSpinorField &field) : ColorSpinorFieldOrder(field),
field(field), volume_4d_h(1), Ls(0)
{
if (field.Ndim() != 5) errorQuda("Error, wrong number of dimensions for this ColorSpinorFieldOrder");
for (int i=0; i<4; i++) volume_4dh = field.x[i];
if (subset == QUDA_FULL_SITE_SUBSET) volume_4d_h /= 2;
Ls = field.x[4];
volume_5d_h=volume_4d_h_Ls;
}
virtual ~QOPDomainWallOrder() { ; }
const Float& operator()(const int &x, const int &s, const int &c, const int &z) const {
int ls;
int x_4d;
if (subset == QUDA_FULL_SITESUBSET) {
// Here x=x1+L1h(x2+L2(x3+L3(x4+L4_(xs+Ls_oddBit)))).
int oddBit=x/volume_5d_h;
int x_5d=x-oddBit_volume_5d_h;
ls=x_5d/volume_4d_h;
x_4d=x_5d-ls_volume_4d_h;
if (oddBit) x_4d += volume_4dh;
} else {
// Here x=x1+L1h(x2+L2(x3+L3(x4+L4_xs))).
ls = x / volume_4d_h;
x_4d = x - ls_volume_4d_h;
if (parity == QUDA_ODD_PARITY) {
x_4d += volume_4d_h;
}
}
unsigned long index_4d = ((x_4d_field.nColor+c)_field.nSpin+s)2+z;
return ((Float*)(field.v))[ls][index_4d];
}
(Bug report and suggested fix from Joel Giedt)
Below is what I did with the ordering. The issue is that the array of 4d spinors is a full parity spinor, with the even parity part in the first half of the 4d spinor and the odd parity part in the second half. If the 5d quda spinor is a full parity spinor then we have to reorder the parity parts so all Ls evens go into the first half of the 5d spinor and all Ls odds go in the second half of the 5d array. That gives one type of mapping for the variable x below. On the other hand, if the 5d quda spinor is a parity spinor, then we have to know which parity so that we know which part of the 4d spinor to use. Also the mapping for x is different. The result of this is that the function needs to know whether it is a full parity spinor or a half, and if half, it needs to know which parity. This would require additions to the cpuColorSpinorField so that the information gets stored in field.
Joel
QOPDomainWallOrder(cpuColorSpinorField &field) : ColorSpinorFieldOrder(field),
field(field), volume_4d_h(1), Ls(0)
{
if (field.Ndim() != 5) errorQuda("Error, wrong number of dimensions for this ColorSpinorFieldOrder");
for (int i=0; i<4; i++) volume_4dh = field.x[i];
if (subset == QUDA_FULL_SITE_SUBSET) volume_4d_h /= 2;
Ls = field.x[4];
volume_5d_h=volume_4d_h_Ls;
}
virtual ~QOPDomainWallOrder() { ; }
const Float& operator()(const int &x, const int &s, const int &c, const int &z) const { int ls; int x_4d; if (subset == QUDA_FULL_SITESUBSET) { // Here x=x1+L1h(x2+L2(x3+L3(x4+L4_(xs+Ls_oddBit)))). int oddBit=x/volume_5d_h; int x_5d=x-oddBit_volume_5d_h; ls=x_5d/volume_4d_h; x_4d=x_5d-ls_volume_4d_h; if (oddBit) x_4d += volume_4dh; } else { // Here x=x1+L1h(x2+L2(x3+L3(x4+L4_xs))). ls = x / volume_4d_h; x_4d = x - ls_volume_4d_h; if (parity == QUDA_ODD_PARITY) { x_4d += volume_4d_h; } } unsigned long index_4d = ((x_4d_field.nColor+c)_field.nSpin+s)2+z; return ((Float*)(field.v))[ls][index_4d]; }