taskset setting updated

Add _fh variants of fderivs, fdderivs, kodis, and lopsided that accept
a caller-provided fh work array instead of allocating one internally.
Declare two shared work arrays in compute_rhs_bssn (fh_work2 for
symmetry_bd(2) callers, fh_work3 for symmetry_bd(3) callers) and pass
them to all ~84 subroutine calls, eliminating ~77 redundant automatic
array allocations (~591 MB churn per RHS call, ~2.3 GB per timestep).

Co-Authored-By: Claude Opus 4.6 <noreply@anthropic.com>

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)

AMSS_NCKU_source/ABE.optrpt

@@ -1,710 +0,0 @@
-Begin optimization report for: __cxx_global_var_init
-=================================================================
-
-Begin optimization report for: std::map<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>, int, std::less<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>, std::allocator<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, int>>>::map()
-=================================================================
-
-Begin optimization report for: std::map<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>, int, std::less<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>, std::allocator<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, int>>>::~map()
-=================================================================
-
-Begin optimization report for: __cxx_global_var_init.1
-=================================================================
-
-Begin optimization report for: std::map<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>, double, std::less<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>, std::allocator<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, double>>>::map()
-=================================================================
-
-Begin optimization report for: std::map<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>, double, std::less<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>, std::allocator<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, double>>>::~map()
-=================================================================
-
-Begin optimization report for: __cxx_global_var_init.2
-=================================================================
-
-Begin optimization report for: std::map<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>, std::less<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>, std::allocator<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>>>::map()
-=================================================================
-
-Begin optimization report for: std::map<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>, std::less<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>, std::allocator<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>>>::~map()
-=================================================================
-
-Begin optimization report for: main
-
-LOOP BEGIN at /usr/lib/gcc/x86_64-redhat-linux/14/../../../../include/c++/14/bits/stl_tree.h (1947, 7)
-LOOP END
-
-LOOP BEGIN at ABE.C (125, 13)
-LOOP END
-=================================================================
-
-Begin optimization report for: std::allocator<char>::allocator()
-=================================================================
-
-Begin optimization report for: std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>::basic_string<std::allocator<char>>(char const*, std::allocator<char> const&)
-=================================================================
-
-Begin optimization report for: std::__new_allocator<char>::~__new_allocator()
-=================================================================
-
-Begin optimization report for: std::map<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>, std::less<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>, std::allocator<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>>>::insert(std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>&&)
-=================================================================
-
-Begin optimization report for: std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>::pair<char const (&) [9], std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>&, true>(char const (&) [9], std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>&)
-=================================================================
-
-Begin optimization report for: std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>::~pair()
-=================================================================
-
-Begin optimization report for: std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>::~basic_string()
-=================================================================
-
-Begin optimization report for: std::_Rb_tree_iterator<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>>::_Rb_tree_iterator()
-=================================================================
-
-Begin optimization report for: std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>::basic_string()
-=================================================================
-
-Begin optimization report for: std::map<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>, std::less<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>, std::allocator<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>>>::find(std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const&)
-
-LOOP BEGIN at /usr/lib/gcc/x86_64-redhat-linux/14/../../../../include/c++/14/bits/stl_tree.h (1947, 7)
-LOOP END
-=================================================================
-
-Begin optimization report for: std::operator!=(std::_Rb_tree_iterator<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>> const&, std::_Rb_tree_iterator<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>> const&)
-=================================================================
-
-Begin optimization report for: std::map<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>, std::less<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>, std::allocator<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>>>::end()
-=================================================================
-
-Begin optimization report for: std::_Rb_tree_iterator<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>>::operator->() const
-=================================================================
-
-Begin optimization report for: std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>::operator=(std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const&)
-=================================================================
-
-Begin optimization report for: std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>::c_str() const
-=================================================================
-
-Begin optimization report for: std::basic_ostream<char, std::char_traits<char>>& std::operator<<<std::char_traits<char>>(std::basic_ostream<char, std::char_traits<char>>&, char const*)
-=================================================================
-
-Begin optimization report for: std::ostream::operator<<(std::ostream& (*)(std::ostream&))
-=================================================================
-
-Begin optimization report for: std::basic_ostream<char, std::char_traits<char>>& std::endl<char, std::char_traits<char>>(std::basic_ostream<char, std::char_traits<char>>&)
-=================================================================
-
-Begin optimization report for: std::basic_ifstream<char, std::char_traits<char>>::basic_ifstream(char const*, std::_Ios_Openmode)
-=================================================================
-
-Begin optimization report for: std::basic_ios<char, std::char_traits<char>>::good() const
-=================================================================
-
-Begin optimization report for: std::istream::getline(char*, long)
-=================================================================
-
-Begin optimization report for: std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>::operator=(char const*)
-=================================================================
-
-Begin optimization report for: std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>::basic_string(std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const&)
-=================================================================
-
-Begin optimization report for: bool std::operator==<char, std::char_traits<char>, std::allocator<char>>(std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const&, char const*)
-=================================================================
-
-Begin optimization report for: atoi
-=================================================================
-
-Begin optimization report for: atof
-=================================================================
-
-Begin optimization report for: std::basic_ifstream<char, std::char_traits<char>>::close()
-=================================================================
-
-Begin optimization report for: std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>::pair<char const (&) [11], std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>&, true>(char const (&) [11], std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>&)
-=================================================================
-
-Begin optimization report for: std::basic_ifstream<char, std::char_traits<char>>::~basic_ifstream()
-=================================================================
-
-Begin optimization report for: std::basic_ofstream<char, std::char_traits<char>>::basic_ofstream()
-=================================================================
-
-Begin optimization report for: std::basic_ofstream<char, std::char_traits<char>>::open(char const*, std::_Ios_Openmode)
-=================================================================
-
-Begin optimization report for: std::basic_ofstream<char, std::char_traits<char>>::close()
-=================================================================
-
-Begin optimization report for: std::basic_ofstream<char, std::char_traits<char>>::~basic_ofstream()
-=================================================================
-
-Begin optimization report for: std::ostream::operator<<(double)
-=================================================================
-
-Begin optimization report for: std::_Rb_tree<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>, std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, int>, std::_Select1st<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, int>>, std::less<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>, std::allocator<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, int>>>::_Rb_tree()
-=================================================================
-
-Begin optimization report for: std::_Rb_tree<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>, std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, int>, std::_Select1st<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, int>>, std::less<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>, std::allocator<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, int>>>::_Rb_tree_impl<std::less<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>, true>::_Rb_tree_impl()
-=================================================================
-
-Begin optimization report for: std::allocator<std::_Rb_tree_node<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, int>>>::allocator()
-=================================================================
-
-Begin optimization report for: std::_Rb_tree_key_compare<std::less<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>>::_Rb_tree_key_compare()
-=================================================================
-
-Begin optimization report for: std::_Rb_tree_header::_Rb_tree_header()
-=================================================================
-
-Begin optimization report for: std::__new_allocator<std::_Rb_tree_node<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, int>>>::__new_allocator()
-=================================================================
-
-Begin optimization report for: std::_Rb_tree_header::_M_reset()
-=================================================================
-
-Begin optimization report for: __clang_call_terminate
-=================================================================
-
-Begin optimization report for: std::_Rb_tree<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>, std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, int>, std::_Select1st<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, int>>, std::less<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>, std::allocator<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, int>>>::~_Rb_tree()
-=================================================================
-
-Begin optimization report for: std::_Rb_tree<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>, std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, int>, std::_Select1st<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, int>>, std::less<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>, std::allocator<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, int>>>::_M_erase(std::_Rb_tree_node<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, int>>*)
-
-LOOP BEGIN at /usr/lib/gcc/x86_64-redhat-linux/14/../../../../include/c++/14/bits/stl_tree.h (1930, 7)
-    remark #15521: Loop was not vectorized: loop control variable was not identified. Explicitly compute the iteration count before executing the loop or try using canonical loop form from OpenMP specification
-    remark #25478: While loop unrolled by 2
-LOOP END
-=================================================================
-
-Begin optimization report for: std::_Rb_tree<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>, std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, int>, std::_Select1st<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, int>>, std::less<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>, std::allocator<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, int>>>::_M_begin()
-=================================================================
-
-Begin optimization report for: std::__new_allocator<std::_Rb_tree_node<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, int>>>::~__new_allocator()
-=================================================================
-
-Begin optimization report for: std::_Rb_tree<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>, std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, int>, std::_Select1st<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, int>>, std::less<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>, std::allocator<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, int>>>::_S_right(std::_Rb_tree_node_base*)
-=================================================================
-
-Begin optimization report for: std::_Rb_tree<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>, std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, int>, std::_Select1st<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, int>>, std::less<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>, std::allocator<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, int>>>::_S_left(std::_Rb_tree_node_base*)
-=================================================================
-
-Begin optimization report for: std::_Rb_tree<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>, std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, int>, std::_Select1st<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, int>>, std::less<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>, std::allocator<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, int>>>::_M_drop_node(std::_Rb_tree_node<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, int>>*)
-=================================================================
-
-Begin optimization report for: std::_Rb_tree<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>, std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, int>, std::_Select1st<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, int>>, std::less<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>, std::allocator<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, int>>>::_M_destroy_node(std::_Rb_tree_node<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, int>>*)
-=================================================================
-
-Begin optimization report for: std::_Rb_tree<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>, std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, int>, std::_Select1st<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, int>>, std::less<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>, std::allocator<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, int>>>::_M_put_node(std::_Rb_tree_node<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, int>>*)
-=================================================================
-
-Begin optimization report for: void std::allocator_traits<std::allocator<std::_Rb_tree_node<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, int>>>>::destroy<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, int>>(std::allocator<std::_Rb_tree_node<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, int>>>&, std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, int>*)
-=================================================================
-
-Begin optimization report for: std::_Rb_tree<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>, std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, int>, std::_Select1st<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, int>>, std::less<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>, std::allocator<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, int>>>::_M_get_Node_allocator()
-=================================================================
-
-Begin optimization report for: std::_Rb_tree_node<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, int>>::_M_valptr()
-=================================================================
-
-Begin optimization report for: void std::__new_allocator<std::_Rb_tree_node<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, int>>>::destroy<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, int>>(std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, int>*)
-=================================================================
-
-Begin optimization report for: std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, int>::~pair()
-=================================================================
-
-Begin optimization report for: __gnu_cxx::__aligned_membuf<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, int>>::_M_ptr()
-=================================================================
-
-Begin optimization report for: __gnu_cxx::__aligned_membuf<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, int>>::_M_addr()
-=================================================================
-
-Begin optimization report for: std::allocator_traits<std::allocator<std::_Rb_tree_node<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, int>>>>::deallocate(std::allocator<std::_Rb_tree_node<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, int>>>&, std::_Rb_tree_node<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, int>>*, unsigned long)
-=================================================================
-
-Begin optimization report for: std::__new_allocator<std::_Rb_tree_node<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, int>>>::deallocate(std::_Rb_tree_node<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, int>>*, unsigned long)
-=================================================================
-
-Begin optimization report for: std::_Rb_tree<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>, std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, int>, std::_Select1st<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, int>>, std::less<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>, std::allocator<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, int>>>::_M_mbegin() const
-=================================================================
-
-Begin optimization report for: std::_Rb_tree<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>, std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, double>, std::_Select1st<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, double>>, std::less<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>, std::allocator<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, double>>>::_Rb_tree()
-=================================================================
-
-Begin optimization report for: std::_Rb_tree<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>, std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, double>, std::_Select1st<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, double>>, std::less<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>, std::allocator<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, double>>>::_Rb_tree_impl<std::less<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>, true>::_Rb_tree_impl()
-=================================================================
-
-Begin optimization report for: std::allocator<std::_Rb_tree_node<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, double>>>::allocator()
-=================================================================
-
-Begin optimization report for: std::__new_allocator<std::_Rb_tree_node<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, double>>>::__new_allocator()
-=================================================================
-
-Begin optimization report for: std::_Rb_tree<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>, std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, double>, std::_Select1st<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, double>>, std::less<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>, std::allocator<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, double>>>::~_Rb_tree()
-=================================================================
-
-Begin optimization report for: std::_Rb_tree<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>, std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, double>, std::_Select1st<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, double>>, std::less<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>, std::allocator<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, double>>>::_M_erase(std::_Rb_tree_node<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, double>>*)
-
-LOOP BEGIN at /usr/lib/gcc/x86_64-redhat-linux/14/../../../../include/c++/14/bits/stl_tree.h (1930, 7)
-    remark #15521: Loop was not vectorized: loop control variable was not identified. Explicitly compute the iteration count before executing the loop or try using canonical loop form from OpenMP specification
-    remark #25478: While loop unrolled by 2
-LOOP END
-=================================================================
-
-Begin optimization report for: std::_Rb_tree<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>, std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, double>, std::_Select1st<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, double>>, std::less<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>, std::allocator<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, double>>>::_M_begin()
-=================================================================
-
-Begin optimization report for: std::__new_allocator<std::_Rb_tree_node<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, double>>>::~__new_allocator()
-=================================================================
-
-Begin optimization report for: std::_Rb_tree<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>, std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, double>, std::_Select1st<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, double>>, std::less<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>, std::allocator<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, double>>>::_S_right(std::_Rb_tree_node_base*)
-=================================================================
-
-Begin optimization report for: std::_Rb_tree<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>, std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, double>, std::_Select1st<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, double>>, std::less<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>, std::allocator<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, double>>>::_S_left(std::_Rb_tree_node_base*)
-=================================================================
-
-Begin optimization report for: std::_Rb_tree<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>, std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, double>, std::_Select1st<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, double>>, std::less<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>, std::allocator<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, double>>>::_M_drop_node(std::_Rb_tree_node<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, double>>*)
-=================================================================
-
-Begin optimization report for: std::_Rb_tree<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>, std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, double>, std::_Select1st<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, double>>, std::less<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>, std::allocator<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, double>>>::_M_destroy_node(std::_Rb_tree_node<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, double>>*)
-=================================================================
-
-Begin optimization report for: std::_Rb_tree<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>, std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, double>, std::_Select1st<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, double>>, std::less<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>, std::allocator<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, double>>>::_M_put_node(std::_Rb_tree_node<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, double>>*)
-=================================================================
-
-Begin optimization report for: void std::allocator_traits<std::allocator<std::_Rb_tree_node<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, double>>>>::destroy<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, double>>(std::allocator<std::_Rb_tree_node<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, double>>>&, std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, double>*)
-=================================================================
-
-Begin optimization report for: std::_Rb_tree<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>, std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, double>, std::_Select1st<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, double>>, std::less<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>, std::allocator<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, double>>>::_M_get_Node_allocator()
-=================================================================
-
-Begin optimization report for: std::_Rb_tree_node<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, double>>::_M_valptr()
-=================================================================
-
-Begin optimization report for: void std::__new_allocator<std::_Rb_tree_node<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, double>>>::destroy<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, double>>(std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, double>*)
-=================================================================
-
-Begin optimization report for: std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, double>::~pair()
-=================================================================
-
-Begin optimization report for: __gnu_cxx::__aligned_membuf<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, double>>::_M_ptr()
-=================================================================
-
-Begin optimization report for: __gnu_cxx::__aligned_membuf<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, double>>::_M_addr()
-=================================================================
-
-Begin optimization report for: std::allocator_traits<std::allocator<std::_Rb_tree_node<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, double>>>>::deallocate(std::allocator<std::_Rb_tree_node<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, double>>>&, std::_Rb_tree_node<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, double>>*, unsigned long)
-=================================================================
-
-Begin optimization report for: std::__new_allocator<std::_Rb_tree_node<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, double>>>::deallocate(std::_Rb_tree_node<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, double>>*, unsigned long)
-=================================================================
-
-Begin optimization report for: std::_Rb_tree<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>, std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, double>, std::_Select1st<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, double>>, std::less<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>, std::allocator<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, double>>>::_M_mbegin() const
-=================================================================
-
-Begin optimization report for: std::_Rb_tree<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>, std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>, std::_Select1st<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>>, std::less<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>, std::allocator<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>>>::_Rb_tree()
-=================================================================
-
-Begin optimization report for: std::_Rb_tree<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>, std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>, std::_Select1st<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>>, std::less<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>, std::allocator<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>>>::_Rb_tree_impl<std::less<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>, true>::_Rb_tree_impl()
-=================================================================
-
-Begin optimization report for: std::allocator<std::_Rb_tree_node<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>>>::allocator()
-=================================================================
-
-Begin optimization report for: std::__new_allocator<std::_Rb_tree_node<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>>>::__new_allocator()
-=================================================================
-
-Begin optimization report for: std::_Rb_tree<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>, std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>, std::_Select1st<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>>, std::less<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>, std::allocator<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>>>::~_Rb_tree()
-=================================================================
-
-Begin optimization report for: std::_Rb_tree<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>, std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>, std::_Select1st<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>>, std::less<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>, std::allocator<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>>>::_M_erase(std::_Rb_tree_node<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>>*)
-
-LOOP BEGIN at /usr/lib/gcc/x86_64-redhat-linux/14/../../../../include/c++/14/bits/stl_tree.h (1930, 7)
-    remark #15521: Loop was not vectorized: loop control variable was not identified. Explicitly compute the iteration count before executing the loop or try using canonical loop form from OpenMP specification
-LOOP END
-=================================================================
-
-Begin optimization report for: std::_Rb_tree<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>, std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>, std::_Select1st<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>>, std::less<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>, std::allocator<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>>>::_M_begin()
-=================================================================
-
-Begin optimization report for: std::__new_allocator<std::_Rb_tree_node<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>>>::~__new_allocator()
-=================================================================
-
-Begin optimization report for: std::_Rb_tree<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>, std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>, std::_Select1st<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>>, std::less<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>, std::allocator<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>>>::_S_right(std::_Rb_tree_node_base*)
-=================================================================
-
-Begin optimization report for: std::_Rb_tree<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>, std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>, std::_Select1st<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>>, std::less<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>, std::allocator<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>>>::_S_left(std::_Rb_tree_node_base*)
-=================================================================
-
-Begin optimization report for: std::_Rb_tree<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>, std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>, std::_Select1st<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>>, std::less<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>, std::allocator<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>>>::_M_drop_node(std::_Rb_tree_node<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>>*)
-=================================================================
-
-Begin optimization report for: std::_Rb_tree<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>, std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>, std::_Select1st<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>>, std::less<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>, std::allocator<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>>>::_M_destroy_node(std::_Rb_tree_node<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>>*)
-=================================================================
-
-Begin optimization report for: std::_Rb_tree<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>, std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>, std::_Select1st<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>>, std::less<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>, std::allocator<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>>>::_M_put_node(std::_Rb_tree_node<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>>*)
-=================================================================
-
-Begin optimization report for: void std::allocator_traits<std::allocator<std::_Rb_tree_node<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>>>>::destroy<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>>(std::allocator<std::_Rb_tree_node<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>>>&, std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>*)
-=================================================================
-
-Begin optimization report for: std::_Rb_tree<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>, std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>, std::_Select1st<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>>, std::less<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>, std::allocator<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>>>::_M_get_Node_allocator()
-=================================================================
-
-Begin optimization report for: std::_Rb_tree_node<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>>::_M_valptr()
-=================================================================
-
-Begin optimization report for: void std::__new_allocator<std::_Rb_tree_node<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>>>::destroy<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>>(std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>*)
-=================================================================
-
-Begin optimization report for: __gnu_cxx::__aligned_membuf<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>>::_M_ptr()
-=================================================================
-
-Begin optimization report for: __gnu_cxx::__aligned_membuf<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>>::_M_addr()
-=================================================================
-
-Begin optimization report for: std::allocator_traits<std::allocator<std::_Rb_tree_node<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>>>>::deallocate(std::allocator<std::_Rb_tree_node<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>>>&, std::_Rb_tree_node<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>>*, unsigned long)
-=================================================================
-
-Begin optimization report for: std::__new_allocator<std::_Rb_tree_node<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>>>::deallocate(std::_Rb_tree_node<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>>*, unsigned long)
-=================================================================
-
-Begin optimization report for: std::_Rb_tree<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>, std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>, std::_Select1st<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>>, std::less<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>, std::allocator<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>>>::_M_mbegin() const
-=================================================================
-
-Begin optimization report for: std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>::_M_data() const
-=================================================================
-
-Begin optimization report for: std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>::_M_local_data()
-=================================================================
-
-Begin optimization report for: std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>::_Alloc_hider::_Alloc_hider(char*, std::allocator<char>&&)
-=================================================================
-
-Begin optimization report for: std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>::_M_init_local_buf()
-=================================================================
-
-Begin optimization report for: std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>::_M_set_length(unsigned long)
-=================================================================
-
-Begin optimization report for: std::__ptr_traits_ptr_to<char*, char, false>::pointer_to(char&)
-=================================================================
-
-Begin optimization report for: std::__new_allocator<char>::__new_allocator(std::__new_allocator<char> const&)
-=================================================================
-
-Begin optimization report for: std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>::_M_length(unsigned long)
-=================================================================
-
-Begin optimization report for: std::char_traits<char>::assign(char&, char const&)
-=================================================================
-
-Begin optimization report for: std::__new_allocator<char>::__new_allocator()
-=================================================================
-
-Begin optimization report for: std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>::_M_dispose()
-=================================================================
-
-Begin optimization report for: std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>::_Alloc_hider::~_Alloc_hider()
-=================================================================
-
-Begin optimization report for: std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>::_M_is_local() const
-=================================================================
-
-Begin optimization report for: std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>::_M_destroy(unsigned long)
-=================================================================
-
-Begin optimization report for: std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>::_M_local_data() const
-=================================================================
-
-Begin optimization report for: std::__ptr_traits_ptr_to<char const*, char const, false>::pointer_to(char const&)
-=================================================================
-
-Begin optimization report for: std::allocator_traits<std::allocator<char>>::deallocate(std::allocator<char>&, char*, unsigned long)
-=================================================================
-
-Begin optimization report for: std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>::_M_get_allocator()
-=================================================================
-
-Begin optimization report for: std::__new_allocator<char>::deallocate(char*, unsigned long)
-=================================================================
-
-Begin optimization report for: std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>::_Alloc_hider::_Alloc_hider(char*, std::allocator<char> const&)
-=================================================================
-
-Begin optimization report for: std::char_traits<char>::length(char const*)
-=================================================================
-
-Begin optimization report for: void std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>::_M_construct<char const*>(char const*, char const*, std::forward_iterator_tag)
-=================================================================
-
-Begin optimization report for: std::iterator_traits<char const*>::difference_type std::distance<char const*>(char const*, char const*)
-=================================================================
-
-Begin optimization report for: std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>::_M_data(char*)
-=================================================================
-
-Begin optimization report for: std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>::_M_capacity(unsigned long)
-=================================================================
-
-Begin optimization report for: void std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>::_M_construct<char const*>(char const*, char const*, std::forward_iterator_tag)::_Guard::_Guard(std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>*)
-=================================================================
-
-Begin optimization report for: std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>::_S_copy_chars(char*, char const*, char const*)
-=================================================================
-
-Begin optimization report for: void std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>::_M_construct<char const*>(char const*, char const*, std::forward_iterator_tag)::_Guard::~_Guard()
-=================================================================
-
-Begin optimization report for: std::iterator_traits<char const*>::difference_type std::__distance<char const*>(char const*, char const*, std::random_access_iterator_tag)
-=================================================================
-
-Begin optimization report for: std::iterator_traits<char const*>::iterator_category std::__iterator_category<char const*>(char const* const&)
-=================================================================
-
-Begin optimization report for: std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>::_S_copy(char*, char const*, unsigned long)
-=================================================================
-
-Begin optimization report for: std::char_traits<char>::copy(char*, char const*, unsigned long)
-=================================================================
-
-Begin optimization report for: __gnu_cxx::__alloc_traits<std::allocator<char>, char>::_S_select_on_copy(std::allocator<char> const&)
-=================================================================
-
-Begin optimization report for: std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>::_M_get_allocator() const
-=================================================================
-
-Begin optimization report for: void std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>::_M_construct<char*>(char*, char*, std::forward_iterator_tag)
-=================================================================
-
-Begin optimization report for: std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>::length() const
-=================================================================
-
-Begin optimization report for: std::allocator_traits<std::allocator<char>>::select_on_container_copy_construction(std::allocator<char> const&)
-=================================================================
-
-Begin optimization report for: std::allocator<char>::allocator(std::allocator<char> const&)
-=================================================================
-
-Begin optimization report for: std::iterator_traits<char*>::difference_type std::distance<char*>(char*, char*)
-=================================================================
-
-Begin optimization report for: void std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>::_M_construct<char*>(char*, char*, std::forward_iterator_tag)::_Guard::_Guard(std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>*)
-=================================================================
-
-Begin optimization report for: std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>::_S_copy_chars(char*, char*, char*)
-=================================================================
-
-Begin optimization report for: void std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>::_M_construct<char*>(char*, char*, std::forward_iterator_tag)::_Guard::~_Guard()
-=================================================================
-
-Begin optimization report for: std::iterator_traits<char*>::difference_type std::__distance<char*>(char*, char*, std::random_access_iterator_tag)
-=================================================================
-
-Begin optimization report for: std::iterator_traits<char*>::iterator_category std::__iterator_category<char*>(char* const&)
-=================================================================
-
-Begin optimization report for: std::pair<std::_Rb_tree_iterator<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>>, bool> std::_Rb_tree<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>, std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>, std::_Select1st<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>>, std::less<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>, std::allocator<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>>>::_M_insert_unique<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>>(std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>&&)
-
-LOOP BEGIN at /usr/lib/gcc/x86_64-redhat-linux/14/../../../../include/c++/14/bits/stl_tree.h (2113, 7)
-LOOP END
-=================================================================
-
-Begin optimization report for: std::_Rb_tree<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>, std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>, std::_Select1st<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>>, std::less<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>, std::allocator<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>>>::_M_get_insert_unique_pos(std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const&)
-
-LOOP BEGIN at /usr/lib/gcc/x86_64-redhat-linux/14/../../../../include/c++/14/bits/stl_tree.h (2113, 7)
-LOOP END
-=================================================================
-
-Begin optimization report for: std::_Select1st<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>>::operator()(std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>&) const
-=================================================================
-
-Begin optimization report for: std::_Rb_tree<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>, std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>, std::_Select1st<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>>, std::less<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>, std::allocator<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>>>::_Alloc_node::_Alloc_node(std::_Rb_tree<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>, std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>, std::_Select1st<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>>, std::less<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>, std::allocator<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>>>&)
-=================================================================
-
-Begin optimization report for: std::_Rb_tree_iterator<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>> std::_Rb_tree<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>, std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>, std::_Select1st<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>>, std::less<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>, std::allocator<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>>>::_M_insert_<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>, std::_Rb_tree<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>, std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>, std::_Select1st<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>>, std::less<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>, std::allocator<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>>>::_Alloc_node>(std::_Rb_tree_node_base*, std::_Rb_tree_node_base*, std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>&&, std::_Rb_tree<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>, std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>, std::_Select1st<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>>, std::less<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>, std::allocator<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>>>::_Alloc_node&)
-=================================================================
-
-Begin optimization report for: std::pair<std::_Rb_tree_iterator<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>>, bool>::pair<std::_Rb_tree_iterator<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>>, bool, true>(std::_Rb_tree_iterator<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>>&&, bool&&)
-=================================================================
-
-Begin optimization report for: std::_Rb_tree_iterator<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>>::_Rb_tree_iterator(std::_Rb_tree_node_base*)
-=================================================================
-
-Begin optimization report for: std::_Rb_tree<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>, std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>, std::_Select1st<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>>, std::less<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>, std::allocator<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>>>::_M_end()
-=================================================================
-
-Begin optimization report for: std::less<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>::operator()(std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const&, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const&) const
-=================================================================
-
-Begin optimization report for: std::_Rb_tree<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>, std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>, std::_Select1st<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>>, std::less<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>, std::allocator<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>>>::_S_key(std::_Rb_tree_node<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>> const*)
-=================================================================
-
-Begin optimization report for: std::operator==(std::_Rb_tree_iterator<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>> const&, std::_Rb_tree_iterator<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>> const&)
-=================================================================
-
-Begin optimization report for: std::_Rb_tree<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>, std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>, std::_Select1st<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>>, std::less<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>, std::allocator<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>>>::begin()
-=================================================================
-
-Begin optimization report for: std::pair<std::_Rb_tree_node_base*, std::_Rb_tree_node_base*>::pair<std::_Rb_tree_node<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>>*&, std::_Rb_tree_node_base*&, true>(std::_Rb_tree_node<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>>*&, std::_Rb_tree_node_base*&)
-=================================================================
-
-Begin optimization report for: std::_Rb_tree_iterator<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>>::operator--()
-=================================================================
-
-Begin optimization report for: std::_Rb_tree<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>, std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>, std::_Select1st<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>>, std::less<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>, std::allocator<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>>>::_S_key(std::_Rb_tree_node_base const*)
-=================================================================
-
-Begin optimization report for: std::pair<std::_Rb_tree_node_base*, std::_Rb_tree_node_base*>::pair<std::_Rb_tree_node_base*, std::_Rb_tree_node_base*, true>(std::_Rb_tree_node_base* const&, std::_Rb_tree_node_base* const&)
-=================================================================
-
-Begin optimization report for: bool std::operator<<char, std::char_traits<char>, std::allocator<char>>(std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const&, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const&)
-=================================================================
-
-Begin optimization report for: std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>::compare(std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const&) const
-=================================================================
-
-Begin optimization report for: std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>::size() const
-=================================================================
-
-Begin optimization report for: unsigned long const& std::min<unsigned long>(unsigned long const&, unsigned long const&)
-=================================================================
-
-Begin optimization report for: std::char_traits<char>::compare(char const*, char const*, unsigned long)
-=================================================================
-
-Begin optimization report for: std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>::data() const
-=================================================================
-
-Begin optimization report for: std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>::_S_compare(unsigned long, unsigned long)
-=================================================================
-
-Begin optimization report for: std::_Select1st<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>>::operator()(std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>> const&) const
-=================================================================
-
-Begin optimization report for: std::_Rb_tree_node<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>>::_M_valptr() const
-=================================================================
-
-Begin optimization report for: __gnu_cxx::__aligned_membuf<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>>::_M_ptr() const
-=================================================================
-
-Begin optimization report for: __gnu_cxx::__aligned_membuf<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>>::_M_addr() const
-=================================================================
-
-Begin optimization report for: std::_Rb_tree_node<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>>* std::_Rb_tree<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>, std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>, std::_Select1st<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>>, std::less<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>, std::allocator<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>>>::_Alloc_node::operator()<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>>(std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>&&) const
-=================================================================
-
-Begin optimization report for: std::_Rb_tree_node<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>>* std::_Rb_tree<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>, std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>, std::_Select1st<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>>, std::less<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>, std::allocator<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>>>::_M_create_node<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>>(std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>&&)
-=================================================================
-
-Begin optimization report for: std::_Rb_tree<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>, std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>, std::_Select1st<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>>, std::less<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>, std::allocator<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>>>::_M_get_node()
-=================================================================
-
-Begin optimization report for: void std::_Rb_tree<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>, std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>, std::_Select1st<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>>, std::less<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>, std::allocator<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>>>::_M_construct_node<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>>(std::_Rb_tree_node<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>>*, std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>&&)
-=================================================================
-
-Begin optimization report for: std::allocator_traits<std::allocator<std::_Rb_tree_node<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>>>>::allocate(std::allocator<std::_Rb_tree_node<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>>>&, unsigned long)
-=================================================================
-
-Begin optimization report for: std::__new_allocator<std::_Rb_tree_node<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>>>::allocate(unsigned long, void const*)
-=================================================================
-
-Begin optimization report for: std::__new_allocator<std::_Rb_tree_node<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>>>::_M_max_size() const
-=================================================================
-
-Begin optimization report for: void std::allocator_traits<std::allocator<std::_Rb_tree_node<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>>>>::construct<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>, std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>>(std::allocator<std::_Rb_tree_node<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>>>&, std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>*, std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>&&)
-=================================================================
-
-Begin optimization report for: void std::__new_allocator<std::_Rb_tree_node<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>>>::construct<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>, std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>>(std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>*, std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>&&)
-=================================================================
-
-Begin optimization report for: std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>::pair(std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>&&)
-=================================================================
-
-Begin optimization report for: std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>::basic_string(std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>&&)
-=================================================================
-
-Begin optimization report for: std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>::_M_use_local_data()
-=================================================================
-
-Begin optimization report for: std::_Rb_tree<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>, std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>, std::_Select1st<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>>, std::less<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>, std::allocator<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>>>::find(std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const&)
-
-LOOP BEGIN at /usr/lib/gcc/x86_64-redhat-linux/14/../../../../include/c++/14/bits/stl_tree.h (1947, 7)
-LOOP END
-=================================================================
-
-Begin optimization report for: std::_Rb_tree<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>, std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>, std::_Select1st<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>>, std::less<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>, std::allocator<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>>>::_M_lower_bound(std::_Rb_tree_node<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>>*, std::_Rb_tree_node_base*, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const&)
-
-LOOP BEGIN at /usr/lib/gcc/x86_64-redhat-linux/14/../../../../include/c++/14/bits/stl_tree.h (1947, 7)
-LOOP END
-=================================================================
-
-Begin optimization report for: std::_Rb_tree<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>, std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>, std::_Select1st<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>>, std::less<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>, std::allocator<std::pair<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>>>>::end()
-=================================================================
-
-Begin optimization report for: std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>::assign(std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>> const&)
-=================================================================
-
-Begin optimization report for: std::basic_ios<char, std::char_traits<char>>::setstate(std::_Ios_Iostate)
-=================================================================
-
-Begin optimization report for: std::operator|(std::_Ios_Iostate, std::_Ios_Iostate)
-=================================================================
-
-Begin optimization report for: std::basic_ios<char, std::char_traits<char>>::rdstate() const
-=================================================================
-
-Begin optimization report for: std::basic_ostream<char, std::char_traits<char>>& std::flush<char, std::char_traits<char>>(std::basic_ostream<char, std::char_traits<char>>&)
-=================================================================
-
-Begin optimization report for: std::basic_ios<char, std::char_traits<char>>::widen(char) const
-=================================================================
-
-Begin optimization report for: std::ctype<char> const& std::__check_facet<std::ctype<char>>(std::ctype<char> const*)
-=================================================================
-
-Begin optimization report for: std::ctype<char>::widen(char) const
-=================================================================
-
-Begin optimization report for: std::basic_ios<char, std::char_traits<char>>::basic_ios()
-=================================================================
-
-Begin optimization report for: std::basic_istream<char, std::char_traits<char>>::basic_istream()
-=================================================================
-
-Begin optimization report for: std::basic_ifstream<char, std::char_traits<char>>::open(char const*, std::_Ios_Openmode)
-=================================================================
-
-Begin optimization report for: std::basic_filebuf<char, std::char_traits<char>>::~basic_filebuf()
-=================================================================
-
-Begin optimization report for: std::basic_istream<char, std::char_traits<char>>::~basic_istream()
-=================================================================
-
-Begin optimization report for: std::operator|(std::_Ios_Openmode, std::_Ios_Openmode)
-=================================================================
-
-Begin optimization report for: std::basic_streambuf<char, std::char_traits<char>>::~basic_streambuf()
-=================================================================
-
-Begin optimization report for: std::basic_ifstream<char, std::char_traits<char>>::~basic_ifstream()
-=================================================================
-
-Begin optimization report for: virtual thunk to std::basic_ifstream<char, std::char_traits<char>>::~basic_ifstream()
-=================================================================
-
-Begin optimization report for: std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char>>::assign(char const*)
-=================================================================
-
-Begin optimization report for: std::basic_ostream<char, std::char_traits<char>>::basic_ostream()
-=================================================================
-
-Begin optimization report for: std::basic_ostream<char, std::char_traits<char>>::~basic_ostream()
-=================================================================
-
-Begin optimization report for: std::basic_ofstream<char, std::char_traits<char>>::~basic_ofstream()
-=================================================================
-
-Begin optimization report for: virtual thunk to std::basic_ofstream<char, std::char_traits<char>>::~basic_ofstream()
-=================================================================
-
-Begin optimization report for: _GLOBAL__sub_I_ABE.C
-=================================================================
-

AMSS_NCKU_source/TwoPunctures.C

@@ -5,6 +5,7 @@
 #include <cstdio>
 #include <cstdlib>
 #include <string>
+#include <cstring>
 #include <iostream>
 #include <iomanip>
 #include <fstream>
@@ -60,13 +61,110 @@ TwoPunctures::TwoPunctures(double mp, double mm, double b,
   F = dvector(0, ntotal - 1);
   allocate_derivs(&u, ntotal);
   allocate_derivs(&v, ntotal);
+
+  // Allocate workspace buffers for hot-path allocation elimination
+  int N = maximum3(n1, n2, n3);
+  int maxn = maximum2(n1, n2);
+
+  // LineRelax_be workspace (sized for n2)
+  ws_diag_be = new double[n2];
+  ws_e_be = new double[n2 - 1];
+  ws_f_be = new double[n2 - 1];
+  ws_b_be = new double[n2];
+  ws_x_be = new double[n2];
+
+  // LineRelax_al workspace (sized for n1)
+  ws_diag_al = new double[n1];
+  ws_e_al = new double[n1 - 1];
+  ws_f_al = new double[n1 - 1];
+  ws_b_al = new double[n1];
+  ws_x_al = new double[n1];
+
+  // ThomasAlgorithm workspace (sized for max(n1,n2))
+  ws_thomas_y = new double[maxn];
+
+  // JFD_times_dv workspace (sized for nvar)
+  ws_jfd_values = dvector(0, nvar - 1);
+  allocate_derivs(&ws_jfd_dU, nvar);
+  allocate_derivs(&ws_jfd_U, nvar);
+
+  // chebft_Zeros workspace (sized for N+1)
+  ws_cheb_c = dvector(0, N);
+
+  // fourft workspace (sized for N/2+1 each)
+  ws_four_a = dvector(0, N / 2);
+  ws_four_b = dvector(0, N / 2);
+
+  // Derivatives_AB3 workspace
+  ws_deriv_p = dvector(0, N);
+  ws_deriv_dp = dvector(0, N);
+  ws_deriv_d2p = dvector(0, N);
+  ws_deriv_q = dvector(0, N);
+  ws_deriv_dq = dvector(0, N);
+  ws_deriv_r = dvector(0, N);
+  ws_deriv_dr = dvector(0, N);
+  ws_deriv_indx = ivector(0, N);
+
+  // F_of_v workspace
+  ws_fov_sources = new double[n1 * n2 * n3];
+  ws_fov_values = dvector(0, nvar - 1);
+  allocate_derivs(&ws_fov_U, nvar);
+
+  // J_times_dv workspace
+  ws_jtdv_values = dvector(0, nvar - 1);
+  allocate_derivs(&ws_jtdv_dU, nvar);
+  allocate_derivs(&ws_jtdv_U, nvar);
 }
 
 TwoPunctures::~TwoPunctures()
 {
+  int const nvar = 1, n1 = npoints_A, n2 = npoints_B, n3 = npoints_phi;
+  int N = maximum3(n1, n2, n3);
+
   free_dvector(F, 0, ntotal - 1);
   free_derivs(&u, ntotal);
   free_derivs(&v, ntotal);
+
+  // Free workspace buffers
+  delete[] ws_diag_be;
+  delete[] ws_e_be;
+  delete[] ws_f_be;
+  delete[] ws_b_be;
+  delete[] ws_x_be;
+
+  delete[] ws_diag_al;
+  delete[] ws_e_al;
+  delete[] ws_f_al;
+  delete[] ws_b_al;
+  delete[] ws_x_al;
+
+  delete[] ws_thomas_y;
+
+  free_dvector(ws_jfd_values, 0, nvar - 1);
+  free_derivs(&ws_jfd_dU, nvar);
+  free_derivs(&ws_jfd_U, nvar);
+
+  free_dvector(ws_cheb_c, 0, N);
+
+  free_dvector(ws_four_a, 0, N / 2);
+  free_dvector(ws_four_b, 0, N / 2);
+
+  free_dvector(ws_deriv_p, 0, N);
+  free_dvector(ws_deriv_dp, 0, N);
+  free_dvector(ws_deriv_d2p, 0, N);
+  free_dvector(ws_deriv_q, 0, N);
+  free_dvector(ws_deriv_dq, 0, N);
+  free_dvector(ws_deriv_r, 0, N);
+  free_dvector(ws_deriv_dr, 0, N);
+  free_ivector(ws_deriv_indx, 0, N);
+
+  delete[] ws_fov_sources;
+  free_dvector(ws_fov_values, 0, nvar - 1);
+  free_derivs(&ws_fov_U, nvar);
+
+  free_dvector(ws_jtdv_values, 0, nvar - 1);
+  free_derivs(&ws_jtdv_dU, nvar);
+  free_derivs(&ws_jtdv_U, nvar);
 }
 
 void TwoPunctures::Solve()
@@ -655,7 +753,7 @@ void TwoPunctures::chebft_Zeros(double u[], int n, int inv)
   int k, j, isignum;
   double fac, sum, Pion, *c;
 
-  c = dvector(0, n);
+  c = ws_cheb_c;
   Pion = Pi / n;
   if (inv == 0)
   {
@@ -686,7 +784,6 @@ void TwoPunctures::chebft_Zeros(double u[], int n, int inv)
   }
   for (j = 0; j < n; j++)
     u[j] = c[j];
-  free_dvector(c, 0, n);
 }
 
 /* --------------------------------------------------------------------------*/
@@ -774,8 +871,8 @@ void TwoPunctures::fourft(double *u, int N, int inv)
   double x, x1, fac, Pi_fac, *a, *b;
 
   M = N / 2;
-  a = dvector(0, M);
+  a = ws_four_a;
-  b = dvector(1, M); /* Actually: b=vector(1,M-1) but this is problematic if M=1*/
+  b = ws_four_b - 1; /* offset to match dvector(1,M) indexing */
   fac = 1. / M;
   Pi_fac = Pi * fac;
   if (inv == 0)
@@ -824,8 +921,6 @@ void TwoPunctures::fourft(double *u, int N, int inv)
       iy = -iy;
     }
   }
-  free_dvector(a, 0, M);
-  free_dvector(b, 1, M);
 }
 
 /* -----------------------------------------*/
@@ -1118,14 +1213,14 @@ void TwoPunctures::Derivatives_AB3(int nvar, int n1, int n2, int n3, derivs v)
   double *p, *dp, *d2p, *q, *dq, *r, *dr;
 
   N = maximum3(n1, n2, n3);
-  p = dvector(0, N);
+  p = ws_deriv_p;
-  dp = dvector(0, N);
+  dp = ws_deriv_dp;
-  d2p = dvector(0, N);
+  d2p = ws_deriv_d2p;
-  q = dvector(0, N);
+  q = ws_deriv_q;
-  dq = dvector(0, N);
+  dq = ws_deriv_dq;
-  r = dvector(0, N);
+  r = ws_deriv_r;
-  dr = dvector(0, N);
+  dr = ws_deriv_dr;
-  indx = ivector(0, N);
+  indx = ws_deriv_indx;
 
   for (ivar = 0; ivar < nvar; ivar++)
   {
@@ -1208,14 +1303,6 @@ void TwoPunctures::Derivatives_AB3(int nvar, int n1, int n2, int n3, derivs v)
       }
     }
   }
-  free_dvector(p, 0, N);
-  free_dvector(dp, 0, N);
-  free_dvector(d2p, 0, N);
-  free_dvector(q, 0, N);
-  free_dvector(dq, 0, N);
-  free_dvector(r, 0, N);
-  free_dvector(dr, 0, N);
-  free_ivector(indx, 0, N);
 }
 /* --------------------------------------------------------------------------*/
 void TwoPunctures::Newton(int const nvar, int const n1, int const n2, int const n3,
@@ -1284,10 +1371,11 @@ void TwoPunctures::F_of_v(int nvar, int n1, int n2, int n3, derivs v, double *F,
   derivs U;
   double *sources;
 
-  values = dvector(0, nvar - 1);
+  values = ws_fov_values;
-  allocate_derivs(&U, nvar);
+  U = ws_fov_U;
 
-  sources = (double *)calloc(n1 * n2 * n3, sizeof(double));
+  sources = ws_fov_sources;
+  memset(sources, 0, n1 * n2 * n3 * sizeof(double));
   if (0)
   {
     double *s_x, *s_y, *s_z;
@@ -1442,9 +1530,6 @@ void TwoPunctures::F_of_v(int nvar, int n1, int n2, int n3, derivs v, double *F,
   {
     fclose(debugfile);
   }
-  free(sources);
-  free_dvector(values, 0, nvar - 1);
-  free_derivs(&U, nvar);
 }
 /* --------------------------------------------------------------------------*/
 double TwoPunctures::norm_inf(double const *F, int const ntotal)
@@ -1850,11 +1935,12 @@ void TwoPunctures::J_times_dv(int nvar, int n1, int n2, int n3, derivs dv, doubl
 
   Derivatives_AB3(nvar, n1, n2, n3, dv);
 
+  values = ws_jtdv_values;
+  dU = ws_jtdv_dU;
+  U = ws_jtdv_U;
+
   for (i = 0; i < n1; i++)
   {
-    values = dvector(0, nvar - 1);
-    allocate_derivs(&dU, nvar);
-    allocate_derivs(&U, nvar);
     for (j = 0; j < n2; j++)
     {
       for (k = 0; k < n3; k++)
@@ -1908,9 +1994,6 @@ void TwoPunctures::J_times_dv(int nvar, int n1, int n2, int n3, derivs dv, doubl
         }
       }
     }
-    free_dvector(values, 0, nvar - 1);
-    free_derivs(&dU, nvar);
-    free_derivs(&U, nvar);
   }
 }
 /* --------------------------------------------------------------------------*/
@@ -1957,17 +2040,11 @@ void TwoPunctures::LineRelax_be(double *dv,
 {
   int j, m, Ic, Ip, Im, col, ivar;
 
-  double *diag = new double[n2];
+  double *diag = ws_diag_be;
-  double *e = new double[n2 - 1]; /* above diagonal */
+  double *e = ws_e_be;     /* above diagonal */
-  double *f = new double[n2 - 1]; /* below diagonal */
+  double *f = ws_f_be;     /* below diagonal */
-  double *b = new double[n2];     /* rhs */
+  double *b = ws_b_be;     /* rhs */
-  double *x = new double[n2];     /* solution vector */
+  double *x = ws_x_be;     /* solution vector */
-
-  //  gsl_vector *diag = gsl_vector_alloc(n2);
-  //  gsl_vector *e = gsl_vector_alloc(n2-1); /* above diagonal */
-  //  gsl_vector *f = gsl_vector_alloc(n2-1); /* below diagonal */
-  //  gsl_vector *b = gsl_vector_alloc(n2);   /* rhs */
-  //  gsl_vector *x = gsl_vector_alloc(n2);   /* solution vector */
 
   for (ivar = 0; ivar < nvar; ivar++)
   {
@@ -1977,62 +2054,35 @@ void TwoPunctures::LineRelax_be(double *dv,
     }
     diag[n2 - 1] = 0;
 
-    //    gsl_vector_set_zero(diag);
-    //    gsl_vector_set_zero(e);
-    //    gsl_vector_set_zero(f);
     for (j = 0; j < n2; j++)
     {
       Ip = Index(ivar, i, j + 1, k, nvar, n1, n2, n3);
       Ic = Index(ivar, i, j, k, nvar, n1, n2, n3);
       Im = Index(ivar, i, j - 1, k, nvar, n1, n2, n3);
       b[j] = rhs[Ic];
-      //      gsl_vector_set(b,j,rhs[Ic]);
       for (m = 0; m < ncols[Ic]; m++)
       {
         col = cols[Ic][m];
         if (col != Ip && col != Ic && col != Im)
           b[j] -= JFD[Ic][m] * dv[col];
-        //          *gsl_vector_ptr(b, j) -= JFD[Ic][m] * dv[col];
         else
         {
           if (col == Im && j > 0)
             f[j - 1] = JFD[Ic][m];
-          //            gsl_vector_set(f,j-1,JFD[Ic][m]);
           if (col == Ic)
             diag[j] = JFD[Ic][m];
-          //            gsl_vector_set(diag,j,JFD[Ic][m]);
           if (col == Ip && j < n2 - 1)
             e[j] = JFD[Ic][m];
-          //            gsl_vector_set(e,j,JFD[Ic][m]);
         }
       }
     }
-    //          A x = b
-    //          A = ( d_0 e_0  0   0  )
-    //              ( f_0 d_1 e_1  0  )
-    //              (  0  f_1 d_2 e_2 )
-    //              (  0   0  f_2 d_3 )
-    //
     ThomasAlgorithm(n2, f, diag, e, x, b);
-    //    gsl_linalg_solve_tridiag(diag, e, f, b, x);
     for (j = 0; j < n2; j++)
     {
       Ic = Index(ivar, i, j, k, nvar, n1, n2, n3);
       dv[Ic] = x[j];
-      //      dv[Ic] = gsl_vector_get(x, j);
     }
   }
-
-  delete[] diag;
-  delete[] e;
-  delete[] f;
-  delete[] b;
-  delete[] x;
-  //  gsl_vector_free(diag);
-  //  gsl_vector_free(e);
-  //  gsl_vector_free(f);
-  //  gsl_vector_free(b);
-  //  gsl_vector_free(x);
 }
 /* --------------------------------------------------------------------------*/
 void TwoPunctures::JFD_times_dv(int i, int j, int k, int nvar, int n1, int n2,
@@ -2049,8 +2099,8 @@ void TwoPunctures::JFD_times_dv(int i, int j, int k, int nvar, int n1, int n2,
       ha, ga, ga2, hb, gb, gb2, hp, gp, gp2, gagb, gagp, gbgp;
   derivs dU, U;
 
-  allocate_derivs(&dU, nvar);
+  dU = ws_jfd_dU;
-  allocate_derivs(&U, nvar);
+  U = ws_jfd_U;
 
   if (k < 0)
     k = k + n3;
@@ -2168,9 +2218,6 @@ void TwoPunctures::JFD_times_dv(int i, int j, int k, int nvar, int n1, int n2,
   LinEquations(A, B, X, R, x, r, phi, y, z, dU, U, values);
   for (ivar = 0; ivar < nvar; ivar++)
     values[ivar] *= FAC;
-
-  free_derivs(&dU, nvar);
-  free_derivs(&U, nvar);
 }
 #undef FAC
 /*-----------------------------------------------------------*/
@@ -2202,17 +2249,11 @@ void TwoPunctures::LineRelax_al(double *dv,
 {
   int i, m, Ic, Ip, Im, col, ivar;
 
-  double *diag = new double[n1];
+  double *diag = ws_diag_al;
-  double *e = new double[n1 - 1]; /* above diagonal */
+  double *e = ws_e_al;     /* above diagonal */
-  double *f = new double[n1 - 1]; /* below diagonal */
+  double *f = ws_f_al;     /* below diagonal */
-  double *b = new double[n1];     /* rhs */
+  double *b = ws_b_al;     /* rhs */
-  double *x = new double[n1];     /* solution vector */
+  double *x = ws_x_al;     /* solution vector */
-
-  //  gsl_vector *diag = gsl_vector_alloc(n1);
-  //  gsl_vector *e = gsl_vector_alloc(n1-1); /* above diagonal */
-  //  gsl_vector *f = gsl_vector_alloc(n1-1); /* below diagonal */
-  //  gsl_vector *b = gsl_vector_alloc(n1);   /* rhs */
-  //  gsl_vector *x = gsl_vector_alloc(n1);   /* solution vector */
 
   for (ivar = 0; ivar < nvar; ivar++)
   {
@@ -2222,57 +2263,35 @@ void TwoPunctures::LineRelax_al(double *dv,
     }
     diag[n1 - 1] = 0;
 
-    //    gsl_vector_set_zero(diag);
-    //    gsl_vector_set_zero(e);
-    //    gsl_vector_set_zero(f);
     for (i = 0; i < n1; i++)
     {
       Ip = Index(ivar, i + 1, j, k, nvar, n1, n2, n3);
       Ic = Index(ivar, i, j, k, nvar, n1, n2, n3);
       Im = Index(ivar, i - 1, j, k, nvar, n1, n2, n3);
       b[i] = rhs[Ic];
-      //      gsl_vector_set(b,i,rhs[Ic]);
       for (m = 0; m < ncols[Ic]; m++)
       {
         col = cols[Ic][m];
         if (col != Ip && col != Ic && col != Im)
           b[i] -= JFD[Ic][m] * dv[col];
-        //          *gsl_vector_ptr(b, i) -= JFD[Ic][m] * dv[col];
         else
         {
           if (col == Im && i > 0)
             f[i - 1] = JFD[Ic][m];
-          //            gsl_vector_set(f,i-1,JFD[Ic][m]);
           if (col == Ic)
             diag[i] = JFD[Ic][m];
-          //            gsl_vector_set(diag,i,JFD[Ic][m]);
           if (col == Ip && i < n1 - 1)
             e[i] = JFD[Ic][m];
-          //            gsl_vector_set(e,i,JFD[Ic][m]);
         }
       }
     }
     ThomasAlgorithm(n1, f, diag, e, x, b);
-    //    gsl_linalg_solve_tridiag(diag, e, f, b, x);
     for (i = 0; i < n1; i++)
     {
       Ic = Index(ivar, i, j, k, nvar, n1, n2, n3);
       dv[Ic] = x[i];
-      //      dv[Ic] = gsl_vector_get(x, i);
     }
   }
-
-  delete[] diag;
-  delete[] e;
-  delete[] f;
-  delete[] b;
-  delete[] x;
-
-  //  gsl_vector_free(diag);
-  //  gsl_vector_free(e);
-  //  gsl_vector_free(f);
-  //  gsl_vector_free(b);
-  //  gsl_vector_free(x);
 }
 /* -------------------------------------------------------------------------*/
 // a[N], b[N-1], c[N-1], x[N], q[N]
@@ -2284,44 +2303,29 @@ void TwoPunctures::LineRelax_al(double *dv,
 //"Parallel Scientific Computing in C++ and MPI" P361
 void TwoPunctures::ThomasAlgorithm(int N, double *b, double *a, double *c, double *x, double *q)
 {
+  // In-place Thomas algorithm: uses a[] as d workspace, b[] as l workspace.
+  // c[] is already u (above-diagonal). ws_thomas_y is pre-allocated workspace.
   int i;
-  double *l, *u, *d, *y;
+  double *y = ws_thomas_y;
-  l = new double[N - 1];
-  u = new double[N - 1];
-  d = new double[N];
-  y = new double[N];
-
-  /* LU Decomposition */
-  d[0] = a[0];
-  u[0] = c[0];
 
+  /* LU Decomposition (in-place: a becomes d, b becomes l) */
   for (i = 0; i < N - 2; i++)
   {
-    l[i] = b[i] / d[i];
+    b[i] = b[i] / a[i];
-    d[i + 1] = a[i + 1] - l[i] * u[i];
+    a[i + 1] = a[i + 1] - b[i] * c[i];
-    u[i + 1] = c[i + 1];
   }
-
+  b[N - 2] = b[N - 2] / a[N - 2];
-  l[N - 2] = b[N - 2] / d[N - 2];
+  a[N - 1] = a[N - 1] - b[N - 2] * c[N - 2];
-  d[N - 1] = a[N - 1] - l[N - 2] * u[N - 2];
 
   /* Forward Substitution [L][y] = [q] */
   y[0] = q[0];
   for (i = 1; i < N; i++)
-    y[i] = q[i] - l[i - 1] * y[i - 1];
+    y[i] = q[i] - b[i - 1] * y[i - 1];
 
   /* Backward Substitution [U][x] = [y] */
-  x[N - 1] = y[N - 1] / d[N - 1];
+  x[N - 1] = y[N - 1] / a[N - 1];
-
   for (i = N - 2; i >= 0; i--)
-    x[i] = (y[i] - u[i] * x[i + 1]) / d[i];
+    x[i] = (y[i] - c[i] * x[i + 1]) / a[i];
-
-  delete[] l;
-  delete[] u;
-  delete[] d;
-  delete[] y;
-
-  return;
 }
 // --------------------------------------------------------------------------*/
 // Calculates the value of v at an arbitrary position (x,y,z) if the spectral coefficients are know*/*/

AMSS_NCKU_source/TwoPunctures.h

@@ -42,6 +42,33 @@ private:
 
        int ntotal;
 
+       // Pre-allocated workspace buffers for hot-path allocation elimination
+       // LineRelax_be workspace (sized for n2)
+       double *ws_diag_be, *ws_e_be, *ws_f_be, *ws_b_be, *ws_x_be;
+       // LineRelax_al workspace (sized for n1)
+       double *ws_diag_al, *ws_e_al, *ws_f_al, *ws_b_al, *ws_x_al;
+       // ThomasAlgorithm workspace (sized for max(n1,n2))
+       double *ws_thomas_y;
+       // JFD_times_dv workspace (sized for nvar)
+       double *ws_jfd_values;
+       derivs ws_jfd_dU, ws_jfd_U;
+       // chebft_Zeros workspace (sized for max(n1,n2,n3)+1)
+       double *ws_cheb_c;
+       // fourft workspace (sized for max(n1,n2,n3)/2+1 each)
+       double *ws_four_a, *ws_four_b;
+       // Derivatives_AB3 workspace
+       double *ws_deriv_p, *ws_deriv_dp, *ws_deriv_d2p;
+       double *ws_deriv_q, *ws_deriv_dq;
+       double *ws_deriv_r, *ws_deriv_dr;
+       int *ws_deriv_indx;
+       // F_of_v workspace
+       double *ws_fov_sources;
+       double *ws_fov_values;
+       derivs ws_fov_U;
+       // J_times_dv workspace
+       double *ws_jtdv_values;
+       derivs ws_jtdv_dU, ws_jtdv_U;
+
        struct parameters
        {
               int nvar, n1, n2, n3;

AMSS_NCKU_source/bssn_rhs.f90

@@ -24,9 +24,6 @@
                         Gmx_Res, Gmy_Res, Gmz_Res,                             &
                Symmetry,Lev,eps,co)  result(gont)
 ! calculate constraint violation when co=0               
-#if (GAUGE == 6 || GAUGE == 7)
-  use mkl_vml
-#endif
   implicit none
 
 !~~~~~~> Input parameters:
@@ -86,6 +83,10 @@
   real*8, dimension(ex(1),ex(2),ex(3)) :: gupxx,gupxy,gupxz
   real*8, dimension(ex(1),ex(2),ex(3)) :: gupyy,gupyz,gupzz
 
+! shared work arrays (memory pool) to avoid repeated allocation in subroutines
+  real*8, dimension(-1:ex(1),-1:ex(2),-1:ex(3)) :: fh_work2  ! for fderivs_fh/fdderivs_fh (ghost_width=2)
+  real*8, dimension(-2:ex(1),-2:ex(2),-2:ex(3)) :: fh_work3  ! for kodis_fh/lopsided_fh (ghost_width=3)
+
   real*8,dimension(3) ::SSS,AAS,ASA,SAA,ASS,SAS,SSA
   real*8            :: dX, dY, dZ, PI
   real*8, parameter :: ZEO = 0.d0,ONE = 1.D0, TWO = 2.D0, FOUR = 4.D0
@@ -100,18 +101,17 @@
 #endif
 
 #if (GAUGE == 6 || GAUGE == 7)
-  integer :: BHN,i,j,k, idx, total_points
+  integer :: BHN,i,j,k
   real*8, dimension(9) :: Porg
   real*8, dimension(3) :: Mass
-  real*8 :: r1,r2,M,A,w1,w2,C1,C2,denom_r
+  real*8 :: r1,r2,M,A,w1,w2,C1,C2
   real*8, dimension(ex(1),ex(2),ex(3)) :: reta
-  real*8, allocatable :: vml_r1(:), vml_r2(:)
-  real*8, allocatable :: vml_res1(:), vml_res2(:)
 
   call getpbh(BHN,Porg,Mass)
 #endif
 
-!!! sanity check
+!!! sanity check (disabled in production builds for performance)
+#ifdef DEBUG
   dX = sum(chi)+sum(trK)+sum(dxx)+sum(gxy)+sum(gxz)+sum(dyy)+sum(gyz)+sum(dzz) &
       +sum(Axx)+sum(Axy)+sum(Axz)+sum(Ayy)+sum(Ayz)+sum(Azz)                   &
       +sum(Gamx)+sum(Gamy)+sum(Gamz)                                           &
@@ -141,6 +141,7 @@
      gont = 1
      return
   endif
+#endif
 
   PI = dacos(-ONE)
 
@@ -154,22 +155,22 @@
   gyy = dyy + ONE
   gzz = dzz + ONE
 
-  call fderivs(ex,betax,betaxx,betaxy,betaxz,X,Y,Z,ANTI, SYM, SYM,Symmetry,Lev)
+  call fderivs_fh(ex,betax,betaxx,betaxy,betaxz,X,Y,Z,ANTI, SYM, SYM,Symmetry,Lev,fh_work2)
-  call fderivs(ex,betay,betayx,betayy,betayz,X,Y,Z, SYM,ANTI, SYM,Symmetry,Lev)
+  call fderivs_fh(ex,betay,betayx,betayy,betayz,X,Y,Z, SYM,ANTI, SYM,Symmetry,Lev,fh_work2)
-  call fderivs(ex,betaz,betazx,betazy,betazz,X,Y,Z, SYM, SYM,ANTI,Symmetry,Lev)
+  call fderivs_fh(ex,betaz,betazx,betazy,betazz,X,Y,Z, SYM, SYM,ANTI,Symmetry,Lev,fh_work2)
   
   div_beta = betaxx + betayy + betazz
  
-  call fderivs(ex,chi,chix,chiy,chiz,X,Y,Z,SYM,SYM,SYM,symmetry,Lev)
+  call fderivs_fh(ex,chi,chix,chiy,chiz,X,Y,Z,SYM,SYM,SYM,symmetry,Lev,fh_work2)
 
   chi_rhs = F2o3 *chin1*( alpn1 * trK - div_beta ) !rhs for chi
 
-  call fderivs(ex,dxx,gxxx,gxxy,gxxz,X,Y,Z,SYM ,SYM ,SYM ,Symmetry,Lev)
+  call fderivs_fh(ex,dxx,gxxx,gxxy,gxxz,X,Y,Z,SYM ,SYM ,SYM ,Symmetry,Lev,fh_work2)
-  call fderivs(ex,gxy,gxyx,gxyy,gxyz,X,Y,Z,ANTI,ANTI,SYM ,Symmetry,Lev)
+  call fderivs_fh(ex,gxy,gxyx,gxyy,gxyz,X,Y,Z,ANTI,ANTI,SYM ,Symmetry,Lev,fh_work2)
-  call fderivs(ex,gxz,gxzx,gxzy,gxzz,X,Y,Z,ANTI,SYM ,ANTI,Symmetry,Lev)
+  call fderivs_fh(ex,gxz,gxzx,gxzy,gxzz,X,Y,Z,ANTI,SYM ,ANTI,Symmetry,Lev,fh_work2)
-  call fderivs(ex,dyy,gyyx,gyyy,gyyz,X,Y,Z,SYM ,SYM ,SYM ,Symmetry,Lev)
+  call fderivs_fh(ex,dyy,gyyx,gyyy,gyyz,X,Y,Z,SYM ,SYM ,SYM ,Symmetry,Lev,fh_work2)
-  call fderivs(ex,gyz,gyzx,gyzy,gyzz,X,Y,Z,SYM ,ANTI,ANTI,Symmetry,Lev)
+  call fderivs_fh(ex,gyz,gyzx,gyzy,gyzz,X,Y,Z,SYM ,ANTI,ANTI,Symmetry,Lev,fh_work2)
-  call fderivs(ex,dzz,gzzx,gzzy,gzzz,X,Y,Z,SYM ,SYM ,SYM ,Symmetry,Lev)
+  call fderivs_fh(ex,dzz,gzzx,gzzy,gzzz,X,Y,Z,SYM ,SYM ,SYM ,Symmetry,Lev,fh_work2)
 
   gxx_rhs = - TWO * alpn1 * Axx    -  F2o3 * gxx * div_beta          + &
               TWO *(  gxx * betaxx +   gxy * betayx +   gxz * betazx)
@@ -287,8 +288,8 @@
           (gupyy * gupzz       + gupyz * gupyz)* Ayz
 
 ! Right hand side for Gam^i without shift terms...
-  call fderivs(ex,Lap,Lapx,Lapy,Lapz,X,Y,Z,SYM,SYM,SYM,Symmetry,Lev)
+  call fderivs_fh(ex,Lap,Lapx,Lapy,Lapz,X,Y,Z,SYM,SYM,SYM,Symmetry,Lev,fh_work2)
-  call fderivs(ex,trK,Kx,Ky,Kz,X,Y,Z,SYM,SYM,SYM,symmetry,Lev)
+  call fderivs_fh(ex,trK,Kx,Ky,Kz,X,Y,Z,SYM,SYM,SYM,symmetry,Lev,fh_work2)
 
    Gamx_rhs = - TWO * (   Lapx * Rxx +   Lapy * Rxy +   Lapz * Rxz ) + &
         TWO * alpn1 * (                                                &
@@ -317,12 +318,12 @@
                         Gamzxx * Rxx + Gamzyy * Ryy + Gamzzz * Rzz   + &
                 TWO * ( Gamzxy * Rxy + Gamzxz * Rxz + Gamzyz * Ryz ) )
 
-  call fdderivs(ex,betax,gxxx,gxyx,gxzx,gyyx,gyzx,gzzx,&
+  call fdderivs_fh(ex,betax,gxxx,gxyx,gxzx,gyyx,gyzx,gzzx,&
-                X,Y,Z,ANTI,SYM, SYM ,Symmetry,Lev)
+                X,Y,Z,ANTI,SYM, SYM ,Symmetry,Lev,fh_work2)
-  call fdderivs(ex,betay,gxxy,gxyy,gxzy,gyyy,gyzy,gzzy,&
+  call fdderivs_fh(ex,betay,gxxy,gxyy,gxzy,gyyy,gyzy,gzzy,&
-                X,Y,Z,SYM ,ANTI,SYM ,Symmetry,Lev)
+                X,Y,Z,SYM ,ANTI,SYM ,Symmetry,Lev,fh_work2)
-  call fdderivs(ex,betaz,gxxz,gxyz,gxzz,gyyz,gyzz,gzzz,&
+  call fdderivs_fh(ex,betaz,gxxz,gxyz,gxzz,gyyz,gyzz,gzzz,&
-                X,Y,Z,SYM ,SYM, ANTI,Symmetry,Lev)
+                X,Y,Z,SYM ,SYM, ANTI,Symmetry,Lev,fh_work2)
 
   fxx = gxxx + gxyy + gxzz
   fxy = gxyx + gyyy + gyzz
@@ -335,9 +336,9 @@
   Gamza =       gupxx * Gamzxx + gupyy * Gamzyy + gupzz * Gamzzz + &
           TWO*( gupxy * Gamzxy + gupxz * Gamzxz + gupyz * Gamzyz )
 
-  call fderivs(ex,Gamx,Gamxx,Gamxy,Gamxz,X,Y,Z,ANTI,SYM ,SYM ,Symmetry,Lev)
+  call fderivs_fh(ex,Gamx,Gamxx,Gamxy,Gamxz,X,Y,Z,ANTI,SYM ,SYM ,Symmetry,Lev,fh_work2)
-  call fderivs(ex,Gamy,Gamyx,Gamyy,Gamyz,X,Y,Z,SYM ,ANTI,SYM ,Symmetry,Lev)
+  call fderivs_fh(ex,Gamy,Gamyx,Gamyy,Gamyz,X,Y,Z,SYM ,ANTI,SYM ,Symmetry,Lev,fh_work2)
-  call fderivs(ex,Gamz,Gamzx,Gamzy,Gamzz,X,Y,Z,SYM ,SYM ,ANTI,Symmetry,Lev)
+  call fderivs_fh(ex,Gamz,Gamzx,Gamzy,Gamzz,X,Y,Z,SYM ,SYM ,ANTI,Symmetry,Lev,fh_work2)
 
   Gamx_rhs =               Gamx_rhs +  F2o3 *  Gamxa * div_beta        - &
                      Gamxa * betaxx - Gamya * betaxy - Gamza * betaxz  + &
@@ -380,27 +381,27 @@
   gzzz = gxz * Gamxzz + gyz * Gamyzz + gzz * Gamzzz
 
 !compute Ricci tensor for tilted metric
-   call fdderivs(ex,dxx,fxx,fxy,fxz,fyy,fyz,fzz,X,Y,Z,SYM ,SYM ,SYM ,symmetry,Lev)
+   call fdderivs_fh(ex,dxx,fxx,fxy,fxz,fyy,fyz,fzz,X,Y,Z,SYM ,SYM ,SYM ,symmetry,Lev,fh_work2)
    Rxx =   gupxx * fxx + gupyy * fyy + gupzz * fzz + &
          ( gupxy * fxy + gupxz * fxz + gupyz * fyz ) * TWO
 
-   call fdderivs(ex,dyy,fxx,fxy,fxz,fyy,fyz,fzz,X,Y,Z,SYM ,SYM ,SYM ,symmetry,Lev)
+   call fdderivs_fh(ex,dyy,fxx,fxy,fxz,fyy,fyz,fzz,X,Y,Z,SYM ,SYM ,SYM ,symmetry,Lev,fh_work2)
    Ryy =   gupxx * fxx + gupyy * fyy + gupzz * fzz + &
          ( gupxy * fxy + gupxz * fxz + gupyz * fyz ) * TWO
 
-   call fdderivs(ex,dzz,fxx,fxy,fxz,fyy,fyz,fzz,X,Y,Z,SYM ,SYM ,SYM ,symmetry,Lev)
+   call fdderivs_fh(ex,dzz,fxx,fxy,fxz,fyy,fyz,fzz,X,Y,Z,SYM ,SYM ,SYM ,symmetry,Lev,fh_work2)
    Rzz =   gupxx * fxx + gupyy * fyy + gupzz * fzz + &
          ( gupxy * fxy + gupxz * fxz + gupyz * fyz ) * TWO
 
-   call fdderivs(ex,gxy,fxx,fxy,fxz,fyy,fyz,fzz,X,Y,Z,ANTI, ANTI,SYM ,symmetry,Lev)
+   call fdderivs_fh(ex,gxy,fxx,fxy,fxz,fyy,fyz,fzz,X,Y,Z,ANTI, ANTI,SYM ,symmetry,Lev,fh_work2)
    Rxy =   gupxx * fxx + gupyy * fyy + gupzz * fzz + &
          ( gupxy * fxy + gupxz * fxz + gupyz * fyz ) * TWO
 
-   call fdderivs(ex,gxz,fxx,fxy,fxz,fyy,fyz,fzz,X,Y,Z,ANTI ,SYM ,ANTI,symmetry,Lev)
+   call fdderivs_fh(ex,gxz,fxx,fxy,fxz,fyy,fyz,fzz,X,Y,Z,ANTI ,SYM ,ANTI,symmetry,Lev,fh_work2)
    Rxz =   gupxx * fxx + gupyy * fyy + gupzz * fzz + &
          ( gupxy * fxy + gupxz * fxz + gupyz * fyz ) * TWO
 
-   call fdderivs(ex,gyz,fxx,fxy,fxz,fyy,fyz,fzz,X,Y,Z,SYM ,ANTI ,ANTI,symmetry,Lev)
+   call fdderivs_fh(ex,gyz,fxx,fxy,fxz,fyy,fyz,fzz,X,Y,Z,SYM ,ANTI ,ANTI,symmetry,Lev,fh_work2)
    Ryz =   gupxx * fxx + gupyy * fyy + gupzz * fzz + &
          ( gupxy * fxy + gupxz * fxz + gupyz * fyz ) * TWO
 
@@ -605,7 +606,7 @@
             Gamxzz * gxzy + Gamyzz * gyzy + Gamzzz * gzzy  + &
             Gamxyz * gzzx + Gamyyz * gzzy + Gamzyz * gzzz )
 !covariant second derivative of chi respect to tilted metric
-  call fdderivs(ex,chi,fxx,fxy,fxz,fyy,fyz,fzz,X,Y,Z,SYM,SYM,SYM,Symmetry,Lev)
+  call fdderivs_fh(ex,chi,fxx,fxy,fxz,fyy,fyz,fzz,X,Y,Z,SYM,SYM,SYM,Symmetry,Lev,fh_work2)
 
   fxx = fxx - Gamxxx * chix - Gamyxx * chiy - Gamzxx * chiz
   fxy = fxy - Gamxxy * chix - Gamyxy * chiy - Gamzxy * chiz
@@ -631,8 +632,8 @@
   Ryz = Ryz + (fyz - chiy*chiz/chin1/TWO + gyz * f)/chin1/TWO
 
 ! covariant second derivatives of the lapse respect to physical metric
-  call fdderivs(ex,Lap,fxx,fxy,fxz,fyy,fyz,fzz,X,Y,Z, &
+  call fdderivs_fh(ex,Lap,fxx,fxy,fxz,fyy,fyz,fzz,X,Y,Z, &
-                SYM,SYM,SYM,symmetry,Lev)
+                SYM,SYM,SYM,symmetry,Lev,fh_work2)
 
   gxxx = (gupxx * chix + gupxy * chiy + gupxz * chiz)/chin1
   gxxy = (gupxy * chix + gupyy * chiy + gupyz * chiz)/chin1
@@ -815,7 +816,7 @@
   betay_rhs = FF*dtSfy
   betaz_rhs = FF*dtSfz
 
-  call fderivs(ex,chi,dtSfx_rhs,dtSfy_rhs,dtSfz_rhs,X,Y,Z,SYM,SYM,SYM,Symmetry,Lev)
+  call fderivs_fh(ex,chi,dtSfx_rhs,dtSfy_rhs,dtSfz_rhs,X,Y,Z,SYM,SYM,SYM,Symmetry,Lev,fh_work2)
   reta = gupxx * dtSfx_rhs * dtSfx_rhs + gupyy * dtSfy_rhs * dtSfy_rhs + gupzz * dtSfz_rhs * dtSfz_rhs + &
        TWO * (gupxy * dtSfx_rhs * dtSfy_rhs + gupxz * dtSfx_rhs * dtSfz_rhs + gupyz * dtSfy_rhs * dtSfz_rhs)
   reta = 1.31d0/2*dsqrt(reta/chin1)/(1-dsqrt(chin1))**2
@@ -827,7 +828,7 @@
   betay_rhs = FF*dtSfy
   betaz_rhs = FF*dtSfz
 
-  call fderivs(ex,chi,dtSfx_rhs,dtSfy_rhs,dtSfz_rhs,X,Y,Z,SYM,SYM,SYM,Symmetry,Lev)
+  call fderivs_fh(ex,chi,dtSfx_rhs,dtSfy_rhs,dtSfz_rhs,X,Y,Z,SYM,SYM,SYM,Symmetry,Lev,fh_work2)
   reta = gupxx * dtSfx_rhs * dtSfx_rhs + gupyy * dtSfy_rhs * dtSfy_rhs + gupzz * dtSfz_rhs * dtSfz_rhs + &
        TWO * (gupxy * dtSfx_rhs * dtSfy_rhs + gupxz * dtSfx_rhs * dtSfz_rhs + gupyz * dtSfy_rhs * dtSfz_rhs)
   reta = 1.31d0/2*dsqrt(reta/chin1)/(1-chin1)**2
@@ -835,7 +836,7 @@
   dtSfy_rhs = Gamy_rhs - reta*dtSfy
   dtSfz_rhs = Gamz_rhs - reta*dtSfz
 #elif (GAUGE == 4)
-  call fderivs(ex,chi,dtSfx_rhs,dtSfy_rhs,dtSfz_rhs,X,Y,Z,SYM,SYM,SYM,Symmetry,Lev)
+  call fderivs_fh(ex,chi,dtSfx_rhs,dtSfy_rhs,dtSfz_rhs,X,Y,Z,SYM,SYM,SYM,Symmetry,Lev,fh_work2)
   reta = gupxx * dtSfx_rhs * dtSfx_rhs + gupyy * dtSfy_rhs * dtSfy_rhs + gupzz * dtSfz_rhs * dtSfz_rhs + &
        TWO * (gupxy * dtSfx_rhs * dtSfy_rhs + gupxz * dtSfx_rhs * dtSfz_rhs + gupyz * dtSfy_rhs * dtSfz_rhs)
   reta = 1.31d0/2*dsqrt(reta/chin1)/(1-dsqrt(chin1))**2
@@ -847,7 +848,7 @@
   dtSfy_rhs = ZEO
   dtSfz_rhs = ZEO
 #elif (GAUGE == 5)
-  call fderivs(ex,chi,dtSfx_rhs,dtSfy_rhs,dtSfz_rhs,X,Y,Z,SYM,SYM,SYM,Symmetry,Lev)
+  call fderivs_fh(ex,chi,dtSfx_rhs,dtSfy_rhs,dtSfz_rhs,X,Y,Z,SYM,SYM,SYM,Symmetry,Lev,fh_work2)
   reta = gupxx * dtSfx_rhs * dtSfx_rhs + gupyy * dtSfy_rhs * dtSfy_rhs + gupzz * dtSfz_rhs * dtSfz_rhs + &
        TWO * (gupxy * dtSfx_rhs * dtSfy_rhs + gupxz * dtSfx_rhs * dtSfz_rhs + gupyz * dtSfy_rhs * dtSfz_rhs)
   reta = 1.31d0/2*dsqrt(reta/chin1)/(1-chin1)**2
@@ -867,41 +868,17 @@
    C1 = 1.d0/Mass(1) - A
    C2 = 1.d0/Mass(2) - A
 
-   denom_r = ((Porg(1)-Porg(4))**2+(Porg(2)-Porg(5))**2+(Porg(3)-Porg(6))**2)
-   total_points = ex(1)*ex(2)*ex(3)
-   allocate(vml_r1(total_points), vml_r2(total_points))
-   allocate(vml_res1(total_points), vml_res2(total_points))
-
-   idx = 0
    do k=1,ex(3)
    do j=1,ex(2)
    do i=1,ex(1)
-     idx = idx + 1
+     r1 = ((Porg(1)-X(i))**2+(Porg(2)-Y(j))**2+(Porg(3)-Z(k))**2)/ &
-     r1 = ((Porg(1)-X(i))**2+(Porg(2)-Y(j))**2+(Porg(3)-Z(k))**2)/denom_r
+          ((Porg(1)-Porg(4))**2+(Porg(2)-Porg(5))**2+(Porg(3)-Porg(6))**2)
-     r2 = ((Porg(4)-X(i))**2+(Porg(5)-Y(j))**2+(Porg(6)-Z(k))**2)/denom_r
+     r2 = ((Porg(4)-X(i))**2+(Porg(5)-Y(j))**2+(Porg(6)-Z(k))**2)/ &
-     
+          ((Porg(1)-Porg(4))**2+(Porg(2)-Porg(5))**2+(Porg(3)-Porg(6))**2)
-     ! Prepare for 1/(1+w*r) -> vdInv
+     reta(i,j,k) = A + C1/(ONE+w1*r1) + C2/(ONE+w2*r2)
-     vml_r1(idx) = ONE + w1*r1
-     vml_r2(idx) = ONE + w2*r2
     enddo
     enddo
     enddo
-    
-    call vdInv(total_points, vml_r1, vml_res1)
-    call vdInv(total_points, vml_r2, vml_res2)
-
-    idx = 0
-    do k=1,ex(3)
-    do j=1,ex(2)
-    do i=1,ex(1)
-       idx = idx + 1
-       reta(i,j,k) = A + C1*vml_res1(idx) + C2*vml_res2(idx)
-    enddo
-    enddo
-    enddo
-
-    deallocate(vml_r1, vml_r2, vml_res1, vml_res2)
-
   else
     write(*,*) "not support BH_num in Jason's form 1",BHN
   endif
@@ -921,41 +898,17 @@
    C1 = 1.d0/Mass(1) - A
    C2 = 1.d0/Mass(2) - A
 
-   denom_r = ((Porg(1)-Porg(4))**2+(Porg(2)-Porg(5))**2+(Porg(3)-Porg(6))**2)
-   total_points = ex(1)*ex(2)*ex(3)
-   allocate(vml_r1(total_points), vml_r2(total_points))
-   allocate(vml_res1(total_points), vml_res2(total_points))
-
-   idx = 0
    do k=1,ex(3)
    do j=1,ex(2)
    do i=1,ex(1)
-     idx = idx + 1
+     r1 = ((Porg(1)-X(i))**2+(Porg(2)-Y(j))**2+(Porg(3)-Z(k))**2)/ &
-     r1 = ((Porg(1)-X(i))**2+(Porg(2)-Y(j))**2+(Porg(3)-Z(k))**2)/denom_r
+          ((Porg(1)-Porg(4))**2+(Porg(2)-Porg(5))**2+(Porg(3)-Porg(6))**2)
-     r2 = ((Porg(4)-X(i))**2+(Porg(5)-Y(j))**2+(Porg(6)-Z(k))**2)/denom_r
+     r2 = ((Porg(4)-X(i))**2+(Porg(5)-Y(j))**2+(Porg(6)-Z(k))**2)/ &
-     
+          ((Porg(1)-Porg(4))**2+(Porg(2)-Porg(5))**2+(Porg(3)-Porg(6))**2)
-     ! Prepare for dexp -> vdExp
+     reta(i,j,k) = A + C1*dexp(-w1*r1) + C2*dexp(-w2*r2)
-     vml_r1(idx) = -w1*r1
-     vml_r2(idx) = -w2*r2
     enddo
     enddo
     enddo
-
-    call vdExp(total_points, vml_r1, vml_res1)
-    call vdExp(total_points, vml_r2, vml_res2)
-
-    idx = 0
-    do k=1,ex(3)
-    do j=1,ex(2)
-    do i=1,ex(1)
-       idx = idx + 1
-       reta(i,j,k) = A + C1*vml_res1(idx) + C2*vml_res2(idx)
-    enddo
-    enddo
-    enddo
-
-    deallocate(vml_r1, vml_r2, vml_res1, vml_res2)
-
   else
     write(*,*) "not support BH_num in Jason's form 2",BHN
   endif
@@ -998,51 +951,51 @@
 
 !!!!!!!!!advection term part
 
-  call lopsided(ex,X,Y,Z,gxx,gxx_rhs,betax,betay,betaz,Symmetry,SSS)
+  call lopsided_fh(ex,X,Y,Z,gxx,gxx_rhs,betax,betay,betaz,Symmetry,SSS,fh_work3)
-  call lopsided(ex,X,Y,Z,gxy,gxy_rhs,betax,betay,betaz,Symmetry,AAS)
+  call lopsided_fh(ex,X,Y,Z,gxy,gxy_rhs,betax,betay,betaz,Symmetry,AAS,fh_work3)
-  call lopsided(ex,X,Y,Z,gxz,gxz_rhs,betax,betay,betaz,Symmetry,ASA)
+  call lopsided_fh(ex,X,Y,Z,gxz,gxz_rhs,betax,betay,betaz,Symmetry,ASA,fh_work3)
-  call lopsided(ex,X,Y,Z,gyy,gyy_rhs,betax,betay,betaz,Symmetry,SSS)
+  call lopsided_fh(ex,X,Y,Z,gyy,gyy_rhs,betax,betay,betaz,Symmetry,SSS,fh_work3)
-  call lopsided(ex,X,Y,Z,gyz,gyz_rhs,betax,betay,betaz,Symmetry,SAA)
+  call lopsided_fh(ex,X,Y,Z,gyz,gyz_rhs,betax,betay,betaz,Symmetry,SAA,fh_work3)
-  call lopsided(ex,X,Y,Z,gzz,gzz_rhs,betax,betay,betaz,Symmetry,SSS)
+  call lopsided_fh(ex,X,Y,Z,gzz,gzz_rhs,betax,betay,betaz,Symmetry,SSS,fh_work3)
 
-  call lopsided(ex,X,Y,Z,Axx,Axx_rhs,betax,betay,betaz,Symmetry,SSS)
+  call lopsided_fh(ex,X,Y,Z,Axx,Axx_rhs,betax,betay,betaz,Symmetry,SSS,fh_work3)
-  call lopsided(ex,X,Y,Z,Axy,Axy_rhs,betax,betay,betaz,Symmetry,AAS)
+  call lopsided_fh(ex,X,Y,Z,Axy,Axy_rhs,betax,betay,betaz,Symmetry,AAS,fh_work3)
-  call lopsided(ex,X,Y,Z,Axz,Axz_rhs,betax,betay,betaz,Symmetry,ASA)
+  call lopsided_fh(ex,X,Y,Z,Axz,Axz_rhs,betax,betay,betaz,Symmetry,ASA,fh_work3)
-  call lopsided(ex,X,Y,Z,Ayy,Ayy_rhs,betax,betay,betaz,Symmetry,SSS)
+  call lopsided_fh(ex,X,Y,Z,Ayy,Ayy_rhs,betax,betay,betaz,Symmetry,SSS,fh_work3)
-  call lopsided(ex,X,Y,Z,Ayz,Ayz_rhs,betax,betay,betaz,Symmetry,SAA)
+  call lopsided_fh(ex,X,Y,Z,Ayz,Ayz_rhs,betax,betay,betaz,Symmetry,SAA,fh_work3)
-  call lopsided(ex,X,Y,Z,Azz,Azz_rhs,betax,betay,betaz,Symmetry,SSS)
+  call lopsided_fh(ex,X,Y,Z,Azz,Azz_rhs,betax,betay,betaz,Symmetry,SSS,fh_work3)
 
-  call lopsided(ex,X,Y,Z,chi,chi_rhs,betax,betay,betaz,Symmetry,SSS)
+  call lopsided_fh(ex,X,Y,Z,chi,chi_rhs,betax,betay,betaz,Symmetry,SSS,fh_work3)
-  call lopsided(ex,X,Y,Z,trK,trK_rhs,betax,betay,betaz,Symmetry,SSS)
+  call lopsided_fh(ex,X,Y,Z,trK,trK_rhs,betax,betay,betaz,Symmetry,SSS,fh_work3)
 
-  call lopsided(ex,X,Y,Z,Gamx,Gamx_rhs,betax,betay,betaz,Symmetry,ASS)
+  call lopsided_fh(ex,X,Y,Z,Gamx,Gamx_rhs,betax,betay,betaz,Symmetry,ASS,fh_work3)
-  call lopsided(ex,X,Y,Z,Gamy,Gamy_rhs,betax,betay,betaz,Symmetry,SAS)
+  call lopsided_fh(ex,X,Y,Z,Gamy,Gamy_rhs,betax,betay,betaz,Symmetry,SAS,fh_work3)
-  call lopsided(ex,X,Y,Z,Gamz,Gamz_rhs,betax,betay,betaz,Symmetry,SSA)
+  call lopsided_fh(ex,X,Y,Z,Gamz,Gamz_rhs,betax,betay,betaz,Symmetry,SSA,fh_work3)
 !!
-  call lopsided(ex,X,Y,Z,Lap,Lap_rhs,betax,betay,betaz,Symmetry,SSS)
+  call lopsided_fh(ex,X,Y,Z,Lap,Lap_rhs,betax,betay,betaz,Symmetry,SSS,fh_work3)
 
 #if (GAUGE == 0 || GAUGE == 1 || GAUGE == 2 || GAUGE == 3 || GAUGE == 4 || GAUGE == 5 || GAUGE == 6 || GAUGE == 7)
-  call lopsided(ex,X,Y,Z,betax,betax_rhs,betax,betay,betaz,Symmetry,ASS)
+  call lopsided_fh(ex,X,Y,Z,betax,betax_rhs,betax,betay,betaz,Symmetry,ASS,fh_work3)
-  call lopsided(ex,X,Y,Z,betay,betay_rhs,betax,betay,betaz,Symmetry,SAS)
+  call lopsided_fh(ex,X,Y,Z,betay,betay_rhs,betax,betay,betaz,Symmetry,SAS,fh_work3)
-  call lopsided(ex,X,Y,Z,betaz,betaz_rhs,betax,betay,betaz,Symmetry,SSA)
+  call lopsided_fh(ex,X,Y,Z,betaz,betaz_rhs,betax,betay,betaz,Symmetry,SSA,fh_work3)
 #endif
 
 #if (GAUGE == 0 || GAUGE == 2 || GAUGE == 3 || GAUGE == 6 || GAUGE == 7)
-  call lopsided(ex,X,Y,Z,dtSfx,dtSfx_rhs,betax,betay,betaz,Symmetry,ASS)
+  call lopsided_fh(ex,X,Y,Z,dtSfx,dtSfx_rhs,betax,betay,betaz,Symmetry,ASS,fh_work3)
-  call lopsided(ex,X,Y,Z,dtSfy,dtSfy_rhs,betax,betay,betaz,Symmetry,SAS)
+  call lopsided_fh(ex,X,Y,Z,dtSfy,dtSfy_rhs,betax,betay,betaz,Symmetry,SAS,fh_work3)
-  call lopsided(ex,X,Y,Z,dtSfz,dtSfz_rhs,betax,betay,betaz,Symmetry,SSA)
+  call lopsided_fh(ex,X,Y,Z,dtSfz,dtSfz_rhs,betax,betay,betaz,Symmetry,SSA,fh_work3)
 #endif
 
   if(eps>0)then 
 ! usual Kreiss-Oliger dissipation      
-  call kodis(ex,X,Y,Z,chi,chi_rhs,SSS,Symmetry,eps)
+  call kodis_fh(ex,X,Y,Z,chi,chi_rhs,SSS,Symmetry,eps,fh_work3)
-  call kodis(ex,X,Y,Z,trK,trK_rhs,SSS,Symmetry,eps)
+  call kodis_fh(ex,X,Y,Z,trK,trK_rhs,SSS,Symmetry,eps,fh_work3)
-  call kodis(ex,X,Y,Z,dxx,gxx_rhs,SSS,Symmetry,eps)
+  call kodis_fh(ex,X,Y,Z,dxx,gxx_rhs,SSS,Symmetry,eps,fh_work3)
-  call kodis(ex,X,Y,Z,gxy,gxy_rhs,AAS,Symmetry,eps)
+  call kodis_fh(ex,X,Y,Z,gxy,gxy_rhs,AAS,Symmetry,eps,fh_work3)
-  call kodis(ex,X,Y,Z,gxz,gxz_rhs,ASA,Symmetry,eps)
+  call kodis_fh(ex,X,Y,Z,gxz,gxz_rhs,ASA,Symmetry,eps,fh_work3)
-  call kodis(ex,X,Y,Z,dyy,gyy_rhs,SSS,Symmetry,eps)
+  call kodis_fh(ex,X,Y,Z,dyy,gyy_rhs,SSS,Symmetry,eps,fh_work3)
-  call kodis(ex,X,Y,Z,gyz,gyz_rhs,SAA,Symmetry,eps)
+  call kodis_fh(ex,X,Y,Z,gyz,gyz_rhs,SAA,Symmetry,eps,fh_work3)
-  call kodis(ex,X,Y,Z,dzz,gzz_rhs,SSS,Symmetry,eps)
+  call kodis_fh(ex,X,Y,Z,dzz,gzz_rhs,SSS,Symmetry,eps,fh_work3)
 #if 0
 #define i 42
 #define j 40
@@ -1056,7 +1009,7 @@ endif
 #undef k
 !!stop
 #endif
-  call kodis(ex,X,Y,Z,Axx,Axx_rhs,SSS,Symmetry,eps)
+  call kodis_fh(ex,X,Y,Z,Axx,Axx_rhs,SSS,Symmetry,eps,fh_work3)
 #if 0
 #define i 42
 #define j 40
@@ -1070,25 +1023,25 @@ endif
 #undef k
 !!stop
 #endif
-  call kodis(ex,X,Y,Z,Axy,Axy_rhs,AAS,Symmetry,eps)
+  call kodis_fh(ex,X,Y,Z,Axy,Axy_rhs,AAS,Symmetry,eps,fh_work3)
-  call kodis(ex,X,Y,Z,Axz,Axz_rhs,ASA,Symmetry,eps)
+  call kodis_fh(ex,X,Y,Z,Axz,Axz_rhs,ASA,Symmetry,eps,fh_work3)
-  call kodis(ex,X,Y,Z,Ayy,Ayy_rhs,SSS,Symmetry,eps)
+  call kodis_fh(ex,X,Y,Z,Ayy,Ayy_rhs,SSS,Symmetry,eps,fh_work3)
-  call kodis(ex,X,Y,Z,Ayz,Ayz_rhs,SAA,Symmetry,eps)
+  call kodis_fh(ex,X,Y,Z,Ayz,Ayz_rhs,SAA,Symmetry,eps,fh_work3)
-  call kodis(ex,X,Y,Z,Azz,Azz_rhs,SSS,Symmetry,eps)
+  call kodis_fh(ex,X,Y,Z,Azz,Azz_rhs,SSS,Symmetry,eps,fh_work3)
-  call kodis(ex,X,Y,Z,Gamx,Gamx_rhs,ASS,Symmetry,eps)
+  call kodis_fh(ex,X,Y,Z,Gamx,Gamx_rhs,ASS,Symmetry,eps,fh_work3)
-  call kodis(ex,X,Y,Z,Gamy,Gamy_rhs,SAS,Symmetry,eps)
+  call kodis_fh(ex,X,Y,Z,Gamy,Gamy_rhs,SAS,Symmetry,eps,fh_work3)
-  call kodis(ex,X,Y,Z,Gamz,Gamz_rhs,SSA,Symmetry,eps)
+  call kodis_fh(ex,X,Y,Z,Gamz,Gamz_rhs,SSA,Symmetry,eps,fh_work3)
 
 #if 1 
 !! bam does not apply dissipation on gauge variables
-  call kodis(ex,X,Y,Z,Lap,Lap_rhs,SSS,Symmetry,eps)
+  call kodis_fh(ex,X,Y,Z,Lap,Lap_rhs,SSS,Symmetry,eps,fh_work3)
-  call kodis(ex,X,Y,Z,betax,betax_rhs,ASS,Symmetry,eps)
+  call kodis_fh(ex,X,Y,Z,betax,betax_rhs,ASS,Symmetry,eps,fh_work3)
-  call kodis(ex,X,Y,Z,betay,betay_rhs,SAS,Symmetry,eps)
+  call kodis_fh(ex,X,Y,Z,betay,betay_rhs,SAS,Symmetry,eps,fh_work3)
-  call kodis(ex,X,Y,Z,betaz,betaz_rhs,SSA,Symmetry,eps)
+  call kodis_fh(ex,X,Y,Z,betaz,betaz_rhs,SSA,Symmetry,eps,fh_work3)
 #if (GAUGE == 0 || GAUGE == 2 || GAUGE == 3 || GAUGE == 6 || GAUGE == 7)
-  call kodis(ex,X,Y,Z,dtSfx,dtSfx_rhs,ASS,Symmetry,eps)
+  call kodis_fh(ex,X,Y,Z,dtSfx,dtSfx_rhs,ASS,Symmetry,eps,fh_work3)
-  call kodis(ex,X,Y,Z,dtSfy,dtSfy_rhs,SAS,Symmetry,eps)
+  call kodis_fh(ex,X,Y,Z,dtSfy,dtSfy_rhs,SAS,Symmetry,eps,fh_work3)
-  call kodis(ex,X,Y,Z,dtSfz,dtSfz_rhs,SSA,Symmetry,eps)
+  call kodis_fh(ex,X,Y,Z,dtSfz,dtSfz_rhs,SSA,Symmetry,eps,fh_work3)
 #endif
 #endif
 
@@ -1129,12 +1082,12 @@ endif
 
 ! mov_Res_j = gupkj*(-1/chi d_k chi*A_ij + D_k A_ij) - 2/3 d_j trK - 8 PI s_j where D respect to physical metric
 ! store D_i A_jk - 1/chi d_i chi*A_jk in gjk_i
-  call fderivs(ex,Axx,gxxx,gxxy,gxxz,X,Y,Z,SYM ,SYM ,SYM ,Symmetry,0)
+  call fderivs_fh(ex,Axx,gxxx,gxxy,gxxz,X,Y,Z,SYM ,SYM ,SYM ,Symmetry,0,fh_work2)
-  call fderivs(ex,Axy,gxyx,gxyy,gxyz,X,Y,Z,ANTI,ANTI,SYM ,Symmetry,0)
+  call fderivs_fh(ex,Axy,gxyx,gxyy,gxyz,X,Y,Z,ANTI,ANTI,SYM ,Symmetry,0,fh_work2)
-  call fderivs(ex,Axz,gxzx,gxzy,gxzz,X,Y,Z,ANTI,SYM ,ANTI,Symmetry,0)
+  call fderivs_fh(ex,Axz,gxzx,gxzy,gxzz,X,Y,Z,ANTI,SYM ,ANTI,Symmetry,0,fh_work2)
-  call fderivs(ex,Ayy,gyyx,gyyy,gyyz,X,Y,Z,SYM ,SYM ,SYM ,Symmetry,0)
+  call fderivs_fh(ex,Ayy,gyyx,gyyy,gyyz,X,Y,Z,SYM ,SYM ,SYM ,Symmetry,0,fh_work2)
-  call fderivs(ex,Ayz,gyzx,gyzy,gyzz,X,Y,Z,SYM ,ANTI,ANTI,Symmetry,0)
+  call fderivs_fh(ex,Ayz,gyzx,gyzy,gyzz,X,Y,Z,SYM ,ANTI,ANTI,Symmetry,0,fh_work2)
-  call fderivs(ex,Azz,gzzx,gzzy,gzzz,X,Y,Z,SYM ,SYM ,SYM ,Symmetry,0)
+  call fderivs_fh(ex,Azz,gzzx,gzzy,gzzz,X,Y,Z,SYM ,SYM ,SYM ,Symmetry,0,fh_work2)
 
   gxxx = gxxx - (  Gamxxx * Axx + Gamyxx * Axy + Gamzxx * Axz &
                  + Gamxxx * Axx + Gamyxx * Axy + Gamzxx * Axz) - chix*Axx/chin1

AMSS_NCKU_source/diff_new.f90

@@ -69,12 +69,8 @@
   fy = ZEO
   fz = ZEO
 
-! Intel compiler directives for aggressive vectorization
-!DIR$ SIMD
-!DIR$ IVDEP
   do k=1,ex(3)-1
   do j=1,ex(2)-1
-!DIR$ SIMD
   do i=1,ex(1)-1
 ! x direction   
         if(i+1 <= imax .and. i-1 >= imin)then
@@ -962,7 +958,6 @@
   real*8,dimension(-1:ex(1),-1:ex(2),-1:ex(3))   :: fh
   real*8, dimension(3) :: SoA
   integer :: imin,jmin,kmin,imax,jmax,kmax,i,j,k
-  integer :: ib_s, ib_e, jb_s, jb_e, kb_s, kb_e
   real*8 :: d12dx,d12dy,d12dz,d2dx,d2dy,d2dz
   integer, parameter :: NO_SYMM = 0, EQ_SYMM = 1, OCTANT = 2
   real*8,  parameter :: ZEO=0.d0,ONE=1.d0, F60=6.d1
@@ -1002,66 +997,105 @@
   fy = ZEO
   fz = ZEO
 
-  ib_s = max(1, imin + 2)
-  ib_e = min(ex(1)-1, imax - 2)
-  jb_s = max(1, jmin + 2)
-  jb_e = min(ex(2)-1, jmax - 2)
-  kb_s = max(1, kmin + 2)
-  kb_e = min(ex(3)-1, kmax - 2)
-
-! Intel compiler directives for aggressive vectorization
-!DIR$ SIMD
-!DIR$ IVDEP
   do k=1,ex(3)-1
   do j=1,ex(2)-1
+  do i=1,ex(1)-1
+#if 0  
+! x direction   
+        if(i+2 <= imax .and. i-2 >= imin)then
+!
+!              f(i-2) - 8 f(i-1) + 8 f(i+1) - f(i+2)
+!  fx(i) = ---------------------------------------------
+!                             12 dx
+      fx(i,j,k)=d12dx*(fh(i-2,j,k)-EIT*fh(i-1,j,k)+EIT*fh(i+1,j,k)-fh(i+2,j,k))
 
-    ! Check if (j, k) are within the safe 4th-order range
+    elseif(i+1 <= imax .and. i-1 >= imin)then
-    if (k >= kb_s .and. k <= kb_e .and. j >= jb_s .and. j <= jb_e) then
+!
-        
+!              - f(i-1) + f(i+1)
-        ! 1. Left Boundary Peel
+!  fx(i) = --------------------------------
-        do i=1, min(ex(1)-1, ib_s-1)
+!                     2 dx
-           if(i+1 <= imax .and. i-1 >= imin .and. &
-              j+1 <= jmax .and. j-1 >= jmin .and. &
-              k+1 <= kmax .and. k-1 >= kmin) then
       fx(i,j,k)=d2dx*(-fh(i-1,j,k)+fh(i+1,j,k))
-              fy(i,j,k)=d2dy*(-fh(i,j-1,k)+fh(i,j+1,k))
-              fz(i,j,k)=d2dz*(-fh(i,j,k-1)+fh(i,j,k+1))
-           endif
-        enddo
 
-        ! 2. BULK CORE (No branches, full SIMD)
+! set imax and imin 0
-        !DIR$ SIMD
+    endif
-        do i=ib_s, ib_e
+! 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))
+
+    elseif(j+1 <= jmax .and. j-1 >= jmin)then
+
+     fy(i,j,k)=d2dy*(-fh(i,j-1,k)+fh(i,j+1,k))
+
+! set jmax and jmin 0
+    endif
+! 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))
+
+    elseif(k+1 <= kmax .and. k-1 >= kmin)then
+
+      fz(i,j,k)=d2dz*(-fh(i,j,k-1)+fh(i,j,k+1))
+
+! set kmax and kmin 0
+    endif
+#elif 0
+! 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)
+!  fx(i) = ---------------------------------------------
+!                             12 dx
+      fx(i,j,k)=d12dx*(fh(i-2,j,k)-EIT*fh(i-1,j,k)+EIT*fh(i+1,j,k)-fh(i+2,j,k))
+
+    elseif(i+3 <= imax .and. i-1 >= imin)then
+      fx(i,j,k)=d12dx*(-3.d0*fh(i-1,j,k)-1.d1*fh(i,j,k)+1.8d1*fh(i+1,j,k)-6.d0*fh(i+2,j,k)+fh(i+3,j,k))
+    elseif(i+1 <= imax .and. i-3 >= imin)then
+      fx(i,j,k)=d12dx*( 3.d0*fh(i+1,j,k)+1.d1*fh(i,j,k)-1.8d1*fh(i-1,j,k)+6.d0*fh(i-2,j,k)-fh(i-3,j,k))
+! set imax and imin 0
+    endif
+! 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))
+
+    elseif(j+3 <= jmax .and. j-1 >= jmin)then
+      fy(i,j,k)=d12dy*(-3.d0*fh(i,j-1,k)-1.d1*fh(i,j,k)+1.8d1*fh(i,j+1,k)-6.d0*fh(i,j+2,k)+fh(i,j+3,k))
+    elseif(j+1 <= jmax .and. j-3 >= jmin)then
+      fy(i,j,k)=d12dy*( 3.d0*fh(i,j+1,k)+1.d1*fh(i,j,k)-1.8d1*fh(i,j-1,k)+6.d0*fh(i,j-2,k)-fh(i,j-3,k))
+
+! set jmax and jmin 0
+    endif
+! 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))
+
+    elseif(k+3 <= kmax .and. k-1 >= kmin)then
+      fz(i,j,k)=d12dz*(-3.d0*fh(i,j,k-1)-1.d1*fh(i,j,k)+1.8d1*fh(i,j,k+1)-6.d0*fh(i,j,k+2)+fh(i,j,k+3))
+    elseif(k+1 <= kmax .and. k-3 >= kmin)then
+      fz(i,j,k)=d12dz*( 3.d0*fh(i,j,k+1)+1.d1*fh(i,j,k)-1.8d1*fh(i,j,k-1)+6.d0*fh(i,j,k-2)-fh(i,j,k-3))
+
+! set kmax and kmin 0
+    endif
+#else
+! for bam comparison
+   if(i+2 <= imax .and. i-2 >= imin .and. &
+      j+2 <= jmax .and. j-2 >= jmin .and. &
+      k+2 <= kmax .and. k-2 >= kmin) then
       fx(i,j,k)=d12dx*(fh(i-2,j,k)-EIT*fh(i-1,j,k)+EIT*fh(i+1,j,k)-fh(i+2,j,k))
       fy(i,j,k)=d12dy*(fh(i,j-2,k)-EIT*fh(i,j-1,k)+EIT*fh(i,j+1,k)-fh(i,j+2,k))
       fz(i,j,k)=d12dz*(fh(i,j,k-2)-EIT*fh(i,j,k-1)+EIT*fh(i,j,k+1)-fh(i,j,k+2))
-        enddo
+   elseif(i+1 <= imax .and. i-1 >= imin .and. &
-
-        ! 3. Right Boundary Peel
-        do i=max(1, ib_e+1), ex(1)-1
-           if(i+1 <= imax .and. i-1 >= imin .and. &
           j+1 <= jmax .and. j-1 >= jmin .and. &
           k+1 <= kmax .and. k-1 >= kmin) then
       fx(i,j,k)=d2dx*(-fh(i-1,j,k)+fh(i+1,j,k))
       fy(i,j,k)=d2dy*(-fh(i,j-1,k)+fh(i,j+1,k))
       fz(i,j,k)=d2dz*(-fh(i,j,k-1)+fh(i,j,k+1))
    endif
+#endif
   enddo
-        
-    else
-        ! We are in a boundary slab (j or k is near edge)
-        ! Fallback to standard loop with checks (max 2nd order possible here)
-        do i=1,ex(1)-1
-           if(i+1 <= imax .and. i-1 >= imin .and. &
-              j+1 <= jmax .and. j-1 >= jmin .and. &
-              k+1 <= kmax .and. k-1 >= kmin) then
-              fx(i,j,k)=d2dx*(-fh(i-1,j,k)+fh(i+1,j,k))
-              fy(i,j,k)=d2dy*(-fh(i,j-1,k)+fh(i,j+1,k))
-              fz(i,j,k)=d2dz*(-fh(i,j,k-1)+fh(i,j,k+1))
-           endif
-        enddo
-    endif
-
   enddo
   enddo
 
@@ -1069,6 +1103,103 @@
 
   end subroutine fderivs
 !-----------------------------------------------------------------------------
+! fderivs variant: reuses caller-provided fh work array (memory pool)
+!-----------------------------------------------------------------------------
+  subroutine fderivs_fh(ex,f,fx,fy,fz,X,Y,Z,SYM1,SYM2,SYM3, &
+                         symmetry,onoff,fh)
+  implicit none
+
+  integer,                               intent(in ):: ex(1:3),symmetry,onoff
+  real*8,  dimension(ex(1),ex(2),ex(3)), intent(in ):: f
+  real*8,  dimension(ex(1),ex(2),ex(3)), intent(out):: fx,fy,fz
+  real*8,                                intent(in) :: X(ex(1)),Y(ex(2)),Z(ex(3))
+  real*8,                                intent(in ):: SYM1,SYM2,SYM3
+  real*8,dimension(-1:ex(1),-1:ex(2),-1:ex(3)),intent(inout):: fh
+
+  real*8 :: dX,dY,dZ
+  real*8, dimension(3) :: SoA
+  integer :: imin,jmin,kmin,imax,jmax,kmax,i,j,k
+  real*8 :: d12dx,d12dy,d12dz,d2dx,d2dy,d2dz
+  integer, parameter :: NO_SYMM = 0, EQ_SYMM = 1, OCTANT = 2
+  real*8,  parameter :: ZEO=0.d0,ONE=1.d0
+  real*8,  parameter :: TWO=2.d0,EIT=8.d0
+  real*8,  parameter :: F12=1.2d1
+
+  dX = X(2)-X(1)
+  dY = Y(2)-Y(1)
+  dZ = Z(2)-Z(1)
+
+  imax = ex(1)
+  jmax = ex(2)
+  kmax = ex(3)
+
+  imin = 1
+  jmin = 1
+  kmin = 1
+  if(Symmetry > NO_SYMM .and. dabs(Z(1)) < dZ) kmin = -1
+  if(Symmetry > EQ_SYMM .and. dabs(X(1)) < dX) imin = -1
+  if(Symmetry > EQ_SYMM .and. dabs(Y(1)) < dY) jmin = -1
+
+  SoA(1) = SYM1
+  SoA(2) = SYM2
+  SoA(3) = SYM3
+
+  call symmetry_bd(2,ex,f,fh,SoA)
+
+  d12dx = ONE/F12/dX
+  d12dy = ONE/F12/dY
+  d12dz = ONE/F12/dZ
+
+  d2dx = ONE/TWO/dX
+  d2dy = ONE/TWO/dY
+  d2dz = ONE/TWO/dZ
+
+  fx = ZEO
+  fy = ZEO
+  fz = ZEO
+
+  do k=1,ex(3)-1
+  do j=1,ex(2)-1
+  do i=1,ex(1)-1
+#if 0
+   if(i+2 <= imax .and. i-2 >= imin)then
+      fx(i,j,k)=d12dx*(fh(i-2,j,k)-EIT*fh(i-1,j,k)+EIT*fh(i+1,j,k)-fh(i+2,j,k))
+    elseif(i+1 <= imax .and. i-1 >= imin)then
+      fx(i,j,k)=d2dx*(-fh(i-1,j,k)+fh(i+1,j,k))
+    endif
+        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))
+    elseif(j+1 <= jmax .and. j-1 >= jmin)then
+     fy(i,j,k)=d2dy*(-fh(i,j-1,k)+fh(i,j+1,k))
+    endif
+        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))
+    elseif(k+1 <= kmax .and. k-1 >= kmin)then
+      fz(i,j,k)=d2dz*(-fh(i,j,k-1)+fh(i,j,k+1))
+    endif
+#else
+   if(i+2 <= imax .and. i-2 >= imin .and. &
+      j+2 <= jmax .and. j-2 >= jmin .and. &
+      k+2 <= kmax .and. k-2 >= kmin) then
+      fx(i,j,k)=d12dx*(fh(i-2,j,k)-EIT*fh(i-1,j,k)+EIT*fh(i+1,j,k)-fh(i+2,j,k))
+      fy(i,j,k)=d12dy*(fh(i,j-2,k)-EIT*fh(i,j-1,k)+EIT*fh(i,j+1,k)-fh(i,j+2,k))
+      fz(i,j,k)=d12dz*(fh(i,j,k-2)-EIT*fh(i,j,k-1)+EIT*fh(i,j,k+1)-fh(i,j,k+2))
+   elseif(i+1 <= imax .and. i-1 >= imin .and. &
+          j+1 <= jmax .and. j-1 >= jmin .and. &
+          k+1 <= kmax .and. k-1 >= kmin) then
+      fx(i,j,k)=d2dx*(-fh(i-1,j,k)+fh(i+1,j,k))
+      fy(i,j,k)=d2dy*(-fh(i,j-1,k)+fh(i,j+1,k))
+      fz(i,j,k)=d2dz*(-fh(i,j,k-1)+fh(i,j,k+1))
+   endif
+#endif
+  enddo
+  enddo
+  enddo
+
+  return
+
+  end subroutine fderivs_fh
+!-----------------------------------------------------------------------------
 !
 ! single derivatives dx
 !
@@ -1906,6 +2037,162 @@
 
   end subroutine fddyz
 
+!-----------------------------------------------------------------------------
+! fdderivs variant: reuses caller-provided fh work array (memory pool)
+!-----------------------------------------------------------------------------
+  subroutine fdderivs_fh(ex,f,fxx,fxy,fxz,fyy,fyz,fzz,X,Y,Z, &
+                          SYM1,SYM2,SYM3,symmetry,onoff,fh)
+  implicit none
+
+  integer,                             intent(in ):: ex(1:3),symmetry,onoff
+  real*8, dimension(ex(1),ex(2),ex(3)),intent(in ):: f
+  real*8, dimension(ex(1),ex(2),ex(3)),intent(out):: fxx,fxy,fxz,fyy,fyz,fzz
+  real*8,                              intent(in ):: X(ex(1)),Y(ex(2)),Z(ex(3))
+  real*8,                              intent(in ):: SYM1,SYM2,SYM3
+  real*8,dimension(-1:ex(1),-1:ex(2),-1:ex(3)),intent(inout):: fh
+
+  real*8 :: dX,dY,dZ
+  real*8, dimension(3) :: SoA
+  integer :: imin,jmin,kmin,imax,jmax,kmax,i,j,k
+  real*8  :: Sdxdx,Sdydy,Sdzdz,Fdxdx,Fdydy,Fdzdz
+  real*8  :: Sdxdy,Sdxdz,Sdydz,Fdxdy,Fdxdz,Fdydz
+  integer, parameter :: NO_SYMM = 0, EQ_SYMM = 1, OCTANT = 2
+  real*8, parameter :: ZEO=0.d0, ONE=1.d0, TWO=2.d0, F1o4=2.5d-1
+  real*8, parameter :: F8=8.d0, F16=1.6d1, F30=3.d1
+  real*8, parameter :: F1o12=ONE/1.2d1, F1o144=ONE/1.44d2
+
+  dX = X(2)-X(1)
+  dY = Y(2)-Y(1)
+  dZ = Z(2)-Z(1)
+
+  imax = ex(1)
+  jmax = ex(2)
+  kmax = ex(3)
+
+  imin = 1
+  jmin = 1
+  kmin = 1
+  if(Symmetry > NO_SYMM .and. dabs(Z(1)) < dZ) kmin = -1
+  if(Symmetry > EQ_SYMM .and. dabs(X(1)) < dX) imin = -1
+  if(Symmetry > EQ_SYMM .and. dabs(Y(1)) < dY) jmin = -1
+
+  SoA(1) = SYM1
+  SoA(2) = SYM2
+  SoA(3) = SYM3
+
+  call symmetry_bd(2,ex,f,fh,SoA)
+
+  Sdxdx =  ONE /( dX * dX )
+  Sdydy =  ONE /( dY * dY )
+  Sdzdz =  ONE /( dZ * dZ )
+
+  Fdxdx = F1o12 /( dX * dX )
+  Fdydy = F1o12 /( dY * dY )
+  Fdzdz = F1o12 /( dZ * dZ )
+
+  Sdxdy = F1o4 /( dX * dY )
+  Sdxdz = F1o4 /( dX * dZ )
+  Sdydz = F1o4 /( dY * dZ )
+
+  Fdxdy = F1o144 /( dX * dY )
+  Fdxdz = F1o144 /( dX * dZ )
+  Fdydz = F1o144 /( dY * dZ )
+
+  fxx = ZEO
+  fyy = ZEO
+  fzz = ZEO
+  fxy = ZEO
+  fxz = ZEO
+  fyz = ZEO
+
+  do k=1,ex(3)-1
+  do j=1,ex(2)-1
+  do i=1,ex(1)-1
+#if 0
+   if(i+2 <= imax .and. i-2 >= imin)then
+   fxx(i,j,k) = Fdxdx*(-fh(i-2,j,k)+F16*fh(i-1,j,k)-F30*fh(i,j,k) &
+                       -fh(i+2,j,k)+F16*fh(i+1,j,k)              )
+   elseif(i+1 <= imax .and. i-1 >= imin)then
+   fxx(i,j,k) = Sdxdx*(fh(i-1,j,k)-TWO*fh(i,j,k)+fh(i+1,j,k))
+   endif
+        if(j+2 <= jmax .and. j-2 >= jmin)then
+   fyy(i,j,k) = Fdydy*(-fh(i,j-2,k)+F16*fh(i,j-1,k)-F30*fh(i,j,k) &
+                       -fh(i,j+2,k)+F16*fh(i,j+1,k)              )
+   elseif(j+1 <= jmax .and. j-1 >= jmin)then
+   fyy(i,j,k) = Sdydy*(fh(i,j-1,k)-TWO*fh(i,j,k)+fh(i,j+1,k))
+   endif
+        if(k+2 <= kmax .and. k-2 >= kmin)then
+   fzz(i,j,k) = Fdzdz*(-fh(i,j,k-2)+F16*fh(i,j,k-1)-F30*fh(i,j,k) &
+                       -fh(i,j,k+2)+F16*fh(i,j,k+1)              )
+   elseif(k+1 <= kmax .and. k-1 >= kmin)then
+   fzz(i,j,k) = Sdzdz*(fh(i,j,k-1)-TWO*fh(i,j,k)+fh(i,j,k+1))
+   endif
+       if(i+2 <= imax .and. i-2 >= imin .and. j+2 <= jmax .and. j-2 >= jmin)then
+   fxy(i,j,k) = Fdxdy*(     (fh(i-2,j-2,k)-F8*fh(i-1,j-2,k)+F8*fh(i+1,j-2,k)-fh(i+2,j-2,k))  &
+                       -F8 *(fh(i-2,j-1,k)-F8*fh(i-1,j-1,k)+F8*fh(i+1,j-1,k)-fh(i+2,j-1,k))  &
+                       +F8 *(fh(i-2,j+1,k)-F8*fh(i-1,j+1,k)+F8*fh(i+1,j+1,k)-fh(i+2,j+1,k))  &
+                       -    (fh(i-2,j+2,k)-F8*fh(i-1,j+2,k)+F8*fh(i+1,j+2,k)-fh(i+2,j+2,k)))
+   elseif(i+1 <= imax .and. i-1 >= imin .and. j+1 <= jmax .and. j-1 >= jmin)then
+   fxy(i,j,k) = Sdxdy*(fh(i-1,j-1,k)-fh(i+1,j-1,k)-fh(i-1,j+1,k)+fh(i+1,j+1,k))
+   endif
+       if(i+2 <= imax .and. i-2 >= imin .and. k+2 <= kmax .and. k-2 >= kmin)then
+   fxz(i,j,k) = Fdxdz*(     (fh(i-2,j,k-2)-F8*fh(i-1,j,k-2)+F8*fh(i+1,j,k-2)-fh(i+2,j,k-2))  &
+                       -F8 *(fh(i-2,j,k-1)-F8*fh(i-1,j,k-1)+F8*fh(i+1,j,k-1)-fh(i+2,j,k-1))  &
+                       +F8 *(fh(i-2,j,k+1)-F8*fh(i-1,j,k+1)+F8*fh(i+1,j,k+1)-fh(i+2,j,k+1))  &
+                       -    (fh(i-2,j,k+2)-F8*fh(i-1,j,k+2)+F8*fh(i+1,j,k+2)-fh(i+2,j,k+2)))
+   elseif(i+1 <= imax .and. i-1 >= imin .and. k+1 <= kmax .and. k-1 >= kmin)then
+   fxz(i,j,k) = Sdxdz*(fh(i-1,j,k-1)-fh(i+1,j,k-1)-fh(i-1,j,k+1)+fh(i+1,j,k+1))
+   endif
+       if(j+2 <= jmax .and. j-2 >= jmin .and. k+2 <= kmax .and. k-2 >= kmin)then
+   fyz(i,j,k) = Fdydz*(     (fh(i,j-2,k-2)-F8*fh(i,j-1,k-2)+F8*fh(i,j+1,k-2)-fh(i,j+2,k-2))  &
+                       -F8 *(fh(i,j-2,k-1)-F8*fh(i,j-1,k-1)+F8*fh(i,j+1,k-1)-fh(i,j+2,k-1))  &
+                       +F8 *(fh(i,j-2,k+1)-F8*fh(i,j-1,k+1)+F8*fh(i,j+1,k+1)-fh(i,j+2,k+1))  &
+                       -    (fh(i,j-2,k+2)-F8*fh(i,j-1,k+2)+F8*fh(i,j+1,k+2)-fh(i,j+2,k+2)))
+   elseif(j+1 <= jmax .and. j-1 >= jmin .and. k+1 <= kmax .and. k-1 >= kmin)then
+   fyz(i,j,k) = Sdydz*(fh(i,j-1,k-1)-fh(i,j+1,k-1)-fh(i,j-1,k+1)+fh(i,j+1,k+1))
+   endif
+#else
+! for bam comparison
+   if(i+2 <= imax .and. i-2 >= imin .and. &
+      j+2 <= jmax .and. j-2 >= jmin .and. &
+      k+2 <= kmax .and. k-2 >= kmin) then
+   fxx(i,j,k) = Fdxdx*(-fh(i-2,j,k)+F16*fh(i-1,j,k)-F30*fh(i,j,k) &
+                       -fh(i+2,j,k)+F16*fh(i+1,j,k)              )
+   fyy(i,j,k) = Fdydy*(-fh(i,j-2,k)+F16*fh(i,j-1,k)-F30*fh(i,j,k) &
+                       -fh(i,j+2,k)+F16*fh(i,j+1,k)              )
+   fzz(i,j,k) = Fdzdz*(-fh(i,j,k-2)+F16*fh(i,j,k-1)-F30*fh(i,j,k) &
+                       -fh(i,j,k+2)+F16*fh(i,j,k+1)              )
+   fxy(i,j,k) = Fdxdy*(     (fh(i-2,j-2,k)-F8*fh(i-1,j-2,k)+F8*fh(i+1,j-2,k)-fh(i+2,j-2,k))  &
+                       -F8 *(fh(i-2,j-1,k)-F8*fh(i-1,j-1,k)+F8*fh(i+1,j-1,k)-fh(i+2,j-1,k))  &
+                       +F8 *(fh(i-2,j+1,k)-F8*fh(i-1,j+1,k)+F8*fh(i+1,j+1,k)-fh(i+2,j+1,k))  &
+                       -    (fh(i-2,j+2,k)-F8*fh(i-1,j+2,k)+F8*fh(i+1,j+2,k)-fh(i+2,j+2,k)))
+   fxz(i,j,k) = Fdxdz*(     (fh(i-2,j,k-2)-F8*fh(i-1,j,k-2)+F8*fh(i+1,j,k-2)-fh(i+2,j,k-2))  &
+                       -F8 *(fh(i-2,j,k-1)-F8*fh(i-1,j,k-1)+F8*fh(i+1,j,k-1)-fh(i+2,j,k-1))  &
+                       +F8 *(fh(i-2,j,k+1)-F8*fh(i-1,j,k+1)+F8*fh(i+1,j,k+1)-fh(i+2,j,k+1))  &
+                       -    (fh(i-2,j,k+2)-F8*fh(i-1,j,k+2)+F8*fh(i+1,j,k+2)-fh(i+2,j,k+2)))
+   fyz(i,j,k) = Fdydz*(     (fh(i,j-2,k-2)-F8*fh(i,j-1,k-2)+F8*fh(i,j+1,k-2)-fh(i,j+2,k-2))  &
+                       -F8 *(fh(i,j-2,k-1)-F8*fh(i,j-1,k-1)+F8*fh(i,j+1,k-1)-fh(i,j+2,k-1))  &
+                       +F8 *(fh(i,j-2,k+1)-F8*fh(i,j-1,k+1)+F8*fh(i,j+1,k+1)-fh(i,j+2,k+1))  &
+                       -    (fh(i,j-2,k+2)-F8*fh(i,j-1,k+2)+F8*fh(i,j+1,k+2)-fh(i,j+2,k+2)))
+   elseif(i+1 <= imax .and. i-1 >= imin .and. &
+          j+1 <= jmax .and. j-1 >= jmin .and. &
+          k+1 <= kmax .and. k-1 >= kmin) then
+   fxx(i,j,k) = Sdxdx*(fh(i-1,j,k)-TWO*fh(i,j,k)+fh(i+1,j,k))
+   fyy(i,j,k) = Sdydy*(fh(i,j-1,k)-TWO*fh(i,j,k)+fh(i,j+1,k))
+   fzz(i,j,k) = Sdzdz*(fh(i,j,k-1)-TWO*fh(i,j,k)+fh(i,j,k+1))
+   fxy(i,j,k) = Sdxdy*(fh(i-1,j-1,k)-fh(i+1,j-1,k)-fh(i-1,j+1,k)+fh(i+1,j+1,k))
+   fxz(i,j,k) = Sdxdz*(fh(i-1,j,k-1)-fh(i+1,j,k-1)-fh(i-1,j,k+1)+fh(i+1,j,k+1))
+   fyz(i,j,k) = Sdydz*(fh(i,j-1,k-1)-fh(i,j+1,k-1)-fh(i,j-1,k+1)+fh(i,j+1,k+1))
+   endif
+#endif
+   enddo
+   enddo
+   enddo
+
+  return
+
+  end subroutine fdderivs_fh
+
 #elif (ghost_width == 4)
 ! sixth order code
 
@@ -1977,12 +2264,8 @@
   fy = ZEO
   fz = ZEO
 
-! Intel compiler directives for aggressive vectorization
-!DIR$ SIMD
-!DIR$ IVDEP
   do k=1,ex(3)-1
   do j=1,ex(2)-1
-!DIR$ SIMD
   do i=1,ex(1)-1
 ! x direction   
         if(i+3 <= imax .and. i-3 >= imin)then
@@ -3050,12 +3333,8 @@
   fy = ZEO
   fz = ZEO
 
-! Intel compiler directives for aggressive vectorization
-!DIR$ SIMD
-!DIR$ IVDEP
   do k=1,ex(3)-1
   do j=1,ex(2)-1
-!DIR$ SIMD
   do i=1,ex(1)-1
 ! x direction
         if(i+4 <= imax .and. i-4 >= imin)then

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,7 +883,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+1,1:extc(2),1:extc(3))*SoA(1)
@@ -912,7 +908,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 +936,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)
@@ -1118,64 +1112,65 @@ end subroutine d2dump
 ! Lagrangian polynomial interpolation
 !------------------------------------------------------------------------------
 
-  subroutine polint(xa,ya,x,y,dy,ordn)
+  subroutine polint(xa, ya, x, y, dy, ordn)
-
   implicit none
 
-!~~~~~~> Input Parameter:
+  integer, intent(in) :: ordn
-  integer,intent(in) :: ordn
+  real*8, dimension(ordn), intent(in) :: xa, ya
-  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
 
-!~~~~~~> Other parameter:
+  integer :: i, m, ns, n_m
+  real*8, dimension(ordn) :: c, d, ho
+  real*8 :: dif, dift, hp, h, den_val
 
-  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))
 
-  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
+
-    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
+
+      if (den_val == 0.0d0) then
         write(*,*) 'failure in polint for point',x
         write(*,*) 'with input points: ',xa
         stop
-    endif
-    den(1:n-m)=(c(2:n-m+1)-d(1:n-m))/den(1:n-m)
-    d(1:n-m)=ho(1+m:n)*den(1:n-m)
-    c(1:n-m)=ho(1:n-m)*den(1:n-m)
-    if (2*ns < n-m) then
-      dy=c(ns+1)
-    else
-      dy=d(ns)
-      ns=ns-1
       end if
-    y=y+dy
+
+      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
 !------------------------------------------------------------------------------
 !
@@ -1183,35 +1178,37 @@ end subroutine d2dump
 !
 !------------------------------------------------------------------------------
   subroutine polin2(x1a,x2a,ya,x1,x2,y,dy,ordn)
-
   implicit none
 
-!~~~~~~> Input parameters:
   integer,intent(in) :: ordn
   real*8, dimension(1:ordn), intent(in) :: x1a,x2a
   real*8, dimension(1:ordn,1:ordn), intent(in) :: ya
   real*8, intent(in) :: x1,x2
   real*8, intent(out) :: y,dy
 
-!~~~~~~> Other parameters:
+#ifdef POLINT_LEGACY_ORDER
-
   integer  :: i,m
   real*8, dimension(ordn) :: ymtmp
   real*8, dimension(ordn) :: yntmp
 
   m=size(x1a)
-  
   do i=1,m
-
     yntmp=ya(i,:)
     call polint(x2a,yntmp,x2,ymtmp(i),dy,ordn)
-
   end do
-
   call polint(x1a,ymtmp,x1,y,dy,ordn)
+#else
+  integer  :: j
+  real*8, dimension(ordn) :: ymtmp
+  real*8 :: dy_temp
+
+  do j=1,ordn
+    call polint(x1a, ya(:,j), x1, ymtmp(j), dy_temp, ordn)
+  end do
+  call polint(x2a, ymtmp, x2, y, dy, ordn)
+#endif
 
   return
-
   end subroutine polin2
 !------------------------------------------------------------------------------
 !
@@ -1219,18 +1216,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 +1233,33 @@ end subroutine d2dump
 
   m=size(x1a)
   n=size(x2a)
-  
   do i=1,m
    do j=1,n
-
     yqtmp=ya(i,j,:)
     call polint(x3a,yqtmp,x3,yatmp(i,j),dy,ordn)
-
    end do
-
     yntmp=yatmp(i,:)
     call polint(x2a,yntmp,x2,ymtmp(i),dy,ordn)
-
   end do
-
   call polint(x1a,ymtmp,x1,y,dy,ordn)
+#else
+  integer  :: j, k
+  real*8, dimension(ordn,ordn) :: yatmp
+  real*8, dimension(ordn) :: ymtmp
+  real*8 :: dy_temp
+
+  do k=1,ordn
+    do j=1,ordn
+      call polint(x1a, ya(:,j,k), x1, yatmp(j,k), dy_temp, ordn)
+    end do
+  end do
+  do k=1,ordn
+    call polint(x2a, yatmp(:,k), x2, ymtmp(k), dy_temp, ordn)
+  end do
+  call polint(x3a, ymtmp, x3, y, dy, ordn)
+#endif
 
   return
-
   end subroutine polin3
 !--------------------------------------------------------------------------------------
 ! calculate L2norm

AMSS_NCKU_source/kodiss.f90

@@ -215,6 +215,99 @@ integer, parameter :: NO_SYMM=0, OCTANT=2
 
   end subroutine kodis
 
+!-----------------------------------------------------------------------------
+! kodis variant: reuses caller-provided fh work array (memory pool)
+!-----------------------------------------------------------------------------
+subroutine kodis_fh(ex,X,Y,Z,f,f_rhs,SoA,Symmetry,eps,fh)
+
+implicit none
+! argument variables
+integer,intent(in) :: Symmetry
+integer,dimension(3),intent(in)::ex
+real*8, dimension(1:3), intent(in) :: SoA
+double precision,intent(in),dimension(ex(1))::X
+double precision,intent(in),dimension(ex(2))::Y
+double precision,intent(in),dimension(ex(3))::Z
+double precision,intent(in),dimension(ex(1),ex(2),ex(3))::f
+double precision,intent(inout),dimension(ex(1),ex(2),ex(3))::f_rhs
+real*8,intent(in) :: eps
+real*8,dimension(-2:ex(1),-2:ex(2),-2:ex(3)),intent(inout):: fh
+! local variables
+integer :: imin,jmin,kmin,imax,jmax,kmax
+integer :: i,j,k
+real*8  :: dX,dY,dZ
+real*8, parameter :: ONE=1.d0,SIX=6.d0,FIT=1.5d1,TWT=2.d1
+real*8,parameter::cof=6.4d1   ! 2^6
+integer, parameter :: NO_SYMM=0, OCTANT=2
+
+  dX = X(2)-X(1)
+  dY = Y(2)-Y(1)
+  dZ = Z(2)-Z(1)
+
+  imax = ex(1)
+  jmax = ex(2)
+  kmax = ex(3)
+
+  imin = 1
+  jmin = 1
+  kmin = 1
+
+  if(Symmetry > NO_SYMM .and. dabs(Z(1)) < dZ) kmin = -2
+  if(Symmetry == OCTANT .and. dabs(X(1)) < dX) imin = -2
+  if(Symmetry == OCTANT .and. dabs(Y(1)) < dY) jmin = -2
+
+  call symmetry_bd(3,ex,f,fh,SoA)
+
+  do k=1,ex(3)
+  do j=1,ex(2)
+  do i=1,ex(1)
+
+  if(i-3 >= imin .and. i+3 <= imax .and. &
+     j-3 >= jmin .and. j+3 <= jmax .and. &
+     k-3 >= kmin .and. k+3 <= kmax) then
+#if 0
+   f_rhs(i,j,k)       = f_rhs(i,j,k) + eps/dX/cof * (     &
+                              (fh(i-3,j,k)+fh(i+3,j,k)) - &
+                          SIX*(fh(i-2,j,k)+fh(i+2,j,k)) + &
+                          FIT*(fh(i-1,j,k)+fh(i+1,j,k)) - &
+                          TWT* fh(i,j,k)            )
+   f_rhs(i,j,k)       = f_rhs(i,j,k) + eps/dY/cof * (     &
+                              (fh(i,j-3,k)+fh(i,j+3,k)) - &
+                          SIX*(fh(i,j-2,k)+fh(i,j+2,k)) + &
+                          FIT*(fh(i,j-1,k)+fh(i,j+1,k)) - &
+                          TWT* fh(i,j,k)            )
+   f_rhs(i,j,k)       = f_rhs(i,j,k) + eps/dZ/cof * (     &
+                              (fh(i,j,k-3)+fh(i,j,k+3)) - &
+                          SIX*(fh(i,j,k-2)+fh(i,j,k+2)) + &
+                          FIT*(fh(i,j,k-1)+fh(i,j,k+1)) - &
+                          TWT* fh(i,j,k)            )
+#else
+   f_rhs(i,j,k)       = f_rhs(i,j,k) + eps/cof *( (     &
+                              (fh(i-3,j,k)+fh(i+3,j,k)) - &
+                          SIX*(fh(i-2,j,k)+fh(i+2,j,k)) + &
+                          FIT*(fh(i-1,j,k)+fh(i+1,j,k)) - &
+                          TWT* fh(i,j,k)            )/dX + &
+                                                  (     &
+                              (fh(i,j-3,k)+fh(i,j+3,k)) - &
+                          SIX*(fh(i,j-2,k)+fh(i,j+2,k)) + &
+                          FIT*(fh(i,j-1,k)+fh(i,j+1,k)) - &
+                          TWT* fh(i,j,k)            )/dY + &
+                                                  (     &
+                              (fh(i,j,k-3)+fh(i,j,k+3)) - &
+                          SIX*(fh(i,j,k-2)+fh(i,j,k+2)) + &
+                          FIT*(fh(i,j,k-1)+fh(i,j,k+1)) - &
+                          TWT* fh(i,j,k)            )/dZ )
+#endif
+  endif
+
+  enddo
+  enddo
+  enddo
+
+  return
+
+  end subroutine kodis_fh
+
 #elif (ghost_width == 4)
 ! sixth order code
 !------------------------------------------------------------------------------------------------------------------------------

AMSS_NCKU_source/lopsidediff.f90

@@ -487,6 +487,160 @@ subroutine lopsided(ex,X,Y,Z,f,f_rhs,Sfx,Sfy,Sfz,Symmetry,SoA)
 
   end subroutine lopsided
 
+!-----------------------------------------------------------------------------
+! lopsided variant: reuses caller-provided fh work array (memory pool)
+!-----------------------------------------------------------------------------
+subroutine lopsided_fh(ex,X,Y,Z,f,f_rhs,Sfx,Sfy,Sfz,Symmetry,SoA,fh)
+  implicit none
+
+!~~~~~~> Input parameters:
+
+  integer, intent(in)  :: ex(1:3),Symmetry
+  real*8,  intent(in)  :: X(1:ex(1)),Y(1:ex(2)),Z(1:ex(3))
+  real*8,dimension(ex(1),ex(2),ex(3)),intent(in)   :: f,Sfx,Sfy,Sfz
+
+  real*8,dimension(ex(1),ex(2),ex(3)),intent(inout):: f_rhs
+  real*8,dimension(3),intent(in) ::SoA
+  real*8,dimension(-2:ex(1),-2:ex(2),-2:ex(3)),intent(inout):: fh
+
+!~~~~~~> local variables:
+  integer :: imin,jmin,kmin,imax,jmax,kmax,i,j,k
+  real*8 :: dX,dY,dZ
+  real*8 :: d12dx,d12dy,d12dz,d2dx,d2dy,d2dz
+  real*8,  parameter :: ZEO=0.d0,ONE=1.d0, F3=3.d0
+  real*8,  parameter :: TWO=2.d0,F6=6.0d0,F18=1.8d1
+  real*8,  parameter :: F12=1.2d1, F10=1.d1,EIT=8.d0
+  integer, parameter :: NO_SYMM = 0, EQ_SYMM = 1, OCTANT = 2
+
+  dX = X(2)-X(1)
+  dY = Y(2)-Y(1)
+  dZ = Z(2)-Z(1)
+
+  d12dx = ONE/F12/dX
+  d12dy = ONE/F12/dY
+  d12dz = ONE/F12/dZ
+
+  d2dx = ONE/TWO/dX
+  d2dy = ONE/TWO/dY
+  d2dz = ONE/TWO/dZ
+
+  imax = ex(1)
+  jmax = ex(2)
+  kmax = ex(3)
+
+  imin = 1
+  jmin = 1
+  kmin = 1
+  if(Symmetry > NO_SYMM .and. dabs(Z(1)) < dZ) kmin = -2
+  if(Symmetry > EQ_SYMM .and. dabs(X(1)) < dX) imin = -2
+  if(Symmetry > EQ_SYMM .and. dabs(Y(1)) < dY) jmin = -2
+
+  call symmetry_bd(3,ex,f,fh,SoA)
+
+! upper bound set ex-1 only for efficiency,
+! the loop body will set ex 0 also
+  do k=1,ex(3)-1
+  do j=1,ex(2)-1
+  do i=1,ex(1)-1
+#if 0
+!! old code - same as original lopsided
+#else
+!! new code, 2012dec27, based on bam
+! x direction
+    if(Sfx(i,j,k) > ZEO)then
+      if(i+3 <= imax)then
+     f_rhs(i,j,k)=f_rhs(i,j,k)+                                                   &
+                  Sfx(i,j,k)*d12dx*(-F3*fh(i-1,j,k)-F10*fh(i,j,k)+F18*fh(i+1,j,k) &
+                                    -F6*fh(i+2,j,k)+    fh(i+3,j,k))
+     elseif(i+2 <= imax)then
+     f_rhs(i,j,k)=f_rhs(i,j,k)+                                                           &
+                  Sfx(i,j,k)*d12dx*(fh(i-2,j,k)-EIT*fh(i-1,j,k)+EIT*fh(i+1,j,k)-fh(i+2,j,k))
+     elseif(i+1 <= imax)then
+     f_rhs(i,j,k)=f_rhs(i,j,k)-                                                   &
+                  Sfx(i,j,k)*d12dx*(-F3*fh(i+1,j,k)-F10*fh(i,j,k)+F18*fh(i-1,j,k) &
+                                    -F6*fh(i-2,j,k)+    fh(i-3,j,k))
+     endif
+   elseif(Sfx(i,j,k) < ZEO)then
+      if(i-3 >= imin)then
+     f_rhs(i,j,k)=f_rhs(i,j,k)-                                                   &
+                  Sfx(i,j,k)*d12dx*(-F3*fh(i+1,j,k)-F10*fh(i,j,k)+F18*fh(i-1,j,k) &
+                                    -F6*fh(i-2,j,k)+    fh(i-3,j,k))
+     elseif(i-2 >= imin)then
+     f_rhs(i,j,k)=f_rhs(i,j,k)+                                                           &
+                  Sfx(i,j,k)*d12dx*(fh(i-2,j,k)-EIT*fh(i-1,j,k)+EIT*fh(i+1,j,k)-fh(i+2,j,k))
+     elseif(i-1 >= imin)then
+     f_rhs(i,j,k)=f_rhs(i,j,k)+                                                   &
+                  Sfx(i,j,k)*d12dx*(-F3*fh(i-1,j,k)-F10*fh(i,j,k)+F18*fh(i+1,j,k) &
+                                    -F6*fh(i+2,j,k)+    fh(i+3,j,k))
+     endif
+   endif
+
+! y direction
+    if(Sfy(i,j,k) > ZEO)then
+      if(j+3 <= jmax)then
+     f_rhs(i,j,k)=f_rhs(i,j,k)+                                                   &
+                  Sfy(i,j,k)*d12dy*(-F3*fh(i,j-1,k)-F10*fh(i,j,k)+F18*fh(i,j+1,k) &
+                                    -F6*fh(i,j+2,k)+    fh(i,j+3,k))
+     elseif(j+2 <= jmax)then
+     f_rhs(i,j,k)=f_rhs(i,j,k)+                                                           &
+                  Sfy(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))
+     elseif(j+1 <= jmax)then
+     f_rhs(i,j,k)=f_rhs(i,j,k)-                                                   &
+                  Sfy(i,j,k)*d12dy*(-F3*fh(i,j+1,k)-F10*fh(i,j,k)+F18*fh(i,j-1,k) &
+                                    -F6*fh(i,j-2,k)+    fh(i,j-3,k))
+     endif
+   elseif(Sfy(i,j,k) < ZEO)then
+      if(j-3 >= jmin)then
+     f_rhs(i,j,k)=f_rhs(i,j,k)-                                                   &
+                  Sfy(i,j,k)*d12dy*(-F3*fh(i,j+1,k)-F10*fh(i,j,k)+F18*fh(i,j-1,k) &
+                                    -F6*fh(i,j-2,k)+    fh(i,j-3,k))
+     elseif(j-2 >= jmin)then
+     f_rhs(i,j,k)=f_rhs(i,j,k)+                                                           &
+                  Sfy(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))
+     elseif(j-1 >= jmin)then
+     f_rhs(i,j,k)=f_rhs(i,j,k)+                                                   &
+                  Sfy(i,j,k)*d12dy*(-F3*fh(i,j-1,k)-F10*fh(i,j,k)+F18*fh(i,j+1,k) &
+                                    -F6*fh(i,j+2,k)+    fh(i,j+3,k))
+     endif
+   endif
+
+! z direction
+    if(Sfz(i,j,k) > ZEO)then
+      if(k+3 <= kmax)then
+     f_rhs(i,j,k)=f_rhs(i,j,k)+                                                   &
+                  Sfz(i,j,k)*d12dz*(-F3*fh(i,j,k-1)-F10*fh(i,j,k)+F18*fh(i,j,k+1) &
+                                    -F6*fh(i,j,k+2)+    fh(i,j,k+3))
+     elseif(k+2 <= kmax)then
+     f_rhs(i,j,k)=f_rhs(i,j,k)+                                                           &
+                  Sfz(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))
+     elseif(k+1 <= kmax)then
+     f_rhs(i,j,k)=f_rhs(i,j,k)-                                                   &
+                  Sfz(i,j,k)*d12dz*(-F3*fh(i,j,k+1)-F10*fh(i,j,k)+F18*fh(i,j,k-1) &
+                                    -F6*fh(i,j,k-2)+    fh(i,j,k-3))
+     endif
+   elseif(Sfz(i,j,k) < ZEO)then
+      if(k-3 >= kmin)then
+     f_rhs(i,j,k)=f_rhs(i,j,k)-                                                   &
+                  Sfz(i,j,k)*d12dz*(-F3*fh(i,j,k+1)-F10*fh(i,j,k)+F18*fh(i,j,k-1) &
+                                    -F6*fh(i,j,k-2)+    fh(i,j,k-3))
+     elseif(k-2 >= kmin)then
+     f_rhs(i,j,k)=f_rhs(i,j,k)+                                                           &
+                  Sfz(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))
+     elseif(k-1 >= kmin)then
+     f_rhs(i,j,k)=f_rhs(i,j,k)+                                                   &
+                  Sfz(i,j,k)*d12dz*(-F3*fh(i,j,k-1)-F10*fh(i,j,k)+F18*fh(i,j,k+1) &
+                                    -F6*fh(i,j,k+2)+    fh(i,j,k+3))
+     endif
+   endif
+#endif
+  enddo
+  enddo
+  enddo
+
+  return
+
+  end subroutine lopsided_fh
+
 #elif (ghost_width == 4)
 ! sixth order code
 ! Compute advection terms in right hand sides of field equations

AMSS_NCKU_source/makefile.inc

@@ -6,34 +6,20 @@
 ## Intel oneAPI version with oneMKL (Optimized for performance)
 filein  = -I/usr/include/ -I${MKLROOT}/include
 
-## Using multi-threaded MKL for better scalability with MPI
+## Using sequential MKL (OpenMP disabled for better single-threaded performance)
-## This allows MKL functions (FFT, BLAS, LAPACK) to use multiple threads internally
+## Added -lifcore for Intel Fortran runtime and -limf for Intel math library
-## while keeping the application code as pure MPI (no OpenMP pragmas in user code)
+LDLIBS  = -L${MKLROOT}/lib -lmkl_intel_lp64 -lmkl_sequential -lmkl_core -lifcore -limf -lpthread -lm -ldl
-LDLIBS  = -L/usr/lib/x86_64-linux-gnu -L/usr/lib64 -lifcore -limf -lmpi \
-          -L${MKLROOT}/lib -lmkl_intel_lp64 -lmkl_intel_thread -lmkl_core \
-          -liomp5 -lpthread -lm -ldl
 
-## Aggressive optimization flags for maximum performance:
+## Aggressive optimization flags:
-## -O3: Maximum optimization level
+## -O3: Maximum optimization
 ## -xHost: Optimize for the host CPU architecture (Intel/AMD compatible)
-## -fp-model fast=2: Aggressive floating-point optimizations (allows reassociation)
+## -fp-model fast=2: Aggressive floating-point optimizations
 ## -fma: Enable fused multiply-add instructions
-## -qopt-report=5: Generate detailed optimization reports
-## -qopt-report-phase=vec,loop: Report vectorization and loop optimizations
-## -march=native: Use all available CPU instructions
-## -mtune=native: Tune for the specific CPU
-## -funroll-loops: Aggressively unroll loops
-## -fno-alias: Assume no pointer aliasing (safe for Fortran arrays)
-## -qopt-prefetch: Enable aggressive prefetching
 ## Note: OpenMP has been disabled (-qopenmp removed) due to performance issues
-CXXAPPFLAGS  = -O3 -xHost -march=native -mtune=native -fp-model fast=2 -fma \
+CXXAPPFLAGS  = -O3 -xHost -fp-model fast=2 -fma -ipo \
-               -qopt-report=5 -qopt-report-phase=vec,loop \
-               -funroll-loops -fno-alias -qopt-prefetch \
                -Dfortran3 -Dnewc -I${MKLROOT}/include
-f90appflags  = -O3 -xHost -march=native -mtune=native -fp-model fast=2 -fma \
+f90appflags  = -O3 -xHost -fp-model fast=2 -fma -ipo \
-               -qopt-report=5 -qopt-report-phase=vec,loop \
+               -align array64byte -fpp -I${MKLROOT}/include
-               -funroll-loops -fno-alias -qopt-prefetch \
-               -fpp -I${MKLROOT}/include
 f90          = ifx
 f77          = ifx
 CXX          = icpx

makefile_and_run.py

@@ -15,19 +15,13 @@ import subprocess
 ## 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 4-55,60-111"
+NUMACTL_CPU_BIND = "taskset -c 16-47,64-95"
+#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
+BUILD_JOBS = 64
-
-## MKL threading configuration for hybrid MPI+threaded-MKL execution
-## Total cores available: 104 (cores 4-55, 60-111)
-## MPI processes: configured in AMSS_NCKU_Input.py (typically 48)
-## MKL threads per process: 104 / MPI_processes ≈ 2
-## This ensures full utilization of available cores
-MKL_NUM_THREADS = max(1, 104 // input_data.MPI_processes)
 
 
 ##################################################################
@@ -119,20 +113,14 @@ def run_ABE():
     print(                                                      )
     print( " Running the AMSS-NCKU executable file ABE/ABEGPU " ) 
     print(                                                      )
-    print( f" MPI processes: {input_data.MPI_processes}" )
-    print( f" MKL threads per process: {MKL_NUM_THREADS}" )
-    print( f" Total threads: {input_data.MPI_processes * MKL_NUM_THREADS}" )
-    print(                                                      )
 
     ## Define the command to run; cast other values to strings as needed
-    ## Set MKL threading environment variables for optimal performance
-    env_vars = f"export MKL_NUM_THREADS={MKL_NUM_THREADS} && export MKL_DYNAMIC=FALSE && "
     
     if (input_data.GPU_Calculation == "no"):
-        mpi_command         = env_vars + NUMACTL_CPU_BIND + " 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         = env_vars + NUMACTL_CPU_BIND + " 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
@@ -170,9 +158,8 @@ def run_TwoPunctureABE():
     print( " Running the AMSS-NCKU executable file TwoPunctureABE " ) 
     print(                                                          )
     
-    ## Define the command to run with MKL threading configuration
+    ## Define the command to run
-    env_vars = f"export MKL_NUM_THREADS={MKL_NUM_THREADS} && export MKL_DYNAMIC=FALSE && "
+    TwoPuncture_command         = NUMACTL_CPU_BIND + " ./TwoPunctureABE"
-    TwoPuncture_command         = env_vars + NUMACTL_CPU_BIND + " ./TwoPunctureABE"
     TwoPuncture_command_outfile = "TwoPunctureABE_out.log"
 
     ## Execute the command with subprocess.Popen and stream output