Trigger-Discipline: parallelize result plotting

.gitignore

@@ -1,2 +1,5 @@
 __pycache__
 GW150914
+GW150914*
+.codex
+docs/

AMSS_NCKU_Program.py

@@ -9,6 +9,16 @@
 ##################################################################
 
 
+##################################################################
+
+## Guard against re-execution by multiprocessing child processes.
+## Without this, using 'spawn' or 'forkserver' context would cause every
+## worker to re-run the entire script.
+if __name__ != '__main__':
+    import sys as _sys
+    _sys.exit(0)
+
+
 ##################################################################
 
 ## Print program introduction
@@ -424,26 +434,31 @@ print(
 
 import plot_xiaoqu
 import plot_GW_strain_amplitude_xiaoqu
+from parallel_plot_helper import run_plot_tasks_parallel
+
+plot_tasks = []
 
 ## Plot black hole trajectory
-plot_xiaoqu.generate_puncture_orbit_plot(   binary_results_directory, figure_directory )
+plot_tasks.append( ( plot_xiaoqu.generate_puncture_orbit_plot,   (binary_results_directory, figure_directory) ) )
-plot_xiaoqu.generate_puncture_orbit_plot3D( binary_results_directory, figure_directory )
+plot_tasks.append( ( plot_xiaoqu.generate_puncture_orbit_plot3D, (binary_results_directory, figure_directory) ) )
 
 ## Plot black hole separation vs. time
-plot_xiaoqu.generate_puncture_distence_plot( binary_results_directory, figure_directory )
+plot_tasks.append( ( plot_xiaoqu.generate_puncture_distence_plot, (binary_results_directory, figure_directory) ) )
 
 ## Plot gravitational waveforms (psi4 and strain amplitude)
 for i in range(input_data.Detector_Number):
-    plot_xiaoqu.generate_gravitational_wave_psi4_plot( binary_results_directory, figure_directory, i )
+    plot_tasks.append( ( plot_xiaoqu.generate_gravitational_wave_psi4_plot, (binary_results_directory, figure_directory, i) ) )
-    plot_GW_strain_amplitude_xiaoqu.generate_gravitational_wave_amplitude_plot( binary_results_directory, figure_directory, i )
+    plot_tasks.append( ( plot_GW_strain_amplitude_xiaoqu.generate_gravitational_wave_amplitude_plot, (binary_results_directory, figure_directory, i) ) )
 
 ## Plot ADM mass evolution
 for i in range(input_data.Detector_Number):
-    plot_xiaoqu.generate_ADMmass_plot( binary_results_directory, figure_directory, i )
+    plot_tasks.append( ( plot_xiaoqu.generate_ADMmass_plot, (binary_results_directory, figure_directory, i) ) )
 
 ## Plot Hamiltonian constraint violation over time
 for i in range(input_data.grid_level):
-    plot_xiaoqu.generate_constraint_check_plot( binary_results_directory, figure_directory, i )
+    plot_tasks.append( ( plot_xiaoqu.generate_constraint_check_plot, (binary_results_directory, figure_directory, i) ) )
+
+run_plot_tasks_parallel(plot_tasks)
 
 ## Plot stored binary data
 plot_xiaoqu.generate_binary_data_plot( binary_results_directory, figure_directory )

AMSS_NCKU_source/TwoPunctures.h

@@ -1,7 +1,8 @@
-
 #ifndef TWO_PUNCTURES_H
 #define TWO_PUNCTURES_H
 
+#include <omp.h>
+
 #define StencilSize 19
 #define N_PlaneRelax 1
 #define NRELAX 200
@@ -42,6 +43,18 @@ private:
 
        int ntotal;
 
+       // ===== Precomputed spectral derivative matrices =====
+       double *D1_A, *D2_A;
+       double *D1_B, *D2_B;
+       double *DF1_phi, *DF2_phi;
+
+       // ===== Pre-allocated workspace for LineRelax (per-thread) =====
+       int max_threads;
+       double **ws_diag_be, **ws_e_be, **ws_f_be, **ws_b_be, **ws_x_be;
+       double **ws_l_be, **ws_u_be, **ws_d_be, **ws_y_be;
+       double **ws_diag_al, **ws_e_al, **ws_f_al, **ws_b_al, **ws_x_al;
+       double **ws_l_al, **ws_u_al, **ws_d_al, **ws_y_al;
+
        struct parameters
        {
               int nvar, n1, n2, n3;
@@ -58,6 +71,28 @@ public:
                     int Newtonmaxit);
        ~TwoPunctures();
 
+       // 02/07: New/modified methods
+       void allocate_workspace();
+       void free_workspace();
+       void precompute_derivative_matrices();
+       void build_cheb_deriv_matrices(int n, double *D1, double *D2);
+       void build_fourier_deriv_matrices(int N, double *DF1, double *DF2);
+       void Derivatives_AB3_MatMul(int nvar, int n1, int n2, int n3, derivs v);
+       void ThomasAlgorithm_ws(int N, double *b, double *a, double *c, double *x, double *q,
+                                double *l, double *u_ws, double *d, double *y);
+       void LineRelax_be_omp(double *dv,
+                         int const i, int const k, int const nvar,
+                         int const n1, int const n2, int const n3,
+                         double const *rhs, int const *ncols, int **cols,
+                         double **JFD, int tid);
+       void LineRelax_al_omp(double *dv,
+                         int const j, int const k, int const nvar,
+                         int const n1, int const n2, int const n3,
+                         double const *rhs, int const *ncols,
+                         int **cols, double **JFD, int tid);
+       void relax_omp(double *dv, int const nvar, int const n1, int const n2, int const n3,
+                  double const *rhs, int const *ncols, int **cols, double **JFD);
+
        void Solve();
        void set_initial_guess(derivs v);
        int index(int i, int j, int k, int l, int a, int b, int c, int d);
@@ -116,23 +151,11 @@ public:
        double BY_KKofxyz(double x, double y, double z);
        void SetMatrix_JFD(int nvar, int n1, int n2, int n3, derivs u, int *ncols, int **cols, double **Matrix);
        void J_times_dv(int nvar, int n1, int n2, int n3, derivs dv, double *Jdv, derivs u);
-       void relax(double *dv, int const nvar, int const n1, int const n2, int const n3,
-                  double const *rhs, int const *ncols, int **cols, double **JFD);
-       void LineRelax_be(double *dv,
-                         int const i, int const k, int const nvar,
-                         int const n1, int const n2, int const n3,
-                         double const *rhs, int const *ncols, int **cols,
-                         double **JFD);
        void JFD_times_dv(int i, int j, int k, int nvar, int n1, int n2,
                          int n3, derivs dv, derivs u, double *values);
        void LinEquations(double A, double B, double X, double R,
                          double x, double r, double phi,
                          double y, double z, derivs dU, derivs U, double *values);
-       void LineRelax_al(double *dv,
-                         int const j, int const k, int const nvar,
-                         int const n1, int const n2, int const n3,
-                         double const *rhs, int const *ncols,
-                         int **cols, double **JFD);
        void ThomasAlgorithm(int N, double *b, double *a, double *c, double *x, double *q);
        void Save(char *fname);
        // provided by Vasileios Paschalidis (vpaschal@illinois.edu)

AMSS_NCKU_source/enforce_algebra.f90

@@ -19,48 +19,60 @@
 
 !~~~~~~~> Local variable:
 
-  real*8, dimension(ex(1),ex(2),ex(3)) :: trA,detg
+  integer :: i,j,k
-  real*8, dimension(ex(1),ex(2),ex(3)) :: gxx,gyy,gzz 
+  real*8 :: lgxx,lgyy,lgzz,ldetg
-  real*8, dimension(ex(1),ex(2),ex(3)) :: gupxx,gupxy,gupxz,gupyy,gupyz,gupzz
+  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
+  do k=1,ex(3)
-  gyy = dyy + ONE
+  do j=1,ex(2)
-  gzz = dzz + ONE
+  do i=1,ex(1)
 
-  detg =  gxx * gyy * gzz + gxy * gyz * gxz + gxz * gxy * gyz - &
+    lgxx = dxx(i,j,k) + ONE
-          gxz * gyy * gxz - gxy * gxy * gzz - gxx * gyz * gyz
+    lgyy = dyy(i,j,k) + ONE
-  gupxx =   ( gyy * gzz - gyz * gyz ) / detg
+    lgzz = dzz(i,j,k) + ONE
-  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
 
-  trA =         gupxx * Axx + gupyy * Ayy + gupzz * Azz &
+    ldetg =  lgxx * lgyy * lgzz &
-       + TWO * (gupxy * Axy + gupxz * Axz + gupyz * Ayz)
+           + 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
+    lgupxx =   ( lgyy * lgzz - gyz(i,j,k) * gyz(i,j,k) ) / ldetg
-  Axy = Axy - F1o3 * gxy * trA
+    lgupxy = - ( gxy(i,j,k) * lgzz - gyz(i,j,k) * gxz(i,j,k) ) / ldetg
-  Axz = Axz - F1o3 * gxz * trA
+    lgupxz =   ( gxy(i,j,k) * gyz(i,j,k) - lgyy * gxz(i,j,k) ) / ldetg
-  Ayy = Ayy - F1o3 * gyy * trA
+    lgupyy =   ( lgxx * lgzz - gxz(i,j,k) * gxz(i,j,k) ) / ldetg
-  Ayz = Ayz - F1o3 * gyz * trA
+    lgupyz = - ( lgxx * gyz(i,j,k) - gxy(i,j,k) * gxz(i,j,k) ) / ldetg
-  Azz = Azz - F1o3 * gzz * trA
+    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
+    Axx(i,j,k) = Axx(i,j,k) - F1o3 * lgxx * ltrA
-  gxy = gxy * detg
+    Axy(i,j,k) = Axy(i,j,k) - F1o3 * gxy(i,j,k) * ltrA
-  gxz = gxz * detg
+    Axz(i,j,k) = Axz(i,j,k) - F1o3 * gxz(i,j,k) * ltrA
-  gyy = gyy * detg
+    Ayy(i,j,k) = Ayy(i,j,k) - F1o3 * lgyy * ltrA
-  gyz = gyz * detg
+    Ayz(i,j,k) = Ayz(i,j,k) - F1o3 * gyz(i,j,k) * ltrA
-  gzz = gzz * detg
+    Azz(i,j,k) = Azz(i,j,k) - F1o3 * lgzz * ltrA
 
-  dxx = gxx - ONE
+    lscale = ONE / ( ldetg ** F1o3 )
-  dyy = gyy - ONE
+
-  dzz = gzz - ONE
+    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
+  integer :: i,j,k
-  real*8, dimension(ex(1),ex(2),ex(3)) :: gxx,gyy,gzz 
+  real*8 :: lgxx,lgyy,lgzz,lscale
-  real*8, dimension(ex(1),ex(2),ex(3)) :: gupxx,gupxy,gupxz,gupyy,gupyz,gupzz
+  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
+  do k=1,ex(3)
-  gyy = dyy + ONE
+  do j=1,ex(2)
-  gzz = dzz + ONE
+  do i=1,ex(1)
-! for g
-  gupzz =  gxx * gyy * gzz + gxy * gyz * gxz + gxz * gxy * gyz - &
-           gxz * gyy * gxz - gxy * gxy * gzz - gxx * gyz * gyz
 
-  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
+    lscale =  lgxx * lgyy * lgzz + lgxy * lgyz * lgxz &
-  gxy = gxy * gupzz
+            + lgxz * lgxy * lgyz - lgxz * lgyy * lgxz &
-  gxz = gxz * gupzz
+            - lgxy * lgxy * lgzz - lgxx * lgyz * lgyz
-  gyy = gyy * gupzz
-  gyz = gyz * gupzz
-  gzz = gzz * gupzz
 
-  dxx = gxx - ONE
+    lscale = ONE / ( lscale ** F1o3 )
-  dyy = gyy - ONE
-  dzz = gzz - ONE
-! for A  
 
-  gupxx =   ( gyy * gzz - gyz * gyz )
+    lgxx = lgxx * lscale
-  gupxy = - ( gxy * gzz - gyz * gxz )
+    lgxy = lgxy * lscale
-  gupxz =   ( gxy * gyz - gyy * gxz )
+    lgxz = lgxz * lscale
-  gupyy =   ( gxx * gzz - gxz * gxz )
+    lgyy = lgyy * lscale
-  gupyz = - ( gxx * gyz - gxy * gxz )
+    lgyz = lgyz * lscale
-  gupzz =   ( gxx * gyy - gxy * gxy )
+    lgzz = lgzz * lscale
 
-  trA =         gupxx * Axx + gupyy * Ayy + gupzz * Azz &
+    dxx(i,j,k) = lgxx - ONE
-       + TWO * (gupxy * Axy + gupxz * Axz + gupyz * Ayz)
+    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
+! for A: trace-free using normalized metric (det=1, no division needed)
-  Axy = Axy - F1o3 * gxy * trA
+    lgupxx =   ( lgyy * lgzz - lgyz * lgyz )
-  Axz = Axz - F1o3 * gxz * trA
+    lgupxy = - ( lgxy * lgzz - lgyz * lgxz )
-  Ayy = Ayy - F1o3 * gyy * trA
+    lgupxz =   ( lgxy * lgyz - lgyy * lgxz )
-  Ayz = Ayz - F1o3 * gyz * trA
+    lgupyy =   ( lgxx * lgzz - lgxz * lgxz )
-  Azz = Azz - F1o3 * gzz * trA
+    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
 

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)
@@ -886,14 +883,17 @@ subroutine symmetry_bd(ord,extc,func,funcc,SoA)
 
   integer::i
 
-  funcc = 0.d0
+!DIR$ SIMD VECTORLENGTHFOR(KNOWN_INTEGER=8)
   funcc(1:extc(1),1:extc(2),1:extc(3)) = func
+!DIR$ SIMD VECTORLENGTHFOR(KNOWN_INTEGER=8)
    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
+!DIR$ SIMD VECTORLENGTHFOR(KNOWN_INTEGER=8)
    do i=0,ord-1
       funcc(:,-i,1:extc(3)) = funcc(:,i+1,1:extc(3))*SoA(2)
    enddo
+!DIR$ SIMD VECTORLENGTHFOR(KNOWN_INTEGER=8)
    do i=0,ord-1
       funcc(:,:,-i) = funcc(:,:,i+1)*SoA(3)
    enddo
@@ -912,7 +912,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 +940,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)
@@ -1117,27 +1115,21 @@ end subroutine d2dump
 !------------------------------------------------------------------------------
 ! Lagrangian polynomial interpolation
 !------------------------------------------------------------------------------
+#ifndef POLINT6_USE_BARYCENTRIC
+#define POLINT6_USE_BARYCENTRIC 1
+#endif
 
-  subroutine polint(xa,ya,x,y,dy,ordn)
+!DIR$ ATTRIBUTES FORCEINLINE :: polint6_neville
-
+  subroutine polint6_neville(xa, ya, x, y, dy)
   implicit none
 
-!~~~~~~> Input Parameter:
+  real*8, dimension(6), intent(in) :: xa, ya
-  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 :: i, m, ns, n_m
-
+  real*8, dimension(6) :: c, d, ho
-  integer :: m,n,ns
+  real*8 :: dif, dift, hp, h, den_val
-  real*8, dimension(ordn) :: c,d,den,ho
-  real*8 :: dif,dift
-
-!~~~~~~>
-
-  n=ordn
-  m=ordn
 
   c = ya
   d = ya
@@ -1145,27 +1137,38 @@ end subroutine d2dump
 
   ns = 1
   dif = abs(x - xa(1))
-  do m=1,n
+
-   dift=abs(x-xa(m))
+  do i = 2, 6
+    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)
+  do m = 1, 5
-    if (any(den(1:n-m) == 0.0))then
+    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(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)
+      den_val = (c(i+1) - d(i)) / den_val
-    c(1:n-m)=ho(1:n-m)*den(1:n-m)
+
-    if (2*ns < n-m) then
+      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,43 +1178,227 @@ end subroutine d2dump
   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)
+  implicit none
+
+  integer, intent(in) :: ordn
+  real*8, dimension(ordn), intent(in) :: xa, ya
+  real*8, intent(in) :: x
+  real*8, intent(out) :: y, dy
+
+  integer :: i, m, ns, n_m
+  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
+
+  ns = 1
+  dif = abs(x - xa(1))
+
+  do i = 2, ordn
+    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, 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_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 polint
 !------------------------------------------------------------------------------
+! 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
 !
 !------------------------------------------------------------------------------
   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 +1406,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 +1423,40 @@ 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  :: i, j, k
+  real*8, dimension(ordn) :: w1, w2
+  real*8, dimension(ordn) :: ymtmp
+  real*8 :: yx_sum, x_sum
 
+  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
 
-  call polint(x1a,ymtmp,x1,y,dy,ordn)
+  call polint(x3a, ymtmp, x3, y, dy, ordn)
+#endif
 
   return
-
   end subroutine polin3
 !--------------------------------------------------------------------------------------
 ! calculate L2norm  
@@ -1276,7 +1476,9 @@ end subroutine d2dump
   real*8            :: dX, dY, dZ
   integer::imin,jmin,kmin
   integer::imax,jmax,kmax
-  integer::i,j,k
+  integer::i,j,k,n_elements
+  real*8, dimension(:), allocatable :: f_flat
+  real*8, external :: DDOT
 
   dX = X(2) - X(1)
   dY = Y(2) - Y(1)
@@ -1300,13 +1502,89 @@ if(dabs(X(1)-xmin) < dX) imin = 1
 if(dabs(Y(1)-ymin) < dY) jmin = 1
 if(dabs(Z(1)-zmin) < dZ) kmin = 1
 
-f_out = sum(f(imin:imax,jmin:jmax,kmin:kmax)*f(imin:imax,jmin:jmax,kmin:kmax))
+  n_elements = (imax-imin+1)*(jmax-jmin+1)*(kmax-kmin+1)
+  allocate(f_flat(n_elements))
+  f_flat = reshape(f(imin:imax,jmin:jmax,kmin:kmax), [n_elements])
+  f_out = DDOT(n_elements, f_flat, 1, f_flat, 1)
+  deallocate(f_flat)
 
 f_out = f_out*dX*dY*dZ
 
   return
 
   end subroutine l2normhelper
+!--------------------------------------------------------------------------------------
+  subroutine l2normhelper7(ex, X, Y, Z,xmin,ymin,zmin,xmax,ymax,zmax,&
+                           f1,f2,f3,f4,f5,f6,f7,f_out,gw)
+
+  implicit none
+!~~~~~~> Input parameters:
+  integer,intent(in ):: ex(1:3)
+  real*8, intent(in ):: X(1:ex(1)),Y(1:ex(2)),Z(1:ex(3)),xmin,ymin,zmin,xmax,ymax,zmax
+  integer,intent(in)::gw
+  real*8, dimension(ex(1),ex(2),ex(3)),intent(in) :: f1,f2,f3,f4,f5,f6,f7
+  real*8, intent(out) :: f_out(7)
+!~~~~~~> Other variables:
+
+  real*8            :: dX, dY, dZ
+  integer::imin,jmin,kmin
+  integer::imax,jmax,kmax
+  integer::i,j,k
+  real*8 :: s1,s2,s3,s4,s5,s6,s7
+
+  dX = X(2) - X(1)
+  dY = Y(2) - Y(1)
+  dZ = Z(2) - Z(1)
+
+   imin = gw+1
+   jmin = gw+1
+   kmin = gw+1
+
+   imax = ex(1) - gw
+   jmax = ex(2) - gw
+   kmax = ex(3) - gw
+
+if(dabs(X(ex(1))-xmax) < dX) imax = ex(1)
+if(dabs(Y(ex(2))-ymax) < dY) jmax = ex(2)
+if(dabs(Z(ex(3))-zmax) < dZ) kmax = ex(3)
+if(dabs(X(1)-xmin) < dX) imin = 1
+if(dabs(Y(1)-ymin) < dY) jmin = 1
+if(dabs(Z(1)-zmin) < dZ) kmin = 1
+
+  s1 = 0.d0
+  s2 = 0.d0
+  s3 = 0.d0
+  s4 = 0.d0
+  s5 = 0.d0
+  s6 = 0.d0
+  s7 = 0.d0
+
+  do k=kmin,kmax
+    do j=jmin,jmax
+!DIR$ SIMD REDUCTION(+:s1,s2,s3,s4,s5,s6,s7)
+      do i=imin,imax
+        s1 = s1 + f1(i,j,k)*f1(i,j,k)
+        s2 = s2 + f2(i,j,k)*f2(i,j,k)
+        s3 = s3 + f3(i,j,k)*f3(i,j,k)
+        s4 = s4 + f4(i,j,k)*f4(i,j,k)
+        s5 = s5 + f5(i,j,k)*f5(i,j,k)
+        s6 = s6 + f6(i,j,k)*f6(i,j,k)
+        s7 = s7 + f7(i,j,k)*f7(i,j,k)
+      enddo
+    enddo
+  enddo
+
+  f_out(1) = s1*dX*dY*dZ
+  f_out(2) = s2*dX*dY*dZ
+  f_out(3) = s3*dX*dY*dZ
+  f_out(4) = s4*dX*dY*dZ
+  f_out(5) = s5*dX*dY*dZ
+  f_out(6) = s6*dX*dY*dZ
+  f_out(7) = s7*dX*dY*dZ
+
+  return
+
+  end subroutine l2normhelper7
 !--------------------------------------------------------------------------------------
 ! calculate L2norm especially for shell Blocks
   subroutine l2normhelper_sh(ex, X, Y, Z,xmin,ymin,zmin,xmax,ymax,zmax,&
@@ -1325,7 +1603,9 @@ f_out = f_out*dX*dY*dZ
   real*8            :: dX, dY, dZ
   integer::imin,jmin,kmin
   integer::imax,jmax,kmax
-  integer::i,j,k
+  integer::i,j,k,n_elements
+  real*8, dimension(:), allocatable :: f_flat
+  real*8, external :: DDOT
 
   real*8 :: PIo4
 
@@ -1388,7 +1668,11 @@ if(Symmetry==2)then
   if(dabs(ymin+gw*dY)<dY.and.Y(1)<0.d0) jmin = gw+1
 endif
 
-f_out = sum(f(imin:imax,jmin:jmax,kmin:kmax)*f(imin:imax,jmin:jmax,kmin:kmax))
+  n_elements = (imax-imin+1)*(jmax-jmin+1)*(kmax-kmin+1)
+  allocate(f_flat(n_elements))
+  f_flat = reshape(f(imin:imax,jmin:jmax,kmin:kmax), [n_elements])
+  f_out = DDOT(n_elements, f_flat, 1, f_flat, 1)
+  deallocate(f_flat)
 
 f_out = f_out*dX*dY*dZ
 
@@ -1416,6 +1700,8 @@ f_out = f_out*dX*dY*dZ
   integer::imin,jmin,kmin
   integer::imax,jmax,kmax
   integer::i,j,k
+  real*8, dimension(:), allocatable :: f_flat
+  real*8, external :: DDOT
 
   real*8 :: PIo4
 
@@ -1478,11 +1764,11 @@ if(Symmetry==2)then
   if(dabs(ymin+gw*dY)<dY.and.Y(1)<0.d0) jmin = gw+1
 endif
 
-f_out = sum(f(imin:imax,jmin:jmax,kmin:kmax)*f(imin:imax,jmin:jmax,kmin:kmax))
-
-f_out = f_out
-
 Nout = (imax-imin+1)*(jmax-jmin+1)*(kmax-kmin+1)
+  allocate(f_flat(Nout))
+  f_flat = reshape(f(imin:imax,jmin:jmax,kmin:kmax), [Nout])
+  f_out = DDOT(Nout, f_flat, 1, f_flat, 1)
+  deallocate(f_flat)
 
   return
 
@@ -1584,8 +1870,11 @@ Nout = (imax-imin+1)*(jmax-jmin+1)*(kmax-kmin+1)
 ! 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
 
@@ -1680,6 +1969,7 @@ Nout = (imax-imin+1)*(jmax-jmin+1)*(kmax-kmin+1)
   real*8, dimension(ORDN,ORDN) :: tmp2
   real*8, dimension(ORDN) :: tmp1
   real*8, dimension(3) :: SoAh
+  real*8, external :: DDOT
 
 ! +1 because c++ gives 0 for first point
   cxB = inds+1  
@@ -1725,10 +2015,7 @@ Nout = (imax-imin+1)*(jmax-jmin+1)*(kmax-kmin+1)
     tmp1 = tmp1 + coef(ORDN+m)*tmp2(:,m)
   enddo
 
-  f_int=0
+  f_int = DDOT(ORDN, coef(1:ORDN), 1, tmp1, 1)
-  do m=1,ORDN
-    f_int = f_int + coef(m)*tmp1(m)
-  enddo
 
   return
 
@@ -1758,6 +2045,7 @@ Nout = (imax-imin+1)*(jmax-jmin+1)*(kmax-kmin+1)
   real*8, dimension(ORDN,ORDN) :: ya
   real*8, dimension(ORDN) :: tmp1
   real*8, dimension(2) :: SoAh
+  real*8, external :: DDOT
 
 ! +1 because c++ gives 0 for first point
   cxB = inds(1:2)+1  
@@ -1792,10 +2080,7 @@ Nout = (imax-imin+1)*(jmax-jmin+1)*(kmax-kmin+1)
     tmp1 = tmp1 + coef(ORDN+m)*ya(:,m)
   enddo
 
-  f_int=0
+  f_int = DDOT(ORDN, coef(1:ORDN), 1, tmp1, 1)
-  do m=1,ORDN
-    f_int = f_int + coef(m)*tmp1(m)
-  enddo
 
   return
 
@@ -1826,6 +2111,7 @@ Nout = (imax-imin+1)*(jmax-jmin+1)*(kmax-kmin+1)
   real*8, dimension(ORDN) :: ya
   real*8 :: SoAh
   integer,dimension(3) :: inds
+  real*8, external :: DDOT
 
 ! +1 because c++ gives 0 for first point
   inds = indsi + 1
@@ -1886,10 +2172,7 @@ Nout = (imax-imin+1)*(jmax-jmin+1)*(kmax-kmin+1)
           write(*,*)"error in global_interpind1d, not recognized dumyd = ",dumyd
   endif
 
-  f_int=0
+  f_int = DDOT(ORDN, coef, 1, ya, 1)
-  do m=1,ORDN
-    f_int = f_int + coef(m)*ya(m)
-  enddo
 
   return
 
@@ -2128,17 +2411,25 @@ Nout = (imax-imin+1)*(jmax-jmin+1)*(kmax-kmin+1)
   real*8 :: gont
 
   integer :: i
+  real*8, parameter, dimension(0:20) :: fact_table = [ &
+    1.d0, 1.d0, 2.d0, 6.d0, 24.d0, 120.d0, 720.d0, 5040.d0, 40320.d0, &
+    362880.d0, 3628800.d0, 39916800.d0, 479001600.d0, 6227020800.d0, &
+    87178291200.d0, 1307674368000.d0, 20922789888000.d0, &
+    355687428096000.d0, 6402373705728000.d0, 121645100408832000.d0, &
+    2432902008176640000.d0 ]
 
 ! sanity check
   if(N < 0)then
      write(*,*) "ffact: error input for factorial"
+     gont = 1.d0
      return
   endif
 
-  gont = 1.d0
+  if(N <= 20)then
-  do i=1,N
+     gont = fact_table(N)
-     gont = gont*i
+  else
-  enddo
+     gont = exp(log_gamma(dble(N+1)))
+  endif
 
   return
 

AMSS_NCKU_source/fmisc.h

@@ -13,6 +13,7 @@
 #define f_global_interpind2d global_interpind2d
 #define f_global_interpind1d global_interpind1d
 #define f_l2normhelper l2normhelper
+#define f_l2normhelper7 l2normhelper7
 #define f_l2normhelper_sh l2normhelper_sh
 #define f_l2normhelper_sh_rms l2normhelper_sh_rms
 #define f_average average
@@ -42,6 +43,7 @@
 #define f_global_interpind2d GLOBAL_INTERPIND2D
 #define f_global_interpind1d GLOBAL_INTERPIND1D
 #define f_l2normhelper L2NORMHELPER
+#define f_l2normhelper7 L2NORMHELPER7
 #define f_l2normhelper_sh L2NORMHELPER_SH
 #define f_l2normhelper_sh_rms L2NORMHELPER_SH_RMS
 #define f_average AVERAGE
@@ -71,6 +73,7 @@
 #define f_global_interpind2d global_interpind2d_
 #define f_global_interpind1d global_interpind1d_
 #define f_l2normhelper l2normhelper_
+#define f_l2normhelper7 l2normhelper7_
 #define f_l2normhelper_sh l2normhelper_sh_
 #define f_l2normhelper_sh_rms l2normhelper_sh_rms_
 #define f_average average_
@@ -164,6 +167,15 @@ extern "C"
 						double *, double &, int &);
 }
 
+extern "C"
+{
+	void f_l2normhelper7(int *, double *, double *, double *,
+						 double &, double &, double &,
+						 double &, double &, double &,
+						 double *, double *, double *, double *,
+						 double *, double *, double *, double *, int &);
+}
+
 extern "C"
 {
 	void f_l2normhelper_sh(int *, double *, double *, double *,

AMSS_NCKU_source/macrodef.h

@@ -2,7 +2,7 @@
 #ifndef MICRODEF_H
 #define MICRODEF_H
 
-#include "microdef.fh"
+#include "macrodef.fh"
 
 // application parameters
 

AMSS_NCKU_source/makefile

@@ -2,6 +2,20 @@
 
 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)
+
+ARCH_OPT = -march=x86-64-v4
+CXXAPPFLAGS = -O3 $(ARCH_OPT) -fp-model fast=2 -fma -ipo \
+              -Dfortran3 -Dnewc -I${MKLROOT}/include
+f90appflags = -O3 $(ARCH_OPT) -fp-model fast=2 -fma -ipo \
+              -align array64byte -fpp -I${MKLROOT}/include $(POLINT6_FLAG)
+TP_OPTFLAGS = -O3 $(ARCH_OPT) -fp-model fast=2 -fma -ipo \
+              -Dfortran3 -Dnewc -I${MKLROOT}/include
+
 .SUFFIXES: .o .f90 .C .for .cu
 
 .f90.o:
@@ -16,6 +30,12 @@ include makefile.inc
 .cu.o:
 	$(Cu) $(CUDA_APP_FLAGS) -c $< -o $@ $(CUDA_LIB_PATH)
 
+TwoPunctures.o: TwoPunctures.C
+	${CXX} $(TP_OPTFLAGS) -qopenmp -c $< -o $@
+
+TwoPunctureABE.o: TwoPunctureABE.C
+	${CXX} $(TP_OPTFLAGS) -qopenmp -c $< -o $@
+
 # Input files
 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\
@@ -96,7 +116,7 @@ ABEGPU: $(C++FILES_GPU) $(F90FILES) $(F77FILES) $(AHFDOBJS) $(CUDAFILES)
 	$(CLINKER) $(CXXAPPFLAGS) -o $@ $(C++FILES_GPU) $(F90FILES) $(F77FILES) $(AHFDOBJS) $(CUDAFILES) $(LDLIBS)
 
 TwoPunctureABE: $(TwoPunctureFILES)
-	$(CLINKER) $(CXXAPPFLAGS) -o $@ $(TwoPunctureFILES) $(LDLIBS)
+	$(CLINKER) $(TP_OPTFLAGS) -qopenmp -o $@ $(TwoPunctureFILES) $(LDLIBS)
 
 clean:
 	rm *.o ABE ABEGPU TwoPunctureABE make.log -f

AMSS_NCKU_source/makefile.inc

@@ -1,19 +1,32 @@
 
+## GCC version (commented out)
 ## filein  = -I/usr/include -I/usr/lib/x86_64-linux-gnu/mpich/include -I/usr/lib/x86_64-linux-gnu/openmpi/lib/ -I/usr/lib/gcc/x86_64-linux-gnu/11/ -I/usr/include/c++/11/
+## filein  = -I/usr/include/ -I/usr/include/openmpi-x86_64/ -I/usr/lib/x86_64-linux-gnu/openmpi/include/ -I/usr/lib/x86_64-linux-gnu/openmpi/lib/ -I/usr/lib/gcc/x86_64-linux-gnu/11/ -I/usr/include/c++/11/
+## LDLIBS  = -L/usr/lib/x86_64-linux-gnu -L/usr/lib64 -L/usr/lib/gcc/x86_64-linux-gnu/11 -lgfortran -lmpi -lgfortran
 
-filein  = -I/usr/include/ -I/usr/include/openmpi-x86_64/ -I/usr/lib/x86_64-linux-gnu/openmpi/include/ -I/usr/lib/x86_64-linux-gnu/openmpi/lib/ -I/usr/lib/gcc/x86_64-linux-gnu/11/ -I/usr/include/c++/11/
+## Intel oneAPI version with oneMKL
+filein  = -I/usr/include/ -I${MKLROOT}/include
 
-## LDLIBS  = -L/usr/lib/x86_64-linux-gnu -lmpich -L/usr/lib64 -L/usr/lib/gcc/x86_64-linux-gnu/11 -lgfortran
+## Use sequential oneMKL to avoid introducing extra OpenMP behavior into ABE.
-LDLIBS  = -L/usr/lib/x86_64-linux-gnu -L/usr/lib64 -L/usr/lib/gcc/x86_64-linux-gnu/11 -lgfortran -lmpi -lgfortran
+LDLIBS  = -L${MKLROOT}/lib -lmkl_intel_lp64 -lmkl_sequential -lmkl_core -lifcore -limf -lpthread -lm -ldl -liomp5
 
-CXXAPPFLAGS  = -O0 -Wno-deprecated -Dfortran3 -Dnewc
+## Optional Intel oneTBB allocator, kept aligned with main's build environment.
-#f90appflags = -O0 -fpp
+USE_TBBMALLOC ?= 1
-f90appflags  = -O0 -x f95-cpp-input
+TBBMALLOC_SO ?= /home/intel/oneapi/2025.3/lib/libtbbmalloc.so
-f90          = gfortran 
+ifneq ($(wildcard $(TBBMALLOC_SO)),)
-f77          = gfortran 
+TBBMALLOC_LIBS = -Wl,--no-as-needed $(TBBMALLOC_SO) -Wl,--as-needed
-CXX          = g++
+else
-CC           = gcc
+TBBMALLOC_LIBS = -Wl,--no-as-needed -ltbbmalloc -Wl,--as-needed
-CLINKER      = mpic++ 
+endif
+ifeq ($(USE_TBBMALLOC),1)
+LDLIBS := $(TBBMALLOC_LIBS) $(LDLIBS)
+endif
+
+f90          = ifx
+f77          = ifx
+CXX          = icpx
+CC           = icx
+CLINKER      = mpiicpx
 
 Cu = nvcc
 CUDA_LIB_PATH = -L/usr/lib/cuda/lib64 -I/usr/include -I/usr/lib/cuda/include

parallel_plot_helper.py

@@ -0,0 +1,12 @@
+import multiprocessing
+
+
+def run_plot_task(task):
+    func, args = task
+    return func(*args)
+
+
+def run_plot_tasks_parallel(plot_tasks):
+    ctx = multiprocessing.get_context('fork')
+    with ctx.Pool() as pool:
+        pool.map(run_plot_task, plot_tasks)

plot_GW_strain_amplitude_xiaoqu.py

@@ -11,6 +11,8 @@
 import numpy                               ## numpy for array operations
 import scipy                               ## scipy for interpolation and signal processing
 import math
+import matplotlib
+matplotlib.use('Agg')                      ## use non-interactive backend for multiprocessing safety
 import matplotlib.pyplot    as     plt     ## matplotlib for plotting
 import os                                  ## os for system/file operations
 

plot_binary_data.py

@@ -8,16 +8,21 @@
 ##
 #################################################
 
+## Restrict OpenMP to one thread per process so that parallel
+## subprocess plotting does not multiply BLAS thread counts.
+import os
+os.environ.setdefault("OMP_NUM_THREADS", "1")
+
 import numpy
 import scipy
+import matplotlib
+matplotlib.use('Agg')                      ## use non-interactive backend for multiprocessing safety
 import matplotlib.pyplot    as     plt
 from   matplotlib.colors    import LogNorm
 from   mpl_toolkits.mplot3d import Axes3D
 ## import torch
 import AMSS_NCKU_Input      as input_data
 
-import os
-
 
 #########################################################################################
 
@@ -192,3 +197,11 @@ def get_data_xy( Rmin, Rmax, n, data0, time, figure_title, figure_outdir ):
 
 ####################################################################################
 
+## Allow standalone subprocess execution for parallel binary-data plotting.
+if __name__ == '__main__':
+    import sys
+    if len(sys.argv) != 4:
+        print(f"Usage: {sys.argv[0]} <filename> <binary_outdir> <figure_outdir>")
+        sys.exit(1)
+    plot_binary_data(sys.argv[1], sys.argv[2], sys.argv[3])
+

plot_xiaoqu.py

@@ -8,6 +8,8 @@
 #################################################
 
 import numpy                               ## numpy for array operations
+import matplotlib
+matplotlib.use('Agg')                      ## use non-interactive backend for multiprocessing safety
 import matplotlib.pyplot    as     plt     ## matplotlib for plotting
 from   mpl_toolkits.mplot3d import Axes3D  ## needed for 3D plots
 import glob
@@ -15,6 +17,9 @@ import os                                  ## operating system utilities
 
 import plot_binary_data
 import AMSS_NCKU_Input as input_data
+import subprocess
+import sys
+import multiprocessing
 
 # plt.rcParams['text.usetex'] = True  ## enable LaTeX fonts in plots
 
@@ -50,10 +55,34 @@ def generate_binary_data_plot( binary_outdir, figure_outdir ):
         file_list.append(x)
         print(x)
 
-    ## Plot each file in the list
+    ## Plot each file in parallel using subprocesses.
+    ## Each subprocess starts with BLAS thread limits in plot_binary_data.py.
+    script = os.path.join( os.path.dirname(__file__), "plot_binary_data.py" )
+    max_workers = min( multiprocessing.cpu_count(), len(file_list) ) if file_list else 0
+
+    running = []
+    failed  = []
     for filename in file_list:
         print(filename)
-        plot_binary_data.plot_binary_data(filename, binary_outdir, figure_outdir)
+        proc = subprocess.Popen(
+            [sys.executable, script, filename, binary_outdir, figure_outdir],
+        )
+        running.append( (proc, filename) )
+        if len(running) >= max_workers:
+            p, fn = running.pop(0)
+            p.wait()
+            if p.returncode != 0:
+                failed.append(fn)
+
+    for p, fn in running:
+        p.wait()
+        if p.returncode != 0:
+            failed.append(fn)
+
+    if failed:
+        print( " WARNING: the following binary data plots failed:" )
+        for fn in failed:
+            print( "   ", fn )
 
     print(                        )
     print( " Binary Data Plot Has been Finished " )