baseline updated

2026-02-05 19:53:55 +08:00
17 changed files with 228 additions and 540 deletions
--- a/.gitignore
+++ b/.gitignore
@@ -1,6 +1,2 @@
 __pycache__
 GW150914
-GW150914-origin
-docs
-*.tmp
-
--- a/AMSS_NCKU_Input.py
+++ b/AMSS_NCKU_Input.py
@@ -16,7 +16,7 @@ import numpy
 File_directory   = "GW150914"                    ## output file directory
 Output_directory = "binary_output"               ## binary data file directory
                                                 ## The file directory name should not be too long
-MPI_processes    = 48                             ## number of mpi processes used in the simulation
+MPI_processes    = 64                             ## number of mpi processes used in the simulation

 GPU_Calculation  = "no"                          ## Use GPU or not 
                                                 ## (prefer "no" in the current version, because the GPU part may have bugs when integrated in this Python interface)
--- a/AMSS_NCKU_Verify_ASC26.py
+++ b/AMSS_NCKU_Verify_ASC26.py
@@ -1,279 +0,0 @@
-#!/usr/bin/env python3
-"""
-AMSS-NCKU GW150914 Simulation Regression Test Script
-
-Verification Requirements:
-1. XY-plane trajectory RMS error < 1% (Optimized vs. baseline, max of BH1 and BH2)
-2. ADM constraint violation < 2 (Grid Level 0)
-
-RMS Calculation Method:
- Computes trajectory deviation on the XY plane independently for BH1 and BH2
- For each black hole: RMS = sqrt((1/M) * sum((Δr_i / r_i^max)^2)) × 100%
- Final RMS = max(RMS_BH1, RMS_BH2)
-
-Usage: python3 AMSS_NCKU_Verify_ASC26.py [output_dir]
-Default: output_dir = GW150914/AMSS_NCKU_output
-Reference: GW150914-origin (baseline simulation)
-"""
-
-import numpy as np
-import sys
-import os
-
-# ANSI Color Codes
-class Color:
-    GREEN = '\033[92m'
-    RED = '\033[91m'
-    YELLOW = '\033[93m'
-    BLUE = '\033[94m'
-    BOLD = '\033[1m'
-    RESET = '\033[0m'
-
-def get_status_text(passed):
-    if passed:
-        return f"{Color.GREEN}{Color.BOLD}PASS{Color.RESET}"
-    else:
-        return f"{Color.RED}{Color.BOLD}FAIL{Color.RESET}"
-
-def load_bh_trajectory(filepath):
-    """Load black hole trajectory data"""
-    data = np.loadtxt(filepath)
-    return {
-        'time': data[:, 0],
-        'x1': data[:, 1], 'y1': data[:, 2], 'z1': data[:, 3],
-        'x2': data[:, 4], 'y2': data[:, 5], 'z2': data[:, 6]
-    }
-
-
-def load_constraint_data(filepath):
-    """Load constraint violation data"""
-    data = []
-    with open(filepath, 'r') as f:
-        for line in f:
-            if line.startswith('#'):
-                continue
-            parts = line.split()
-            if len(parts) >= 8:
-                data.append([float(x) for x in parts[:8]])
-    return np.array(data)
-
-
-def calculate_rms_error(bh_data_ref, bh_data_target):
-    """
-    Calculate trajectory-based RMS error on the XY plane between baseline and optimized simulations.
-
-    This function computes the RMS error independently for BH1 and BH2 trajectories,
-    then returns the maximum of the two as the final RMS error metric.
-
-    For each black hole, the RMS is calculated as:
-        RMS = sqrt( (1/M) * sum( (Δr_i / r_i^max)^2 ) ) × 100%
-
-    where:
-        Δr_i = sqrt((x_ref,i - x_new,i)^2 + (y_ref,i - y_new,i)^2)
-        r_i^max = max(sqrt(x_ref,i^2 + y_ref,i^2), sqrt(x_new,i^2 + y_new,i^2))
-
-    Args:
-        bh_data_ref: Reference (baseline) trajectory data
-        bh_data_target: Target (optimized) trajectory data
-
-    Returns:
-        rms_value: Final RMS error as a percentage (max of BH1 and BH2)
-        error: Error message if any
-    """
-    # Align data: truncate to the length of the shorter dataset
-    M = min(len(bh_data_ref['time']), len(bh_data_target['time']))
-
-    if M < 10:
-        return None, "Insufficient data points for comparison"
-
-    # Extract XY coordinates for both black holes
-    x1_ref = bh_data_ref['x1'][:M]
-    y1_ref = bh_data_ref['y1'][:M]
-    x2_ref = bh_data_ref['x2'][:M]
-    y2_ref = bh_data_ref['y2'][:M]
-
-    x1_new = bh_data_target['x1'][:M]
-    y1_new = bh_data_target['y1'][:M]
-    x2_new = bh_data_target['x2'][:M]
-    y2_new = bh_data_target['y2'][:M]
-
-    # Calculate RMS for BH1
-    delta_r1 = np.sqrt((x1_ref - x1_new)**2 + (y1_ref - y1_new)**2)
-    r1_ref = np.sqrt(x1_ref**2 + y1_ref**2)
-    r1_new = np.sqrt(x1_new**2 + y1_new**2)
-    r1_max = np.maximum(r1_ref, r1_new)
-
-    # Calculate RMS for BH2
-    delta_r2 = np.sqrt((x2_ref - x2_new)**2 + (y2_ref - y2_new)**2)
-    r2_ref = np.sqrt(x2_ref**2 + y2_ref**2)
-    r2_new = np.sqrt(x2_new**2 + y2_new**2)
-    r2_max = np.maximum(r2_ref, r2_new)
-
-    # Avoid division by zero for BH1
-    valid_mask1 = r1_max > 1e-15
-    if np.sum(valid_mask1) < 10:
-        return None, "Insufficient valid data points for BH1"
-
-    terms1 = (delta_r1[valid_mask1] / r1_max[valid_mask1])**2
-    rms_bh1 = np.sqrt(np.mean(terms1)) * 100
-
-    # Avoid division by zero for BH2
-    valid_mask2 = r2_max > 1e-15
-    if np.sum(valid_mask2) < 10:
-        return None, "Insufficient valid data points for BH2"
-
-    terms2 = (delta_r2[valid_mask2] / r2_max[valid_mask2])**2
-    rms_bh2 = np.sqrt(np.mean(terms2)) * 100
-
-    # Final RMS is the maximum of BH1 and BH2
-    rms_final = max(rms_bh1, rms_bh2)
-
-    return rms_final, None
-
-
-def analyze_constraint_violation(constraint_data, n_levels=9):
-    """
-    Analyze ADM constraint violation
-    Return maximum constraint violation for Grid Level 0
-    """
-    # Extract Grid Level 0 data (first entry for each time step)
-    level0_data = constraint_data[::n_levels]
-
-    # Calculate maximum absolute value for each constraint
-    results = {
-        'Ham': np.max(np.abs(level0_data[:, 1])),
-        'Px': np.max(np.abs(level0_data[:, 2])),
-        'Py': np.max(np.abs(level0_data[:, 3])),
-        'Pz': np.max(np.abs(level0_data[:, 4])),
-        'Gx': np.max(np.abs(level0_data[:, 5])),
-        'Gy': np.max(np.abs(level0_data[:, 6])),
-        'Gz': np.max(np.abs(level0_data[:, 7]))
-    }
-
-    results['max_violation'] = max(results.values())
-    return results
-
-
-def print_header():
-    """Print report header"""
-    print("\n" + Color.BLUE + Color.BOLD + "=" * 65 + Color.RESET)
-    print(Color.BOLD + "   AMSS-NCKU GW150914 Simulation Regression Test Report" + Color.RESET)
-    print(Color.BLUE + Color.BOLD + "=" * 65 + Color.RESET)
-
-
-def print_rms_results(rms_rel, error, threshold=1.0):
-    """Print RMS error results"""
-    print(f"\n{Color.BOLD}1. RMS Error Analysis (Baseline vs Optimized){Color.RESET}")
-    print("-" * 45)
-
-    if error:
-        print(f"   {Color.RED}Error: {error}{Color.RESET}")
-        return False
-
-    passed = rms_rel < threshold
-
-    print(f"   RMS relative error: {rms_rel:.4f}%")
-    print(f"   Requirement:        < {threshold}%")
-    print(f"   Status:             {get_status_text(passed)}")
-
-    return passed
-
-
-def print_constraint_results(results, threshold=2.0):
-    """Print constraint violation results"""
-    print(f"\n{Color.BOLD}2. ADM Constraint Violation Analysis (Grid Level 0){Color.RESET}")
-    print("-" * 45)
-
-    names = ['Ham', 'Px', 'Py', 'Pz', 'Gx', 'Gy', 'Gz']
-    for i, name in enumerate(names):
-        print(f"   Max |{name:3}|: {results[name]:.6f}", end="   ")
-        if (i + 1) % 2 == 0: print()
-    if len(names) % 2 != 0: print()
-
-    passed = results['max_violation'] < threshold
-    
-    print(f"\n   Maximum violation:  {results['max_violation']:.6f}")
-    print(f"   Requirement:        < {threshold}")
-    print(f"   Status:             {get_status_text(passed)}")
-
-    return passed
-
-
-def print_summary(rms_passed, constraint_passed):
-    """Print summary"""
-    print("\n" + Color.BLUE + Color.BOLD + "=" * 65 + Color.RESET)
-    print(Color.BOLD + "Verification Summary" + Color.RESET)
-    print(Color.BLUE + Color.BOLD + "=" * 65 + Color.RESET)
-
-    all_passed = rms_passed and constraint_passed
-    
-    res_rms = get_status_text(rms_passed)
-    res_con = get_status_text(constraint_passed)
-
-    print(f"   [1] RMS trajectory check:         {res_rms}")
-    print(f"   [2] ADM constraint check:         {res_con}")
-    
-    final_status = f"{Color.GREEN}{Color.BOLD}ALL CHECKS PASSED{Color.RESET}" if all_passed else f"{Color.RED}{Color.BOLD}SOME CHECKS FAILED{Color.RESET}"
-    print(f"\n   Overall result: {final_status}")
-    print(Color.BLUE + Color.BOLD + "=" * 65 + Color.RESET + "\n")
-
-    return all_passed
-
-
-def main():
-    # Determine target (optimized) output directory
-    if len(sys.argv) > 1:
-        target_dir = sys.argv[1]
-    else:
-        script_dir = os.path.dirname(os.path.abspath(__file__))
-        target_dir = os.path.join(script_dir, "GW150914/AMSS_NCKU_output")
-
-    # Determine reference (baseline) directory
-    script_dir = os.path.dirname(os.path.abspath(__file__))
-    reference_dir = os.path.join(script_dir, "GW150914-origin/AMSS_NCKU_output")
-
-    # Data file paths
-    bh_file_ref = os.path.join(reference_dir, "bssn_BH.dat")
-    bh_file_target = os.path.join(target_dir, "bssn_BH.dat")
-    constraint_file = os.path.join(target_dir, "bssn_constraint.dat")
-
-    # Check if files exist
-    if not os.path.exists(bh_file_ref):
-        print(f"{Color.RED}{Color.BOLD}Error:{Color.RESET} Baseline trajectory file not found: {bh_file_ref}")
-        sys.exit(1)
-
-    if not os.path.exists(bh_file_target):
-        print(f"{Color.RED}{Color.BOLD}Error:{Color.RESET} Target trajectory file not found: {bh_file_target}")
-        sys.exit(1)
-
-    if not os.path.exists(constraint_file):
-        print(f"{Color.RED}{Color.BOLD}Error:{Color.RESET} Constraint data file not found: {constraint_file}")
-        sys.exit(1)
-
-    # Print header
-    print_header()
-    print(f"\n{Color.BOLD}Reference (Baseline):{Color.RESET} {Color.BLUE}{reference_dir}{Color.RESET}")
-    print(f"{Color.BOLD}Target (Optimized):  {Color.RESET} {Color.BLUE}{target_dir}{Color.RESET}")
-
-    # Load data
-    bh_data_ref = load_bh_trajectory(bh_file_ref)
-    bh_data_target = load_bh_trajectory(bh_file_target)
-    constraint_data = load_constraint_data(constraint_file)
-
-    # Calculate RMS error
-    rms_rel, error = calculate_rms_error(bh_data_ref, bh_data_target)
-    rms_passed = print_rms_results(rms_rel, error)
-
-    # Analyze constraint violation
-    constraint_results = analyze_constraint_violation(constraint_data)
-    constraint_passed = print_constraint_results(constraint_results)
-
-    # Print summary
-    all_passed = print_summary(rms_passed, constraint_passed)
-
-    # Return exit code
-    sys.exit(0 if all_passed else 1)
-
-
-if __name__ == "__main__":
-    main()
--- a/AMSS_NCKU_source/FFT.f90
+++ b/AMSS_NCKU_source/FFT.f90
@@ -37,51 +37,57 @@ 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, intent(in) :: isign, nn
-DOUBLE PRECISION, dimension(2*nn), intent(inout) :: dataa
-
-type(DFTI_DESCRIPTOR), pointer :: desc
-integer :: status
-
-! Create DFTI descriptor for 1D complex-to-complex transform
-status = DftiCreateDescriptor(desc, DFTI_DOUBLE, DFTI_COMPLEX, 1, nn)
-if (status /= 0) return
-
-! Set input/output storage as interleaved complex (default)
-status = DftiSetValue(desc, DFTI_PLACEMENT, DFTI_INPLACE)
-if (status /= 0) then
-   status = DftiFreeDescriptor(desc)
-   return
+INTEGER::isign,nn
+double precision,dimension(2*nn)::dataa
+INTEGER::i,istep,j,m,mmax,n
+double precision::tempi,tempr
+DOUBLE PRECISION::theta,wi,wpi,wpr,wr,wtemp
+n=2*nn
+j=1
+do i=1,n,2
+  if(j.gt.i)then
+     tempr=dataa(j)
+     tempi=dataa(j+1)
+     dataa(j)=dataa(i)
+     dataa(j+1)=dataa(i+1)
+     dataa(i)=tempr
+     dataa(i+1)=tempi
+  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
--- a/AMSS_NCKU_source/TwoPunctures.C
+++ b/AMSS_NCKU_source/TwoPunctures.C
@@ -27,7 +27,6 @@ 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,
@@ -892,17 +891,25 @@ double TwoPunctures::norm1(double *v, int n)
 /* -------------------------------------------------------------------------*/
 double TwoPunctures::norm2(double *v, int n)
 {
-  // Optimized with oneMKL BLAS DNRM2
-  // Computes: sqrt(sum(v[i]^2))
-  return cblas_dnrm2(n, v, 1);
+  int i;
+  double result = 0;
+
+  for (i = 0; i < n; i++)
+    result += v[i] * v[i];
+
+  return sqrt(result);
 }

 /* -------------------------------------------------------------------------*/
 double TwoPunctures::scalarproduct(double *v, double *w, int n)
 {
-  // Optimized with oneMKL BLAS DDOT
-  // Computes: sum(v[i] * w[i])
-  return cblas_ddot(n, v, 1, w, 1);
+  int i;
+  double result = 0;
+
+  for (i = 0; i < n; i++)
+    result += v[i] * w[i];
+
+  return result;
 }

 /* -------------------------------------------------------------------------*/
--- a/AMSS_NCKU_source/diff_new.f90
+++ b/AMSS_NCKU_source/diff_new.f90
@@ -997,11 +997,10 @@
  fy = ZEO
  fz = ZEO

-  !$omp parallel do collapse(3) private(i,j,k) if(ex(1)*ex(2)*ex(3) > 4096)
  do k=1,ex(3)-1
  do j=1,ex(2)-1
  do i=1,ex(1)-1
-#if 0
+#if 0  
 ! x direction   
        if(i+2 <= imax .and. i-2 >= imin)then
 !
@@ -1152,11 +1151,10 @@

  fx = ZEO

-  !$omp parallel do collapse(3) private(i,j,k) if(ex(1)*ex(2)*ex(3) > 4096)
  do k=1,ex(3)-1
  do j=1,ex(2)-1
  do i=1,ex(1)-1
-! x direction
+! x direction   
        if(i+2 <= imax .and. i-2 >= imin)then
 !
 !              f(i-2) - 8 f(i-1) + 8 f(i+1) - f(i+2)
@@ -1229,11 +1227,10 @@

  fy = ZEO

-  !$omp parallel do collapse(3) private(i,j,k) if(ex(1)*ex(2)*ex(3) > 4096)
  do k=1,ex(3)-1
  do j=1,ex(2)-1
  do i=1,ex(1)-1
-! y direction
+! y direction   
        if(j+2 <= jmax .and. j-2 >= jmin)then

      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))
@@ -1300,11 +1297,10 @@

  fz = ZEO

-  !$omp parallel do collapse(3) private(i,j,k) if(ex(1)*ex(2)*ex(3) > 4096)
  do k=1,ex(3)-1
  do j=1,ex(2)-1
  do i=1,ex(1)-1
-! z direction
+! z direction   
        if(k+2 <= kmax .and. k-2 >= kmin)then

      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))
@@ -1405,11 +1401,10 @@
  fxz = ZEO
  fyz = ZEO

-  !$omp parallel do collapse(3) private(i,j,k) if(ex(1)*ex(2)*ex(3) > 4096)
  do k=1,ex(3)-1
  do j=1,ex(2)-1
  do i=1,ex(1)-1
-#if 0
+#if 0  
 !~~~~~~ fxx
        if(i+2 <= imax .and. i-2 >= imin)then
 !
@@ -1581,7 +1576,6 @@

  fxx = ZEO

-  !$omp parallel do collapse(3) private(i,j,k) if(ex(1)*ex(2)*ex(3) > 4096)
  do k=1,ex(3)-1
  do j=1,ex(2)-1
  do i=1,ex(1)-1
@@ -1649,7 +1643,6 @@

  fyy = ZEO

-  !$omp parallel do collapse(3) private(i,j,k) if(ex(1)*ex(2)*ex(3) > 4096)
  do k=1,ex(3)-1
  do j=1,ex(2)-1
  do i=1,ex(1)-1
@@ -1719,7 +1712,6 @@

  fzz = ZEO

-  !$omp parallel do collapse(3) private(i,j,k) if(ex(1)*ex(2)*ex(3) > 4096)
  do k=1,ex(3)-1
  do j=1,ex(2)-1
  do i=1,ex(1)-1
@@ -1789,7 +1781,6 @@

  fxy = ZEO

-  !$omp parallel do collapse(3) private(i,j,k) if(ex(1)*ex(2)*ex(3) > 4096)
  do k=1,ex(3)-1
  do j=1,ex(2)-1
  do i=1,ex(1)-1
@@ -1860,7 +1851,6 @@

  fxz = ZEO

-  !$omp parallel do collapse(3) private(i,j,k) if(ex(1)*ex(2)*ex(3) > 4096)
  do k=1,ex(3)-1
  do j=1,ex(2)-1
  do i=1,ex(1)-1
@@ -1929,7 +1919,6 @@

  fyz = ZEO

-  !$omp parallel do collapse(3) private(i,j,k) if(ex(1)*ex(2)*ex(3) > 4096)
  do k=1,ex(3)-1
  do j=1,ex(2)-1
  do i=1,ex(1)-1
--- a/AMSS_NCKU_source/diff_new_sh.f90
+++ b/AMSS_NCKU_source/diff_new_sh.f90
@@ -1019,11 +1019,10 @@
  fy = ZEO
  fz = ZEO

-  !$omp parallel do collapse(3) private(i,j,k) if(ex(1)*ex(2)*ex(3) > 4096)
  do k=1,ex(3)
  do j=1,ex(2)
  do i=1,ex(1)
-! x direction
+! x direction   
        if(i+2 <= imax .and. i-2 >= imin)then
 !
 !              f(i-2) - 8 f(i-1) + 8 f(i+1) - f(i+2)
@@ -1135,11 +1134,10 @@

  fx = ZEO

-  !$omp parallel do collapse(3) private(i,j,k) if(ex(1)*ex(2)*ex(3) > 4096)
  do k=1,ex(3)
  do j=1,ex(2)
  do i=1,ex(1)
-! x direction
+! x direction   
        if(i+2 <= imax .and. i-2 >= imin)then
 !
 !              f(i-2) - 8 f(i-1) + 8 f(i+1) - f(i+2)
@@ -1229,11 +1227,10 @@

  fy = ZEO

-  !$omp parallel do collapse(3) private(i,j,k) if(ex(1)*ex(2)*ex(3) > 4096)
  do k=1,ex(3)
  do j=1,ex(2)
  do i=1,ex(1)
-! y direction
+! y direction   
        if(j+2 <= jmax .and. j-2 >= jmin)then

      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))
@@ -1317,11 +1314,10 @@

  fz = ZEO

-  !$omp parallel do collapse(3) private(i,j,k) if(ex(1)*ex(2)*ex(3) > 4096)
  do k=1,ex(3)
  do j=1,ex(2)
  do i=1,ex(1)
-! z direction
+! z direction   
        if(k+2 <= kmax .and. k-2 >= kmin)then

      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))
@@ -1434,7 +1430,6 @@
  fxz = ZEO
  fyz = ZEO

-  !$omp parallel do collapse(3) private(i,j,k) if(ex(1)*ex(2)*ex(3) > 4096)
  do k=1,ex(3)
  do j=1,ex(2)
  do i=1,ex(1)
@@ -1585,7 +1580,6 @@

  fxx = ZEO

-  !$omp parallel do collapse(3) private(i,j,k) if(ex(1)*ex(2)*ex(3) > 4096)
  do k=1,ex(3)
  do j=1,ex(2)
  do i=1,ex(1)
@@ -1665,7 +1659,6 @@

  fyy = ZEO

-  !$omp parallel do collapse(3) private(i,j,k) if(ex(1)*ex(2)*ex(3) > 4096)
  do k=1,ex(3)
  do j=1,ex(2)
  do i=1,ex(1)
@@ -1747,7 +1740,6 @@

  fzz = ZEO

-  !$omp parallel do collapse(3) private(i,j,k) if(ex(1)*ex(2)*ex(3) > 4096)
  do k=1,ex(3)
  do j=1,ex(2)
  do i=1,ex(1)
@@ -1829,7 +1821,6 @@

  fxy = ZEO

-  !$omp parallel do collapse(3) private(i,j,k) if(ex(1)*ex(2)*ex(3) > 4096)
  do k=1,ex(3)
  do j=1,ex(2)
  do i=1,ex(1)
@@ -1912,7 +1903,6 @@

  fxz = ZEO

-  !$omp parallel do collapse(3) private(i,j,k) if(ex(1)*ex(2)*ex(3) > 4096)
  do k=1,ex(3)
  do j=1,ex(2)
  do i=1,ex(1)
@@ -1993,7 +1983,6 @@

  fyz = ZEO

-  !$omp parallel do collapse(3) private(i,j,k) if(ex(1)*ex(2)*ex(3) > 4096)
  do k=1,ex(3)
  do j=1,ex(2)
  do i=1,ex(1)
--- a/AMSS_NCKU_source/diff_newwb.f90
+++ b/AMSS_NCKU_source/diff_newwb.f90
@@ -1186,11 +1186,10 @@
  fy = ZEO
  fz = ZEO

-  !$omp parallel do collapse(3) private(i,j,k) if(ex(1)*ex(2)*ex(3) > 4096)
  do k=1,ex(3)
  do j=1,ex(2)
  do i=1,ex(1)
-! x direction
+! x direction   
        if(i+2 <= imax .and. i-2 >= imin)then
 !
 !              f(i-2) - 8 f(i-1) + 8 f(i+1) - f(i+2)
@@ -1301,11 +1300,10 @@

  fx = ZEO

-  !$omp parallel do collapse(3) private(i,j,k) if(ex(1)*ex(2)*ex(3) > 4096)
  do k=1,ex(3)
  do j=1,ex(2)
  do i=1,ex(1)
-! x direction
+! x direction   
        if(i+2 <= imax .and. i-2 >= imin)then
 !
 !              f(i-2) - 8 f(i-1) + 8 f(i+1) - f(i+2)
@@ -1383,11 +1381,10 @@

  fy = ZEO

-  !$omp parallel do collapse(3) private(i,j,k) if(ex(1)*ex(2)*ex(3) > 4096)
  do k=1,ex(3)
  do j=1,ex(2)
  do i=1,ex(1)
-! y direction
+! y direction   
        if(j+2 <= jmax .and. j-2 >= jmin)then

      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))
@@ -1459,11 +1456,10 @@

  fz = ZEO

-  !$omp parallel do collapse(3) private(i,j,k) if(ex(1)*ex(2)*ex(3) > 4096)
  do k=1,ex(3)
  do j=1,ex(2)
  do i=1,ex(1)
-! z direction
+! z direction   
        if(k+2 <= kmax .and. k-2 >= kmin)then

      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))
@@ -1569,7 +1565,6 @@
  fxz = ZEO
  fyz = ZEO

-  !$omp parallel do collapse(3) private(i,j,k) if(ex(1)*ex(2)*ex(3) > 4096)
  do k=1,ex(3)
  do j=1,ex(2)
  do i=1,ex(1)
@@ -1786,7 +1781,6 @@

  fxx = ZEO

-  !$omp parallel do collapse(3) private(i,j,k) if(ex(1)*ex(2)*ex(3) > 4096)
  do k=1,ex(3)
  do j=1,ex(2)
  do i=1,ex(1)
@@ -1862,7 +1856,6 @@

  fyy = ZEO

-  !$omp parallel do collapse(3) private(i,j,k) if(ex(1)*ex(2)*ex(3) > 4096)
  do k=1,ex(3)
  do j=1,ex(2)
  do i=1,ex(1)
@@ -1940,7 +1933,6 @@

  fzz = ZEO

-  !$omp parallel do collapse(3) private(i,j,k) if(ex(1)*ex(2)*ex(3) > 4096)
  do k=1,ex(3)
  do j=1,ex(2)
  do i=1,ex(1)
@@ -2018,7 +2010,6 @@

  fxy = ZEO

-  !$omp parallel do collapse(3) private(i,j,k) if(ex(1)*ex(2)*ex(3) > 4096)
  do k=1,ex(3)
  do j=1,ex(2)
  do i=1,ex(1)
@@ -2107,7 +2098,6 @@

  fxz = ZEO

-  !$omp parallel do collapse(3) private(i,j,k) if(ex(1)*ex(2)*ex(3) > 4096)
  do k=1,ex(3)
  do j=1,ex(2)
  do i=1,ex(1)
@@ -2194,7 +2184,6 @@

  fyz = ZEO

-  !$omp parallel do collapse(3) private(i,j,k) if(ex(1)*ex(2)*ex(3) > 4096)
  do k=1,ex(3)
  do j=1,ex(2)
  do i=1,ex(1)
--- a/AMSS_NCKU_source/fmisc.f90
+++ b/AMSS_NCKU_source/fmisc.f90
@@ -1259,7 +1259,7 @@ end subroutine d2dump

  end subroutine polin3
 !--------------------------------------------------------------------------------------
-! calculate L2norm
+! calculate L2norm  
  subroutine l2normhelper(ex, X, Y, Z,xmin,ymin,zmin,xmax,ymax,zmax,&
                          f,f_out,gw)

@@ -1276,9 +1276,7 @@ end subroutine d2dump
  real*8            :: dX, dY, dZ
  integer::imin,jmin,kmin
  integer::imax,jmax,kmax
-  integer::i,j,k,n_elements
-  real*8, dimension(:), allocatable :: f_flat
-  real*8, external :: DDOT
+  integer::i,j,k

  dX = X(2) - X(1)
  dY = Y(2) - Y(1)
@@ -1302,12 +1300,7 @@ if(dabs(X(1)-xmin) < dX) imin = 1
 if(dabs(Y(1)-ymin) < dY) jmin = 1
 if(dabs(Z(1)-zmin) < dZ) kmin = 1

-! 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 = sum(f(imin:imax,jmin:jmax,kmin:kmax)*f(imin:imax,jmin:jmax,kmin:kmax))

 f_out = f_out*dX*dY*dZ

@@ -1332,9 +1325,7 @@ f_out = f_out*dX*dY*dZ
  real*8            :: dX, dY, dZ
  integer::imin,jmin,kmin
  integer::imax,jmax,kmax
-  integer::i,j,k,n_elements
-  real*8, dimension(:), allocatable :: f_flat
-  real*8, external :: DDOT
+  integer::i,j,k

  real*8 :: PIo4

@@ -1397,12 +1388,7 @@ if(Symmetry==2)then
  if(dabs(ymin+gw*dY)<dY.and.Y(1)<0.d0) jmin = gw+1
 endif

-! 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 = sum(f(imin:imax,jmin:jmax,kmin:kmax)*f(imin:imax,jmin:jmax,kmin:kmax))

 f_out = f_out*dX*dY*dZ

@@ -1430,8 +1416,6 @@ 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

@@ -1494,12 +1478,11 @@ if(Symmetry==2)then
  if(dabs(ymin+gw*dY)<dY.and.Y(1)<0.d0) jmin = gw+1
 endif

-! Optimized with oneMKL BLAS DDOT for dot product
+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

@@ -1697,7 +1680,6 @@ deallocate(f_flat)
  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  
@@ -1733,21 +1715,20 @@ deallocate(f_flat)
     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

-  ! Third dimension: x-direction weighted sum using BLAS DDOT
-  f_int = DDOT(ORDN, coef(1:ORDN), 1, tmp1, 1)
+  f_int=0
+  do m=1,ORDN
+    f_int = f_int + coef(m)*tmp1(m)
+  enddo

  return

@@ -1777,7 +1758,6 @@ deallocate(f_flat)
  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  
@@ -1807,14 +1787,15 @@ deallocate(f_flat)
     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

-  ! Use BLAS DDOT for final weighted sum
-  f_int = DDOT(ORDN, coef(1:ORDN), 1, tmp1, 1)
+  f_int=0
+  do m=1,ORDN
+    f_int = f_int + coef(m)*tmp1(m)
+  enddo

  return

@@ -1845,7 +1826,6 @@ deallocate(f_flat)
  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
@@ -1906,8 +1886,10 @@ deallocate(f_flat)
          write(*,*)"error in global_interpind1d, not recognized dumyd = ",dumyd
  endif

-  ! Optimized with BLAS DDOT for weighted sum
-  f_int = DDOT(ORDN, coef, 1, ya, 1)
+  f_int=0
+  do m=1,ORDN
+    f_int = f_int + coef(m)*ya(m)
+  enddo

  return

@@ -2139,38 +2121,24 @@ deallocate(f_flat)

  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 ]
+  integer :: i

 ! sanity check
  if(N < 0)then
     write(*,*) "ffact: error input for factorial"
-     gont = 1.d0
     return
  endif

-  ! Use lookup table for small N (fast path)
-  if(N <= 20)then
-     gont = fact_table(N)
-  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
+  gont = 1.d0
+  do i=1,N
+     gont = gont*i
+  enddo

  return

--- a/AMSS_NCKU_source/gaussj.C
+++ b/AMSS_NCKU_source/gaussj.C
@@ -16,66 +16,115 @@ using namespace std;
 #include <string.h>
 #include <math.h>
 #endif
-
-// Intel oneMKL LAPACK interface
-#include <mkl_lapacke.h>
-/* Linear equation solution using Intel oneMKL LAPACK.
+/* Linear equation solution by Gauss-Jordan elimination.
 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.
-
-Mathematical equivalence:
-  Solves: A * x = b  =>  x = A^(-1) * b
-  Original Gauss-Jordan and LAPACK dgesv/dgetri produce identical results
-  within numerical precision. */
+corresponding set of solution vectors */

 int gaussj(double *a, double *b, int n)
 {
-  // Allocate pivot array and workspace
-  lapack_int *ipiv = new lapack_int[n];
-  lapack_int info;
+  double swap;

-  // Make a copy of matrix a for solving (dgesv modifies it to LU form)
-  double *a_copy = new double[n * n];
-  for (int i = 0; i < n * n; i++) {
-    a_copy[i] = a[i];
+  int *indxc, *indxr, *ipiv;
+  indxc = new int[n];
+  indxr = new int[n];
+  ipiv = new int[n];
+
+  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;
+    if (irow != icol)
+    {
+      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];
+      b[irow] = b[icol];
+      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;
+      }
  }

-  // Step 1: Solve linear system A*x = b using LU decomposition
-  // 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);
-
-  if (info != 0) {
-    cout << "gaussj: Singular Matrix (dgesv info=" << info << ")" << endl;
-    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;
+  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 0;
 }
--- a/AMSS_NCKU_source/ilucg.f90
+++ b/AMSS_NCKU_source/ilucg.f90
@@ -512,10 +512,11 @@
      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)

-      DOUBLE PRECISION, EXTERNAL :: DDOT
-      DGVV = DDOT(N, V, 1, W, 1)
+      SUM = 0.0D0
+            DO 10 I = 1,N
+            SUM = SUM + V(I)*W(I)
+10          CONTINUE
+      DGVV = SUM
      RETURN
      END
--- a/AMSS_NCKU_source/kodiss.f90
+++ b/AMSS_NCKU_source/kodiss.f90
@@ -159,7 +159,6 @@ integer, parameter :: NO_SYMM=0, OCTANT=2

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

-  !$omp parallel do collapse(3) private(i,j,k) if(ex(1)*ex(2)*ex(3) > 4096)
  do k=1,ex(3)
  do j=1,ex(2)
  do i=1,ex(1)
--- a/AMSS_NCKU_source/kodiss_sh.f90
+++ b/AMSS_NCKU_source/kodiss_sh.f90
@@ -369,12 +369,11 @@ integer, parameter :: NO_SYMM=0, EQ_SYMM=1, OCTANT=2

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

-  !$omp parallel do collapse(3) private(i,j,k) if(ex(1)*ex(2)*ex(3) > 4096)
  do k=1,ex(3)
  do j=1,ex(2)
  do i=1,ex(1)

-#if 1
+#if 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
--- a/AMSS_NCKU_source/lopsidediff.f90
+++ b/AMSS_NCKU_source/lopsidediff.f90
@@ -231,9 +231,8 @@ subroutine lopsided(ex,X,Y,Z,f,f_rhs,Sfx,Sfy,Sfz,Symmetry,SoA)

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

-! upper bound set ex-1 only for efficiency,
+! upper bound set ex-1 only for efficiency, 
 ! the loop body will set ex 0 also
-  !$omp parallel do collapse(3) private(i,j,k) if(ex(1)*ex(2)*ex(3) > 4096)
  do k=1,ex(3)-1
  do j=1,ex(2)-1
  do i=1,ex(1)-1
--- a/AMSS_NCKU_source/macrodef.h
+++ b/AMSS_NCKU_source/macrodef.h
@@ -2,7 +2,7 @@
 #ifndef MICRODEF_H
 #define MICRODEF_H

-#include "macrodef.fh"
+#include "microdef.fh"

 // application parameters

--- a/AMSS_NCKU_source/makefile.inc
+++ b/AMSS_NCKU_source/makefile.inc
@@ -1,30 +1,19 @@
-## 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

-## Intel oneAPI version with oneMKL (Optimized for performance)
-filein  = -I/usr/include/ -I${MKLROOT}/include
+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/

-## 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 -qopenmp
+## LDLIBS  = -L/usr/lib/x86_64-linux-gnu -lmpich -L/usr/lib64 -L/usr/lib/gcc/x86_64-linux-gnu/11 -lgfortran
+LDLIBS  = -L/usr/lib/x86_64-linux-gnu -L/usr/lib64 -L/usr/lib/gcc/x86_64-linux-gnu/11 -lgfortran -lmpi -lgfortran

-## Aggressive optimization flags:
-## -O3: Maximum optimization
-## -xHost: Optimize for the host CPU architecture (Intel/AMD compatible)
-## -fp-model fast=2: Aggressive floating-point optimizations
-## -fma: Enable fused multiply-add instructions
-## OpenMP re-enabled for MPI+OpenMP hybrid parallelism (MKL stays sequential to avoid nested parallelism)
-CXXAPPFLAGS  = -O3 -xHost -fp-model fast=2 -fma -qopenmp \
-               -Dfortran3 -Dnewc -I${MKLROOT}/include
-f90appflags  = -O3 -xHost -fp-model fast=2 -fma -qopenmp \
-               -fpp -I${MKLROOT}/include
-f90          = ifx
-f77          = ifx
-CXX          = icpx
-CC           = icx
-CLINKER      = mpiicpx 
+CXXAPPFLAGS  = -O0 -Wno-deprecated -Dfortran3 -Dnewc
+#f90appflags = -O0 -fpp
+f90appflags  = -O0 -x f95-cpp-input
+f90          = gfortran 
+f77          = gfortran 
+CXX          = g++
+CC           = gcc
+CLINKER      = mpic++ 

 Cu = nvcc
 CUDA_LIB_PATH = -L/usr/lib/cuda/lib64 -I/usr/include -I/usr/lib/cuda/include
--- a/makefile_and_run.py
+++ b/makefile_and_run.py
@@ -11,13 +11,6 @@
 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
-NUMACTL_CPU_BIND = "taskset -c 0-111"
-
-## Build parallelism configuration
-BUILD_JOBS = 104
-

 ##################################################################

@@ -33,11 +26,11 @@ def makefile_ABE():
    print( " Compiling the AMSS-NCKU executable file ABE/ABEGPU " ) 
    print(                                                        )

-    ## Build command with CPU binding to nohz_full cores
+    ## Build command
    if (input_data.GPU_Calculation == "no"):
-        makefile_command  = f"{NUMACTL_CPU_BIND} make -j{BUILD_JOBS} ABE"
+        makefile_command  = "make -j96" + " ABE"
    elif (input_data.GPU_Calculation == "yes"):
-        makefile_command  = f"{NUMACTL_CPU_BIND} make -j{BUILD_JOBS} ABEGPU"
+        makefile_command  = "make -j4" + " ABEGPU"
    else:
        print( " CPU/GPU numerical calculation setting is wrong " )
        print(                                                    )
@@ -74,8 +67,8 @@ def makefile_TwoPunctureABE():
    print( " Compiling the AMSS-NCKU executable file TwoPunctureABE " )
    print(                                                            )
    
-    ## Build command with CPU binding to nohz_full cores
-    makefile_command = f"{NUMACTL_CPU_BIND} make -j{BUILD_JOBS} TwoPunctureABE"
+    ## Build command
+    makefile_command = "make" + " 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) 
@@ -110,18 +103,12 @@ def run_ABE():
    print(                                                      )

    ## Define the command to run; cast other values to strings as needed
-    ## MPI+OpenMP hybrid: compute threads per rank from total cores / MPI ranks
-    omp_threads = max(1, 96 // input_data.MPI_processes)
-    omp_env = (f" -genv OMP_NUM_THREADS={omp_threads}"
-               f" -genv OMP_PROC_BIND=close"
-               f" -genv OMP_PLACES=cores"
-               f" -genv I_MPI_PIN_DOMAIN=omp")
-
+    
    if (input_data.GPU_Calculation == "no"):
-        mpi_command         = NUMACTL_CPU_BIND + " mpirun -np " + str(input_data.MPI_processes) + omp_env + " ./ABE"
+        mpi_command         = "mpirun -np " + str(input_data.MPI_processes) + " ./ABE"
        mpi_command_outfile = "ABE_out.log"
    elif (input_data.GPU_Calculation == "yes"):
-        mpi_command         = NUMACTL_CPU_BIND + " mpirun -np " + str(input_data.MPI_processes) + omp_env + " ./ABEGPU"
+        mpi_command         = "mpirun -np " + str(input_data.MPI_processes) + " ./ABEGPU"
        mpi_command_outfile = "ABEGPU_out.log"
 
    ## Execute the MPI command and stream output
@@ -160,7 +147,7 @@ def run_TwoPunctureABE():
    print(                                                          )
    
    ## Define the command to run
-    TwoPuncture_command         = NUMACTL_CPU_BIND + " ./TwoPunctureABE"
+    TwoPuncture_command         = "./TwoPunctureABE"
    TwoPuncture_command_outfile = "TwoPunctureABE_out.log"

    ## Execute the command with subprocess.Popen and stream output