diff --git a/cpp/field.cpp b/cpp/field.cpp
index a5779752ee4ca232f4883042c3171ec4c45ab18a..c4cb192a8a543733d015d43c793104c7b5969868 100644
--- a/cpp/field.cpp
+++ b/cpp/field.cpp
@@ -1266,23 +1266,17 @@ void field<rnumber, be, fc>::normalize()
template <typename rnumber,
field_backend be,
field_components fc>
-void field<rnumber, be, fc>::symmetrize()
+void field<rnumber, be, fc>::Hermitian_reflect()
{
- TIMEZONE("field::symmetrize");
+ TIMEZONE("field::Hermitian_reflect");
assert(!this->real_space_representation);
- // for debugging, just use FFTW
- //this->ift();
- //this->dft();
- //this->normalize();
- //return;
typename fftw_interface<rnumber>::complex *cdata = this->get_cdata();
- // symmetrize kx = 0 plane, line by line, for ky != 0
+ // reflect kx = 0 plane, line by line, for ky != 0
MPI_Status *mpistatus = new MPI_Status;
// bufferp will hold data from "plus", i.e. iy
// bufferm will hold data from "minus", i.e. ny - iy
- typename fftw_interface<rnumber>::complex *bufferp = new typename fftw_interface<rnumber>::complex[ncomp(fc)*this->clayout->sizes[1]];
+ typename fftw_interface<rnumber>::complex *bufferp = new typename fftw_interface<rnumber>::complex[ncomp(fc)*this->clayout->sizes[1]];
typename fftw_interface<rnumber>::complex *bufferm = new typename fftw_interface<rnumber>::complex[ncomp(fc)*this->clayout->sizes[1]];
- /*ptrdiff_t tindex;*/
int rankp, rankm;
// for each ky slice except ky=0
for (ptrdiff_t iy = 1; iy < ptrdiff_t(this->clayout->sizes[0]/2); iy++)
@@ -1360,15 +1354,15 @@ void field<rnumber, be, fc>::symmetrize()
ptrdiff_t cindex = this->get_cindex(0, iyy, izz);
for (int cc = 0; cc < int(ncomp(fc)); cc++)
{
- (*(cdata + ncomp(fc)*cindex + cc))[0] = ((*(bufferp + ncomp(fc)*iz+cc))[0] + (*(bufferm + ncomp(fc)*iz+cc))[0])/2;
- (*(cdata + ncomp(fc)*cindex + cc))[1] = ((*(bufferp + ncomp(fc)*iz+cc))[1] - (*(bufferp + ncomp(fc)*iz+cc))[1])/2;
+ (*(cdata + ncomp(fc)*cindex + cc))[0] = (*(bufferm + ncomp(fc)*iz+cc))[0];
+ (*(cdata + ncomp(fc)*cindex + cc))[1] = -(*(bufferm + ncomp(fc)*iz+cc))[1];
}
}
ptrdiff_t cindex = this->get_cindex(0, iyy, 0);
for (int cc = 0; cc < int(ncomp(fc)); cc++)
{
- (*(cdata + cc + ncomp(fc)*cindex))[0] = ((*(bufferp + cc))[0] + (*(bufferm + cc))[0])/2;
- (*(cdata + cc + ncomp(fc)*cindex))[1] = ((*(bufferp + cc))[1] - (*(bufferm + cc))[1])/2;
+ (*(cdata + cc + ncomp(fc)*cindex))[0] = (*(bufferm + cc))[0];
+ (*(cdata + cc + ncomp(fc)*cindex))[1] = -(*(bufferm + cc))[1];
}
}
// if I my rank is either plus or minus, I should update my slice
@@ -1381,15 +1375,15 @@ void field<rnumber, be, fc>::symmetrize()
ptrdiff_t cindex = this->get_cindex(0, iyy, izz);
for (int cc = 0; cc < int(ncomp(fc)); cc++)
{
- (*(cdata + ncomp(fc)*cindex + cc))[0] = ((*(bufferp + ncomp(fc)*iz+cc))[0] + (*(bufferm + ncomp(fc)*iz+cc))[0])/2;
- (*(cdata + ncomp(fc)*cindex + cc))[1] = -((*(bufferp + ncomp(fc)*iz+cc))[1] - (*(bufferp + ncomp(fc)*iz+cc))[1])/2;
+ (*(cdata + ncomp(fc)*cindex + cc))[0] = (*(bufferp + ncomp(fc)*iz+cc))[0];
+ (*(cdata + ncomp(fc)*cindex + cc))[1] = -(*(bufferp + ncomp(fc)*iz+cc))[1];
}
}
ptrdiff_t cindex = this->get_cindex(0, iyy, 0);
for (int cc = 0; cc < int(ncomp(fc)); cc++)
{
- (*(cdata + cc + ncomp(fc)*cindex))[0] = ((*(bufferp + cc))[0] + (*(bufferm + cc))[0])/2;
- (*(cdata + cc + ncomp(fc)*cindex))[1] = -((*(bufferp + cc))[1] - (*(bufferm + cc))[1])/2;
+ (*(cdata + cc + ncomp(fc)*cindex))[0] = (*(bufferp + cc))[0];
+ (*(cdata + cc + ncomp(fc)*cindex))[1] = -(*(bufferp + cc))[1];
}
}
}
@@ -1406,12 +1400,14 @@ void field<rnumber, be, fc>::symmetrize()
ptrdiff_t cindex1 = this->get_cindex(0, 0, this->clayout->sizes[1] - iz);
for (int cc = 0; cc < int(ncomp(fc)); cc++)
{
- double re = ((*(cdata + cc + ncomp(fc)*cindex0))[0] + (*(cdata + cc + ncomp(fc)*cindex1))[0])/2;
- double im = ((*(cdata + cc + ncomp(fc)*cindex0))[1] - (*(cdata + cc + ncomp(fc)*cindex1))[1])/2;
- (*(cdata + cc + ncomp(fc)*cindex0))[0] = re;
- (*(cdata + cc + ncomp(fc)*cindex0))[1] = im;
- (*(cdata + cc + ncomp(fc)*cindex1))[0] = re;
- (*(cdata + cc + ncomp(fc)*cindex1))[1] = -im;
+ double rep = (*(cdata + cc + ncomp(fc)*cindex0))[0];
+ double imp = (*(cdata + cc + ncomp(fc)*cindex0))[1];
+ double rem = (*(cdata + cc + ncomp(fc)*cindex1))[0];
+ double imm = (*(cdata + cc + ncomp(fc)*cindex1))[1];
+ (*(cdata + cc + ncomp(fc)*cindex0))[0] = rem;
+ (*(cdata + cc + ncomp(fc)*cindex0))[1] = -imm;
+ (*(cdata + cc + ncomp(fc)*cindex1))[0] = rep;
+ (*(cdata + cc + ncomp(fc)*cindex1))[1] = -imp;
}
}
}
@@ -1456,6 +1452,200 @@ void field<rnumber, be, fc>::symmetrize()
}
}
+template <typename rnumber,
+ field_backend be,
+ field_components fc>
+void field<rnumber, be, fc>::symmetrize()
+{
+ TIMEZONE("field::symmetrize");
+ assert(!this->real_space_representation);
+ typename fftw_interface<rnumber>::complex *cdata = this->get_cdata();
+
+ // make 0 mode real
+ if (this->myrank == this->clayout->rank[0][0])
+ {
+ for (ptrdiff_t cc = 0; cc < ncomp(fc); cc++)
+ cdata[cc][1] = 0.0;
+ }
+ // put kx = nx/2 modes to 0
+ for (ptrdiff_t iy = 0; iy < ptrdiff_t(this->clayout->subsizes[0]); iy++)
+ for (ptrdiff_t iz = 0; iz < ptrdiff_t(this->clayout->subsizes[1]); iz++)
+ {
+ ptrdiff_t cindex = this->get_cindex(this->clayout->sizes[2]-1, iy, iz);
+ for (int cc = 0; cc < int(ncomp(fc)); cc++) {
+ (*(cdata + cc + ncomp(fc)*cindex))[0] = 0.0;
+ (*(cdata + cc + ncomp(fc)*cindex))[1] = 0.0;
+ }
+ }
+ // put ky = ny/2 modes to 0
+ if (this->clayout->myrank == this->clayout->rank[0][this->clayout->sizes[0]/2])
+ {
+ for (ptrdiff_t iz = 0; iz < ptrdiff_t(this->clayout->subsizes[1]); iz++)
+ for (ptrdiff_t ix = 0; ix < ptrdiff_t(this->clayout->subsizes[2]); ix++)
+ {
+ ptrdiff_t cindex = this->get_cindex(ix, this->clayout->sizes[0]/2-this->clayout->starts[0], iz);
+ for (int cc = 0; cc < int(ncomp(fc)); cc++) {
+ (*(cdata + cc + ncomp(fc)*cindex))[0] = 0.0;
+ (*(cdata + cc + ncomp(fc)*cindex))[1] = 0.0;
+ }
+ }
+ }
+ // put kz = nz/2 modes to 0
+ for (ptrdiff_t iy = 0; iy < ptrdiff_t(this->clayout->subsizes[0]); iy++)
+ for (ptrdiff_t ix = 0; ix < ptrdiff_t(this->clayout->subsizes[2]); ix++)
+ {
+ ptrdiff_t cindex = this->get_cindex(ix, iy, this->clayout->sizes[1]/2);
+ for (int cc = 0; cc < int(ncomp(fc)); cc++) {
+ (*(cdata + cc + ncomp(fc)*cindex))[0] = 0.0;
+ (*(cdata + cc + ncomp(fc)*cindex))[1] = 0.0;
+ }
+ }
+
+ // for debugging, just use FFTW
+ this->ift();
+ this->dft();
+ this->normalize();
+ return;
+ // symmetrize kx = 0 plane, line by line, for ky != 0
+ MPI_Status *mpistatus = new MPI_Status;
+ // bufferp will hold data from "plus", i.e. iy
+ // bufferm will hold data from "minus", i.e. ny - iy
+ typename fftw_interface<rnumber>::complex *bufferp = new typename fftw_interface<rnumber>::complex[ncomp(fc)*this->clayout->sizes[1]];
+ typename fftw_interface<rnumber>::complex *bufferm = new typename fftw_interface<rnumber>::complex[ncomp(fc)*this->clayout->sizes[1]];
+ int rankp, rankm;
+ // for each ky slice except ky=0
+ for (ptrdiff_t iy = 1; iy < ptrdiff_t(this->clayout->sizes[0]/2); iy++)
+ {
+ // read rank plus and rank minus
+ rankp = this->clayout->rank[0][iy];
+ rankm = this->clayout->rank[0][this->clayout->sizes[0] - iy];
+ // if my rank is rank plus, I should fill out the plus buffer
+ if (this->clayout->myrank == rankp)
+ {
+ ptrdiff_t iyy = iy - this->clayout->starts[0];
+ for (ptrdiff_t iz = 0; iz < ptrdiff_t(this->clayout->sizes[1]); iz++)
+ {
+ ptrdiff_t cindex = this->get_cindex(0, iyy, iz);
+ for (int cc = 0; cc < int(ncomp(fc)); cc++)
+ for (int imag_comp=0; imag_comp<2; imag_comp++)
+ (*(bufferp + ncomp(fc)*iz+cc))[imag_comp] =
+ (*(cdata + ncomp(fc)*cindex + cc))[imag_comp];
+ }
+ }
+ // if my rank is rank minus, I should fill out the minus buffer
+ if (this->clayout->myrank == rankm)
+ {
+ ptrdiff_t iyy = (this->clayout->sizes[0] - iy) - this->clayout->starts[0];
+ for (ptrdiff_t iz = 0; iz < ptrdiff_t(this->clayout->sizes[1]); iz++)
+ {
+ ptrdiff_t cindex = this->get_cindex(0, iyy, iz);
+ for (int cc = 0; cc < int(ncomp(fc)); cc++)
+ for (int imag_comp=0; imag_comp<2; imag_comp++)
+ (*(bufferm + ncomp(fc)*iz+cc))[imag_comp] =
+ (*(cdata + ncomp(fc)*cindex + cc))[imag_comp];
+ }
+ }
+ // if ranks are different, send and receive
+ if (rankp != rankm)
+ {
+ // if my rank is rank plus, send buffer plus, receive buffer minus
+ if (this->clayout->myrank == rankp)
+ {
+ MPI_Send((void*)bufferp,
+ ncomp(fc)*this->clayout->sizes[1],
+ mpi_real_type<rnumber>::complex(),
+ rankm, 2*iy+0,
+ this->clayout->comm);
+ MPI_Recv((void*)bufferm,
+ ncomp(fc)*this->clayout->sizes[1],
+ mpi_real_type<rnumber>::complex(),
+ rankm, 2*iy+1,
+ this->clayout->comm,
+ mpistatus);
+ }
+ // if my rank is rank minus, receive buffer plus, send buffer minus
+ if (this->clayout->myrank == rankm)
+ {
+ MPI_Recv((void*)bufferp,
+ ncomp(fc)*this->clayout->sizes[1],
+ mpi_real_type<rnumber>::complex(),
+ rankp, 2*iy+0,
+ this->clayout->comm,
+ mpistatus);
+ MPI_Send((void*)bufferm,
+ ncomp(fc)*this->clayout->sizes[1],
+ mpi_real_type<rnumber>::complex(),
+ rankp, 2*iy+1,
+ this->clayout->comm);
+ }
+ }
+ // if I my rank is either plus or minus, I should update my slice
+ if (this->clayout->myrank == rankp)
+ {
+ ptrdiff_t iyy = iy - this->clayout->starts[0];
+ for (ptrdiff_t iz = 1; iz < ptrdiff_t(this->clayout->sizes[1]); iz++)
+ {
+ ptrdiff_t izz = (this->clayout->sizes[1] - iz);
+ ptrdiff_t cindex = this->get_cindex(0, iyy, izz);
+ for (int cc = 0; cc < int(ncomp(fc)); cc++)
+ {
+ (*(cdata + ncomp(fc)*cindex + cc))[0] = ((*(bufferp + ncomp(fc)*iz+cc))[0] + (*(bufferm + ncomp(fc)*iz+cc))[0])/2;
+ (*(cdata + ncomp(fc)*cindex + cc))[1] = ((*(bufferp + ncomp(fc)*iz+cc))[1] - (*(bufferm + ncomp(fc)*iz+cc))[1])/2;
+ }
+ }
+ ptrdiff_t cindex = this->get_cindex(0, iyy, 0);
+ for (int cc = 0; cc < int(ncomp(fc)); cc++)
+ {
+ (*(cdata + cc + ncomp(fc)*cindex))[0] = ((*(bufferp + cc))[0] + (*(bufferm + cc))[0])/2;
+ (*(cdata + cc + ncomp(fc)*cindex))[1] = ((*(bufferp + cc))[1] - (*(bufferm + cc))[1])/2;
+ }
+ }
+ // if I my rank is either plus or minus, I should update my slice
+ if (this->clayout->myrank == rankm)
+ {
+ ptrdiff_t iyy = (this->clayout->sizes[0] - iy) - this->clayout->starts[0];
+ for (ptrdiff_t iz = 1; iz < ptrdiff_t(this->clayout->sizes[1]); iz++)
+ {
+ ptrdiff_t izz = (this->clayout->sizes[1] - iz);
+ ptrdiff_t cindex = this->get_cindex(0, iyy, izz);
+ for (int cc = 0; cc < int(ncomp(fc)); cc++)
+ {
+ (*(cdata + ncomp(fc)*cindex + cc))[0] = ((*(bufferp + ncomp(fc)*iz+cc))[0] + (*(bufferm + ncomp(fc)*iz+cc))[0])/2;
+ (*(cdata + ncomp(fc)*cindex + cc))[1] = -((*(bufferp + ncomp(fc)*iz+cc))[1] - (*(bufferm + ncomp(fc)*iz+cc))[1])/2;
+ }
+ }
+ ptrdiff_t cindex = this->get_cindex(0, iyy, 0);
+ for (int cc = 0; cc < int(ncomp(fc)); cc++)
+ {
+ (*(cdata + cc + ncomp(fc)*cindex))[0] = ((*(bufferp + cc))[0] + (*(bufferm + cc))[0])/2;
+ (*(cdata + cc + ncomp(fc)*cindex))[1] = -((*(bufferp + cc))[1] - (*(bufferm + cc))[1])/2;
+ }
+ }
+ }
+ //fftw_interface<rnumber>::free(buffer);
+ delete[] bufferp;
+ delete[] bufferm;
+ delete mpistatus;
+ // symmetrize kx = 0, ky = 0 line
+ if (this->clayout->myrank == this->clayout->rank[0][0])
+ {
+ for (ptrdiff_t iz = 1; iz < ptrdiff_t(this->clayout->sizes[1]/2); iz++)
+ {
+ ptrdiff_t cindex0 = this->get_cindex(0, 0, iz);
+ ptrdiff_t cindex1 = this->get_cindex(0, 0, this->clayout->sizes[1] - iz);
+ for (int cc = 0; cc < int(ncomp(fc)); cc++)
+ {
+ double re = ((*(cdata + cc + ncomp(fc)*cindex0))[0] + (*(cdata + cc + ncomp(fc)*cindex1))[0])/2;
+ double im = ((*(cdata + cc + ncomp(fc)*cindex0))[1] - (*(cdata + cc + ncomp(fc)*cindex1))[1])/2;
+ (*(cdata + cc + ncomp(fc)*cindex0))[0] = re;
+ (*(cdata + cc + ncomp(fc)*cindex0))[1] = im;
+ (*(cdata + cc + ncomp(fc)*cindex1))[0] = re;
+ (*(cdata + cc + ncomp(fc)*cindex1))[1] = -im;
+ }
+ }
+ }
+}
+
template <typename rnumber,
field_backend be,
field_components fc>
diff --git a/cpp/field.hpp b/cpp/field.hpp
index feb19c804766470197e37f66f8405d8b3bfb9c9f..cd7a06d74fd4bc670920121976bd28008b40003c 100644
--- a/cpp/field.hpp
+++ b/cpp/field.hpp
@@ -122,6 +122,7 @@ class field
void ift();
void normalize();
void symmetrize();
+ void Hermitian_reflect();
/* stats */
void compute_rspace_xincrement_stats(