dpchid Function

public function dpchid(n, x, f, d, incfd, skip, ia, ib, ierr) result(value)

Piecewise Cubic Hermite Integrator, Data limits

Evaluates the definite integral of the cubic Hermite function defined by N, X, F, D over the interval [X(IA), X(IB)].

To provide compatibility with dpchim and dpchic, includes an increment between successive values of the F- and D-arrays.

Programming notes

  1. This routine uses a special formula that is valid only for integrals whose limits coincide with data values. This is mathematically equivalent to, but much more efficient than, calls to dchfie.

Arguments

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

number of data points. (Error return if N<2).
real(kind=wp), intent(in) :: x(*)

array of independent variable values. The elements of X must be strictly increasing: X(I-1) < X(I), I = 2(1)N. (Error return if not.)
real(kind=wp), intent(in) :: f(incfd,*)

array of function values. F(1+(I-1)*INCFD) is the value corresponding to X(I).
real(kind=wp), intent(in) :: d(incfd,*)

array of derivative values. D(1+(I-1)*INCFD) is the value corresponding to X(I).
integer, intent(in) :: incfd

increment between successive values in F and D. (Error return if INCFD<1).
logical, intent(inout) :: skip

logical variable which should be set to .TRUE. if the user wishes to skip checks for validity of preceding parameters, or to .FALSE. otherwise. This will save time in case these checks have already been performed (say, in dpchim or dpchic). SKIP will be set to .TRUE. on return with IERR = 0 or -4.
integer, intent(in) :: ia

indices in X-array for the limits of integration. both must be in the range [1,N]. (Error return if not.) No restrictions on their relative values.
integer, intent(in) :: ib

indices in X-array for the limits of integration. both must be in the range [1,N]. (Error return if not.) No restrictions on their relative values.
integer, intent(out) :: ierr

error flag.

  • Normal return:

  • IERR = 0 (no errors).

  • "Recoverable" errors (VALUE will be zero in any of these cases.):

  • IERR = -1 if N<2 .
  • IERR = -2 if INCFD<1 .
  • IERR = -3 if the X-array is not strictly increasing.
  • IERR = -4 if IA or IB is out of range.

NOTE: The above errors are checked in the order listed, and following arguments have NOT been validated.

Return Value real(kind=wp)

value of the requested integral.


Calls

proc~~dpchid~~CallsGraph proc~dpchid pchip_module::dpchid proc~xermsg pchip_module::xermsg proc~dpchid->proc~xermsg

Called by

proc~~dpchid~~CalledByGraph proc~dpchid pchip_module::dpchid proc~dpchia pchip_module::dpchia proc~dpchia->proc~dpchid

Source Code

    function dpchid (n, x, f, d, incfd, skip, ia, ib, ierr) result(value)

    integer,intent(in)    :: n          !! number of data points.  (Error return if N<2).
    integer,intent(in)    :: incfd      !! increment between successive values in F and D.
                                        !! (Error return if  INCFD<1).
    integer,intent(in)    :: ia         !! indices in X-array for the limits of integration.
                                        !! both must be in the range [1,N].  (Error return if not.)
                                        !! No restrictions on their relative values.
    integer,intent(in)    :: ib         !! indices in X-array for the limits of integration.
                                        !! both must be in the range [1,N].  (Error return if not.)
                                        !! No restrictions on their relative values.
    integer,intent(out)   :: ierr       !! error flag.
                                        !!
                                        !! * Normal return:
                                        !!    * IERR = 0  (no errors).
                                        !!
                                        !! * "Recoverable" errors (VALUE will be zero in any of these cases.):
                                        !!    * IERR = -1  if N<2 .
                                        !!    * IERR = -2  if INCFD<1 .
                                        !!    * IERR = -3  if the X-array is not strictly increasing.
                                        !!    * IERR = -4  if IA or IB is out of range.
                                        !!
                                        !! NOTE: The above errors are checked in the order listed,
                                        !! and following arguments have **NOT** been validated.
    real(wp),intent(in)   :: x(*)       !! array of independent variable values.  The
                                        !! elements of X must be strictly increasing:
                                        !! `X(I-1) < X(I),  I = 2(1)N`.
                                        !! (Error return if not.)
    real(wp),intent(in)   :: f(incfd,*) !! array of function values.  `F(1+(I-1)*INCFD)` is
                                        !! the value corresponding to `X(I)`.
    real(wp),intent(in)   :: d(incfd,*) !! array of derivative values. ` D(1+(I-1)*INCFD)`
                                        !! is the value corresponding to `X(I)`.
    logical,intent(inout) :: skip       !! logical variable which should be set to
                                        !! .TRUE. if the user wishes to skip checks for validity of
                                        !! preceding parameters, or to .FALSE. otherwise.
                                        !! This will save time in case these checks have already
                                        !! been performed (say, in [[dpchim]] or [[dpchic]]).
                                        !! SKIP will be set to .TRUE. on return with IERR = 0 or -4.
    real(wp) :: value   !! value of the requested integral.

    integer :: i, iup, low
    real(wp) :: h, sum

    value = zero

    ! validity-check arguments.
    if (.not. skip) then
        if ( n<2 ) then
            ierr = -1
            call xermsg ('PCHIP', 'dpchid', 'number of data points less than two', ierr, 1)
            return
        end if
        if ( incfd<1 ) then
            ierr = -2
            call xermsg ('PCHIP', 'dpchid', 'increment less than one', ierr,1)
            return
        end if
        do i = 2, n
            if ( x(i)<=x(i-1) ) then
                ! x-array not strictly increasing.
                ierr = -3
                call xermsg ('PCHIP', 'dpchid', 'x-array not strictly increasing', ierr, 1)
                return
            end if
        end do
    end if

    ! function definition is ok, go on.
    skip = .true.
    if ((ia<1) .or. (ia>n) .or. (ib<1) .or. (ib>n)) then
        ! ia or ib out of range return.
        ierr = -4
        call xermsg ('PCHIP', 'dpchid', 'ia or ib out of range', ierr, 1)
        return
    end if

    ierr = 0

    ! compute integral value.

    if (ia /= ib) then
        low = min(ia, ib)
        iup = max(ia, ib) - 1
        sum = zero
        do i = low, iup
            h = x(i+1) - x(i)
            sum = sum + h*( (f(1,i) + f(1,i+1)) + &
                            (d(1,i) - d(1,i+1))*(h/six) )
        end do
        value = half * sum
        if (ia > ib)  value = -value
    endif

    end function dpchid