From 3ece1425a749b1ff0a1e7e2fec3731656a035d30 Mon Sep 17 00:00:00 2001
From: eric bylaska <eric.bylaska@pnnl.gov>
Date: Mon, 20 May 2024 13:58:28 -0700
Subject: [PATCH] memory leak bug fix...EJB

---
 Nwpw/nwpwlib/C3dB/CBalance.cpp   |   23 +-
 Nwpw/nwpwlib/C3dB/CGrid.cpp      |  796 +++++++++++----------
 Nwpw/nwpwlib/C3dB/c3db.cpp       |   62 +-
 Nwpw/nwpwlib/D3dB/d3db.cpp       |  180 ++---
 Nwpw/nwpwlib/device/gdevice2.cpp |   10 +-
 Nwpw/nwpwlib/lattice/Balance.cpp |   50 +-
 Nwpw/nwpwlib/lattice/PGrid.cpp   | 1127 +++++++++++++++++++-----------
 Nwpw/nwpwlib/utilities/util.hpp  |    6 +
 8 files changed, 1348 insertions(+), 906 deletions(-)

diff --git a/Nwpw/nwpwlib/C3dB/CBalance.cpp b/Nwpw/nwpwlib/C3dB/CBalance.cpp
index 345e69df..1b0a29e3 100644
--- a/Nwpw/nwpwlib/C3dB/CBalance.cpp
+++ b/Nwpw/nwpwlib/C3dB/CBalance.cpp
@@ -295,7 +295,7 @@ CBalance::~CBalance()
 
 /***********************************
  *                                 *
- *    CBalance::c_unbalance_start   *
+ *   CBalance::c_unbalance_start   *
  *                                 *
  ***********************************/
 /**
@@ -311,24 +311,23 @@ CBalance::~CBalance()
  */
 void CBalance::c_unbalance_start(const int nffts, const int nb, double *a, const int request_indx, const int msgtype) 
 {
-   int j, pto, pfrom, msglen, indx;
- 
+   //int j, pto, pfrom, msglen, indx;
    if (sender_list[nb])
-      for (j = 0; j < npacket_list[nb]; ++j) 
+      for (auto j=0; j<npacket_list[nb]; ++j) 
       {
-         pfrom = proc_to_list[nb][j];
-         msglen = 2*packet_size_list[nb][j];
-         indx = 2*indx_start_list[nb][j];
+         int pfrom  = proc_to_list[nb][j];
+         int msglen = 2*packet_size_list[nb][j];
+         int indx   = 2*indx_start_list[nb][j];
          if (msglen > 0)
             parall->adreceive(request_indx, msgtype, pfrom, msglen, a + indx);
       }
  
    if (receiver_list[nb])
-      for (j = 0; j < npacket_list[nb]; ++j) 
+      for (auto j=0; j<npacket_list[nb]; ++j) 
       {
-         pto = proc_from_list[nb][j];
-         msglen = 2*packet_size_list[nb][j];
-         indx = 2*indx_start_list[nb][j];
+         int pto    = proc_from_list[nb][j];
+         int msglen = 2*packet_size_list[nb][j];
+         int indx   = 2*indx_start_list[nb][j];
          if (msglen > 0)
             parall->adsend(request_indx, msgtype, pto, msglen, a + indx);
       }
@@ -357,7 +356,7 @@ void CBalance::c_unbalance_end(const int nffts, const int nb, double *a, const i
 
 /********************************
  *                              *
- *    CBalance::c_unbalance      *
+ *    CBalance::c_unbalance     *
  *                              *
  ********************************/
 /**
diff --git a/Nwpw/nwpwlib/C3dB/CGrid.cpp b/Nwpw/nwpwlib/C3dB/CGrid.cpp
index 3c9362fb..3f8b8779 100644
--- a/Nwpw/nwpwlib/C3dB/CGrid.cpp
+++ b/Nwpw/nwpwlib/C3dB/CGrid.cpp
@@ -1383,7 +1383,8 @@ void CGrid::c_unpack_start(const int nffts, const int nb, double *tmp1, double *
                            const int request_indx, const int msgtype) 
 {
    if (balanced)
-      mybalance->c_unbalance_start(nffts, nb, tmp1, request_indx, msgtype);
+      for (auto s=0; s<nffts; ++s)
+         mybalance->c_unbalance_start(nffts, nb, tmp1+s*n2ft3d, request_indx, msgtype);
 }
 
 /********************************
@@ -1397,10 +1398,11 @@ void CGrid::c_unpack_mid(const int nffts, const int nb, double *tmp1, double *tm
    if (balanced)
      mybalance->c_unbalance_end(nffts, nb, tmp1, request_indx);
  
-   std::memcpy(tmp2, tmp1, 2 * (nidb2[nb]) * sizeof(double));
-   std::memset(tmp1, 0, 2*nfft3d * sizeof(double));
+   std::memcpy(tmp2, tmp1, nffts*2*(nidb2[nb])*sizeof(double));
+   std::memset(tmp1, 0, nffts*2*nfft3d*sizeof(double));
  
-   c_bindexcopy((nidb2[nb]), packarray[nb], tmp2, tmp1);
+   for (auto s=0; s<nffts; ++s)
+      c_bindexcopy((nidb2[nb]), packarray[nb], tmp2+s*2*nfft3d, tmp1+s*2*nfft3d);
 }
 
 /********************************
@@ -1432,29 +1434,32 @@ void CGrid::pfftbz(const int nffts, const int nb, double *tmp1, double *tmp2, in
        ***     A(kx,nz,ky) <- fft1d^(-1)[A(kx,kz,ky)]  ***
        ***************************************************/
       int indx0 = 0;
-      int indx2 = 0;
       int nn = 0;
-      for (auto q = 0; q < nq; ++q) 
+      for (auto s=0; s<nffts; ++s)
       {
-         for (auto i = 0; i<nx; ++i) 
+         int indx2 = 0;
+         for (auto q=0; q<nq; ++q) 
          {
-            if (!zero_row3[nb][indx2]) 
+            for (auto i=0; i<nx; ++i) 
             {
-               auto kk = 0;
-               auto indx3 = 2*i + indx0;
-               auto shift = 2*nz*nn;
-               for (auto k=0; k<nz; ++k) 
+               if (!zero_row3[nb][indx2]) 
                {
-                  tmp2[kk   + shift] = tmp1[indx3];
-                  tmp2[kk+1 + shift] = tmp1[indx3 + 1];
-                  kk += 2;
-                  indx3 += nx;
+                  auto kk = 0;
+                  auto indx3 = 2*i + indx0;
+                  auto shift = 2*nz*nn;
+                  for (auto k=0; k<nz; ++k) 
+                  {
+                     tmp2[kk   + shift] = tmp1[indx3];
+                     tmp2[kk+1 + shift] = tmp1[indx3 + 1];
+                     kk += 2;
+                     indx3 += nx;
+                  }
+                  nn += 1;
                }
-               nn += 1;
+               ++indx2;
             }
-            ++indx2;
+            indx0 += nxz;
          }
-         indx0 += nxz;
       }
      
       //c3db::mygdevice.batch_cfftz_tmpz(c3db::fft_tag,false, nz, nn, 2*nfft3d, tmp2, c3db::tmpz);
@@ -1463,29 +1468,32 @@ void CGrid::pfftbz(const int nffts, const int nb, double *tmp1, double *tmp2, in
       //    dcfftb_(&nz,tmp2+2*nz*i,c3db::tmpz);
      
       indx0 = 0;
-      indx2 = 0;
       nn = 0;
-      for (auto q=0; q<nq; ++q) 
+      for (auto s=0; s<nffts; ++s)
       {
-         for (auto i=0; i<nx; ++i) 
+         int indx2 = 0;
+         for (auto q=0; q<nq; ++q) 
          {
-            if (!zero_row3[nb][indx2]) 
+            for (auto i=0; i<nx; ++i) 
             {
-               auto kk = 0;
-               auto indx3 = 2 * i + indx0;
-               auto shift = 2 * nz * nn;
-               for (auto k=0; k<nz; ++k) 
+               if (!zero_row3[nb][indx2]) 
                {
-                  tmp1[indx3]   = tmp2[kk   + shift];
-                  tmp1[indx3+1] = tmp2[kk+1 + shift];
-                  kk += 2;
-                  indx3 += nx;
+                  auto kk = 0;
+                  auto indx3 = 2 * i + indx0;
+                  auto shift = 2 * nz * nn;
+                  for (auto k=0; k<nz; ++k) 
+                  {
+                     tmp1[indx3]   = tmp2[kk   + shift];
+                     tmp1[indx3+1] = tmp2[kk+1 + shift];
+                     kk += 2;
+                     indx3 += nx;
+                  }
+                  nn += 1;
                }
-               nn += 1;
+               ++indx2;
             }
-            ++indx2;
+            indx0 += nxz;
          }
-         indx0 += nxz;
       }
      
      
@@ -1505,9 +1513,9 @@ void CGrid::pfftbz(const int nffts, const int nb, double *tmp1, double *tmp2, in
        ***   A(nz,kx,ky) <- fft1d^(-1)[A(kz,kx,ky)] ***
        ************************************************/
       //c3db::mygdevice.batch_cfftz_tmpz_zero(c3db::fft_tag,false,nz,nq3,2*nfft3d,tmp1,c3db::tmpz,zero_row3[nb]);
-      c3db::mygdevice.batch_cfft_zero(c3db::fft_tag,false,nz,nq3,nz,tmp1,c3db::backward_z,c3db::tmpz,zero_row3[nb],2);
+      c3db::mygdevice.batch_cfft_zero(c3db::fft_tag,false,nz,nffts*nq3,nz,tmp1,c3db::backward_z,c3db::tmpz,zero_row3[nb],2);
      
-      c3db::c_ptranspose_ijk_start(nffts, nb, 2, tmp1, tmp2, tmp1, request_indx, 45);
+      c3db::c_ptranspose_ijk_start(nffts,nb,2,tmp1,tmp2,tmp1,request_indx, 45);
       // c3db::c_ptranspose_ijk(nb,2,tmp1,tmp2,tmp1);
    }
 }
@@ -1537,29 +1545,32 @@ void CGrid::pfftby(const int nffts, const int nb, double *tmp1, double *tmp2, in
        ***   A(kx,ky,nz) <- fft1d[A(kx,ny,nz)]  ***
        ********************************************/
       int indx0 = 0;
-      int indx2 = 0;
       int nn = 0;
-      for (auto q=0; q<nq; ++q) 
+      for (auto s=0; s<nffts; ++s)
       {
-         for (auto i=0; i<nx; ++i) 
+         int indx2 = 0;
+         for (auto q=0; q<nq; ++q) 
          {
-            if (!zero_row2[nb][indx2]) 
+            for (auto i=0; i<nx; ++i) 
             {
-               auto jj = 0;
-               auto indx3 = 2*i + indx0;
-               auto shift = 2*ny*nn;
-               for (auto j=0; j<ny; ++j) 
+               if (!zero_row2[nb][indx2]) 
                {
-                  tmp1[jj   + shift] = tmp2[indx3];
-                  tmp1[jj+1 + shift] = tmp2[indx3 + 1];
-                  jj += 2;
-                  indx3 += nx;
+                  auto jj = 0;
+                  auto indx3 = 2*i + indx0;
+                  auto shift = 2*ny*nn;
+                  for (auto j=0; j<ny; ++j) 
+                  {
+                     tmp1[jj   + shift] = tmp2[indx3];
+                     tmp1[jj+1 + shift] = tmp2[indx3 + 1];
+                     jj += 2;
+                     indx3 += nx;
+                  }
+                  nn += 1;
                }
-               nn += 1;
+               ++indx2;
             }
-            ++indx2;
+            indx0 += nxy;
          }
-         indx0 += nxy;
       }
      
       //c3db::mygdevice.batch_cffty_tmpy(c3db::fft_tag,false, ny, nn, 2*nfft3d, tmp1, c3db::tmpy);
@@ -1568,31 +1579,36 @@ void CGrid::pfftby(const int nffts, const int nb, double *tmp1, double *tmp2, in
       //    dcfftb_(&ny,tmp1+2*ny*i,c3db::tmpy);
      
       indx0 = 0;
-      indx2 = 0;
       nn = 0;
-      for (auto q=0; q<nq; ++q) 
+      for (auto s=0; s<nffts; ++s)
       {
-         for (auto i=0; i<nx; ++i) 
+         int indx2 = 0;
+         for (auto q=0; q<nq; ++q) 
          {
-            if (!zero_row2[nb][indx2]) 
+            for (auto i=0; i<nx; ++i) 
             {
-               auto jj = 0;
-               auto indx3 = 2*i + indx0;
-               auto shift = 2*ny*nn;
-               for (auto j = 0; j < ny; ++j) 
+               if (!zero_row2[nb][indx2]) 
                {
-                  tmp2[indx3]   = tmp1[jj   + shift];
-                  tmp2[indx3+1] = tmp1[jj+1 + shift];
-                  jj += 2;
-                  indx3 += nx;
+                  auto jj = 0;
+                  auto indx3 = 2*i + indx0;
+                  auto shift = 2*ny*nn;
+                  for (auto j = 0; j < ny; ++j) 
+                  {
+                     tmp2[indx3]   = tmp1[jj   + shift];
+                     tmp2[indx3+1] = tmp1[jj+1 + shift];
+                     jj += 2;
+                     indx3 += nx;
+                  }
+                  nn += 1;
                }
-               nn += 1;
+               ++indx2;
             }
-            ++indx2;
+            indx0 += nxy;
          }
-         indx0 += nxy;
       }
+
    }
+
    /*************************
     **** hilbert mapping ****
     *************************/
@@ -1605,9 +1621,9 @@ void CGrid::pfftby(const int nffts, const int nb, double *tmp1, double *tmp2, in
        ***   A(ky,nz,kx) <- fft1d[A(ny,nz,kx)]  ***
        ********************************************/
       //c3db::mygdevice.batch_cffty_tmpy_zero(c3db::fft_tag,false,ny,nq2,2*nfft3d,tmp2,c3db::tmpy,zero_row2[nb]);
-      c3db::mygdevice.batch_cfft_zero(c3db::fft_tag,false,ny,nq2,ny,tmp2,c3db::backward_y,c3db::tmpy,zero_row2[nb],1);
+      c3db::mygdevice.batch_cfft_zero(c3db::fft_tag,false,ny,nffts*nq2,ny,tmp2,c3db::backward_y,c3db::tmpy,zero_row2[nb],1);
      
-      c3db::c_ptranspose_ijk_start(nffts, nb, 3, tmp2, tmp1, tmp2, request_indx, 46);
+      c3db::c_ptranspose_ijk_start(nffts,nb,3,tmp2,tmp1,tmp2,request_indx, 46);
       // c3db::c_ptranspose_ijk(nb,3,tmp2,tmp1,tmp2);
    }
 }
@@ -1630,24 +1646,25 @@ void CGrid::pfftbx(const int nffts, const int nb, double *tmp1, double *tmp2, in
        ***   A(nx,ny,nz) <- fft1d^(-1)[A(kx,ny,nz)] ***
        ************************************************/
       //c3db::mygdevice.batch_cfftx_tmpx(c3db::fft_tag,false, nx, ny * nq, 2*nfft3d, tmp2, c3db::tmpx);
-      c3db::mygdevice.batch_cfft(c3db::fft_tag,false,nx,ny*nq,nx,tmp2,c3db::backward_x,c3db::tmpx,0);
-      std::memcpy(tmp1, tmp2, 2*nfft3d * sizeof(double));
+      c3db::mygdevice.batch_cfft(c3db::fft_tag,false,nx,nffts*ny*nq,nx,tmp2,c3db::backward_x,c3db::tmpx,0);
+      std::memcpy(tmp1,tmp2,nffts*2*nfft3d*sizeof(double));
    }
    /*************************
     **** hilbert mapping ****
     *************************/
    else 
    {
-      c3db::c_ptranspose_ijk_end(nffts, nb, 3, tmp1, tmp2, request_indx);
+      c3db::c_ptranspose_ijk_end(nffts,nb,3,tmp1,tmp2,request_indx);
      
       /************************************************
        ***     do fft along kx dimension            ***
        ***   A(nx,ny,nz) <- fft1d^(-1)[A(kx,ny,nz)] ***
        ************************************************/
       //c3db::mygdevice.batch_cfftx_tmpx(c3db::fft_tag,false, nx, nq1, 2*nfft3d, tmp1, c3db::tmpx);
-      c3db::mygdevice.batch_cfft(c3db::fft_tag,false,nx,nq1,nx,tmp1,c3db::backward_x,c3db::tmpx,0);
+      c3db::mygdevice.batch_cfft(c3db::fft_tag,false,nx,nffts*nq1,nx,tmp1,c3db::backward_x,c3db::tmpx,0);
       if (2*nfft3d_map < 2*nfft3d)
-         std::memset(tmp1 + 2*nfft3d_map, 0, (2*nfft3d - 2*nfft3d_map) * sizeof(double));
+         for (auto s=0; s<nffts; ++s)
+            std::memset(tmp1+s*2*nfft3d + 2*nfft3d_map, 0, (2*nfft3d-2*nfft3d_map)*sizeof(double));
    }
 }
 
@@ -1663,7 +1680,8 @@ void CGrid::pfftb_step(const int step, const int nffts, const int nb, double *a,
    {
       // c3db::parall->astart(request_indx,parall->np_i());
       // unpack start, tmp1-->tmp1
-      std::memcpy(tmp1, a, 2 * (nidb[nb]) * sizeof(double));
+      for (auto s=0; s<nffts; ++s)
+      std::memcpy(tmp1, a, nffts*2*(nidb[nb])*sizeof(double));
       this->c_unpack_start(nffts, nb, tmp1, tmp2, request_indx, 47);
    } 
    else if (step == 1) 
@@ -1715,36 +1733,39 @@ void CGrid::pfftbz_start(const int nffts, const int nb, double *tmp1, double *tm
       ***     A(kx,nz,ky) <- fft1d^(-1)[A(kx,kz,ky)]  ***
       ***************************************************/
      int indx0 = 0;
-     int indx2 = 0;
      int nn = 0;
-     for (auto q=0; q<nq; ++q) 
+     for (auto s=0; s<nffts; ++s)
      {
-        for (auto i=0; i<nx; ++i) 
+        int indx2 = 0;
+        for (auto q=0; q<nq; ++q) 
         {
-           if (!zero_row3[nb][indx2]) 
+           for (auto i=0; i<nx; ++i) 
            {
-              auto kk = 0;
-              auto indx3 = 2*i + indx0;
-              auto shift = 2*nz*nn;
-              for (auto k=0; k<nz; ++k) 
+              if (!zero_row3[nb][indx2]) 
               {
-                 tmp2[kk   + shift] = tmp1[indx3];
-                 tmp2[kk+1 + shift] = tmp1[indx3+1];
-                 kk += 2;
-                 indx3 += nx;
+                 auto kk = 0;
+                 auto indx3 = 2*i + indx0;
+                 auto shift = 2*nz*nn;
+                 for (auto k=0; k<nz; ++k) 
+                 {
+                    tmp2[kk   + shift] = tmp1[indx3];
+                    tmp2[kk+1 + shift] = tmp1[indx3+1];
+                    kk += 2;
+                    indx3 += nx;
+                 }
+                 nn += 1;
               }
-              nn += 1;
+              ++indx2;
            }
-           ++indx2;
+           indx0 += nxz;
         }
-        indx0 += nxz;
      }
     
      c3db::mygdevice.batch_cfftz_stages_tmpz(0,c3db::fft_tag,false, nz, nn, 2*nfft3d, tmp2, c3db::tmpz,da_indx);
      // for (auto i=0; i<nn; ++i)
      //    dcfftb_(&nz,tmp2+2*nz*i,c3db::tmpz);
-
   }
+
   /*************************
    **** hilbert mapping ****
    *************************/
@@ -1754,7 +1775,7 @@ void CGrid::pfftbz_start(const int nffts, const int nb, double *tmp1, double *tm
       ***     do fft along kz dimension            ***
       ***   A(nz,kx,ky) <- fft1d^(-1)[A(kz,kx,ky)] ***
       ************************************************/
-     c3db::mygdevice.batch_cfftz_stages_tmpz_zero(0,c3db::fft_tag,false, nz, nq3, 2*nfft3d, tmp1, c3db::tmpz, zero_row3[nb], da_indx);
+     c3db::mygdevice.batch_cfftz_stages_tmpz_zero(0,c3db::fft_tag,false, nz, nffts*nq3, 2*nfft3d, tmp1, c3db::tmpz, zero_row3[nb], da_indx);
   }
 
 }
@@ -1778,17 +1799,20 @@ void CGrid::pfftbz_compute(const int nffts, const int nb, double *tmp1, double *
       ***     A(kx,nz,ky) <- fft1d^(-1)[A(kx,kz,ky)]  ***
       ***************************************************/
      int indx0 = 0;
-     int indx2 = 0;
      int nn = 0;
-     for (auto q=0; q<nq; ++q) 
+     for (auto s=0; s<nffts; ++s)
      {
-        for (auto i=0; i<nx; ++i) 
+        int indx2 = 0;
+        for (auto q=0; q<nq; ++q) 
         {
-           if (!zero_row3[nb][indx2]) 
+           for (auto i=0; i<nx; ++i) 
            {
-              nn += 1;
+              if (!zero_row3[nb][indx2]) 
+              {
+                 nn += 1;
+              }
+              ++indx2;
            }
-           ++indx2;
         }
      }
     
@@ -1806,7 +1830,7 @@ void CGrid::pfftbz_compute(const int nffts, const int nb, double *tmp1, double *
        ***     do fft along kz dimension            ***
        ***   A(nz,kx,ky) <- fft1d^(-1)[A(kz,kx,ky)] ***
        ************************************************/
-      c3db::mygdevice.batch_cfftz_stages_tmpz_zero(1,c3db::fft_tag,false, nz, nq3, 2*nfft3d, tmp1, c3db::tmpz,zero_row3[nb],da_indx);
+      c3db::mygdevice.batch_cfftz_stages_tmpz_zero(1,c3db::fft_tag,false, nz, nffts*nq3, 2*nfft3d, tmp1, c3db::tmpz,zero_row3[nb],da_indx);
   }
 
 }
@@ -1830,17 +1854,20 @@ void CGrid::pfftbz_end(const int nffts, const int nb, double *tmp1, double *tmp2
        ***     A(kx,nz,ky) <- fft1d^(-1)[A(kx,kz,ky)]  ***
        ***************************************************/
       int indx0 = 0;
-      int indx2 = 0;
       int nn = 0;
-      for (auto q=0; q<nq; ++q) 
+      for (auto s=0; s<nffts; ++s)
       {
-         for (auto i=0; i<nx; ++i) 
+         int indx2 = 0;
+         for (auto q=0; q<nq; ++q) 
          {
-            if (!zero_row3[nb][indx2]) 
+            for (auto i=0; i<nx; ++i) 
             {
-               nn += 1;
+               if (!zero_row3[nb][indx2]) 
+               {
+                  nn += 1;
+               }
+               ++indx2;
             }
-            ++indx2;
          }
       }
      
@@ -1849,29 +1876,32 @@ void CGrid::pfftbz_end(const int nffts, const int nb, double *tmp1, double *tmp2
       //    dcfftb_(&nz,tmp2+2*nz*i,c3db::tmpz);
      
       indx0 = 0;
-      indx2 = 0;
       nn = 0;
-      for (auto q=0; q<nq; ++q) 
+      for (auto s=0; s<nffts; ++s)
       {
-         for (auto i=0; i<nx; ++i) 
+         int indx2 = 0;
+         for (auto q=0; q<nq; ++q) 
          {
-            if (!zero_row3[nb][indx2]) 
+            for (auto i=0; i<nx; ++i) 
             {
-               auto kk = 0;
-               auto indx3 = 2 * i + indx0;
-               auto shift = 2 * nz * nn;
-               for (auto k=0; k<nz; ++k) 
+               if (!zero_row3[nb][indx2]) 
                {
-                  tmp1[indx3]   = tmp2[kk   + shift];
-                  tmp1[indx3+1] = tmp2[kk+1 + shift];
-                  kk += 2;
-                  indx3 += nx;
+                  auto kk = 0;
+                  auto indx3 = 2 * i + indx0;
+                  auto shift = 2 * nz * nn;
+                  for (auto k=0; k<nz; ++k) 
+                  {
+                     tmp1[indx3]   = tmp2[kk   + shift];
+                     tmp1[indx3+1] = tmp2[kk+1 + shift];
+                     kk += 2;
+                     indx3 += nx;
+                  }
+                  nn += 1;
                }
-               nn += 1;
+               ++indx2;
             }
-            ++indx2;
+            indx0 += nxz;
          }
-         indx0 += nxz;
       }
      
       /***********************************************
@@ -1890,15 +1920,13 @@ void CGrid::pfftbz_end(const int nffts, const int nb, double *tmp1, double *tmp2
       ***     do fft along kz dimension            ***
       ***   A(nz,kx,ky) <- fft1d^(-1)[A(kz,kx,ky)] ***
       ************************************************/
-     c3db::mygdevice.batch_cfftz_stages_tmpz_zero(2,c3db::fft_tag,false, nz, nq3, 2*nfft3d, tmp1, c3db::tmpz, zero_row3[nb], da_indx);
+     c3db::mygdevice.batch_cfftz_stages_tmpz_zero(2,c3db::fft_tag,false,nz,nffts*nq3,2*nfft3d,tmp1,c3db::tmpz,zero_row3[nb],da_indx);
  
-     c3db::c_ptranspose_ijk_start(nffts, nb, 2, tmp1, tmp2, tmp1, request_indx, 45);
+     c3db::c_ptranspose_ijk_start(nffts,nb,2,tmp1,tmp2,tmp1,request_indx, 45);
      // c3db::c_ptranspose_ijk(nb,2,tmp1,tmp2,tmp1);
    }
 }
 
-
-
 /********************************
  *                              *
  *      CGrid::pfftby_start     *
@@ -1924,48 +1952,51 @@ void CGrid::pfftby_start(const int nffts, const int nb, double *tmp1, double *tm
        ***   A(kx,ky,nz) <- fft1d[A(kx,ny,nz)]  ***
        ********************************************/
       int indx0 = 0;
-      int indx2 = 0;
       int nn = 0;
-      for (auto q=0; q<nq; ++q) 
+      for (auto s=0; s<nffts; ++s)
       {
-         for (auto i=0; i<nx; ++i) 
+         int indx2 = 0;
+         for (auto q=0; q<nq; ++q) 
          {
-            if (!zero_row2[nb][indx2]) 
+            for (auto i=0; i<nx; ++i) 
             {
-               auto jj = 0;
-               auto indx3 = 2*i + indx0;
-               auto shift = 2*ny*nn;
-               for (auto j = 0; j < ny; ++j) 
+               if (!zero_row2[nb][indx2]) 
                {
-                  tmp1[jj   + shift] = tmp2[indx3];
-                  tmp1[jj+1 + shift] = tmp2[indx3 + 1];
-                  jj += 2;
-                  indx3 += nx;
+                  auto jj = 0;
+                  auto indx3 = 2*i + indx0;
+                  auto shift = 2*ny*nn;
+                  for (auto j = 0; j < ny; ++j) 
+                  {
+                     tmp1[jj   + shift] = tmp2[indx3];
+                     tmp1[jj+1 + shift] = tmp2[indx3 + 1];
+                     jj += 2;
+                     indx3 += nx;
+                  }
+                  nn += 1;
                }
-               nn += 1;
+               ++indx2;
             }
-            ++indx2;
+            indx0 += nxy;
          }
-         indx0 += nxy;
       }
      
       c3db::mygdevice.batch_cffty_stages_tmpy(0,c3db::fft_tag,false,ny,nn,2*nfft3d,tmp1,c3db::tmpy,da_indx);
       // for (auto i=0; i<nn; ++i)
       //    dcfftb_(&ny,tmp1+2*ny*i,c3db::tmpy);
- 
    }
+
    /*************************
     **** hilbert mapping ****
     *************************/
    else 
    {
-      c3db::c_ptranspose_ijk_end(nffts, nb, 2, tmp2, tmp1, request_indx);
+      c3db::c_ptranspose_ijk_end(nffts,nb,2,tmp2,tmp1,request_indx);
      
       /********************************************
        ***     do fft along ny dimension        ***
        ***   A(ky,nz,kx) <- fft1d[A(ny,nz,kx)]  ***
        ********************************************/
-      c3db::mygdevice.batch_cffty_stages_tmpy_zero(0,c3db::fft_tag,false,ny,nq2,2*nfft3d,tmp2,c3db::tmpy,zero_row2[nb],da_indx);
+      c3db::mygdevice.batch_cffty_stages_tmpy_zero(0,c3db::fft_tag,false,ny,nffts*nq2,2*nfft3d,tmp2,c3db::tmpy,zero_row2[nb],da_indx);
    }
 }
 
@@ -1989,17 +2020,20 @@ void CGrid::pfftby_compute(const int nffts, const int nb, double *tmp1, double *
        ***   A(kx,ky,nz) <- fft1d[A(kx,ny,nz)]  ***
        ********************************************/
       int indx0 = 0;
-      int indx2 = 0;
       int nn = 0;
-      for (auto q=0; q<nq; ++q) 
+      for (auto s=0; s<nffts; ++s)
       {
-         for (auto i=0; i<nx; ++i) 
+         int indx2 = 0;
+         for (auto q=0; q<nq; ++q) 
          {
-            if (!zero_row2[nb][indx2]) 
+            for (auto i=0; i<nx; ++i) 
             {
-               nn += 1;
+               if (!zero_row2[nb][indx2]) 
+               {
+                  nn += 1;
+               }
+               ++indx2;
             }
-            ++indx2;
          }
       }
      
@@ -2007,6 +2041,7 @@ void CGrid::pfftby_compute(const int nffts, const int nb, double *tmp1, double *
       // for (auto i=0; i<nn; ++i)
       //    dcfftb_(&ny,tmp1+2*ny*i,c3db::tmpy);
    }
+
    /*************************
     **** hilbert mapping ****
     *************************/
@@ -2040,17 +2075,20 @@ void CGrid::pfftby_end(const int nffts, const int nb, double *tmp1, double *tmp2
        ***   A(kx,ky,nz) <- fft1d[A(kx,ny,nz)]  ***
        ********************************************/
       int indx0 = 0;
-      int indx2 = 0;
       int nn = 0;
-      for (auto q=0; q<nq; ++q) 
+      for (auto s=0; s<nffts; ++s)
       {
-         for (auto i=0; i<nx; ++i) 
+         int indx2 = 0;
+         for (auto q=0; q<nq; ++q) 
          {
-            if (!zero_row2[nb][indx2]) 
+            for (auto i=0; i<nx; ++i) 
             {
-               nn += 1;
+               if (!zero_row2[nb][indx2]) 
+               {
+                  nn += 1;
+               }
+               ++indx2;
             }
-            ++indx2;
          }
       }
      
@@ -2059,31 +2097,36 @@ void CGrid::pfftby_end(const int nffts, const int nb, double *tmp1, double *tmp2
       //    dcfftb_(&ny,tmp1+2*ny*i,c3db::tmpy);
      
       indx0 = 0;
-      indx2 = 0;
       nn = 0;
-      for (auto q=0; q<nq; ++q) 
+      for (auto s=0; s<nffts; ++s)
       {
-         for (auto i=0; i<nx; ++i) 
+         int indx2 = 0;
+         for (auto q=0; q<nq; ++q) 
          {
-            if (!zero_row2[nb][indx2]) 
+            for (auto i=0; i<nx; ++i) 
             {
-               auto jj = 0;
-               auto indx3 = 2*i + indx0;
-               auto shift = 2*ny*nn;
-               for (auto j=0; j<ny; ++j) 
+               if (!zero_row2[nb][indx2]) 
                {
-                  tmp2[indx3]   = tmp1[jj   + shift];
-                  tmp2[indx3+1] = tmp1[jj+1 + shift];
-                  jj += 2;
-                  indx3 += nx;
+                  auto jj = 0;
+                  auto indx3 = 2*i + indx0;
+                  auto shift = 2*ny*nn;
+                  for (auto j=0; j<ny; ++j) 
+                  {
+                     tmp2[indx3]   = tmp1[jj   + shift];
+                     tmp2[indx3+1] = tmp1[jj+1 + shift];
+                     jj += 2;
+                     indx3 += nx;
+                  }
+                  nn += 1;
                }
-               nn += 1;
+               ++indx2;
             }
-            ++indx2;
+            indx0 += nxy;
          }
-         indx0 += nxy;
       }
+
    }
+
    /*************************
     **** hilbert mapping ****
     *************************/
@@ -2093,7 +2136,7 @@ void CGrid::pfftby_end(const int nffts, const int nb, double *tmp1, double *tmp2
       ***     do fft along ny dimension        ***
       ***   A(ky,nz,kx) <- fft1d[A(ny,nz,kx)]  ***
       ********************************************/
-      c3db::mygdevice.batch_cffty_stages_tmpy_zero(2,c3db::fft_tag,false,ny,nq2,2*nfft3d,tmp2,c3db::tmpy,zero_row2[nb],da_indx);
+      c3db::mygdevice.batch_cffty_stages_tmpy_zero(2,c3db::fft_tag,false,ny,nffts*nq2,2*nfft3d,tmp2,c3db::tmpy,zero_row2[nb],da_indx);
  
       c3db::c_ptranspose_ijk_start(nffts, nb, 3, tmp2, tmp1, tmp2, request_indx, 46);
       // c3db::c_ptranspose_ijk(nb,3,tmp2,tmp1,tmp2);
@@ -2119,7 +2162,7 @@ void CGrid::pfftbx_start(const int nffts, const int nb, double *tmp1, double *tm
        ***     do fft along kx dimension            ***
        ***   A(nx,ny,nz) <- fft1d^(-1)[A(kx,ny,nz)] ***
        ************************************************/
-      c3db::mygdevice.batch_cfftx_stages_tmpx(0,c3db::fft_tag,false, nx, ny * nq, 2*nfft3d, tmp2, c3db::tmpx,da_indx);
+      c3db::mygdevice.batch_cfftx_stages_tmpx(0,c3db::fft_tag,false,nx,nffts*ny*nq,2*nfft3d,tmp2,c3db::tmpx,da_indx);
    }
    /*************************
     **** hilbert mapping ****
@@ -2132,7 +2175,7 @@ void CGrid::pfftbx_start(const int nffts, const int nb, double *tmp1, double *tm
        ***     do fft along kx dimension            ***
        ***   A(nx,ny,nz) <- fft1d^(-1)[A(kx,ny,nz)] ***
        ************************************************/
-      c3db::mygdevice.batch_cfftx_stages_tmpx(0,c3db::fft_tag,false, nx, nq1, 2*nfft3d, tmp1, c3db::tmpx,da_indx);
+      c3db::mygdevice.batch_cfftx_stages_tmpx(0,c3db::fft_tag,false,nx,nffts*nq1,2*nfft3d,tmp1,c3db::tmpx,da_indx);
    }
 }
 
@@ -2153,7 +2196,7 @@ void CGrid::pfftbx_compute(const int nffts, const int nb, double *tmp1, double *
        ***     do fft along kx dimension            ***
        ***   A(nx,ny,nz) <- fft1d^(-1)[A(kx,ny,nz)] ***
        ************************************************/
-      c3db::mygdevice.batch_cfftx_stages_tmpx(1,c3db::fft_tag,false, nx, ny * nq, 2*nfft3d, tmp2, c3db::tmpx,da_indx);
+      c3db::mygdevice.batch_cfftx_stages_tmpx(1,c3db::fft_tag,false,nx,nffts*ny*nq,2*nfft3d,tmp2,c3db::tmpx,da_indx);
    }
 
    /*************************
@@ -2165,7 +2208,7 @@ void CGrid::pfftbx_compute(const int nffts, const int nb, double *tmp1, double *
        ***     do fft along kx dimension            ***
        ***   A(nx,ny,nz) <- fft1d^(-1)[A(kx,ny,nz)] ***
        ************************************************/
-      c3db::mygdevice.batch_cfftx_stages_tmpx(1,c3db::fft_tag,false, nx, nq1, 2*nfft3d, tmp1, c3db::tmpx,da_indx);
+      c3db::mygdevice.batch_cfftx_stages_tmpx(1,c3db::fft_tag,false,nx,nffts*nq1,2*nfft3d,tmp1,c3db::tmpx,da_indx);
    }
 }
 
@@ -2187,9 +2230,10 @@ void CGrid::pfftbx_end(const int nffts, const int nb, double *tmp1, double *tmp2
        ***     do fft along kx dimension            ***
        ***   A(nx,ny,nz) <- fft1d^(-1)[A(kx,ny,nz)] ***
        ************************************************/
-      c3db::mygdevice.batch_cfftx_stages_tmpx(2,c3db::fft_tag,false, nx, ny * nq, 2*nfft3d, tmp2, c3db::tmpx,da_indx);
-      std::memcpy(tmp1, tmp2, 2*nfft3d * sizeof(double));
+      c3db::mygdevice.batch_cfftx_stages_tmpx(2,c3db::fft_tag,false, nx, nffts*ny*nq,2*nfft3d,tmp2,c3db::tmpx,da_indx);
+      std::memcpy(tmp1,tmp2,nffts*2*nfft3d*sizeof(double));
    }
+
    /*************************
     **** hilbert mapping ****
     *************************/
@@ -2199,7 +2243,7 @@ void CGrid::pfftbx_end(const int nffts, const int nb, double *tmp1, double *tmp2
        ***     do fft along kx dimension            ***
        ***   A(nx,ny,nz) <- fft1d^(-1)[A(kx,ny,nz)] ***
        ************************************************/
-      c3db::mygdevice.batch_cfftx_stages_tmpx(2,c3db::fft_tag,false, nx, nq1, 2*nfft3d, tmp1, c3db::tmpx,da_indx);
+      c3db::mygdevice.batch_cfftx_stages_tmpx(2,c3db::fft_tag,false,nx,nffts*nq1,2*nfft3d,tmp1,c3db::tmpx,da_indx);
       if (2*nfft3d_map < 2*nfft3d)
          std::memset(tmp1 + 2*nfft3d_map, 0, (2*nfft3d - 2*nfft3d_map) * sizeof(double));
    }
@@ -2220,7 +2264,8 @@ void CGrid::pfftb_step12(const int step, const int nffts, const int nb, double *
    {
       // c3db::parall->astart(request_indx,parall->np_i());
       // unpack start, tmp1-->tmp1
-      std::memcpy(tmp1, a, 2 * (nidb[nb]) * sizeof(double));
+      for (auto s=0; s<nffts; ++s)
+         std::memcpy(tmp1 + s*n2ft3d, a + s*2*nidb[nb], 2*(nidb[nb])*sizeof(double));
       this->c_unpack_start(nffts, nb, tmp1, tmp2, request_indx, 47);
    } 
    else if (step == 1) 
@@ -2374,8 +2419,8 @@ void CGrid::pfftfx(const int nffts, const int nb, double *a, double *tmp1, doubl
    {
       // do fft along nx dimension
       //c3db::mygdevice.batch_cfftx_tmpx(c3db::fft_tag,true, nx, ny*nq, 2*nfft3d, a, c3db::tmpx);
-      c3db::mygdevice.batch_cfft(c3db::fft_tag,true,nx,ny*nq,nx,a,c3db::forward_x,c3db::tmpx,0);
-      std::memcpy(tmp1, a, 2*nfft3d * sizeof(double));
+      c3db::mygdevice.batch_cfft(c3db::fft_tag,true,nx,nffts*ny*nq,nx,a,c3db::forward_x,c3db::tmpx,0);
+      std::memcpy(tmp1, a, nffts*2*nfft3d * sizeof(double));
    }
    /**** hilbert mapping ****/
    else 
@@ -2383,8 +2428,8 @@ void CGrid::pfftfx(const int nffts, const int nb, double *a, double *tmp1, doubl
       // do fft along nx dimension
       // A(kx,ny,nz) <- fft1d[A(nx,ny,nz)]
       //c3db::mygdevice.batch_cfftx_tmpx(c3db::fft_tag,true, nx, nq1, 2*nfft3d, a, c3db::tmpx);
-      c3db::mygdevice.batch_cfft(c3db::fft_tag,true,nx,nq1,nx,a,c3db::forward_x,c3db::tmpx,0);
-      c3db::c_ptranspose_ijk_start(nffts, nb, 0, a, tmp1, tmp2, request_indx, 40);
+      c3db::mygdevice.batch_cfft(c3db::fft_tag,true,nx,nffts*nq1,nx,a,c3db::forward_x,c3db::tmpx,0);
+      c3db::c_ptranspose_ijk_start(nffts,nb,0,a,tmp1,tmp2,request_indx, 40);
    }
 }
 
@@ -2403,58 +2448,64 @@ void CGrid::pfftfy(const int nffts, const int nb, double *tmp1, double *tmp2, in
       // do fft along ny dimension
       // A(kx,ky,nz) <- fft1d[A(kx,ny,nz)]
       int indx0 = 0;
-      int indx2 = 0;
       int nn = 0;
-      for (auto q=0; q<nq; ++q) 
+      for (auto s=0; s<nffts; ++s)
       {
-         for (auto i=0; i<nx; ++i) 
+         int indx2 = 0;
+         for (auto q=0; q<nq; ++q) 
          {
-            if (!zero_row2[nb][indx2]) 
+            for (auto i=0; i<nx; ++i) 
             {
-               auto jj = 0;
-               auto indx3 = 2*i + indx0;
-               auto shift = 2*ny*nn;
-               for (auto j=0; j<ny; ++j) 
+               if (!zero_row2[nb][indx2]) 
                {
-                  tmp2[jj   + shift] = tmp1[indx3];
-                  tmp2[jj+1 + shift] = tmp1[indx3 + 1];
-                  jj += 2;
-                  indx3 += nx;
+                  auto jj = 0;
+                  auto indx3 = 2*i + indx0;
+                  auto shift = 2*ny*nn;
+                  for (auto j=0; j<ny; ++j) 
+                  {
+                     tmp2[jj   + shift] = tmp1[indx3];
+                     tmp2[jj+1 + shift] = tmp1[indx3 + 1];
+                     jj += 2;
+                     indx3 += nx;
+                  }
+                  ++nn;
                }
-               ++nn;
+               ++indx2;
             }
-            ++indx2;
+            indx0 += nxy;
          }
-         indx0 += nxy;
       }
      
       //c3db::mygdevice.batch_cffty_tmpy(c3db::fft_tag,true, ny, nn, 2*nfft3d, tmp2, c3db::tmpy);
       c3db::mygdevice.batch_cfft(c3db::fft_tag,true,ny,nn,ny,tmp2,c3db::forward_y,c3db::tmpy,1);
      
       indx0 = 0;
-      indx2 = 0;
       nn = 0;
-      for (auto q=0; q<nq; ++q) 
+      for (auto s=0; s<nffts; ++s)
       {
-         for (auto i=0; i<nx; ++i) 
+         int indx2 = 0;
+         for (auto q=0; q<nq; ++q) 
          {
-            if (!zero_row2[nb][indx2]) 
+            for (auto i=0; i<nx; ++i) 
             {
-               auto jj = 0;
-               auto indx3 = 2*i + indx0;
-               auto shift = 2*ny*nn;
-               for (auto j=0; j<ny; ++j) 
+               if (!zero_row2[nb][indx2]) 
                {
-                  tmp1[indx3] = tmp2[jj + shift];
-                  tmp1[indx3 + 1] = tmp2[jj + 1 + shift];
-                  jj += 2;
-                  indx3 += nx;
+                  auto jj = 0;
+                  auto indx3 = 2*i + indx0;
+                  auto shift = 2*ny*nn;
+                  for (auto j=0; j<ny; ++j) 
+                  {
+                     tmp1[indx3] = tmp2[jj + shift];
+                     tmp1[indx3 + 1] = tmp2[jj + 1 + shift];
+                     jj += 2;
+                     indx3 += nx;
+                  }
+                  ++nn;
                }
-               ++nn;
+               ++indx2;
             }
-            ++indx2;
+            indx0 += nxy;
          }
-         indx0 += nxy;
       }
      
       // Do a transpose of A
@@ -2471,10 +2522,10 @@ void CGrid::pfftfy(const int nffts, const int nb, double *tmp1, double *tmp2, in
       // do fft along ny dimension
       // A(ky,nz,kx) <- fft1d[A(ny,nz,kx)]
       //c3db::mygdevice.batch_cffty_tmpy_zero(c3db::fft_tag,true,ny,nq2,2*nfft3d,tmp1,c3db::tmpy,zero_row2[nb]);
-      c3db::mygdevice.batch_cfft_zero(c3db::fft_tag,true,ny,nq2,ny,tmp1,c3db::forward_y,c3db::tmpy,zero_row2[nb],1);
+      c3db::mygdevice.batch_cfft_zero(c3db::fft_tag,true,ny,nffts*nq2,ny,tmp1,c3db::forward_y,c3db::tmpy,zero_row2[nb],1);
      
       // in=tmp2, out=tmp2
-      c3db::c_ptranspose_ijk_start(nffts, nb, 1, tmp1, tmp2, tmp1, request_indx, 42);
+      c3db::c_ptranspose_ijk_start(nffts,nb,1,tmp1,tmp2,tmp1,request_indx, 42);
    }
 }
 
@@ -2496,59 +2547,66 @@ void CGrid::pfftfz(const int nffts, const int nb, double *tmp1, double *tmp2, in
       // do fft along nz dimension
       // A(kx,kz,ky) <- fft1d[A(kx,nz,ky)]
       int indx0 = 0;
-      int indx2 = 0;
       int nn = 0;
-      for (auto q=0; q<nq; ++q) 
+      for (auto s=0; s<nffts; ++s)
       {
-         for (auto i=0; i<nx; ++i) 
+         int indx2 = 0;
+         for (auto q=0; q<nq; ++q) 
          {
-            if (!zero_row3[nb][indx2]) 
+            for (auto i=0; i<nx; ++i) 
             {
-               auto kk = 0;
-               auto indx3 = 2*i + indx0;
-               auto shift = 2*nz*nn;
-               for (auto k=0; k<nz; ++k) 
+               if (!zero_row3[nb][indx2]) 
                {
-                  tmp1[kk   + shift] = tmp2[indx3];
-                  tmp1[kk+1 + shift] = tmp2[indx3 + 1];
-                  kk += 2;
-                  indx3 += nx;
+                  auto kk = 0;
+                  auto indx3 = 2*i + indx0;
+                  auto shift = 2*nz*nn;
+                  for (auto k=0; k<nz; ++k) 
+                  {
+                     tmp1[kk   + shift] = tmp2[indx3];
+                     tmp1[kk+1 + shift] = tmp2[indx3 + 1];
+                     kk += 2;
+                     indx3 += nx;
+                  }
+                  ++nn;
                }
-               ++nn;
+               ++indx2;
             }
-            ++indx2;
+            indx0 += nxz;
          }
-         indx0 += nxz;
       }
      
       //c3db::mygdevice.batch_cfftz_tmpz(c3db::fft_tag,true, nz, nn, 2*nfft3d, tmp1, c3db::tmpz);
       c3db::mygdevice.batch_cfft(c3db::fft_tag,true,nz,nn,nz,tmp1,c3db::forward_z,c3db::tmpz,2);
      
       indx0 = 0;
-      indx2 = 0;
       nn = 0;
-      for (auto q=0; q<nq; ++q) 
+      for (auto s=0; s<nffts; ++s)
       {
-         for (auto i=0; i<nx; ++i) 
+         int indx2 = 0;
+         for (auto q=0; q<nq; ++q) 
          {
-            if (!zero_row3[nb][indx2]) 
+            for (auto i=0; i<nx; ++i) 
             {
-               auto kk = 0;
-               auto indx3 = 2 * i + indx0;
-               auto shift = 2 * nz * nn;
-               for (auto k=0; k<nz; ++k) 
+               if (!zero_row3[nb][indx2]) 
                {
-                  tmp2[indx3]   = tmp1[kk   + shift];
-                  tmp2[indx3+1] = tmp1[kk+1 + shift];
-                  kk += 2;
-                  indx3 += nx;
+                  auto kk = 0;
+                  auto indx3 = 2 * i + indx0;
+                  auto shift = 2 * nz * nn;
+                  for (auto k=0; k<nz; ++k) 
+                  {
+                     tmp2[indx3]   = tmp1[kk   + shift];
+                     tmp2[indx3+1] = tmp1[kk+1 + shift];
+                     kk += 2;
+                     indx3 += nx;
+                  }
+                  ++nn;
                }
-               ++nn;
+               ++indx2;
             }
-            ++indx2;
+            indx0 += nxz;
          }
-         indx0 += nxz;
       }
+
    }
  
    /**** hilbert mapping ****/
@@ -2560,7 +2618,7 @@ void CGrid::pfftfz(const int nffts, const int nb, double *tmp1, double *tmp2, in
       // do fft along nz dimension
       // A(kz,kx,ky) <- fft1d[A(nz,kx,ky)]
       //c3db::mygdevice.batch_cfftz_tmpz_zero(c3db::fft_tag,true, nz, nq3, 2*nfft3d, tmp2, c3db::tmpz,zero_row3[nb]);
-      c3db::mygdevice.batch_cfft_zero(c3db::fft_tag,true,nz,nq3,nz,tmp2,c3db::forward_z,c3db::tmpz,zero_row3[nb],2);
+      c3db::mygdevice.batch_cfft_zero(c3db::fft_tag,true,nz,nffts*nq3,nz,tmp2,c3db::forward_z,c3db::tmpz,zero_row3[nb],2);
    }
 }
 
@@ -2605,16 +2663,16 @@ void CGrid::pfftfx_start(const int nffts, const int nb, double *a, double *tmp1,
    if (maptype == 1) 
    {
       // do fft along nx dimension
-      std::memcpy(tmp2, a, 2*nfft3d * sizeof(double));
-      c3db::mygdevice.batch_cfftx_stages_tmpx(0,c3db::fft_tag,true, nx, ny*nq, 2*nfft3d, tmp2, c3db::tmpx,da_indx);
+      std::memcpy(tmp2,a,nffts*2*nfft3d*sizeof(double));
+      c3db::mygdevice.batch_cfftx_stages_tmpx(0,c3db::fft_tag,true,nx,nffts*ny*nq,2*nfft3d,tmp2,c3db::tmpx,da_indx);
    }
    /**** hilbert mapping ****/
    else 
    {
       // do fft along nx dimension
       // A(kx,ny,nz) <- fft1d[A(nx,ny,nz)]
-      std::memcpy(tmp2, a, 2*nfft3d * sizeof(double));
-      c3db::mygdevice.batch_cfftx_stages_tmpx(0,c3db::fft_tag,true, nx, nq1, 2*nfft3d, tmp2, c3db::tmpx,da_indx);
+      std::memcpy(tmp2,a,nffts*2*nfft3d*sizeof(double));
+      c3db::mygdevice.batch_cfftx_stages_tmpx(0,c3db::fft_tag,true,nx,nffts*nq1,2*nfft3d,tmp2,c3db::tmpx,da_indx);
    }
 }
 
@@ -2629,14 +2687,15 @@ void CGrid::pfftfx_compute(const int nffts, const int nb, double *a, double *tmp
    if (maptype == 1) 
    {
       // do fft along nx dimension
-      c3db::mygdevice.batch_cfftx_stages_tmpx(1,c3db::fft_tag,true, nx, ny*nq, 2*nfft3d, tmp2, c3db::tmpx,da_indx);
+      c3db::mygdevice.batch_cfftx_stages_tmpx(1,c3db::fft_tag,true,nx,nffts*ny*nq,2*nfft3d,tmp2,c3db::tmpx,da_indx);
    }
+
    /**** hilbert mapping ****/
    else 
    {
       // do fft along nx dimension
       // A(kx,ny,nz) <- fft1d[A(nx,ny,nz)]
-      c3db::mygdevice.batch_cfftx_stages_tmpx(1,c3db::fft_tag,true, nx, nq1, 2*nfft3d, tmp2, c3db::tmpx,da_indx);
+      c3db::mygdevice.batch_cfftx_stages_tmpx(1,c3db::fft_tag,true,nx,nffts*nq1,2*nfft3d,tmp2,c3db::tmpx,da_indx);
    }
 }
 
@@ -2652,23 +2711,21 @@ void CGrid::pfftfx_end(const int nffts, const int nb, double *a, double *tmp1, d
    if (maptype == 1) 
    {
       // do fft along nx dimension
-      c3db::mygdevice.batch_cfftx_stages_tmpx(2,c3db::fft_tag,true, nx, ny*nq, 2*nfft3d, tmp2, c3db::tmpx,da_indx);
-      std::memcpy(tmp1, tmp2, 2*nfft3d * sizeof(double));
+      c3db::mygdevice.batch_cfftx_stages_tmpx(2,c3db::fft_tag,true,nx,nffts*ny*nq,2*nfft3d,tmp2,c3db::tmpx,da_indx);
+      std::memcpy(tmp1, tmp2, nffts*2*nfft3d*sizeof(double));
    }
    /**** hilbert mapping ****/
    else 
    {
       // do fft along nx dimension
       // A(kx,ny,nz) <- fft1d[A(nx,ny,nz)]
-      c3db::mygdevice.batch_cfftx_stages_tmpx(2,c3db::fft_tag,true, nx, nq1, 2*nfft3d, tmp2, c3db::tmpx,da_indx);
+      c3db::mygdevice.batch_cfftx_stages_tmpx(2,c3db::fft_tag,true,nx,nffts*nq1,2*nfft3d,tmp2,c3db::tmpx,da_indx);
  
-      c3db::c_ptranspose_ijk_start(nffts, nb, 0, tmp2, tmp1, tmp2, request_indx, 40);
+      c3db::c_ptranspose_ijk_start(nffts,nb,0,tmp2,tmp1,tmp2,request_indx, 40);
  
    }
 }
 
-
-
 /********************************
  *                              *
  *      CGrid::pfftfy_start     *
@@ -2684,29 +2741,32 @@ void CGrid::pfftfy_start(const int nffts, const int nb, double *tmp1, double *tm
       // do fft along ny dimension
       // A(kx,ky,nz) <- fft1d[A(kx,ny,nz)]
       int indx0 = 0;
-      int indx2 = 0;
       int nn = 0;
-      for (auto q=0; q<nq; ++q) 
+      for (auto s=0; s<nffts; ++s)
       {
-         for (auto i=0; i<nx; ++i) 
+         int indx2 = 0;
+         for (auto q=0; q<nq; ++q) 
          {
-            if (!zero_row2[nb][indx2]) 
+            for (auto i=0; i<nx; ++i) 
             {
-               auto jj = 0;
-               auto indx3 = 2*i + indx0;
-               auto shift = 2*ny*nn;
-               for (auto j=0; j<ny; ++j) 
+               if (!zero_row2[nb][indx2]) 
                {
-                  tmp2[jj   + shift] = tmp1[indx3];
-                  tmp2[jj+1 + shift] = tmp1[indx3 + 1];
-                  jj += 2;
-                  indx3 += nx;
+                  auto jj = 0;
+                  auto indx3 = 2*i + indx0;
+                  auto shift = 2*ny*nn;
+                  for (auto j=0; j<ny; ++j) 
+                  {
+                     tmp2[jj   + shift] = tmp1[indx3];
+                     tmp2[jj+1 + shift] = tmp1[indx3 + 1];
+                     jj += 2;
+                     indx3 += nx;
+                  }
+                  ++nn;
                }
-               ++nn;
+               ++indx2;
             }
-            ++indx2;
+            indx0 += nxy;
          }
-         indx0 += nxy;
       }
      
       c3db::mygdevice.batch_cffty_stages_tmpy(0,c3db::fft_tag,true, ny, nn, 2*nfft3d, tmp2, c3db::tmpy,da_indx);
@@ -2716,11 +2776,11 @@ void CGrid::pfftfy_start(const int nffts, const int nb, double *tmp1, double *tm
    else 
    {
       // in=tmp1, out=tmp2
-      c3db::c_ptranspose_ijk_end(nffts, nb, 0, tmp1, tmp2, request_indx);
+      c3db::c_ptranspose_ijk_end(nffts,nb,0,tmp1,tmp2,request_indx);
      
       // do fft along ny dimension
       // A(ky,nz,kx) <- fft1d[A(ny,nz,kx)]
-      c3db::mygdevice.batch_cffty_stages_tmpy_zero(0,c3db::fft_tag,true,ny,nq2,2*nfft3d,tmp1,c3db::tmpy,zero_row2[nb],da_indx);
+      c3db::mygdevice.batch_cffty_stages_tmpy_zero(0,c3db::fft_tag,true,ny,nffts*nq2,2*nfft3d,tmp1,c3db::tmpy,zero_row2[nb],da_indx);
    }
 }
 
@@ -2739,17 +2799,20 @@ void CGrid::pfftfy_compute(const int nffts, const int nb, double *tmp1, double *
       // do fft along ny dimension
       // A(kx,ky,nz) <- fft1d[A(kx,ny,nz)]
       int indx0 = 0;
-      int indx2 = 0;
       int nn = 0;
-      for (auto q=0; q<nq; ++q) 
+      for (auto s=0; s<nffts; ++s)
       {
-         for (auto i=0; i<nx; ++i) 
+         int indx2 = 0;
+         for (auto q=0; q<nq; ++q) 
          {
-            if (!zero_row2[nb][indx2]) 
+            for (auto i=0; i<nx; ++i) 
             {
-               ++nn;
+               if (!zero_row2[nb][indx2]) 
+               {
+                  ++nn;
+               }
+               ++indx2;
             }
-            ++indx2;
          }
       }
      
@@ -2762,13 +2825,10 @@ void CGrid::pfftfy_compute(const int nffts, const int nb, double *tmp1, double *
       // in=tmp1, out=tmp2
       // do fft along ny dimension
       // A(ky,nz,kx) <- fft1d[A(ny,nz,kx)]
-      c3db::mygdevice.batch_cffty_stages_tmpy_zero(1,c3db::fft_tag,true,ny,nq2,2*nfft3d,tmp1,c3db::tmpy,zero_row2[nb],da_indx);
+      c3db::mygdevice.batch_cffty_stages_tmpy_zero(1,c3db::fft_tag,true,ny,nffts*nq2,2*nfft3d,tmp1,c3db::tmpy,zero_row2[nb],da_indx);
    }
 }
 
-
-
-
 /********************************
  *                              *
  *      CGrid::pfftfy_end       *
@@ -2784,61 +2844,68 @@ void CGrid::pfftfy_end(const int nffts, const int nb, double *tmp1, double *tmp2
       // do fft along ny dimension
       // A(kx,ky,nz) <- fft1d[A(kx,ny,nz)]
       int indx0 = 0;
-      int indx2 = 0;
       int nn = 0;
-      for (auto q=0; q<nq; ++q) 
+      for (auto s=0; s<nffts; ++s)
       {
-         for (auto i=0; i<nx; ++i) 
+         int indx2 = 0;
+         for (auto q=0; q<nq; ++q) 
          {
-            if (!zero_row2[nb][indx2]) 
+            for (auto i=0; i<nx; ++i) 
             {
-               ++nn;
+               if (!zero_row2[nb][indx2]) 
+               {
+                  ++nn;
+               }
+               ++indx2;
             }
-            ++indx2;
          }
       }
      
       c3db::mygdevice.batch_cffty_stages_tmpy(2,c3db::fft_tag,true, ny, nn, 2*nfft3d, tmp2, c3db::tmpy,da_indx);
      
       indx0 = 0;
-      indx2 = 0;
       nn = 0;
-      for (auto q=0; q<nq; ++q) 
+      for (auto s=0; s<nffts; ++s)
       {
-         for (auto i=0; i<nx; ++i) 
+         int indx2 = 0;
+         for (auto q=0; q<nq; ++q) 
          {
-            if (!zero_row2[nb][indx2]) 
+            for (auto i=0; i<nx; ++i) 
             {
-               auto jj = 0;
-               auto indx3 = 2*i + indx0;
-               auto shift = 2*ny*nn;
-               for (auto j=0; j<ny; ++j) 
+               if (!zero_row2[nb][indx2]) 
                {
-                 tmp1[indx3]   = tmp2[jj   + shift];
-                 tmp1[indx3+1] = tmp2[jj+1 + shift];
-                 jj += 2;
-                 indx3 += nx;
+                  auto jj = 0;
+                  auto indx3 = 2*i + indx0;
+                  auto shift = 2*ny*nn;
+                  for (auto j=0; j<ny; ++j) 
+                  {
+                    tmp1[indx3]   = tmp2[jj   + shift];
+                    tmp1[indx3+1] = tmp2[jj+1 + shift];
+                    jj += 2;
+                    indx3 += nx;
+                  }
+                  ++nn;
                }
-               ++nn;
+               ++indx2;
             }
-            ++indx2;
+            indx0 += nxy;
          }
-         indx0 += nxy;
       }
      
       // Do a transpose of A
       // A(ky,nz,ky) <- A(kx,ky,nz)
       c3db::c_ptranspose2_jk_start(nffts, nb, tmp1, tmp2, tmp1, request_indx, 41);
    }
+
    /**** hilbert mapping ****/
    else 
    {
       // do fft along ny dimension
       // A(ky,nz,kx) <- fft1d[A(ny,nz,kx)]
-      c3db::mygdevice.batch_cffty_stages_tmpy_zero(2,c3db::fft_tag,true,ny,nq2,2*nfft3d,tmp1,c3db::tmpy,zero_row2[nb],da_indx);
+      c3db::mygdevice.batch_cffty_stages_tmpy_zero(2,c3db::fft_tag,true,ny,nffts*nq2,2*nfft3d,tmp1,c3db::tmpy,zero_row2[nb],da_indx);
      
       // in=tmp2, out=tmp2
-      c3db::c_ptranspose_ijk_start(nffts, nb, 1, tmp1, tmp2, tmp1, request_indx, 42);
+      c3db::c_ptranspose_ijk_start(nffts,nb,1,tmp1,tmp2,tmp1,request_indx, 42);
    }
 }
 
@@ -2861,33 +2928,37 @@ void CGrid::pfftfz_start(const int nffts, const int nb, double *tmp1, double *tm
       // do fft along nz dimension
       // A(kx,kz,ky) <- fft1d[A(kx,nz,ky)]
       int indx0 = 0;
-      int indx2 = 0;
       int nn = 0;
-      for (auto q=0; q<nq; ++q) 
+      for (auto s=0; s<nffts; ++s)
       {
-         for (auto i=0; i<nx; ++i) 
+         int indx2 = 0;
+         for (auto q=0; q<nq; ++q) 
          {
-            if (!zero_row3[nb][indx2]) 
+            for (auto i=0; i<nx; ++i) 
             {
-               auto kk = 0;
-               auto indx3 = 2*i + indx0;
-               auto shift = 2*nz*nn;
-               for (auto k=0; k<nz; ++k) 
+               if (!zero_row3[nb][indx2]) 
                {
-                  tmp1[kk   + shift] = tmp2[indx3];
-                  tmp1[kk+1 + shift] = tmp2[indx3 + 1];
-                  kk += 2;
-                  indx3 += nx;
+                  auto kk = 0;
+                  auto indx3 = 2*i + indx0;
+                  auto shift = 2*nz*nn;
+                  for (auto k=0; k<nz; ++k) 
+                  {
+                     tmp1[kk   + shift] = tmp2[indx3];
+                     tmp1[kk+1 + shift] = tmp2[indx3 + 1];
+                     kk += 2;
+                     indx3 += nx;
+                  }
+                  ++nn;
                }
-               ++nn;
+               ++indx2;
             }
-            ++indx2;
+            indx0 += nxz;
          }
-         indx0 += nxz;
       }
      
-      c3db::mygdevice.batch_cfftz_stages_tmpz(0,c3db::fft_tag,true, nz, nn, 2*nfft3d, tmp1, c3db::tmpz,da_indx);
+      c3db::mygdevice.batch_cfftz_stages_tmpz(0,c3db::fft_tag,true,nz,nn,2*nfft3d,tmp1,c3db::tmpz,da_indx);
    }
+
    /**** hilbert mapping ****/
    else 
    {
@@ -2896,7 +2967,7 @@ void CGrid::pfftfz_start(const int nffts, const int nb, double *tmp1, double *tm
      
       // do fft along nz dimension
       // A(kz,kx,ky) <- fft1d[A(nz,kx,ky)]
-      c3db::mygdevice.batch_cfftz_stages_tmpz_zero(0,c3db::fft_tag,true, nz, nq3, 2*nfft3d, tmp2, c3db::tmpz, zero_row3[nb],da_indx);
+      c3db::mygdevice.batch_cfftz_stages_tmpz_zero(0,c3db::fft_tag,true,nz,nffts*nq3,2*nfft3d,tmp2,c3db::tmpz,zero_row3[nb],da_indx);
    }
 }
 
@@ -2916,27 +2987,32 @@ void CGrid::pfftfz_compute(const int nffts, const int nb, double *tmp1, double *
       // do fft along nz dimension
       // A(kx,kz,ky) <- fft1d[A(kx,nz,ky)]
       int indx0 = 0;
-      int indx2 = 0;
       int nn = 0;
-      for (auto q=0; q<nq; ++q) 
+      for (auto s=0; s<nffts; ++s)
       {
-         for (auto i=0; i<nx; ++i) 
+         int indx2 = 0;
+         for (auto q=0; q<nq; ++q) 
          {
-            if (!zero_row3[nb][indx2]) 
+            for (auto i=0; i<nx; ++i) 
             {
-               ++nn;
+               if (!zero_row3[nb][indx2]) 
+               {
+                  ++nn;
+               }
+               ++indx2;
             }
-            ++indx2;
          }
       }
-      c3db::mygdevice.batch_cfftz_stages_tmpz(1,c3db::fft_tag,true, nz, nn, 2*nfft3d, tmp1, c3db::tmpz,da_indx);
+
+      c3db::mygdevice.batch_cfftz_stages_tmpz(1,c3db::fft_tag,true,nz,nn,2*nfft3d,tmp1,c3db::tmpz,da_indx);
    }
+
    /**** hilbert mapping ****/
    else 
    {
       // do fft along nz dimension
       // A(kz,kx,ky) <- fft1d[A(nz,kx,ky)]
-      c3db::mygdevice.batch_cfftz_stages_tmpz_zero(1,c3db::fft_tag,true, nz, nq3, 2*nfft3d, tmp2, c3db::tmpz, zero_row3[nb],da_indx);
+      c3db::mygdevice.batch_cfftz_stages_tmpz_zero(1,c3db::fft_tag,true,nz,nffts*nq3,2*nfft3d,tmp2,c3db::tmpz,zero_row3[nb],da_indx);
    }
 }
 
@@ -2955,47 +3031,54 @@ void CGrid::pfftfz_end(const int nffts, const int nb, double *tmp1, double *tmp2
       // do fft along nz dimension
       // A(kx,kz,ky) <- fft1d[A(kx,nz,ky)]
       int indx0 = 0;
-      int indx2 = 0;
       int nn = 0;
-      for (auto q=0; q<nq; ++q) 
+      for (auto s=0; s<nffts; ++s)
       {
-         for (auto i=0; i<nx; ++i) 
+         int indx2 = 0;
+         for (auto q=0; q<nq; ++q) 
          {
-            if (!zero_row3[nb][indx2]) 
+            for (auto i=0; i<nx; ++i) 
             {
-               ++nn;
+               if (!zero_row3[nb][indx2]) 
+               {
+                  ++nn;
+               }
+               ++indx2;
             }
-            ++indx2;
          }
       }
      
       c3db::mygdevice.batch_cfftz_stages_tmpz(2,c3db::fft_tag,true, nz, nn, 2*nfft3d, tmp1, c3db::tmpz,da_indx);
      
       indx0 = 0;
-      indx2 = 0;
       nn = 0;
-      for (auto q=0; q<nq; ++q) 
+      for (auto s=0; s<nffts; ++s)
       {
-         for (auto i=0; i<nx; ++i) 
+         int indx2 = 0;
+         for (auto q=0; q<nq; ++q) 
          {
-            if (!zero_row3[nb][indx2]) 
+            for (auto i=0; i<nx; ++i) 
             {
-               auto kk = 0;
-               auto indx3 = 2*i + indx0;
-               auto shift = 2*nz*nn;
-               for (auto k=0; k<nz; ++k) 
+               if (!zero_row3[nb][indx2]) 
                {
-                  tmp2[indx3]   = tmp1[kk   + shift];
-                  tmp2[indx3+1] = tmp1[kk+1 + shift];
-                  kk += 2;
-                  indx3 += nx;
+                  auto kk = 0;
+                  auto indx3 = 2*i + indx0;
+                  auto shift = 2*nz*nn;
+                  for (auto k=0; k<nz; ++k) 
+                  {
+                     tmp2[indx3]   = tmp1[kk   + shift];
+                     tmp2[indx3+1] = tmp1[kk+1 + shift];
+                     kk += 2;
+                     indx3 += nx;
+                  }
+                  ++nn;
                }
-               ++nn;
+               ++indx2;
             }
-            ++indx2;
+            indx0 += nxz;
          }
-         indx0 += nxz;
       }
+
    }
  
    /**** hilbert mapping ****/
@@ -3080,13 +3163,14 @@ void CGrid::c_pack_start(const int nffts, const int nb, double *a, double *tmp1,
   // int one=1;
 
   // DCOPY_PWDFT(2*nfft3d,a,one,tmp,one);
-  std::memcpy(tmp1, a, 2*nfft3d * sizeof(double));
-  std::memset(a, 0, 2*nfft3d * sizeof(double));
+  std::memcpy(tmp1, a, nffts*2*nfft3d * sizeof(double));
+  std::memset(a, 0, nffts*2*nfft3d * sizeof(double));
 
-  c_aindexcopy(nidb2[nb], packarray[nb], tmp1, a);
+  for (auto s=0; s<nffts; ++s)
+     c_aindexcopy(nidb2[nb], packarray[nb], tmp1 + s*n2ft3d, a + s*n2ft3d);
 
   if (balanced)
-    mybalance->c_balance_start(nffts, nb, a, request_indx, msgtype);
+     mybalance->c_balance_start(nffts, nb, a, request_indx, msgtype);
 
   return;
 }
diff --git a/Nwpw/nwpwlib/C3dB/c3db.cpp b/Nwpw/nwpwlib/C3dB/c3db.cpp
index a29bd4da..9c93dcf2 100644
--- a/Nwpw/nwpwlib/C3dB/c3db.cpp
+++ b/Nwpw/nwpwlib/C3dB/c3db.cpp
@@ -4181,10 +4181,11 @@ void c3db::c_ptranspose1_jk_start(const int nffts, const int nb, double *a, doub
    int msglen;
    int n1 = p_i1_start[nb][0][np];
  
-   c_aindexcopy(n1, p_iq_to_i1[nb][0], a, tmp1);
+   for (auto s=0; s<nffts; ++s)
+      c_aindexcopy_stride(nffts,n1, p_iq_to_i1[nb][0], a + s*n2ft3d, tmp1 + 2*s);
  
    /* it = 0, transpose data on same thread */
-   msglen = 2*(p_i2_start[nb][0][1] - p_i2_start[nb][0][0]);
+   msglen = nffts*2*(p_i2_start[nb][0][1] - p_i2_start[nb][0][0]);
    std::memcpy(tmp2 + 2*p_i2_start[nb][0][0], tmp1 + 2*p_i1_start[nb][0][0], msglen*sizeof(double));
  
    /* receive packed array data */
@@ -4192,7 +4193,7 @@ void c3db::c_ptranspose1_jk_start(const int nffts, const int nb, double *a, doub
    {
       /* synchronous receive of tmp */
       int proc_from = (taskid - it + np) % np;
-      msglen = 2*(p_i2_start[nb][0][it+1] - p_i2_start[nb][0][it]);
+      msglen = nffts*2*(p_i2_start[nb][0][it+1] - p_i2_start[nb][0][it]);
       if (msglen > 0)
          parall->adreceive(request_indx, msgtype, proc_from, msglen, tmp2 + 2*p_i2_start[nb][0][it]);
      // parall->adreceive(request_indx,msgtype,proc_from,msglen,&tmp2[2*p_i2_start[nb][0][it]]);
@@ -4200,7 +4201,7 @@ void c3db::c_ptranspose1_jk_start(const int nffts, const int nb, double *a, doub
    for (auto it=1; it<np; ++it) 
    {
       int proc_to = (taskid + it) % np;
-      msglen = 2*(p_i1_start[nb][0][it+1] - p_i1_start[nb][0][it]);
+      msglen = nffts*2*(p_i1_start[nb][0][it+1] - p_i1_start[nb][0][it]);
       if (msglen > 0)
          parall->adsend(request_indx, msgtype, proc_to, msglen, tmp1 + 2*p_i1_start[nb][0][it]);
      // parall->adsend(request_indx,msgtype,proc_to,msglen,&tmp1[2*p_i1_start[nb][0][it]]);
@@ -4217,8 +4218,11 @@ void c3db::c_ptranspose1_jk_end(const int nffts, const int nb, double *a, double
    parall->awaitall(request_indx);
 
    int n2 = p_i2_start[nb][0][np];
-   c_bindexcopy(n2, p_iq_to_i2[nb][0], tmp2, a);
-   c_bindexzero(nfft3d - n2, p_iz_to_i2[nb][0], a);
+   for (auto s=0; s<nffts; ++s)
+   {
+      c_bindexcopy_stride(nffts,n2, p_iq_to_i2[nb][0], tmp2 + 2*s, a + s*n2ft3d);
+      c_bindexzero(nfft3d - n2, p_iz_to_i2[nb][0], a + s*n2ft3d);
+   }
 }
 
 /**************************************
@@ -4233,26 +4237,26 @@ void c3db::c_ptranspose2_jk_start(const int nffts, const int nb, double *a, doub
    int msglen;
    int n1 = p_j1_start[nb][0][np];
  
-   c_aindexcopy(n1, p_jq_to_i1[nb][0], a, tmp1);
+   for (auto s=0; s<nffts; ++s)
+      c_aindexcopy_stride(nffts,n1, p_jq_to_i1[nb][0], a+s*n2ft3d, tmp1+2*s);
  
    /* it = 0, transpose data on same thread */
-   msglen = 2*(p_j2_start[nb][0][1] - p_j2_start[nb][0][0]);
-   std::memcpy(tmp2 + 2 * p_j2_start[nb][0][0], tmp1 + 2 * p_j1_start[nb][0][0],
-               msglen * sizeof(double));
+   msglen = nffts*2*(p_j2_start[nb][0][1] - p_j2_start[nb][0][0]);
+   std::memcpy(tmp2 + 2*p_j2_start[nb][0][0], tmp1 + 2*p_j1_start[nb][0][0], msglen*sizeof(double));
  
    /* receive packed array data */
    for (auto it=1; it<np; ++it) 
    {
       /* synchronous receive of tmp */
       int proc_from = (taskid - it + np) % np;
-      msglen = 2 * (p_j2_start[nb][0][it+1] - p_j2_start[nb][0][it]);
+      msglen = nffts*2*(p_j2_start[nb][0][it+1] - p_j2_start[nb][0][it]);
       if (msglen > 0)
          parall->adreceive(request_indx, msgtype, proc_from, msglen, tmp2 + 2*p_j2_start[nb][0][it]);
    }
    for (auto it = 1; it < np; ++it) 
    {
       int proc_to = (taskid + it) % np;
-      msglen = 2 * (p_j1_start[nb][0][it + 1] - p_j1_start[nb][0][it]);
+      msglen = nffts*2*(p_j1_start[nb][0][it + 1] - p_j1_start[nb][0][it]);
       if (msglen > 0)
          parall->adsend(request_indx, msgtype, proc_to, msglen, tmp1 + 2*p_j1_start[nb][0][it]);
    }
@@ -4269,8 +4273,11 @@ void c3db::c_ptranspose2_jk_end(const int nffts, const int nb, double *a, double
    parall->awaitall(request_indx);
 
    int n2 = p_j2_start[nb][0][np];
-   c_bindexcopy(n2, p_jq_to_i2[nb][0], tmp2, a);
-   c_bindexzero(nfft3d - n2, p_jz_to_i2[nb][0], a);
+   for (auto s=0; s<nffts; ++s)
+   {
+      c_bindexcopy_stride(nffts,n2, p_jq_to_i2[nb][0], tmp2+2*s, a+s*n2ft3d);
+      c_bindexzero(nfft3d - n2, p_jz_to_i2[nb][0], a + s*n2ft3d);
+   }
 }
 
 /**************************************
@@ -4287,31 +4294,32 @@ void c3db::c_ptranspose_ijk_start(const int nffts, const int nb, const int op, d
    int n1 = p_i1_start[nb][op][np];
  
    /* pack a array - tmp1->tmp2 */
-   c_aindexcopy(n1, p_iq_to_i1[nb][op], a, tmp1);
+   for (auto s=0; s<nffts; ++s)
+      c_aindexcopy_stride(nffts,n1, p_iq_to_i1[nb][op], a + s*n2ft3d, tmp1 + 2*s);
  
    /* it = 0, transpose data on same thread  - tmp2->tmp1*/
-   msglen = 2*(p_i2_start[nb][op][1] - p_i2_start[nb][op][0]);
-   std::memcpy(tmp2 + 2 * p_i2_start[nb][op][0],
-               tmp1 + 2 * p_i1_start[nb][op][0], msglen * sizeof(double));
+   msglen = nffts*2*(p_i2_start[nb][op][1] - p_i2_start[nb][op][0]);
+   std::memcpy(tmp2 + 2*p_i2_start[nb][op][0],
+               tmp1 + 2*p_i1_start[nb][op][0], msglen*sizeof(double));
  
    /* receive packed array data */
    for (auto it=1; it<np; ++it) 
    {
       /* synchronous receive of tmp */
       int proc_from = (taskid - it + np) % np;
-      msglen = 2*(p_i2_start[nb][op][it+1] - p_i2_start[nb][op][it]);
+      msglen = nffts*2*(p_i2_start[nb][op][it+1] - p_i2_start[nb][op][it]);
       if (msglen > 0)
-        parall->adreceive(request_indx, msgtype, proc_from, msglen,
-                          tmp2 + 2 * p_i2_start[nb][op][it]);
+         parall->adreceive(request_indx, msgtype, proc_from, msglen,
+                           tmp2 + 2*p_i2_start[nb][op][it]);
       // parall->adreceive(request_indx,msgtype,proc_from,msglen,&tmp2[2*p_i2_start[nb][op][it]]);
    }
    for (auto it=1; it<np; ++it) 
    {
       int proc_to = (taskid + it) % np;
-      msglen = 2*(p_i1_start[nb][op][it+1] - p_i1_start[nb][op][it]);
+      msglen = nffts*2*(p_i1_start[nb][op][it+1] - p_i1_start[nb][op][it]);
       if (msglen > 0)
-        parall->adsend(request_indx, msgtype, proc_to, msglen,
-                       tmp1 + 2 * p_i1_start[nb][op][it]);
+         parall->adsend(request_indx, msgtype, proc_to, msglen,
+                        tmp1 + 2*p_i1_start[nb][op][it]);
       // parall->adsend(request_indx,msgtype,proc_to,msglen,&tmp1[2*p_i1_start[nb][op][it]]);
    }
 }
@@ -4332,10 +4340,10 @@ void c3db::c_ptranspose_ijk_end(const int nffts, const int nb, const int op, dou
    parall->awaitall(request_indx);
 
    /* unpack a array */
-   for (auto i=0; i<nffts; ++i)
+   for (auto s=0; s<nffts; ++s)
    {
-      c_bindexcopy(n2, p_iq_to_i2[nb][op], tmp2 + i*n2ft3d, a + i*n2ft3d);
-      c_bindexzero(n3, p_iz_to_i2[nb][op], a + i*n2ft3d);
+      c_bindexcopy_stride(nffts,n2, p_iq_to_i2[nb][op], tmp2 + 2*s, a + s*n2ft3d);
+      c_bindexzero(n3, p_iz_to_i2[nb][op], a + s*n2ft3d);
    }
 }
 
diff --git a/Nwpw/nwpwlib/D3dB/d3db.cpp b/Nwpw/nwpwlib/D3dB/d3db.cpp
index 554732f8..95f76cda 100644
--- a/Nwpw/nwpwlib/D3dB/d3db.cpp
+++ b/Nwpw/nwpwlib/D3dB/d3db.cpp
@@ -3773,34 +3773,33 @@ void d3db::c_ptranspose1_jk_start(const int nffts, const int nb, double *a, doub
    int msglen;
  
    int n1 = p_i1_start[nb][0][np];
- 
-   c_aindexcopy(n1, p_iq_to_i1[nb][0], a, tmp1);
+   for (auto s=0; s<nffts; ++s) 
+      c_aindexcopy_stride(nffts, n1, p_iq_to_i1[nb][0], a + s*n2ft3d, tmp1 + 2*s);
  
    /* it = 0, transpose data on same thread */
-   msglen = 2 * (p_i2_start[nb][0][1] - p_i2_start[nb][0][0]);
+   msglen = nffts*2*(p_i2_start[nb][0][1] - p_i2_start[nb][0][0]);
    // int one=1;
    // DCOPY_PWDFT(msglen,&(tmp1[2*p_i1_start[nb][0][0]]),one,&(tmp2[2*p_i2_start[nb][0][0]]),one);
-   std::memcpy(tmp2 + 2 * p_i2_start[nb][0][0], tmp1 + 2 * p_i1_start[nb][0][0],
-               msglen * sizeof(double));
+   std::memcpy(tmp2 + 2*p_i2_start[nb][0][0], tmp1 + 2*p_i1_start[nb][0][0], msglen*sizeof(double));
  
    /* receive packed array data */
-   for (it = 1; it < np; ++it) 
+   for (it=1; it<np; ++it) 
    {
       /* synchronous receive of tmp */
       proc_from = (taskid - it + np) % np;
-      msglen = 2 * (p_i2_start[nb][0][it + 1] - p_i2_start[nb][0][it]);
+      msglen = nffts*2*(p_i2_start[nb][0][it + 1] - p_i2_start[nb][0][it]);
       if (msglen > 0)
-        parall->adreceive(request_indx, msgtype, proc_from, msglen,
-                          tmp2 + 2 * p_i2_start[nb][0][it]);
+         parall->adreceive(request_indx, msgtype, proc_from, msglen,
+                           tmp2 + 2 * p_i2_start[nb][0][it]);
       // parall->adreceive(request_indx,msgtype,proc_from,msglen,&tmp2[2*p_i2_start[nb][0][it]]);
    }
-   for (it = 1; it < np; ++it) 
+   for (it=1; it<np; ++it) 
    {
       proc_to = (taskid + it) % np;
-      msglen = 2 * (p_i1_start[nb][0][it + 1] - p_i1_start[nb][0][it]);
+      msglen = nffts*2*(p_i1_start[nb][0][it + 1] - p_i1_start[nb][0][it]);
       if (msglen > 0)
-        parall->adsend(request_indx, msgtype, proc_to, msglen,
-                       tmp1 + 2 * p_i1_start[nb][0][it]);
+         parall->adsend(request_indx, msgtype, proc_to, msglen,
+                        tmp1 + 2 * p_i1_start[nb][0][it]);
       // parall->adsend(request_indx,msgtype,proc_to,msglen,&tmp1[2*p_i1_start[nb][0][it]]);
    }
 }
@@ -3815,8 +3814,11 @@ void d3db::c_ptranspose1_jk_end(const int nffts, const int nb, double *a, double
    parall->awaitall(request_indx);
 
    int n2 = p_i2_start[nb][0][np];
-   c_bindexcopy(n2, p_iq_to_i2[nb][0], tmp2, a);
-   c_bindexzero(nfft3d - n2, p_iz_to_i2[nb][0], a);
+   for (auto s=0; s<nffts; ++s)
+   {
+      c_bindexcopy_stride(nffts,n2, p_iq_to_i2[nb][0], tmp2 + 2*s, a+2*n2ft3d);
+      c_bindexzero(nfft3d - n2, p_iz_to_i2[nb][0], a+2*n2ft3d);
+   }
 }
 
 /**************************************
@@ -3826,37 +3828,40 @@ void d3db::c_ptranspose1_jk_end(const int nffts, const int nb, double *a, double
  **************************************/
 void d3db::c_ptranspose2_jk_start(const int nffts, const int nb, double *a, double *tmp1,
                                   double *tmp2, const int request_indx,
-                                  const int msgtype) {
-  int it, proc_from, proc_to;
-  int msglen;
-
-  int n1 = p_j1_start[nb][0][np];
-
-  c_aindexcopy(n1, p_jq_to_i1[nb][0], a, tmp1);
-
-  /* it = 0, transpose data on same thread */
-  msglen = 2 * (p_j2_start[nb][0][1] - p_j2_start[nb][0][0]);
-  // int one=1;
-  // DCOPY_PWDFT(msglen,&(tmp1[2*p_j1_start[nb][0][0]]),one,&(tmp2[2*p_j2_start[nb][0][0]]),one);
-  std::memcpy(tmp2 + 2 * p_j2_start[nb][0][0], tmp1 + 2 * p_j1_start[nb][0][0],
-              msglen * sizeof(double));
-
-  /* receive packed array data */
-  for (it = 1; it < np; ++it) {
-    /* synchronous receive of tmp */
-    proc_from = (taskid - it + np) % np;
-    msglen = 2 * (p_j2_start[nb][0][it + 1] - p_j2_start[nb][0][it]);
-    if (msglen > 0)
-      parall->adreceive(request_indx, msgtype, proc_from, msglen,
-                        &tmp2[2 * p_j2_start[nb][0][it]]);
-  }
-  for (it = 1; it < np; ++it) {
-    proc_to = (taskid + it) % np;
-    msglen = 2 * (p_j1_start[nb][0][it + 1] - p_j1_start[nb][0][it]);
-    if (msglen > 0)
-      parall->adsend(request_indx, msgtype, proc_to, msglen,
-                     &tmp1[2 * p_j1_start[nb][0][it]]);
-  }
+                                  const int msgtype) 
+{
+   int it, proc_from, proc_to;
+   int msglen;
+ 
+   int n1 = p_j1_start[nb][0][np];
+ 
+   for (auto s=0; s<nffts; ++s)
+      c_aindexcopy_stride(nffts, n1, p_jq_to_i1[nb][0], a + s*n2ft3d, tmp1 + 2*s);
+ 
+   /* it = 0, transpose data on same thread */
+   msglen = nffts*2*(p_j2_start[nb][0][1] - p_j2_start[nb][0][0]);
+   // int one=1;
+   // DCOPY_PWDFT(msglen,&(tmp1[2*p_j1_start[nb][0][0]]),one,&(tmp2[2*p_j2_start[nb][0][0]]),one);
+   std::memcpy(tmp2 + 2*p_j2_start[nb][0][0], tmp1 + 2*p_j1_start[nb][0][0], msglen*sizeof(double));
+ 
+   /* receive packed array data */
+   for (it=1; it<np; ++it) 
+   {
+      /* synchronous receive of tmp */
+      proc_from = (taskid - it + np) % np;
+      msglen = nffts*2*(p_j2_start[nb][0][it+1] - p_j2_start[nb][0][it]);
+      if (msglen > 0)
+         parall->adreceive(request_indx, msgtype, proc_from, msglen,
+                           tmp2+2*p_j2_start[nb][0][it]);
+   }
+   for (it=1; it<np; ++it) 
+   {
+      proc_to = (taskid + it) % np;
+      msglen = nffts*2*(p_j1_start[nb][0][it + 1] - p_j1_start[nb][0][it]);
+      if (msglen>0)
+         parall->adsend(request_indx, msgtype, proc_to, msglen,
+                        tmp1 + 2*p_j1_start[nb][0][it]);
+   }
 }
 
 /**************************************
@@ -3865,12 +3870,16 @@ void d3db::c_ptranspose2_jk_start(const int nffts, const int nb, double *a, doub
  *                                    *
  **************************************/
 void d3db::c_ptranspose2_jk_end(const int nffts, const int nb, double *a, double *tmp2,
-                                const int request_indx) {
-  parall->awaitall(request_indx);
+                                const int request_indx) 
+{
+   parall->awaitall(request_indx);
 
-  int n2 = p_j2_start[nb][0][np];
-  c_bindexcopy(n2, p_jq_to_i2[nb][0], tmp2, a);
-  c_bindexzero(nfft3d - n2, p_jz_to_i2[nb][0], a);
+   int n2 = p_j2_start[nb][0][np];
+   for (auto s=0; s<nffts; ++s)
+   {
+      c_bindexcopy_stride(nffts,n2, p_jq_to_i2[nb][0], tmp2 + 2*s, a + s*n2ft3d);
+      c_bindexzero(nfft3d - n2, p_jz_to_i2[nb][0], a + s*n2ft3d);
+   }
 }
 
 /**************************************
@@ -3888,34 +3897,35 @@ void d3db::c_ptranspose_ijk_start(const int nffts, const int nb, const int op, d
    int n1 = p_i1_start[nb][op][np];
  
    /* pack a array - tmp1->tmp2 */
-   c_aindexcopy(n1, p_iq_to_i1[nb][op], a, tmp1);
+   for (auto s=0; s<nffts; ++s)
+      c_aindexcopy_stride(nffts,n1, p_iq_to_i1[nb][op], a + s*n2ft3d, tmp1 + 2*s);
  
    /* it = 0, transpose data on same thread  - tmp2->tmp1*/
-   msglen = 2 * (p_i2_start[nb][op][1] - p_i2_start[nb][op][0]);
+   msglen = nffts*2*(p_i2_start[nb][op][1] - p_i2_start[nb][op][0]);
    // int one=1;
    // DCOPY_PWDFT(msglen,&(tmp1[2*p_i1_start[nb][op][0]]),one,&(tmp2[2*p_i2_start[nb][op][0]]),one);
    // std::memcpy(&(tmp2[2*p_i2_start[nb][op][0]]),&(tmp1[2*p_i1_start[nb][op][0]]),msglen*sizeof(double));
-   std::memcpy(tmp2 + 2 * p_i2_start[nb][op][0],
-               tmp1 + 2 * p_i1_start[nb][op][0], msglen * sizeof(double));
+   std::memcpy(tmp2 + 2*p_i2_start[nb][op][0],
+               tmp1 + 2*p_i1_start[nb][op][0], msglen*sizeof(double));
  
    /* receive packed array data */
-   for (it = 1; it < np; ++it) 
+   for (it=1; it<np; ++it) 
    {
       /* synchronous receive of tmp */
       proc_from = (taskid - it + np) % np;
-      msglen = 2 * (p_i2_start[nb][op][it + 1] - p_i2_start[nb][op][it]);
+      msglen = nffts*2*(p_i2_start[nb][op][it + 1] - p_i2_start[nb][op][it]);
       if (msglen > 0)
-        parall->adreceive(request_indx, msgtype, proc_from, msglen,
-                          tmp2 + 2 * p_i2_start[nb][op][it]);
+         parall->adreceive(request_indx, msgtype, proc_from, msglen,
+                           tmp2 + 2 * p_i2_start[nb][op][it]);
       // parall->adreceive(request_indx,msgtype,proc_from,msglen,&tmp2[2*p_i2_start[nb][op][it]]);
    }
-   for (it = 1; it < np; ++it) 
+   for (it=1; it<np; ++it) 
    {
       proc_to = (taskid + it) % np;
-      msglen = 2 * (p_i1_start[nb][op][it + 1] - p_i1_start[nb][op][it]);
+      msglen = nffts*2*(p_i1_start[nb][op][it + 1] - p_i1_start[nb][op][it]);
       if (msglen > 0)
-        parall->adsend(request_indx, msgtype, proc_to, msglen,
-                       tmp1 + 2 * p_i1_start[nb][op][it]);
+         parall->adsend(request_indx, msgtype, proc_to, msglen,
+                        tmp1 + 2 * p_i1_start[nb][op][it]);
       // parall->adsend(request_indx,msgtype,proc_to,msglen,&tmp1[2*p_i1_start[nb][op][it]]);
    }
 }
@@ -3935,8 +3945,11 @@ void d3db::c_ptranspose_ijk_end(const int nffts, const int nb, const int op, dou
    parall->awaitall(request_indx);
  
    /* unpack a array */
-   c_bindexcopy(n2, p_iq_to_i2[nb][op], tmp2, a);
-   c_bindexzero(n3, p_iz_to_i2[nb][op], a);
+   for (auto s=0; s<nffts; ++s)
+   {
+      c_bindexcopy_stride(nffts,n2, p_iq_to_i2[nb][op], tmp2 + 2*s, a + s*n2ft3d);
+      c_bindexzero(n3, p_iz_to_i2[nb][op], a + s*n2ft3d);
+   }
 }
 
 /********************************
@@ -4507,30 +4520,30 @@ void d3db::c_timereverse_start(const int nffts, double *a, double *tmp1, double
  
    indx = t_i1_start[0];
    nnfft3d = (t_i1_start[np] - t_i1_start[0] + 0);
-   c_aindexcopy(nnfft3d, t_iq_to_i1 + indx, a, tmp1 + indx);
+   for (auto s=0; s<nffts; ++s)
+      c_aindexcopy_stride(nffts,nnfft3d,t_iq_to_i1+indx, a + s*n2ft3d, tmp1 + indx + 2*s);
  
    /* it = 0, transpose data on same thread */
-   msglen = 2 * (t_i2_start[1] - t_i2_start[0]);
+   msglen = nffts*2*(t_i2_start[1] - t_i2_start[0]);
    // int one=1;
    // DCOPY_PWDFT(msglen,&(tmp1[2*t_i1_start[0]]),one,&(tmp2[2*t_i2_start[0]]),one);
-   std::memcpy(tmp2 + 2 * t_i2_start[0], tmp1 + 2 * t_i1_start[0],
-               msglen * sizeof(double));
+   std::memcpy(tmp2 + 2*t_i2_start[0], tmp1 + 2*t_i1_start[0], msglen*sizeof(double));
  
    /* receive packed array data */
-   for (it = 1; it < np; ++it) {
-     /* synchronous receive of tmp */
-     proc_from = (taskid - it + np) % np;
-     msglen = 2 * (t_i2_start[it + 1] - t_i2_start[it]);
-     if (msglen > 0)
-       parall->adreceive(request_indx, msgtype, proc_from, msglen,
-                         &tmp2[2 * t_i2_start[it]]);
+   for (it = 1; it < np; ++it) 
+   {
+      /* synchronous receive of tmp */
+      proc_from = (taskid - it + np) % np;
+      msglen = nffts*2*(t_i2_start[it + 1] - t_i2_start[it]);
+      if (msglen > 0)
+         parall->adreceive(request_indx, msgtype, proc_from, msglen, tmp2 + 2*t_i2_start[it]);
    }
-   for (it = 1; it < np; ++it) {
-     proc_to = (taskid + it) % np;
-     msglen = 2 * (t_i1_start[it + 1] - t_i1_start[it]);
-     if (msglen > 0)
-       parall->adsend(request_indx, msgtype, proc_to, msglen,
-                      &tmp1[2 * t_i1_start[it]]);
+   for (it = 1; it < np; ++it) 
+   {
+      proc_to = (taskid + it) % np;
+      msglen = nffts*2*(t_i1_start[it + 1] - t_i1_start[it]);
+      if (msglen > 0)
+         parall->adsend(request_indx, msgtype, proc_to, msglen, tmp1 + 2*t_i1_start[it]);
    }
 }
 
@@ -4566,7 +4579,8 @@ void d3db::c_timereverse_end(const int nffts, double *a, double *tmp1, double *t
    int nnfft3d = (t_i2_start[np] - t_i2_start[0] + 0);
 
    parall->awaitall(request_indx);
-   c_bindexcopy_conjg(nnfft3d, t_iq_to_i2 + indx, tmp2 + indx, a);
+   for (auto s=0; s<nffts; ++s)
+      c_bindexcopy_conjg_stride(nffts, nnfft3d, t_iq_to_i2+indx, tmp2 + indx + 2*s, a + s*n2ft3d);
 }
 
 
diff --git a/Nwpw/nwpwlib/device/gdevice2.cpp b/Nwpw/nwpwlib/device/gdevice2.cpp
index 78e9964c..0e1c248b 100644
--- a/Nwpw/nwpwlib/device/gdevice2.cpp
+++ b/Nwpw/nwpwlib/device/gdevice2.cpp
@@ -299,7 +299,7 @@ void gdevice2::batch_cffty_stages_tmpy(const int stage,const int tag,bool forwar
                                        double *a, double *tmpy, int da) {
 #if defined(NWPW_CUDA) || defined(NWPW_HIP)
    if (mygdevice2->hasgpu)
-      mygdevice2->batch_cffty_stages(stage,tag,forward, ny, nq, n2ft3d, a,da);
+      mygdevice2->batch_cffty_stages(stage,tag,forward,ny,nq,n2ft3d,a,da);
 #endif
 }
 
@@ -319,9 +319,9 @@ void gdevice2::batch_cfftz_tmpz(const int tag, bool forward, int nz, int nq, int
                               double *a, double *tmpz) {
 #if defined(NWPW_CUDA) || defined(NWPW_HIP)
    if (mygdevice2->hasgpu)
-      mygdevice2->batch_cfftz(tag, forward, nz, nq, n2ft3d, a);
+      mygdevice2->batch_cfftz(tag,forward,nz,nq,n2ft3d,a);
 #else
-   mygdevice2->batch_cfftz_tmpz(forward, nz, nq, n2ft3d, a, tmpz);
+   mygdevice2->batch_cfftz_tmpz(forward,nz,nq,n2ft3d,a,tmpz);
 #endif
 }
 
@@ -329,9 +329,9 @@ void gdevice2::batch_cfftz_tmpz_zero(const int tag, bool forward, int nz, int nq
                                      double *a, double *tmpz, bool *zero) {
 #if defined(NWPW_CUDA) || defined(NWPW_HIP)
    if (mygdevice2->hasgpu)
-      mygdevice2->batch_cfftz(tag,forward, nz, nq, n2ft3d, a);
+      mygdevice2->batch_cfftz(tag,forward,nz,nq,n2ft3d,a);
 #else
-   mygdevice2->batch_cfftz_tmpz_zero(forward, nz, nq, n2ft3d, a, tmpz, zero);
+   mygdevice2->batch_cfftz_tmpz_zero(forward,nz,nq,n2ft3d,a,tmpz,zero);
 #endif
 }
 
diff --git a/Nwpw/nwpwlib/lattice/Balance.cpp b/Nwpw/nwpwlib/lattice/Balance.cpp
index 815a72dd..2118dc6d 100644
--- a/Nwpw/nwpwlib/lattice/Balance.cpp
+++ b/Nwpw/nwpwlib/lattice/Balance.cpp
@@ -265,7 +265,7 @@ Balance::Balance(Parallel *inparall, const int maxsize0, const int *nidb, int *n
  */
 Balance::~Balance() 
 {
-   for (int nb = 0; nb < maxsize; ++nb) 
+   for (int nb=0; nb<maxsize; ++nb) 
    {
       delete[] proc_to_list[nb];
       delete[] proc_from_list[nb];
@@ -293,24 +293,22 @@ Balance::~Balance()
  */
 void Balance::c_unbalance_start(const int nffts, const int nb, double *a, const int request_indx, const int msgtype) 
 {
-   int j, pto, pfrom, msglen, indx;
- 
    if (sender_list[nb])
-     for (j = 0; j < npacket_list[nb]; ++j) 
+     for (auto j=0; j<npacket_list[nb]; ++j) 
      {
-        pfrom = proc_to_list[nb][j];
-        msglen = 2 * packet_size_list[nb][j];
-        indx = 2 * indx_start_list[nb][j];
+        int pfrom  = proc_to_list[nb][j];
+        int msglen = 2*packet_size_list[nb][j];
+        int indx   = 2*indx_start_list[nb][j];
         if (msglen > 0)
            parall->adreceive(request_indx, msgtype, pfrom, msglen, a + indx);
      }
  
    if (receiver_list[nb])
-     for (j = 0; j < npacket_list[nb]; ++j) 
+     for (auto j=0; j<npacket_list[nb]; ++j) 
      {
-        pto = proc_from_list[nb][j];
-        msglen = 2 * packet_size_list[nb][j];
-        indx = 2 * indx_start_list[nb][j];
+        int pto    = proc_from_list[nb][j];
+        int msglen = 2*packet_size_list[nb][j];
+        int indx   = 2*indx_start_list[nb][j];
         if (msglen > 0)
            parall->adsend(request_indx, msgtype, pto, msglen, a + indx);
      }
@@ -355,22 +353,24 @@ void Balance::c_unbalance(const int nb, double *a) {
   int j, pto, pfrom, msglen, indx;
 
   if (sender_list[nb])
-    for (j = 0; j < npacket_list[nb]; ++j) {
-      pfrom = proc_to_list[nb][j];
-      msglen = 2 * packet_size_list[nb][j];
-      indx = 2 * indx_start_list[nb][j];
-      if (msglen > 0)
-        parall->dreceive(1, 9, pfrom, msglen, &a[indx]);
-    }
+     for (auto j=0; j<npacket_list[nb]; ++j) 
+     {
+        int pfrom  = proc_to_list[nb][j];
+        int msglen = 2*packet_size_list[nb][j];
+        int indx   = 2*indx_start_list[nb][j];
+        if (msglen > 0)
+           parall->dreceive(1, 9, pfrom, msglen, a + indx);
+     }
 
   if (receiver_list[nb])
-    for (j = 0; j < npacket_list[nb]; ++j) {
-      pto = proc_from_list[nb][j];
-      msglen = 2 * packet_size_list[nb][j];
-      indx = 2 * indx_start_list[nb][j];
-      if (msglen > 0)
-        parall->dsend(1, 9, pto, msglen, &a[indx]);
-    }
+     for (auto j=0; j<npacket_list[nb]; ++j) 
+     {
+        int pto    = proc_from_list[nb][j];
+        int msglen = 2*packet_size_list[nb][j];
+        int indx   = 2*indx_start_list[nb][j];
+        if (msglen > 0)
+           parall->dsend(1, 9, pto, msglen, a + indx);
+     }
 }
 
 /********************************
diff --git a/Nwpw/nwpwlib/lattice/PGrid.cpp b/Nwpw/nwpwlib/lattice/PGrid.cpp
index d842a613..486c5218 100644
--- a/Nwpw/nwpwlib/lattice/PGrid.cpp
+++ b/Nwpw/nwpwlib/lattice/PGrid.cpp
@@ -1505,7 +1505,10 @@ void PGrid::c_unpack_start(const int nffts, const int nb, double *tmp1, double *
                            const int request_indx, const int msgtype) 
 {
    if (balanced)
-      mybalance->c_unbalance_start(nffts, nb, tmp1, request_indx, msgtype);
+   {
+      for (auto s=0; s<nffts; ++s)
+         mybalance->c_unbalance_start(nffts, nb, tmp1+s*n2ft3d, request_indx, msgtype);
+   }
 }
 
 /********************************
@@ -1519,10 +1522,11 @@ void PGrid::c_unpack_mid(const int nffts, const int nb, double *tmp1, double *tm
    if (balanced)
       mybalance->c_unbalance_end(nffts, nb, tmp1, request_indx);
  
-   std::memcpy(tmp2, tmp1, 2 * (nida[nb] + nidb2[nb]) * sizeof(double));
-   std::memset(tmp1, 0, n2ft3d * sizeof(double));
+   std::memcpy(tmp2, tmp1, nffts*2*(nida[nb]+nidb2[nb])*sizeof(double));
+   std::memset(tmp1, 0, nffts*n2ft3d*sizeof(double));
  
-   c_bindexcopy((nida[nb] + nidb2[nb]), packarray[nb], tmp2, tmp1);
+   for (auto s=0; s<nffts; ++s)
+      c_bindexcopy((nida[nb]+nidb2[nb]),packarray[nb], tmp2 + s*n2ft3d, tmp1 + s*n2ft3d);
    // c_bindexcopy(nida[nb]+nidb[nb],packarray[nb],tmp2,tmp1);
  
    d3db::c_timereverse_start(nffts, tmp1, zplane_tmp1, zplane_tmp2, request_indx, msgtype);
@@ -1560,26 +1564,32 @@ void PGrid::pfftbz(const int nffts, const int nb, double *tmp1, double *tmp2, in
        ***     A(kx,nz,ky) <- fft1d^(-1)[A(kx,kz,ky)]  ***
        ***************************************************/
       int indx0 = 0;
-      int indx2 = 0;
       int nn = 0;
-      for (auto q = 0; q < nq; ++q) {
-        for (auto i = 0; i < nxh; ++i) {
-          if (!zero_row3[nb][indx2]) {
-            auto kk = 0;
-            auto indx3 = 2 * i + indx0;
-            auto shift = 2 * nz * nn;
-            for (auto k = 0; k < nz; ++k) 
+      for (auto s=0; s<nffts; ++s)
+      {
+         int indx2 = 0;
+         for (auto q=0; q<nq; ++q) 
+         {
+            for (auto i=0; i<nxh; ++i) 
             {
-               tmp2[kk   + shift] = tmp1[indx3];
-               tmp2[kk+1 + shift] = tmp1[indx3 + 1];
-               kk += 2;
-               indx3 += nxh2;
+               if (!zero_row3[nb][indx2]) 
+               {
+                  auto kk = 0;
+                  auto indx3 = 2*i + indx0;
+                  auto shift = 2*nz*nn;
+                  for (auto k=0; k<nz; ++k) 
+                  {
+                     tmp2[kk   + shift] = tmp1[indx3];
+                     tmp2[kk+1 + shift] = tmp1[indx3 + 1];
+                     kk += 2;
+                     indx3 += nxh2;
+                  }
+                  nn += 1;
+               }
+               ++indx2;
             }
-            nn += 1;
-          }
-          ++indx2;
-        }
-        indx0 += nxhz2;
+            indx0 += nxhz2;
+         }
       }
      
       d3db::mygdevice.batch_cfftz_tmpz(d3db::fft_tag,false, nz, nn, n2ft3d, tmp2, d3db::tmpz);
@@ -1587,35 +1597,41 @@ void PGrid::pfftbz(const int nffts, const int nb, double *tmp1, double *tmp2, in
       //    dcfftb_(&nz,tmp2+2*nz*i,d3db::tmpz);
      
       indx0 = 0;
-      indx2 = 0;
       nn = 0;
-      for (auto q = 0; q < nq; ++q) {
-        for (auto i = 0; i < nxh; ++i) {
-          if (!zero_row3[nb][indx2]) {
-            auto kk = 0;
-            auto indx3 = 2 * i + indx0;
-            auto shift = 2 * nz * nn;
-            for (auto k = 0; k < nz; ++k) 
+      for (auto s=0; s<nffts; ++s)
+      {
+         int indx2 = 0;
+         for (auto q=0; q<nq; ++q) 
+         {
+            for (auto i=0; i<nxh; ++i) 
             {
-               tmp1[indx3]   = tmp2[kk   + shift];
-               tmp1[indx3+1] = tmp2[kk+1 + shift];
-               kk += 2;
-               indx3 += nxh2;
+               if (!zero_row3[nb][indx2]) 
+               {
+                  auto kk = 0;
+                  auto indx3 = 2 * i + indx0;
+                  auto shift = 2 * nz * nn;
+                  for (auto k = 0; k < nz; ++k) 
+                  {
+                     tmp1[indx3]   = tmp2[kk   + shift];
+                     tmp1[indx3+1] = tmp2[kk+1 + shift];
+                     kk += 2;
+                     indx3 += nxh2;
+                  }
+                  nn += 1;
+               }
+               ++indx2;
             }
-            nn += 1;
-          }
-          ++indx2;
-        }
-        indx0 += nxhz2;
+            indx0 += nxhz2;
+         }
       }
      
-     
       /***********************************************
        ***         Do a ptranspose of A            ***
        ***       A(kx,ky,nz) <- A(kx,nz,ky)        ***
        ************************************************/
       d3db::c_ptranspose1_jk_start(nffts, nb, tmp1, tmp2, tmp1, request_indx, 44);
    }
+
    /*************************
     **** hilbert mapping ****
     *************************/
@@ -1625,7 +1641,7 @@ void PGrid::pfftbz(const int nffts, const int nb, double *tmp1, double *tmp2, in
        ***     do fft along kz dimension            ***
        ***   A(nz,kx,ky) <- fft1d^(-1)[A(kz,kx,ky)] ***
        ************************************************/
-      d3db::mygdevice.batch_cfftz_tmpz_zero(d3db::fft_tag,false, nz, nq3, n2ft3d, tmp1, d3db::tmpz,
+      d3db::mygdevice.batch_cfftz_tmpz_zero(d3db::fft_tag,false, nz, nffts*nq3, n2ft3d, tmp1, d3db::tmpz,
                                     zero_row3[nb]);
      
       d3db::c_ptranspose_ijk_start(nffts, nb, 2, tmp1, tmp2, tmp1, request_indx, 45);
@@ -1645,7 +1661,7 @@ void PGrid::pfftby(const int nffts, const int nb, double *tmp1, double *tmp2, in
     **********************/
    if (maptype == 1) 
    {
-      auto nxh = nx / 2 + 1;
+      auto nxh = nx/2 + 1;
       auto nxh2 = nx + 2;
       auto nxhy2 = nxh2 * ny;
      
@@ -1660,27 +1676,32 @@ void PGrid::pfftby(const int nffts, const int nb, double *tmp1, double *tmp2, in
        ***   A(kx,ky,nz) <- fft1d[A(kx,ny,nz)]  ***
        ********************************************/
       int indx0 = 0;
-      int indx2 = 0;
       int nn = 0;
-      for (auto q=0; q<nq; ++q) {
-        for (auto i=0; i<nxh; ++i) {
-          if (!zero_row2[nb][indx2]) 
-          {
-             auto jj = 0;
-             auto indx3 = 2 * i + indx0;
-             auto shift = 2 * ny * nn;
-             for (auto j=0; j<ny; ++j) 
-             {
-                tmp1[jj   + shift] = tmp2[indx3];
-                tmp1[jj+1 + shift] = tmp2[indx3 + 1];
-                jj += 2;
-                indx3 += nxh2;
-             }
-             nn += 1;
-          }
-          ++indx2;
-        }
-        indx0 += nxhy2;
+      for (auto s=0; s<nffts; ++s)
+      {
+         int indx2 = 0;
+         for (auto q=0; q<nq; ++q) 
+         {
+            for (auto i=0; i<nxh; ++i) 
+            {
+               if (!zero_row2[nb][indx2]) 
+               {
+                  auto jj = 0;
+                  auto indx3 = 2 * i + indx0;
+                  auto shift = 2 * ny * nn;
+                  for (auto j=0; j<ny; ++j) 
+                  {
+                     tmp1[jj   + shift] = tmp2[indx3];
+                     tmp1[jj+1 + shift] = tmp2[indx3 + 1];
+                     jj += 2;
+                     indx3 += nxh2;
+                  }
+                  nn += 1;
+               }
+               ++indx2;
+            }
+            indx0 += nxhy2;
+         }
       }
      
       d3db::mygdevice.batch_cffty_tmpy(d3db::fft_tag,false, ny, nn, n2ft3d, tmp1, d3db::tmpy);
@@ -1688,29 +1709,35 @@ void PGrid::pfftby(const int nffts, const int nb, double *tmp1, double *tmp2, in
       //    dcfftb_(&ny,tmp1+2*ny*i,d3db::tmpy);
      
       indx0 = 0;
-      indx2 = 0;
       nn = 0;
-      for (auto q=0; q<nq; ++q) {
-        for (auto i=0; i<nxh; ++i) {
-          if (!zero_row2[nb][indx2]) 
-          {
-             auto jj = 0;
-             auto indx3 = 2 * i + indx0;
-             auto shift = 2 * ny * nn;
-             for (auto j = 0; j < ny; ++j) 
-             {
-                tmp2[indx3]   = tmp1[jj   + shift];
-                tmp2[indx3+1] = tmp1[jj+1 + shift];
-                jj += 2;
-                indx3 += nxh2;
-             }
-             nn += 1;
-          }
-          ++indx2;
-        }
-        indx0 += nxhy2;
+      for (auto s=0; s<nffts; ++s)
+      {
+         int indx2 = 0;
+         for (auto q=0; q<nq; ++q) 
+         {
+            for (auto i=0; i<nxh; ++i) 
+            {
+               if (!zero_row2[nb][indx2]) 
+               {
+                  auto jj = 0;
+                  auto indx3 = 2 * i + indx0;
+                  auto shift = 2 * ny * nn;
+                  for (auto j = 0; j < ny; ++j) 
+                  {
+                     tmp2[indx3]   = tmp1[jj   + shift];
+                     tmp2[indx3+1] = tmp1[jj+1 + shift];
+                     jj += 2;
+                     indx3 += nxh2;
+                  }
+                  nn += 1;
+               }
+               ++indx2;
+            }
+            indx0 += nxhy2;
+         }
       }
    }
+
    /*************************
     **** hilbert mapping ****
     *************************/
@@ -1722,7 +1749,7 @@ void PGrid::pfftby(const int nffts, const int nb, double *tmp1, double *tmp2, in
        ***     do fft along ny dimension        ***
        ***   A(ky,nz,kx) <- fft1d[A(ny,nz,kx)]  ***
        ********************************************/
-      d3db::mygdevice.batch_cffty_tmpy_zero(d3db::fft_tag,false,ny,nq2,n2ft3d,tmp2,d3db::tmpy,zero_row2[nb]);
+      d3db::mygdevice.batch_cffty_tmpy_zero(d3db::fft_tag,false,ny,nffts*nq2,n2ft3d,tmp2,d3db::tmpy,zero_row2[nb]);
      
       d3db::c_ptranspose_ijk_start(nffts, nb, 3, tmp2, tmp1, tmp2, request_indx, 46);
       // d3db::c_ptranspose_ijk(nb,3,tmp2,tmp1,tmp2);
@@ -1736,7 +1763,6 @@ void PGrid::pfftby(const int nffts, const int nb, double *tmp1, double *tmp2, in
  ********************************/
 void PGrid::pfftbx(const int nffts, const int nb, double *tmp1, double *tmp2, int request_indx) 
 {
- 
    /**********************
     **** slab mapping ****
     **********************/
@@ -1746,10 +1772,12 @@ void PGrid::pfftbx(const int nffts, const int nb, double *tmp1, double *tmp2, in
        ***     do fft along kx dimension            ***
        ***   A(nx,ny,nz) <- fft1d^(-1)[A(kx,ny,nz)] ***
        ************************************************/
-      d3db::mygdevice.batch_rfftx_tmpx(d3db::fft_tag,false, nx, ny * nq, n2ft3d, tmp2, d3db::tmpx);
-      d3db::zeroend_fftb(nx, ny, nq, 1, tmp2);
-      std::memcpy(tmp1, tmp2, n2ft3d * sizeof(double));
+      d3db::mygdevice.batch_rfftx_tmpx(d3db::fft_tag,false,nx,nffts*ny*nq,n2ft3d,tmp2,d3db::tmpx);
+
+      d3db::zeroend_fftb(nx,ny,nq,nffts,tmp2);
+      std::memcpy(tmp1,tmp2,nffts*n2ft3d*sizeof(double));
    }
+
    /*************************
     **** hilbert mapping ****
     *************************/
@@ -1761,10 +1789,11 @@ void PGrid::pfftbx(const int nffts, const int nb, double *tmp1, double *tmp2, in
        ***     do fft along kx dimension            ***
        ***   A(nx,ny,nz) <- fft1d^(-1)[A(kx,ny,nz)] ***
        ************************************************/
-      d3db::mygdevice.batch_rfftx_tmpx(d3db::fft_tag,false, nx, nq1, n2ft3d, tmp1, d3db::tmpx);
-      d3db::zeroend_fftb(nx, nq1, 1, 1, tmp1);
+      d3db::mygdevice.batch_rfftx_tmpx(d3db::fft_tag,false, nx, nffts*nq1, n2ft3d, tmp1, d3db::tmpx);
+      d3db::zeroend_fftb(nx, nq1, 1, nffts, tmp1);
       if (n2ft3d_map < n2ft3d)
-         std::memset(tmp1 + n2ft3d_map, 0, (n2ft3d - n2ft3d_map) * sizeof(double));
+         for (auto s=0; s<nffts; ++s)
+            std::memset(tmp1+s*n2ft3d + n2ft3d_map,0,(n2ft3d-n2ft3d_map)*sizeof(double));
    }
 }
 
@@ -1780,7 +1809,8 @@ void PGrid::pfftb_step(const int step, const int nffts, const int nb, double *a,
    {
       // parall->astart(request_indx,parall->np_i());
       // unpack start, tmp1-->tmp1
-      std::memcpy(tmp1, a, 2 * (nida[nb] + nidb[nb]) * sizeof(double));
+      for (auto s=0; s<nffts; ++s)
+         std::memcpy(tmp1 + s*n2ft3d, a + s*2*(nida[nb]+nidb[nb]), 2*(nida[nb]+nidb[nb])*sizeof(double));
       this->c_unpack_start(nffts, nb, tmp1, tmp2, request_indx, 47);
    } 
    else if (step == 1) 
@@ -1834,32 +1864,35 @@ void PGrid::pfftbz_start(const int nffts, const int nb, double *tmp1, double *tm
       ***     A(kx,nz,ky) <- fft1d^(-1)[A(kx,kz,ky)]  ***
       ***************************************************/
      int indx0 = 0;
-     int indx2 = 0;
      int nn = 0;
-     for (auto q=0; q<nq; ++q) 
+     for (auto s=0; s<nffts; ++s)
      {
-        for (auto i=0; i<nxh; ++i) 
+        int indx2 = 0;
+        for (auto q=0; q<nq; ++q) 
         {
-           if (!zero_row3[nb][indx2]) 
+           for (auto i=0; i<nxh; ++i) 
            {
-              auto kk = 0;
-              auto indx3 = 2 * i + indx0;
-              auto shift = 2 * nz * nn;
-              for (auto k=0; k<nz; ++k) 
+              if (!zero_row3[nb][indx2]) 
               {
-                 tmp2[kk   + shift] = tmp1[indx3];
-                 tmp2[kk+1 + shift] = tmp1[indx3+1];
-                 kk += 2;
-                 indx3 += nxh2;
+                 auto kk = 0;
+                 auto indx3 = 2 * i + indx0;
+                 auto shift = 2 * nz * nn;
+                 for (auto k=0; k<nz; ++k) 
+                 {
+                    tmp2[kk   + shift] = tmp1[indx3];
+                    tmp2[kk+1 + shift] = tmp1[indx3+1];
+                    kk += 2;
+                    indx3 += nxh2;
+                 }
+                 nn += 1;
               }
-              nn += 1;
+              ++indx2;
            }
-           ++indx2;
+           indx0 += nxhz2;
         }
-        indx0 += nxhz2;
      }
     
-     d3db::mygdevice.batch_cfftz_stages_tmpz(0,d3db::fft_tag,false, nz, nn, n2ft3d, tmp2, d3db::tmpz,da_indx);
+     d3db::mygdevice.batch_cfftz_stages_tmpz(0,d3db::fft_tag,false, nz, nffts*nn, n2ft3d, tmp2, d3db::tmpz,da_indx);
      // for (auto i=0; i<nn; ++i)
      //    dcfftb_(&nz,tmp2+2*nz*i,d3db::tmpz);
 
@@ -1873,7 +1906,7 @@ void PGrid::pfftbz_start(const int nffts, const int nb, double *tmp1, double *tm
       ***     do fft along kz dimension            ***
       ***   A(nz,kx,ky) <- fft1d^(-1)[A(kz,kx,ky)] ***
       ************************************************/
-     d3db::mygdevice.batch_cfftz_stages_tmpz_zero(0,d3db::fft_tag,false, nz, nq3, n2ft3d, tmp1, d3db::tmpz, zero_row3[nb], da_indx);
+     d3db::mygdevice.batch_cfftz_stages_tmpz_zero(0,d3db::fft_tag,false, nz, nffts*nq3, n2ft3d, tmp1, d3db::tmpz, zero_row3[nb], da_indx);
   }
 
 }
@@ -1899,17 +1932,20 @@ void PGrid::pfftbz_compute(const int nffts, const int nb, double *tmp1, double *
        ***     A(kx,nz,ky) <- fft1d^(-1)[A(kx,kz,ky)]  ***
        ***************************************************/
       int indx0 = 0;
-      int indx2 = 0;
       int nn = 0;
-      for (auto q=0; q<nq; ++q) 
+      for (auto s=0; s<nffts; ++s)
       {
-         for (auto i = 0; i < nxh; ++i) 
+         int indx2 = 0;
+         for (auto q=0; q<nq; ++q) 
          {
-            if (!zero_row3[nb][indx2]) 
+            for (auto i = 0; i < nxh; ++i) 
             {
-               nn += 1;
+               if (!zero_row3[nb][indx2]) 
+               {
+                  nn += 1;
+               }
+               ++indx2;
             }
-            ++indx2;
          }
       }
      
@@ -1927,7 +1963,7 @@ void PGrid::pfftbz_compute(const int nffts, const int nb, double *tmp1, double *
        ***     do fft along kz dimension            ***
        ***   A(nz,kx,ky) <- fft1d^(-1)[A(kz,kx,ky)] ***
        ************************************************/
-      d3db::mygdevice.batch_cfftz_stages_tmpz_zero(1,d3db::fft_tag,false, nz, nq3, n2ft3d, tmp1, d3db::tmpz,zero_row3[nb],da_indx);
+      d3db::mygdevice.batch_cfftz_stages_tmpz_zero(1,d3db::fft_tag,false, nz, nffts*nq3, n2ft3d, tmp1, d3db::tmpz,zero_row3[nb],da_indx);
    }
 }
 
@@ -1952,17 +1988,20 @@ void PGrid::pfftbz_end(const int nffts, const int nb, double *tmp1, double *tmp2
        ***     A(kx,nz,ky) <- fft1d^(-1)[A(kx,kz,ky)]  ***
        ***************************************************/
       int indx0 = 0;
-      int indx2 = 0;
       int nn = 0;
-      for (auto q=0; q<nq; ++q) 
+      for (auto s=0; s<nffts; ++s)
       {
-         for (auto i=0; i<nxh; ++i) 
+         int indx2 = 0;
+         for (auto q=0; q<nq; ++q) 
          {
-            if (!zero_row3[nb][indx2]) 
+            for (auto i=0; i<nxh; ++i) 
             {
-               nn += 1;
+               if (!zero_row3[nb][indx2]) 
+               {
+                  nn += 1;
+               }
+               ++indx2;
             }
-            ++indx2;
          }
       }
      
@@ -1971,34 +2010,39 @@ void PGrid::pfftbz_end(const int nffts, const int nb, double *tmp1, double *tmp2
       //    dcfftb_(&nz,tmp2+2*nz*i,d3db::tmpz);
      
       indx0 = 0;
-      indx2 = 0;
       nn = 0;
-      for (auto q=0; q<nq; ++q) {
-        for (auto i=0; i<nxh; ++i) {
-          if (!zero_row3[nb][indx2]) 
-          {
-             auto kk = 0;
-             auto indx3 = 2 * i + indx0;
-             auto shift = 2 * nz * nn;
-             for (auto k=0; k<nz; ++k) 
-             {
-                tmp1[indx3]   = tmp2[kk   + shift];
-                tmp1[indx3+1] = tmp2[kk+1 + shift];
-                kk += 2;
-                indx3 += nxh2;
-             }
-             nn += 1;
-          }
-          ++indx2;
-        }
-        indx0 += nxhz2;
+      for (auto s=0; s<nffts; ++s)
+      {
+         int indx2 = 0;
+         for (auto q=0; q<nq; ++q) 
+         {
+            for (auto i=0; i<nxh; ++i) 
+            {
+               if (!zero_row3[nb][indx2]) 
+               {
+                  auto kk = 0;
+                  auto indx3 = 2 * i + indx0;
+                  auto shift = 2 * nz * nn;
+                  for (auto k=0; k<nz; ++k) 
+                  {
+                     tmp1[indx3]   = tmp2[kk   + shift];
+                     tmp1[indx3+1] = tmp2[kk+1 + shift];
+                     kk += 2;
+                     indx3 += nxh2;
+                  }
+                  nn += 1;
+               }
+               ++indx2;
+            }
+            indx0 += nxhz2;
+         }
       }
      
       /***********************************************
        ***         Do a ptranspose of A            ***
        ***       A(kx,ky,nz) <- A(kx,nz,ky)        ***
        ************************************************/
-      d3db::c_ptranspose1_jk_start(nffts, nb, tmp1, tmp2, tmp1, request_indx, 44);
+      d3db::c_ptranspose1_jk_start(nffts,nb,tmp1,tmp2,tmp1,request_indx, 44);
    }
    /*************************
     **** hilbert mapping ****
@@ -2009,9 +2053,9 @@ void PGrid::pfftbz_end(const int nffts, const int nb, double *tmp1, double *tmp2
        ***     do fft along kz dimension            ***
        ***   A(nz,kx,ky) <- fft1d^(-1)[A(kz,kx,ky)] ***
        ************************************************/
-      d3db::mygdevice.batch_cfftz_stages_tmpz_zero(2,d3db::fft_tag,false, nz, nq3, n2ft3d, tmp1, d3db::tmpz, zero_row3[nb], da_indx);
+      d3db::mygdevice.batch_cfftz_stages_tmpz_zero(2,d3db::fft_tag,false, nz, nffts*nq3, n2ft3d, tmp1, d3db::tmpz, zero_row3[nb], da_indx);
      
-      d3db::c_ptranspose_ijk_start(nffts, nb, 2, tmp1, tmp2, tmp1, request_indx, 45);
+      d3db::c_ptranspose_ijk_start(nffts,nb,2,tmp1,tmp2,tmp1,request_indx, 45);
       // d3db::c_ptranspose_ijk(nb,2,tmp1,tmp2,tmp1);
    }
 }
@@ -2030,7 +2074,7 @@ void PGrid::pfftby_start(const int nffts, const int nb, double *tmp1, double *tm
     **********************/
    if (maptype == 1) 
    {
-      auto nxh = nx / 2 + 1;
+      auto nxh = nx/2 + 1;
       auto nxh2 = nx + 2;
       auto nxhy2 = nxh2 * ny;
      
@@ -2045,27 +2089,32 @@ void PGrid::pfftby_start(const int nffts, const int nb, double *tmp1, double *tm
        ***   A(kx,ky,nz) <- fft1d[A(kx,ny,nz)]  ***
        ********************************************/
       int indx0 = 0;
-      int indx2 = 0;
       int nn = 0;
-      for (auto q=0; q<nq; ++q) {
-        for (auto i=0; i<nxh; ++i) {
-          if (!zero_row2[nb][indx2]) 
-          {
-             auto jj = 0;
-             auto indx3 = 2 * i + indx0;
-             auto shift = 2 * ny * nn;
-             for (auto j = 0; j < ny; ++j) 
-             {
-                tmp1[jj   + shift] = tmp2[indx3];
-                tmp1[jj+1 + shift] = tmp2[indx3 + 1];
-                jj += 2;
-                indx3 += nxh2;
-             }
-             nn += 1;
-          }
-          ++indx2;
-        }
-        indx0 += nxhy2;
+      for (auto s=0; s<nffts; ++s)
+      {
+         int indx2 = 0;
+         for (auto q=0; q<nq; ++q) 
+         {
+            for (auto i=0; i<nxh; ++i) 
+            {
+               if (!zero_row2[nb][indx2]) 
+               {
+                  auto jj = 0;
+                  auto indx3 = 2 * i + indx0;
+                  auto shift = 2 * ny * nn;
+                  for (auto j = 0; j < ny; ++j) 
+                  {
+                     tmp1[jj   + shift] = tmp2[indx3];
+                     tmp1[jj+1 + shift] = tmp2[indx3 + 1];
+                     jj += 2;
+                     indx3 += nxh2;
+                  }
+                  nn += 1;
+               }
+               ++indx2;
+            }
+            indx0 += nxhy2;
+         }
       }
      
       d3db::mygdevice.batch_cffty_stages_tmpy(0,d3db::fft_tag,false,ny,nn,n2ft3d,tmp1,d3db::tmpy,da_indx);
@@ -2084,7 +2133,7 @@ void PGrid::pfftby_start(const int nffts, const int nb, double *tmp1, double *tm
        ***     do fft along ny dimension        ***
        ***   A(ky,nz,kx) <- fft1d[A(ny,nz,kx)]  ***
        ********************************************/
-      d3db::mygdevice.batch_cffty_stages_tmpy_zero(0,d3db::fft_tag,false,ny,nq2,n2ft3d,tmp2,d3db::tmpy,zero_row2[nb],da_indx);
+      d3db::mygdevice.batch_cffty_stages_tmpy_zero(0,d3db::fft_tag,false,ny,nffts*nq2,n2ft3d,tmp2,d3db::tmpy,zero_row2[nb],da_indx);
    }
 }
 
@@ -2101,7 +2150,7 @@ void PGrid::pfftby_compute(const int nffts, const int nb, double *tmp1, double *
     **********************/
    if (maptype == 1) 
    {
-      auto nxh = nx / 2 + 1;
+      auto nxh = nx/2 + 1;
       auto nxh2 = nx + 2;
       auto nxhy2 = nxh2 * ny;
      
@@ -2110,16 +2159,21 @@ void PGrid::pfftby_compute(const int nffts, const int nb, double *tmp1, double *
        ***   A(kx,ky,nz) <- fft1d[A(kx,ny,nz)]  ***
        ********************************************/
       int indx0 = 0;
-      int indx2 = 0;
       int nn = 0;
-      for (auto q=0; q<nq; ++q) {
-        for (auto i=0; i<nxh; ++i) {
-          if (!zero_row2[nb][indx2]) 
-          {
-             nn += 1;
-          }
-          ++indx2;
-        }
+      for (auto s=0; s<nffts; ++s)
+      {
+         int indx2 = 0;
+         for (auto q=0; q<nq; ++q) 
+         {
+            for (auto i=0; i<nxh; ++i) 
+            {
+               if (!zero_row2[nb][indx2]) 
+               {
+                  nn += 1;
+               }
+               ++indx2;
+            }
+         }
       }
      
       d3db::mygdevice.batch_cffty_stages_tmpy(1,d3db::fft_tag,false, ny, nn, n2ft3d, tmp1, d3db::tmpy,da_indx);
@@ -2135,7 +2189,7 @@ void PGrid::pfftby_compute(const int nffts, const int nb, double *tmp1, double *
        ***     do fft along ny dimension        ***
        ***   A(ky,nz,kx) <- fft1d[A(ny,nz,kx)]  ***
        ********************************************/
-      d3db::mygdevice.batch_cffty_stages_tmpy_zero(1,d3db::fft_tag,false,ny,nq2,n2ft3d,tmp2,d3db::tmpy,zero_row2[nb],da_indx);
+      d3db::mygdevice.batch_cffty_stages_tmpy_zero(1,d3db::fft_tag,false,ny,nffts*nq2,n2ft3d,tmp2,d3db::tmpy,zero_row2[nb],da_indx);
    }
 }
 
@@ -2161,16 +2215,21 @@ void PGrid::pfftby_end(const int nffts, const int nb, double *tmp1, double *tmp2
        ***   A(kx,ky,nz) <- fft1d[A(kx,ny,nz)]  ***
        ********************************************/
       int indx0 = 0;
-      int indx2 = 0;
       int nn = 0;
-      for (auto q=0; q<nq; ++q) {
-        for (auto i=0; i<nxh; ++i) {
-          if (!zero_row2[nb][indx2]) 
-          {
-             nn += 1;
-          }
-          ++indx2;
-        }
+      for (auto s=0; s<nffts; ++s)
+      {
+         int indx2 = 0;
+         for (auto q=0; q<nq; ++q) 
+         {
+            for (auto i=0; i<nxh; ++i) 
+            {
+               if (!zero_row2[nb][indx2]) 
+               {
+                  nn += 1;
+               }
+               ++indx2;
+            }
+         }
       }
      
       d3db::mygdevice.batch_cffty_stages_tmpy(2,d3db::fft_tag,false, ny, nn, n2ft3d, tmp1, d3db::tmpy,da_indx);
@@ -2178,27 +2237,32 @@ void PGrid::pfftby_end(const int nffts, const int nb, double *tmp1, double *tmp2
       //    dcfftb_(&ny,tmp1+2*ny*i,d3db::tmpy);
      
       indx0 = 0;
-      indx2 = 0;
       nn = 0;
-      for (auto q=0; q<nq; ++q) {
-        for (auto i=0; i<nxh; ++i) {
-          if (!zero_row2[nb][indx2]) 
-          {
-             auto jj = 0;
-             auto indx3 = 2 * i + indx0;
-             auto shift = 2 * ny * nn;
-             for (auto j = 0; j < ny; ++j) 
-             {
-                tmp2[indx3]   = tmp1[jj   + shift];
-                tmp2[indx3+1] = tmp1[jj+1 + shift];
-                jj += 2;
-                indx3 += nxh2;
-             }
-             nn += 1;
-          }
-          ++indx2;
-        }
-        indx0 += nxhy2;
+      for (auto s=0; s<nffts; ++s)
+      {
+         int indx2 = 0;
+         for (auto q=0; q<nq; ++q) 
+         {
+            for (auto i=0; i<nxh; ++i) 
+            {
+               if (!zero_row2[nb][indx2]) 
+               {
+                  auto jj = 0;
+                  auto indx3 = 2 * i + indx0;
+                  auto shift = 2 * ny * nn;
+                  for (auto j = 0; j < ny; ++j) 
+                  {
+                     tmp2[indx3]   = tmp1[jj   + shift];
+                     tmp2[indx3+1] = tmp1[jj+1 + shift];
+                     jj += 2;
+                     indx3 += nxh2;
+                  }
+                  nn += 1;
+               }
+               ++indx2;
+            }
+            indx0 += nxhy2;
+         }
       }
    }
    /*************************
@@ -2210,7 +2274,7 @@ void PGrid::pfftby_end(const int nffts, const int nb, double *tmp1, double *tmp2
       ***     do fft along ny dimension        ***
       ***   A(ky,nz,kx) <- fft1d[A(ny,nz,kx)]  ***
       ********************************************/
-      d3db::mygdevice.batch_cffty_stages_tmpy_zero(2,d3db::fft_tag,false,ny,nq2,n2ft3d,tmp2,d3db::tmpy,zero_row2[nb],da_indx);
+      d3db::mygdevice.batch_cffty_stages_tmpy_zero(2,d3db::fft_tag,false,ny,nffts*nq2,n2ft3d,tmp2,d3db::tmpy,zero_row2[nb],da_indx);
  
       d3db::c_ptranspose_ijk_start(nffts, nb, 3, tmp2, tmp1, tmp2, request_indx, 46);
       // d3db::c_ptranspose_ijk(nb,3,tmp2,tmp1,tmp2);
@@ -2236,7 +2300,7 @@ void PGrid::pfftbx_start(const int nffts, const int nb, double *tmp1, double *tm
        ***     do fft along kx dimension            ***
        ***   A(nx,ny,nz) <- fft1d^(-1)[A(kx,ny,nz)] ***
        ************************************************/
-      d3db::mygdevice.batch_rfftx_stages_tmpx(0,d3db::fft_tag,false, nx, ny * nq, n2ft3d, tmp2, d3db::tmpx,da_indx);
+      d3db::mygdevice.batch_rfftx_stages_tmpx(0,d3db::fft_tag,false, nx, nffts*ny*nq, n2ft3d, tmp2, d3db::tmpx,da_indx);
    }
    /*************************
     **** hilbert mapping ****
@@ -2249,7 +2313,7 @@ void PGrid::pfftbx_start(const int nffts, const int nb, double *tmp1, double *tm
        ***     do fft along kx dimension            ***
        ***   A(nx,ny,nz) <- fft1d^(-1)[A(kx,ny,nz)] ***
        ************************************************/
-      d3db::mygdevice.batch_rfftx_stages_tmpx(0,d3db::fft_tag,false, nx, nq1, n2ft3d, tmp1, d3db::tmpx,da_indx);
+      d3db::mygdevice.batch_rfftx_stages_tmpx(0,d3db::fft_tag,false, nx, nffts*nq1, n2ft3d, tmp1, d3db::tmpx,da_indx);
    }
 }
 
@@ -2270,7 +2334,7 @@ void PGrid::pfftbx_compute(const int nffts, const int nb, double *tmp1, double *
        ***     do fft along kx dimension            ***
        ***   A(nx,ny,nz) <- fft1d^(-1)[A(kx,ny,nz)] ***
        ************************************************/
-      d3db::mygdevice.batch_rfftx_stages_tmpx(1,d3db::fft_tag,false, nx, ny * nq, n2ft3d, tmp2, d3db::tmpx,da_indx);
+      d3db::mygdevice.batch_rfftx_stages_tmpx(1,d3db::fft_tag,false, nx, nffts*ny*nq, n2ft3d, tmp2, d3db::tmpx,da_indx);
    }
 
    /*************************
@@ -2282,7 +2346,7 @@ void PGrid::pfftbx_compute(const int nffts, const int nb, double *tmp1, double *
        ***     do fft along kx dimension            ***
        ***   A(nx,ny,nz) <- fft1d^(-1)[A(kx,ny,nz)] ***
        ************************************************/
-      d3db::mygdevice.batch_rfftx_stages_tmpx(1,d3db::fft_tag,false, nx, nq1, n2ft3d, tmp1, d3db::tmpx,da_indx);
+      d3db::mygdevice.batch_rfftx_stages_tmpx(1,d3db::fft_tag,false, nx, nffts*nq1, n2ft3d, tmp1, d3db::tmpx,da_indx);
    }
 }
 
@@ -2304,10 +2368,11 @@ void PGrid::pfftbx_end(const int nffts, const int nb, double *tmp1, double *tmp2
        ***     do fft along kx dimension            ***
        ***   A(nx,ny,nz) <- fft1d^(-1)[A(kx,ny,nz)] ***
        ************************************************/
-      d3db::mygdevice.batch_rfftx_stages_tmpx(2,d3db::fft_tag,false, nx, ny * nq, n2ft3d, tmp2, d3db::tmpx,da_indx);
-      d3db::zeroend_fftb(nx, ny, nq, 1, tmp2);
-      std::memcpy(tmp1, tmp2, n2ft3d * sizeof(double));
+      d3db::mygdevice.batch_rfftx_stages_tmpx(2,d3db::fft_tag,false, nx, nffts*ny*nq, n2ft3d, tmp2, d3db::tmpx,da_indx);
+      d3db::zeroend_fftb(nx, ny, nq, nffts, tmp2);
+      std::memcpy(tmp1,tmp2,nffts*n2ft3d*sizeof(double));
    }
+
    /*************************
     **** hilbert mapping ****
     *************************/
@@ -2317,10 +2382,11 @@ void PGrid::pfftbx_end(const int nffts, const int nb, double *tmp1, double *tmp2
        ***     do fft along kx dimension            ***
        ***   A(nx,ny,nz) <- fft1d^(-1)[A(kx,ny,nz)] ***
        ************************************************/
-      d3db::mygdevice.batch_rfftx_stages_tmpx(2,d3db::fft_tag,false, nx, nq1, n2ft3d, tmp1, d3db::tmpx,da_indx);
-      d3db::zeroend_fftb(nx, nq1, 1, 1, tmp1);
+      d3db::mygdevice.batch_rfftx_stages_tmpx(2,d3db::fft_tag,false, nx, nffts*nq1, n2ft3d, tmp1, d3db::tmpx,da_indx);
+      d3db::zeroend_fftb(nx, nq1, 1, nffts, tmp1);
       if (n2ft3d_map < n2ft3d)
-         std::memset(tmp1 + n2ft3d_map, 0, (n2ft3d - n2ft3d_map) * sizeof(double));
+         for (auto s=0; s<nffts; ++s)
+            std::memset(tmp1+s*n2ft3d + n2ft3d_map, 0, (n2ft3d-n2ft3d_map)*sizeof(double));
    }
 }
 
@@ -2339,7 +2405,8 @@ void PGrid::pfftb_step12(const int step, const int nffts, const int nb, double *
    {
       // parall->astart(request_indx,parall->np_i());
       // unpack start, tmp1-->tmp1
-      std::memcpy(tmp1, a, 2 * (nida[nb] + nidb[nb]) * sizeof(double));
+      for (auto s=0; s<nffts; ++s)
+         std::memcpy(tmp1+s*n2ft3d, a+s*n2ft3d, 2*(nida[nb]+nidb[nb])*sizeof(double));
       this->c_unpack_start(nffts, nb, tmp1, tmp2, request_indx, 47);
    } 
    else if (step == 1) 
@@ -2419,6 +2486,11 @@ void PGrid::cr_pfft3b_queuein(const int nb, const int nffts_in, double *a)
    //int nffts_in = 1;
    int shift1, shift2;
    int np = parall->np_i();
+   //std::cout << "cr queuein HERA" << std::endl;
+   //std::cout << "cr_queuein nb=" << nb << " nffts_in=" << nffts_in << std::endl;
+   //std::cout << "       aqsize=" << aqsize << " nffts_max=" << nffts_max << std::endl;
+   //std::cout << "       alast_index=" << alast_index << std::endl;
+   //std::cout << "src ptr:" << a << std::endl;
  
    for (auto q=0; q<aqsize; ++q) 
    {
@@ -2508,15 +2580,15 @@ void PGrid::pfftfx(const int nffts, const int nb, double *a, double *tmp1, doubl
    if (maptype == 1) 
    {
       // do fft along nx dimension
-      d3db::mygdevice.batch_rfftx_tmpx(d3db::fft_tag,true, nx, ny*nq, n2ft3d, a, d3db::tmpx);
-      std::memcpy(tmp1, a, n2ft3d * sizeof(double));
+      d3db::mygdevice.batch_rfftx_tmpx(d3db::fft_tag,true, nx, nffts*ny*nq, n2ft3d, a, d3db::tmpx);
+      std::memcpy(tmp1, a, nffts*n2ft3d*sizeof(double));
    }
    /**** hilbert mapping ****/
    else 
    {
       // do fft along nx dimension
       // A(kx,ny,nz) <- fft1d[A(nx,ny,nz)]
-      d3db::mygdevice.batch_rfftx_tmpx(d3db::fft_tag,true, nx, nq1, n2ft3d, a, d3db::tmpx);
+      d3db::mygdevice.batch_rfftx_tmpx(d3db::fft_tag,true, nx, nffts*nq1, n2ft3d, a, d3db::tmpx);
       d3db::c_ptranspose_ijk_start(nffts, nb, 0, a, tmp1, tmp2, request_indx, 40);
    }
 }
@@ -2526,65 +2598,103 @@ void PGrid::pfftfx(const int nffts, const int nb, double *a, double *tmp1, doubl
  *        PGrid::pfftfy         *
  *                              *
  ********************************/
+/**
+ * Performs the 3D Fast Fourier Transform (FFT) pipeline along the ny dimension of the dataset as part of
+ * a multi-dimensional FFT processing sequence. This function handles FFT computations differently
+ * based on the specified mapping strategy (slab or Hilbert mapping).
+ *
+ * @param nffts        Number of 3D FFTs to perform in this batch, affecting the scale of transformation.
+ * @param nb           Fermi sphere packing number, indicating energy cutoffs:
+ *                     - '1' for lower cutoff energy used in wavefunctions,
+ *                     - '0' for higher cutoff energy used in density grids.
+ * @param tmp1         Temporary storage array used as input data for the FFT or for data transposition.
+ * @param tmp2         Temporary storage array where the FFT output or transposed data is stored.
+ * @param request_indx Identifier for this specific FFT operation sequence, used to manage and track
+ *                     FFT operations across different stages and configurations.
+ *
+ * Processing Details:
+ * - Slab Mapping (maptype == 1):
+ *   - Transposes data for alignment along the ny dimension before performing FFT.
+ *   - Data is reorganized based on zero row filtering, followed by FFT along ny.
+ *   - The function uses optimized data handling to efficiently perform transformations and memory operations.
+ * - Hilbert Mapping (maptype != 1):
+ *   - Data is transposed from tmp1 to tmp2 for subsequent FFT along nz, using configurations that adapt
+ *     the process to Hilbert mapping strategies.
+ *
+ * Note:
+ * - 'maptype' determines the FFT configuration and influences how the function processes the data.
+ * - This function is critical in a high-performance computing context, designed to handle complex,
+ *   large-scale data arrays across multiple dimensions.
+ */
 void PGrid::pfftfy(const int nffts, const int nb, double *tmp1, double *tmp2, int request_indx) 
 {
    /**** slab mapping ****/
    if (maptype == 1) 
    {
-      auto nxh = nx / 2 + 1;
+      auto nxh = nx/2 + 1;
       auto nxh2 = nx + 2;
       auto nxhy2 = nxh2 * ny;
      
       // do fft along ny dimension
       // A(kx,ky,nz) <- fft1d[A(kx,ny,nz)]
       int indx0 = 0;
-      int indx2 = 0;
       int nn = 0;
-      for (auto q = 0; q < nq; ++q) {
-        for (auto i = 0; i < nxh; ++i) {
-          if (!zero_row2[nb][indx2]) 
-          {
-             auto jj = 0;
-             auto indx3 = 2 * i + indx0;
-             auto shift = 2 * ny * nn;
-             for (auto j = 0; j < ny; ++j) 
-             {
-                tmp2[jj   + shift] = tmp1[indx3];
-                tmp2[jj+1 + shift] = tmp1[indx3 + 1];
-                jj += 2;
-                indx3 += nxh2;
-             }
-             ++nn;
-          }
-          ++indx2;
-        }
-        indx0 += nxhy2;
+      for (auto s=0; s<nffts; ++s)
+      {
+         int indx2 = 0;
+         for (auto q=0; q<nq; ++q) 
+         {
+            for (auto i = 0; i < nxh; ++i) 
+            {
+               if (!zero_row2[nb][indx2]) 
+               {
+                  auto jj = 0;
+                  auto indx3 = 2 * i + indx0;
+                  auto shift = 2 * ny * nn;
+                  for (auto j = 0; j < ny; ++j) 
+                  {
+                     tmp2[jj   + shift] = tmp1[indx3];
+                     tmp2[jj+1 + shift] = tmp1[indx3 + 1];
+                     jj += 2;
+                     indx3 += nxh2;
+                  }
+                  ++nn;
+               }
+               ++indx2;
+            }
+            indx0 += nxhy2;
+         }
       }
      
       d3db::mygdevice.batch_cffty_tmpy(d3db::fft_tag,true, ny, nn, n2ft3d, tmp2, d3db::tmpy);
      
       indx0 = 0;
-      indx2 = 0;
       nn = 0;
-      for (auto q=0; q<nq; ++q) {
-        for (auto i=0; i<nxh; ++i) {
-          if (!zero_row2[nb][indx2]) 
-          {
-             auto jj = 0;
-             auto indx3 = 2 * i + indx0;
-             auto shift = 2 * ny * nn;
-             for (auto j = 0; j < ny; ++j) 
-             {
-                tmp1[indx3] = tmp2[jj + shift];
-                tmp1[indx3 + 1] = tmp2[jj + 1 + shift];
-                jj += 2;
-                indx3 += nxh2;
-             }
-             ++nn;
-          }
-          ++indx2;
-        }
-        indx0 += nxhy2;
+      for (auto s=0; s<nffts; ++s)
+      {
+         int indx2 = 0;
+         for (auto q=0; q<nq; ++q) 
+         {
+            for (auto i=0; i<nxh; ++i) 
+            {
+               if (!zero_row2[nb][indx2]) 
+               {
+                  auto jj = 0;
+                  auto indx3 = 2*i + indx0;
+                  auto shift = 2*ny*nn;
+                  for (auto j = 0; j < ny; ++j) 
+                  {
+                     tmp1[indx3] = tmp2[jj + shift];
+                     tmp1[indx3 + 1] = tmp2[jj + 1 + shift];
+                     jj += 2;
+                     indx3 += nxh2;
+                  }
+                  ++nn;
+               }
+               ++indx2;
+            }
+            indx0 += nxhy2;
+         }
       }
      
       // Do a transpose of A
@@ -2600,7 +2710,7 @@ void PGrid::pfftfy(const int nffts, const int nb, double *tmp1, double *tmp2, in
      
       // do fft along ny dimension
       // A(ky,nz,kx) <- fft1d[A(ny,nz,kx)]
-      d3db::mygdevice.batch_cffty_tmpy_zero(d3db::fft_tag,true,ny,nq2,n2ft3d,tmp1,d3db::tmpy,zero_row2[nb]);
+      d3db::mygdevice.batch_cffty_tmpy_zero(d3db::fft_tag,true,ny,nffts*nq2,n2ft3d,tmp1,d3db::tmpy,zero_row2[nb]);
      
       // in=tmp2, out=tmp2
       d3db::c_ptranspose_ijk_start(nffts, nb, 1, tmp1, tmp2, tmp1, request_indx, 42);
@@ -2612,6 +2722,34 @@ void PGrid::pfftfy(const int nffts, const int nb, double *tmp1, double *tmp2, in
  *        PGrid::pfftfz         *
  *                              *
  ********************************/
+/**
+ * Performs the 3D Fast Fourier Transform (FFT) pipeline along the nz dimension of the dataset as part of
+ * a multi-dimensional FFT processing sequence. This function handles FFT computations differently
+ * based on the specified mapping strategy (slab or Hilbert mapping).
+ *
+ * @param nffts        Number of 3D FFTs to perform in this batch, affecting the scale of transformation.
+ * @param nb           Fermi sphere packing number, indicating energy cutoffs:
+ *                     - '1' for lower cutoff energy used in wavefunctions,
+ *                     - '0' for higher cutoff energy used in density grids.
+ * @param tmp1         Temporary storage array used as input data for the FFT or for data transposition.
+ * @param tmp2         Temporary storage array where the FFT output or transposed data is stored.
+ * @param request_indx Identifier for this specific FFT operation sequence, used to manage and track
+ *                     FFT operations across different stages and configurations.
+ *
+ * Processing Details:
+ * - Slab Mapping (maptype == 1):
+ *   - Transposes data for alignment along the nz dimension before performing FFT.
+ *   - Data is reorganized based on zero row filtering, followed by FFT along nz.
+ *   - The function uses optimized data handling to efficiently perform transformations and memory operations.
+ * - Hilbert Mapping (maptype != 1):
+ *   - Data is transposed from tmp1 to tmp2 for subsequent FFT along nz, using configurations that adapt
+ *     the process to Hilbert mapping strategies.
+ *
+ * Note:
+ * - 'maptype' determines the FFT configuration and influences how the function processes the data.
+ * - This function is critical in a high-performance computing context, designed to handle complex,
+ *   large-scale data arrays across multiple dimensions.
+ */
 void PGrid::pfftfz(const int nffts, const int nb, double *tmp1, double *tmp2, int request_indx) 
 {
 
@@ -2627,52 +2765,61 @@ void PGrid::pfftfz(const int nffts, const int nb, double *tmp1, double *tmp2, in
       // do fft along nz dimension
       // A(kx,kz,ky) <- fft1d[A(kx,nz,ky)]
       int indx0 = 0;
-      int indx2 = 0;
       int nn = 0;
-      for (auto q=0; q<nq; ++q) {
-        for (auto i=0; i<nxh; ++i) {
-          if (!zero_row3[nb][indx2]) 
-          {
-             auto kk = 0;
-             auto indx3 = 2 * i + indx0;
-             auto shift = 2 * nz * nn;
-             for (auto k=0; k<nz; ++k) 
+      for (auto s=0; s<nffts; ++s)
+      {
+         int indx2 = 0;
+         for (auto q=0; q<nq; ++q) {
+           for (auto i=0; i<nxh; ++i) {
+             if (!zero_row3[nb][indx2]) 
              {
-                tmp1[kk   + shift] = tmp2[indx3];
-                tmp1[kk+1 + shift] = tmp2[indx3 + 1];
-                kk += 2;
-                indx3 += nxh2;
+                auto kk = 0;
+                auto indx3 = 2 * i + indx0;
+                auto shift = 2 * nz * nn;
+                for (auto k=0; k<nz; ++k) 
+                {
+                   tmp1[kk   + shift] = tmp2[indx3];
+                   tmp1[kk+1 + shift] = tmp2[indx3 + 1];
+                   kk += 2;
+                   indx3 += nxh2;
+                }
+                ++nn;
              }
-             ++nn;
-          }
-          ++indx2;
-        }
-        indx0 += nxhz2;
+             ++indx2;
+           }
+           indx0 += nxhz2;
+         }
       }
      
       d3db::mygdevice.batch_cfftz_tmpz(d3db::fft_tag,true, nz, nn, n2ft3d, tmp1, d3db::tmpz);
      
       indx0 = 0;
-      indx2 = 0;
       nn = 0;
-      for (auto q=0; q<nq; ++q) {
-        for (auto i=0; i<nxh; ++i) {
-          if (!zero_row3[nb][indx2]) 
-          {
-             auto kk = 0;
-             auto indx3 = 2 * i + indx0;
-             auto shift = 2 * nz * nn;
-             for (auto k = 0; k < nz; ++k) {
-               tmp2[indx3]   = tmp1[kk   + shift];
-               tmp2[indx3+1] = tmp1[kk+1 + shift];
-               kk += 2;
-               indx3 += nxh2;
-             }
-             ++nn;
-          }
-          ++indx2;
-        }
-        indx0 += nxhz2;
+      for (auto s=0; s<nffts; ++s)
+      {
+         int indx2 = 0;
+         for (auto q=0; q<nq; ++q) 
+         {
+            for (auto i=0; i<nxh; ++i) 
+            {
+               if (!zero_row3[nb][indx2]) 
+               {
+                  auto kk = 0;
+                  auto indx3 = 2*i + indx0;
+                  auto shift = 2*nz*nn;
+                  for (auto k=0; k<nz; ++k) 
+                  {
+                     tmp2[indx3]   = tmp1[kk   + shift];
+                     tmp2[indx3+1] = tmp1[kk+1 + shift];
+                     kk += 2;
+                     indx3 += nxh2;
+                  }
+                  ++nn;
+               }
+               ++indx2;
+            }
+            indx0 += nxhz2;
+         }
       }
    }
  
@@ -2684,7 +2831,7 @@ void PGrid::pfftfz(const int nffts, const int nb, double *tmp1, double *tmp2, in
      
       // do fft along nz dimension
       // A(kz,kx,ky) <- fft1d[A(nz,kx,ky)]
-      d3db::mygdevice.batch_cfftz_tmpz_zero(d3db::fft_tag,true, nz, nq3, n2ft3d, tmp2, d3db::tmpz,zero_row3[nb]);
+      d3db::mygdevice.batch_cfftz_tmpz_zero(d3db::fft_tag,true,nz,nffts*nq3,n2ft3d,tmp2,d3db::tmpz,zero_row3[nb]);
    }
 }
 
@@ -2693,6 +2840,38 @@ void PGrid::pfftfz(const int nffts, const int nb, double *tmp1, double *tmp2, in
  *       PGrid::pfftf_step      *
  *                              *
  ********************************/
+/**
+ * Controls the execution of multi-dimensional Fast Fourier Transform (FFT) pipeline operations 
+ * on a dataset by orchestrating steps across different dimensions (x, y, z). This function
+ * sequentially triggers FFT operations or packing steps based on the 'step' parameter.
+ *
+ * @param step         The specific stage of the FFT process to execute, dictating which FFT
+ *                     or packing function is called:
+ *                     - 0: Initial transfer of data for FFT along x-axis.
+ *                     - 1: Perform FFT along y-axis.
+ *                     - 2: Finalize FFT along z-axis and start packing data.
+ *                     - 3: Complete packing of data.
+ * @param nffts        Number of 3D FFTs to perform, affects the scope of transformation.
+ * @param nb           Fermi sphere packing number, influencing configurations like:
+ *                     - '1' for lower energy cutoff used in wavefunctions.
+ *                     - '0' for higher energy cutoff used in density grids.
+ * @param a            Input data array for the initial FFT step.
+ * @param tmp1         Temporary storage array used as an intermediary in FFT transformations.
+ * @param tmp2         Temporary storage array used to hold the output of FFT operations.
+ * @param request_indx Identifier for this specific FFT operation sequence, used to track
+ *                     and manage the request through various transformation stages.
+ *
+ * Processing Flow:
+ * - step 0: Moves memory to device and prepares for FFT along the x-dimension.
+ * - step 1: Processes FFT along the y-dimension.
+ * - step 2: Completes FFT along the z-dimension and begins data packing.
+ * - step 3: Finalizes data packing and prepares it for further processing or output.
+ *
+ * Note:
+ * - Each step corresponds to a specific operation and must be executed in sequence for correct
+ *   processing. Misalignment in step ordering can lead to incorrect results or data corruption.
+ * - Designed for high-performance computing contexts dealing with large-scale data sets in multi-dimensional spaces.
+ */
 void PGrid::pfftf_step(const int step, const int nffts,  const int nb, double *a, double *tmp1, double *tmp2, int request_indx)
 {
    if (step==0)
@@ -2723,22 +2902,49 @@ void PGrid::pfftf_step(const int step, const int nffts,  const int nb, double *a
  *    PGrid::pfftfx_start       *
  *                              *
  ********************************/
+/**
+ * Initiates batch processing of 3D Fast Fourier Transforms (FFT) pipeline  on a given dataset, setting up data and
+ * performing the initial stages of FFT. This function supports two types of FFT mappings, each affecting
+ * the execution along the nx dimension.
+ *
+ * @param nffts        Number of 3D FFTs to perform in this batch, affecting the transformation scale.
+ * @param nb           Fermi sphere packing number, with possible values:
+ *                     - '1' for lower cutoff energy used in wavefunctions,
+ *                     - '0' for higher cutoff energy used in density grids.
+ * @param a            Pointer to the input data array.
+ * @param tmp1         Pointer to a temporary storage array (unused in this function).
+ * @param tmp2         Pointer to a temporary storage array, used for storing data before FFT processing.
+ * @param request_indx Index or identifier for this specific FFT request (unused directly in this function).
+ * @param da_indx      Dataset identifier used in GPU operations to specify a subset of data.
+ *
+ * Processing Details:
+ * - Slab Mapping (maptype == 1):
+ *   - Copies data from 'a' to 'tmp2'.
+ *   - Performs FFT along the nx dimension using the copied data.
+ * - Hilbert Mapping (maptype != 1):
+ *   - Copies data from 'a' to 'tmp2'.
+ *   - Executes FFT along the nx dimension, tailored for the Hilbert mapping strategy.
+ *
+ * Note:
+ * - 'maptype' needs to be set externally and determines the FFT configuration.
+ * - Designed for high-performance computing contexts for processing extensive multidimensional datasets.
+ */
 void PGrid::pfftfx_start(const int nffts, const int nb, double *a, double *tmp1, double *tmp2, int request_indx, int da_indx)
 {
    /**** slab mapping ****/
    if (maptype == 1) 
    {
       // do fft along nx dimension
-      std::memcpy(tmp2, a, n2ft3d * sizeof(double));
-      d3db::mygdevice.batch_rfftx_stages_tmpx(0,d3db::fft_tag,true, nx, ny*nq, n2ft3d, tmp2, d3db::tmpx,da_indx);
+      std::memcpy(tmp2, a, nffts*n2ft3d*sizeof(double));
+      d3db::mygdevice.batch_rfftx_stages_tmpx(0,d3db::fft_tag,true,nx,nffts*ny*nq, n2ft3d,tmp2,d3db::tmpx,da_indx);
    }
    /**** hilbert mapping ****/
    else 
    {
       // do fft along nx dimension
       // A(kx,ny,nz) <- fft1d[A(nx,ny,nz)]
-      std::memcpy(tmp2, a, n2ft3d * sizeof(double));
-      d3db::mygdevice.batch_rfftx_stages_tmpx(0,d3db::fft_tag,true, nx, nq1, n2ft3d, tmp2, d3db::tmpx,da_indx);
+      std::memcpy(tmp2, a, nffts*n2ft3d*sizeof(double));
+      d3db::mygdevice.batch_rfftx_stages_tmpx(0,d3db::fft_tag,true,nx,nffts*nq1,n2ft3d,tmp2,d3db::tmpx,da_indx);
    }
 }
 
@@ -2747,103 +2953,184 @@ void PGrid::pfftfx_start(const int nffts, const int nb, double *a, double *tmp1,
  *    PGrid::pfftfx_compute     *
  *                              *
  ********************************/
+/**
+ * Executes the main computation step of the 3D Fast Fourier Transforms (FFT) pipeline  for a given dataset,
+ * which is part of a batch processing sequence initiated by `pfftfx_start`. This function specifically
+ * handles the computation stages of FFT based on the selected mapping strategy.
+ *
+ * @param nffts        Number of 3D FFTs to perform in this batch, determines the scale of transformation.
+ * @param nb           Fermi sphere packing number, indicating energy cutoffs:
+ *                     - '1' for lower cutoff energy used in wavefunctions,
+ *                     - '0' for higher cutoff energy used in density grids.
+ * @param a            Pointer to the input data array (not directly used in this function).
+ * @param tmp1         Pointer to a temporary storage array (unused in this function).
+ * @param tmp2         Pointer to the temporary storage array, used for holding data during FFT processing.
+ * @param request_indx Index or identifier for this specific FFT request (not directly used in this function).
+ * @param da_indx      Dataset identifier used in GPU operations to specify the data subset being processed.
+ *
+ * Processing Details:
+ * - Slab Mapping (maptype == 1):
+ *   - Performs FFT computation along the nx dimension using the data in 'tmp2'.
+ * - Hilbert Mapping (maptype != 1):
+ *   - Executes FFT computation along the nx dimension, tailored for the Hilbert mapping strategy.
+ *
+ * Note:
+ * - 'maptype' is externally defined and selects the FFT configuration.
+ * - This function is crucial within a high-performance computing framework intended for extensive multidimensional data arrays.
+ */
 void PGrid::pfftfx_compute(const int nffts, const int nb, double *a, double *tmp1, double *tmp2, int request_indx, int da_indx)
 {
    /**** slab mapping ****/
    if (maptype == 1) 
    {
       // do fft along nx dimension
-      d3db::mygdevice.batch_rfftx_stages_tmpx(1,d3db::fft_tag,true, nx, ny*nq, n2ft3d, tmp2, d3db::tmpx,da_indx);
+      d3db::mygdevice.batch_rfftx_stages_tmpx(1,d3db::fft_tag,true,nx,nffts*ny*nq,n2ft3d,tmp2,d3db::tmpx,da_indx);
    }
    /**** hilbert mapping ****/
    else 
    {
       // do fft along nx dimension
       // A(kx,ny,nz) <- fft1d[A(nx,ny,nz)]
-      d3db::mygdevice.batch_rfftx_stages_tmpx(1,d3db::fft_tag,true, nx, nq1, n2ft3d, tmp2, d3db::tmpx,da_indx);
+      d3db::mygdevice.batch_rfftx_stages_tmpx(1,d3db::fft_tag,true,nx,nffts*nq1,n2ft3d,tmp2,d3db::tmpx,da_indx);
    }
 }
 
-
 /********************************
  *                              *
  *      PGrid::pfftfx_end       *
  *                              *
  ********************************/
+/**
+ * Finalizes the batch of 3D Fast Fourier Transform (FFT) pipeline operations initiated by `pfftfx_start`.
+ * This function performs the last stages of FFT and handles the necessary post-processing for the data.
+ * It supports different processing strategies based on the 'maptype' setting.
+ *
+ * @param nffts        Number of 3D FFTs to be processed in this batch.
+ * @param nb           Fermi sphere packing number, indicating energy cutoffs:
+ *                     - '1' for lower cutoff energy used in wavefunctions.
+ *                     - '0' for higher cutoff energy used in density grids.
+ * @param a            Pointer to the initial data array (unused directly in this function).
+ * @param tmp1         Pointer to the temporary array where the final processed data is stored.
+ * @param tmp2         Pointer to the temporary array for intermediate data processing.
+ * @param request_indx Identifier for this FFT request, used for data transposition in Hilbert mapping.
+ * @param da_indx      Dataset identifier used in GPU operations.
+ *
+ * Processing Details:
+ * - Slab Mapping (maptype == 1):
+ *   - Completes FFT along the nx dimension.
+ *   - Copies processed data from 'tmp2' to 'tmp1'.
+ * - Hilbert Mapping (maptype != 1):
+ *   - Completes FFT along the nx dimension.
+ *   - Invokes 'c_ptranspose_ijk_start' for complex data transposition.
+ *
+ * Note:
+ * - Assumes external setting of 'maptype' to distinguish between Slab and Hilbert mappings.
+ * - Integral to a high-performance computing framework for multidimensional FFT operations.
+ */
 void PGrid::pfftfx_end(const int nffts, const int nb, double *a, double *tmp1, double *tmp2, int request_indx, int da_indx)
 {
    /**** slab mapping ****/
    if (maptype == 1) 
    {
       // do fft along nx dimension
-      d3db::mygdevice.batch_rfftx_stages_tmpx(2,d3db::fft_tag,true, nx, ny*nq, n2ft3d, tmp2, d3db::tmpx,da_indx);
-      std::memcpy(tmp1, tmp2, n2ft3d * sizeof(double));
+      d3db::mygdevice.batch_rfftx_stages_tmpx(2,d3db::fft_tag,true,nx,nffts*ny*nq,n2ft3d,tmp2,d3db::tmpx,da_indx);
+      std::memcpy(tmp1,tmp2,nffts*n2ft3d*sizeof(double));
    }
    /**** hilbert mapping ****/
    else 
    {
       // do fft along nx dimension
       // A(kx,ny,nz) <- fft1d[A(nx,ny,nz)]
-      d3db::mygdevice.batch_rfftx_stages_tmpx(2,d3db::fft_tag,true, nx, nq1, n2ft3d, tmp2, d3db::tmpx,da_indx);
+      d3db::mygdevice.batch_rfftx_stages_tmpx(2,d3db::fft_tag,true, nx, nffts*nq1, n2ft3d, tmp2, d3db::tmpx,da_indx);
  
-      d3db::c_ptranspose_ijk_start(nffts, nb, 0, tmp2, tmp1, tmp2, request_indx, 40);
+      d3db::c_ptranspose_ijk_start(nffts,nb,0,tmp2,tmp1,tmp2,request_indx,40);
  
    }
 }
 
-
-
 /********************************
  *                              *
  *      PGrid::pfftfy_start     *
  *                              *
  ********************************/
+/**
+ * Initiates the batch processing of 3D Fast Fourier Transforms (FFT) pipeline  along the ny dimension,
+ * catering to different FFT mapping strategies: slab and Hilbert. This function sets up and
+ * executes initial stages of FFT tailored for large-scale data arrays.
+ *
+ * @param nffts        Number of 3D FFTs to perform in this batch, affecting the transformation scale.
+ * @param nb           Fermi sphere packing number, with possible values:
+ *                     - '1' for lower cutoff energy used in wavefunctions,
+ *                     - '0' for higher cutoff energy used in density grids.
+ * @param tmp1         Pointer to the input array of data.
+ * @param tmp2         Pointer to the output array where transformed data is stored.
+ * @param request_indx Index or identifier for this specific FFT request, used in transposition in Hilbert mapping.
+ * @param da_indx      Dataset identifier used in GPU operations to specify a subset of data.
+ *
+ * Processing Details:
+ * - Slab Mapping (maptype == 1):
+ *   - Arranges and copies data based on zero row filtering before performing FFT along ny.
+ *   - Data preparation involves optimized memory handling to set up for batched FFT operations.
+ * - Hilbert Mapping (maptype != 1):
+ *   - Data transposition is performed from tmp1 to tmp2 using pre-defined configurations.
+ *   - Post-transposition, executes FFT along ny dimension, designed for the Hilbert mapping approach.
+ *
+ * Note:
+ * - 'maptype' is externally defined, guiding the FFT configuration and execution approach.
+ * - This function is integral to a high-performance computing framework focused on multidimensional FFT operations.
+ */
 void PGrid::pfftfy_start(const int nffts, const int nb, double *tmp1, double *tmp2, int request_indx, int da_indx) 
 {
    /**** slab mapping ****/
    if (maptype == 1) 
    {
-      auto nxh = nx / 2 + 1;
+      auto nxh = nx/2 + 1;
       auto nxh2 = nx + 2;
       auto nxhy2 = nxh2 * ny;
      
       // do fft along ny dimension
       // A(kx,ky,nz) <- fft1d[A(kx,ny,nz)]
       int indx0 = 0;
-      int indx2 = 0;
       int nn = 0;
-      for (auto q=0; q<nq; ++q) {
-        for (auto i=0; i<nxh; ++i) {
-          if (!zero_row2[nb][indx2]) {
-            auto jj = 0;
-            auto indx3 = 2 * i + indx0;
-            auto shift = 2 * ny * nn;
-            for (auto j=0; j<ny; ++j) 
+      for (auto s=0; s<nffts; ++s)
+      {
+         int indx2 = 0;
+         for (auto q=0; q<(nffts*nq); ++q) 
+         {
+            for (auto i=0; i<nxh; ++i) 
             {
-               tmp2[jj   + shift] = tmp1[indx3];
-               tmp2[jj+1 + shift] = tmp1[indx3 + 1];
-               jj += 2;
-               indx3 += nxh2;
+               if (!zero_row2[nb][indx2]) 
+               {
+                  auto jj = 0;
+                  auto indx3 = 2*i + indx0;
+                  auto shift = 2*ny*nn;
+                  for (auto j=0; j<ny; ++j) 
+                  {
+                     tmp2[jj   + shift] = tmp1[indx3];
+                     tmp2[jj+1 + shift] = tmp1[indx3 + 1];
+                     jj += 2;
+                     indx3 += nxh2;
+                  }
+                  ++nn;
+               }
+               ++indx2;
             }
-            ++nn;
-          }
-          ++indx2;
-        }
-        indx0 += nxhy2;
+            indx0 += nxhy2;
+         }
       }
      
-      d3db::mygdevice.batch_cffty_stages_tmpy(0,d3db::fft_tag,true, ny, nn, n2ft3d, tmp2, d3db::tmpy,da_indx);
+      d3db::mygdevice.batch_cffty_stages_tmpy(0,d3db::fft_tag,true,ny,nn,n2ft3d,tmp2,d3db::tmpy,da_indx);
    }
  
    /**** hilbert mapping ****/
    else 
    {
       // in=tmp1, out=tmp2
-      d3db::c_ptranspose_ijk_end(nffts, nb, 0, tmp1, tmp2, request_indx);
+      d3db::c_ptranspose_ijk_end(nffts,nb,0,tmp1,tmp2,request_indx);
      
       // do fft along ny dimension
       // A(ky,nz,kx) <- fft1d[A(ny,nz,kx)]
-      d3db::mygdevice.batch_cffty_stages_tmpy_zero(0,d3db::fft_tag,true,ny,nq2,n2ft3d,tmp1,d3db::tmpy,zero_row2[nb],da_indx);
+      d3db::mygdevice.batch_cffty_stages_tmpy_zero(0,d3db::fft_tag,true,ny,nffts*nq2,n2ft3d,tmp1,d3db::tmpy,zero_row2[nb],da_indx);
    }
 }
 
@@ -2852,6 +3139,31 @@ void PGrid::pfftfy_start(const int nffts, const int nb, double *tmp1, double *tm
  *      PGrid::pfftfy_compute   *
  *                              *
  ********************************/
+/**
+ * Executes the main computation step of the 3D Fast Fourier Transforms (FFT) pipeline along the ny dimension
+ * for a given dataset. This function is part of a batch processing sequence that deals with large-scale
+ * data arrays, applying FFT based on selected mapping strategies.
+ *
+ * @param nffts        Number of 3D FFTs to perform in this batch, determines the scale of transformation.
+ * @param nb           Fermi sphere packing number, with possible values:
+ *                     - '1' for lower cutoff energy used in wavefunctions,
+ *                     - '0' for higher cutoff energy used in density grids.
+ * @param tmp1         Pointer to the input array of data used in Hilbert mapping.
+ * @param tmp2         Pointer to the output array where transformed data is stored.
+ * @param request_indx Index or identifier for this specific FFT request, used in configurations for Hilbert mapping.
+ * @param da_indx      Dataset identifier used in GPU operations to specify the subset of data.
+ *
+ * Processing Details:
+ * - Slab Mapping (maptype == 1):
+ *   - Arranges data based on zero row filtering before performing FFT along ny.
+ *   - Efficiently handles memory and data flow to optimize FFT computations.
+ * - Hilbert Mapping (maptype != 1):
+ *   - Directly uses pre-configured data from tmp1 for FFT computations along ny, adapted for Hilbert mapping strategy.
+ *
+ * Note:
+ * - 'maptype' must be externally set and determines the FFT configuration and execution approach.
+ * - This function is essential in high-performance computing frameworks focused on extensive multidimensional FFT operations.
+ */
 void PGrid::pfftfy_compute(const int nffts, const int nb, double *tmp1, double *tmp2, int request_indx, int da_indx) 
 {
    /**** slab mapping ****/
@@ -2864,19 +3176,24 @@ void PGrid::pfftfy_compute(const int nffts, const int nb, double *tmp1, double *
       // do fft along ny dimension
       // A(kx,ky,nz) <- fft1d[A(kx,ny,nz)]
       int indx0 = 0;
-      int indx2 = 0;
       int nn = 0;
-      for (auto q=0; q<nq; ++q) {
-        for (auto i=0; i<nxh; ++i) {
-          if (!zero_row2[nb][indx2]) 
-          {
-             ++nn;
-          }
-          ++indx2;
-        }
+      for (auto s=0; s<nffts; ++s)
+      {
+         int indx2 = 0;
+         for (auto q=0; q<nq; ++q)
+         {
+            for (auto i=0; i<nxh; ++i) 
+            {
+               if (!zero_row2[nb][indx2]) 
+               {
+                  ++nn;
+               }
+               ++indx2;
+            }
+         }
       }
      
-      d3db::mygdevice.batch_cffty_stages_tmpy(1,d3db::fft_tag,true, ny, nn, n2ft3d, tmp2, d3db::tmpy,da_indx);
+      d3db::mygdevice.batch_cffty_stages_tmpy(1,d3db::fft_tag,true,ny,nn,n2ft3d,tmp2,d3db::tmpy,da_indx);
    }
  
    /**** hilbert mapping ****/
@@ -2885,7 +3202,7 @@ void PGrid::pfftfy_compute(const int nffts, const int nb, double *tmp1, double *
       // in=tmp1, out=tmp2
       // do fft along ny dimension
       // A(ky,nz,kx) <- fft1d[A(ny,nz,kx)]
-      d3db::mygdevice.batch_cffty_stages_tmpy_zero(1,d3db::fft_tag,true,ny,nq2,n2ft3d,tmp1,d3db::tmpy,zero_row2[nb],da_indx);
+      d3db::mygdevice.batch_cffty_stages_tmpy_zero(1,d3db::fft_tag,true,ny,nffts*nq2,n2ft3d,tmp1,d3db::tmpy,zero_row2[nb],da_indx);
    }
 }
 
@@ -2909,54 +3226,64 @@ void PGrid::pfftfy_end(const int nffts, const int nb, double *tmp1, double *tmp2
       // do fft along ny dimension
       // A(kx,ky,nz) <- fft1d[A(kx,ny,nz)]
       int indx0 = 0;
-      int indx2 = 0;
       int nn = 0;
-      for (auto q=0; q<nq; ++q) {
-        for (auto i=0; i<nxh; ++i) {
-          if (!zero_row2[nb][indx2]) 
-          {
-             ++nn;
-          }
-          ++indx2;
-        }
+      for (auto s=0; s<nffts; ++s)
+      {
+         int indx2 = 0;
+         for (auto q=0; q<nq; ++q) 
+         {
+            for (auto i=0; i<nxh; ++i) 
+            {
+               if (!zero_row2[nb][indx2]) 
+               {
+                  ++nn;
+               }
+               ++indx2;
+            }
+         }
       }
      
-      d3db::mygdevice.batch_cffty_stages_tmpy(2,d3db::fft_tag,true, ny, nn, n2ft3d, tmp2, d3db::tmpy,da_indx);
+      d3db::mygdevice.batch_cffty_stages_tmpy(2,d3db::fft_tag,true,ny,nn,n2ft3d,tmp2,d3db::tmpy,da_indx);
      
       indx0 = 0;
-      indx2 = 0;
       nn = 0;
-      for (auto q=0; q<nq; ++q) {
-        for (auto i=0; i<nxh; ++i) {
-          if (!zero_row2[nb][indx2]) 
-          {
-             auto jj = 0;
-             auto indx3 = 2 * i + indx0;
-             auto shift = 2 * ny * nn;
-             for (auto j=0; j<ny; ++j) 
-             {
-               tmp1[indx3]   = tmp2[jj   + shift];
-               tmp1[indx3+1] = tmp2[jj+1 + shift];
-               jj += 2;
-               indx3 += nxh2;
-             }
-             ++nn;
-          }
-          ++indx2;
-        }
-        indx0 += nxhy2;
+      for (auto s=0; s<nffts; ++s)
+      {
+         int indx2 = 0;
+         for (auto q=0; q<nq; ++q) 
+         {
+            for (auto i=0; i<nxh; ++i) 
+            {
+               if (!zero_row2[nb][indx2]) 
+               {
+                  auto jj = 0;
+                  auto indx3 = 2 * i + indx0;
+                  auto shift = 2 * ny * nn;
+                  for (auto j=0; j<ny; ++j) 
+                  {
+                     tmp1[indx3]   = tmp2[jj   + shift];
+                     tmp1[indx3+1] = tmp2[jj+1 + shift];
+                     jj += 2;
+                     indx3 += nxh2;
+                  }
+                  ++nn;
+               }
+               ++indx2;
+            }
+            indx0 += nxhy2;
+         }
       }
      
       // Do a transpose of A
       // A(ky,nz,ky) <- A(kx,ky,nz)
-      d3db::c_ptranspose2_jk_start(nffts, nb, tmp1, tmp2, tmp1, request_indx, 41);
+      d3db::c_ptranspose2_jk_start(nffts,nb,tmp1,tmp2,tmp1,request_indx,41);
    }
    /**** hilbert mapping ****/
    else 
    {
       // do fft along ny dimension
       // A(ky,nz,kx) <- fft1d[A(ny,nz,kx)]
-      d3db::mygdevice.batch_cffty_stages_tmpy_zero(2,d3db::fft_tag,true,ny,nq2,n2ft3d,tmp1,d3db::tmpy,zero_row2[nb],da_indx);
+      d3db::mygdevice.batch_cffty_stages_tmpy_zero(2,d3db::fft_tag,true,ny,nffts*nq2,n2ft3d,tmp1,d3db::tmpy,zero_row2[nb],da_indx);
      
       // in=tmp2, out=tmp2
       d3db::c_ptranspose_ijk_start(nffts, nb, 1, tmp1, tmp2, tmp1, request_indx, 42);
@@ -3017,7 +3344,7 @@ void PGrid::pfftfz_start(const int nffts, const int nb, double *tmp1, double *tm
      
       // do fft along nz dimension
       // A(kz,kx,ky) <- fft1d[A(nz,kx,ky)]
-      d3db::mygdevice.batch_cfftz_stages_tmpz_zero(0,d3db::fft_tag,true, nz, nq3, n2ft3d, tmp2, d3db::tmpz, zero_row3[nb],da_indx);
+      d3db::mygdevice.batch_cfftz_stages_tmpz_zero(0,d3db::fft_tag,true, nz, nffts*nq3, n2ft3d, tmp2, d3db::tmpz, zero_row3[nb],da_indx);
    }
 }
 
@@ -3057,7 +3384,7 @@ void PGrid::pfftfz_compute(const int nffts, const int nb, double *tmp1, double *
    {
       // do fft along nz dimension
       // A(kz,kx,ky) <- fft1d[A(nz,kx,ky)]
-      d3db::mygdevice.batch_cfftz_stages_tmpz_zero(1,d3db::fft_tag,true, nz, nq3, n2ft3d, tmp2, d3db::tmpz, zero_row3[nb],da_indx);
+      d3db::mygdevice.batch_cfftz_stages_tmpz_zero(1,d3db::fft_tag,true, nz, nffts*nq3, n2ft3d, tmp2, d3db::tmpz, zero_row3[nb],da_indx);
    }
 }
 
@@ -3123,7 +3450,7 @@ void PGrid::pfftfz_end(const int nffts, const int nb, double *tmp1, double *tmp2
    {
       // do fft along nz dimension
       // A(kz,kx,ky) <- fft1d[A(nz,kx,ky)]
-      d3db::mygdevice.batch_cfftz_stages_tmpz_zero(2,d3db::fft_tag,true, nz, nq3, n2ft3d, tmp2, d3db::tmpz, zero_row3[nb],da_indx);
+      d3db::mygdevice.batch_cfftz_stages_tmpz_zero(2,d3db::fft_tag,true, nz, nffts*nq3, n2ft3d, tmp2, d3db::tmpz, zero_row3[nb],da_indx);
    }
 }
 
@@ -3196,19 +3523,22 @@ void PGrid::pfftf_step10(const int step, const int nffts, const int nb, double *
  *                              *
  ********************************/
 void PGrid::c_pack_start(const int nffts, const int nb, double *a, double *tmp1,
-                         const int request_indx, const int msgtype) {
-  // int one=1;
-
-  // DCOPY_PWDFT(n2ft3d,a,one,tmp,one);
-  std::memcpy(tmp1, a, n2ft3d * sizeof(double));
-  std::memset(a, 0, n2ft3d * sizeof(double));
-
-  c_aindexcopy(nida[nb] + nidb2[nb], packarray[nb], tmp1, a);
-
-  if (balanced)
-    mybalance->c_balance_start(nffts, nb, a, request_indx, msgtype);
-
-  return;
+                         const int request_indx, const int msgtype) 
+{
+   // int one=1;
+ 
+   // DCOPY_PWDFT(n2ft3d,a,one,tmp,one);
+   std::memcpy(tmp1, a, nffts*n2ft3d * sizeof(double));
+   std::memset(a, 0, nffts*n2ft3d * sizeof(double));
+ 
+   for (auto s=0; s<nffts; ++s)
+      c_aindexcopy((nida[nb] + nidb2[nb]), packarray[nb], tmp1+s*n2ft3d, a+s*n2ft3d);
+ 
+   if (balanced)
+      for (auto s=0; s<nffts; ++s)
+         mybalance->c_balance_start(nffts, nb, a+s*n2ft3d, request_indx, msgtype);
+ 
+   return;
 }
 
 /********************************
@@ -3216,12 +3546,13 @@ void PGrid::c_pack_start(const int nffts, const int nb, double *a, double *tmp1,
  *       PGrid:c_pack_end       *
  *                              *
  ********************************/
-void PGrid::c_pack_end(const int nffts, const int nb, double *tmp1, const int request_indx) {
-
-  if (balanced)
-    mybalance->c_balance_end(nffts, nb, tmp1, request_indx);
+void PGrid::c_pack_end(const int nffts, const int nb, double *tmp1, const int request_indx) 
+{
+   if (balanced)
+      for (auto s=0; s<nffts; ++s)
+         mybalance->c_balance_end(nffts, nb, tmp1+s*n2ft3d, request_indx);
 
-  return;
+   return;
 }
 
 /********************************
diff --git a/Nwpw/nwpwlib/utilities/util.hpp b/Nwpw/nwpwlib/utilities/util.hpp
index 1557e082..c09de50c 100644
--- a/Nwpw/nwpwlib/utilities/util.hpp
+++ b/Nwpw/nwpwlib/utilities/util.hpp
@@ -15,8 +15,14 @@
 namespace pwdft {
 
 extern void c_aindexcopy(const int, const int *, double *, double *);
+extern void c_aindexcopy_stride(const int, const int, const int *, double *, double *);
+
 extern void c_bindexcopy(const int, const int *, double *, double *);
+extern void c_bindexcopy_stride(const int, const int, const int *, double *, double *);
+
 extern void c_bindexcopy_conjg(const int, const int *, double *, double *);
+extern void c_bindexcopy_conjg_stride(const int, const int, const int *, double *, double *);
+
 extern void c_bindexzero(const int, const int *, double *);
 
 extern void t_aindexcopy(const int, const int *, double *, double *);