diff --git a/cpp/field.cpp b/cpp/field.cpp index c33e6296787fbd9db143778d9293d7e1ed114f68..f523daf2a95e23855fc41c21dab34ad68374dd8f 100644 --- a/cpp/field.cpp +++ b/cpp/field.cpp @@ -1351,11 +1351,11 @@ void field<rnumber, be, fc>::Hermitian_reflect() 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); + ptrdiff_t cindex = this->get_cindex(0, iyy, iz); for (int cc = 0; cc < int(ncomp(fc)); cc++) { - (*(cdata + ncomp(fc)*cindex + cc))[0] = (*(bufferm + ncomp(fc)*iz+cc))[0]; - (*(cdata + ncomp(fc)*cindex + cc))[1] = -(*(bufferm + ncomp(fc)*iz+cc))[1]; + (*(cdata + ncomp(fc)*cindex + cc))[0] = (*(bufferm + ncomp(fc)*izz+cc))[0]; + (*(cdata + ncomp(fc)*cindex + cc))[1] = -(*(bufferm + ncomp(fc)*izz+cc))[1]; } } ptrdiff_t cindex = this->get_cindex(0, iyy, 0); @@ -1391,7 +1391,7 @@ void field<rnumber, be, fc>::Hermitian_reflect() delete[] bufferp; delete[] bufferm; delete mpistatus; - // symmetrize kx = 0, ky = 0 line + // reflect 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++) @@ -1569,106 +1569,111 @@ void field<rnumber, be, fc>::symmetrize() // 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) + // if my rank is plus or minus, then I should do actual work + if (this->clayout->myrank == rankp || + this->clayout->myrank == rankm) { - 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 my rank is rank plus, I should fill out the plus buffer 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); + 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, receive buffer plus, send buffer minus + // if my rank is rank minus, I should fill out the minus buffer 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, iz); - for (int cc = 0; cc < int(ncomp(fc)); cc++) + 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++) { - (*(cdata + ncomp(fc)*cindex + cc))[0] = ((*(bufferp + ncomp(fc)*iz+cc))[0] + (*(bufferm + ncomp(fc)*izz+cc))[0])/2; - (*(cdata + ncomp(fc)*cindex + cc))[1] = ((*(bufferp + ncomp(fc)*iz+cc))[1] - (*(bufferm + ncomp(fc)*izz+cc))[1])/2; + 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]; } } - ptrdiff_t cindex = this->get_cindex(0, iyy, 0); - for (int cc = 0; cc < int(ncomp(fc)); cc++) + // if ranks are different, send and receive + if (rankp != rankm) { - (*(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 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 == 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++) + // if I my rank is either plus or minus, I should update my slice + if (this->clayout->myrank == rankp) { - ptrdiff_t izz = (this->clayout->sizes[1] - iz); - ptrdiff_t cindex = this->get_cindex(0, iyy, izz); + 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, iz); + for (int cc = 0; cc < int(ncomp(fc)); cc++) + { + (*(cdata + ncomp(fc)*cindex + cc))[0] = ((*(bufferp + ncomp(fc)*iz+cc))[0] + (*(bufferm + ncomp(fc)*izz+cc))[0])/2; + (*(cdata + ncomp(fc)*cindex + cc))[1] = ((*(bufferp + ncomp(fc)*iz+cc))[1] - (*(bufferm + ncomp(fc)*izz+cc))[1])/2; + } + } + ptrdiff_t cindex = this->get_cindex(0, iyy, 0); for (int cc = 0; cc < int(ncomp(fc)); cc++) { - (*(cdata + ncomp(fc)*cindex + cc))[0] = ((*(bufferp + ncomp(fc)*iz+cc))[0] + (*(bufferm + ncomp(fc)*izz+cc))[0])/2; - (*(cdata + ncomp(fc)*cindex + cc))[1] = -((*(bufferp + ncomp(fc)*iz+cc))[1] - (*(bufferm + ncomp(fc)*izz+cc))[1])/2; + (*(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; } } - ptrdiff_t cindex = this->get_cindex(0, iyy, 0); - for (int cc = 0; cc < int(ncomp(fc)); cc++) + // if I my rank is either plus or minus, I should update my slice + if (this->clayout->myrank == rankm) { - (*(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; + 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)*izz+cc))[0])/2; + (*(cdata + ncomp(fc)*cindex + cc))[1] = -((*(bufferp + ncomp(fc)*iz+cc))[1] - (*(bufferm + ncomp(fc)*izz+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; + } } } }