prolong3 改为先算实际 stencil 窗口；只有窗口触及对称边界时才走全域 symmetry_bd，否则只复制必需窗口。restrict3 同样改成窗口判定，无触边时仅填 ii/jj/kk 必需窗口。

- compute fi_min/fi_max from output i-range and derive ii_lo/ii_hi
 - replace full ii sweep (-1:extf(1)) with windowed sweep in Z/Y precompute passes
 - keep stencil math unchanged; add bounds sanity check for ii window

AMSS_NCKU_source/MPatch.C

@@ -7,6 +7,7 @@
 #include <string>
 #include <cmath>
 #include <new>
+#include <vector>
 using namespace std;
 
 #include "misc.h"
@@ -17,6 +18,168 @@ using namespace std;
 #include "interp_lb_profile.h"
 #endif
 
+namespace
+{
+struct InterpBlockView
+{
+  Block *bp;
+  double llb[dim];
+  double uub[dim];
+};
+
+struct BlockBinIndex
+{
+  int bins[dim];
+  double lo[dim];
+  double inv[dim];
+  vector<InterpBlockView> views;
+  vector<vector<int>> bin_to_blocks;
+  bool valid;
+
+  BlockBinIndex() : valid(false)
+  {
+    for (int i = 0; i < dim; i++)
+    {
+      bins[i] = 1;
+      lo[i] = 0.0;
+      inv[i] = 0.0;
+    }
+  }
+};
+
+inline int clamp_int(int v, int lo, int hi)
+{
+  return (v < lo) ? lo : ((v > hi) ? hi : v);
+}
+
+inline int coord_to_bin(double x, double lo, double inv, int nb)
+{
+  if (nb <= 1 || inv <= 0.0)
+    return 0;
+  int b = int(floor((x - lo) * inv));
+  return clamp_int(b, 0, nb - 1);
+}
+
+inline int bin_loc(const BlockBinIndex &index, int b0, int b1, int b2)
+{
+  return b0 + index.bins[0] * (b1 + index.bins[1] * b2);
+}
+
+inline bool point_in_block_view(const InterpBlockView &view, const double *pox, const double *DH)
+{
+  for (int i = 0; i < dim; i++)
+  {
+    if (pox[i] - view.llb[i] < -DH[i] / 2 || pox[i] - view.uub[i] > DH[i] / 2)
+      return false;
+  }
+  return true;
+}
+
+void build_block_bin_index(Patch *patch, const double *DH, BlockBinIndex &index)
+{
+  index = BlockBinIndex();
+
+  MyList<Block> *Bp = patch->blb;
+  while (Bp)
+  {
+    Block *BP = Bp->data;
+    InterpBlockView view;
+    view.bp = BP;
+    for (int i = 0; i < dim; i++)
+    {
+#ifdef Vertex
+#ifdef Cell
+#error Both Cell and Vertex are defined
+#endif
+      view.llb[i] = (feq(BP->bbox[i], patch->bbox[i], DH[i] / 2)) ? BP->bbox[i] + patch->lli[i] * DH[i] : BP->bbox[i] + (ghost_width - 0.5) * DH[i];
+      view.uub[i] = (feq(BP->bbox[dim + i], patch->bbox[dim + i], DH[i] / 2)) ? BP->bbox[dim + i] - patch->uui[i] * DH[i] : BP->bbox[dim + i] - (ghost_width - 0.5) * DH[i];
+#else
+#ifdef Cell
+      view.llb[i] = (feq(BP->bbox[i], patch->bbox[i], DH[i] / 2)) ? BP->bbox[i] + patch->lli[i] * DH[i] : BP->bbox[i] + ghost_width * DH[i];
+      view.uub[i] = (feq(BP->bbox[dim + i], patch->bbox[dim + i], DH[i] / 2)) ? BP->bbox[dim + i] - patch->uui[i] * DH[i] : BP->bbox[dim + i] - ghost_width * DH[i];
+#else
+#error Not define Vertex nor Cell
+#endif
+#endif
+    }
+    index.views.push_back(view);
+    if (Bp == patch->ble)
+      break;
+    Bp = Bp->next;
+  }
+
+  const int nblocks = int(index.views.size());
+  if (nblocks <= 0)
+    return;
+
+  int bins_1d = int(ceil(pow(double(nblocks), 1.0 / 3.0)));
+  bins_1d = clamp_int(bins_1d, 1, 32);
+  for (int i = 0; i < dim; i++)
+  {
+    index.bins[i] = bins_1d;
+    index.lo[i] = patch->bbox[i] + patch->lli[i] * DH[i];
+    const double hi = patch->bbox[dim + i] - patch->uui[i] * DH[i];
+    if (hi > index.lo[i] && bins_1d > 1)
+      index.inv[i] = bins_1d / (hi - index.lo[i]);
+    else
+      index.inv[i] = 0.0;
+  }
+
+  index.bin_to_blocks.resize(index.bins[0] * index.bins[1] * index.bins[2]);
+
+  for (int bi = 0; bi < nblocks; bi++)
+  {
+    const InterpBlockView &view = index.views[bi];
+    int bmin[dim], bmax[dim];
+    for (int d = 0; d < dim; d++)
+    {
+      const double low = view.llb[d] - DH[d] / 2;
+      const double up = view.uub[d] + DH[d] / 2;
+      bmin[d] = coord_to_bin(low, index.lo[d], index.inv[d], index.bins[d]);
+      bmax[d] = coord_to_bin(up, index.lo[d], index.inv[d], index.bins[d]);
+      if (bmax[d] < bmin[d])
+      {
+        int t = bmin[d];
+        bmin[d] = bmax[d];
+        bmax[d] = t;
+      }
+    }
+
+    for (int bz = bmin[2]; bz <= bmax[2]; bz++)
+      for (int by = bmin[1]; by <= bmax[1]; by++)
+        for (int bx = bmin[0]; bx <= bmax[0]; bx++)
+          index.bin_to_blocks[bin_loc(index, bx, by, bz)].push_back(bi);
+  }
+
+  index.valid = true;
+}
+
+int find_block_index_for_point(const BlockBinIndex &index, const double *pox, const double *DH)
+{
+  if (!index.valid)
+    return -1;
+
+  const int bx = coord_to_bin(pox[0], index.lo[0], index.inv[0], index.bins[0]);
+  const int by = coord_to_bin(pox[1], index.lo[1], index.inv[1], index.bins[1]);
+  const int bz = coord_to_bin(pox[2], index.lo[2], index.inv[2], index.bins[2]);
+  const vector<int> &cand = index.bin_to_blocks[bin_loc(index, bx, by, bz)];
+
+  for (size_t ci = 0; ci < cand.size(); ci++)
+  {
+    const int bi = cand[ci];
+    if (point_in_block_view(index.views[bi], pox, DH))
+      return bi;
+  }
+
+  // Fallback to full scan for numerical edge cases around bin boundaries.
+  for (size_t bi = 0; bi < index.views.size(); bi++)
+    if (point_in_block_view(index.views[bi], pox, DH))
+      return int(bi);
+
+  return -1;
+}
+} // namespace
+
 Patch::Patch(int DIM, int *shapei, double *bboxi, int levi, bool buflog, int Symmetry) : lev(levi)
 {
 
@@ -367,9 +530,11 @@ void Patch::Interp_Points(MyList<var> *VarList,
   for (int j = 0; j < NN; j++)
     owner_rank[j] = -1;
 
-  double DH[dim], llb[dim], uub[dim];
+  double DH[dim];
   for (int i = 0; i < dim; i++)
     DH[i] = getdX(i);
+  BlockBinIndex block_index;
+  build_block_bin_index(this, DH, block_index);
 
   for (int j = 0; j < NN; j++) // run along points
   {
@@ -392,57 +557,24 @@ void Patch::Interp_Points(MyList<var> *VarList,
       }
     }
 
-    MyList<Block> *Bp = blb;
-    bool notfind = true;
-    while (notfind && Bp) // run along Blocks
+    const int block_i = find_block_index_for_point(block_index, pox, DH);
+    if (block_i >= 0)
     {
-      Block *BP = Bp->data;
-
-      bool flag = true;
-      for (int i = 0; i < dim; i++)
+      Block *BP = block_index.views[block_i].bp;
+      owner_rank[j] = BP->rank;
+      if (myrank == BP->rank)
       {
-#ifdef Vertex
-#ifdef Cell
-#error Both Cell and Vertex are defined
-#endif
-        llb[i] = (feq(BP->bbox[i], bbox[i], DH[i] / 2)) ? BP->bbox[i] + lli[i] * DH[i] : BP->bbox[i] + (ghost_width - 0.5) * DH[i];
-        uub[i] = (feq(BP->bbox[dim + i], bbox[dim + i], DH[i] / 2)) ? BP->bbox[dim + i] - uui[i] * DH[i] : BP->bbox[dim + i] - (ghost_width - 0.5) * DH[i];
-#else
-#ifdef Cell
-        llb[i] = (feq(BP->bbox[i], bbox[i], DH[i] / 2)) ? BP->bbox[i] + lli[i] * DH[i] : BP->bbox[i] + ghost_width * DH[i];
-        uub[i] = (feq(BP->bbox[dim + i], bbox[dim + i], DH[i] / 2)) ? BP->bbox[dim + i] - uui[i] * DH[i] : BP->bbox[dim + i] - ghost_width * DH[i];
-#else
-#error Not define Vertex nor Cell
-#endif
-#endif
-        if (XX[i][j] - llb[i] < -DH[i] / 2 || XX[i][j] - uub[i] > DH[i] / 2)
+        //---> interpolation
+        varl = VarList;
+        int k = 0;
+        while (varl) // run along variables
         {
-          flag = false;
-          break;
+          f_global_interp(BP->shape, BP->X[0], BP->X[1], BP->X[2], BP->fgfs[varl->data->sgfn], Shellf[j * num_var + k],
+                          pox[0], pox[1], pox[2], ordn, varl->data->SoA, Symmetry);
+          varl = varl->next;
+          k++;
         }
       }
-
-      if (flag)
-      {
-        notfind = false;
-        owner_rank[j] = BP->rank;
-        if (myrank == BP->rank)
-        {
-          //---> interpolation
-          varl = VarList;
-          int k = 0;
-          while (varl) // run along variables
-          {
-            f_global_interp(BP->shape, BP->X[0], BP->X[1], BP->X[2], BP->fgfs[varl->data->sgfn], Shellf[j * num_var + k],
-                            pox[0], pox[1], pox[2], ordn, varl->data->SoA, Symmetry);
-            varl = varl->next;
-            k++;
-          }
-        }
-      }
-      if (Bp == ble)
-        break;
-      Bp = Bp->next;
     }
   }
 
@@ -535,9 +667,11 @@ void Patch::Interp_Points(MyList<var> *VarList,
   for (int j = 0; j < NN; j++)
     owner_rank[j] = -1;
 
-  double DH[dim], llb[dim], uub[dim];
+  double DH[dim];
   for (int i = 0; i < dim; i++)
     DH[i] = getdX(i);
+  BlockBinIndex block_index;
+  build_block_bin_index(this, DH, block_index);
 
   // --- Interpolation phase (identical to original) ---
   for (int j = 0; j < NN; j++)
@@ -561,56 +695,23 @@ void Patch::Interp_Points(MyList<var> *VarList,
       }
     }
 
-    MyList<Block> *Bp = blb;
-    bool notfind = true;
-    while (notfind && Bp)
+    const int block_i = find_block_index_for_point(block_index, pox, DH);
+    if (block_i >= 0)
     {
-      Block *BP = Bp->data;
-
-      bool flag = true;
-      for (int i = 0; i < dim; i++)
+      Block *BP = block_index.views[block_i].bp;
+      owner_rank[j] = BP->rank;
+      if (myrank == BP->rank)
       {
-#ifdef Vertex
-#ifdef Cell
-#error Both Cell and Vertex are defined
-#endif
-        llb[i] = (feq(BP->bbox[i], bbox[i], DH[i] / 2)) ? BP->bbox[i] + lli[i] * DH[i] : BP->bbox[i] + (ghost_width - 0.5) * DH[i];
-        uub[i] = (feq(BP->bbox[dim + i], bbox[dim + i], DH[i] / 2)) ? BP->bbox[dim + i] - uui[i] * DH[i] : BP->bbox[dim + i] - (ghost_width - 0.5) * DH[i];
-#else
-#ifdef Cell
-        llb[i] = (feq(BP->bbox[i], bbox[i], DH[i] / 2)) ? BP->bbox[i] + lli[i] * DH[i] : BP->bbox[i] + ghost_width * DH[i];
-        uub[i] = (feq(BP->bbox[dim + i], bbox[dim + i], DH[i] / 2)) ? BP->bbox[dim + i] - uui[i] * DH[i] : BP->bbox[dim + i] - ghost_width * DH[i];
-#else
-#error Not define Vertex nor Cell
-#endif
-#endif
-        if (XX[i][j] - llb[i] < -DH[i] / 2 || XX[i][j] - uub[i] > DH[i] / 2)
+        varl = VarList;
+        int k = 0;
+        while (varl)
         {
-          flag = false;
-          break;
+          f_global_interp(BP->shape, BP->X[0], BP->X[1], BP->X[2], BP->fgfs[varl->data->sgfn], Shellf[j * num_var + k],
+                          pox[0], pox[1], pox[2], ordn, varl->data->SoA, Symmetry);
+          varl = varl->next;
+          k++;
         }
       }
-
-      if (flag)
-      {
-        notfind = false;
-        owner_rank[j] = BP->rank;
-        if (myrank == BP->rank)
-        {
-          varl = VarList;
-          int k = 0;
-          while (varl)
-          {
-            f_global_interp(BP->shape, BP->X[0], BP->X[1], BP->X[2], BP->fgfs[varl->data->sgfn], Shellf[j * num_var + k],
-                            pox[0], pox[1], pox[2], ordn, varl->data->SoA, Symmetry);
-            varl = varl->next;
-            k++;
-          }
-        }
-      }
-      if (Bp == ble)
-        break;
-      Bp = Bp->next;
     }
   }
 
@@ -833,9 +934,11 @@ void Patch::Interp_Points(MyList<var> *VarList,
   MPI_Comm_group(MPI_COMM_WORLD, &world_group);
   MPI_Comm_group(Comm_here, &local_group);
 
-  double DH[dim], llb[dim], uub[dim];
+  double DH[dim];
   for (int i = 0; i < dim; i++)
     DH[i] = getdX(i);
+  BlockBinIndex block_index;
+  build_block_bin_index(this, DH, block_index);
 
   for (int j = 0; j < NN; j++) // run along points
   {
@@ -858,57 +961,24 @@ void Patch::Interp_Points(MyList<var> *VarList,
       }
     }
 
-    MyList<Block> *Bp = blb;
-    bool notfind = true;
-    while (notfind && Bp) // run along Blocks
+    const int block_i = find_block_index_for_point(block_index, pox, DH);
+    if (block_i >= 0)
     {
-      Block *BP = Bp->data;
-
-      bool flag = true;
-      for (int i = 0; i < dim; i++)
+      Block *BP = block_index.views[block_i].bp;
+      owner_rank[j] = BP->rank;
+      if (myrank == BP->rank)
       {
-#ifdef Vertex
-#ifdef Cell
-#error Both Cell and Vertex are defined
-#endif
-        llb[i] = (feq(BP->bbox[i], bbox[i], DH[i] / 2)) ? BP->bbox[i] + lli[i] * DH[i] : BP->bbox[i] + (ghost_width - 0.5) * DH[i];
-        uub[i] = (feq(BP->bbox[dim + i], bbox[dim + i], DH[i] / 2)) ? BP->bbox[dim + i] - uui[i] * DH[i] : BP->bbox[dim + i] - (ghost_width - 0.5) * DH[i];
-#else
-#ifdef Cell
-        llb[i] = (feq(BP->bbox[i], bbox[i], DH[i] / 2)) ? BP->bbox[i] + lli[i] * DH[i] : BP->bbox[i] + ghost_width * DH[i];
-        uub[i] = (feq(BP->bbox[dim + i], bbox[dim + i], DH[i] / 2)) ? BP->bbox[dim + i] - uui[i] * DH[i] : BP->bbox[dim + i] - ghost_width * DH[i];
-#else
-#error Not define Vertex nor Cell
-#endif
-#endif
-        if (XX[i][j] - llb[i] < -DH[i] / 2 || XX[i][j] - uub[i] > DH[i] / 2)
+        //---> interpolation
+        varl = VarList;
+        int k = 0;
+        while (varl) // run along variables
         {
-          flag = false;
-          break;
+          f_global_interp(BP->shape, BP->X[0], BP->X[1], BP->X[2], BP->fgfs[varl->data->sgfn], Shellf[j * num_var + k],
+                          pox[0], pox[1], pox[2], ordn, varl->data->SoA, Symmetry);
+          varl = varl->next;
+          k++;
         }
       }
-
-      if (flag)
-      {
-        notfind = false;
-        owner_rank[j] = BP->rank;
-        if (myrank == BP->rank)
-        {
-          //---> interpolation
-          varl = VarList;
-          int k = 0;
-          while (varl) // run along variables
-          {
-            f_global_interp(BP->shape, BP->X[0], BP->X[1], BP->X[2], BP->fgfs[varl->data->sgfn], Shellf[j * num_var + k],
-                            pox[0], pox[1], pox[2], ordn, varl->data->SoA, Symmetry);
-            varl = varl->next;
-            k++;
-          }
-        }
-      }
-      if (Bp == ble)
-        break;
-      Bp = Bp->next;
     }
   }
 

AMSS_NCKU_source/Parallel.C

@@ -3893,66 +3893,105 @@ void Parallel::transfer(MyList<Parallel::gridseg> **src, MyList<Parallel::gridse
 
   int node;
 
-  MPI_Request *reqs;
-  MPI_Status *stats;
-  reqs = new MPI_Request[2 * cpusize];
-  stats = new MPI_Status[2 * cpusize];
+  MPI_Request *reqs = new MPI_Request[2 * cpusize];
+  MPI_Status *stats = new MPI_Status[2 * cpusize];
+  int *req_node = new int[2 * cpusize];
+  int *req_is_recv = new int[2 * cpusize];
+  int *completed = new int[2 * cpusize];
   int req_no = 0;
+  int pending_recv = 0;
 
-  double **send_data, **rec_data;
-  send_data = new double *[cpusize];
-  rec_data = new double *[cpusize];
-  int length;
+  double **send_data = new double *[cpusize];
+  double **rec_data = new double *[cpusize];
+  int *send_lengths = new int[cpusize];
+  int *recv_lengths = new int[cpusize];
 
   for (node = 0; node < cpusize; node++)
   {
     send_data[node] = rec_data[node] = 0;
-    if (node == myrank)
+    send_lengths[node] = recv_lengths[node] = 0;
+  }
+
+  // Post receives first so peers can progress rendezvous early.
+  for (node = 0; node < cpusize; node++)
+  {
+    if (node == myrank) continue;
+
+    recv_lengths[node] = data_packer(0, src[node], dst[node], node, UNPACK, VarList1, VarList2, Symmetry);
+    if (recv_lengths[node] > 0)
     {
-      if (length = data_packer(0, src[myrank], dst[myrank], node, PACK, VarList1, VarList2, Symmetry))
+      rec_data[node] = new double[recv_lengths[node]];
+      if (!rec_data[node])
       {
-        rec_data[node] = new double[length];
-        if (!rec_data[node])
-        {
-          cout << "out of memory when new in short transfer, place 1" << endl;
-          MPI_Abort(MPI_COMM_WORLD, 1);
-        }
-        data_packer(rec_data[node], src[myrank], dst[myrank], node, PACK, VarList1, VarList2, Symmetry);
+        cout << "out of memory when new in short transfer, place 1" << endl;
+        MPI_Abort(MPI_COMM_WORLD, 1);
       }
+      MPI_Irecv((void *)rec_data[node], recv_lengths[node], MPI_DOUBLE, node, 1, MPI_COMM_WORLD, reqs + req_no);
+      req_node[req_no] = node;
+      req_is_recv[req_no] = 1;
+      req_no++;
+      pending_recv++;
     }
-    else
+  }
+
+  // Local transfer on this rank.
+  recv_lengths[myrank] = data_packer(0, src[myrank], dst[myrank], myrank, PACK, VarList1, VarList2, Symmetry);
+  if (recv_lengths[myrank] > 0)
+  {
+    rec_data[myrank] = new double[recv_lengths[myrank]];
+    if (!rec_data[myrank])
     {
-      // send from this cpu to cpu#node
-      if (length = data_packer(0, src[myrank], dst[myrank], node, PACK, VarList1, VarList2, Symmetry))
+      cout << "out of memory when new in short transfer, place 2" << endl;
+      MPI_Abort(MPI_COMM_WORLD, 1);
+    }
+    data_packer(rec_data[myrank], src[myrank], dst[myrank], myrank, PACK, VarList1, VarList2, Symmetry);
+  }
+
+  // Pack and post sends.
+  for (node = 0; node < cpusize; node++)
+  {
+    if (node == myrank) continue;
+
+    send_lengths[node] = data_packer(0, src[myrank], dst[myrank], node, PACK, VarList1, VarList2, Symmetry);
+    if (send_lengths[node] > 0)
+    {
+      send_data[node] = new double[send_lengths[node]];
+      if (!send_data[node])
       {
-        send_data[node] = new double[length];
-        if (!send_data[node])
-        {
-          cout << "out of memory when new in short transfer, place 2" << endl;
-          MPI_Abort(MPI_COMM_WORLD, 1);
-        }
-        data_packer(send_data[node], src[myrank], dst[myrank], node, PACK, VarList1, VarList2, Symmetry);
-        MPI_Isend((void *)send_data[node], length, MPI_DOUBLE, node, 1, MPI_COMM_WORLD, reqs + req_no++);
+        cout << "out of memory when new in short transfer, place 3" << endl;
+        MPI_Abort(MPI_COMM_WORLD, 1);
       }
-      // receive from cpu#node to this cpu
-      if (length = data_packer(0, src[node], dst[node], node, UNPACK, VarList1, VarList2, Symmetry))
+      data_packer(send_data[node], src[myrank], dst[myrank], node, PACK, VarList1, VarList2, Symmetry);
+      MPI_Isend((void *)send_data[node], send_lengths[node], MPI_DOUBLE, node, 1, MPI_COMM_WORLD, reqs + req_no);
+      req_node[req_no] = node;
+      req_is_recv[req_no] = 0;
+      req_no++;
+    }
+  }
+
+  // Unpack as soon as receive completes to reduce pure wait time.
+  while (pending_recv > 0)
+  {
+    int outcount = 0;
+    MPI_Waitsome(req_no, reqs, &outcount, completed, stats);
+    if (outcount == MPI_UNDEFINED) break;
+
+    for (int i = 0; i < outcount; i++)
+    {
+      int idx = completed[i];
+      if (idx >= 0 && req_is_recv[idx])
       {
-        rec_data[node] = new double[length];
-        if (!rec_data[node])
-        {
-          cout << "out of memory when new in short transfer, place 3" << endl;
-          MPI_Abort(MPI_COMM_WORLD, 1);
-        }
-        MPI_Irecv((void *)rec_data[node], length, MPI_DOUBLE, node, 1, MPI_COMM_WORLD, reqs + req_no++);
+        int recv_node = req_node[idx];
+        data_packer(rec_data[recv_node], src[recv_node], dst[recv_node], recv_node, UNPACK, VarList1, VarList2, Symmetry);
+        pending_recv--;
       }
     }
   }
-  // wait for all requests to complete
-  MPI_Waitall(req_no, reqs, stats);
 
-  for (node = 0; node < cpusize; node++)
-    if (rec_data[node])
-      data_packer(rec_data[node], src[node], dst[node], node, UNPACK, VarList1, VarList2, Symmetry);
+  if (req_no > 0) MPI_Waitall(req_no, reqs, stats);
+
+  if (rec_data[myrank])
+    data_packer(rec_data[myrank], src[myrank], dst[myrank], myrank, UNPACK, VarList1, VarList2, Symmetry);
 
   for (node = 0; node < cpusize; node++)
   {
@@ -3964,8 +4003,13 @@ void Parallel::transfer(MyList<Parallel::gridseg> **src, MyList<Parallel::gridse
 
   delete[] reqs;
   delete[] stats;
+  delete[] req_node;
+  delete[] req_is_recv;
+  delete[] completed;
   delete[] send_data;
   delete[] rec_data;
+  delete[] send_lengths;
+  delete[] recv_lengths;
 }
 //
 void Parallel::transfermix(MyList<Parallel::gridseg> **src, MyList<Parallel::gridseg> **dst,
@@ -3978,66 +4022,105 @@ void Parallel::transfermix(MyList<Parallel::gridseg> **src, MyList<Parallel::gri
 
   int node;
 
-  MPI_Request *reqs;
-  MPI_Status *stats;
-  reqs = new MPI_Request[2 * cpusize];
-  stats = new MPI_Status[2 * cpusize];
+  MPI_Request *reqs = new MPI_Request[2 * cpusize];
+  MPI_Status *stats = new MPI_Status[2 * cpusize];
+  int *req_node = new int[2 * cpusize];
+  int *req_is_recv = new int[2 * cpusize];
+  int *completed = new int[2 * cpusize];
   int req_no = 0;
+  int pending_recv = 0;
 
-  double **send_data, **rec_data;
-  send_data = new double *[cpusize];
-  rec_data = new double *[cpusize];
-  int length;
+  double **send_data = new double *[cpusize];
+  double **rec_data = new double *[cpusize];
+  int *send_lengths = new int[cpusize];
+  int *recv_lengths = new int[cpusize];
 
   for (node = 0; node < cpusize; node++)
   {
     send_data[node] = rec_data[node] = 0;
-    if (node == myrank)
+    send_lengths[node] = recv_lengths[node] = 0;
+  }
+
+  // Post receives first so peers can progress rendezvous early.
+  for (node = 0; node < cpusize; node++)
+  {
+    if (node == myrank) continue;
+
+    recv_lengths[node] = data_packermix(0, src[node], dst[node], node, UNPACK, VarList1, VarList2, Symmetry);
+    if (recv_lengths[node] > 0)
     {
-      if (length = data_packermix(0, src[myrank], dst[myrank], node, PACK, VarList1, VarList2, Symmetry))
+      rec_data[node] = new double[recv_lengths[node]];
+      if (!rec_data[node])
       {
-        rec_data[node] = new double[length];
-        if (!rec_data[node])
-        {
-          cout << "out of memory when new in short transfer, place 1" << endl;
-          MPI_Abort(MPI_COMM_WORLD, 1);
-        }
-        data_packermix(rec_data[node], src[myrank], dst[myrank], node, PACK, VarList1, VarList2, Symmetry);
+        cout << "out of memory when new in short transfer, place 1" << endl;
+        MPI_Abort(MPI_COMM_WORLD, 1);
       }
+      MPI_Irecv((void *)rec_data[node], recv_lengths[node], MPI_DOUBLE, node, 1, MPI_COMM_WORLD, reqs + req_no);
+      req_node[req_no] = node;
+      req_is_recv[req_no] = 1;
+      req_no++;
+      pending_recv++;
     }
-    else
+  }
+
+  // Local transfer on this rank.
+  recv_lengths[myrank] = data_packermix(0, src[myrank], dst[myrank], myrank, PACK, VarList1, VarList2, Symmetry);
+  if (recv_lengths[myrank] > 0)
+  {
+    rec_data[myrank] = new double[recv_lengths[myrank]];
+    if (!rec_data[myrank])
     {
-      // send from this cpu to cpu#node
-      if (length = data_packermix(0, src[myrank], dst[myrank], node, PACK, VarList1, VarList2, Symmetry))
+      cout << "out of memory when new in short transfer, place 2" << endl;
+      MPI_Abort(MPI_COMM_WORLD, 1);
+    }
+    data_packermix(rec_data[myrank], src[myrank], dst[myrank], myrank, PACK, VarList1, VarList2, Symmetry);
+  }
+
+  // Pack and post sends.
+  for (node = 0; node < cpusize; node++)
+  {
+    if (node == myrank) continue;
+
+    send_lengths[node] = data_packermix(0, src[myrank], dst[myrank], node, PACK, VarList1, VarList2, Symmetry);
+    if (send_lengths[node] > 0)
+    {
+      send_data[node] = new double[send_lengths[node]];
+      if (!send_data[node])
       {
-        send_data[node] = new double[length];
-        if (!send_data[node])
-        {
-          cout << "out of memory when new in short transfer, place 2" << endl;
-          MPI_Abort(MPI_COMM_WORLD, 1);
-        }
-        data_packermix(send_data[node], src[myrank], dst[myrank], node, PACK, VarList1, VarList2, Symmetry);
-        MPI_Isend((void *)send_data[node], length, MPI_DOUBLE, node, 1, MPI_COMM_WORLD, reqs + req_no++);
+        cout << "out of memory when new in short transfer, place 3" << endl;
+        MPI_Abort(MPI_COMM_WORLD, 1);
       }
-      // receive from cpu#node to this cpu
-      if (length = data_packermix(0, src[node], dst[node], node, UNPACK, VarList1, VarList2, Symmetry))
+      data_packermix(send_data[node], src[myrank], dst[myrank], node, PACK, VarList1, VarList2, Symmetry);
+      MPI_Isend((void *)send_data[node], send_lengths[node], MPI_DOUBLE, node, 1, MPI_COMM_WORLD, reqs + req_no);
+      req_node[req_no] = node;
+      req_is_recv[req_no] = 0;
+      req_no++;
+    }
+  }
+
+  // Unpack as soon as receive completes to reduce pure wait time.
+  while (pending_recv > 0)
+  {
+    int outcount = 0;
+    MPI_Waitsome(req_no, reqs, &outcount, completed, stats);
+    if (outcount == MPI_UNDEFINED) break;
+
+    for (int i = 0; i < outcount; i++)
+    {
+      int idx = completed[i];
+      if (idx >= 0 && req_is_recv[idx])
       {
-        rec_data[node] = new double[length];
-        if (!rec_data[node])
-        {
-          cout << "out of memory when new in short transfer, place 3" << endl;
-          MPI_Abort(MPI_COMM_WORLD, 1);
-        }
-        MPI_Irecv((void *)rec_data[node], length, MPI_DOUBLE, node, 1, MPI_COMM_WORLD, reqs + req_no++);
+        int recv_node = req_node[idx];
+        data_packermix(rec_data[recv_node], src[recv_node], dst[recv_node], recv_node, UNPACK, VarList1, VarList2, Symmetry);
+        pending_recv--;
       }
     }
   }
-  // wait for all requests to complete
-  MPI_Waitall(req_no, reqs, stats);
 
-  for (node = 0; node < cpusize; node++)
-    if (rec_data[node])
-      data_packermix(rec_data[node], src[node], dst[node], node, UNPACK, VarList1, VarList2, Symmetry);
+  if (req_no > 0) MPI_Waitall(req_no, reqs, stats);
+
+  if (rec_data[myrank])
+    data_packermix(rec_data[myrank], src[myrank], dst[myrank], myrank, UNPACK, VarList1, VarList2, Symmetry);
 
   for (node = 0; node < cpusize; node++)
   {
@@ -4049,8 +4132,13 @@ void Parallel::transfermix(MyList<Parallel::gridseg> **src, MyList<Parallel::gri
 
   delete[] reqs;
   delete[] stats;
+  delete[] req_node;
+  delete[] req_is_recv;
+  delete[] completed;
   delete[] send_data;
   delete[] rec_data;
+  delete[] send_lengths;
+  delete[] recv_lengths;
 }
 void Parallel::Sync(Patch *Pat, MyList<var> *VarList, int Symmetry)
 {
@@ -4289,62 +4377,99 @@ void Parallel::transfer_cached(MyList<Parallel::gridseg> **src, MyList<Parallel:
   int cpusize = cache.cpusize;
 
   int req_no = 0;
+  int pending_recv = 0;
   int node;
+  int *req_node = new int[cache.max_reqs];
+  int *req_is_recv = new int[cache.max_reqs];
+  int *completed = new int[cache.max_reqs];
 
+  // Post receives first so peers can progress rendezvous early.
   for (node = 0; node < cpusize; node++)
   {
-    if (node == myrank)
+    if (node == myrank) continue;
+
+    int rlength = data_packer(0, src[node], dst[node], node, UNPACK, VarList1, VarList2, Symmetry);
+    cache.recv_lengths[node] = rlength;
+    if (rlength > 0)
     {
-      int length = data_packer(0, src[myrank], dst[myrank], node, PACK, VarList1, VarList2, Symmetry);
-      cache.recv_lengths[node] = length;
-      if (length > 0)
+      if (rlength > cache.recv_buf_caps[node])
       {
-        if (length > cache.recv_buf_caps[node])
-        {
-          if (cache.recv_bufs[node]) delete[] cache.recv_bufs[node];
-          cache.recv_bufs[node] = new double[length];
-          cache.recv_buf_caps[node] = length;
-        }
-        data_packer(cache.recv_bufs[node], src[myrank], dst[myrank], node, PACK, VarList1, VarList2, Symmetry);
+        if (cache.recv_bufs[node]) delete[] cache.recv_bufs[node];
+        cache.recv_bufs[node] = new double[rlength];
+        cache.recv_buf_caps[node] = rlength;
       }
+      MPI_Irecv((void *)cache.recv_bufs[node], rlength, MPI_DOUBLE, node, 1, MPI_COMM_WORLD, cache.reqs + req_no);
+      req_node[req_no] = node;
+      req_is_recv[req_no] = 1;
+      req_no++;
+      pending_recv++;
     }
-    else
+  }
+
+  // Local transfer on this rank.
+  int self_len = data_packer(0, src[myrank], dst[myrank], myrank, PACK, VarList1, VarList2, Symmetry);
+  cache.recv_lengths[myrank] = self_len;
+  if (self_len > 0)
+  {
+    if (self_len > cache.recv_buf_caps[myrank])
     {
-      // send
-      int slength = data_packer(0, src[myrank], dst[myrank], node, PACK, VarList1, VarList2, Symmetry);
-      cache.send_lengths[node] = slength;
-      if (slength > 0)
+      if (cache.recv_bufs[myrank]) delete[] cache.recv_bufs[myrank];
+      cache.recv_bufs[myrank] = new double[self_len];
+      cache.recv_buf_caps[myrank] = self_len;
+    }
+    data_packer(cache.recv_bufs[myrank], src[myrank], dst[myrank], myrank, PACK, VarList1, VarList2, Symmetry);
+  }
+
+  // Pack and post sends.
+  for (node = 0; node < cpusize; node++)
+  {
+    if (node == myrank) continue;
+
+    int slength = data_packer(0, src[myrank], dst[myrank], node, PACK, VarList1, VarList2, Symmetry);
+    cache.send_lengths[node] = slength;
+    if (slength > 0)
+    {
+      if (slength > cache.send_buf_caps[node])
       {
-        if (slength > cache.send_buf_caps[node])
-        {
-          if (cache.send_bufs[node]) delete[] cache.send_bufs[node];
-          cache.send_bufs[node] = new double[slength];
-          cache.send_buf_caps[node] = slength;
-        }
-        data_packer(cache.send_bufs[node], src[myrank], dst[myrank], node, PACK, VarList1, VarList2, Symmetry);
-        MPI_Isend((void *)cache.send_bufs[node], slength, MPI_DOUBLE, node, 1, MPI_COMM_WORLD, cache.reqs + req_no++);
+        if (cache.send_bufs[node]) delete[] cache.send_bufs[node];
+        cache.send_bufs[node] = new double[slength];
+        cache.send_buf_caps[node] = slength;
       }
-      // recv
-      int rlength = data_packer(0, src[node], dst[node], node, UNPACK, VarList1, VarList2, Symmetry);
-      cache.recv_lengths[node] = rlength;
-      if (rlength > 0)
+      data_packer(cache.send_bufs[node], src[myrank], dst[myrank], node, PACK, VarList1, VarList2, Symmetry);
+      MPI_Isend((void *)cache.send_bufs[node], slength, MPI_DOUBLE, node, 1, MPI_COMM_WORLD, cache.reqs + req_no);
+      req_node[req_no] = node;
+      req_is_recv[req_no] = 0;
+      req_no++;
+    }
+  }
+
+  // Unpack as soon as receive completes to reduce pure wait time.
+  while (pending_recv > 0)
+  {
+    int outcount = 0;
+    MPI_Waitsome(req_no, cache.reqs, &outcount, completed, cache.stats);
+    if (outcount == MPI_UNDEFINED) break;
+
+    for (int i = 0; i < outcount; i++)
+    {
+      int idx = completed[i];
+      if (idx >= 0 && req_is_recv[idx])
       {
-        if (rlength > cache.recv_buf_caps[node])
-        {
-          if (cache.recv_bufs[node]) delete[] cache.recv_bufs[node];
-          cache.recv_bufs[node] = new double[rlength];
-          cache.recv_buf_caps[node] = rlength;
-        }
-        MPI_Irecv((void *)cache.recv_bufs[node], rlength, MPI_DOUBLE, node, 1, MPI_COMM_WORLD, cache.reqs + req_no++);
+        int recv_node_i = req_node[idx];
+        data_packer(cache.recv_bufs[recv_node_i], src[recv_node_i], dst[recv_node_i], recv_node_i, UNPACK, VarList1, VarList2, Symmetry);
+        pending_recv--;
       }
     }
   }
 
-  MPI_Waitall(req_no, cache.reqs, cache.stats);
+  if (req_no > 0) MPI_Waitall(req_no, cache.reqs, cache.stats);
 
-  for (node = 0; node < cpusize; node++)
-    if (cache.recv_bufs[node] && cache.recv_lengths[node] > 0)
-      data_packer(cache.recv_bufs[node], src[node], dst[node], node, UNPACK, VarList1, VarList2, Symmetry);
+  if (self_len > 0)
+    data_packer(cache.recv_bufs[myrank], src[myrank], dst[myrank], myrank, UNPACK, VarList1, VarList2, Symmetry);
+
+  delete[] req_node;
+  delete[] req_is_recv;
+  delete[] completed;
 }
 // Sync_cached: build grid segment lists on first call, reuse on subsequent calls
 void Parallel::Sync_cached(MyList<Patch> *PatL, MyList<var> *VarList, int Symmetry, SyncCache &cache)
@@ -5807,58 +5932,98 @@ void Parallel::OutBdLow2Himix_cached(MyList<Patch> *PatcL, MyList<Patch> *PatfL,
   int cpusize = cache.cpusize;
 
   int req_no = 0;
+  int pending_recv = 0;
+  int *req_node = new int[cache.max_reqs];
+  int *req_is_recv = new int[cache.max_reqs];
+  int *completed = new int[cache.max_reqs];
+
+  // Post receives first so peers can progress rendezvous early.
   for (int node = 0; node < cpusize; node++)
   {
-    if (node == myrank)
+    if (node == myrank) continue;
+
+    int rlength = data_packermix(0, cache.combined_src[node], cache.combined_dst[node], node, UNPACK, VarList1, VarList2, Symmetry);
+    cache.recv_lengths[node] = rlength;
+    if (rlength > 0)
     {
-      int length = data_packermix(0, cache.combined_src[myrank], cache.combined_dst[myrank], node, PACK, VarList1, VarList2, Symmetry);
-      cache.recv_lengths[node] = length;
-      if (length > 0)
+      if (rlength > cache.recv_buf_caps[node])
       {
-        if (length > cache.recv_buf_caps[node])
-        {
-          if (cache.recv_bufs[node]) delete[] cache.recv_bufs[node];
-          cache.recv_bufs[node] = new double[length];
-          cache.recv_buf_caps[node] = length;
-        }
-        data_packermix(cache.recv_bufs[node], cache.combined_src[myrank], cache.combined_dst[myrank], node, PACK, VarList1, VarList2, Symmetry);
+        if (cache.recv_bufs[node]) delete[] cache.recv_bufs[node];
+        cache.recv_bufs[node] = new double[rlength];
+        cache.recv_buf_caps[node] = rlength;
       }
+      MPI_Irecv((void *)cache.recv_bufs[node], rlength, MPI_DOUBLE, node, 1, MPI_COMM_WORLD, cache.reqs + req_no);
+      req_node[req_no] = node;
+      req_is_recv[req_no] = 1;
+      req_no++;
+      pending_recv++;
     }
-    else
+  }
+
+  // Local transfer on this rank.
+  int self_len = data_packermix(0, cache.combined_src[myrank], cache.combined_dst[myrank], myrank, PACK, VarList1, VarList2, Symmetry);
+  cache.recv_lengths[myrank] = self_len;
+  if (self_len > 0)
+  {
+    if (self_len > cache.recv_buf_caps[myrank])
     {
-      int slength = data_packermix(0, cache.combined_src[myrank], cache.combined_dst[myrank], node, PACK, VarList1, VarList2, Symmetry);
-      cache.send_lengths[node] = slength;
-      if (slength > 0)
+      if (cache.recv_bufs[myrank]) delete[] cache.recv_bufs[myrank];
+      cache.recv_bufs[myrank] = new double[self_len];
+      cache.recv_buf_caps[myrank] = self_len;
+    }
+    data_packermix(cache.recv_bufs[myrank], cache.combined_src[myrank], cache.combined_dst[myrank], myrank, PACK, VarList1, VarList2, Symmetry);
+  }
+
+  // Pack and post sends.
+  for (int node = 0; node < cpusize; node++)
+  {
+    if (node == myrank) continue;
+
+    int slength = data_packermix(0, cache.combined_src[myrank], cache.combined_dst[myrank], node, PACK, VarList1, VarList2, Symmetry);
+    cache.send_lengths[node] = slength;
+    if (slength > 0)
+    {
+      if (slength > cache.send_buf_caps[node])
       {
-        if (slength > cache.send_buf_caps[node])
-        {
-          if (cache.send_bufs[node]) delete[] cache.send_bufs[node];
-          cache.send_bufs[node] = new double[slength];
-          cache.send_buf_caps[node] = slength;
-        }
-        data_packermix(cache.send_bufs[node], cache.combined_src[myrank], cache.combined_dst[myrank], node, PACK, VarList1, VarList2, Symmetry);
-        MPI_Isend((void *)cache.send_bufs[node], slength, MPI_DOUBLE, node, 1, MPI_COMM_WORLD, cache.reqs + req_no++);
+        if (cache.send_bufs[node]) delete[] cache.send_bufs[node];
+        cache.send_bufs[node] = new double[slength];
+        cache.send_buf_caps[node] = slength;
       }
-      int rlength = data_packermix(0, cache.combined_src[node], cache.combined_dst[node], node, UNPACK, VarList1, VarList2, Symmetry);
-      cache.recv_lengths[node] = rlength;
-      if (rlength > 0)
+      data_packermix(cache.send_bufs[node], cache.combined_src[myrank], cache.combined_dst[myrank], node, PACK, VarList1, VarList2, Symmetry);
+      MPI_Isend((void *)cache.send_bufs[node], slength, MPI_DOUBLE, node, 1, MPI_COMM_WORLD, cache.reqs + req_no);
+      req_node[req_no] = node;
+      req_is_recv[req_no] = 0;
+      req_no++;
+    }
+  }
+
+  // Unpack as soon as receive completes to reduce pure wait time.
+  while (pending_recv > 0)
+  {
+    int outcount = 0;
+    MPI_Waitsome(req_no, cache.reqs, &outcount, completed, cache.stats);
+    if (outcount == MPI_UNDEFINED) break;
+
+    for (int i = 0; i < outcount; i++)
+    {
+      int idx = completed[i];
+      if (idx >= 0 && req_is_recv[idx])
       {
-        if (rlength > cache.recv_buf_caps[node])
-        {
-          if (cache.recv_bufs[node]) delete[] cache.recv_bufs[node];
-          cache.recv_bufs[node] = new double[rlength];
-          cache.recv_buf_caps[node] = rlength;
-        }
-        MPI_Irecv((void *)cache.recv_bufs[node], rlength, MPI_DOUBLE, node, 1, MPI_COMM_WORLD, cache.reqs + req_no++);
+        int recv_node_i = req_node[idx];
+        data_packermix(cache.recv_bufs[recv_node_i], cache.combined_src[recv_node_i], cache.combined_dst[recv_node_i], recv_node_i, UNPACK, VarList1, VarList2, Symmetry);
+        pending_recv--;
       }
     }
   }
 
-  MPI_Waitall(req_no, cache.reqs, cache.stats);
+  if (req_no > 0) MPI_Waitall(req_no, cache.reqs, cache.stats);
 
-  for (int node = 0; node < cpusize; node++)
-    if (cache.recv_bufs[node] && cache.recv_lengths[node] > 0)
-      data_packermix(cache.recv_bufs[node], cache.combined_src[node], cache.combined_dst[node], node, UNPACK, VarList1, VarList2, Symmetry);
+  if (self_len > 0)
+    data_packermix(cache.recv_bufs[myrank], cache.combined_src[myrank], cache.combined_dst[myrank], myrank, UNPACK, VarList1, VarList2, Symmetry);
+
+  delete[] req_node;
+  delete[] req_is_recv;
+  delete[] completed;
 }
 
 // collect all buffer grid segments or blocks for given patch

AMSS_NCKU_source/bssn_rhs_cuda.h

@@ -0,0 +1,36 @@
+#ifndef BSSN_RHS_CUDA_H
+#define BSSN_RHS_CUDA_H
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+int f_compute_rhs_bssn(int *ex, double &T,
+                       double *X, double *Y, double *Z,
+                       double *chi, double *trK,
+                       double *dxx, double *gxy, double *gxz, double *dyy, double *gyz, double *dzz,
+                       double *Axx, double *Axy, double *Axz, double *Ayy, double *Ayz, double *Azz,
+                       double *Gamx, double *Gamy, double *Gamz,
+                       double *Lap, double *betax, double *betay, double *betaz,
+                       double *dtSfx, double *dtSfy, double *dtSfz,
+                       double *chi_rhs, double *trK_rhs,
+                       double *gxx_rhs, double *gxy_rhs, double *gxz_rhs, double *gyy_rhs, double *gyz_rhs, double *gzz_rhs,
+                       double *Axx_rhs, double *Axy_rhs, double *Axz_rhs, double *Ayy_rhs, double *Ayz_rhs, double *Azz_rhs,
+                       double *Gamx_rhs, double *Gamy_rhs, double *Gamz_rhs,
+                       double *Lap_rhs, double *betax_rhs, double *betay_rhs, double *betaz_rhs,
+                       double *dtSfx_rhs, double *dtSfy_rhs, double *dtSfz_rhs,
+                       double *rho, double *Sx, double *Sy, double *Sz,
+                       double *Sxx, double *Sxy, double *Sxz, double *Syy, double *Syz, double *Szz,
+                       double *Gamxxx, double *Gamxxy, double *Gamxxz, double *Gamxyy, double *Gamxyz, double *Gamxzz,
+                       double *Gamyxx, double *Gamyxy, double *Gamyxz, double *Gamyyy, double *Gamyyz, double *Gamyzz,
+                       double *Gamzxx, double *Gamzxy, double *Gamzxz, double *Gamzyy, double *Gamzyz, double *Gamzzz,
+                       double *Rxx, double *Rxy, double *Rxz, double *Ryy, double *Ryz, double *Rzz,
+                       double *ham_Res, double *movx_Res, double *movy_Res, double *movz_Res,
+                       double *Gmx_Res, double *Gmy_Res, double *Gmz_Res,
+                       int &Symmetry, int &Lev, double &eps, int &co);
+
+#ifdef __cplusplus
+}
+#endif
+
+#endif

AMSS_NCKU_source/fdderivs_c.C

@@ -71,149 +71,99 @@ void fdderivs(const int ex[3],
     const double Fdxdz = F1o144 / (dX * dZ);
     const double Fdydz = F1o144 / (dY * dZ);
 
-    /* 输出清零：fxx,fyy,fzz,fxy,fxz,fyz = 0 */
-    const size_t all = (size_t)ex1 * (size_t)ex2 * (size_t)ex3;
-    for (size_t p = 0; p < all; ++p) {
-        fxx[p] = ZEO; fyy[p] = ZEO; fzz[p] = ZEO;
-        fxy[p] = ZEO; fxz[p] = ZEO; fyz[p] = ZEO;
+    /* 只清零不被主循环覆盖的边界面 */
+    {
+        /* 高边界：k0=ex3-1 */
+        for (int j0 = 0; j0 < ex2; ++j0)
+            for (int i0 = 0; i0 < ex1; ++i0) {
+                const size_t p = idx_ex(i0, j0, ex3 - 1, ex);
+                fxx[p]=ZEO; fyy[p]=ZEO; fzz[p]=ZEO;
+                fxy[p]=ZEO; fxz[p]=ZEO; fyz[p]=ZEO;
+            }
+        /* 高边界：j0=ex2-1 */
+        for (int k0 = 0; k0 < ex3 - 1; ++k0)
+            for (int i0 = 0; i0 < ex1; ++i0) {
+                const size_t p = idx_ex(i0, ex2 - 1, k0, ex);
+                fxx[p]=ZEO; fyy[p]=ZEO; fzz[p]=ZEO;
+                fxy[p]=ZEO; fxz[p]=ZEO; fyz[p]=ZEO;
+            }
+        /* 高边界：i0=ex1-1 */
+        for (int k0 = 0; k0 < ex3 - 1; ++k0)
+            for (int j0 = 0; j0 < ex2 - 1; ++j0) {
+                const size_t p = idx_ex(ex1 - 1, j0, k0, ex);
+                fxx[p]=ZEO; fyy[p]=ZEO; fzz[p]=ZEO;
+                fxy[p]=ZEO; fxz[p]=ZEO; fyz[p]=ZEO;
+            }
+
+        /* 低边界：当二阶模板也不可用时，对应 i0/j0/k0=0 面 */
+        if (kminF == 1) {
+            for (int j0 = 0; j0 < ex2; ++j0)
+                for (int i0 = 0; i0 < ex1; ++i0) {
+                    const size_t p = idx_ex(i0, j0, 0, ex);
+                    fxx[p]=ZEO; fyy[p]=ZEO; fzz[p]=ZEO;
+                    fxy[p]=ZEO; fxz[p]=ZEO; fyz[p]=ZEO;
+                }
+        }
+        if (jminF == 1) {
+            for (int k0 = 0; k0 < ex3; ++k0)
+                for (int i0 = 0; i0 < ex1; ++i0) {
+                    const size_t p = idx_ex(i0, 0, k0, ex);
+                    fxx[p]=ZEO; fyy[p]=ZEO; fzz[p]=ZEO;
+                    fxy[p]=ZEO; fxz[p]=ZEO; fyz[p]=ZEO;
+                }
+        }
+        if (iminF == 1) {
+            for (int k0 = 0; k0 < ex3; ++k0)
+                for (int j0 = 0; j0 < ex2; ++j0) {
+                    const size_t p = idx_ex(0, j0, k0, ex);
+                    fxx[p]=ZEO; fyy[p]=ZEO; fzz[p]=ZEO;
+                    fxy[p]=ZEO; fxz[p]=ZEO; fyz[p]=ZEO;
+                }
+        }
     }
 
     /*
-     * Fortran:
-     * do k=1,ex3-1
-     * do j=1,ex2-1
-     * do i=1,ex1-1
+     * 两段式：
+     * 1) 二阶可用区域先计算二阶模板
+     * 2) 高阶可用区域再覆盖四阶模板
      */
+    const int i2_lo = (iminF > 0) ? iminF : 0;
+    const int j2_lo = (jminF > 0) ? jminF : 0;
+    const int k2_lo = (kminF > 0) ? kminF : 0;
+    const int i2_hi = ex1 - 2;
+    const int j2_hi = ex2 - 2;
+    const int k2_hi = ex3 - 2;
 
-    for (int k0 = 0; k0 <= ex3 - 2; ++k0) {
-        const int kF = k0 + 1;
-        for (int j0 = 0; j0 <= ex2 - 2; ++j0) {
-            const int jF = j0 + 1;
-            for (int i0 = 0; i0 <= ex1 - 2; ++i0) {
-                const int iF = i0 + 1;
-                const size_t p = idx_ex(i0, j0, k0, ex);
+    const int i4_lo = (iminF + 1 > 0) ? (iminF + 1) : 0;
+    const int j4_lo = (jminF + 1 > 0) ? (jminF + 1) : 0;
+    const int k4_lo = (kminF + 1 > 0) ? (kminF + 1) : 0;
+    const int i4_hi = ex1 - 3;
+    const int j4_hi = ex2 - 3;
+    const int k4_hi = ex3 - 3;
 
-                /* 高阶分支：i±2,j±2,k±2 都在范围内 */
-                if ((iF + 2) <= imaxF && (iF - 2) >= iminF &&
-                    (jF + 2) <= jmaxF && (jF - 2) >= jminF &&
-                    (kF + 2) <= kmaxF && (kF - 2) >= kminF)
-                {
-                    fxx[p] = Fdxdx * (
-                        -fh[idx_fh_F_ord2(iF - 2, jF,     kF,     ex)] +
-                        F16 * fh[idx_fh_F_ord2(iF - 1, jF,     kF,     ex)] -
-                        F30 * fh[idx_fh_F_ord2(iF,     jF,     kF,     ex)] -
-                        fh[idx_fh_F_ord2(iF + 2, jF,     kF,     ex)] +
-                        F16 * fh[idx_fh_F_ord2(iF + 1, jF,     kF,     ex)]
-                    );
+    /*
+     * Strategy A:
+     * Avoid redundant work in overlap of 2nd/4th-order regions.
+     * Only compute 2nd-order on shell points that are NOT overwritten by
+     * the 4th-order pass.
+     */
+    const int has4 = (i4_lo <= i4_hi && j4_lo <= j4_hi && k4_lo <= k4_hi);
 
-                    fyy[p] = Fdydy * (
-                        -fh[idx_fh_F_ord2(iF,     jF - 2, kF,     ex)] +
-                        F16 * fh[idx_fh_F_ord2(iF,     jF - 1, kF,     ex)] -
-                        F30 * fh[idx_fh_F_ord2(iF,     jF,     kF,     ex)] -
-                        fh[idx_fh_F_ord2(iF,     jF + 2, kF,     ex)] +
-                        F16 * fh[idx_fh_F_ord2(iF,     jF + 1, kF,     ex)]
-                    );
-
-                    fzz[p] = Fdzdz * (
-                        -fh[idx_fh_F_ord2(iF,     jF,     kF - 2, ex)] +
-                        F16 * fh[idx_fh_F_ord2(iF,     jF,     kF - 1, ex)] -
-                        F30 * fh[idx_fh_F_ord2(iF,     jF,     kF,     ex)] -
-                        fh[idx_fh_F_ord2(iF,     jF,     kF + 2, ex)] +
-                        F16 * fh[idx_fh_F_ord2(iF,     jF,     kF + 1, ex)]
-                    );
-
-                    /* fxy 高阶：完全照搬 Fortran 的括号结构 */
-                    {
-                        const double t_jm2 =
-                            ( fh[idx_fh_F_ord2(iF - 2, jF - 2, kF, ex)]
-                             -F8*fh[idx_fh_F_ord2(iF - 1, jF - 2, kF, ex)]
-                             +F8*fh[idx_fh_F_ord2(iF + 1, jF - 2, kF, ex)]
-                             -    fh[idx_fh_F_ord2(iF + 2, jF - 2, kF, ex)] );
-
-                        const double t_jm1 =
-                            ( fh[idx_fh_F_ord2(iF - 2, jF - 1, kF, ex)]
-                             -F8*fh[idx_fh_F_ord2(iF - 1, jF - 1, kF, ex)]
-                             +F8*fh[idx_fh_F_ord2(iF + 1, jF - 1, kF, ex)]
-                             -    fh[idx_fh_F_ord2(iF + 2, jF - 1, kF, ex)] );
-
-                        const double t_jp1 =
-                            ( fh[idx_fh_F_ord2(iF - 2, jF + 1, kF, ex)]
-                             -F8*fh[idx_fh_F_ord2(iF - 1, jF + 1, kF, ex)]
-                             +F8*fh[idx_fh_F_ord2(iF + 1, jF + 1, kF, ex)]
-                             -    fh[idx_fh_F_ord2(iF + 2, jF + 1, kF, ex)] );
-
-                        const double t_jp2 =
-                            ( fh[idx_fh_F_ord2(iF - 2, jF + 2, kF, ex)]
-                             -F8*fh[idx_fh_F_ord2(iF - 1, jF + 2, kF, ex)]
-                             +F8*fh[idx_fh_F_ord2(iF + 1, jF + 2, kF, ex)]
-                             -    fh[idx_fh_F_ord2(iF + 2, jF + 2, kF, ex)] );
-
-                        fxy[p] = Fdxdy * ( t_jm2 - F8 * t_jm1 + F8 * t_jp1 - t_jp2 );
+    if (i2_lo <= i2_hi && j2_lo <= j2_hi && k2_lo <= k2_hi) {
+        for (int k0 = k2_lo; k0 <= k2_hi; ++k0) {
+            const int kF = k0 + 1;
+            for (int j0 = j2_lo; j0 <= j2_hi; ++j0) {
+                const int jF = j0 + 1;
+                for (int i0 = i2_lo; i0 <= i2_hi; ++i0) {
+                    if (has4 &&
+                        i0 >= i4_lo && i0 <= i4_hi &&
+                        j0 >= j4_lo && j0 <= j4_hi &&
+                        k0 >= k4_lo && k0 <= k4_hi) {
+                        continue;
                     }
+                    const int iF = i0 + 1;
+                    const size_t p = idx_ex(i0, j0, k0, ex);
 
-                    /* fxz 高阶 */
-                    {
-                        const double t_km2 =
-                            ( fh[idx_fh_F_ord2(iF - 2, jF, kF - 2, ex)]
-                             -F8*fh[idx_fh_F_ord2(iF - 1, jF, kF - 2, ex)]
-                             +F8*fh[idx_fh_F_ord2(iF + 1, jF, kF - 2, ex)]
-                             -    fh[idx_fh_F_ord2(iF + 2, jF, kF - 2, ex)] );
-
-                        const double t_km1 =
-                            ( fh[idx_fh_F_ord2(iF - 2, jF, kF - 1, ex)]
-                             -F8*fh[idx_fh_F_ord2(iF - 1, jF, kF - 1, ex)]
-                             +F8*fh[idx_fh_F_ord2(iF + 1, jF, kF - 1, ex)]
-                             -    fh[idx_fh_F_ord2(iF + 2, jF, kF - 1, ex)] );
-
-                        const double t_kp1 =
-                            ( fh[idx_fh_F_ord2(iF - 2, jF, kF + 1, ex)]
-                             -F8*fh[idx_fh_F_ord2(iF - 1, jF, kF + 1, ex)]
-                             +F8*fh[idx_fh_F_ord2(iF + 1, jF, kF + 1, ex)]
-                             -    fh[idx_fh_F_ord2(iF + 2, jF, kF + 1, ex)] );
-
-                        const double t_kp2 =
-                            ( fh[idx_fh_F_ord2(iF - 2, jF, kF + 2, ex)]
-                             -F8*fh[idx_fh_F_ord2(iF - 1, jF, kF + 2, ex)]
-                             +F8*fh[idx_fh_F_ord2(iF + 1, jF, kF + 2, ex)]
-                             -    fh[idx_fh_F_ord2(iF + 2, jF, kF + 2, ex)] );
-
-                        fxz[p] = Fdxdz * ( t_km2 - F8 * t_km1 + F8 * t_kp1 - t_kp2 );
-                    }
-
-                    /* fyz 高阶 */
-                    {
-                        const double t_km2 =
-                            ( fh[idx_fh_F_ord2(iF, jF - 2, kF - 2, ex)]
-                             -F8*fh[idx_fh_F_ord2(iF, jF - 1, kF - 2, ex)]
-                             +F8*fh[idx_fh_F_ord2(iF, jF + 1, kF - 2, ex)]
-                             -    fh[idx_fh_F_ord2(iF, jF + 2, kF - 2, ex)] );
-
-                        const double t_km1 =
-                            ( fh[idx_fh_F_ord2(iF, jF - 2, kF - 1, ex)]
-                             -F8*fh[idx_fh_F_ord2(iF, jF - 1, kF - 1, ex)]
-                             +F8*fh[idx_fh_F_ord2(iF, jF + 1, kF - 1, ex)]
-                             -    fh[idx_fh_F_ord2(iF, jF + 2, kF - 1, ex)] );
-
-                        const double t_kp1 =
-                            ( fh[idx_fh_F_ord2(iF, jF - 2, kF + 1, ex)]
-                             -F8*fh[idx_fh_F_ord2(iF, jF - 1, kF + 1, ex)]
-                             +F8*fh[idx_fh_F_ord2(iF, jF + 1, kF + 1, ex)]
-                             -    fh[idx_fh_F_ord2(iF, jF + 2, kF + 1, ex)] );
-
-                        const double t_kp2 =
-                            ( fh[idx_fh_F_ord2(iF, jF - 2, kF + 2, ex)]
-                             -F8*fh[idx_fh_F_ord2(iF, jF - 1, kF + 2, ex)]
-                             +F8*fh[idx_fh_F_ord2(iF, jF + 1, kF + 2, ex)]
-                             -    fh[idx_fh_F_ord2(iF, jF + 2, kF + 2, ex)] );
-
-                        fyz[p] = Fdydz * ( t_km2 - F8 * t_km1 + F8 * t_kp1 - t_kp2 );
-                    }
-                }
-                /* 二阶分支：i±1,j±1,k±1 在范围内 */
-                else if ((iF + 1) <= imaxF && (iF - 1) >= iminF &&
-                         (jF + 1) <= jmaxF && (jF - 1) >= jminF &&
-                         (kF + 1) <= kmaxF && (kF - 1) >= kminF)
-                {
                     fxx[p] = Sdxdx * (
                         fh[idx_fh_F_ord2(iF - 1, jF,     kF,     ex)] -
                         TWO * fh[idx_fh_F_ord2(iF,     jF,     kF,     ex)] +
@@ -252,13 +202,127 @@ void fdderivs(const int ex[3],
                         fh[idx_fh_F_ord2(iF, jF - 1, kF + 1, ex)] +
                         fh[idx_fh_F_ord2(iF, jF + 1, kF + 1, ex)]
                     );
-                }else{
-                    fxx[p] = 0.0;
-                    fyy[p] = 0.0;
-                    fzz[p] = 0.0;
-                    fxy[p] = 0.0;
-                    fxz[p] = 0.0;
-                    fyz[p] = 0.0;
+                }
+            }
+        }
+    }
+
+    if (has4) {
+        for (int k0 = k4_lo; k0 <= k4_hi; ++k0) {
+            const int kF = k0 + 1;
+            for (int j0 = j4_lo; j0 <= j4_hi; ++j0) {
+                const int jF = j0 + 1;
+                for (int i0 = i4_lo; i0 <= i4_hi; ++i0) {
+                    const int iF = i0 + 1;
+                    const size_t p = idx_ex(i0, j0, k0, ex);
+
+                    fxx[p] = Fdxdx * (
+                        -fh[idx_fh_F_ord2(iF - 2, jF,     kF,     ex)] +
+                        F16 * fh[idx_fh_F_ord2(iF - 1, jF,     kF,     ex)] -
+                        F30 * fh[idx_fh_F_ord2(iF,     jF,     kF,     ex)] -
+                        fh[idx_fh_F_ord2(iF + 2, jF,     kF,     ex)] +
+                        F16 * fh[idx_fh_F_ord2(iF + 1, jF,     kF,     ex)]
+                    );
+
+                    fyy[p] = Fdydy * (
+                        -fh[idx_fh_F_ord2(iF,     jF - 2, kF,     ex)] +
+                        F16 * fh[idx_fh_F_ord2(iF,     jF - 1, kF,     ex)] -
+                        F30 * fh[idx_fh_F_ord2(iF,     jF,     kF,     ex)] -
+                        fh[idx_fh_F_ord2(iF,     jF + 2, kF,     ex)] +
+                        F16 * fh[idx_fh_F_ord2(iF,     jF + 1, kF,     ex)]
+                    );
+
+                    fzz[p] = Fdzdz * (
+                        -fh[idx_fh_F_ord2(iF,     jF,     kF - 2, ex)] +
+                        F16 * fh[idx_fh_F_ord2(iF,     jF,     kF - 1, ex)] -
+                        F30 * fh[idx_fh_F_ord2(iF,     jF,     kF,     ex)] -
+                        fh[idx_fh_F_ord2(iF,     jF,     kF + 2, ex)] +
+                        F16 * fh[idx_fh_F_ord2(iF,     jF,     kF + 1, ex)]
+                    );
+
+                    {
+                        const double t_jm2 =
+                            ( fh[idx_fh_F_ord2(iF - 2, jF - 2, kF, ex)]
+                             -F8*fh[idx_fh_F_ord2(iF - 1, jF - 2, kF, ex)]
+                             +F8*fh[idx_fh_F_ord2(iF + 1, jF - 2, kF, ex)]
+                             -    fh[idx_fh_F_ord2(iF + 2, jF - 2, kF, ex)] );
+
+                        const double t_jm1 =
+                            ( fh[idx_fh_F_ord2(iF - 2, jF - 1, kF, ex)]
+                             -F8*fh[idx_fh_F_ord2(iF - 1, jF - 1, kF, ex)]
+                             +F8*fh[idx_fh_F_ord2(iF + 1, jF - 1, kF, ex)]
+                             -    fh[idx_fh_F_ord2(iF + 2, jF - 1, kF, ex)] );
+
+                        const double t_jp1 =
+                            ( fh[idx_fh_F_ord2(iF - 2, jF + 1, kF, ex)]
+                             -F8*fh[idx_fh_F_ord2(iF - 1, jF + 1, kF, ex)]
+                             +F8*fh[idx_fh_F_ord2(iF + 1, jF + 1, kF, ex)]
+                             -    fh[idx_fh_F_ord2(iF + 2, jF + 1, kF, ex)] );
+
+                        const double t_jp2 =
+                            ( fh[idx_fh_F_ord2(iF - 2, jF + 2, kF, ex)]
+                             -F8*fh[idx_fh_F_ord2(iF - 1, jF + 2, kF, ex)]
+                             +F8*fh[idx_fh_F_ord2(iF + 1, jF + 2, kF, ex)]
+                             -    fh[idx_fh_F_ord2(iF + 2, jF + 2, kF, ex)] );
+
+                        fxy[p] = Fdxdy * ( t_jm2 - F8 * t_jm1 + F8 * t_jp1 - t_jp2 );
+                    }
+
+                    {
+                        const double t_km2 =
+                            ( fh[idx_fh_F_ord2(iF - 2, jF, kF - 2, ex)]
+                             -F8*fh[idx_fh_F_ord2(iF - 1, jF, kF - 2, ex)]
+                             +F8*fh[idx_fh_F_ord2(iF + 1, jF, kF - 2, ex)]
+                             -    fh[idx_fh_F_ord2(iF + 2, jF, kF - 2, ex)] );
+
+                        const double t_km1 =
+                            ( fh[idx_fh_F_ord2(iF - 2, jF, kF - 1, ex)]
+                             -F8*fh[idx_fh_F_ord2(iF - 1, jF, kF - 1, ex)]
+                             +F8*fh[idx_fh_F_ord2(iF + 1, jF, kF - 1, ex)]
+                             -    fh[idx_fh_F_ord2(iF + 2, jF, kF - 1, ex)] );
+
+                        const double t_kp1 =
+                            ( fh[idx_fh_F_ord2(iF - 2, jF, kF + 1, ex)]
+                             -F8*fh[idx_fh_F_ord2(iF - 1, jF, kF + 1, ex)]
+                             +F8*fh[idx_fh_F_ord2(iF + 1, jF, kF + 1, ex)]
+                             -    fh[idx_fh_F_ord2(iF + 2, jF, kF + 1, ex)] );
+
+                        const double t_kp2 =
+                            ( fh[idx_fh_F_ord2(iF - 2, jF, kF + 2, ex)]
+                             -F8*fh[idx_fh_F_ord2(iF - 1, jF, kF + 2, ex)]
+                             +F8*fh[idx_fh_F_ord2(iF + 1, jF, kF + 2, ex)]
+                             -    fh[idx_fh_F_ord2(iF + 2, jF, kF + 2, ex)] );
+
+                        fxz[p] = Fdxdz * ( t_km2 - F8 * t_km1 + F8 * t_kp1 - t_kp2 );
+                    }
+
+                    {
+                        const double t_km2 =
+                            ( fh[idx_fh_F_ord2(iF, jF - 2, kF - 2, ex)]
+                             -F8*fh[idx_fh_F_ord2(iF, jF - 1, kF - 2, ex)]
+                             +F8*fh[idx_fh_F_ord2(iF, jF + 1, kF - 2, ex)]
+                             -    fh[idx_fh_F_ord2(iF, jF + 2, kF - 2, ex)] );
+
+                        const double t_km1 =
+                            ( fh[idx_fh_F_ord2(iF, jF - 2, kF - 1, ex)]
+                             -F8*fh[idx_fh_F_ord2(iF, jF - 1, kF - 1, ex)]
+                             +F8*fh[idx_fh_F_ord2(iF, jF + 1, kF - 1, ex)]
+                             -    fh[idx_fh_F_ord2(iF, jF + 2, kF - 1, ex)] );
+
+                        const double t_kp1 =
+                            ( fh[idx_fh_F_ord2(iF, jF - 2, kF + 1, ex)]
+                             -F8*fh[idx_fh_F_ord2(iF, jF - 1, kF + 1, ex)]
+                             +F8*fh[idx_fh_F_ord2(iF, jF + 1, kF + 1, ex)]
+                             -    fh[idx_fh_F_ord2(iF, jF + 2, kF + 1, ex)] );
+
+                        const double t_kp2 =
+                            ( fh[idx_fh_F_ord2(iF, jF - 2, kF + 2, ex)]
+                             -F8*fh[idx_fh_F_ord2(iF, jF - 1, kF + 2, ex)]
+                             +F8*fh[idx_fh_F_ord2(iF, jF + 1, kF + 2, ex)]
+                             -    fh[idx_fh_F_ord2(iF, jF + 2, kF + 2, ex)] );
+
+                        fyz[p] = Fdydz * ( t_km2 - F8 * t_km1 + F8 * t_kp1 - t_kp2 );
+                    }
                 }
             }
         }

AMSS_NCKU_source/fmisc.f90

@@ -1115,6 +1115,147 @@ end subroutine d2dump
 !------------------------------------------------------------------------------
 ! Lagrangian polynomial interpolation
 !------------------------------------------------------------------------------
+#ifndef POLINT6_USE_BARYCENTRIC
+#define POLINT6_USE_BARYCENTRIC 1
+#endif
+
+!DIR$ ATTRIBUTES FORCEINLINE :: polint6_neville
+  subroutine polint6_neville(xa, ya, x, y, dy)
+  implicit none
+
+  real*8, dimension(6), intent(in) :: xa, ya
+  real*8, intent(in) :: x
+  real*8, intent(out) :: y, dy
+
+  integer :: i, m, ns, n_m
+  real*8, dimension(6) :: c, d, ho
+  real*8 :: dif, dift, hp, h, den_val
+
+  c = ya
+  d = ya
+  ho = xa - x
+
+  ns = 1
+  dif = abs(x - xa(1))
+
+  do i = 2, 6
+    dift = abs(x - xa(i))
+    if (dift < dif) then
+      ns = i
+      dif = dift
+    end if
+  end do
+
+  y = ya(ns)
+  ns = ns - 1
+
+  do m = 1, 5
+    n_m = 6 - m
+    do i = 1, n_m
+      hp = ho(i)
+      h  = ho(i+m)
+      den_val = hp - h
+
+      if (den_val == 0.0d0) then
+        write(*,*) 'failure in polint for point',x
+        write(*,*) 'with input points: ',xa
+        stop
+      end if
+
+      den_val = (c(i+1) - d(i)) / den_val
+
+      d(i) = h * den_val
+      c(i) = hp * den_val
+    end do
+
+    if (2 * ns < n_m) then
+      dy = c(ns + 1)
+    else
+      dy = d(ns)
+      ns = ns - 1
+    end if
+    y = y + dy
+  end do
+
+  return
+  end subroutine polint6_neville
+
+!DIR$ ATTRIBUTES FORCEINLINE :: polint6_barycentric
+  subroutine polint6_barycentric(xa, ya, x, y, dy)
+  implicit none
+
+  real*8, dimension(6), intent(in) :: xa, ya
+  real*8, intent(in) :: x
+  real*8, intent(out) :: y, dy
+
+  integer :: i, j
+  logical :: is_uniform
+  real*8, dimension(6) :: lambda
+  real*8 :: dx, den_i, term, num, den, step, tol
+  real*8, parameter :: c_uniform(6) = (/ -1.d0, 5.d0, -10.d0, 10.d0, -5.d0, 1.d0 /)
+
+  do i = 1, 6
+    if (x == xa(i)) then
+      y = ya(i)
+      dy = 0.d0
+      return
+    end if
+  end do
+
+  step = xa(2) - xa(1)
+  is_uniform = (step /= 0.d0)
+  if (is_uniform) then
+    tol = 64.d0 * epsilon(1.d0) * max(1.d0, abs(step))
+    do i = 3, 6
+      if (abs((xa(i) - xa(i-1)) - step) > tol) then
+        is_uniform = .false.
+        exit
+      end if
+    end do
+  end if
+
+  if (is_uniform) then
+    num = 0.d0
+    den = 0.d0
+    do i = 1, 6
+      term = c_uniform(i) / (x - xa(i))
+      num = num + term * ya(i)
+      den = den + term
+    end do
+    y = num / den
+    dy = 0.d0
+    return
+  end if
+
+  do i = 1, 6
+    den_i = 1.d0
+    do j = 1, 6
+      if (j /= i) then
+        dx = xa(i) - xa(j)
+        if (dx == 0.0d0) then
+          write(*,*) 'failure in polint for point',x
+          write(*,*) 'with input points: ',xa
+          stop
+        end if
+        den_i = den_i * dx
+      end if
+    end do
+    lambda(i) = 1.d0 / den_i
+  end do
+
+  num = 0.d0
+  den = 0.d0
+  do i = 1, 6
+    term = lambda(i) / (x - xa(i))
+    num = num + term * ya(i)
+    den = den + term
+  end do
+
+  y = num / den
+  dy = 0.d0
+
+  return
+  end subroutine polint6_barycentric
 
 !DIR$ ATTRIBUTES FORCEINLINE :: polint
   subroutine polint(xa, ya, x, y, dy, ordn)
@@ -1129,6 +1270,15 @@ end subroutine d2dump
   real*8, dimension(ordn) :: c, d, ho
   real*8 :: dif, dift, hp, h, den_val
 
+  if (ordn == 6) then
+#if POLINT6_USE_BARYCENTRIC
+    call polint6_barycentric(xa, ya, x, y, dy)
+#else
+    call polint6_neville(xa, ya, x, y, dy)
+#endif
+    return
+  end if
+
   c = ya
   d = ya
   ho = xa - x
@@ -1177,6 +1327,70 @@ end subroutine d2dump
 
   return
   end subroutine polint
+
+  subroutine polint0(xa, ya, y, ordn)
+  ! Lagrange interpolation at x=0, O(n) direct formula
+  implicit none
+  integer, intent(in) :: ordn
+  real*8, dimension(ordn), intent(in) :: xa, ya
+  real*8, intent(out) :: y
+
+  integer :: j, k
+  real*8 :: wj
+
+  y = 0.d0
+  do j = 1, ordn
+    wj = 1.d0
+    do k = 1, ordn
+      if (k .ne. j) then
+        wj = wj * xa(k) / (xa(k) - xa(j))
+      endif
+    enddo
+    y = y + wj * ya(j)
+  enddo
+
+  return
+  end subroutine polint0
+!------------------------------------------------------------------------------
+!
+! interpolation in 2 dimensions, follow yx order
+!
+!------------------------------------------------------------------------------
+!------------------------------------------------------------------------------
+! Compute Lagrange interpolation basis weights for one target point.
+!------------------------------------------------------------------------------
+!DIR$ ATTRIBUTES FORCEINLINE :: polint_lagrange_weights
+  subroutine polint_lagrange_weights(xa, x, w, ordn)
+  implicit none
+
+  integer, intent(in) :: ordn
+  real*8, dimension(1:ordn), intent(in) :: xa
+  real*8, intent(in) :: x
+  real*8, dimension(1:ordn), intent(out) :: w
+
+  integer :: i, j
+  real*8 :: num, den, dx
+
+  do i = 1, ordn
+    num = 1.d0
+    den = 1.d0
+    do j = 1, ordn
+      if (j /= i) then
+        dx = xa(i) - xa(j)
+        if (dx == 0.0d0) then
+          write(*,*) 'failure in polint for point',x
+          write(*,*) 'with input points: ',xa
+          stop
+        end if
+        num = num * (x - xa(j))
+        den = den * dx
+      end if
+    end do
+    w(i) = num / den
+  end do
+
+  return
+  end subroutine polint_lagrange_weights
 !------------------------------------------------------------------------------
 !
 ! interpolation in 2 dimensions, follow yx order
@@ -1248,19 +1462,26 @@ end subroutine d2dump
   end do
   call polint(x1a,ymtmp,x1,y,dy,ordn)
 #else
-  integer  :: j, k
-  real*8, dimension(ordn,ordn) :: yatmp
+  integer  :: i, j, k
+  real*8, dimension(ordn) :: w1, w2
   real*8, dimension(ordn) :: ymtmp
-  real*8 :: dy_temp
+  real*8 :: yx_sum, x_sum
 
-  do k=1,ordn
-    do j=1,ordn
-      call polint(x1a, ya(:,j,k), x1, yatmp(j,k), dy_temp, ordn)
+  call polint_lagrange_weights(x1a, x1, w1, ordn)
+  call polint_lagrange_weights(x2a, x2, w2, ordn)
+
+  do k = 1, ordn
+    yx_sum = 0.d0
+    do j = 1, ordn
+      x_sum = 0.d0
+      do i = 1, ordn
+        x_sum = x_sum + w1(i) * ya(i,j,k)
+      end do
+      yx_sum = yx_sum + w2(j) * x_sum
     end do
+    ymtmp(k) = yx_sum
   end do
-  do k=1,ordn
-    call polint(x2a, yatmp(:,k), x2, ymtmp(k), dy_temp, ordn)
-  end do
+
   call polint(x3a, ymtmp, x3, y, dy, ordn)
 #endif
 
@@ -1609,8 +1830,11 @@ deallocate(f_flat)
 ! f=3/8*f_1 + 3/4*f_2 - 1/8*f_3
 
   real*8,parameter::C1=3.d0/8.d0,C2=3.d0/4.d0,C3=-1.d0/8.d0
+  integer :: i,j,k
 
-  fout = C1*f1+C2*f2+C3*f3
+  do concurrent (k=1:ext(3), j=1:ext(2), i=1:ext(1))
+    fout(i,j,k) = C1*f1(i,j,k)+C2*f2(i,j,k)+C3*f3(i,j,k)
+  end do
 
   return
 

AMSS_NCKU_source/makefile

@@ -2,6 +2,12 @@
 
 include makefile.inc
 
+## polint(ordn=6) kernel selector:
+##   1 (default): barycentric fast path
+##   0          : fallback to Neville path
+POLINT6_USE_BARY ?= 1
+POLINT6_FLAG = -DPOLINT6_USE_BARYCENTRIC=$(POLINT6_USE_BARY)
+
 ## ABE build flags selected by PGO_MODE (set in makefile.inc, default: opt)
 ##   make                        -> opt  (PGO-guided, maximum performance)
 ##   make PGO_MODE=instrument    -> instrument (Phase 1: collect fresh profile data)
@@ -12,7 +18,7 @@ ifeq ($(PGO_MODE),instrument)
 CXXAPPFLAGS = -O3 -xHost -fma -fprofile-instr-generate -ipo \
               -Dfortran3 -Dnewc -I${MKLROOT}/include $(INTERP_LB_FLAGS)
 f90appflags = -O3 -xHost -fma -fprofile-instr-generate -ipo \
-              -align array64byte -fpp -I${MKLROOT}/include
+              -align array64byte -fpp -I${MKLROOT}/include $(POLINT6_FLAG)
 else
 ## opt (default): maximum performance with PGO profile data
 CXXAPPFLAGS = -O3 -xHost -fp-model fast=2 -fma -ipo \
@@ -20,7 +26,7 @@ CXXAPPFLAGS = -O3 -xHost -fp-model fast=2 -fma -ipo \
               -Dfortran3 -Dnewc -I${MKLROOT}/include $(INTERP_LB_FLAGS)
 f90appflags = -O3 -xHost -fp-model fast=2 -fma -ipo \
               -fprofile-instr-use=$(PROFDATA) \
-              -align array64byte -fpp -I${MKLROOT}/include
+              -align array64byte -fpp -I${MKLROOT}/include $(POLINT6_FLAG)
 endif
 
 .SUFFIXES: .o .f90 .C .for .cu
@@ -37,6 +43,10 @@ endif
 .cu.o:
 	$(Cu) $(CUDA_APP_FLAGS) -c $< -o $@ $(CUDA_LIB_PATH)
 
+# CUDA rewrite of BSSN RHS (drop-in replacement for bssn_rhs_c + stencil helpers)
+bssn_rhs_cuda.o: bssn_rhs_cuda.cu macrodef.h
+	$(Cu) $(CUDA_APP_FLAGS) -c $< -o $@ $(CUDA_LIB_PATH)
+
 # C rewrite of BSSN RHS kernel and helpers
 bssn_rhs_c.o: bssn_rhs_c.C
 	${CXX} $(CXXAPPFLAGS) -c $< $(filein) -o $@
@@ -79,6 +89,18 @@ else
 CFILES = bssn_rhs_c.o fderivs_c.o fdderivs_c.o kodiss_c.o lopsided_c.o
 endif
 
+# CUDA rewrite: bssn_rhs_cuda.o replaces all CFILES (stencils are built-in)
+CFILES_CUDA = bssn_rhs_cuda.o
+
+## RK4 kernel switch (independent from USE_CXX_KERNELS)
+ifeq ($(USE_CXX_RK4),1)
+CFILES += rungekutta4_rout_c.o
+CFILES_CUDA += rungekutta4_rout_c.o
+RK4_F90_OBJ =
+else
+RK4_F90_OBJ = rungekutta4_rout.o
+endif
+
 C++FILES = ABE.o Ansorg.o Block.o misc.o monitor.o Parallel.o MPatch.o var.o\
            cgh.o bssn_class.o surface_integral.o ShellPatch.o\
 	   bssnEScalar_class.o perf.o Z4c_class.o NullShellPatch.o\
@@ -96,7 +118,7 @@ C++FILES_GPU = ABE.o Ansorg.o Block.o misc.o monitor.o Parallel.o MPatch.o var.o
 
 F90FILES_BASE = enforce_algebra.o fmisc.o initial_puncture.o prolongrestrict.o\
 	   prolongrestrict_cell.o prolongrestrict_vertex.o\
-	   rungekutta4_rout.o diff_new.o kodiss.o kodiss_sh.o\
+	   $(RK4_F90_OBJ) diff_new.o kodiss.o kodiss_sh.o\
 	   lopsidediff.o sommerfeld_rout.o getnp4.o diff_new_sh.o\
 	   shellfunctions.o bssn_rhs_ss.o Set_Rho_ADM.o\
            getnp4EScalar.o bssnEScalar_rhs.o bssn_constraint.o ricci_gamma.o\
@@ -159,8 +181,11 @@ $(CUDAFILES): bssn_gpu.h gpu_mem.h gpu_rhsSS_mem.h
 misc.o : zbesh.o
 
 # projects
-ABE: $(C++FILES) $(CFILES) $(F90FILES) $(F77FILES) $(AHFDOBJS)
-	$(CLINKER) $(CXXAPPFLAGS) -o $@ $(C++FILES) $(CFILES) $(F90FILES) $(F77FILES) $(AHFDOBJS) $(LDLIBS)
+ABE: $(C++FILES) $(CFILES_CUDA) $(F90FILES) $(F77FILES) $(AHFDOBJS)
+	$(CLINKER) $(CXXAPPFLAGS) -o $@ $(C++FILES) $(CFILES_CUDA) $(F90FILES) $(F77FILES) $(AHFDOBJS) $(LDLIBS) -lcudart $(CUDA_LIB_PATH)
+
+ABE_CUDA: $(C++FILES) $(CFILES_CUDA) $(F90FILES) $(F77FILES) $(AHFDOBJS)
+	$(CLINKER) $(CXXAPPFLAGS) -o $@ $(C++FILES) $(CFILES_CUDA) $(F90FILES) $(F77FILES) $(AHFDOBJS) $(LDLIBS) -lcudart $(CUDA_LIB_PATH)
 
 ABEGPU: $(C++FILES_GPU) $(CFILES) $(F90FILES) $(F77FILES) $(AHFDOBJS) $(CUDAFILES)
 	$(CLINKER) $(CXXAPPFLAGS) -o $@ $(C++FILES_GPU) $(CFILES) $(F90FILES) $(F77FILES) $(AHFDOBJS) $(CUDAFILES) $(LDLIBS)
@@ -169,4 +194,4 @@ TwoPunctureABE: $(TwoPunctureFILES)
 	$(CLINKER) $(TP_OPTFLAGS) -qopenmp -o $@ $(TwoPunctureFILES) $(LDLIBS)
 
 clean:
-	rm *.o ABE ABEGPU TwoPunctureABE make.log -f
+	rm *.o ABE ABE_CUDA ABEGPU TwoPunctureABE make.log -f

AMSS_NCKU_source/makefile.inc

@@ -10,6 +10,20 @@ filein  = -I/usr/include/ -I${MKLROOT}/include
 ## Added -lifcore for Intel Fortran runtime and -limf for Intel math library
 LDLIBS  = -L${MKLROOT}/lib -lmkl_intel_lp64 -lmkl_sequential -lmkl_core -lifcore -limf -lpthread -lm -ldl -liomp5
 
+## Memory allocator switch
+##   1 (default) : link Intel oneTBB allocator (libtbbmalloc)
+##   0           : use system default allocator (ptmalloc)
+USE_TBBMALLOC ?= 1
+TBBMALLOC_SO ?= /home/intel/oneapi/2025.3/lib/libtbbmalloc.so
+ifneq ($(wildcard $(TBBMALLOC_SO)),)
+TBBMALLOC_LIBS = -Wl,--no-as-needed $(TBBMALLOC_SO) -Wl,--as-needed
+else
+TBBMALLOC_LIBS = -Wl,--no-as-needed -ltbbmalloc -Wl,--as-needed
+endif
+ifeq ($(USE_TBBMALLOC),1)
+LDLIBS := $(TBBMALLOC_LIBS) $(LDLIBS)
+endif
+
 ## PGO build mode switch (ABE only; TwoPunctureABE always uses opt flags)
 ##   opt        : (default) maximum performance with PGO profile-guided optimization
 ##   instrument : PGO Phase 1 instrumentation to collect fresh profile data
@@ -33,6 +47,12 @@ endif
 ##   1 (default) : use C++ rewrite of bssn_rhs and helper kernels (faster)
 ##   0           : fall back to original Fortran kernels
 USE_CXX_KERNELS ?= 1
+
+## RK4 kernel implementation switch
+##   1 (default) : use C/C++ rewrite of rungekutta4_rout (for optimization experiments)
+##   0           : use original Fortran rungekutta4_rout.o
+USE_CXX_RK4 ?= 1
+
 f90          = ifx
 f77          = ifx
 CXX          = icpx
@@ -42,4 +62,4 @@ CLINKER      = mpiicpx
 Cu = nvcc
 CUDA_LIB_PATH = -L/usr/lib/cuda/lib64 -I/usr/include -I/usr/lib/cuda/include
 #CUDA_APP_FLAGS = -c -g -O3 --ptxas-options=-v -arch compute_13 -code compute_13,sm_13 -Dfortran3 -Dnewc
-CUDA_APP_FLAGS = -c -g -O3 --ptxas-options=-v -Dfortran3 -Dnewc
+CUDA_APP_FLAGS = -c -g -O3 --ptxas-options=-v -Dfortran3 -Dnewc -arch=sm_80

AMSS_NCKU_source/prolongrestrict_cell.f90

@@ -217,7 +217,6 @@
   real*8,dimension(2*ghost_width) :: X,Y,Z
   real*8, dimension(2*ghost_width,2*ghost_width) :: tmp2
   real*8, dimension(2*ghost_width) :: tmp1
-  real*8 :: ddy
   real*8,dimension(3) :: ccp
 
 #if (ghost_width == 2)
@@ -580,7 +579,7 @@
             tmp1(ghost_width-cxI(1)+cxB(1)  :ghost_width-cxI(1)+cxT(1)  ) = funf(cxB(1):cxT(1),j,k)
          endif
 
-         call polint(X,tmp1,0.d0,funf(i,j,k),ddy,2*ghost_width)
+         call polint0(X,tmp1,funf(i,j,k),2*ghost_width)
 
 ! for y direction
        elseif(sum(fg).eq.2.and.fg(2) .eq. 0.and. &
@@ -690,7 +689,7 @@
             tmp1(ghost_width-cxI(2)+cxB(2)  :ghost_width-cxI(2)+cxT(2)  ) = funf(i,cxB(2):cxT(2),k)
          endif
 
-         call polint(Y,tmp1,0.d0,funf(i,j,k),ddy,2*ghost_width)
+         call polint0(Y,tmp1,funf(i,j,k),2*ghost_width)
 
 ! for z direction
        elseif(sum(fg).eq.2.and.fg(3) .eq. 0.and. &
@@ -802,7 +801,7 @@
             tmp1(ghost_width-cxI(3)+cxB(3)  :ghost_width-cxI(3)+cxT(3)  ) = funf(i,j,cxB(3):cxT(3))
          endif
 
-         call polint(Z,tmp1,0.d0,funf(i,j,k),ddy,2*ghost_width)
+         call polint0(Z,tmp1,funf(i,j,k),2*ghost_width)
 
 #else
 
@@ -1934,18 +1933,35 @@
 ! when if=1 -> ic=0, this is different to vertex center grid 
   real*8, dimension(-2:extc(1),-2:extc(2),-2:extc(3))   :: funcc
   integer,dimension(3) :: cxI
-  integer :: i,j,k,ii,jj,kk
+  integer :: i,j,k,ii,jj,kk,px,py,pz
   real*8, dimension(6,6) :: tmp2
   real*8, dimension(6) :: tmp1
+  integer, dimension(extf(1)) :: cix
+  integer, dimension(extf(2)) :: ciy
+  integer, dimension(extf(3)) :: ciz
+  integer, dimension(extf(1)) :: pix
+  integer, dimension(extf(2)) :: piy
+  integer, dimension(extf(3)) :: piz
 
   real*8, parameter :: C1=7.7d1/8.192d3,C2=-6.93d2/8.192d3,C3=3.465d3/4.096d3
   real*8, parameter :: C6=6.3d1/8.192d3,C5=-4.95d2/8.192d3,C4=1.155d3/4.096d3
+  real*8, dimension(6,2), parameter :: WC = reshape((/&
+      C1,C2,C3,C4,C5,C6,&
+      C6,C5,C4,C3,C2,C1/), (/6,2/))
 
   integer::imini,imaxi,jmini,jmaxi,kmini,kmaxi
   integer::imino,imaxo,jmino,jmaxo,kmino,kmaxo
+  integer::maxcx,maxcy,maxcz
 
   real*8,dimension(3) :: CD,FD
-  
+  real*8 :: tmp_yz(extc(1), 6)      ! 存储整条 X 线上 6 个 Y 轴偏置的 Z 向插值结果
+  real*8 :: tmp_xyz_line(-2:extc(1))   ! 包含 X 向 6 点模板访问所需下界
+  real*8 :: v1, v2, v3, v4, v5, v6
+  integer :: ic, jc, kc, ix_offset,ix,iy,iz,jc_min,jc_max,ic_min,ic_max,kc_min,kc_max
+  integer :: i_lo, i_hi, j_lo, j_hi, k_lo, k_hi
+  logical :: need_full_symmetry
+  real*8 :: res_line
+  real*8 :: tmp_z_slab(-2:extc(1), -2:extc(2))  ! 包含 Y/X 向模板访问所需下界
   if(wei.ne.3)then
      write(*,*)"prolongrestrict.f90::prolong3: this routine only surport 3 dimension"
      write(*,*)"dim = ",wei
@@ -2020,145 +2036,140 @@
           return
   endif
 
-  call symmetry_bd(3,extc,func,funcc,SoA)
-     
-!~~~~~~> prolongation start...
+  do i = imino,imaxo
+     ii = i + lbf(1) - 1
+     cix(i) = ii/2 - lbc(1) + 1
+     if(ii/2*2 == ii)then
+        pix(i) = 1
+     else
+        pix(i) = 2
+     endif
+  enddo
+  do j = jmino,jmaxo
+     jj = j + lbf(2) - 1
+     ciy(j) = jj/2 - lbc(2) + 1
+     if(jj/2*2 == jj)then
+        piy(j) = 1
+     else
+        piy(j) = 2
+     endif
+  enddo
   do k = kmino,kmaxo
-   do j = jmino,jmaxo
-    do i = imino,imaxo
-       cxI(1) = i
-       cxI(2) = j
-       cxI(3) = k
-! change to coarse level reference
-!|---*--- ---*--- ---*--- ---*--- ---*--- ---*--- ---*--- ---*---| 
-!|=======x===============x===============x===============x=======|
-       cxI = (cxI+lbf-1)/2
-! change to array index      
-       cxI = cxI - lbc + 1
-
-       if(any(cxI+3 > extc)) write(*,*)"error in prolong"
-       ii=i+lbf(1)-1
-       jj=j+lbf(2)-1
-       kk=k+lbf(3)-1
-#if 0
-       if(ii/2*2==ii)then
-         if(jj/2*2==jj)then
-           if(kk/2*2==kk)then
-             tmp2= C1*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)-2)+&
-                   C2*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)-1)+&
-                   C3*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)  )+&
-                   C4*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)+1)+&
-                   C5*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)+2)+&
-                   C6*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)+3)
-             tmp1= C1*tmp2(:,1)+C2*tmp2(:,2)+C3*tmp2(:,3)+C4*tmp2(:,4)+C5*tmp2(:,5)+C6*tmp2(:,6)
-             funf(i,j,k)= C1*tmp1(1)+C2*tmp1(2)+C3*tmp1(3)+C4*tmp1(4)+C5*tmp1(5)+C6*tmp1(6)
-           else
-             tmp2= C6*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)-2)+&
-                   C5*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)-1)+&
-                   C4*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)  )+&
-                   C3*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)+1)+&
-                   C2*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)+2)+&
-                   C1*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)+3)
-             tmp1= C1*tmp2(:,1)+C2*tmp2(:,2)+C3*tmp2(:,3)+C4*tmp2(:,4)+C5*tmp2(:,5)+C6*tmp2(:,6)
-             funf(i,j,k)=  C1*tmp1(1)+C2*tmp1(2)+C3*tmp1(3)+C4*tmp1(4)+C5*tmp1(5)+C6*tmp1(6)
-           endif
-         else
-           if(kk/2*2==kk)then
-             tmp2= C1*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)-2)+&
-                   C2*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)-1)+&
-                   C3*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)  )+&
-                   C4*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)+1)+&
-                   C5*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)+2)+&
-                   C6*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)+3)
-             tmp1= C6*tmp2(:,1)+C5*tmp2(:,2)+C4*tmp2(:,3)+C3*tmp2(:,4)+C2*tmp2(:,5)+C1*tmp2(:,6)
-             funf(i,j,k)= C1*tmp1(1)+C2*tmp1(2)+C3*tmp1(3)+C4*tmp1(4)+C5*tmp1(5)+C6*tmp1(6)
-           else
-             tmp2= C6*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)-2)+&
-                   C5*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)-1)+&
-                   C4*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)  )+&
-                   C3*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)+1)+&
-                   C2*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)+2)+&
-                   C1*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)+3)
-             tmp1= C6*tmp2(:,1)+C5*tmp2(:,2)+C4*tmp2(:,3)+C3*tmp2(:,4)+C2*tmp2(:,5)+C1*tmp2(:,6)
-             funf(i,j,k)=  C1*tmp1(1)+C2*tmp1(2)+C3*tmp1(3)+C4*tmp1(4)+C5*tmp1(5)+C6*tmp1(6)
-           endif
-         endif
-       else
-         if(jj/2*2==jj)then
-           if(kk/2*2==kk)then               
-             tmp2= C1*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)-2)+&
-                   C2*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)-1)+&
-                   C3*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)  )+&
-                   C4*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)+1)+&
-                   C5*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)+2)+&
-                   C6*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)+3)
-             tmp1= C1*tmp2(:,1)+C2*tmp2(:,2)+C3*tmp2(:,3)+C4*tmp2(:,4)+C5*tmp2(:,5)+C6*tmp2(:,6)
-             funf(i,j,k)= C6*tmp1(1)+C5*tmp1(2)+C4*tmp1(3)+C3*tmp1(4)+C2*tmp1(5)+C1*tmp1(6)
-           else
-             tmp2= C6*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)-2)+&
-                   C5*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)-1)+&
-                   C4*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)  )+&
-                   C3*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)+1)+&
-                   C2*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)+2)+&
-                   C1*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)+3)
-             tmp1= C1*tmp2(:,1)+C2*tmp2(:,2)+C3*tmp2(:,3)+C4*tmp2(:,4)+C5*tmp2(:,5)+C6*tmp2(:,6)
-             funf(i,j,k)=  C6*tmp1(1)+C5*tmp1(2)+C4*tmp1(3)+C3*tmp1(4)+C2*tmp1(5)+C1*tmp1(6)
-           endif
-         else
-           if(kk/2*2==kk)then
-             tmp2= C1*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)-2)+&
-                   C2*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)-1)+&
-                   C3*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)  )+&
-                   C4*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)+1)+&
-                   C5*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)+2)+&
-                   C6*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)+3)
-             tmp1= C6*tmp2(:,1)+C5*tmp2(:,2)+C4*tmp2(:,3)+C3*tmp2(:,4)+C2*tmp2(:,5)+C1*tmp2(:,6)
-             funf(i,j,k)= C6*tmp1(1)+C5*tmp1(2)+C4*tmp1(3)+C3*tmp1(4)+C2*tmp1(5)+C1*tmp1(6)
-           else
-             tmp2= C6*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)-2)+&
-                   C5*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)-1)+&
-                   C4*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)  )+&
-                   C3*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)+1)+&
-                   C2*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)+2)+&
-                   C1*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)+3)
-             tmp1= C6*tmp2(:,1)+C5*tmp2(:,2)+C4*tmp2(:,3)+C3*tmp2(:,4)+C2*tmp2(:,5)+C1*tmp2(:,6)
-             funf(i,j,k)=  C6*tmp1(1)+C5*tmp1(2)+C4*tmp1(3)+C3*tmp1(4)+C2*tmp1(5)+C1*tmp1(6)
-           endif
-         endif
-       endif
-#else 
-       if(kk/2*2==kk)then
-             tmp2= C1*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)-2)+&
-                   C2*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)-1)+&
-                   C3*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)  )+&
-                   C4*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)+1)+&
-                   C5*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)+2)+&
-                   C6*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)+3)
-       else
-             tmp2= C6*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)-2)+&
-                   C5*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)-1)+&
-                   C4*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)  )+&
-                   C3*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)+1)+&
-                   C2*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)+2)+&
-                   C1*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)+3)
-       endif
-
-       if(jj/2*2==jj)then
-             tmp1= C1*tmp2(:,1)+C2*tmp2(:,2)+C3*tmp2(:,3)+C4*tmp2(:,4)+C5*tmp2(:,5)+C6*tmp2(:,6)
-       else
-             tmp1= C6*tmp2(:,1)+C5*tmp2(:,2)+C4*tmp2(:,3)+C3*tmp2(:,4)+C2*tmp2(:,5)+C1*tmp2(:,6)
-       endif
-
-       if(ii/2*2==ii)then
-             funf(i,j,k)= C1*tmp1(1)+C2*tmp1(2)+C3*tmp1(3)+C4*tmp1(4)+C5*tmp1(5)+C6*tmp1(6)
-       else
-             funf(i,j,k)= C6*tmp1(1)+C5*tmp1(2)+C4*tmp1(3)+C3*tmp1(4)+C2*tmp1(5)+C1*tmp1(6)
-       endif
-#endif
-    enddo
-   enddo
+     kk = k + lbf(3) - 1
+     ciz(k) = kk/2 - lbc(3) + 1
+     if(kk/2*2 == kk)then
+        piz(k) = 1
+     else
+        piz(k) = 2
+     endif
   enddo
 
+  ic_min = minval(cix(imino:imaxo))
+  ic_max = maxval(cix(imino:imaxo))
+  jc_min = minval(ciy(jmino:jmaxo))
+  jc_max = maxval(ciy(jmino:jmaxo))
+  kc_min = minval(ciz(kmino:kmaxo))
+  kc_max = maxval(ciz(kmino:kmaxo))
+
+  maxcx = ic_max
+  maxcy = jc_max
+  maxcz = kc_max
+  if(maxcx+3 > extc(1) .or. maxcy+3 > extc(2) .or. maxcz+3 > extc(3))then
+     write(*,*)"error in prolong"
+     return
+  endif
+
+  i_lo = ic_min - 2
+  i_hi = ic_max + 3
+  j_lo = jc_min - 2
+  j_hi = jc_max + 3
+  k_lo = kc_min - 2
+  k_hi = kc_max + 3
+  need_full_symmetry = (i_lo < 1) .or. (j_lo < 1) .or. (k_lo < 1)
+  if(need_full_symmetry)then
+     call symmetry_bd(3,extc,func,funcc,SoA)
+  else
+     funcc(i_lo:i_hi,j_lo:j_hi,k_lo:k_hi) = func(i_lo:i_hi,j_lo:j_hi,k_lo:k_hi)
+  endif
+
+     ! 对每个 k（pz, kc 固定）预计算 Z 向插值的 2D 切片
+
+do k = kmino, kmaxo
+    pz = piz(k); kc = ciz(k)
+    ! --- Pass 1: Z 方向，只算一次 ---
+    do iy = jc_min-2, jc_max+3   ! 仅需的 iy 范围（对应 jc-2:jc+3）
+        do ii = ic_min-2, ic_max+3  ! 仅需的 ii 范围（对应 cix-2:cix+3）
+            tmp_z_slab(ii, iy) = sum(WC(:,pz) * funcc(ii, iy, kc-2:kc+3))
+        end do
+    end do
+
+    do j = jmino, jmaxo
+        py = piy(j); jc = ciy(j)
+        ! --- Pass 2: Y 方向 ---
+        do ii = ic_min-2, ic_max+3
+            tmp_xyz_line(ii) = sum(WC(:,py) * tmp_z_slab(ii, jc-2:jc+3))
+        end do
+        ! --- Pass 3: X 方向 ---
+        do i = imino, imaxo
+            funf(i,j,k) = sum(WC(:,pix(i)) * tmp_xyz_line(cix(i)-2:cix(i)+3))
+        end do
+    end do
+end do
+
+!~~~~~~> prolongation start...
+#if 0
+ do k = kmino, kmaxo
+     pz = piz(k)
+     kc = ciz(k)
+
+     do j = jmino, jmaxo
+        py = piy(j)
+        jc = ciy(j)
+
+! --- 步骤 1 & 2 融合：分段处理 X 轴，提升 Cache 命中率 ---
+        ! 我们将 ii 循环逻辑重组，减少对 funcc 的跨行重复访问
+        do ii = 1, extc(1)
+           ! 1. 先做 Z 方向的 6 条线插值（针对当前的 ii 和当前的 6 个 iy）
+           ! 我们直接在这里把 Y 方向的加权也做了，省去 tmp_yz 数组
+           ! 这样 funcc 的数据读进来后立即完成所有维度的贡献，不再写回内存
+           
+           res_line = 0.0d0
+           do jj = 1, 6
+              iy = jc - 3 + jj
+              ! 这一行代码是核心：一次性完成 Z 插值并加上 Y 的权重
+              ! 编译器会把 WC(jj, py) 存在寄存器里
+              res_line = res_line + WC(jj, py) * ( &
+                         WC(1, pz) * funcc(ii, iy, kc-2) + &
+                         WC(2, pz) * funcc(ii, iy, kc-1) + &
+                         WC(3, pz) * funcc(ii, iy, kc  ) + &
+                         WC(4, pz) * funcc(ii, iy, kc+1) + &
+                         WC(5, pz) * funcc(ii, iy, kc+2) + &
+                         WC(6, pz) * funcc(ii, iy, kc+3) )
+           end do
+           tmp_xyz_line(ii) = res_line
+        end do
+
+
+
+
+        ! 3. 【降维：X 向】最后在最内层只处理 X 方向的 6 点加权
+        ! 此时每个点的计算量从原来的 200+ 次乘法降到了仅 6 次
+        do i = imino, imaxo
+           px = pix(i)
+           ic = cix(i)
+           
+           ! 直接从预计算好的 line 中读取连续的 6 个点
+           ! ic-2 到 ic+3 对应原始 6 点算子
+           funf(i,j,k) = WC(1,px)*tmp_xyz_line(ic-2) + &
+                         WC(2,px)*tmp_xyz_line(ic-1) + &
+                         WC(3,px)*tmp_xyz_line(ic  ) + &
+                         WC(4,px)*tmp_xyz_line(ic+1) + &
+                         WC(5,px)*tmp_xyz_line(ic+2) + &
+                         WC(6,px)*tmp_xyz_line(ic+3)
+        end do
+     end do
+  end do
+#endif
   return
 
   end subroutine prolong3
@@ -2358,6 +2369,13 @@
 
   real*8,dimension(3) :: CD,FD
 
+  real*8 :: tmp_xz_plane(-1:extf(1), 6)
+  real*8 :: tmp_x_line(-1:extf(1))
+  integer :: fi, fj, fk, ii, jj, kk
+  integer :: fi_min, fi_max, ii_lo, ii_hi
+  integer :: fj_min, fj_max, fk_min, fk_max, jj_lo, jj_hi, kk_lo, kk_hi
+  logical :: need_full_symmetry
+
   if(wei.ne.3)then
      write(*,*)"prolongrestrict.f90::restrict3: this routine only surport 3 dimension"
      write(*,*)"dim = ",wei
@@ -2436,9 +2454,86 @@
           stop
   endif
 
-  call symmetry_bd(2,extf,funf,funff,SoA)
+  ! 仅计算 X 向最终写回所需的窗口：
+  ! func(i,j,k) 只访问 tmp_x_line(fi-2:fi+3)
+  fi_min = 2*(imino + lbc(1) - 1) - 1 - lbf(1) + 1
+  fi_max = 2*(imaxo + lbc(1) - 1) - 1 - lbf(1) + 1
+  fj_min = 2*(jmino + lbc(2) - 1) - 1 - lbf(2) + 1
+  fj_max = 2*(jmaxo + lbc(2) - 1) - 1 - lbf(2) + 1
+  fk_min = 2*(kmino + lbc(3) - 1) - 1 - lbf(3) + 1
+  fk_max = 2*(kmaxo + lbc(3) - 1) - 1 - lbf(3) + 1
+  ii_lo = fi_min - 2
+  ii_hi = fi_max + 3
+  jj_lo = fj_min - 2
+  jj_hi = fj_max + 3
+  kk_lo = fk_min - 2
+  kk_hi = fk_max + 3
+  if(ii_lo < -1 .or. ii_hi > extf(1) .or. &
+     jj_lo < -1 .or. jj_hi > extf(2) .or. &
+     kk_lo < -1 .or. kk_hi > extf(3))then
+      write(*,*)"restrict3: invalid stencil window"
+      write(*,*)"ii=",ii_lo,ii_hi," jj=",jj_lo,jj_hi," kk=",kk_lo,kk_hi
+      write(*,*)"extf=",extf
+      stop
+  endif
+  need_full_symmetry = (ii_lo < 1) .or. (jj_lo < 1) .or. (kk_lo < 1)
+  if(need_full_symmetry)then
+      call symmetry_bd(2,extf,funf,funff,SoA)
+  else
+      funff(ii_lo:ii_hi,jj_lo:jj_hi,kk_lo:kk_hi) = funf(ii_lo:ii_hi,jj_lo:jj_hi,kk_lo:kk_hi)
+  endif
 
 !~~~~~~> restriction start...
+do k = kmino, kmaxo
+    fk = 2*(k + lbc(3) - 1) - 1 - lbf(3) + 1
+
+    do j = jmino, jmaxo
+        fj = 2*(j + lbc(2) - 1) - 1 - lbf(2) + 1
+        
+        ! 优化点 1: 显式展开 Z 方向计算，减少循环开销
+        ! 确保 ii 循环是最内层且连续访问
+        !DIR$ VECTOR ALWAYS
+        do ii = ii_lo, ii_hi
+            ! 预计算当前 j 对应的 6 行在 Z 方向的压缩结果
+            ! 这里直接硬编码 jj 的偏移，彻底消除一层循环
+            tmp_xz_plane(ii, 1) = C1*(funff(ii,fj-2,fk-2)+funff(ii,fj-2,fk+3)) + &
+                                  C2*(funff(ii,fj-2,fk-1)+funff(ii,fj-2,fk+2)) + &
+                                  C3*(funff(ii,fj-2,fk  )+funff(ii,fj-2,fk+1))
+            tmp_xz_plane(ii, 2) = C1*(funff(ii,fj-1,fk-2)+funff(ii,fj-1,fk+3)) + &
+                                  C2*(funff(ii,fj-1,fk-1)+funff(ii,fj-1,fk+2)) + &
+                                  C3*(funff(ii,fj-1,fk  )+funff(ii,fj-1,fk+1))
+            tmp_xz_plane(ii, 3) = C1*(funff(ii,fj  ,fk-2)+funff(ii,fj  ,fk+3)) + &
+                                  C2*(funff(ii,fj  ,fk-1)+funff(ii,fj  ,fk+2)) + &
+                                  C3*(funff(ii,fj  ,fk  )+funff(ii,fj  ,fk+1))
+            tmp_xz_plane(ii, 4) = C1*(funff(ii,fj+1,fk-2)+funff(ii,fj+1,fk+3)) + &
+                                  C2*(funff(ii,fj+1,fk-1)+funff(ii,fj+1,fk+2)) + &
+                                  C3*(funff(ii,fj+1,fk  )+funff(ii,fj+1,fk+1))
+            tmp_xz_plane(ii, 5) = C1*(funff(ii,fj+2,fk-2)+funff(ii,fj+2,fk+3)) + &
+                                  C2*(funff(ii,fj+2,fk-1)+funff(ii,fj+2,fk+2)) + &
+                                  C3*(funff(ii,fj+2,fk  )+funff(ii,fj+2,fk+1))
+            tmp_xz_plane(ii, 6) = C1*(funff(ii,fj+3,fk-2)+funff(ii,fj+3,fk+3)) + &
+                                  C2*(funff(ii,fj+3,fk-1)+funff(ii,fj+3,fk+2)) + &
+                                  C3*(funff(ii,fj+3,fk  )+funff(ii,fj+3,fk+1))
+        end do
+
+        ! 优化点 2: 同样向量化 Y 方向压缩
+        !DIR$ VECTOR ALWAYS
+        do ii = ii_lo, ii_hi
+            tmp_x_line(ii) = C1*(tmp_xz_plane(ii, 1) + tmp_xz_plane(ii, 6)) + &
+                            C2*(tmp_xz_plane(ii, 2) + tmp_xz_plane(ii, 5)) + &
+                            C3*(tmp_xz_plane(ii, 3) + tmp_xz_plane(ii, 4))
+        end do
+
+        ! 优化点 3: 最终写入，利用已经缓存在 tmp_x_line 的数据
+        do i = imino, imaxo
+            fi = 2*(i + lbc(1) - 1) - 1 - lbf(1) + 1
+            func(i, j, k) = C1*(tmp_x_line(fi-2) + tmp_x_line(fi+3)) + &
+                            C2*(tmp_x_line(fi-1) + tmp_x_line(fi+2)) + &
+                            C3*(tmp_x_line(fi  ) + tmp_x_line(fi+1))
+        end do
+    end do
+end do
+#if 0
   do k = kmino,kmaxo
    do j = jmino,jmaxo
     do i = imino,imaxo
@@ -2462,7 +2557,7 @@
     enddo
    enddo
   enddo
-  
+#endif
   return
 
   end subroutine restrict3

AMSS_NCKU_source/prolongrestrict_vertex.f90

@@ -217,7 +217,6 @@
   real*8,dimension(2*ghost_width) :: X,Y,Z
   real*8, dimension(2*ghost_width,2*ghost_width) :: tmp2
   real*8, dimension(2*ghost_width) :: tmp1
-  real*8 :: ddy
 
 #if (ghost_width == 2)
   real*8, parameter :: C1=-1.d0/16,C2=9.d0/16
@@ -470,7 +469,7 @@
 
          tmp1(cxB(1)+ghost_width-i+1:cxT(1)+ghost_width-i+1) = fh(cxB(1):cxT(1),j,k)
 
-         call polint(X,tmp1,0.d0,funf(i,j,k),ddy,2*ghost_width)
+         call polint0(X,tmp1,funf(i,j,k),2*ghost_width)
 
 ! for y direction
        elseif (fg(2) .eq. 0)then
@@ -529,7 +528,7 @@
 
          tmp1(cxB(2)+ghost_width-j+1:cxT(2)+ghost_width-j+1) = fh(i,cxB(2):cxT(2),k)
 
-         call polint(Y,tmp1,0.d0,funf(i,j,k),ddy,2*ghost_width)
+         call polint0(Y,tmp1,funf(i,j,k),2*ghost_width)
 
 ! for z direction
        else
@@ -588,7 +587,7 @@
 
          tmp1(cxB(3)+ghost_width-k+1:cxT(3)+ghost_width-k+1) = fh(i,j,cxB(3):cxT(3))
 
-         call polint(Z,tmp1,0.d0,funf(i,j,k),ddy,2*ghost_width)
+         call polint0(Z,tmp1,funf(i,j,k),2*ghost_width)
 
        endif
       

AMSS_NCKU_source/rungekutta4_rout_c.C

@@ -0,0 +1,155 @@
+#include "rungekutta4_rout.h"
+#include <cstdio>
+#include <cstdlib>
+#include <cstddef>
+#include <immintrin.h>
+
+namespace {
+
+inline void rk4_stage0(std::size_t n,
+                       const double *__restrict f0,
+                       const double *__restrict frhs,
+                       double *__restrict f1,
+                       double c) {
+    std::size_t i = 0;
+#if defined(__AVX512F__)
+    const __m512d vc = _mm512_set1_pd(c);
+    for (; i + 7 < n; i += 8) {
+        const __m512d v0 = _mm512_loadu_pd(f0 + i);
+        const __m512d vr = _mm512_loadu_pd(frhs + i);
+        _mm512_storeu_pd(f1 + i, _mm512_fmadd_pd(vc, vr, v0));
+    }
+#elif defined(__AVX2__)
+    const __m256d vc = _mm256_set1_pd(c);
+    for (; i + 3 < n; i += 4) {
+        const __m256d v0 = _mm256_loadu_pd(f0 + i);
+        const __m256d vr = _mm256_loadu_pd(frhs + i);
+        _mm256_storeu_pd(f1 + i, _mm256_fmadd_pd(vc, vr, v0));
+    }
+#endif
+#pragma ivdep
+    for (; i < n; ++i) {
+        f1[i] = f0[i] + c * frhs[i];
+    }
+}
+
+inline void rk4_rhs_accum(std::size_t n,
+                          const double *__restrict f1,
+                          double *__restrict frhs) {
+    std::size_t i = 0;
+#if defined(__AVX512F__)
+    const __m512d v2 = _mm512_set1_pd(2.0);
+    for (; i + 7 < n; i += 8) {
+        const __m512d v1 = _mm512_loadu_pd(f1 + i);
+        const __m512d vrhs = _mm512_loadu_pd(frhs + i);
+        _mm512_storeu_pd(frhs + i, _mm512_fmadd_pd(v2, v1, vrhs));
+    }
+#elif defined(__AVX2__)
+    const __m256d v2 = _mm256_set1_pd(2.0);
+    for (; i + 3 < n; i += 4) {
+        const __m256d v1 = _mm256_loadu_pd(f1 + i);
+        const __m256d vrhs = _mm256_loadu_pd(frhs + i);
+        _mm256_storeu_pd(frhs + i, _mm256_fmadd_pd(v2, v1, vrhs));
+    }
+#endif
+#pragma ivdep
+    for (; i < n; ++i) {
+        frhs[i] = frhs[i] + 2.0 * f1[i];
+    }
+}
+
+inline void rk4_f1_from_f0_f1(std::size_t n,
+                              const double *__restrict f0,
+                              double *__restrict f1,
+                              double c) {
+    std::size_t i = 0;
+#if defined(__AVX512F__)
+    const __m512d vc = _mm512_set1_pd(c);
+    for (; i + 7 < n; i += 8) {
+        const __m512d v0 = _mm512_loadu_pd(f0 + i);
+        const __m512d v1 = _mm512_loadu_pd(f1 + i);
+        _mm512_storeu_pd(f1 + i, _mm512_fmadd_pd(vc, v1, v0));
+    }
+#elif defined(__AVX2__)
+    const __m256d vc = _mm256_set1_pd(c);
+    for (; i + 3 < n; i += 4) {
+        const __m256d v0 = _mm256_loadu_pd(f0 + i);
+        const __m256d v1 = _mm256_loadu_pd(f1 + i);
+        _mm256_storeu_pd(f1 + i, _mm256_fmadd_pd(vc, v1, v0));
+    }
+#endif
+#pragma ivdep
+    for (; i < n; ++i) {
+        f1[i] = f0[i] + c * f1[i];
+    }
+}
+
+inline void rk4_stage3(std::size_t n,
+                       const double *__restrict f0,
+                       double *__restrict f1,
+                       const double *__restrict frhs,
+                       double c) {
+    std::size_t i = 0;
+#if defined(__AVX512F__)
+    const __m512d vc = _mm512_set1_pd(c);
+    for (; i + 7 < n; i += 8) {
+        const __m512d v0 = _mm512_loadu_pd(f0 + i);
+        const __m512d v1 = _mm512_loadu_pd(f1 + i);
+        const __m512d vr = _mm512_loadu_pd(frhs + i);
+        _mm512_storeu_pd(f1 + i, _mm512_fmadd_pd(vc, _mm512_add_pd(v1, vr), v0));
+    }
+#elif defined(__AVX2__)
+    const __m256d vc = _mm256_set1_pd(c);
+    for (; i + 3 < n; i += 4) {
+        const __m256d v0 = _mm256_loadu_pd(f0 + i);
+        const __m256d v1 = _mm256_loadu_pd(f1 + i);
+        const __m256d vr = _mm256_loadu_pd(frhs + i);
+        _mm256_storeu_pd(f1 + i, _mm256_fmadd_pd(vc, _mm256_add_pd(v1, vr), v0));
+    }
+#endif
+#pragma ivdep
+    for (; i < n; ++i) {
+        f1[i] = f0[i] + c * (f1[i] + frhs[i]);
+    }
+}
+
+} // namespace
+
+extern "C" {
+
+int f_rungekutta4_rout(int *ex, double &dT,
+                       double *f0, double *f1, double *f_rhs,
+                       int &RK4) {
+    const std::size_t n = static_cast<std::size_t>(ex[0]) *
+                          static_cast<std::size_t>(ex[1]) *
+                          static_cast<std::size_t>(ex[2]);
+    const double *const __restrict f0r = f0;
+    double *const __restrict f1r = f1;
+    double *const __restrict frhs = f_rhs;
+
+    if (__builtin_expect(static_cast<unsigned>(RK4) > 3u, 0)) {
+        std::fprintf(stderr, "rungekutta4_rout_c: invalid RK4 stage %d\n", RK4);
+        std::abort();
+    }
+
+    switch (RK4) {
+    case 0:
+        rk4_stage0(n, f0r, frhs, f1r, 0.5 * dT);
+        break;
+    case 1:
+        rk4_rhs_accum(n, f1r, frhs);
+        rk4_f1_from_f0_f1(n, f0r, f1r, 0.5 * dT);
+        break;
+    case 2:
+        rk4_rhs_accum(n, f1r, frhs);
+        rk4_f1_from_f0_f1(n, f0r, f1r, dT);
+        break;
+    default:
+        rk4_stage3(n, f0r, f1r, frhs, (1.0 / 6.0) * dT);
+        break;
+    }
+
+    return 0;
+}
+
+} // extern "C"

AMSS_NCKU_source/share_func.h

@@ -5,6 +5,7 @@
 #include <stddef.h>
 #include <math.h>
 #include <stdio.h>
+#include <string.h>
 /* 主网格：0-based -> 1D */
 static inline size_t idx_ex(int i0, int j0, int k0, const int ex[3]) {
     const int ex1 = ex[0], ex2 = ex[1];
@@ -87,60 +88,159 @@ static inline size_t idx_funcc_F(int iF, int jF, int kF, int ord, const int extc
  *   funcc(:,:,-i) = funcc(:,:,i+1)*SoA(3)
  * enddo
  */
+static inline void symmetry_bd_impl(int ord,
+                 int shift,
+                 const int extc[3],
+                 const double *__restrict func,
+                 double *__restrict funcc,
+                 const double SoA[3])
+{
+    const int extc1 = extc[0], extc2 = extc[1], extc3 = extc[2];
+    const int nx = extc1 + ord;
+    const int ny = extc2 + ord;
+
+    const size_t snx = (size_t)nx;
+    const size_t splane = (size_t)nx * (size_t)ny;
+    const size_t interior_i = (size_t)shift + 1u;          /* iF = 1 */
+    const size_t interior_j = ((size_t)shift + 1u) * snx;  /* jF = 1 */
+    const size_t interior_k = ((size_t)shift + 1u) * splane; /* kF = 1 */
+    const size_t interior0 = interior_k + interior_j + interior_i;
+
+    /* 1) funcc(1:extc1,1:extc2,1:extc3) = func */
+    for (int k0 = 0; k0 < extc3; ++k0) {
+        const double *src_k = func + (size_t)k0 * (size_t)extc2 * (size_t)extc1;
+        const size_t dst_k0 = interior0 + (size_t)k0 * splane;
+        for (int j0 = 0; j0 < extc2; ++j0) {
+            const double *src = src_k + (size_t)j0 * (size_t)extc1;
+            double *dst = funcc + dst_k0 + (size_t)j0 * snx;
+            memcpy(dst, src, (size_t)extc1 * sizeof(double));
+        }
+    }
+
+    /* 2) funcc(-i,1:extc2,1:extc3) = funcc(i+1,1:extc2,1:extc3)*SoA(1) */
+    const double s1 = SoA[0];
+    if (s1 == 1.0) {
+        for (int ii = 0; ii < ord; ++ii) {
+            const size_t dst_i = (size_t)(shift - ii);
+            const size_t src_i = (size_t)(shift + ii + 1);
+            for (int k0 = 0; k0 < extc3; ++k0) {
+                const size_t kbase = interior_k + (size_t)k0 * splane + interior_j;
+                for (int j0 = 0; j0 < extc2; ++j0) {
+                    const size_t off = kbase + (size_t)j0 * snx;
+                    funcc[off + dst_i] = funcc[off + src_i];
+                }
+            }
+        }
+    } else if (s1 == -1.0) {
+        for (int ii = 0; ii < ord; ++ii) {
+            const size_t dst_i = (size_t)(shift - ii);
+            const size_t src_i = (size_t)(shift + ii + 1);
+            for (int k0 = 0; k0 < extc3; ++k0) {
+                const size_t kbase = interior_k + (size_t)k0 * splane + interior_j;
+                for (int j0 = 0; j0 < extc2; ++j0) {
+                    const size_t off = kbase + (size_t)j0 * snx;
+                    funcc[off + dst_i] = -funcc[off + src_i];
+                }
+            }
+        }
+    } else {
+        for (int ii = 0; ii < ord; ++ii) {
+            const size_t dst_i = (size_t)(shift - ii);
+            const size_t src_i = (size_t)(shift + ii + 1);
+            for (int k0 = 0; k0 < extc3; ++k0) {
+                const size_t kbase = interior_k + (size_t)k0 * splane + interior_j;
+                for (int j0 = 0; j0 < extc2; ++j0) {
+                    const size_t off = kbase + (size_t)j0 * snx;
+                    funcc[off + dst_i] = funcc[off + src_i] * s1;
+                }
+            }
+        }
+    }
+
+    /* 3) funcc(:,-j,1:extc3) = funcc(:,j+1,1:extc3)*SoA(2) */
+    const double s2 = SoA[1];
+    if (s2 == 1.0) {
+        for (int jj = 0; jj < ord; ++jj) {
+            const size_t dst_j = (size_t)(shift - jj) * snx;
+            const size_t src_j = (size_t)(shift + jj + 1) * snx;
+            for (int k0 = 0; k0 < extc3; ++k0) {
+                const size_t kbase = interior_k + (size_t)k0 * splane;
+                double *dst = funcc + kbase + dst_j;
+                const double *src = funcc + kbase + src_j;
+                for (int i = 0; i < nx; ++i) dst[i] = src[i];
+            }
+        }
+    } else if (s2 == -1.0) {
+        for (int jj = 0; jj < ord; ++jj) {
+            const size_t dst_j = (size_t)(shift - jj) * snx;
+            const size_t src_j = (size_t)(shift + jj + 1) * snx;
+            for (int k0 = 0; k0 < extc3; ++k0) {
+                const size_t kbase = interior_k + (size_t)k0 * splane;
+                double *dst = funcc + kbase + dst_j;
+                const double *src = funcc + kbase + src_j;
+                for (int i = 0; i < nx; ++i) dst[i] = -src[i];
+            }
+        }
+    } else {
+        for (int jj = 0; jj < ord; ++jj) {
+            const size_t dst_j = (size_t)(shift - jj) * snx;
+            const size_t src_j = (size_t)(shift + jj + 1) * snx;
+            for (int k0 = 0; k0 < extc3; ++k0) {
+                const size_t kbase = interior_k + (size_t)k0 * splane;
+                double *dst = funcc + kbase + dst_j;
+                const double *src = funcc + kbase + src_j;
+                for (int i = 0; i < nx; ++i) dst[i] = src[i] * s2;
+            }
+        }
+    }
+
+    /* 4) funcc(:,:,-k) = funcc(:,:,k+1)*SoA(3) */
+    const double s3 = SoA[2];
+    if (s3 == 1.0) {
+        for (int kk = 0; kk < ord; ++kk) {
+            const size_t dst_k = (size_t)(shift - kk) * splane;
+            const size_t src_k = (size_t)(shift + kk + 1) * splane;
+            double *dst = funcc + dst_k;
+            const double *src = funcc + src_k;
+            for (size_t p = 0; p < splane; ++p) dst[p] = src[p];
+        }
+    } else if (s3 == -1.0) {
+        for (int kk = 0; kk < ord; ++kk) {
+            const size_t dst_k = (size_t)(shift - kk) * splane;
+            const size_t src_k = (size_t)(shift + kk + 1) * splane;
+            double *dst = funcc + dst_k;
+            const double *src = funcc + src_k;
+            for (size_t p = 0; p < splane; ++p) dst[p] = -src[p];
+        }
+    } else {
+        for (int kk = 0; kk < ord; ++kk) {
+            const size_t dst_k = (size_t)(shift - kk) * splane;
+            const size_t src_k = (size_t)(shift + kk + 1) * splane;
+            double *dst = funcc + dst_k;
+            const double *src = funcc + src_k;
+            for (size_t p = 0; p < splane; ++p) dst[p] = src[p] * s3;
+        }
+    }
+}
+
 static inline void symmetry_bd(int ord,
                  const int extc[3],
                  const double *func,
                  double *funcc,
                  const double SoA[3])
 {
-    const int extc1 = extc[0], extc2 = extc[1], extc3 = extc[2];
+    if (ord <= 0) return;
 
-    // 1) funcc(1:extc1,1:extc2,1:extc3) = func
-    // Fortran 的 (iF=1..extc1) 对应 C 的 func(i0=0..extc1-1)
-    for (int k0 = 0; k0 < extc3; ++k0) {
-        for (int j0 = 0; j0 < extc2; ++j0) {
-            for (int i0 = 0; i0 < extc1; ++i0) {
-                const int iF = i0 + 1, jF = j0 + 1, kF = k0 + 1;
-                funcc[idx_funcc_F(iF, jF, kF, ord, extc)] = func[idx_func0(i0, j0, k0, extc)];
-            }
-        }
+    /* Fast paths used by current C kernels: ord=2 (derivs), ord=3 (lopsided/KO). */
+    if (ord == 2) {
+        symmetry_bd_impl(2, 1, extc, func, funcc, SoA);
+        return;
+    }
+    if (ord == 3) {
+        symmetry_bd_impl(3, 2, extc, func, funcc, SoA);
+        return;
     }
 
-    // 2) do i=0..ord-1: funcc(-i, 1:extc2, 1:extc3) = funcc(i+1, ...)*SoA(1)
-    for (int ii = 0; ii <= ord - 1; ++ii) {
-        const int iF_dst = -ii;       // 0, -1, -2, ...
-        const int iF_src = ii + 1;    // 1, 2, 3, ...
-        for (int kF = 1; kF <= extc3; ++kF) {
-            for (int jF = 1; jF <= extc2; ++jF) {
-                funcc[idx_funcc_F(iF_dst, jF, kF, ord, extc)] =
-                    funcc[idx_funcc_F(iF_src, jF, kF, ord, extc)] * SoA[0];
-            }
-        }
-    }
-
-    // 3) do i=0..ord-1: funcc(:,-i, 1:extc3) = funcc(:, i+1, 1:extc3)*SoA(2)
-    // 注意 Fortran 这里的 ":" 表示 iF 从 (-ord+1..extc1) 全覆盖
-    for (int jj = 0; jj <= ord - 1; ++jj) {
-        const int jF_dst = -jj;
-        const int jF_src = jj + 1;
-        for (int kF = 1; kF <= extc3; ++kF) {
-            for (int iF = -ord + 1; iF <= extc1; ++iF) {
-                funcc[idx_funcc_F(iF, jF_dst, kF, ord, extc)] =
-                    funcc[idx_funcc_F(iF, jF_src, kF, ord, extc)] * SoA[1];
-            }
-        }
-    }
-
-    // 4) do i=0..ord-1: funcc(:,:,-i) = funcc(:,:, i+1)*SoA(3)
-    for (int kk = 0; kk <= ord - 1; ++kk) {
-        const int kF_dst = -kk;
-        const int kF_src = kk + 1;
-        for (int jF = -ord + 1; jF <= extc2; ++jF) {
-            for (int iF = -ord + 1; iF <= extc1; ++iF) {
-                funcc[idx_funcc_F(iF, jF, kF_dst, ord, extc)] =
-                    funcc[idx_funcc_F(iF, jF, kF_src, ord, extc)] * SoA[2];
-            }
-        }
-    }
+    symmetry_bd_impl(ord, ord - 1, extc, func, funcc, SoA);
 }
 #endif

makefile_and_run.py

@@ -43,7 +43,8 @@ def get_last_n_cores_per_socket(n=32):
     cpu_str = ",".join(segments)
     total = len(segments) * n
     print(f" CPU binding: taskset -c {cpu_str}  ({total} cores, last {n} per socket)")
-    return f"taskset -c {cpu_str}"
+    #return f"taskset -c {cpu_str}"
+    return f""
 
 
 ## CPU core binding: dynamically select the last 32 cores of each socket (64 cores total)