22 Jun 2012 6:38
#10396
This module crashes:
Error 14, Attempt to alter an actual argument that is a constant, an expression, an INTENT(IN) argument, or a DO variable. .
It appears to work correctly in GNU Fortran. Any help would be appreciated!! 😃
module pchip_module
implicit none
!
! This module smooths the transition between the INTOS and SLATEC interpolators.
!
! 1. INTOS takes X values in either ascending or descending order. SLATEC only ascending.
! This is easy to fix by just flipping the sign of the X array if it is descending.
!
! 2. For each XBAR, INTOS needs to search the X array and then generate the
! parabola that fits those points. SLATEC would prefer to generate all of
! the interpolating polynomials at once and then reuse these polynomials on subsequent
! calls. Not doing this might be very time consuming.
!
! The data type 'polynomials' stores a pointer to the a given set of Y values.
! Y values come in an array, YTB(IDIM1, IDIM2), so that several sets of Y values
! could be stored in one 2D array. INDX selects which set of Y values to use.
!
! If the polynomials have already been generated for the the incoming set of Y values
! the derivatives, 'D', can be retrieved and reused. If not, a new entry is made and
! space is allocated to save those derivatives.
!
! The 'filled' attribute keeps track of which sets of Y values from a given
! 2D YTB array have already been processed and their derivatives stored in 'D'.
!
INTEGER,PARAMETER :: MAXIMUM_POLYNOMIALS = 20
type polynomials
private
REAL, POINTER :: YTBP
REAL, DIMENSION(:,:), POINTER :: D
LOGICAL, DIMENSION(:), POINTER :: filled
end type polynomials
integer, save, private :: numberOfPolynomials = 0
type(polynomials), save, private :: knownPolynomials(MAXIMUM_POLYNOMIALS)
!
! Return the derivatives that have been stored away or need to be generated
! and stored for the given X and YTB(INDX,*).
!
private findDerivatives
contains
function findDerivatives(IDIM1, IDIM2, INDX, NUM2, X, YTB, IERR)
REAL, DIMENSION(:,:), POINTER :: findDerivatives
LOGICAL, DIMENSION(:), POINTER :: filledArray
LOGICAL :: filled
INTEGER :: polynomial
INTEGER, INTENT(IN) :: IDIM1, IDIM2, INDX, NUM2
REAL, DIMENSION(:), INTENT(IN) :: X
REAL, DIMENSION(:,:), INTENT(IN), TARGET :: YTB
INTEGER, INTENT(OUT) :: IERR
REAL, POINTER :: YTBP
INTEGER :: i, allocateStatus
YTBP => YTB(1, 1)
nullify(findDerivatives)
filled = .FALSE.
polynomial = 0
do i = 1, numberOfPolynomials
if ( ASSOCIATED( YTBP, knownPolynomials(i)%YTBP ) ) then
findDerivatives => knownPolynomials(i)%D
filled = knownPolynomials(i)%filled(INDX)
polynomial = i
exit
end if
end do
if (polynomial .EQ. 0) then
if (numberOfPolynomials .LT. MAXIMUM_POLYNOMIALS) then
ALLOCATE (findDerivatives(IDIM1, IDIM2), &
STAT = allocateStatus)
if (allocateStatus /= 0) then
STOP '*** Not enough memory ***'
end if
ALLOCATE(filledArray(IDIM1), STAT = allocateStatus)
if (allocateStatus /= 0) then
STOP '*** Not enough memory ***'
end if
numberOfPolynomials = numberOfPolynomials + 1
!!! write (6,54) numberOfPolynomials
!!!54 format('Number of Polynomials ', I4)
polynomial = numberOfPolynomials
knownPolynomials(polynomial)%YTBP => YTBP
knownPolynomials(polynomial)%D => findDerivatives
knownPolynomials(polynomial)%filled => filledArray
do i = 1, IDIM1
knownPolynomials(polynomial)%filled(i) = .FALSE.
end do
else
STOP '**** Maximum number of Polynomials exceeded ****'
end if
end if
if (.NOT. filled) then
!!! write(6,55) polynomial, INDX
!!!55 format('Filling polynomial/derivative', I4, I4)
call pchim(NUM2, X, YTB(INDX,1), &
findDerivatives(INDX,1), IDIM1, IERR)
knownPolynomials(polynomial)%filled(INDX) = .TRUE.
end if
end function findDerivatives
subroutine interpolate(IDIM1, IDIM2, INDX, NUM2, XBAR, &
XTB, YTB, YBAR, IERR)
!C * * * * * * * * * * * * * * * * * * * * * * *
!C
!C THIS SUBROUTINE PERFORMS PIECEWISE CUBIC HERMITE INTERPOLATION.
!C (CSU = 110805-0096-00)
!C
!C INPUT VARIABLES
!C IDIM1 - FIRST ARRAY DIMENSION (1 = RANK 1)
!C IDIM2 - SECOND ARRAY DIMENSION
!C INDX - DATA LOCATION FOR IDIM1
!C NUM2 - NUMBER OF DATA POINTS
!C XBAR - INDEPENDENT VARIABLE
!C XTB - INDEPENDENT VARIABLE ARRAY
!C YTB - DEPENDENT VARIABLE ARRAY
!C
!C OUTPUT VARIABLES
!C YBAR - INTERPOLATED VALUE
!C IERR - ERROR FLAG (1, IF IDIM2 <= IDEG)
!C
!C * * * * * * * * * * * * * * * * * * * * * * *
!C
INTEGER,INTENT(IN) :: IDIM1, IDIM2, INDX, NUM2
REAL,INTENT(IN) :: XBAR, XTB(IDIM1,IDIM2), YTB(IDIM1,IDIM2)
REAL,INTENT(OUT) :: YBAR
INTEGER,INTENT(OUT) :: IERR
REAL, POINTER :: D(:, :)
REAL, DIMENSION(:) :: X(NUM2)
REAL :: xdiff
REAL :: flip = 1.0
LOGICAL :: SKIP = .TRUE.
INTEGER :: i
INTEGER,PARAMETER :: NE = 1
!
! SLATEC can only support increasing X
!
XDIFF = XTB(INDX, NUM2) - XTB(INDX, 1)
FLIP = SIGN(FLIP, XDIFF)
do i = 1, NUM2
X(i) = FLIP * XTB(INDX, i)
end do
D => findDerivatives(IDIM1, IDIM2, INDX, NUM2, X, YTB, IERR)
call pchfe(NUM2, X, YTB(INDX,1), D(INDX,1), &
IDIM1, SKIP, NE, FLIP * XBAR, YBAR, IERR)
!!! if (YBAR .GT. 201) then
!!! do i = 1, num2
!!! write(6,67) i, indx, XBAR, X(i), ytb(indx, i), D(indx, i)
!!!67 format(2I4, 4F12.4)
!!! end do
!!! end if
end subroutine interpolate
end module pchip_module
Here's the test driver.
program test
implicit none
call backAndForth()
end program test
subroutine backAndForth()
INTEGER, PARAMETER :: IDIM1 = 2
INTEGER, PARAMETER :: IDIM2 = 10
INTEGER, PARAMETER :: NUM2 = 10
INTEGER, PARAMETER :: IDEG = 2
INTEGER, PARAMETER :: INDX = 1
INTEGER, PARAMETER :: INDX2 = 2
REAL, PARAMETER :: TOL = 1.0E-5
INTEGER :: I, IFIRST, IERR
REAL, DIMENSION(IDIM1,IDIM2) :: XTB, YTB
REAL :: XBAR, YBAR, YBAR2
do i = 1, IDIM2
XTB(INDX, i) = i - 1
XTB(INDX2, i) = 10 - i
if (i .GT. 5) THEN
YTB(INDX, i) = 0
YTB(INDX2, i) = 1
else
YTB(INDX, i) = 1
YTB(INDX2, i) = 0
end if
end do
do I = 1, NUM2
write(6,1050) I, XTB(INDX, I), YTB(INDX, I), XTB(INDX2, i), YTB(INDX2, I)
1050 format(I3, 4F12.4)
end do
xbar = 0.0
do i = 1, 2*NUM2
call PCHIP (IDEG, IDIM1, IDIM2, INDX, NUM2, XBAR, &
XTB, YTB, TOL, IFIRST, YBAR, IERR)
call PCHIP (IDEG, IDIM1, IDIM2, INDX2, NUM2, XBAR, &
XTB, YTB, TOL, IFIRST, YBAR2, IERR)
write (6, 100) i, XBAR, YBAR, YBAR2
100 format(i4, 3F14.4)
xbar = xbar + .5
end do
end subroutine backAndForth
A wrapper subroutine....
SUBROUTINE PCHIP (IDEG, IDIM1, IDIM2, INDX, NUM2, XBAR, &
XTB, YTB, TOL, IFIRST, YBAR, IERR)
use pchip_module
!C * * * * * * * * * * * * * * * * * * * * * * *
!C
!C THIS SUBROUTINE PERFORMS PIECEWISE CUBIC HERMITE INTERPOLATION.
!C (CSU = 110805-0096-00)
!C
!C INPUT VARIABLES
!C IDEG - DEGREE OF INTERPOLATION
!C IDIM1 - FIRST ARRAY DIMENSION (1 = RANK 1)
!C IDIM2 - SECOND ARRAY DIMENSION
!C INDX - DATA LOCATION FOR IDIM1
!C NUM2 - NUMBER OF DATA POINTS
!C XBAR - INDEPENDENT VARIABLE
!C XTB - INDEPENDENT VARIABLE ARRAY
!C YTB - DEPENDENT VARIABLE ARRAY
!C TOL - TOLERANCE
!C
!C INPUT/OUTPUT VARIABLE
!C IFIRST - 0, DETERMINE INITIAL INTERPOLATION POINT
!C >0, BEGIN INTERPOLATION AT SUBSCRIPT IFIRST
!C
!C OUTPUT VARIABLES
!C YBAR - INTERPOLATED VALUE
!C IERR - ERROR FLAG (1, IF IDIM2 <= IDEG)
!C
!C
!C THERE ARE NO CALLS TO OTHER ROUTINES
!C
!C THIS ROUTINE IS CALLED BY DCALC, PSCALC, RADGEN, RNGALT
!C
!C * * * * * * * * * * * * * * * * * * * * * * *
!C
INTEGER,INTENT(IN) :: IDIM1, IDIM2, INDX, NUM2
REAL,INTENT(IN) :: XBAR
REAL,INTENT(IN), DIMENSION (IDIM1, IDIM2) :: XTB, YTB
REAL,INTENT(OUT) :: YBAR
INTEGER,INTENT(OUT) :: IERR
REAL :: X0, X1
REAL :: xrange, dToX0, dToX1
! write(6,1000) IDEG, INDX, NUM2, XBAR
!1000 format(3I4, F12.4)
! do 1100 I = 1, NUM2
! write(6,1050) I, XTB(INDX, I), YTB(INDX, I)
!1050 format(I3, F12.4, F12.4)
!1100 continue
!C
!C * * CHECK FOR SUFFICIENT DATA VALUES
IERR = 0
X0 = XTB(INDX,1)
X1 = XTB(INDX, NUM2)
dToX0 = abs(XBAR - X0)
dToX1 = abs(XBAR - X1)
xrange = dToX0 + dToX1
if (xrange .GT. abs( X0 - X1 ) + TOL) then
write(6,666) XBAR, X0, X1
666 format ('Interpolation ', F12.4, ' outside of range ', 2F12.4)
! write(6,667) dToX0, dToX1, abs(X0-X1)
!667 format(3F12.4)
!C
!C * *
!C Return the Y value associated with the closest end
!C * *
!C
if (dToX0 .LT. dToX1) then
YBAR = YTB(INDX, 1)
else
YBAR = YTB(INDX, NUM2)
end if
return
end if
!C
!C * * PROCEED WITH INTERPOLATION
!C
call interpolate(IDIM1, IDIM2, INDX, NUM2, XBAR, XTB, YTB, YBAR, IERR)
!C*** write(6,2000) IFIRST, X0, XBAR, X1, Y0, YBAR, Y1, slope0, slope1
!C***2000 format(I2, 8F12.4)
!C
END SUBROUTINE PCHIP
Two stubs.
SUBROUTINE PCHFE (N, X, F, D, INCFD, SKIP, NE, XE, FE, IERR)
C
C DECLARE ARGUMENTS.
C
INTEGER N, INCFD, NE, IERR
REAL X(*), F(INCFD,*), D(INCFD,*), XE(*), FE(*)
LOGICAL SKIP
RETURN
C------------- LAST LINE OF PCHFE FOLLOWS ------------------------------
END
SUBROUTINE PCHIM (N, X, F, D, INCFD, IERR)
C
C DECLARE ARGUMENTS.
C
INTEGER N, INCFD, IERR
REAL X(*), F(INCFD,*), D(INCFD,*)
RETURN
C------------- LAST LINE OF PCHIM FOLLOWS ------------------------------
END