64-BitBrainstorm_2026/AMSS-NCKU
8
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Compiler-level and hot-path optimizations for GW150914

Browse Source 
- makefile.inc: add -ipo (interprocedural optimization) and
  -align array64byte (64-byte array alignment for vectorization)
- fmisc.f90: remove redundant funcc=0.d0 zeroing from symmetry_bd,
  symmetry_tbd, symmetry_stbd (~328+ full-array memsets eliminated
  per timestep)
- enforce_algebra.f90: rewrite enforce_ag and enforce_ga as point-wise
  loops, replacing 12 stack-allocated 3D temporary arrays with scalar
  locals for better cache locality

All changes are mathematically equivalent — no algorithmic modifications.

Co-Authored-By: Claude Opus 4.6 <noreply@anthropic.com>
This commit is contained in:
CGH0S7
2026-02-06 17:13:39 +08:00
parent 223ec17a54
commit 6738854a9d
3 changed files with 105 additions and 79 deletions
Download Patch File Download Diff File Expand all files Collapse all files

164
AMSS_NCKU_source/enforce_algebra.f90
View File

@@ -19,48 +19,60 @@
!~~~~~~~> Local variable: !~~~~~~~> 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 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 return
@@ -83,50 +95,70 @@
!~~~~~~~> Local variable: !~~~~~~~> 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 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 return

6
AMSS_NCKU_source/fmisc.f90
View File

@@ -324,7 +324,6 @@ subroutine symmetry_bd(ord,extc,func,funcc,SoA)
integer::i integer::i
funcc = 0.d0
funcc(1:extc(1),1:extc(2),1:extc(3)) = func funcc(1:extc(1),1:extc(2),1:extc(3)) = func
do i=0,ord-1 do i=0,ord-1
funcc(-i,1:extc(2),1:extc(3)) = funcc(i+2,1:extc(2),1:extc(3))*SoA(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 integer::i
funcc = 0.d0
funcc(1:extc(1),1:extc(2),1:extc(3)) = func funcc(1:extc(1),1:extc(2),1:extc(3)) = func
do i=0,ord-1 do i=0,ord-1
funcc(-i,1:extc(2),1:extc(3)) = funcc(i+2,1:extc(2),1:extc(3))*SoA(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 integer::i
funcc = 0.d0
funcc(1:extc(1),1:extc(2),1:extc(3)) = func funcc(1:extc(1),1:extc(2),1:extc(3)) = func
do i=0,ord-1 do i=0,ord-1
funcc(-i,1:extc(2),1:extc(3)) = funcc(i+2,1:extc(2),1:extc(3))*SoA(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 integer::i
funcc = 0.d0
funcc(1:extc(1),1:extc(2),1:extc(3)) = func funcc(1:extc(1),1:extc(2),1:extc(3)) = func
do i=0,ord-1 do i=0,ord-1
funcc(-i,1:extc(2),1:extc(3)) = funcc(i+1,1:extc(2),1:extc(3))*SoA(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 integer::i
funcc = 0.d0
funcc(1:extc(1),1:extc(2),1:extc(3)) = func funcc(1:extc(1),1:extc(2),1:extc(3)) = func
do i=0,ord-1 do i=0,ord-1
funcc(-i,1:extc(2),1:extc(3)) = funcc(i+1,1:extc(2),1:extc(3))*SoA(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 integer::i
funcc = 0.d0
funcc(1:extc(1),1:extc(2),1:extc(3)) = func funcc(1:extc(1),1:extc(2),1:extc(3)) = func
do i=0,ord-1 do i=0,ord-1
funcc(-i,1:extc(2),1:extc(3)) = funcc(i+1,1:extc(2),1:extc(3))*SoA(1) funcc(-i,1:extc(2),1:extc(3)) = funcc(i+1,1:extc(2),1:extc(3))*SoA(1)

6
AMSS_NCKU_source/makefile.inc
View File

@@ -16,10 +16,10 @@ LDLIBS = -L${MKLROOT}/lib -lmkl_intel_lp64 -lmkl_sequential -lmkl_core -lifcore
## -fp-model fast=2: Aggressive floating-point optimizations ## -fp-model fast=2: Aggressive floating-point optimizations
## -fma: Enable fused multiply-add instructions ## -fma: Enable fused multiply-add instructions
## Note: OpenMP has been disabled (-qopenmp removed) due to performance issues ## Note: OpenMP has been disabled (-qopenmp removed) due to performance issues
CXXAPPFLAGS = -O3 -xHost -fp-model fast=2 -fma \ CXXAPPFLAGS = -O3 -xHost -fp-model fast=2 -fma -ipo \
-Dfortran3 -Dnewc -I${MKLROOT}/include -Dfortran3 -Dnewc -I${MKLROOT}/include
f90appflags = -O3 -xHost -fp-model fast=2 -fma \ f90appflags = -O3 -xHost -fp-model fast=2 -fma -ipo \
-fpp -I${MKLROOT}/include -align array64byte -fpp -I${MKLROOT}/include
f90 = ifx f90 = ifx
f77 = ifx f77 = ifx
CXX = icpx CXX = icpx