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

Trigger-Discipline: port conservative build and fmisc optimizations

Browse Source
This commit is contained in:
CGH0S7
2026-04-24 09:09:50 +08:00
parent 79af79d471
commit 34fe3e6aa5
6 changed files with 718 additions and 350 deletions
Download Patch File Download Diff File Expand all files Collapse all files

715
AMSS_NCKU_source/fmisc.f90
View File

@@ -324,8 +324,7 @@ subroutine symmetry_bd(ord,extc,func,funcc,SoA)
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
funcc(-i,1:extc(2),1:extc(3)) = funcc(i+2,1:extc(2),1:extc(3))*SoA(1)
enddo
@@ -350,8 +349,7 @@ subroutine symmetry_tbd(ord,extc,func,funcc,SoA)
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
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)
@@ -379,8 +377,7 @@ subroutine symmetry_stbd(ord,extc,func,funcc,SoA)
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
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)
@@ -886,17 +883,20 @@ 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)
enddo
do i=0,ord-1
funcc(:,-i,1:extc(3)) = funcc(:,i+1,1:extc(3))*SoA(2)
enddo
do i=0,ord-1
funcc(:,:,-i) = funcc(:,:,i+1)*SoA(3)
enddo
!DIR$ SIMD VECTORLENGTHFOR(KNOWN_INTEGER=8)
funcc(1:extc(1),1:extc(2),1:extc(3)) = func
!DIR$ SIMD VECTORLENGTHFOR(KNOWN_INTEGER=8)
do i=0,ord-1
funcc(-i,1:extc(2),1:extc(3)) = funcc(i+1,1:extc(2),1:extc(3))*SoA(1)
enddo
!DIR$ SIMD VECTORLENGTHFOR(KNOWN_INTEGER=8)
do i=0,ord-1
funcc(:,-i,1:extc(3)) = funcc(:,i+1,1:extc(3))*SoA(2)
enddo
!DIR$ SIMD VECTORLENGTHFOR(KNOWN_INTEGER=8)
do i=0,ord-1
funcc(:,:,-i) = funcc(:,:,i+1)*SoA(3)
enddo
end subroutine symmetry_bd
@@ -912,8 +912,7 @@ subroutine symmetry_tbd(ord,extc,func,funcc,SoA)
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
funcc(-i,1:extc(2),1:extc(3)) = funcc(i+1,1:extc(2),1:extc(3))*SoA(1)
funcc(extc(1)+1+i,1:extc(2),1:extc(3)) = funcc(extc(1)-i,1:extc(2),1:extc(3))*SoA(1)
@@ -941,8 +940,7 @@ subroutine symmetry_stbd(ord,extc,func,funcc,SoA)
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
funcc(-i,1:extc(2),1:extc(3)) = funcc(i+1,1:extc(2),1:extc(3))*SoA(1)
funcc(extc(1)+1+i,1:extc(2),1:extc(3)) = funcc(extc(1)-i,1:extc(2),1:extc(3))*SoA(1)
@@ -1113,151 +1111,353 @@ end subroutine d2dump
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! common code for cell and vertex
!------------------------------------------------------------------------------
! Lagrangian polynomial interpolation
!------------------------------------------------------------------------------
subroutine polint(xa,ya,x,y,dy,ordn)
implicit none
!~~~~~~> Input Parameter:
integer,intent(in) :: ordn
real*8, dimension(ordn), intent(in) :: xa,ya
real*8, intent(in) :: x
real*8, intent(out) :: y,dy
!~~~~~~> Other parameter:
integer :: m,n,ns
real*8, dimension(ordn) :: c,d,den,ho
real*8 :: dif,dift
!~~~~~~>
n=ordn
m=ordn
c=ya
d=ya
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)
ns=ns-1
do m=1,n-1
den(1:n-m)=ho(1:n-m)-ho(1+m:n)
if (any(den(1:n-m) == 0.0))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
end do
return
end subroutine polint
!------------------------------------------------------------------------------
!
! interpolation in 2 dimensions, follow yx order
!
!------------------------------------------------------------------------------
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:
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)
return
end subroutine polin2
!------------------------------------------------------------------------------
!
! interpolation in 3 dimensions, follow zyx order
!
!------------------------------------------------------------------------------
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:
integer :: i,j,m,n
real*8, dimension(ordn,ordn) :: yatmp
real*8, dimension(ordn) :: ymtmp
real*8, dimension(ordn) :: yntmp
real*8, dimension(ordn) :: yqtmp
m=size(x1a)
n=size(x2a)
do i=1,m
do j=1,n
yqtmp=ya(i,j,:)
call polint(x3a,yqtmp,x3,yatmp(i,j),dy,ordn)
end do
yntmp=yatmp(i,:)
call polint(x2a,yntmp,x2,ymtmp(i),dy,ordn)
end do
call polint(x1a,ymtmp,x1,y,dy,ordn)
return
end subroutine polin3
! common code for cell and vertex
!------------------------------------------------------------------------------
! Lagrangian polynomial interpolation
!------------------------------------------------------------------------------
#ifndef POLINT6_USE_BARYCENTRIC
#define POLINT6_USE_BARYCENTRIC 1
#endif
!DIR$ ATTRIBUTES FORCEINLINE :: polint6_neville
subroutine polint6_neville(xa, ya, x, y, dy)
implicit none
real*8, dimension(6), intent(in) :: xa, ya
real*8, intent(in) :: x
real*8, intent(out) :: y, dy
integer :: i, m, ns, n_m
real*8, dimension(6) :: c, d, ho
real*8 :: dif, dift, hp, h, den_val
c = ya
d = ya
ho = xa - x
ns = 1
dif = abs(x - xa(1))
do i = 2, 6
dift = abs(x - xa(i))
if (dift < dif) then
ns = i
dif = dift
end if
end do
y = ya(ns)
ns = ns - 1
do m = 1, 5
n_m = 6 - 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
end if
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 polint6_neville
!DIR$ ATTRIBUTES FORCEINLINE :: polint6_barycentric
subroutine polint6_barycentric(xa, ya, x, y, dy)
implicit none
real*8, dimension(6), intent(in) :: xa, ya
real*8, intent(in) :: x
real*8, intent(out) :: y, dy
integer :: i, j
logical :: is_uniform
real*8, dimension(6) :: lambda
real*8 :: dx, den_i, term, num, den, step, tol
real*8, parameter :: c_uniform(6) = (/ -1.d0, 5.d0, -10.d0, 10.d0, -5.d0, 1.d0 /)
do i = 1, 6
if (x == xa(i)) then
y = ya(i)
dy = 0.d0
return
end if
end do
step = xa(2) - xa(1)
is_uniform = (step /= 0.d0)
if (is_uniform) then
tol = 64.d0 * epsilon(1.d0) * max(1.d0, abs(step))
do i = 3, 6
if (abs((xa(i) - xa(i-1)) - step) > tol) then
is_uniform = .false.
exit
end if
end do
end if
if (is_uniform) then
num = 0.d0
den = 0.d0
do i = 1, 6
term = c_uniform(i) / (x - xa(i))
num = num + term * ya(i)
den = den + term
end do
y = num / den
dy = 0.d0
return
end if
do i = 1, 6
den_i = 1.d0
do j = 1, 6
if (j /= i) then
dx = xa(i) - xa(j)
if (dx == 0.0d0) then
write(*,*) 'failure in polint for point',x
write(*,*) 'with input points: ',xa
stop
end if
den_i = den_i * dx
end if
end do
lambda(i) = 1.d0 / den_i
end do
num = 0.d0
den = 0.d0
do i = 1, 6
term = lambda(i) / (x - xa(i))
num = num + term * ya(i)
den = den + term
end do
y = num / den
dy = 0.d0
return
end subroutine polint6_barycentric
!DIR$ ATTRIBUTES FORCEINLINE :: polint
subroutine polint(xa, ya, x, y, dy, ordn)
implicit none
integer, intent(in) :: ordn
real*8, dimension(ordn), intent(in) :: xa, ya
real*8, intent(in) :: x
real*8, intent(out) :: y, dy
integer :: i, m, ns, n_m
real*8, dimension(ordn) :: c, d, ho
real*8 :: dif, dift, hp, h, den_val
if (ordn == 6) then
#if POLINT6_USE_BARYCENTRIC
call polint6_barycentric(xa, ya, x, y, dy)
#else
call polint6_neville(xa, ya, x, y, dy)
#endif
return
end if
c = ya
d = ya
ho = xa - x
ns = 1
dif = abs(x - xa(1))
do i = 2, ordn
dift = abs(x - xa(i))
if (dift < dif) then
ns = i
dif = dift
end if
end do
y = ya(ns)
ns = ns - 1
do m = 1, ordn - 1
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
end if
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
!------------------------------------------------------------------------------
! Compute Lagrange interpolation basis weights for one target point.
!------------------------------------------------------------------------------
!DIR$ ATTRIBUTES FORCEINLINE :: polint_lagrange_weights
subroutine polint_lagrange_weights(xa, x, w, ordn)
implicit none
integer, intent(in) :: ordn
real*8, dimension(1:ordn), intent(in) :: xa
real*8, intent(in) :: x
real*8, dimension(1:ordn), intent(out) :: w
integer :: i, j
real*8 :: num, den, dx
do i = 1, ordn
num = 1.d0
den = 1.d0
do j = 1, ordn
if (j /= i) then
dx = xa(i) - xa(j)
if (dx == 0.0d0) then
write(*,*) 'failure in polint for point',x
write(*,*) 'with input points: ',xa
stop
end if
num = num * (x - xa(j))
den = den * dx
end if
end do
w(i) = num / den
end do
return
end subroutine polint_lagrange_weights
!------------------------------------------------------------------------------
!
! interpolation in 2 dimensions, follow yx order
!
!------------------------------------------------------------------------------
subroutine polin2(x1a,x2a,ya,x1,x2,y,dy,ordn)
implicit none
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
#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
!------------------------------------------------------------------------------
!
! interpolation in 3 dimensions, follow zyx order
!
!------------------------------------------------------------------------------
subroutine polin3(x1a,x2a,x3a,ya,x1,x2,x3,y,dy,ordn)
implicit none
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
#ifdef POLINT_LEGACY_ORDER
integer :: i,j,m,n
real*8, dimension(ordn,ordn) :: yatmp
real*8, dimension(ordn) :: ymtmp
real*8, dimension(ordn) :: yntmp
real*8, dimension(ordn) :: yqtmp
m=size(x1a)
n=size(x2a)
do i=1,m
do j=1,n
yqtmp=ya(i,j,:)
call polint(x3a,yqtmp,x3,yatmp(i,j),dy,ordn)
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 :: i, j, k
real*8, dimension(ordn) :: w1, w2
real*8, dimension(ordn) :: ymtmp
real*8 :: yx_sum, x_sum
call polint_lagrange_weights(x1a, x1, w1, ordn)
call polint_lagrange_weights(x2a, x2, w2, ordn)
do k = 1, ordn
yx_sum = 0.d0
do j = 1, ordn
x_sum = 0.d0
do i = 1, ordn
x_sum = x_sum + w1(i) * ya(i,j,k)
end do
yx_sum = yx_sum + w2(j) * x_sum
end do
ymtmp(k) = yx_sum
end do
call polint(x3a, ymtmp, x3, y, dy, ordn)
#endif
return
end subroutine polin3
!--------------------------------------------------------------------------------------
! calculate L2norm
subroutine l2normhelper(ex, X, Y, Z,xmin,ymin,zmin,xmax,ymax,zmax,&
@@ -1276,7 +1476,9 @@ end subroutine d2dump
real*8 :: dX, dY, dZ
integer::imin,jmin,kmin
integer::imax,jmax,kmax
integer::i,j,k
integer::i,j,k,n_elements
real*8, dimension(:), allocatable :: f_flat
real*8, external :: DDOT
dX = X(2) - X(1)
dY = Y(2) - Y(1)
@@ -1300,15 +1502,91 @@ if(dabs(X(1)-xmin) < dX) imin = 1
if(dabs(Y(1)-ymin) < dY) jmin = 1
if(dabs(Z(1)-zmin) < dZ) kmin = 1
f_out = sum(f(imin:imax,jmin:jmax,kmin:kmax)*f(imin:imax,jmin:jmax,kmin:kmax))
n_elements = (imax-imin+1)*(jmax-jmin+1)*(kmax-kmin+1)
allocate(f_flat(n_elements))
f_flat = reshape(f(imin:imax,jmin:jmax,kmin:kmax), [n_elements])
f_out = DDOT(n_elements, f_flat, 1, f_flat, 1)
deallocate(f_flat)
f_out = f_out*dX*dY*dZ
return
end subroutine l2normhelper
!--------------------------------------------------------------------------------------
! calculate L2norm especially for shell Blocks
!--------------------------------------------------------------------------------------
subroutine l2normhelper7(ex, X, Y, Z,xmin,ymin,zmin,xmax,ymax,zmax,&
f1,f2,f3,f4,f5,f6,f7,f_out,gw)
implicit none
!~~~~~~> Input parameters:
integer,intent(in ):: ex(1:3)
real*8, intent(in ):: X(1:ex(1)),Y(1:ex(2)),Z(1:ex(3)),xmin,ymin,zmin,xmax,ymax,zmax
integer,intent(in)::gw
real*8, dimension(ex(1),ex(2),ex(3)),intent(in) :: f1,f2,f3,f4,f5,f6,f7
real*8, intent(out) :: f_out(7)
!~~~~~~> Other variables:
real*8 :: dX, dY, dZ
integer::imin,jmin,kmin
integer::imax,jmax,kmax
integer::i,j,k
real*8 :: s1,s2,s3,s4,s5,s6,s7
dX = X(2) - X(1)
dY = Y(2) - Y(1)
dZ = Z(2) - Z(1)
imin = gw+1
jmin = gw+1
kmin = gw+1
imax = ex(1) - gw
jmax = ex(2) - gw
kmax = ex(3) - gw
if(dabs(X(ex(1))-xmax) < dX) imax = ex(1)
if(dabs(Y(ex(2))-ymax) < dY) jmax = ex(2)
if(dabs(Z(ex(3))-zmax) < dZ) kmax = ex(3)
if(dabs(X(1)-xmin) < dX) imin = 1
if(dabs(Y(1)-ymin) < dY) jmin = 1
if(dabs(Z(1)-zmin) < dZ) kmin = 1
s1 = 0.d0
s2 = 0.d0
s3 = 0.d0
s4 = 0.d0
s5 = 0.d0
s6 = 0.d0
s7 = 0.d0
do k=kmin,kmax
do j=jmin,jmax
!DIR$ SIMD REDUCTION(+:s1,s2,s3,s4,s5,s6,s7)
do i=imin,imax
s1 = s1 + f1(i,j,k)*f1(i,j,k)
s2 = s2 + f2(i,j,k)*f2(i,j,k)
s3 = s3 + f3(i,j,k)*f3(i,j,k)
s4 = s4 + f4(i,j,k)*f4(i,j,k)
s5 = s5 + f5(i,j,k)*f5(i,j,k)
s6 = s6 + f6(i,j,k)*f6(i,j,k)
s7 = s7 + f7(i,j,k)*f7(i,j,k)
enddo
enddo
enddo
f_out(1) = s1*dX*dY*dZ
f_out(2) = s2*dX*dY*dZ
f_out(3) = s3*dX*dY*dZ
f_out(4) = s4*dX*dY*dZ
f_out(5) = s5*dX*dY*dZ
f_out(6) = s6*dX*dY*dZ
f_out(7) = s7*dX*dY*dZ
return
end subroutine l2normhelper7
!--------------------------------------------------------------------------------------
! calculate L2norm especially for shell Blocks
subroutine l2normhelper_sh(ex, X, Y, Z,xmin,ymin,zmin,xmax,ymax,zmax,&
f,f_out,gw,ogw,Symmetry)
@@ -1325,7 +1603,9 @@ f_out = f_out*dX*dY*dZ
real*8 :: dX, dY, dZ
integer::imin,jmin,kmin
integer::imax,jmax,kmax
integer::i,j,k
integer::i,j,k,n_elements
real*8, dimension(:), allocatable :: f_flat
real*8, external :: DDOT
real*8 :: PIo4
@@ -1388,7 +1668,11 @@ if(Symmetry==2)then
if(dabs(ymin+gw*dY)<dY.and.Y(1)<0.d0) jmin = gw+1
endif
f_out = sum(f(imin:imax,jmin:jmax,kmin:kmax)*f(imin:imax,jmin:jmax,kmin:kmax))
n_elements = (imax-imin+1)*(jmax-jmin+1)*(kmax-kmin+1)
allocate(f_flat(n_elements))
f_flat = reshape(f(imin:imax,jmin:jmax,kmin:kmax), [n_elements])
f_out = DDOT(n_elements, f_flat, 1, f_flat, 1)
deallocate(f_flat)
f_out = f_out*dX*dY*dZ
@@ -1415,7 +1699,9 @@ f_out = f_out*dX*dY*dZ
real*8 :: dX, dY, dZ
integer::imin,jmin,kmin
integer::imax,jmax,kmax
integer::i,j,k
integer::i,j,k
real*8, dimension(:), allocatable :: f_flat
real*8, external :: DDOT
real*8 :: PIo4
@@ -1478,11 +1764,11 @@ if(Symmetry==2)then
if(dabs(ymin+gw*dY)<dY.and.Y(1)<0.d0) jmin = gw+1
endif
f_out = sum(f(imin:imax,jmin:jmax,kmin:kmax)*f(imin:imax,jmin:jmax,kmin:kmax))
f_out = f_out
Nout = (imax-imin+1)*(jmax-jmin+1)*(kmax-kmin+1)
Nout = (imax-imin+1)*(jmax-jmin+1)*(kmax-kmin+1)
allocate(f_flat(Nout))
f_flat = reshape(f(imin:imax,jmin:jmax,kmin:kmax), [Nout])
f_out = DDOT(Nout, f_flat, 1, f_flat, 1)
deallocate(f_flat)
return
@@ -1583,9 +1869,12 @@ Nout = (imax-imin+1)*(jmax-jmin+1)*(kmax-kmin+1)
! ^
! f=3/8*f_1 + 3/4*f_2 - 1/8*f_3
real*8,parameter::C1=3.d0/8.d0,C2=3.d0/4.d0,C3=-1.d0/8.d0
fout = C1*f1+C2*f2+C3*f3
real*8,parameter::C1=3.d0/8.d0,C2=3.d0/4.d0,C3=-1.d0/8.d0
integer :: i,j,k
do concurrent (k=1:ext(3), j=1:ext(2), i=1:ext(1))
fout(i,j,k) = C1*f1(i,j,k)+C2*f2(i,j,k)+C3*f3(i,j,k)
end do
return
@@ -1679,7 +1968,8 @@ Nout = (imax-imin+1)*(jmax-jmin+1)*(kmax-kmin+1)
real*8, dimension(ORDN,ORDN,ORDN) :: ya
real*8, dimension(ORDN,ORDN) :: tmp2
real*8, dimension(ORDN) :: tmp1
real*8, dimension(3) :: SoAh
real*8, dimension(3) :: SoAh
real*8, external :: DDOT
! +1 because c++ gives 0 for first point
cxB = inds+1
@@ -1725,10 +2015,7 @@ Nout = (imax-imin+1)*(jmax-jmin+1)*(kmax-kmin+1)
tmp1 = tmp1 + coef(ORDN+m)*tmp2(:,m)
enddo
f_int=0
do m=1,ORDN
f_int = f_int + coef(m)*tmp1(m)
enddo
f_int = DDOT(ORDN, coef(1:ORDN), 1, tmp1, 1)
return
@@ -1757,7 +2044,8 @@ Nout = (imax-imin+1)*(jmax-jmin+1)*(kmax-kmin+1)
integer,dimension(2) :: cxB,cxT
real*8, dimension(ORDN,ORDN) :: ya
real*8, dimension(ORDN) :: tmp1
real*8, dimension(2) :: SoAh
real*8, dimension(2) :: SoAh
real*8, external :: DDOT
! +1 because c++ gives 0 for first point
cxB = inds(1:2)+1
@@ -1792,10 +2080,7 @@ Nout = (imax-imin+1)*(jmax-jmin+1)*(kmax-kmin+1)
tmp1 = tmp1 + coef(ORDN+m)*ya(:,m)
enddo
f_int=0
do m=1,ORDN
f_int = f_int + coef(m)*tmp1(m)
enddo
f_int = DDOT(ORDN, coef(1:ORDN), 1, tmp1, 1)
return
@@ -1821,11 +2106,12 @@ Nout = (imax-imin+1)*(jmax-jmin+1)*(kmax-kmin+1)
!~~~~~~> Other parameters:
real*8, dimension(-ORDN+1:ex(1)+ORDN,-ORDN+1:ex(2)+ORDN,ex(3)) :: fh
integer :: m
integer :: cxB,cxT
real*8, dimension(ORDN) :: ya
real*8 :: SoAh
integer,dimension(3) :: inds
integer :: m
integer :: cxB,cxT
real*8, dimension(ORDN) :: ya
real*8 :: SoAh
integer,dimension(3) :: inds
real*8, external :: DDOT
! +1 because c++ gives 0 for first point
inds = indsi + 1
@@ -1886,10 +2172,7 @@ Nout = (imax-imin+1)*(jmax-jmin+1)*(kmax-kmin+1)
write(*,*)"error in global_interpind1d, not recognized dumyd = ",dumyd
endif
f_int=0
do m=1,ORDN
f_int = f_int + coef(m)*ya(m)
enddo
f_int = DDOT(ORDN, coef, 1, ya, 1)
return
@@ -2125,20 +2408,28 @@ Nout = (imax-imin+1)*(jmax-jmin+1)*(kmax-kmin+1)
implicit none
integer,intent(in) :: N
real*8 :: gont
integer :: i
real*8 :: gont
integer :: i
real*8, parameter, dimension(0:20) :: fact_table = [ &
1.d0, 1.d0, 2.d0, 6.d0, 24.d0, 120.d0, 720.d0, 5040.d0, 40320.d0, &
362880.d0, 3628800.d0, 39916800.d0, 479001600.d0, 6227020800.d0, &
87178291200.d0, 1307674368000.d0, 20922789888000.d0, &
355687428096000.d0, 6402373705728000.d0, 121645100408832000.d0, &
2432902008176640000.d0 ]
! sanity check
if(N < 0)then
write(*,*) "ffact: error input for factorial"
return
endif
gont = 1.d0
do i=1,N
gont = gont*i
enddo
if(N < 0)then
write(*,*) "ffact: error input for factorial"
gont = 1.d0
return
endif
if(N <= 20)then
gont = fact_table(N)
else
gont = exp(log_gamma(dble(N+1)))
endif
return