dpchci Subroutine

private subroutine dpchci(n, h, slope, d, incfd)

Set interior derivatives for DPCHIC

Called by dpchic to set derivatives needed to determine a monotone piecewise cubic Hermite interpolant to the data.

Default boundary conditions are provided which are compatible with monotonicity. If the data are only piecewise monotonic, the interpolant will have an extremum at each point where monotonicity switches direction.

To facilitate two-dimensional applications, includes an increment between successive values of the D-array.

The resulting piecewise cubic Hermite function should be identical (within roundoff error) to that produced by dpchim.

Programming notes

  1. The function dpchst(ARG1,ARG2) is assumed to return zero if either argument is zero, +1 if they are of the same sign, and -1 if they are of opposite sign.

Warning

This routine does no validity-checking of arguments.

Arguments

Type IntentOptional Attributes Name
integer, intent(in) :: n

number of data points. If N=2, simply does linear interpolation.
real(kind=wp), intent(in) :: h(*)

array of interval lengths.
real(kind=wp), intent(in) :: slope(*)

array of data slopes. If the data are (X(I),Y(I)), I=1(1)N, then these inputs are:

  • H(I) = X(I+1)-X(I),
  • SLOPE(I) = (Y(I+1)-Y(I))/H(I), I=1(1)N-1.
real(kind=wp), intent(out) :: d(incfd,*)

array of derivative values at data points. If the data are monotonic, these values will determine a a monotone cubic Hermite function. The value corresponding to X(I) is stored in D(1+(I-1)*INCFD), I=1(1)N. No other entries in D are changed.
integer, intent(in) :: incfd

increment between successive values in D. This argument is provided primarily for 2-D applications.

Calls

proc~~dpchci~~CallsGraph proc~dpchci pchip_module::dpchci proc~dpchst pchip_module::dpchst proc~dpchci->proc~dpchst

Called by

proc~~dpchci~~CalledByGraph proc~dpchci pchip_module::dpchci proc~dpchic pchip_module::dpchic proc~dpchic->proc~dpchci

Source Code

    subroutine dpchci (n, h, slope, d, incfd)

    integer,intent(in)   :: n           !! number of data points.
                                        !! If N=2, simply does linear interpolation.
    integer,intent(in)   :: incfd       !! increment between successive values in D.
                                        !! This argument is provided primarily for 2-D applications.
    real(wp),intent(in)  :: h(*)        !! array of interval lengths.
    real(wp),intent(in)  :: slope(*)    !! array of data slopes.
                                        !! If the data are (X(I),Y(I)), I=1(1)N, then these inputs are:
                                        !!
                                        !! * H(I) =  X(I+1)-X(I),
                                        !! * SLOPE(I) = (Y(I+1)-Y(I))/H(I),  I=1(1)N-1.
    real(wp),intent(out) :: d(incfd,*)  !! array of derivative values at data points.
                                        !! If the data are monotonic, these values will determine a
                                        !! a monotone cubic Hermite function.
                                        !! The value corresponding to X(I) is stored in
                                        !!      `D(1+(I-1)*INCFD),  I=1(1)N`.
                                        !! No other entries in D are changed.

    integer :: i, nless1
    real(wp) :: del1, del2, dmax, dmin, drat1, drat2, hsum, hsumt3, w1, w2

    nless1 = n - 1
    del1 = slope(1)

    if (nless1 <= 1)  then
        ! special case n=2 -- use linear interpolation.

        d(1,1) = del1
        d(1,n) = del1

    else
        ! normal case  (n >= 3).

        del2 = slope(2)

        ! set d(1) via non-centered three-point formula, adjusted to be
        ! shape-preserving.

        hsum = h(1) + h(2)
        w1 = (h(1) + hsum)/hsum
        w2 = -h(1)/hsum
        d(1,1) = w1*del1 + w2*del2
        if ( dpchst(d(1,1),del1) <= zero) then
            d(1,1) = zero
        else if ( dpchst(del1,del2) < zero) then
            ! need do this check only if monotonicity switches.
            dmax = three*del1
            if (abs(d(1,1)) > abs(dmax))  d(1,1) = dmax
        endif

        ! loop through interior points.

        do i = 2, nless1

            if (i /= 2) then
                hsum = h(i-1) + h(i)
                del1 = del2
                del2 = slope(i)
            end if

            ! set d(i)=0 unless data are strictly monotonic.
            d(1,i) = zero
            if ( dpchst(del1,del2) > zero) then
                ! use brodlie modification of butland formula.
                hsumt3 = hsum+hsum+hsum
                w1 = (hsum + h(i-1))/hsumt3
                w2 = (hsum + h(i)  )/hsumt3
                dmax = max( abs(del1), abs(del2) )
                dmin = min( abs(del1), abs(del2) )
                drat1 = del1/dmax
                drat2 = del2/dmax
                d(1,i) = dmin/(w1*drat1 + w2*drat2)
            end if

        end do

        ! set d(n) via non-centered three-point formula, adjusted to be
        ! shape-preserving.

        w1 = -h(n-1)/hsum
        w2 = (h(n-1) + hsum)/hsum
        d(1,n) = w1*del1 + w2*del2
        if ( dpchst(d(1,n),del2) <= zero) then
            d(1,n) = zero
        else if ( dpchst(del1,del2) < zero) then
            ! need do this check only if monotonicity switches.
            dmax = three*del2
            if (abs(d(1,n)) > abs(dmax))  d(1,n) = dmax
        endif

    end if

    end subroutine dpchci