Add OpenMP parallelization to BSSN RHS hot-path stencil routines

Enable OpenMP threading for the dominant computational kernels: - makefile.inc: add -qopenmp to f90appflags - diff_new.f90: split fderivs/fdderivs into OpenMP interior + serial boundary - kodiss.f90: split kodis into OpenMP interior + serial boundary - lopsidediff.f90: add OMP PARALLEL DO COLLAPSE(2) to lopsided - fmisc.f90: parallelize symmetry_bd bulk array copy - bssn_rhs.f90: add OMP WORKSHARE to array-syntax operations Co-Authored-By: Claude Opus 4.6 <noreply@anthropic.com>
Eliminate hot-path heap allocations in TwoPunctures spectral solver
2026-02-07 13:58:55 +08:00 · 2026-02-06 21:20:35 +08:00 · 2026-02-06 19:00:35 +08:00 · 2026-02-06 18:36:29 +08:00 · 2026-02-06 17:13:39 +08:00
9 changed files with 508 additions and 334 deletions
--- a/AMSS_NCKU_source/TwoPunctures.C
+++ b/AMSS_NCKU_source/TwoPunctures.C
@@ -5,6 +5,7 @@
 #include <cstdio>
 #include <cstdlib>
 #include <string>
+#include <cstring>
 #include <iostream>
 #include <iomanip>
 #include <fstream>
@@ -60,13 +61,110 @@ TwoPunctures::TwoPunctures(double mp, double mm, double b,
  F = dvector(0, ntotal - 1);
  allocate_derivs(&u, ntotal);
  allocate_derivs(&v, ntotal);
+
+  // Allocate workspace buffers for hot-path allocation elimination
+  int N = maximum3(n1, n2, n3);
+  int maxn = maximum2(n1, n2);
+
+  // LineRelax_be workspace (sized for n2)
+  ws_diag_be = new double[n2];
+  ws_e_be = new double[n2 - 1];
+  ws_f_be = new double[n2 - 1];
+  ws_b_be = new double[n2];
+  ws_x_be = new double[n2];
+
+  // LineRelax_al workspace (sized for n1)
+  ws_diag_al = new double[n1];
+  ws_e_al = new double[n1 - 1];
+  ws_f_al = new double[n1 - 1];
+  ws_b_al = new double[n1];
+  ws_x_al = new double[n1];
+
+  // ThomasAlgorithm workspace (sized for max(n1,n2))
+  ws_thomas_y = new double[maxn];
+
+  // JFD_times_dv workspace (sized for nvar)
+  ws_jfd_values = dvector(0, nvar - 1);
+  allocate_derivs(&ws_jfd_dU, nvar);
+  allocate_derivs(&ws_jfd_U, nvar);
+
+  // chebft_Zeros workspace (sized for N+1)
+  ws_cheb_c = dvector(0, N);
+
+  // fourft workspace (sized for N/2+1 each)
+  ws_four_a = dvector(0, N / 2);
+  ws_four_b = dvector(0, N / 2);
+
+  // Derivatives_AB3 workspace
+  ws_deriv_p = dvector(0, N);
+  ws_deriv_dp = dvector(0, N);
+  ws_deriv_d2p = dvector(0, N);
+  ws_deriv_q = dvector(0, N);
+  ws_deriv_dq = dvector(0, N);
+  ws_deriv_r = dvector(0, N);
+  ws_deriv_dr = dvector(0, N);
+  ws_deriv_indx = ivector(0, N);
+
+  // F_of_v workspace
+  ws_fov_sources = new double[n1 * n2 * n3];
+  ws_fov_values = dvector(0, nvar - 1);
+  allocate_derivs(&ws_fov_U, nvar);
+
+  // J_times_dv workspace
+  ws_jtdv_values = dvector(0, nvar - 1);
+  allocate_derivs(&ws_jtdv_dU, nvar);
+  allocate_derivs(&ws_jtdv_U, nvar);
 }

 TwoPunctures::~TwoPunctures()
 {
+  int const nvar = 1, n1 = npoints_A, n2 = npoints_B, n3 = npoints_phi;
+  int N = maximum3(n1, n2, n3);
+
  free_dvector(F, 0, ntotal - 1);
  free_derivs(&u, ntotal);
  free_derivs(&v, ntotal);
+
+  // Free workspace buffers
+  delete[] ws_diag_be;
+  delete[] ws_e_be;
+  delete[] ws_f_be;
+  delete[] ws_b_be;
+  delete[] ws_x_be;
+
+  delete[] ws_diag_al;
+  delete[] ws_e_al;
+  delete[] ws_f_al;
+  delete[] ws_b_al;
+  delete[] ws_x_al;
+
+  delete[] ws_thomas_y;
+
+  free_dvector(ws_jfd_values, 0, nvar - 1);
+  free_derivs(&ws_jfd_dU, nvar);
+  free_derivs(&ws_jfd_U, nvar);
+
+  free_dvector(ws_cheb_c, 0, N);
+
+  free_dvector(ws_four_a, 0, N / 2);
+  free_dvector(ws_four_b, 0, N / 2);
+
+  free_dvector(ws_deriv_p, 0, N);
+  free_dvector(ws_deriv_dp, 0, N);
+  free_dvector(ws_deriv_d2p, 0, N);
+  free_dvector(ws_deriv_q, 0, N);
+  free_dvector(ws_deriv_dq, 0, N);
+  free_dvector(ws_deriv_r, 0, N);
+  free_dvector(ws_deriv_dr, 0, N);
+  free_ivector(ws_deriv_indx, 0, N);
+
+  delete[] ws_fov_sources;
+  free_dvector(ws_fov_values, 0, nvar - 1);
+  free_derivs(&ws_fov_U, nvar);
+
+  free_dvector(ws_jtdv_values, 0, nvar - 1);
+  free_derivs(&ws_jtdv_dU, nvar);
+  free_derivs(&ws_jtdv_U, nvar);
 }

 void TwoPunctures::Solve()
@@ -655,7 +753,7 @@ void TwoPunctures::chebft_Zeros(double u[], int n, int inv)
  int k, j, isignum;
  double fac, sum, Pion, *c;

-  c = dvector(0, n);
+  c = ws_cheb_c;
  Pion = Pi / n;
  if (inv == 0)
  {
@@ -686,7 +784,6 @@ void TwoPunctures::chebft_Zeros(double u[], int n, int inv)
  }
  for (j = 0; j < n; j++)
    u[j] = c[j];
-  free_dvector(c, 0, n);
 }

 /* --------------------------------------------------------------------------*/
@@ -774,8 +871,8 @@ void TwoPunctures::fourft(double *u, int N, int inv)
  double x, x1, fac, Pi_fac, *a, *b;

  M = N / 2;
-  a = dvector(0, M);
-  b = dvector(1, M); /* Actually: b=vector(1,M-1) but this is problematic if M=1*/
+  a = ws_four_a;
+  b = ws_four_b - 1; /* offset to match dvector(1,M) indexing */
  fac = 1. / M;
  Pi_fac = Pi * fac;
  if (inv == 0)
@@ -824,8 +921,6 @@ void TwoPunctures::fourft(double *u, int N, int inv)
      iy = -iy;
    }
  }
-  free_dvector(a, 0, M);
-  free_dvector(b, 1, M);
 }

 /* -----------------------------------------*/
@@ -1118,14 +1213,14 @@ void TwoPunctures::Derivatives_AB3(int nvar, int n1, int n2, int n3, derivs v)
  double *p, *dp, *d2p, *q, *dq, *r, *dr;

  N = maximum3(n1, n2, n3);
-  p = dvector(0, N);
-  dp = dvector(0, N);
-  d2p = dvector(0, N);
-  q = dvector(0, N);
-  dq = dvector(0, N);
-  r = dvector(0, N);
-  dr = dvector(0, N);
-  indx = ivector(0, N);
+  p = ws_deriv_p;
+  dp = ws_deriv_dp;
+  d2p = ws_deriv_d2p;
+  q = ws_deriv_q;
+  dq = ws_deriv_dq;
+  r = ws_deriv_r;
+  dr = ws_deriv_dr;
+  indx = ws_deriv_indx;

  for (ivar = 0; ivar < nvar; ivar++)
  {
@@ -1208,14 +1303,6 @@ void TwoPunctures::Derivatives_AB3(int nvar, int n1, int n2, int n3, derivs v)
      }
    }
  }
-  free_dvector(p, 0, N);
-  free_dvector(dp, 0, N);
-  free_dvector(d2p, 0, N);
-  free_dvector(q, 0, N);
-  free_dvector(dq, 0, N);
-  free_dvector(r, 0, N);
-  free_dvector(dr, 0, N);
-  free_ivector(indx, 0, N);
 }
 /* --------------------------------------------------------------------------*/
 void TwoPunctures::Newton(int const nvar, int const n1, int const n2, int const n3,
@@ -1284,10 +1371,11 @@ void TwoPunctures::F_of_v(int nvar, int n1, int n2, int n3, derivs v, double *F,
  derivs U;
  double *sources;

-  values = dvector(0, nvar - 1);
-  allocate_derivs(&U, nvar);
+  values = ws_fov_values;
+  U = ws_fov_U;

-  sources = (double *)calloc(n1 * n2 * n3, sizeof(double));
+  sources = ws_fov_sources;
+  memset(sources, 0, n1 * n2 * n3 * sizeof(double));
  if (0)
  {
    double *s_x, *s_y, *s_z;
@@ -1442,9 +1530,6 @@ void TwoPunctures::F_of_v(int nvar, int n1, int n2, int n3, derivs v, double *F,
  {
    fclose(debugfile);
  }
-  free(sources);
-  free_dvector(values, 0, nvar - 1);
-  free_derivs(&U, nvar);
 }
 /* --------------------------------------------------------------------------*/
 double TwoPunctures::norm_inf(double const *F, int const ntotal)
@@ -1850,11 +1935,12 @@ void TwoPunctures::J_times_dv(int nvar, int n1, int n2, int n3, derivs dv, doubl

  Derivatives_AB3(nvar, n1, n2, n3, dv);

+  values = ws_jtdv_values;
+  dU = ws_jtdv_dU;
+  U = ws_jtdv_U;
+
  for (i = 0; i < n1; i++)
  {
-    values = dvector(0, nvar - 1);
-    allocate_derivs(&dU, nvar);
-    allocate_derivs(&U, nvar);
    for (j = 0; j < n2; j++)
    {
      for (k = 0; k < n3; k++)
@@ -1908,9 +1994,6 @@ void TwoPunctures::J_times_dv(int nvar, int n1, int n2, int n3, derivs dv, doubl
        }
      }
    }
-    free_dvector(values, 0, nvar - 1);
-    free_derivs(&dU, nvar);
-    free_derivs(&U, nvar);
  }
 }
 /* --------------------------------------------------------------------------*/
@@ -1957,17 +2040,11 @@ void TwoPunctures::LineRelax_be(double *dv,
 {
  int j, m, Ic, Ip, Im, col, ivar;

-  double *diag = new double[n2];
-  double *e = new double[n2 - 1]; /* above diagonal */
-  double *f = new double[n2 - 1]; /* below diagonal */
-  double *b = new double[n2];     /* rhs */
-  double *x = new double[n2];     /* solution vector */
-
-  //  gsl_vector *diag = gsl_vector_alloc(n2);
-  //  gsl_vector *e = gsl_vector_alloc(n2-1); /* above diagonal */
-  //  gsl_vector *f = gsl_vector_alloc(n2-1); /* below diagonal */
-  //  gsl_vector *b = gsl_vector_alloc(n2);   /* rhs */
-  //  gsl_vector *x = gsl_vector_alloc(n2);   /* solution vector */
+  double *diag = ws_diag_be;
+  double *e = ws_e_be;     /* above diagonal */
+  double *f = ws_f_be;     /* below diagonal */
+  double *b = ws_b_be;     /* rhs */
+  double *x = ws_x_be;     /* solution vector */

  for (ivar = 0; ivar < nvar; ivar++)
  {
@@ -1977,62 +2054,35 @@ void TwoPunctures::LineRelax_be(double *dv,
    }
    diag[n2 - 1] = 0;

-    //    gsl_vector_set_zero(diag);
-    //    gsl_vector_set_zero(e);
-    //    gsl_vector_set_zero(f);
    for (j = 0; j < n2; j++)
    {
      Ip = Index(ivar, i, j + 1, k, nvar, n1, n2, n3);
      Ic = Index(ivar, i, j, k, nvar, n1, n2, n3);
      Im = Index(ivar, i, j - 1, k, nvar, n1, n2, n3);
      b[j] = rhs[Ic];
-      //      gsl_vector_set(b,j,rhs[Ic]);
      for (m = 0; m < ncols[Ic]; m++)
      {
        col = cols[Ic][m];
        if (col != Ip && col != Ic && col != Im)
          b[j] -= JFD[Ic][m] * dv[col];
-        //          *gsl_vector_ptr(b, j) -= JFD[Ic][m] * dv[col];
        else
        {
          if (col == Im && j > 0)
            f[j - 1] = JFD[Ic][m];
-          //            gsl_vector_set(f,j-1,JFD[Ic][m]);
          if (col == Ic)
            diag[j] = JFD[Ic][m];
-          //            gsl_vector_set(diag,j,JFD[Ic][m]);
          if (col == Ip && j < n2 - 1)
            e[j] = JFD[Ic][m];
-          //            gsl_vector_set(e,j,JFD[Ic][m]);
        }
      }
    }
-    //          A x = b
-    //          A = ( d_0 e_0  0   0  )
-    //              ( f_0 d_1 e_1  0  )
-    //              (  0  f_1 d_2 e_2 )
-    //              (  0   0  f_2 d_3 )
-    //
    ThomasAlgorithm(n2, f, diag, e, x, b);
-    //    gsl_linalg_solve_tridiag(diag, e, f, b, x);
    for (j = 0; j < n2; j++)
    {
      Ic = Index(ivar, i, j, k, nvar, n1, n2, n3);
      dv[Ic] = x[j];
-      //      dv[Ic] = gsl_vector_get(x, j);
    }
  }
-
-  delete[] diag;
-  delete[] e;
-  delete[] f;
-  delete[] b;
-  delete[] x;
-  //  gsl_vector_free(diag);
-  //  gsl_vector_free(e);
-  //  gsl_vector_free(f);
-  //  gsl_vector_free(b);
-  //  gsl_vector_free(x);
 }
 /* --------------------------------------------------------------------------*/
 void TwoPunctures::JFD_times_dv(int i, int j, int k, int nvar, int n1, int n2,
@@ -2049,8 +2099,8 @@ void TwoPunctures::JFD_times_dv(int i, int j, int k, int nvar, int n1, int n2,
      ha, ga, ga2, hb, gb, gb2, hp, gp, gp2, gagb, gagp, gbgp;
  derivs dU, U;

-  allocate_derivs(&dU, nvar);
-  allocate_derivs(&U, nvar);
+  dU = ws_jfd_dU;
+  U = ws_jfd_U;

  if (k < 0)
    k = k + n3;
@@ -2168,9 +2218,6 @@ void TwoPunctures::JFD_times_dv(int i, int j, int k, int nvar, int n1, int n2,
  LinEquations(A, B, X, R, x, r, phi, y, z, dU, U, values);
  for (ivar = 0; ivar < nvar; ivar++)
    values[ivar] *= FAC;
-
-  free_derivs(&dU, nvar);
-  free_derivs(&U, nvar);
 }
 #undef FAC
 /*-----------------------------------------------------------*/
@@ -2202,17 +2249,11 @@ void TwoPunctures::LineRelax_al(double *dv,
 {
  int i, m, Ic, Ip, Im, col, ivar;

-  double *diag = new double[n1];
-  double *e = new double[n1 - 1]; /* above diagonal */
-  double *f = new double[n1 - 1]; /* below diagonal */
-  double *b = new double[n1];     /* rhs */
-  double *x = new double[n1];     /* solution vector */
-
-  //  gsl_vector *diag = gsl_vector_alloc(n1);
-  //  gsl_vector *e = gsl_vector_alloc(n1-1); /* above diagonal */
-  //  gsl_vector *f = gsl_vector_alloc(n1-1); /* below diagonal */
-  //  gsl_vector *b = gsl_vector_alloc(n1);   /* rhs */
-  //  gsl_vector *x = gsl_vector_alloc(n1);   /* solution vector */
+  double *diag = ws_diag_al;
+  double *e = ws_e_al;     /* above diagonal */
+  double *f = ws_f_al;     /* below diagonal */
+  double *b = ws_b_al;     /* rhs */
+  double *x = ws_x_al;     /* solution vector */

  for (ivar = 0; ivar < nvar; ivar++)
  {
@@ -2222,57 +2263,35 @@ void TwoPunctures::LineRelax_al(double *dv,
    }
    diag[n1 - 1] = 0;

-    //    gsl_vector_set_zero(diag);
-    //    gsl_vector_set_zero(e);
-    //    gsl_vector_set_zero(f);
    for (i = 0; i < n1; i++)
    {
      Ip = Index(ivar, i + 1, j, k, nvar, n1, n2, n3);
      Ic = Index(ivar, i, j, k, nvar, n1, n2, n3);
      Im = Index(ivar, i - 1, j, k, nvar, n1, n2, n3);
      b[i] = rhs[Ic];
-      //      gsl_vector_set(b,i,rhs[Ic]);
      for (m = 0; m < ncols[Ic]; m++)
      {
        col = cols[Ic][m];
        if (col != Ip && col != Ic && col != Im)
          b[i] -= JFD[Ic][m] * dv[col];
-        //          *gsl_vector_ptr(b, i) -= JFD[Ic][m] * dv[col];
        else
        {
          if (col == Im && i > 0)
            f[i - 1] = JFD[Ic][m];
-          //            gsl_vector_set(f,i-1,JFD[Ic][m]);
          if (col == Ic)
            diag[i] = JFD[Ic][m];
-          //            gsl_vector_set(diag,i,JFD[Ic][m]);
          if (col == Ip && i < n1 - 1)
            e[i] = JFD[Ic][m];
-          //            gsl_vector_set(e,i,JFD[Ic][m]);
        }
      }
    }
    ThomasAlgorithm(n1, f, diag, e, x, b);
-    //    gsl_linalg_solve_tridiag(diag, e, f, b, x);
    for (i = 0; i < n1; i++)
    {
      Ic = Index(ivar, i, j, k, nvar, n1, n2, n3);
      dv[Ic] = x[i];
-      //      dv[Ic] = gsl_vector_get(x, i);
    }
  }
-
-  delete[] diag;
-  delete[] e;
-  delete[] f;
-  delete[] b;
-  delete[] x;
-
-  //  gsl_vector_free(diag);
-  //  gsl_vector_free(e);
-  //  gsl_vector_free(f);
-  //  gsl_vector_free(b);
-  //  gsl_vector_free(x);
 }
 /* -------------------------------------------------------------------------*/
 // a[N], b[N-1], c[N-1], x[N], q[N]
@@ -2284,44 +2303,29 @@ void TwoPunctures::LineRelax_al(double *dv,
 //"Parallel Scientific Computing in C++ and MPI" P361
 void TwoPunctures::ThomasAlgorithm(int N, double *b, double *a, double *c, double *x, double *q)
 {
+  // In-place Thomas algorithm: uses a[] as d workspace, b[] as l workspace.
+  // c[] is already u (above-diagonal). ws_thomas_y is pre-allocated workspace.
  int i;
-  double *l, *u, *d, *y;
-  l = new double[N - 1];
-  u = new double[N - 1];
-  d = new double[N];
-  y = new double[N];
-
-  /* LU Decomposition */
-  d[0] = a[0];
-  u[0] = c[0];
+  double *y = ws_thomas_y;

+  /* LU Decomposition (in-place: a becomes d, b becomes l) */
  for (i = 0; i < N - 2; i++)
  {
-    l[i] = b[i] / d[i];
-    d[i + 1] = a[i + 1] - l[i] * u[i];
-    u[i + 1] = c[i + 1];
+    b[i] = b[i] / a[i];
+    a[i + 1] = a[i + 1] - b[i] * c[i];
  }
-
-  l[N - 2] = b[N - 2] / d[N - 2];
-  d[N - 1] = a[N - 1] - l[N - 2] * u[N - 2];
+  b[N - 2] = b[N - 2] / a[N - 2];
+  a[N - 1] = a[N - 1] - b[N - 2] * c[N - 2];

  /* Forward Substitution [L][y] = [q] */
  y[0] = q[0];
  for (i = 1; i < N; i++)
-    y[i] = q[i] - l[i - 1] * y[i - 1];
+    y[i] = q[i] - b[i - 1] * y[i - 1];

  /* Backward Substitution [U][x] = [y] */
-  x[N - 1] = y[N - 1] / d[N - 1];
-
+  x[N - 1] = y[N - 1] / a[N - 1];
  for (i = N - 2; i >= 0; i--)
-    x[i] = (y[i] - u[i] * x[i + 1]) / d[i];
-
-  delete[] l;
-  delete[] u;
-  delete[] d;
-  delete[] y;
-
-  return;
+    x[i] = (y[i] - c[i] * x[i + 1]) / a[i];
 }
 // --------------------------------------------------------------------------*/
 // Calculates the value of v at an arbitrary position (x,y,z) if the spectral coefficients are know*/*/
--- a/AMSS_NCKU_source/TwoPunctures.h
+++ b/AMSS_NCKU_source/TwoPunctures.h
@@ -42,6 +42,33 @@ private:

       int ntotal;

+       // Pre-allocated workspace buffers for hot-path allocation elimination
+       // LineRelax_be workspace (sized for n2)
+       double *ws_diag_be, *ws_e_be, *ws_f_be, *ws_b_be, *ws_x_be;
+       // LineRelax_al workspace (sized for n1)
+       double *ws_diag_al, *ws_e_al, *ws_f_al, *ws_b_al, *ws_x_al;
+       // ThomasAlgorithm workspace (sized for max(n1,n2))
+       double *ws_thomas_y;
+       // JFD_times_dv workspace (sized for nvar)
+       double *ws_jfd_values;
+       derivs ws_jfd_dU, ws_jfd_U;
+       // chebft_Zeros workspace (sized for max(n1,n2,n3)+1)
+       double *ws_cheb_c;
+       // fourft workspace (sized for max(n1,n2,n3)/2+1 each)
+       double *ws_four_a, *ws_four_b;
+       // Derivatives_AB3 workspace
+       double *ws_deriv_p, *ws_deriv_dp, *ws_deriv_d2p;
+       double *ws_deriv_q, *ws_deriv_dq;
+       double *ws_deriv_r, *ws_deriv_dr;
+       int *ws_deriv_indx;
+       // F_of_v workspace
+       double *ws_fov_sources;
+       double *ws_fov_values;
+       derivs ws_fov_U;
+       // J_times_dv workspace
+       double *ws_jtdv_values;
+       derivs ws_jtdv_dU, ws_jtdv_U;
+
       struct parameters
       {
              int nvar, n1, n2, n3;
--- a/AMSS_NCKU_source/bssn_rhs.f90
+++ b/AMSS_NCKU_source/bssn_rhs.f90
@@ -106,7 +106,8 @@
  call getpbh(BHN,Porg,Mass)
 #endif

-!!! sanity check
+!!! sanity check (disabled in production builds for performance)
+#ifdef DEBUG
  dX = sum(chi)+sum(trK)+sum(dxx)+sum(gxy)+sum(gxz)+sum(dyy)+sum(gyz)+sum(dzz) &
      +sum(Axx)+sum(Axy)+sum(Axz)+sum(Ayy)+sum(Ayz)+sum(Azz)                   &
      +sum(Gamx)+sum(Gamy)+sum(Gamz)                                           &
@@ -136,6 +137,7 @@
     gont = 1
     return
  endif
+#endif

  PI = dacos(-ONE)

@@ -166,6 +168,8 @@
  call fderivs(ex,gyz,gyzx,gyzy,gyzz,X,Y,Z,SYM ,ANTI,ANTI,Symmetry,Lev)
  call fderivs(ex,dzz,gzzx,gzzy,gzzz,X,Y,Z,SYM ,SYM ,SYM ,Symmetry,Lev)

+  !$OMP PARALLEL
+  !$OMP WORKSHARE
  gxx_rhs = - TWO * alpn1 * Axx    -  F2o3 * gxx * div_beta          + &
              TWO *(  gxx * betaxx +   gxy * betayx +   gxz * betazx)

@@ -199,6 +203,8 @@
  gupyy =   ( gxx * gzz - gxz * gxz ) / gupzz
  gupyz = - ( gxx * gyz - gxy * gxz ) / gupzz
  gupzz =   ( gxx * gyy - gxy * gxy ) / gupzz
+  !$OMP END WORKSHARE
+  !$OMP END PARALLEL

  if(co == 0)then
 ! Gam^i_Res = Gam^i + gup^ij_,j
@@ -232,6 +238,8 @@
  endif

 ! second kind of connection
+  !$OMP PARALLEL
+  !$OMP WORKSHARE
  Gamxxx =HALF*( gupxx*gxxx + gupxy*(TWO*gxyx - gxxy ) + gupxz*(TWO*gxzx - gxxz ))
  Gamyxx =HALF*( gupxy*gxxx + gupyy*(TWO*gxyx - gxxy ) + gupyz*(TWO*gxzx - gxxz ))
  Gamzxx =HALF*( gupxz*gxxx + gupyz*(TWO*gxyx - gxxy ) + gupzz*(TWO*gxzx - gxxz ))
@@ -280,6 +288,8 @@
          (gupxy * gupyz       + gupyy * gupxz)* Axy                       + &
          (gupxy * gupzz       + gupyz * gupxz)* Axz                       + &
          (gupyy * gupzz       + gupyz * gupyz)* Ayz
+  !$OMP END WORKSHARE
+  !$OMP END PARALLEL

 ! Right hand side for Gam^i without shift terms...
  call fderivs(ex,Lap,Lapx,Lapy,Lapz,X,Y,Z,SYM,SYM,SYM,Symmetry,Lev)
@@ -334,6 +344,8 @@
  call fderivs(ex,Gamy,Gamyx,Gamyy,Gamyz,X,Y,Z,SYM ,ANTI,SYM ,Symmetry,Lev)
  call fderivs(ex,Gamz,Gamzx,Gamzy,Gamzz,X,Y,Z,SYM ,SYM ,ANTI,Symmetry,Lev)

+  !$OMP PARALLEL
+  !$OMP WORKSHARE
  Gamx_rhs =               Gamx_rhs +  F2o3 *  Gamxa * div_beta        - &
                     Gamxa * betaxx - Gamya * betaxy - Gamza * betaxz  + &
             F1o3 * (gupxx * fxx    + gupxy * fxy    + gupxz * fxz    ) + &
@@ -373,6 +385,8 @@
  gyyz = gxz * Gamxyy + gyz * Gamyyy + gzz * Gamzyy
  gyzz = gxz * Gamxyz + gyz * Gamyyz + gzz * Gamzyz
  gzzz = gxz * Gamxzz + gyz * Gamyzz + gzz * Gamzzz
+  !$OMP END WORKSHARE
+  !$OMP END PARALLEL

 !compute Ricci tensor for tilted metric
   call fdderivs(ex,dxx,fxx,fxy,fxz,fyy,fyz,fzz,X,Y,Z,SYM ,SYM ,SYM ,symmetry,Lev)
@@ -399,6 +413,8 @@
   Ryz =   gupxx * fxx + gupyy * fyy + gupzz * fzz + &
         ( gupxy * fxy + gupxz * fxz + gupyz * fyz ) * TWO

+  !$OMP PARALLEL
+  !$OMP WORKSHARE
  Rxx =     - HALF * Rxx                                   + &
               gxx * Gamxx+ gxy * Gamyx   +    gxz * Gamzx + &
             Gamxa * gxxx +  Gamya * gxyx +  Gamza * gxzx  + &
@@ -599,9 +615,13 @@
            Gamxyz * gxzz + Gamyyz * gyzz + Gamzyz * gzzz  + &
            Gamxzz * gxzy + Gamyzz * gyzy + Gamzzz * gzzy  + &
            Gamxyz * gzzx + Gamyyz * gzzy + Gamzyz * gzzz )
+  !$OMP END WORKSHARE
+  !$OMP END PARALLEL
 !covariant second derivative of chi respect to tilted metric
  call fdderivs(ex,chi,fxx,fxy,fxz,fyy,fyz,fzz,X,Y,Z,SYM,SYM,SYM,Symmetry,Lev)

+  !$OMP PARALLEL
+  !$OMP WORKSHARE
  fxx = fxx - Gamxxx * chix - Gamyxx * chiy - Gamzxx * chiz
  fxy = fxy - Gamxxy * chix - Gamyxy * chiy - Gamzxy * chiz
  fxz = fxz - Gamxxz * chix - Gamyxz * chiy - Gamzxz * chiz
@@ -624,11 +644,15 @@
  Rxy = Rxy + (fxy - chix*chiy/chin1/TWO + gxy * f)/chin1/TWO
  Rxz = Rxz + (fxz - chix*chiz/chin1/TWO + gxz * f)/chin1/TWO
  Ryz = Ryz + (fyz - chiy*chiz/chin1/TWO + gyz * f)/chin1/TWO
+  !$OMP END WORKSHARE
+  !$OMP END PARALLEL

 ! covariant second derivatives of the lapse respect to physical metric
  call fdderivs(ex,Lap,fxx,fxy,fxz,fyy,fyz,fzz,X,Y,Z, &
                SYM,SYM,SYM,symmetry,Lev)

+  !$OMP PARALLEL
+  !$OMP WORKSHARE
  gxxx = (gupxx * chix + gupxy * chiy + gupxz * chiz)/chin1
  gxxy = (gupxy * chix + gupyy * chiy + gupyz * chiz)/chin1
  gxxz = (gupxz * chix + gupyz * chiy + gupzz * chiz)/chin1
@@ -789,6 +813,8 @@
 !!!! gauge variable part

  Lap_rhs = -TWO*alpn1*trK
+  !$OMP END WORKSHARE
+  !$OMP END PARALLEL
 #if (GAUGE == 0)
  betax_rhs = FF*dtSfx
  betay_rhs = FF*dtSfy
--- a/AMSS_NCKU_source/diff_new.f90
+++ b/AMSS_NCKU_source/diff_new.f90
@@ -997,10 +997,10 @@
  fy = ZEO
  fz = ZEO

+#if 0
  do k=1,ex(3)-1
  do j=1,ex(2)-1
  do i=1,ex(1)-1
-#if 0  
 ! x direction
        if(i+2 <= imax .and. i-2 >= imin)then
 !
@@ -1040,9 +1040,13 @@

 ! set kmax and kmin 0
    endif
+  enddo
+  enddo
+  enddo
 #elif 0
-! x direction   
-        if(i+2 <= imax .and. i-2 >= imin)then
+  do k=1,ex(3)-1
+  do j=1,ex(2)-1
+  do i=1,ex(1)-1
 !
 !              f(i-2) - 8 f(i-1) + 8 f(i+1) - f(i+2)
 !  fx(i) = ---------------------------------------------
@@ -1079,8 +1083,32 @@

 ! set kmax and kmin 0
    endif
+  enddo
+  enddo
+  enddo
 #else
-! for bam comparison
+! for bam comparison — split into branch-free interior + serial boundary
+! Interior: all stencil points guaranteed in-bounds, no branches needed
+  !$OMP PARALLEL DO COLLAPSE(2) SCHEDULE(static) PRIVATE(i,j,k)
+  do k=max(3,1),min(ex(3)-1,kmax-2)
+  do j=max(3,1),min(ex(2)-1,jmax-2)
+  !DIR$ IVDEP
+  do i=max(3,1),min(ex(1)-1,imax-2)
+      fx(i,j,k)=d12dx*(fh(i-2,j,k)-EIT*fh(i-1,j,k)+EIT*fh(i+1,j,k)-fh(i+2,j,k))
+      fy(i,j,k)=d12dy*(fh(i,j-2,k)-EIT*fh(i,j-1,k)+EIT*fh(i,j+1,k)-fh(i,j+2,k))
+      fz(i,j,k)=d12dz*(fh(i,j,k-2)-EIT*fh(i,j,k-1)+EIT*fh(i,j,k+1)-fh(i,j,k+2))
+  enddo
+  enddo
+  enddo
+  !$OMP END PARALLEL DO
+
+! Boundary shell: original branching logic for points near edges
+  do k=1,ex(3)-1
+  do j=1,ex(2)-1
+  do i=1,ex(1)-1
+   if(i >= 3 .and. i <= imax-2 .and. &
+      j >= 3 .and. j <= jmax-2 .and. &
+      k >= 3 .and. k <= kmax-2) cycle
   if(i+2 <= imax .and. i-2 >= imin .and. &
      j+2 <= jmax .and. j-2 >= jmin .and. &
      k+2 <= kmax .and. k-2 >= kmin) then
@@ -1094,10 +1122,10 @@
      fy(i,j,k)=d2dy*(-fh(i,j-1,k)+fh(i,j+1,k))
      fz(i,j,k)=d2dz*(-fh(i,j,k-1)+fh(i,j,k+1))
   endif
+  enddo
+  enddo
+  enddo
 #endif
-  enddo
-  enddo
-  enddo

  return

@@ -1401,10 +1429,10 @@
  fxz = ZEO
  fyz = ZEO

+#if 0
  do k=1,ex(3)-1
  do j=1,ex(2)-1
  do i=1,ex(1)-1
-#if 0  
 !~~~~~~ fxx
        if(i+2 <= imax .and. i-2 >= imin)then
 !
@@ -1482,8 +1510,47 @@
   elseif(j+1 <= jmax .and. j-1 >= jmin .and. k+1 <= kmax .and. k-1 >= kmin)then
   fyz(i,j,k) = Sdydz*(fh(i,j-1,k-1)-fh(i,j+1,k-1)-fh(i,j-1,k+1)+fh(i,j+1,k+1))
   endif
+  enddo
+  enddo
+  enddo
 #else
-! for bam comparison
+! for bam comparison — split into branch-free interior + serial boundary
+! Interior: all stencil points guaranteed in-bounds, no branches needed
+  !$OMP PARALLEL DO COLLAPSE(2) SCHEDULE(static) PRIVATE(i,j,k)
+  do k=max(3,1),min(ex(3)-1,kmax-2)
+  do j=max(3,1),min(ex(2)-1,jmax-2)
+  !DIR$ IVDEP
+  do i=max(3,1),min(ex(1)-1,imax-2)
+   fxx(i,j,k) = Fdxdx*(-fh(i-2,j,k)+F16*fh(i-1,j,k)-F30*fh(i,j,k) &
+                       -fh(i+2,j,k)+F16*fh(i+1,j,k)              )
+   fyy(i,j,k) = Fdydy*(-fh(i,j-2,k)+F16*fh(i,j-1,k)-F30*fh(i,j,k) &
+                       -fh(i,j+2,k)+F16*fh(i,j+1,k)              )
+   fzz(i,j,k) = Fdzdz*(-fh(i,j,k-2)+F16*fh(i,j,k-1)-F30*fh(i,j,k) &
+                       -fh(i,j,k+2)+F16*fh(i,j,k+1)              )
+   fxy(i,j,k) = Fdxdy*(     (fh(i-2,j-2,k)-F8*fh(i-1,j-2,k)+F8*fh(i+1,j-2,k)-fh(i+2,j-2,k))  &
+                       -F8 *(fh(i-2,j-1,k)-F8*fh(i-1,j-1,k)+F8*fh(i+1,j-1,k)-fh(i+2,j-1,k))  &
+                       +F8 *(fh(i-2,j+1,k)-F8*fh(i-1,j+1,k)+F8*fh(i+1,j+1,k)-fh(i+2,j+1,k))  &
+                       -    (fh(i-2,j+2,k)-F8*fh(i-1,j+2,k)+F8*fh(i+1,j+2,k)-fh(i+2,j+2,k)))
+   fxz(i,j,k) = Fdxdz*(     (fh(i-2,j,k-2)-F8*fh(i-1,j,k-2)+F8*fh(i+1,j,k-2)-fh(i+2,j,k-2))  &
+                       -F8 *(fh(i-2,j,k-1)-F8*fh(i-1,j,k-1)+F8*fh(i+1,j,k-1)-fh(i+2,j,k-1))  &
+                       +F8 *(fh(i-2,j,k+1)-F8*fh(i-1,j,k+1)+F8*fh(i+1,j,k+1)-fh(i+2,j,k+1))  &
+                       -    (fh(i-2,j,k+2)-F8*fh(i-1,j,k+2)+F8*fh(i+1,j,k+2)-fh(i+2,j,k+2)))
+   fyz(i,j,k) = Fdydz*(     (fh(i,j-2,k-2)-F8*fh(i,j-1,k-2)+F8*fh(i,j+1,k-2)-fh(i,j+2,k-2))  &
+                       -F8 *(fh(i,j-2,k-1)-F8*fh(i,j-1,k-1)+F8*fh(i,j+1,k-1)-fh(i,j+2,k-1))  &
+                       +F8 *(fh(i,j-2,k+1)-F8*fh(i,j-1,k+1)+F8*fh(i,j+1,k+1)-fh(i,j+2,k+1))  &
+                       -    (fh(i,j-2,k+2)-F8*fh(i,j-1,k+2)+F8*fh(i,j+1,k+2)-fh(i,j+2,k+2)))
+  enddo
+  enddo
+  enddo
+  !$OMP END PARALLEL DO
+
+! Boundary shell: original branching logic for points near edges
+  do k=1,ex(3)-1
+  do j=1,ex(2)-1
+  do i=1,ex(1)-1
+   if(i >= 3 .and. i <= imax-2 .and. &
+      j >= 3 .and. j <= jmax-2 .and. &
+      k >= 3 .and. k <= kmax-2) cycle
   if(i+2 <= imax .and. i-2 >= imin .and. &
      j+2 <= jmax .and. j-2 >= jmin .and. &
      k+2 <= kmax .and. k-2 >= kmin) then
@@ -1518,10 +1585,10 @@
   fxz(i,j,k) = Sdxdz*(fh(i-1,j,k-1)-fh(i+1,j,k-1)-fh(i-1,j,k+1)+fh(i+1,j,k+1))
   fyz(i,j,k) = Sdydz*(fh(i,j-1,k-1)-fh(i,j+1,k-1)-fh(i,j-1,k+1)+fh(i,j+1,k+1))
   endif
+  enddo
+  enddo
+  enddo
 #endif
-   enddo
-   enddo
-   enddo

  return

--- a/AMSS_NCKU_source/enforce_algebra.f90
+++ b/AMSS_NCKU_source/enforce_algebra.f90
@@ -19,48 +19,60 @@

 !~~~~~~~> Local variable:

-  real*8, dimension(ex(1),ex(2),ex(3)) :: trA,detg
-  real*8, dimension(ex(1),ex(2),ex(3)) :: gxx,gyy,gzz 
-  real*8, dimension(ex(1),ex(2),ex(3)) :: gupxx,gupxy,gupxz,gupyy,gupyz,gupzz
+  integer :: i,j,k
+  real*8 :: lgxx,lgyy,lgzz,ldetg
+  real*8 :: lgupxx,lgupxy,lgupxz,lgupyy,lgupyz,lgupzz
+  real*8 :: ltrA,lscale
  real*8, parameter :: F1o3 = 1.D0 / 3.D0, ONE = 1.D0, TWO = 2.D0

 !~~~~~~>

-  gxx = dxx + ONE
-  gyy = dyy + ONE
-  gzz = dzz + ONE
+  do k=1,ex(3)
+  do j=1,ex(2)
+  do i=1,ex(1)

-  detg =  gxx * gyy * gzz + gxy * gyz * gxz + gxz * gxy * gyz - &
-          gxz * gyy * gxz - gxy * gxy * gzz - gxx * gyz * gyz
-  gupxx =   ( gyy * gzz - gyz * gyz ) / detg
-  gupxy = - ( gxy * gzz - gyz * gxz ) / detg
-  gupxz =   ( gxy * gyz - gyy * gxz ) / detg
-  gupyy =   ( gxx * gzz - gxz * gxz ) / detg
-  gupyz = - ( gxx * gyz - gxy * gxz ) / detg
-  gupzz =   ( gxx * gyy - gxy * gxy ) / detg
+    lgxx = dxx(i,j,k) + ONE
+    lgyy = dyy(i,j,k) + ONE
+    lgzz = dzz(i,j,k) + ONE

-  trA =         gupxx * Axx + gupyy * Ayy + gupzz * Azz &
-       + TWO * (gupxy * Axy + gupxz * Axz + gupyz * Ayz)
+    ldetg =  lgxx * lgyy * lgzz &
+           + gxy(i,j,k) * gyz(i,j,k) * gxz(i,j,k) &
+           + gxz(i,j,k) * gxy(i,j,k) * gyz(i,j,k) &
+           - gxz(i,j,k) * lgyy * gxz(i,j,k) &
+           - gxy(i,j,k) * gxy(i,j,k) * lgzz &
+           - lgxx * gyz(i,j,k) * gyz(i,j,k)

-  Axx = Axx - F1o3 * gxx * trA
-  Axy = Axy - F1o3 * gxy * trA
-  Axz = Axz - F1o3 * gxz * trA
-  Ayy = Ayy - F1o3 * gyy * trA
-  Ayz = Ayz - F1o3 * gyz * trA
-  Azz = Azz - F1o3 * gzz * trA
+    lgupxx =   ( lgyy * lgzz - gyz(i,j,k) * gyz(i,j,k) ) / ldetg
+    lgupxy = - ( gxy(i,j,k) * lgzz - gyz(i,j,k) * gxz(i,j,k) ) / ldetg
+    lgupxz =   ( gxy(i,j,k) * gyz(i,j,k) - lgyy * gxz(i,j,k) ) / ldetg
+    lgupyy =   ( lgxx * lgzz - gxz(i,j,k) * gxz(i,j,k) ) / ldetg
+    lgupyz = - ( lgxx * gyz(i,j,k) - gxy(i,j,k) * gxz(i,j,k) ) / ldetg
+    lgupzz =   ( lgxx * lgyy - gxy(i,j,k) * gxy(i,j,k) ) / ldetg

-  detg = ONE / ( detg ** F1o3 ) 
+    ltrA =         lgupxx * Axx(i,j,k) + lgupyy * Ayy(i,j,k) &
+                 + lgupzz * Azz(i,j,k) &
+         + TWO * (lgupxy * Axy(i,j,k) + lgupxz * Axz(i,j,k) &
+                 + lgupyz * Ayz(i,j,k))

-  gxx = gxx * detg
-  gxy = gxy * detg
-  gxz = gxz * detg
-  gyy = gyy * detg
-  gyz = gyz * detg
-  gzz = gzz * detg
+    Axx(i,j,k) = Axx(i,j,k) - F1o3 * lgxx * ltrA
+    Axy(i,j,k) = Axy(i,j,k) - F1o3 * gxy(i,j,k) * ltrA
+    Axz(i,j,k) = Axz(i,j,k) - F1o3 * gxz(i,j,k) * ltrA
+    Ayy(i,j,k) = Ayy(i,j,k) - F1o3 * lgyy * ltrA
+    Ayz(i,j,k) = Ayz(i,j,k) - F1o3 * gyz(i,j,k) * ltrA
+    Azz(i,j,k) = Azz(i,j,k) - F1o3 * lgzz * ltrA

-  dxx = gxx - ONE
-  dyy = gyy - ONE
-  dzz = gzz - ONE
+    lscale = ONE / ( ldetg ** F1o3 )
+
+    dxx(i,j,k) = lgxx * lscale - ONE
+    gxy(i,j,k) = gxy(i,j,k) * lscale
+    gxz(i,j,k) = gxz(i,j,k) * lscale
+    dyy(i,j,k) = lgyy * lscale - ONE
+    gyz(i,j,k) = gyz(i,j,k) * lscale
+    dzz(i,j,k) = lgzz * lscale - ONE
+
+  enddo
+  enddo
+  enddo

  return

@@ -83,50 +95,70 @@

 !~~~~~~~> Local variable:

-  real*8, dimension(ex(1),ex(2),ex(3)) :: trA
-  real*8, dimension(ex(1),ex(2),ex(3)) :: gxx,gyy,gzz 
-  real*8, dimension(ex(1),ex(2),ex(3)) :: gupxx,gupxy,gupxz,gupyy,gupyz,gupzz
+  integer :: i,j,k
+  real*8 :: lgxx,lgyy,lgzz,lscale
+  real*8 :: lgxy,lgxz,lgyz
+  real*8 :: lgupxx,lgupxy,lgupxz,lgupyy,lgupyz,lgupzz
+  real*8 :: ltrA
  real*8, parameter :: F1o3 = 1.D0 / 3.D0, ONE = 1.D0, TWO = 2.D0

 !~~~~~~>

-  gxx = dxx + ONE
-  gyy = dyy + ONE
-  gzz = dzz + ONE
-! for g
-  gupzz =  gxx * gyy * gzz + gxy * gyz * gxz + gxz * gxy * gyz - &
-           gxz * gyy * gxz - gxy * gxy * gzz - gxx * gyz * gyz
+  do k=1,ex(3)
+  do j=1,ex(2)
+  do i=1,ex(1)

-  gupzz = ONE / ( gupzz ** F1o3 ) 
+! for g: normalize determinant first
+    lgxx = dxx(i,j,k) + ONE
+    lgyy = dyy(i,j,k) + ONE
+    lgzz = dzz(i,j,k) + ONE
+    lgxy = gxy(i,j,k)
+    lgxz = gxz(i,j,k)
+    lgyz = gyz(i,j,k)

-  gxx = gxx * gupzz
-  gxy = gxy * gupzz
-  gxz = gxz * gupzz
-  gyy = gyy * gupzz
-  gyz = gyz * gupzz
-  gzz = gzz * gupzz
+    lscale =  lgxx * lgyy * lgzz + lgxy * lgyz * lgxz &
+            + lgxz * lgxy * lgyz - lgxz * lgyy * lgxz &
+            - lgxy * lgxy * lgzz - lgxx * lgyz * lgyz

-  dxx = gxx - ONE
-  dyy = gyy - ONE
-  dzz = gzz - ONE
-! for A  
+    lscale = ONE / ( lscale ** F1o3 )

-  gupxx =   ( gyy * gzz - gyz * gyz )
-  gupxy = - ( gxy * gzz - gyz * gxz )
-  gupxz =   ( gxy * gyz - gyy * gxz )
-  gupyy =   ( gxx * gzz - gxz * gxz )
-  gupyz = - ( gxx * gyz - gxy * gxz )
-  gupzz =   ( gxx * gyy - gxy * gxy )
+    lgxx = lgxx * lscale
+    lgxy = lgxy * lscale
+    lgxz = lgxz * lscale
+    lgyy = lgyy * lscale
+    lgyz = lgyz * lscale
+    lgzz = lgzz * lscale

-  trA =         gupxx * Axx + gupyy * Ayy + gupzz * Azz &
-       + TWO * (gupxy * Axy + gupxz * Axz + gupyz * Ayz)
+    dxx(i,j,k) = lgxx - ONE
+    gxy(i,j,k) = lgxy
+    gxz(i,j,k) = lgxz
+    dyy(i,j,k) = lgyy - ONE
+    gyz(i,j,k) = lgyz
+    dzz(i,j,k) = lgzz - ONE

-  Axx = Axx - F1o3 * gxx * trA
-  Axy = Axy - F1o3 * gxy * trA
-  Axz = Axz - F1o3 * gxz * trA
-  Ayy = Ayy - F1o3 * gyy * trA
-  Ayz = Ayz - F1o3 * gyz * trA
-  Azz = Azz - F1o3 * gzz * trA
+! for A: trace-free using normalized metric (det=1, no division needed)
+    lgupxx =   ( lgyy * lgzz - lgyz * lgyz )
+    lgupxy = - ( lgxy * lgzz - lgyz * lgxz )
+    lgupxz =   ( lgxy * lgyz - lgyy * lgxz )
+    lgupyy =   ( lgxx * lgzz - lgxz * lgxz )
+    lgupyz = - ( lgxx * lgyz - lgxy * lgxz )
+    lgupzz =   ( lgxx * lgyy - lgxy * lgxy )
+
+    ltrA =         lgupxx * Axx(i,j,k) + lgupyy * Ayy(i,j,k) &
+                 + lgupzz * Azz(i,j,k) &
+         + TWO * (lgupxy * Axy(i,j,k) + lgupxz * Axz(i,j,k) &
+                 + lgupyz * Ayz(i,j,k))
+
+    Axx(i,j,k) = Axx(i,j,k) - F1o3 * lgxx * ltrA
+    Axy(i,j,k) = Axy(i,j,k) - F1o3 * lgxy * ltrA
+    Axz(i,j,k) = Axz(i,j,k) - F1o3 * lgxz * ltrA
+    Ayy(i,j,k) = Ayy(i,j,k) - F1o3 * lgyy * ltrA
+    Ayz(i,j,k) = Ayz(i,j,k) - F1o3 * lgyz * ltrA
+    Azz(i,j,k) = Azz(i,j,k) - F1o3 * lgzz * ltrA
+
+  enddo
+  enddo
+  enddo

  return

--- a/AMSS_NCKU_source/fmisc.f90
+++ b/AMSS_NCKU_source/fmisc.f90
@@ -324,7 +324,6 @@ subroutine symmetry_bd(ord,extc,func,funcc,SoA)

  integer::i

-  funcc = 0.d0
  funcc(1:extc(1),1:extc(2),1:extc(3)) = func
   do i=0,ord-1
      funcc(-i,1:extc(2),1:extc(3)) = funcc(i+2,1:extc(2),1:extc(3))*SoA(1)
@@ -350,7 +349,6 @@ subroutine symmetry_tbd(ord,extc,func,funcc,SoA)

  integer::i

-  funcc = 0.d0
  funcc(1:extc(1),1:extc(2),1:extc(3)) = func
   do i=0,ord-1
      funcc(-i,1:extc(2),1:extc(3)) = funcc(i+2,1:extc(2),1:extc(3))*SoA(1)
@@ -379,7 +377,6 @@ subroutine symmetry_stbd(ord,extc,func,funcc,SoA)

  integer::i

-  funcc = 0.d0
  funcc(1:extc(1),1:extc(2),1:extc(3)) = func
   do i=0,ord-1
      funcc(-i,1:extc(2),1:extc(3)) = funcc(i+2,1:extc(2),1:extc(3))*SoA(1)
@@ -884,10 +881,18 @@ subroutine symmetry_bd(ord,extc,func,funcc,SoA)
  real*8, dimension(-ord+1:extc(1),-ord+1:extc(2),-ord+1:extc(3)),intent(out):: funcc
  real*8, dimension(1:3), intent(in) :: SoA

-  integer::i
+  integer::i,j,k
+
+  !$OMP PARALLEL DO COLLAPSE(2) SCHEDULE(static) PRIVATE(i,j,k)
+  do k=1,extc(3)
+  do j=1,extc(2)
+  do i=1,extc(1)
+     funcc(i,j,k) = func(i,j,k)
+  enddo
+  enddo
+  enddo
+  !$OMP END PARALLEL DO

-  funcc = 0.d0
-  funcc(1:extc(1),1:extc(2),1:extc(3)) = func
   do i=0,ord-1
      funcc(-i,1:extc(2),1:extc(3)) = funcc(i+1,1:extc(2),1:extc(3))*SoA(1)
   enddo
@@ -912,7 +917,6 @@ subroutine symmetry_tbd(ord,extc,func,funcc,SoA)

  integer::i

-  funcc = 0.d0
  funcc(1:extc(1),1:extc(2),1:extc(3)) = func
   do i=0,ord-1
      funcc(-i,1:extc(2),1:extc(3)) = funcc(i+1,1:extc(2),1:extc(3))*SoA(1)
@@ -941,7 +945,6 @@ subroutine symmetry_stbd(ord,extc,func,funcc,SoA)

  integer::i

-  funcc = 0.d0
  funcc(1:extc(1),1:extc(2),1:extc(3)) = func
   do i=0,ord-1
      funcc(-i,1:extc(2),1:extc(3)) = funcc(i+1,1:extc(2),1:extc(3))*SoA(1)
@@ -1119,25 +1122,16 @@ end subroutine d2dump
 !------------------------------------------------------------------------------

  subroutine polint(xa, ya, x, y, dy, ordn)
-
  implicit none

-!~~~~~~> Input Parameter:
  integer, intent(in) :: ordn
  real*8, dimension(ordn), intent(in) :: xa, ya
  real*8, intent(in) :: x
  real*8, intent(out) :: y, dy

-!~~~~~~> Other parameter:
-
-  integer :: m,n,ns
-  real*8, dimension(ordn) :: c,d,den,ho
-  real*8 :: dif,dift
-
-!~~~~~~>
-
-  n=ordn
-  m=ordn
+  integer :: i, m, ns, n_m
+  real*8, dimension(ordn) :: c, d, ho
+  real*8 :: dif, dift, hp, h, den_val

  c = ya
  d = ya
@@ -1145,27 +1139,38 @@ end subroutine d2dump

  ns = 1
  dif = abs(x - xa(1))
-  do m=1,n
-   dift=abs(x-xa(m))
+
+  do i = 2, ordn
+    dift = abs(x - xa(i))
    if (dift < dif) then
-    ns=m
+      ns = i
      dif = dift
    end if
  end do

  y = ya(ns)
  ns = ns - 1
-  do m=1,n-1
-    den(1:n-m)=ho(1:n-m)-ho(1+m:n)
-    if (any(den(1:n-m) == 0.0))then
+
+  do m = 1, ordn - 1
+    n_m = ordn - 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(1:n-m)=(c(2:n-m+1)-d(1:n-m))/den(1:n-m)
-    d(1:n-m)=ho(1+m:n)*den(1:n-m)
-    c(1:n-m)=ho(1:n-m)*den(1:n-m)
-    if (2*ns < n-m) then
+
+      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)
@@ -1175,7 +1180,6 @@ end subroutine d2dump
  end do

  return
-
  end subroutine polint
 !------------------------------------------------------------------------------
 !
@@ -1183,35 +1187,37 @@ end subroutine d2dump
 !
 !------------------------------------------------------------------------------
  subroutine polin2(x1a,x2a,ya,x1,x2,y,dy,ordn)
-
  implicit none

-!~~~~~~> Input parameters:
  integer,intent(in) :: ordn
  real*8, dimension(1:ordn), intent(in) :: x1a,x2a
  real*8, dimension(1:ordn,1:ordn), intent(in) :: ya
  real*8, intent(in) :: x1,x2
  real*8, intent(out) :: y,dy

-!~~~~~~> Other parameters:
-
+#ifdef POLINT_LEGACY_ORDER
  integer  :: i,m
  real*8, dimension(ordn) :: ymtmp
  real*8, dimension(ordn) :: yntmp

  m=size(x1a)
-  
  do i=1,m
-
    yntmp=ya(i,:)
    call polint(x2a,yntmp,x2,ymtmp(i),dy,ordn)
-
  end do
-
  call polint(x1a,ymtmp,x1,y,dy,ordn)
+#else
+  integer  :: j
+  real*8, dimension(ordn) :: ymtmp
+  real*8 :: dy_temp
+
+  do j=1,ordn
+    call polint(x1a, ya(:,j), x1, ymtmp(j), dy_temp, ordn)
+  end do
+  call polint(x2a, ymtmp, x2, y, dy, ordn)
+#endif

  return
-
  end subroutine polin2
 !------------------------------------------------------------------------------
 !
@@ -1219,18 +1225,15 @@ end subroutine d2dump
 !
 !------------------------------------------------------------------------------
  subroutine polin3(x1a,x2a,x3a,ya,x1,x2,x3,y,dy,ordn)
-
  implicit none

-!~~~~~~> Input parameters:
  integer,intent(in) :: ordn
  real*8, dimension(1:ordn), intent(in) :: x1a,x2a,x3a
  real*8, dimension(1:ordn,1:ordn,1:ordn), intent(in) :: ya
  real*8, intent(in) :: x1,x2,x3
  real*8, intent(out) :: y,dy

-!~~~~~~> Other parameters:
-
+#ifdef POLINT_LEGACY_ORDER
  integer  :: i,j,m,n
  real*8, dimension(ordn,ordn) :: yatmp
  real*8, dimension(ordn) :: ymtmp
@@ -1239,24 +1242,33 @@ end subroutine d2dump

  m=size(x1a)
  n=size(x2a)
-  
  do i=1,m
   do j=1,n
-
    yqtmp=ya(i,j,:)
    call polint(x3a,yqtmp,x3,yatmp(i,j),dy,ordn)
-
   end do
-
    yntmp=yatmp(i,:)
    call polint(x2a,yntmp,x2,ymtmp(i),dy,ordn)
-
  end do
-
  call polint(x1a,ymtmp,x1,y,dy,ordn)
+#else
+  integer  :: j, k
+  real*8, dimension(ordn,ordn) :: yatmp
+  real*8, dimension(ordn) :: ymtmp
+  real*8 :: dy_temp
+
+  do k=1,ordn
+    do j=1,ordn
+      call polint(x1a, ya(:,j,k), x1, yatmp(j,k), dy_temp, ordn)
+    end do
+  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

  return
-
  end subroutine polin3
 !--------------------------------------------------------------------------------------
 ! calculate L2norm
--- a/AMSS_NCKU_source/kodiss.f90
+++ b/AMSS_NCKU_source/kodiss.f90
@@ -159,36 +159,12 @@ integer, parameter :: NO_SYMM=0, OCTANT=2

  call symmetry_bd(3,ex,f,fh,SoA)

-  do k=1,ex(3)
-  do j=1,ex(2)
-  do i=1,ex(1)
-
-  if(i-3 >= imin .and. i+3 <= imax .and. &
-     j-3 >= jmin .and. j+3 <= jmax .and. &
-     k-3 >= kmin .and. k+3 <= kmax) then
-#if 0     
-! x direction
-   f_rhs(i,j,k)       = f_rhs(i,j,k) + eps/dX/cof * (     &
-                              (fh(i-3,j,k)+fh(i+3,j,k)) - &
-                          SIX*(fh(i-2,j,k)+fh(i+2,j,k)) + &
-                          FIT*(fh(i-1,j,k)+fh(i+1,j,k)) - &
-                          TWT* fh(i,j,k)            )
-! y direction
-
-   f_rhs(i,j,k)       = f_rhs(i,j,k) + eps/dY/cof * (     &
-                              (fh(i,j-3,k)+fh(i,j+3,k)) - &
-                          SIX*(fh(i,j-2,k)+fh(i,j+2,k)) + &
-                          FIT*(fh(i,j-1,k)+fh(i,j+1,k)) - &
-                          TWT* fh(i,j,k)            )
-! z direction
-
-   f_rhs(i,j,k)       = f_rhs(i,j,k) + eps/dZ/cof * (     &
-                              (fh(i,j,k-3)+fh(i,j,k+3)) - &
-                          SIX*(fh(i,j,k-2)+fh(i,j,k+2)) + &
-                          FIT*(fh(i,j,k-1)+fh(i,j,k+1)) - &
-                          TWT* fh(i,j,k)            )
-#else
-! calculation order if important ?
+! Interior: all stencil points guaranteed in-bounds
+  !$OMP PARALLEL DO COLLAPSE(2) SCHEDULE(static) PRIVATE(i,j,k)
+  do k=4,ex(3)-3
+  do j=4,ex(2)-3
+  !DIR$ IVDEP
+  do i=4,ex(1)-3
   f_rhs(i,j,k)       = f_rhs(i,j,k) + eps/cof *( (     &
                              (fh(i-3,j,k)+fh(i+3,j,k)) - &
                          SIX*(fh(i-2,j,k)+fh(i+2,j,k)) + &
@@ -204,9 +180,37 @@ integer, parameter :: NO_SYMM=0, OCTANT=2
                          SIX*(fh(i,j,k-2)+fh(i,j,k+2)) + &
                          FIT*(fh(i,j,k-1)+fh(i,j,k+1)) - &
                          TWT* fh(i,j,k)            )/dZ )
-#endif
-  endif
+  enddo
+  enddo
+  enddo
+  !$OMP END PARALLEL DO

+! Boundary shell: original branching logic for points near edges
+  do k=1,ex(3)
+  do j=1,ex(2)
+  do i=1,ex(1)
+  if(i >= 4 .and. i <= ex(1)-3 .and. &
+     j >= 4 .and. j <= ex(2)-3 .and. &
+     k >= 4 .and. k <= ex(3)-3) cycle
+  if(i-3 >= imin .and. i+3 <= imax .and. &
+     j-3 >= jmin .and. j+3 <= jmax .and. &
+     k-3 >= kmin .and. k+3 <= kmax) then
+   f_rhs(i,j,k)       = f_rhs(i,j,k) + eps/cof *( (     &
+                              (fh(i-3,j,k)+fh(i+3,j,k)) - &
+                          SIX*(fh(i-2,j,k)+fh(i+2,j,k)) + &
+                          FIT*(fh(i-1,j,k)+fh(i+1,j,k)) - &
+                          TWT* fh(i,j,k)            )/dX + &
+                                                  (     &
+                              (fh(i,j-3,k)+fh(i,j+3,k)) - &
+                          SIX*(fh(i,j-2,k)+fh(i,j+2,k)) + &
+                          FIT*(fh(i,j-1,k)+fh(i,j+1,k)) - &
+                          TWT* fh(i,j,k)            )/dY + &
+                                                  (     &
+                              (fh(i,j,k-3)+fh(i,j,k+3)) - &
+                          SIX*(fh(i,j,k-2)+fh(i,j,k+2)) + &
+                          FIT*(fh(i,j,k-1)+fh(i,j,k+1)) - &
+                          TWT* fh(i,j,k)            )/dZ )
+  endif
  enddo
  enddo
  enddo
--- a/AMSS_NCKU_source/lopsidediff.f90
+++ b/AMSS_NCKU_source/lopsidediff.f90
@@ -233,6 +233,7 @@ subroutine lopsided(ex,X,Y,Z,f,f_rhs,Sfx,Sfy,Sfz,Symmetry,SoA)

 ! upper bound set ex-1 only for efficiency,
 ! the loop body will set ex 0 also
+  !$OMP PARALLEL DO COLLAPSE(2) SCHEDULE(static) PRIVATE(i,j,k)
  do k=1,ex(3)-1
  do j=1,ex(2)-1
  do i=1,ex(1)-1
@@ -482,6 +483,7 @@ subroutine lopsided(ex,X,Y,Z,f,f_rhs,Sfx,Sfy,Sfz,Symmetry,SoA)
  enddo
  enddo
  enddo
+  !$OMP END PARALLEL DO

  return

--- a/AMSS_NCKU_source/makefile.inc
+++ b/AMSS_NCKU_source/makefile.inc
@@ -8,7 +8,7 @@ filein  = -I/usr/include/ -I${MKLROOT}/include

 ## Using sequential MKL (OpenMP disabled for better single-threaded performance)
 ## 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
+LDLIBS  = -L${MKLROOT}/lib -lmkl_intel_lp64 -lmkl_intel_thread -lmkl_core -lifcore -limf -lpthread -lm -ldl -qopenmp

 ## Aggressive optimization flags:
 ## -O3: Maximum optimization
@@ -16,10 +16,10 @@ LDLIBS  = -L${MKLROOT}/lib -lmkl_intel_lp64 -lmkl_sequential -lmkl_core -lifcore
 ## -fp-model fast=2: Aggressive floating-point optimizations
 ## -fma: Enable fused multiply-add instructions
 ## Note: OpenMP has been disabled (-qopenmp removed) due to performance issues
-CXXAPPFLAGS  = -O3 -xHost -fp-model fast=2 -fma \
+CXXAPPFLAGS  = -O3 -xHost -fp-model fast=2 -fma -ipo -qopenmp \
               -Dfortran3 -Dnewc -I${MKLROOT}/include
-f90appflags  = -O3 -xHost -fp-model fast=2 -fma \
-               -fpp -I${MKLROOT}/include
+f90appflags  = -O3 -xHost -fp-model fast=2 -fma -ipo -qopenmp \
+               -align array64byte -fpp -I${MKLROOT}/include
 f90          = ifx
 f77          = ifx
 CXX          = icpx
Author	SHA1	Message	Date
CGH0S7	c6e4d4ab71	Add OpenMP parallelization to BSSN RHS hot-path stencil routines Enable OpenMP threading for the dominant computational kernels: - makefile.inc: add -qopenmp to f90appflags - diff_new.f90: split fderivs/fdderivs into OpenMP interior + serial boundary - kodiss.f90: split kodis into OpenMP interior + serial boundary - lopsidediff.f90: add OMP PARALLEL DO COLLAPSE(2) to lopsided - fmisc.f90: parallelize symmetry_bd bulk array copy - bssn_rhs.f90: add OMP WORKSHARE to array-syntax operations Co-Authored-By: Claude Opus 4.6 <noreply@anthropic.com>	2026-02-07 13:58:55 +08:00
CGH0S7	09ffdb553d	Eliminate hot-path heap allocations in TwoPunctures spectral solver Pre-allocate workspace buffers as class members to remove ~8M malloc/free pairs per Newton iteration from LineRelax, ThomasAlgorithm, JFD_times_dv, J_times_dv, chebft_Zeros, fourft, Derivatives_AB3, and F_of_v. Rewrite ThomasAlgorithm to operate in-place on input arrays. Results are bit-identical; no algorithmic changes. Co-Authored-By: Claude Opus 4.6 <noreply@anthropic.com>	2026-02-06 21:20:35 +08:00
CGH0S7	699e443c7a	Optimize polint/polin2/polin3 interpolation for cache locality Changes: - polint: Rewrite Neville algorithm from array-slice operations to scalar loop. Mathematically identical, avoids temporary array allocations for den(1:n-m) slices. (Credit: yx-fmisc branch) - polin2: Swap interpolation order so inner loop accesses ya(:,j) (contiguous in Fortran column-major) instead of ya(i,:) (strided). Tensor product interpolation is commutative; all call sites pass identical coordinate arrays for all dimensions. - polin3: Swap interpolation order to process contiguous first dimension first: ya(:,j,k) -> yatmp(:,k) -> ymtmp(:). Same commutativity argument as polin2. Compile-time safety switch: -DPOLINT_LEGACY_ORDER restores original dimension ordering Default (no flag): uses optimized contiguous-memory ordering Usage: # Production (optimized order): make clean && make -j ABE # Fallback if results differ (original order): Add -DPOLINT_LEGACY_ORDER to f90appflags in makefile.inc Co-Authored-By: Claude Opus 4.6 <noreply@anthropic.com>	2026-02-06 19:00:35 +08:00
CGH0S7	24bfa44911	Disable NaN sanity check in bssn_rhs.f90 for production builds Wrap the NaN sanity check (21 sum() full-array traversals per RHS call) with #ifdef DEBUG so it is compiled out in production builds. This eliminates 84 redundant full-array scans per timestep (21 per RHS call × 4 RK4 substages) that serve no purpose when input data is valid. Usage: - Production build (default): NaN check is disabled, no changes needed. - Debug build: add -DDEBUG to f90appflags in makefile.inc, e.g. f90appflags = -O3 ... -DDEBUG -fpp ... to re-enable the NaN sanity check. Co-Authored-By: Claude Opus 4.6 <noreply@anthropic.com>	2026-02-06 18:36:29 +08:00
CGH0S7	6738854a9d	Compiler-level and hot-path optimizations for GW150914 - makefile.inc: add -ipo (interprocedural optimization) and -align array64byte (64-byte array alignment for vectorization) - fmisc.f90: remove redundant funcc=0.d0 zeroing from symmetry_bd, symmetry_tbd, symmetry_stbd (~328+ full-array memsets eliminated per timestep) - enforce_algebra.f90: rewrite enforce_ag and enforce_ga as point-wise loops, replacing 12 stack-allocated 3D temporary arrays with scalar locals for better cache locality All changes are mathematically equivalent — no algorithmic modifications. Co-Authored-By: Claude Opus 4.6 <noreply@anthropic.com>	2026-02-06 17:13:39 +08:00