dpchce Subroutine

private subroutine dpchce(ic, vc, n, x, h, slope, d, incfd, ierr)

Set boundary conditions for dpchic

Called by dpchic to set end derivatives as requested by the user. It must be called after interior derivative values have been set.

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

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.
  2. One could reduce the number of arguments and amount of local storage, at the expense of reduced code clarity, by passing in the array WK (rather than splitting it into H and SLOPE) and increasing its length enough to incorporate STEMP and XTEMP.
  3. The two monotonicity checks only use the sufficient conditions. Thus, it is possible (but unlikely) for a boundary condition to be changed, even though the original interpolant was monotonic. (At least the result is a continuous function of the data.)

Warning

This routine does no validity-checking of arguments.

Arguments

Type IntentOptional Attributes Name
integer, intent(in) :: ic(2)

array of length 2 specifying desired boundary conditions:

  • IC(1) = IBEG, desired condition at beginning of data.
  • IC(2) = IEND, desired condition at end of data.

( see prologue to dpchic for details. )
real(kind=wp), intent(in) :: vc(2)

array of length 2 specifying desired boundary values.

  • VC(1) need be set only if IC(1) = 2 or 3 .
  • VC(2) need be set only if IC(2) = 2 or 3 .
integer, intent(in) :: n

number of data points. (assumes N>=2)
real(kind=wp), intent(in) :: x(*)

array of independent variable values. (the elements of X are assumed to be strictly increasing.)
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(inout) :: d(incfd,*)

(input) real8 array of derivative values at the data points. The value corresponding to X(I) must be stored in D(1+(I-1)INCFD), I=1(1)N.

(output) the value of D at X(1) and/or X(N) is changed, if necessary, to produce the requested boundary conditions. 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.
integer, intent(out) :: ierr

error flag.

  • Normal return:
  • IERR = 0 (no errors).
  • Warning errors:
  • IERR = 1 if IBEG<0 and D(1) had to be adjusted for monotonicity.
  • IERR = 2 if IEND<0 and D(1+(N-1)*INCFD) had to be adjusted for monotonicity.
  • IERR = 3 if both of the above are true.

Calls

proc~~dpchce~~CallsGraph proc~dpchce pchip_module::dpchce proc~dpchdf pchip_module::dpchdf proc~dpchce->proc~dpchdf proc~dpchst pchip_module::dpchst proc~dpchce->proc~dpchst proc~xermsg pchip_module::xermsg proc~dpchce->proc~xermsg proc~dpchdf->proc~xermsg

Called by

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

Source Code

    subroutine dpchce (ic, vc, n, x, h, slope, d, incfd, ierr)

    integer,intent(in)  :: ic(2)    !! array of length 2 specifying desired
                                    !! boundary conditions:
                                    !!
                                    !! * IC(1) = IBEG, desired condition at beginning of data.
                                    !! * IC(2) = IEND, desired condition at end of data.
                                    !!
                                    !! ( see prologue to [[dpchic]] for details. )
    integer,intent(in)     :: n     !! number of data points.  (assumes N>=2)
    integer,intent(in)     :: incfd !! increment between successive values in D.
                                    !! This argument is provided primarily for 2-D applications.
    integer,intent(out)    :: ierr  !! error flag.
                                    !!
                                    !! * Normal return:
                                    !!    * IERR = 0  (no errors).
                                    !! * Warning errors:
                                    !!    * IERR = 1  if IBEG<0 and D(1) had to be adjusted for
                                    !!              monotonicity.
                                    !!    * IERR = 2  if IEND<0 and D(1+(N-1)*INCFD) had to be
                                    !!              adjusted for monotonicity.
                                    !!    * IERR = 3  if both of the above are true.
    real(wp),intent(in)    :: vc(2) !! array of length 2 specifying desired boundary
                                    !! values.
                                    !!
                                    !! * VC(1) need be set only if IC(1) = 2 or 3 .
                                    !! * VC(2) need be set only if IC(2) = 2 or 3 .
    real(wp),intent(in)    :: x(*)  !! array of independent variable values.  (the
                                    !! elements of X are assumed to be strictly increasing.)
    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(inout) :: d(incfd,*) !! (input) real*8 array of derivative values at the data points.
                                         !!  The value corresponding to X(I) must be stored in
                                         !!  D(1+(I-1)*INCFD),  I=1(1)N.
                                         !!
                                         !! (output) the value of D at X(1) and/or X(N) is changed, if
                                         !!  necessary, to produce the requested boundary conditions.
                                         !!  no other entries in D are changed.

    integer :: ibeg, iend, ierf, index, j, k
    real(wp) :: stemp(3), xtemp(4)

    ibeg = ic(1)
    iend = ic(2)
    ierr = 0

    ! set to default boundary conditions if n is too small.

    if ( abs(ibeg)>n )  ibeg = 0
    if ( abs(iend)>n )  iend = 0

    ! treat beginning boundary condition.

    if (ibeg /= 0) then
        k = abs(ibeg)
        if (k == 1) then
            ! boundary value provided.
            d(1,1) = vc(1)
        else if (k == 2) then
            ! boundary second derivative provided.
            d(1,1) = half*( (three*slope(1) - d(1,2)) - half*vc(1)*h(1) )
        else if (k < 5) then
            ! use k-point derivative formula.
            ! pick up first k points, in reverse order.
            do j = 1, k
                index = k-j+1
                ! index runs from k down to 1.
                xtemp(j) = x(index)
                if (j < k)  stemp(j) = slope(index-1)
            end do
            ! -----------------------------
            d(1,1) = dpchdf (k, xtemp, stemp, ierf)
            ! -----------------------------
            if (ierf /= 0) then
                ! *** this case should never occur ***
                ierr = -1
                call xermsg ('PCHIP', 'dpchce', 'error return from dpchdf', ierr, 1)
                return
            end if
        else
            ! use 'not a knot' condition.
            d(1,1) = ( three*(h(1)*slope(2) + h(2)*slope(1)) &
                       - two*(h(1)+h(2))*d(1,2) - h(1)*d(1,3) ) / h(2)
        endif

        if (ibeg <= 0) then

            ! check d(1,1) for compatibility with monotonicity.

            if (slope(1) == zero) then
                if (d(1,1) /= zero) then
                    d(1,1) = zero
                    ierr = ierr + 1
                endif
            else if ( dpchst(d(1,1),slope(1)) < zero) then
                d(1,1) = zero
                ierr = ierr + 1
            else if ( abs(d(1,1)) > three*abs(slope(1)) ) then
                d(1,1) = three*slope(1)
                ierr = ierr + 1
            endif

        end if

    end if

    ! treat end boundary condition.

    if (iend == 0) return
    k = abs(iend)
    if (k == 1) then
        ! boundary value provided.
        d(1,n) = vc(2)
    else if (k == 2) then
        ! boundary second derivative provided.
        d(1,n) = half*( (three*slope(n-1) - d(1,n-1)) + half*vc(2)*h(n-1) )
    else if (k < 5) then
        ! use k-point derivative formula.
        ! pick up last k points.
        do j = 1, k
            index = n-k+j
            ! index runs from n+1-k up to n.
            xtemp(j) = x(index)
            if (j < k)  stemp(j) = slope(index)
        end do
        ! -----------------------------
        d(1,n) = dpchdf (k, xtemp, stemp, ierf)
        ! -----------------------------
        if (ierf /= 0) then
            ! *** this case should never occur ***
            ierr = -1
            call xermsg ('PCHIP', 'dpchce', 'error return from dpchdf', ierr, 1)
            return
        end if
    else
        ! use 'not a knot' condition.
        d(1,n) = ( three*(h(n-1)*slope(n-2) + h(n-2)*slope(n-1)) &
                   - two*(h(n-1)+h(n-2))*d(1,n-1) - h(n-1)*d(1,n-2) ) / h(n-2)
    endif

    if (iend > 0)  return

    ! check d(1,n) for compatibility with monotonicity.

    if (slope(n-1) == zero) then
        if (d(1,n) /= zero) then
            d(1,n) = zero
            ierr = ierr + 2
        endif
    else if ( dpchst(d(1,n),slope(n-1)) < zero) then
        d(1,n) = zero
        ierr = ierr + 2
    else if ( abs(d(1,n)) > three*abs(slope(n-1)) ) then
        d(1,n) = three*slope(n-1)
        ierr = ierr + 2
    endif

    end subroutine dpchce