From 6f62e75245d6d504790cad600186ccaffa12b6a0 Mon Sep 17 00:00:00 2001 From: CGH0S7 <776459475@qq.com> Date: Sun, 17 May 2026 00:35:17 +0800 Subject: [PATCH] Add conservative fmisc safe mode --- AMSS_NCKU_source/fmisc.f90 | 175 +++++++++++++++++++++++++--------- AMSS_NCKU_source/makefile | 9 +- AMSS_NCKU_source/makefile.inc | 5 + 3 files changed, 140 insertions(+), 49 deletions(-) diff --git a/AMSS_NCKU_source/fmisc.f90 b/AMSS_NCKU_source/fmisc.f90 index e55f3de..851b79a 100644 --- a/AMSS_NCKU_source/fmisc.f90 +++ b/AMSS_NCKU_source/fmisc.f90 @@ -312,8 +312,8 @@ end function decide3d !--------------------------------------------------------------------------------------- -subroutine symmetry_bd(ord,extc,func,funcc,SoA) - implicit none +subroutine symmetry_bd(ord,extc,func,funcc,SoA) + implicit none !~~~~~~> input arguments integer,intent(in) :: ord @@ -322,9 +322,12 @@ subroutine symmetry_bd(ord,extc,func,funcc,SoA) real*8, dimension(-ord+1:extc(1),-ord+1:extc(2),-ord+1:extc(3)),intent(out):: funcc real*8, dimension(1:3), intent(in) :: SoA - integer::i - - funcc(1:extc(1),1:extc(2),1:extc(3)) = func + integer::i + +#if USE_FMISC_SAFE_MODE + funcc = 0.d0 +#endif + 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) enddo @@ -337,8 +340,8 @@ subroutine symmetry_bd(ord,extc,func,funcc,SoA) end subroutine symmetry_bd -subroutine symmetry_tbd(ord,extc,func,funcc,SoA) - implicit none +subroutine symmetry_tbd(ord,extc,func,funcc,SoA) + implicit none !~~~~~~> input arguments integer,intent(in) :: ord @@ -347,9 +350,12 @@ subroutine symmetry_tbd(ord,extc,func,funcc,SoA) real*8, dimension(-ord+1:extc(1)+ord,-ord+1:extc(2)+ord,-ord+1:extc(3)+ord),intent(out):: funcc real*8, dimension(1:3), intent(in) :: SoA - integer::i - - funcc(1:extc(1),1:extc(2),1:extc(3)) = func + integer::i + +#if USE_FMISC_SAFE_MODE + funcc = 0.d0 +#endif + 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) funcc(extc(1)+1+i,1:extc(2),1:extc(3)) = funcc(extc(1)-1-i,1:extc(2),1:extc(3))*SoA(1) @@ -365,8 +371,8 @@ subroutine symmetry_tbd(ord,extc,func,funcc,SoA) end subroutine symmetry_tbd -subroutine symmetry_stbd(ord,extc,func,funcc,SoA) - implicit none +subroutine symmetry_stbd(ord,extc,func,funcc,SoA) + implicit none !~~~~~~> input arguments integer,intent(in) :: ord @@ -375,9 +381,12 @@ subroutine symmetry_stbd(ord,extc,func,funcc,SoA) real*8, dimension(-ord+1:extc(1)+ord,-ord+1:extc(2)+ord,extc(3)),intent(out):: funcc real*8, dimension(2), intent(in) :: SoA - integer::i - - funcc(1:extc(1),1:extc(2),1:extc(3)) = func + integer::i + +#if USE_FMISC_SAFE_MODE + funcc = 0.d0 +#endif + 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) funcc(extc(1)+1+i,1:extc(2),1:extc(3)) = funcc(extc(1)-1-i,1:extc(2),1:extc(3))*SoA(1) @@ -1119,6 +1128,10 @@ end subroutine d2dump #define POLINT6_USE_BARYCENTRIC 1 #endif +#ifndef USE_FMISC_SAFE_MODE +#define USE_FMISC_SAFE_MODE 0 +#endif + !DIR$ ATTRIBUTES FORCEINLINE :: polint6_neville subroutine polint6_neville(xa, ya, x, y, dy) implicit none @@ -1271,7 +1284,9 @@ end subroutine d2dump real*8 :: dif, dift, hp, h, den_val if (ordn == 6) then -#if POLINT6_USE_BARYCENTRIC +#if USE_FMISC_SAFE_MODE + call polint6_neville(xa, ya, x, y, dy) +#elif POLINT6_USE_BARYCENTRIC call polint6_barycentric(xa, ya, x, y, dy) #else call polint6_neville(xa, ya, x, y, dy) @@ -1376,10 +1391,10 @@ end subroutine d2dump real*8, intent(in) :: x1,x2 real*8, intent(out) :: y,dy -#ifdef POLINT_LEGACY_ORDER - integer :: i,m - real*8, dimension(ordn) :: ymtmp - real*8, dimension(ordn) :: yntmp +#if USE_FMISC_SAFE_MODE || defined(POLINT_LEGACY_ORDER) + integer :: i,m + real*8, dimension(ordn) :: ymtmp + real*8, dimension(ordn) :: yntmp m=size(x1a) do i=1,m @@ -1414,7 +1429,7 @@ end subroutine d2dump real*8, intent(in) :: x1,x2,x3 real*8, intent(out) :: y,dy -#ifdef POLINT_LEGACY_ORDER +#if USE_FMISC_SAFE_MODE || defined(POLINT_LEGACY_ORDER) integer :: i,j,m,n real*8, dimension(ordn,ordn) :: yatmp real*8, dimension(ordn) :: ymtmp @@ -1502,12 +1517,23 @@ if(dabs(X(1)-xmin) < dX) imin = 1 if(dabs(Y(1)-ymin) < dY) jmin = 1 if(dabs(Z(1)-zmin) < dZ) kmin = 1 -! Optimized with oneMKL BLAS DDOT for dot product -n_elements = (imax-imin+1)*(jmax-jmin+1)*(kmax-kmin+1) -allocate(f_flat(n_elements)) -f_flat = reshape(f(imin:imax,jmin:jmax,kmin:kmax), [n_elements]) -f_out = DDOT(n_elements, f_flat, 1, f_flat, 1) -deallocate(f_flat) +#if USE_FMISC_SAFE_MODE + f_out = 0.d0 + do k = kmin, kmax + do j = jmin, jmax + do i = imin, imax + f_out = f_out + f(i,j,k)*f(i,j,k) + end do + end do + end do +#else +! Optimized with oneMKL BLAS DDOT for dot product + n_elements = (imax-imin+1)*(jmax-jmin+1)*(kmax-kmin+1) + allocate(f_flat(n_elements)) + f_flat = reshape(f(imin:imax,jmin:jmax,kmin:kmax), [n_elements]) + f_out = DDOT(n_elements, f_flat, 1, f_flat, 1) + deallocate(f_flat) +#endif f_out = f_out*dX*dY*dZ @@ -1565,7 +1591,9 @@ if(dabs(Z(1)-zmin) < dZ) kmin = 1 do k=kmin,kmax do j=jmin,jmax +#if !USE_FMISC_SAFE_MODE !DIR$ SIMD REDUCTION(+:s1,s2,s3,s4,s5,s6,s7) +#endif do i=imin,imax s1 = s1 + f1(i,j,k)*f1(i,j,k) s2 = s2 + f2(i,j,k)*f2(i,j,k) @@ -1672,12 +1700,23 @@ if(Symmetry==2)then if(dabs(ymin+gw*dY)