pchipMod.f90 Source File

!! Piecewise Cubic Hermite Interpolating Polynomial calculations

MODULE pchipmod
  !! Subroutines and functions related to the calculation of the Piecewise Cubic Hermite Interpolating Polynomial (PCHIP)

  USE splinec

  IMPLICIT NONE
  PUBLIC
  SAVE

  REAL (KIND=8), PRIVATE                 :: h
  REAL (KIND=8), PRIVATE                 :: fprime_r, fprime_i

CONTAINS

  SUBROUTINE PCHIP( x, y, N, PolyCoef, csWork )

    ! This implements the monotone piecewise cubic Hermite interpolating
    ! polynomial (PCHIP) algorithm. This is a new variant of monotone PCHIP
    ! (paper submitted to JASA). Also see;
    !
    ! F. N. Fritsch and J. Butland, "A Method for Constructing Local Monotone
    ! Piecewise Cubic Interpolants", SIAM Journal on Scientific and Statistical
    ! Computing, 5(2):300-304, (1984) https://doi.org/10.1137/0905021
    !
    ! F. N. Fritsch and R. E. Carlson. "Monotone Piecewise Cubic Interpolation",
    ! SIAM Journal on Numerical Analysis, 17(2):238-246, (1980) https://doi.org/10.1137/0717021
    !
    ! N is the number of nodes
    ! x is a vector of the abscissa values
    ! y is a vector of the associated ordinate values
    ! PolyCoef are the coefficients of the standard polynomial
    ! csWork is a temporary work space for the cubic spline

    INTEGER,          INTENT( IN  )   :: N
    REAL    (KIND=8), INTENT( IN  )   :: x(:)
    COMPLEX (KIND=8), INTENT( IN  )   :: y(:)
    COMPLEX (KIND=8), INTENT( INOUT ) :: PolyCoef( 4, N ), csWork( 4, N )

    INTEGER           :: ix, iBCBeg, iBCEnd
    REAL     (KIND=8) :: h1, h2
    COMPLEX  (KIND=8) :: del1, del2, f1, f2, f1prime, f2prime, fprimeT

    !  Precompute estimates of the derivatives at the nodes
    !  The vector PolyCoef(1,*) holds the ordinate values at the nodes
    !  The vector PolyCoef(2,*) holds the ordinate derivatives at the nodes

    IF ( N == 2 ) THEN

    ! handle special case of two data points seperately (linear interpolation)

    PolyCoef( 1, 1 ) = y( 1 )
    PolyCoef( 2, 1 ) = ( y( 2 ) - y( 1 ) ) / ( x( 2 ) - x( 1 ) )
    PolyCoef( 3, 1 ) = 0.0D0
    PolyCoef( 4, 1 ) = 0.0D0

    ELSE

    ! general case of more than two data points

    PolyCoef( 1, 1 : N ) = y( 1 : N )

    ! left endpoint (non-centered 3-point difference formula)

    CALL h_del( x, y, 2, h1, h2, del1, del2 )
    fprimeT = ( ( 2.0D0 * h1 + h2 ) * del1 - h1 * del2 ) / ( h1 + h2 )
    PolyCoef( 2, 1 ) = fprime_left_end_Cmplx( del1, del2, fprimeT )

    ! right endpoint (non-centered 3-point difference formula)

    CALL h_del( x, y, N - 1, h1, h2, del1, del2 )
    fprimeT = ( -h2 * del1 + ( h1 + 2.0D0 * h2 ) * del2 ) / ( h1 + h2 )
    PolyCoef( 2, N ) = fprime_right_end_Cmplx( del1, del2, fprimeT )

    ! compute coefficients of the cubic spline interpolating polynomial

    iBCBeg = 1   ! specified derivatives at the end points
    iBCEnd = 1
    csWork( 1, 1 : N ) = PolyCoef( 1, 1 : N )
    csWork( 2,     1 ) = PolyCoef( 2, 1 )
    csWork( 2,     N ) = PolyCoef( 2, N )
    CALL CSpline( x, csWork( 1, 1 ), N, iBCBeg, iBCEnd, N )

    ! interior nodes (use derivatives from the cubic spline as initial estimate)

    DO ix = 2, N - 1
       CALL h_del( x, y, ix, h1, h2, del1, del2 )
       ! check if the derivative from the cubic spline satisfies monotonicity
       PolyCoef( 2, ix ) = fprime_interior_Cmplx( del1, del2, csWork( 2, ix ) )
    END DO

    !                                                               2      3
    ! compute coefficients of std cubic polynomial: c0 + c1*x + c2*x + c3*x
    !

    DO ix = 1, N - 1
       h  = x( ix + 1 ) - x( ix )

       f1 = PolyCoef( 1, ix )
       f2 = PolyCoef( 1, ix + 1 )

       f1prime = PolyCoef( 2, ix )
       f2prime = PolyCoef( 2, ix + 1 )

       PolyCoef( 3, ix ) = ( 3.0D0 * ( f2 - f1 )  - h * ( 2.0D0 * f1prime + f2prime ) ) / h**2
       PolyCoef( 4, ix ) = ( h * ( f1prime + f2prime ) - 2.0D0 * ( f2 - f1 ) ) / h**3
    END DO

    END IF

    RETURN
  END SUBROUTINE PCHIP

  ! **********************************************************************!

  SUBROUTINE h_del( x, y, ix, h1, h2, del1, del2 )

    INTEGER,          INTENT( IN  ) :: ix   ! index of the center point
    REAL    (KIND=8), INTENT( IN  ) :: x(:)
    COMPLEX (KIND=8), INTENT( IN  ) :: y(:)
    REAL    (KIND=8), INTENT( OUT ) :: h1, h2
    COMPLEX (KIND=8), INTENT( OUT ) :: del1, del2

    h1   =   x( ix     ) - x( ix - 1 )
    h2   =   x( ix + 1 ) - x( ix     )

    del1 = ( y( ix     ) - y( ix - 1 ) ) / h1
    del2 = ( y( ix + 1 ) - y( ix     ) ) / h2

    RETURN
  END SUBROUTINE h_del

  ! **********************************************************************!

  FUNCTION fprime_interior_Cmplx( del1, del2, fprime )

    COMPLEX (KIND=8), INTENT( IN ) :: del1, del2, fprime
    COMPLEX (KIND=8)               :: fprime_interior_Cmplx

    fprime_r = fprime_interior( REAL(  del1 ), REAL(  del2 ), REAL(  fprime ) )
    fprime_i = fprime_interior( AIMAG( del1 ), AIMAG( del2 ), AIMAG( fprime ) )

    fprime_interior_Cmplx = CMPLX( fprime_r, fprime_i, KIND=8 )

  END FUNCTION fprime_interior_Cmplx

  ! **********************************************************************!

  FUNCTION fprime_left_end_Cmplx( del1, del2, fprime )

    COMPLEX (KIND=8), INTENT( IN ) :: del1, del2, fprime
    COMPLEX (KIND=8)               :: fprime_left_end_Cmplx

    fprime_r = fprime_left_end( REAL(  del1 ), REAL(  del2 ), REAL(  fprime ) )
    fprime_i = fprime_left_end( AIMAG( del1 ), AIMAG( del2 ), AIMAG( fprime ) )

    fprime_left_end_Cmplx = CMPLX( fprime_r, fprime_i, KIND=8 )

  END FUNCTION fprime_left_end_Cmplx

  ! **********************************************************************!

  FUNCTION fprime_right_end_Cmplx( del1, del2, fprime )

    COMPLEX (KIND=8), INTENT( IN ) :: del1, del2, fprime
    COMPLEX (KIND=8)               :: fprime_right_end_Cmplx

    fprime_r = fprime_right_end( REAL(  del1 ), REAL(  del2 ), REAL(  fprime ) )
    fprime_i = fprime_right_end( AIMAG( del1 ), AIMAG( del2 ), AIMAG( fprime ) )

    fprime_right_end_Cmplx = CMPLX( fprime_r, fprime_i, KIND=8 )

  END FUNCTION fprime_right_end_Cmplx

 ! **********************************************************************!

  FUNCTION fprime_interior( del1, del2, fprime )

    REAL (KIND=8), INTENT( IN ) :: del1, del2, fprime
    REAL (KIND=8)               :: fprime_interior

    ! check if derivative is within the trust region, project into it if not

    IF ( del1 * del2 > 0.0 ) THEN
      ! adjacent secant slopes have the same sign, enforce monotonicity
      IF ( del1 > 0.0 ) THEN
        fprime_interior = MIN( MAX(fprime, 0.0D0), 3.0D0 * MIN(del1, del2) )
      ELSE
        fprime_interior = MAX( MIN(fprime, 0.0D0), 3.0D0 * MAX(del1, del2) )
      END IF
    ELSE
      ! force the interpolant to have an extrema here
      fprime_interior = 0.0D0;
    END IF

  END FUNCTION fprime_interior

  ! **********************************************************************!

  FUNCTION fprime_left_end( del1, del2, fprime )

    REAL (KIND=8), INTENT( IN ) :: del1, del2, fprime
    REAL (KIND=8)               :: fprime_left_end

    fprime_left_end = fprime

    IF ( del1 * fprime <= 0.0D0 ) THEN
       ! set derivative to zero if the sign differs from sign of secant slope
       fprime_left_end = 0.0;
    ELSE IF ( ( del1 * del2 <= 0.0D0 ) .AND. ( ABS( fprime ) > ABS( 3.0D0 * del1 ) ) ) THEN
       ! adjust derivative value to enforce monotonicity
       fprime_left_end = 3.0D0 * del1;
    END IF

  END FUNCTION fprime_left_end

  ! **********************************************************************!

  FUNCTION fprime_right_end( del1, del2, fprime )

    ! This is essentially the same as fprime_left_end( del2, del1, fprime )
    ! Written separately for clarity

    REAL (KIND=8), INTENT( IN ) :: del1, del2, fprime
    REAL (KIND=8)               :: fprime_right_end

    fprime_right_end = fprime

    IF ( del2 * fprime <= 0.0D0 ) THEN
       ! set derivative to zero if the sign differs from sign of secant slope
       fprime_right_end = 0.0;
    ELSE IF ( ( del1 * del2 <= 0.0D0 ) .AND. ( ABS( fprime ) > ABS( 3.0D0 * del2 ) ) ) THEN
       ! adjust derivative value to enforce monotonicity
       fprime_right_end = 3.0D0 * del2;
    END IF

  END FUNCTION fprime_right_end

END MODULE pchipmod