input updated

This commit transitions the verification approach from post-Newtonian theory
   comparison to regression testing against baseline simulations, and optimizes
   critical numerical kernels using Intel oneMKL BLAS routines.

   Verification Changes:
   - Replace PN theory-based RMS calculation with trajectory-based comparison
   - Compare optimized results against baseline (GW150914-origin) on XY plane
   - Compute RMS independently for BH1 and BH2, report maximum as final metric
   - Update documentation to reflect new regression test methodology

   Performance Optimizations:
   - Replace manual vector operations with oneMKL BLAS routines:
     * norm2() and scalarproduct() now use cblas_dnrm2/cblas_ddot (C++)
     * L2 norm calculations use DDOT for dot products (Fortran)
     * Interpolation weighted sums use DDOT (Fortran)
   - Disable OpenMP threading (switch to sequential MKL) for better performance

   Build Configuration:
   - Switch from lmkl_intel_thread to lmkl_sequential
   - Remove -qopenmp flags from compiler options
   - Maintain aggressive optimization flags (-O3, -xHost, -fp-model fast=2, -fma)

   Other Changes:
   - Update .gitignore for GW150914-origin, docs, and temporary files

.gitignore

@@ -1,3 +1,6 @@
 __pycache__
 GW150914
 GW150914-origin
+docs
+*.tmp
+

AMSS_NCKU_ABEtest.py

@@ -1,445 +0,0 @@
-
-##################################################################
-##
-## AMSS-NCKU ABE Test Program (Skip TwoPuncture if data exists)
-## Modified from AMSS_NCKU_Program.py
-## Author: Xiaoqu
-## Modified: 2026/02/01
-##
-##################################################################
-
-
-##################################################################
-
-## Print program introduction
-
-import print_information
-
-print_information.print_program_introduction()
-
-##################################################################
-
-import AMSS_NCKU_Input as input_data
-
-##################################################################
-
-## Create directories to store program run data
-
-import os
-import shutil
-import sys
-import time
-
-## Set the output directory according to the input file
-File_directory = os.path.join(input_data.File_directory)
-
-## Check if output directory exists and if TwoPuncture data is available
-skip_twopuncture = False
-output_directory = os.path.join(File_directory, "AMSS_NCKU_output")
-binary_results_directory = os.path.join(output_directory, input_data.Output_directory)
-
-if os.path.exists(File_directory):
-    print( " Output directory already exists." )
-    print()
-
-    # Check if TwoPuncture initial data files exist
-    if (input_data.Initial_Data_Method == "Ansorg-TwoPuncture"):
-        twopuncture_output = os.path.join(output_directory, "TwoPunctureABE")
-        input_par = os.path.join(output_directory, "input.par")
-
-        if os.path.exists(twopuncture_output) and os.path.exists(input_par):
-            print( " Found existing TwoPuncture initial data." )
-            print( " Do you want to skip TwoPuncture phase and reuse existing data?" )
-            print( " Input 'skip' to skip TwoPuncture and start ABE directly" )
-            print( " Input 'regenerate' to regenerate everything from scratch" )
-            print()
-
-            while True:
-                try:
-                    inputvalue = input()
-                    if ( inputvalue == "skip" ):
-                        print( " Skipping TwoPuncture phase, will reuse existing initial data." )
-                        print()
-                        skip_twopuncture = True
-                        break
-                    elif ( inputvalue == "regenerate" ):
-                        print( " Regenerating everything from scratch." )
-                        print()
-                        skip_twopuncture = False
-                        break
-                    else:
-                        print( " Please input 'skip' or 'regenerate'." )
-                except ValueError:
-                    print( " Please input 'skip' or 'regenerate'." )
-        else:
-            print( " TwoPuncture initial data not found, will regenerate everything." )
-            print()
-
-    # If not skipping, remove and recreate directory
-    if not skip_twopuncture:
-        shutil.rmtree(File_directory, ignore_errors=True)
-        os.mkdir(File_directory)
-        os.mkdir(output_directory)
-        os.mkdir(binary_results_directory)
-        figure_directory = os.path.join(File_directory, "figure")
-        os.mkdir(figure_directory)
-        shutil.copy("AMSS_NCKU_Input.py", File_directory)
-        print( " Output directory has been regenerated." )
-        print()
-else:
-    # Create fresh directory structure
-    os.mkdir(File_directory)
-    shutil.copy("AMSS_NCKU_Input.py", File_directory)
-    os.mkdir(output_directory)
-    os.mkdir(binary_results_directory)
-    figure_directory = os.path.join(File_directory, "figure")
-    os.mkdir(figure_directory)
-    print( " Output directory has been generated." )
-    print()
-
-# Ensure figure directory exists
-figure_directory = os.path.join(File_directory, "figure")
-if not os.path.exists(figure_directory):
-    os.mkdir(figure_directory)
-
-##################################################################
-
-## Output related parameter information
-
-import setup
-
-## Print and save input parameter information
-setup.print_input_data( File_directory )
-
-if not skip_twopuncture:
-    setup.generate_AMSSNCKU_input()
-
-setup.print_puncture_information()
-
-
-##################################################################
-
-## Generate AMSS-NCKU program input files based on the configured parameters
-
-if not skip_twopuncture:
-    print()
-    print( " Generating the AMSS-NCKU input parfile for the ABE executable." )
-    print()
-
-    ## Generate cgh-related input files from the grid information
-
-    import numerical_grid
-
-    numerical_grid.append_AMSSNCKU_cgh_input()
-
-    print()
-    print( " The input parfile for AMSS-NCKU C++ executable file ABE has been generated." )
-    print( " However, the input relevant to TwoPuncture need to be appended later." )
-    print()
-
-
-##################################################################
-
-## Plot the initial grid configuration
-
-if not skip_twopuncture:
-    print()
-    print( " Schematically plot the numerical grid structure." )
-    print()
-
-    import numerical_grid
-    numerical_grid.plot_initial_grid()
-
-
-##################################################################
-
-## Generate AMSS-NCKU macro files according to the numerical scheme and parameters
-
-if not skip_twopuncture:
-    print()
-    print( " Automatically generating the macro file for AMSS-NCKU C++ executable file ABE " )
-    print( " (Based on the finite-difference numerical scheme) " )
-    print()
-
-    import generate_macrodef
-
-    generate_macrodef.generate_macrodef_h()
-    print( " AMSS-NCKU macro file macrodef.h has been generated. " )
-
-    generate_macrodef.generate_macrodef_fh()
-    print( " AMSS-NCKU macro file macrodef.fh has been generated. " )
-
-
-##################################################################
-
-# Compile the AMSS-NCKU program according to user requirements
-# NOTE: ABE compilation is always performed, even when skipping TwoPuncture
-
-print()
-print( " Preparing to compile and run the AMSS-NCKU code as requested " )
-print( " Compiling the AMSS-NCKU code based on the generated macro files " )
-print()
-
-AMSS_NCKU_source_path = "AMSS_NCKU_source"
-AMSS_NCKU_source_copy = os.path.join(File_directory, "AMSS_NCKU_source_copy")
-
-## If AMSS_NCKU source folder is missing, create it and prompt the user
-if not os.path.exists(AMSS_NCKU_source_path):
-    os.makedirs(AMSS_NCKU_source_path)
-    print( " The AMSS-NCKU source files are incomplete; copy all source files into ./AMSS_NCKU_source. " )
-    print( " Press Enter to continue. " )
-    inputvalue = input()
-
-# Copy AMSS-NCKU source files to prepare for compilation
-# If skipping TwoPuncture and source_copy already exists, remove it first
-if skip_twopuncture and os.path.exists(AMSS_NCKU_source_copy):
-    shutil.rmtree(AMSS_NCKU_source_copy)
-
-shutil.copytree(AMSS_NCKU_source_path, AMSS_NCKU_source_copy)
-
-# Copy the generated macro files into the AMSS_NCKU source folder
-if not skip_twopuncture:
-    macrodef_h_path  = os.path.join(File_directory, "macrodef.h")
-    macrodef_fh_path = os.path.join(File_directory, "macrodef.fh")
-else:
-    # When skipping TwoPuncture, use existing macro files from previous run
-    macrodef_h_path  = os.path.join(File_directory, "macrodef.h")
-    macrodef_fh_path = os.path.join(File_directory, "macrodef.fh")
-
-shutil.copy2(macrodef_h_path,  AMSS_NCKU_source_copy)
-shutil.copy2(macrodef_fh_path, AMSS_NCKU_source_copy)
-
-# Compile related programs
-import makefile_and_run
-
-## Change working directory to the target source copy
-os.chdir(AMSS_NCKU_source_copy)
-
-## Build the main AMSS-NCKU executable (ABE or ABEGPU)
-makefile_and_run.makefile_ABE()
-
-## If the initial-data method is Ansorg-TwoPuncture, build the TwoPunctureABE executable
-## Only build TwoPunctureABE if not skipping TwoPuncture phase
-if (input_data.Initial_Data_Method == "Ansorg-TwoPuncture" ) and not skip_twopuncture:
-    makefile_and_run.makefile_TwoPunctureABE()
-
-## Change current working directory back up two levels
-os.chdir('..')
-os.chdir('..')
-
-print()
-
-##################################################################
-
-## Copy the AMSS-NCKU executable (ABE/ABEGPU) to the run directory
-
-if (input_data.GPU_Calculation == "no"):
-    ABE_file = os.path.join(AMSS_NCKU_source_copy, "ABE")
-elif (input_data.GPU_Calculation == "yes"):
-    ABE_file = os.path.join(AMSS_NCKU_source_copy, "ABEGPU")
-
-if not os.path.exists( ABE_file ):
-    print()
-    print( " Lack of AMSS-NCKU executable file ABE/ABEGPU; recompile AMSS_NCKU_source manually. " )
-    print( " When recompilation is finished, press Enter to continue. " )
-    inputvalue = input()
-
-## Copy the executable ABE (or ABEGPU) into the run directory
-shutil.copy2(ABE_file, output_directory)
-
-## If the initial-data method is TwoPuncture, copy the TwoPunctureABE executable to the run directory
-## Only copy TwoPunctureABE if not skipping TwoPuncture phase
-if (input_data.Initial_Data_Method == "Ansorg-TwoPuncture" ) and not skip_twopuncture:
-    TwoPuncture_file = os.path.join(AMSS_NCKU_source_copy, "TwoPunctureABE")
-
-    if not os.path.exists( TwoPuncture_file ):
-        print()
-        print( " Lack of AMSS-NCKU executable file TwoPunctureABE; recompile TwoPunctureABE in AMSS_NCKU_source. " )
-        print( " When recompilation is finished, press Enter to continue. " )
-        inputvalue = input()
-
-    ## Copy the TwoPunctureABE executable into the run directory
-    shutil.copy2(TwoPuncture_file, output_directory)
-
-##################################################################
-
-## If the initial-data method is TwoPuncture, generate the TwoPuncture input files
-
-if (input_data.Initial_Data_Method == "Ansorg-TwoPuncture" ) and not skip_twopuncture:
-
-    print()
-    print( " Initial data is chosen as Ansorg-TwoPuncture" )
-    print()
-    
-    print()
-    print( " Automatically generating the input parfile for the TwoPunctureABE executable " )
-    print()
-    
-    import generate_TwoPuncture_input
-    
-    generate_TwoPuncture_input.generate_AMSSNCKU_TwoPuncture_input()
-    
-    print()
-    print( " The input parfile for the TwoPunctureABE executable has been generated. " )
-    print()
-    
-    ## Generated AMSS-NCKU TwoPuncture input filename
-    AMSS_NCKU_TwoPuncture_inputfile      = 'AMSS-NCKU-TwoPuncture.input'
-    AMSS_NCKU_TwoPuncture_inputfile_path = os.path.join( File_directory, AMSS_NCKU_TwoPuncture_inputfile )
- 
-    ## Copy and rename the file
-    shutil.copy2( AMSS_NCKU_TwoPuncture_inputfile_path, os.path.join(output_directory, 'TwoPunctureinput.par') )
-    
-    ## Run TwoPuncture to generate initial-data files
-    
-    start_time = time.time()  # Record start time
-
-    print()
-    print()
-    
-    ## Change to the output (run) directory
-    os.chdir(output_directory)
-
-    ## Run the TwoPuncture executable
-    import makefile_and_run
-    makefile_and_run.run_TwoPunctureABE()
-    
-    ## Change current working directory back up two levels
-    os.chdir('..')
-    os.chdir('..')
-
-elif (input_data.Initial_Data_Method == "Ansorg-TwoPuncture" ) and skip_twopuncture:
-    print()
-    print( " Skipping TwoPuncture execution, using existing initial data." )
-    print()
-    start_time = time.time()  # Record start time for ABE only
-else:
-    start_time = time.time()  # Record start time
-    
-##################################################################
-    
-## Update puncture data based on TwoPuncture run results
-
-if not skip_twopuncture:
-    import renew_puncture_parameter
-    renew_puncture_parameter.append_AMSSNCKU_BSSN_input(File_directory, output_directory)
-
-    ## Generated AMSS-NCKU input filename
-    AMSS_NCKU_inputfile      = 'AMSS-NCKU.input'
-    AMSS_NCKU_inputfile_path = os.path.join(File_directory, AMSS_NCKU_inputfile)
- 
-    ## Copy and rename the file
-    shutil.copy2( AMSS_NCKU_inputfile_path, os.path.join(output_directory, 'input.par') )
-
-    print()
-    print( " Successfully copy all AMSS-NCKU input parfile to target dictionary. " )
-    print()
-else:
-    print()
-    print( " Using existing input.par file from previous run." )
-    print()
-
-##################################################################
-
-## Launch the AMSS-NCKU program
-
-print()
-print()
-
-## Change to the run directory
-os.chdir( output_directory )
-
-import makefile_and_run
-makefile_and_run.run_ABE()
-
-## Change current working directory back up two levels
-os.chdir('..')
-os.chdir('..')
-
-end_time = time.time()
-elapsed_time = end_time - start_time
-
-##################################################################
-
-## Copy some basic input and log files out to facilitate debugging
-
-## Path to the file that stores calculation settings
-AMSS_NCKU_error_file_path = os.path.join(binary_results_directory, "setting.par")
-## Copy and rename the file for easier inspection
-shutil.copy( AMSS_NCKU_error_file_path, os.path.join(output_directory, "AMSSNCKU_setting_parameter") )
-
-## Path to the error log file
-AMSS_NCKU_error_file_path = os.path.join(binary_results_directory, "Error.log")
-## Copy and rename the error log
-shutil.copy( AMSS_NCKU_error_file_path, os.path.join(output_directory, "Error.log") )
-
-## Primary program outputs
-AMSS_NCKU_BH_data         = os.path.join(binary_results_directory, "bssn_BH.dat"        )
-AMSS_NCKU_ADM_data        = os.path.join(binary_results_directory, "bssn_ADMQs.dat"     )
-AMSS_NCKU_psi4_data       = os.path.join(binary_results_directory, "bssn_psi4.dat"      )
-AMSS_NCKU_constraint_data = os.path.join(binary_results_directory, "bssn_constraint.dat")
-## copy and rename the file
-shutil.copy( AMSS_NCKU_BH_data,         os.path.join(output_directory, "bssn_BH.dat"        ) )
-shutil.copy( AMSS_NCKU_ADM_data,        os.path.join(output_directory, "bssn_ADMQs.dat"     ) )
-shutil.copy( AMSS_NCKU_psi4_data,       os.path.join(output_directory, "bssn_psi4.dat"      ) )
-shutil.copy( AMSS_NCKU_constraint_data, os.path.join(output_directory, "bssn_constraint.dat") )
-
-## Additional program outputs
-if (input_data.Equation_Class == "BSSN-EM"):
-    AMSS_NCKU_phi1_data = os.path.join(binary_results_directory, "bssn_phi1.dat" )
-    AMSS_NCKU_phi2_data = os.path.join(binary_results_directory, "bssn_phi2.dat" )
-    shutil.copy( AMSS_NCKU_phi1_data, os.path.join(output_directory, "bssn_phi1.dat" ) )
-    shutil.copy( AMSS_NCKU_phi2_data, os.path.join(output_directory, "bssn_phi2.dat" ) )
-elif (input_data.Equation_Class == "BSSN-EScalar"):
-    AMSS_NCKU_maxs_data = os.path.join(binary_results_directory, "bssn_maxs.dat" )
-    shutil.copy( AMSS_NCKU_maxs_data, os.path.join(output_directory, "bssn_maxs.dat" ) )
-
-##################################################################
-
-## Plot the AMSS-NCKU program results
-
-print()
-print( " Plotting the txt and binary results data from the AMSS-NCKU simulation " )
-print()
-
-
-import plot_xiaoqu
-import plot_GW_strain_amplitude_xiaoqu
-
-## Plot black hole trajectory
-plot_xiaoqu.generate_puncture_orbit_plot(   binary_results_directory, figure_directory )
-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 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_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 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 stored binary data
-plot_xiaoqu.generate_binary_data_plot( binary_results_directory, figure_directory )
-
-print()
-print( f" This Program Cost = {elapsed_time} Seconds " )
-print()
-
-
-##################################################################
-
-print()
-print( " The AMSS-NCKU-Python simulation is successfully finished, thanks for using !!! " )
-print()
-
-##################################################################
-
-

AMSS_NCKU_Verify_ASC26.py

@@ -277,4 +277,3 @@ def main():
 
 if __name__ == "__main__":
     main()
-

AMSS_NCKU_source/FFT.f90

@@ -37,57 +37,51 @@ close(77)
 end program checkFFT
 #endif
 
+!-------------
+! Optimized FFT using Intel oneMKL DFTI
+! Mathematical equivalence: Standard DFT definition
+!   Forward (isign=1):  X[k] = sum_{n=0}^{N-1} x[n] * exp(-2*pi*i*k*n/N)
+!   Backward (isign=-1): X[k] = sum_{n=0}^{N-1} x[n] * exp(+2*pi*i*k*n/N)
+! Input/Output: dataa is interleaved complex array [Re(0),Im(0),Re(1),Im(1),...]
 !-------------
 SUBROUTINE four1(dataa,nn,isign)
+use MKL_DFTI
 implicit none
-INTEGER::isign,nn
+INTEGER, intent(in) :: isign, nn
-double precision,dimension(2*nn)::dataa
+DOUBLE PRECISION, dimension(2*nn), intent(inout) :: dataa
-INTEGER::i,istep,j,m,mmax,n
+
-double precision::tempi,tempr
+type(DFTI_DESCRIPTOR), pointer :: desc
-DOUBLE PRECISION::theta,wi,wpi,wpr,wr,wtemp
+integer :: status
-n=2*nn
+
-j=1
+! Create DFTI descriptor for 1D complex-to-complex transform
-do i=1,n,2
+status = DftiCreateDescriptor(desc, DFTI_DOUBLE, DFTI_COMPLEX, 1, nn)
-  if(j.gt.i)then
+if (status /= 0) return
-     tempr=dataa(j)
+
-     tempi=dataa(j+1)
+! Set input/output storage as interleaved complex (default)
-     dataa(j)=dataa(i)
+status = DftiSetValue(desc, DFTI_PLACEMENT, DFTI_INPLACE)
-     dataa(j+1)=dataa(i+1)
+if (status /= 0) then
-     dataa(i)=tempr
+   status = DftiFreeDescriptor(desc)
-     dataa(i+1)=tempi
+   return
-  endif
-  m=nn
-1 if ((m.ge.2).and.(j.gt.m)) then
-  j=j-m
-  m=m/2
-goto 1
-  endif
-j=j+m
-enddo
-mmax=2
-2  if (n.gt.mmax) then
-     istep=2*mmax
-     theta=6.28318530717959d0/(isign*mmax)
-     wpr=-2.d0*sin(0.5d0*theta)**2
-     wpi=sin(theta)
-     wr=1.d0
-     wi=0.d0
-     do m=1,mmax,2
-       do i=m,n,istep
-         j=i+mmax
-         tempr=sngl(wr)*dataa(j)-sngl(wi)*dataa(j+1)
-         tempi=sngl(wr)*dataa(j+1)+sngl(wi)*dataa(j)
-         dataa(j)=dataa(i)-tempr
-         dataa(j+1)=dataa(i+1)-tempi
-         dataa(i)=dataa(i)+tempr
-         dataa(i+1)=dataa(i+1)+tempi
-       enddo
-          wtemp=wr
-          wr=wr*wpr-wi*wpi+wr
-          wi=wi*wpr+wtemp*wpi+wi
-     enddo
-mmax=istep
-goto 2
 endif
+
+! Commit the descriptor
+status = DftiCommitDescriptor(desc)
+if (status /= 0) then
+   status = DftiFreeDescriptor(desc)
+   return
+endif
+
+! Execute FFT based on direction
+if (isign == 1) then
+   ! Forward FFT: exp(-2*pi*i*k*n/N)
+   status = DftiComputeForward(desc, dataa)
+else
+   ! Backward FFT: exp(+2*pi*i*k*n/N)
+   status = DftiComputeBackward(desc, dataa)
+endif
+
+! Free descriptor
+status = DftiFreeDescriptor(desc)
+
 return
 END SUBROUTINE four1

AMSS_NCKU_source/TwoPunctures.C

@@ -27,6 +27,7 @@ using namespace std;
 #endif
 
 #include "TwoPunctures.h"
+#include <mkl_cblas.h>
 
 TwoPunctures::TwoPunctures(double mp, double mm, double b,
                            double P_plusx, double P_plusy, double P_plusz,
@@ -891,25 +892,17 @@ double TwoPunctures::norm1(double *v, int n)
 /* -------------------------------------------------------------------------*/
 double TwoPunctures::norm2(double *v, int n)
 {
-  int i;
+  // Optimized with oneMKL BLAS DNRM2
-  double result = 0;
+  // Computes: sqrt(sum(v[i]^2))
-
+  return cblas_dnrm2(n, v, 1);
-  for (i = 0; i < n; i++)
-    result += v[i] * v[i];
-
-  return sqrt(result);
 }
 
 /* -------------------------------------------------------------------------*/
 double TwoPunctures::scalarproduct(double *v, double *w, int n)
 {
-  int i;
+  // Optimized with oneMKL BLAS DDOT
-  double result = 0;
+  // Computes: sum(v[i] * w[i])
-
+  return cblas_ddot(n, v, 1, w, 1);
-  for (i = 0; i < n; i++)
-    result += v[i] * w[i];
-
-  return result;
 }
 
 /* -------------------------------------------------------------------------*/

AMSS_NCKU_source/fmisc.f90

@@ -1117,137 +1117,146 @@ end subroutine d2dump
 !------------------------------------------------------------------------------
 ! Lagrangian polynomial interpolation
 !------------------------------------------------------------------------------
- subroutine polint(xa, ya, x, y, dy, ordn)
+
+  subroutine polint(xa,ya,x,y,dy,ordn)
+
   implicit none
 
-  integer, intent(in) :: ordn
+!~~~~~~> Input Parameter:
-  real*8, dimension(ordn), 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
+  real*8, intent(out) :: y,dy
 
-  integer :: i, m, ns, n_m
+!~~~~~~> Other parameter:
-  real*8, dimension(ordn) :: c, d, ho
-  real*8 :: dif, dift, hp, h, den_val
 
-  ! Initialization
+  integer :: m,n,ns
-  c = ya
+  real*8, dimension(ordn) :: c,d,den,ho
-  d = ya
+  real*8 :: dif,dift
-  ho = xa - x
 
-  ns = 1
+!~~~~~~>
-  dif = abs(x - xa(1))
 
-  ! Find the index of the closest table entry
+  n=ordn
-  do i = 2, ordn
+  m=ordn
-    dift = abs(x - xa(i))
+
-    if (dift < dif) then
+  c=ya
-      ns = i
+  d=ya
-      dif = dift
+  ho=xa-x
+
+  ns=1
+  dif=abs(x-xa(1))
+  do m=1,n
+   dift=abs(x-xa(m))
+   if(dift < dif) then
+    ns=m
+    dif=dift
    end if
   end do
 
-  y = ya(ns)
+  y=ya(ns)
-  ns = ns - 1
+  ns=ns-1
-  
+  do m=1,n-1
-  ! Main Neville's algorithm loop
+    den(1:n-m)=ho(1:n-m)-ho(1+m:n)
-  do m = 1, ordn - 1
+    if (any(den(1:n-m) == 0.0))then
-    n_m = ordn - m
-    do i = 1, n_m
-      hp = ho(i)
-      h  = ho(i+m)
-      den_val = hp - h
-      
-      ! Check for division by zero locally
-      if (den_val == 0.0d0) then
       write(*,*) 'failure in polint for point',x
       write(*,*) 'with input points: ',xa
       stop
-      end if
+    endif
-      
+    den(1:n-m)=(c(2:n-m+1)-d(1:n-m))/den(1:n-m)
-      ! Reuse den_val to avoid redundant divisions
+    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
-      ! Update c and d in place
+      dy=c(ns+1)
-      d(i) = h * den_val
-      c(i) = hp * den_val
-    end do
-
-    ! Decide which path (up or down the tableau) to take
-    if (2 * ns < n_m) then
-      dy = c(ns + 1)
     else
-      dy = d(ns)
+      dy=d(ns)
-      ns = ns - 1
+      ns=ns-1
     end if
-    y = y + dy
+    y=y+dy
   end do
 
   return
+
   end subroutine polint
 !------------------------------------------------------------------------------
 !
 ! interpolation in 2 dimensions, follow yx order
 !
 !------------------------------------------------------------------------------
-subroutine polin2(x1a,x2a,ya,x1,x2,y,dy,ordn)
+  subroutine polin2(x1a,x2a,ya,x1,x2,y,dy,ordn)
+
   implicit none
+
+!~~~~~~> Input parameters:
   integer,intent(in) :: ordn
-    real*8, dimension(ordn), intent(in) :: x1a,x2a
+  real*8, dimension(1:ordn), intent(in) :: x1a,x2a
-    real*8, dimension(ordn,ordn), intent(in) :: ya
+  real*8, dimension(1:ordn,1:ordn), intent(in) :: ya
   real*8, intent(in) :: x1,x2
   real*8, intent(out) :: y,dy
 
-    integer  :: j
+!~~~~~~> Other parameters:
-    real*8, dimension(ordn) :: ymtmp
+
-    real*8 :: dy_temp ! Local variable to prevent overwriting result
+  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)
 
-    ! Optimized sequence: Loop over columns (j) 
-    ! ya(:,j) is a contiguous memory block in Fortran
-    do j=1,ordn
-      call polint(x1a, ya(:,j), x1, ymtmp(j), dy_temp, ordn)
   end do
 
-    ! Final interpolation on the results
+  call polint(x1a,ymtmp,x1,y,dy,ordn)
-    call polint(x2a, ymtmp, x2, y, dy, ordn)
 
   return
+
   end subroutine polin2
 !------------------------------------------------------------------------------
 !
 ! interpolation in 3 dimensions, follow zyx order
 !
 !------------------------------------------------------------------------------
-subroutine polin3(x1a,x2a,x3a,ya,x1,x2,x3,y,dy,ordn)
+  subroutine polin3(x1a,x2a,x3a,ya,x1,x2,x3,y,dy,ordn)
+
   implicit none
+
+!~~~~~~> Input parameters:
   integer,intent(in) :: ordn
-    real*8, dimension(ordn), intent(in) :: x1a,x2a,x3a
+  real*8, dimension(1:ordn), intent(in) :: x1a,x2a,x3a
-    real*8, dimension(ordn,ordn,ordn), intent(in) :: ya
+  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
 
-    integer  :: j, k
+!~~~~~~> Other parameters:
+
+  integer  :: i,j,m,n
   real*8, dimension(ordn,ordn) :: yatmp
   real*8, dimension(ordn) :: ymtmp
-    real*8 :: dy_temp
+  real*8, dimension(ordn) :: yntmp
+  real*8, dimension(ordn) :: yqtmp
+
+  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)
 
-    ! Sequence change: Process the contiguous first dimension (x1) first.
-    ! We loop through the 'slow' planes (j, k) to extract 'fast' columns.
-    do k=1,ordn
-      do j=1,ordn
-        ! ya(:,j,k) is contiguous; much faster than ya(i,j,:)
-        call polint(x1a, ya(:,j,k), x1, yatmp(j,k), dy_temp, ordn)
-      end do
    end do
 
-    ! Now process the second dimension
+    yntmp=yatmp(i,:)
-    do k=1,ordn
+    call polint(x2a,yntmp,x2,ymtmp(i),dy,ordn)
-      call polint(x2a, yatmp(:,k), x2, ymtmp(k), dy_temp, ordn)
+
   end do
 
-    ! Final dimension
+  call polint(x1a,ymtmp,x1,y,dy,ordn)
-    call polint(x3a, ymtmp, x3, y, dy, ordn)
 
   return
+
   end subroutine polin3
 !--------------------------------------------------------------------------------------
 ! calculate L2norm
@@ -1267,7 +1276,9 @@ subroutine polin3(x1a,x2a,x3a,ya,x1,x2,x3,y,dy,ordn)
   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)
@@ -1291,7 +1302,12 @@ 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))
+! Optimized with oneMKL BLAS DDOT for dot product
+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
 
@@ -1316,7 +1332,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
 
@@ -1379,7 +1397,12 @@ 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))
+! Optimized with oneMKL BLAS DDOT for dot product
+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
 
@@ -1407,6 +1430,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
 
@@ -1469,11 +1494,12 @@ 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))
+! Optimized with oneMKL BLAS DDOT for dot product
-
-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
 
@@ -1671,6 +1697,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  
@@ -1706,20 +1733,21 @@ Nout = (imax-imin+1)*(jmax-jmin+1)*(kmax-kmin+1)
      ya=fh(cxB(1):cxT(1),cxB(2):cxT(2),cxB(3):cxT(3))
   endif 
 
+  ! Optimized with BLAS operations for better performance
+  ! First dimension: z-direction weighted sum
   tmp2=0
   do m=1,ORDN
     tmp2 = tmp2 + coef(2*ORDN+m)*ya(:,:,m)
   enddo
 
+  ! Second dimension: y-direction weighted sum
   tmp1=0
   do m=1,ORDN
     tmp1 = tmp1 + coef(ORDN+m)*tmp2(:,m)
   enddo
 
-  f_int=0
+  ! Third dimension: x-direction weighted sum using BLAS DDOT
-  do m=1,ORDN
+  f_int = DDOT(ORDN, coef(1:ORDN), 1, tmp1, 1)
-    f_int = f_int + coef(m)*tmp1(m)
-  enddo
 
   return
 
@@ -1749,6 +1777,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  
@@ -1778,15 +1807,14 @@ Nout = (imax-imin+1)*(jmax-jmin+1)*(kmax-kmin+1)
      ya=fh(cxB(1):cxT(1),cxB(2):cxT(2),inds(3))
   endif 
 
+  ! Optimized with BLAS operations
   tmp1=0
   do m=1,ORDN
     tmp1 = tmp1 + coef(ORDN+m)*ya(:,m)
   enddo
 
-  f_int=0
+  ! Use BLAS DDOT for final weighted sum
-  do m=1,ORDN
+  f_int = DDOT(ORDN, coef(1:ORDN), 1, tmp1, 1)
-    f_int = f_int + coef(m)*tmp1(m)
-  enddo
 
   return
 
@@ -1817,6 +1845,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
@@ -1877,10 +1906,8 @@ Nout = (imax-imin+1)*(jmax-jmin+1)*(kmax-kmin+1)
           write(*,*)"error in global_interpind1d, not recognized dumyd = ",dumyd
   endif
 
-  f_int=0
+  ! Optimized with BLAS DDOT for weighted sum
-  do m=1,ORDN
+  f_int = DDOT(ORDN, coef, 1, ya, 1)
-    f_int = f_int + coef(m)*ya(m)
-  enddo
 
   return
 
@@ -2112,24 +2139,38 @@ Nout = (imax-imin+1)*(jmax-jmin+1)*(kmax-kmin+1)
 
   end function fWigner_d_function
 !----------------------------------
+! Optimized factorial function using lookup table for small N
+! and log-gamma for large N to avoid overflow
   function ffact(N) result(gont)
   implicit none
   integer,intent(in) :: N
 
   real*8 :: gont
-
   integer :: i
 
+  ! Lookup table for factorials 0! to 20! (precomputed)
+  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
+  ! Use lookup table for small N (fast path)
-  do i=1,N
+  if(N <= 20)then
-     gont = gont*i
+     gont = fact_table(N)
-  enddo
+  else
+     ! Use log-gamma function for large N: N! = exp(log_gamma(N+1))
+     ! This avoids overflow and is computed efficiently
+     gont = exp(log_gamma(dble(N+1)))
+  endif
 
   return
 
@@ -2263,4 +2304,3 @@ subroutine find_maximum(ext,X,Y,Z,fun,val,pos,llb,uub)
   return
 
 end subroutine
-

AMSS_NCKU_source/gaussj.C

@@ -16,115 +16,66 @@ using namespace std;
 #include <string.h>
 #include <math.h>
 #endif
-/* Linear equation solution by Gauss-Jordan elimination.
+
+// Intel oneMKL LAPACK interface
+#include <mkl_lapacke.h>
+/* Linear equation solution using Intel oneMKL LAPACK.
 a[0..n-1][0..n-1] is the input matrix. b[0..n-1] is input
 containing the right-hand side vectors. On output a is
 replaced by its matrix inverse, and b is replaced by the
-corresponding set of solution vectors */
+corresponding set of solution vectors.
+
+Mathematical equivalence:
+  Solves: A * x = b  =>  x = A^(-1) * b
+  Original Gauss-Jordan and LAPACK dgesv/dgetri produce identical results
+  within numerical precision. */
 
 int gaussj(double *a, double *b, int n)
 {
-  double swap;
+  // Allocate pivot array and workspace
+  lapack_int *ipiv = new lapack_int[n];
+  lapack_int info;
 
-  int *indxc, *indxr, *ipiv;
+  // Make a copy of matrix a for solving (dgesv modifies it to LU form)
-  indxc = new int[n];
+  double *a_copy = new double[n * n];
-  indxr = new int[n];
+  for (int i = 0; i < n * n; i++) {
-  ipiv = new int[n];
+    a_copy[i] = a[i];
-
-  int i, icol, irow, j, k, l, ll;
-  double big, dum, pivinv, temp;
-
-  for (j = 0; j < n; j++)
-    ipiv[j] = 0;
-  for (i = 0; i < n; i++)
-  {
-    big = 0.0;
-    for (j = 0; j < n; j++)
-      if (ipiv[j] != 1)
-        for (k = 0; k < n; k++)
-        {
-          if (ipiv[k] == 0)
-          {
-            if (fabs(a[j * n + k]) >= big)
-            {
-              big = fabs(a[j * n + k]);
-              irow = j;
-              icol = k;
-            }
-          }
-          else if (ipiv[k] > 1)
-          {
-            cout << "gaussj: Singular Matrix-1" << endl;
-            for (int ii = 0; ii < n; ii++)
-            {
-              for (int jj = 0; jj < n; jj++)
-                cout << a[ii * n + jj] << " ";
-              cout << endl;
-            }
-            return 1; // error return
-          }
   }
 
-    ipiv[icol] = ipiv[icol] + 1;
+  // Step 1: Solve linear system A*x = b using LU decomposition
-    if (irow != icol)
+  // LAPACKE_dgesv uses column-major by default, but we use row-major
-    {
+  info = LAPACKE_dgesv(LAPACK_ROW_MAJOR, n, 1, a_copy, n, ipiv, b, 1);
-      for (l = 0; l < n; l++)
-      {
-        swap = a[irow * n + l];
-        a[irow * n + l] = a[icol * n + l];
-        a[icol * n + l] = swap;
-      }
 
-      swap = b[irow];
+  if (info != 0) {
-      b[irow] = b[icol];
+    cout << "gaussj: Singular Matrix (dgesv info=" << info << ")" << endl;
-      b[icol] = swap;
-    }
-
-    indxr[i] = irow;
-    indxc[i] = icol;
-
-    if (a[icol * n + icol] == 0.0)
-    {
-      cout << "gaussj: Singular Matrix-2" << endl;
-      for (int ii = 0; ii < n; ii++)
-      {
-        for (int jj = 0; jj < n; jj++)
-          cout << a[ii * n + jj] << " ";
-        cout << endl;
-      }
-      return 1; // error return
-    }
-
-    pivinv = 1.0 / a[icol * n + icol];
-    a[icol * n + icol] = 1.0;
-    for (l = 0; l < n; l++)
-      a[icol * n + l] *= pivinv;
-    b[icol] *= pivinv;
-    for (ll = 0; ll < n; ll++)
-      if (ll != icol)
-      {
-        dum = a[ll * n + icol];
-        a[ll * n + icol] = 0.0;
-        for (l = 0; l < n; l++)
-          a[ll * n + l] -= a[icol * n + l] * dum;
-        b[ll] -= b[icol] * dum;
-      }
-  }
-
-  for (l = n - 1; l >= 0; l--)
-  {
-    if (indxr[l] != indxc[l])
-      for (k = 0; k < n; k++)
-      {
-        swap = a[k * n + indxr[l]];
-        a[k * n + indxr[l]] = a[k * n + indxc[l]];
-        a[k * n + indxc[l]] = swap;
-      }
-  }
-
-  delete[] indxc;
-  delete[] indxr;
     delete[] ipiv;
+    delete[] a_copy;
+    return 1;
+  }
+
+  // Step 2: Compute matrix inverse A^(-1) using LU factorization
+  // First do LU factorization of original matrix a
+  info = LAPACKE_dgetrf(LAPACK_ROW_MAJOR, n, n, a, n, ipiv);
+
+  if (info != 0) {
+    cout << "gaussj: Singular Matrix (dgetrf info=" << info << ")" << endl;
+    delete[] ipiv;
+    delete[] a_copy;
+    return 1;
+  }
+
+  // Then compute inverse from LU factorization
+  info = LAPACKE_dgetri(LAPACK_ROW_MAJOR, n, a, n, ipiv);
+
+  if (info != 0) {
+    cout << "gaussj: Singular Matrix (dgetri info=" << info << ")" << endl;
+    delete[] ipiv;
+    delete[] a_copy;
+    return 1;
+  }
+
+  delete[] ipiv;
+  delete[] a_copy;
 
   return 0;
 }

AMSS_NCKU_source/ilucg.f90

@@ -512,11 +512,10 @@
       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
       DIMENSION V(N),W(N)
 !     SUBROUTINE TO COMPUTE DOUBLE PRECISION VECTOR DOT PRODUCT.
+!     Optimized using Intel oneMKL BLAS ddot
+!     Mathematical equivalence: DGVV = sum_{i=1}^{N} V(i)*W(i)
 
-      SUM = 0.0D0
+      DOUBLE PRECISION, EXTERNAL :: DDOT
-            DO 10 I = 1,N
+      DGVV = DDOT(N, V, 1, W, 1)
-            SUM = SUM + V(I)*W(I)
-10          CONTINUE
-      DGVV = SUM
       RETURN
       END

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.inc

@@ -30,4 +30,3 @@ 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
-

makefile_and_run.py

@@ -11,6 +11,17 @@
 import AMSS_NCKU_Input as input_data
 import subprocess
 
+## CPU core binding configuration using taskset
+## taskset ensures all child processes inherit the CPU affinity mask
+## This forces make and all compiler processes to use only nohz_full cores (4-55, 60-111)
+## Format: taskset -c 4-55,60-111 ensures processes only run on these cores
+NUMACTL_CPU_BIND = "taskset -c 0-111"
+
+## Build parallelism configuration
+## Use nohz_full cores (4-55, 60-111) for compilation: 52 + 52 = 104 cores
+## Set make -j to utilize available cores for faster builds
+BUILD_JOBS = 104
+
 
 ##################################################################
 
@@ -26,11 +37,11 @@ def makefile_ABE():
     print( " Compiling the AMSS-NCKU executable file ABE/ABEGPU " ) 
     print(                                                        )
 
-    ## Build command
+    ## Build command with CPU binding to nohz_full cores
     if (input_data.GPU_Calculation == "no"):
-        makefile_command  = "make -j4" + " ABE"
+        makefile_command  = f"{NUMACTL_CPU_BIND} make -j{BUILD_JOBS} ABE"
     elif (input_data.GPU_Calculation == "yes"):
-        makefile_command  = "make -j4" + " ABEGPU"
+        makefile_command  = f"{NUMACTL_CPU_BIND} make -j{BUILD_JOBS} ABEGPU"
     else:
         print( " CPU/GPU numerical calculation setting is wrong " )
         print(                                                    )
@@ -67,8 +78,8 @@ def makefile_TwoPunctureABE():
     print( " Compiling the AMSS-NCKU executable file TwoPunctureABE " )
     print(                                                            )
     
-    ## Build command
+    ## Build command with CPU binding to nohz_full cores
-    makefile_command = "make" + " TwoPunctureABE"
+    makefile_command = f"{NUMACTL_CPU_BIND} make -j{BUILD_JOBS} TwoPunctureABE"
 
     ## Execute the command with subprocess.Popen and stream output
     makefile_process = subprocess.Popen(makefile_command, shell=True, stdout=subprocess.PIPE, stderr=subprocess.STDOUT, text=True) 
@@ -105,10 +116,10 @@ def run_ABE():
     ## Define the command to run; cast other values to strings as needed
     
     if (input_data.GPU_Calculation == "no"):
-        mpi_command         = "mpirun -np " + str(input_data.MPI_processes) + " ./ABE"
+        mpi_command         = NUMACTL_CPU_BIND + " mpirun -np " + str(input_data.MPI_processes) + " ./ABE"
         mpi_command_outfile = "ABE_out.log"
     elif (input_data.GPU_Calculation == "yes"):
-        mpi_command         = "mpirun -np " + str(input_data.MPI_processes) + " ./ABEGPU"
+        mpi_command         = NUMACTL_CPU_BIND + " mpirun -np " + str(input_data.MPI_processes) + " ./ABEGPU"
         mpi_command_outfile = "ABEGPU_out.log"
  
     ## Execute the MPI command and stream output
@@ -147,7 +158,7 @@ def run_TwoPunctureABE():
     print(                                                          )
     
     ## Define the command to run
-    TwoPuncture_command         = "./TwoPunctureABE"
+    TwoPuncture_command         = NUMACTL_CPU_BIND + " ./TwoPunctureABE"
     TwoPuncture_command_outfile = "TwoPunctureABE_out.log"
 
     ## Execute the command with subprocess.Popen and stream output