对fmisc.f90的polint修改
This commit is contained in:
@@ -1117,146 +1117,137 @@ 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
|
||||
real*8, dimension(ordn), intent(in) :: xa,ya
|
||||
integer, intent(in) :: ordn
|
||||
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
|
||||
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
|
||||
! Initialization
|
||||
c = ya
|
||||
d = ya
|
||||
ho = xa - x
|
||||
|
||||
ns = 1
|
||||
dif = abs(x - xa(1))
|
||||
|
||||
! Find the index of the closest table entry
|
||||
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,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)
|
||||
y = ya(ns)
|
||||
ns = ns - 1
|
||||
|
||||
! Main Neville's algorithm loop
|
||||
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
|
||||
|
||||
! Check for division by zero locally
|
||||
if (den_val == 0.0d0) then
|
||||
write(*,*) 'failure in polint for point',x
|
||||
write(*,*) 'with input points: ',xa
|
||||
stop
|
||||
end if
|
||||
|
||||
! Reuse den_val to avoid redundant divisions
|
||||
den_val = (c(i+1) - d(i)) / den_val
|
||||
|
||||
! Update c and d in place
|
||||
d(i) = h * den_val
|
||||
c(i) = hp * den_val
|
||||
end do
|
||||
|
||||
! Decide which path (up or down the tableau) to take
|
||||
if (2 * ns < n_m) then
|
||||
dy = c(ns + 1)
|
||||
else
|
||||
dy=d(ns)
|
||||
ns=ns-1
|
||||
dy = d(ns)
|
||||
ns = ns - 1
|
||||
end if
|
||||
y=y+dy
|
||||
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)
|
||||
subroutine polin2(x1a,x2a,ya,x1,x2,y,dy,ordn)
|
||||
implicit none
|
||||
integer,intent(in) :: ordn
|
||||
real*8, dimension(ordn), intent(in) :: x1a,x2a
|
||||
real*8, dimension(ordn,ordn), intent(in) :: ya
|
||||
real*8, intent(in) :: x1,x2
|
||||
real*8, intent(out) :: y,dy
|
||||
|
||||
implicit none
|
||||
integer :: j
|
||||
real*8, dimension(ordn) :: ymtmp
|
||||
real*8 :: dy_temp ! Local variable to prevent overwriting result
|
||||
|
||||
!~~~~~~> 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
|
||||
! Optimized sequence: Loop over columns (j)
|
||||
! ya(:,j) is a contiguous memory block in Fortran
|
||||
do j=1,ordn
|
||||
call polint(x1a, ya(:,j), x1, ymtmp(j), dy_temp, ordn)
|
||||
end do
|
||||
|
||||
!~~~~~~> 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
|
||||
! Final interpolation on the results
|
||||
call polint(x2a, ymtmp, x2, 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)
|
||||
subroutine polin3(x1a,x2a,x3a,ya,x1,x2,x3,y,dy,ordn)
|
||||
implicit none
|
||||
integer,intent(in) :: ordn
|
||||
real*8, dimension(ordn), intent(in) :: x1a,x2a,x3a
|
||||
real*8, dimension(ordn,ordn,ordn), intent(in) :: ya
|
||||
real*8, intent(in) :: x1,x2,x3
|
||||
real*8, intent(out) :: y,dy
|
||||
|
||||
implicit none
|
||||
integer :: j, k
|
||||
real*8, dimension(ordn,ordn) :: yatmp
|
||||
real*8, dimension(ordn) :: ymtmp
|
||||
real*8 :: dy_temp
|
||||
|
||||
!~~~~~~> 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
|
||||
! Sequence change: Process the contiguous first dimension (x1) first.
|
||||
! We loop through the 'slow' planes (j, k) to extract 'fast' columns.
|
||||
do k=1,ordn
|
||||
do j=1,ordn
|
||||
! ya(:,j,k) is contiguous; much faster than ya(i,j,:)
|
||||
call polint(x1a, ya(:,j,k), x1, yatmp(j,k), dy_temp, ordn)
|
||||
end do
|
||||
end do
|
||||
|
||||
!~~~~~~> Other parameters:
|
||||
! Now process the second dimension
|
||||
do k=1,ordn
|
||||
call polint(x2a, yatmp(:,k), x2, ymtmp(k), dy_temp, ordn)
|
||||
end do
|
||||
|
||||
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
|
||||
! Final dimension
|
||||
call polint(x3a, ymtmp, x3, y, dy, ordn)
|
||||
|
||||
return
|
||||
end subroutine polin3
|
||||
!--------------------------------------------------------------------------------------
|
||||
! calculate L2norm
|
||||
@@ -2272,3 +2263,4 @@ subroutine find_maximum(ext,X,Y,Z,fun,val,pos,llb,uub)
|
||||
return
|
||||
|
||||
end subroutine
|
||||
|
||||
|
||||
Reference in New Issue
Block a user