dpchia Function

public function dpchia(n, x, f, d, incfd, skip, a, b, ierr) result(value)

Piecewise Cubic Hermite Integrator, Arbitrary limits

Evaluates the definite integral of the cubic Hermite function defined by N, X, F, D over the interval [A, B].

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

Programming notes

  1. The error flag from dpchid is tested, because a logic flaw could conceivably result in IERD=-4, which should be reported.

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 .
real(kind=wp), intent(in) :: a

the limits of integration. NOTE: There is no requirement that [A,B] be contained in [X(1),X(N)]. However, the resulting integral value will be highly suspect, if not.
real(kind=wp), intent(in) :: b

the limits of integration. NOTE: There is no requirement that [A,B] be contained in [X(1),X(N)]. However, the resulting integral value will be highly suspect, if not.
integer, intent(inout) :: ierr

error flag.

  • Normal return:
    • IERR = 0 (no errors).
  • Warning errors:
    • IERR = 1 if A is outside the interval [X(1),X(N)].
    • IERR = 2 if B is outside the interval [X(1),X(N)].
    • IERR = 3 if both of the above are true. (Note that this means that either [A,B] contains data interval or the intervals do not intersect at all.)
  • "Recoverable" errors:
    • IERR = -1 if N<2 .
    • IERR = -2 if INCFD<1 .
    • IERR = -3 if the X-array is not strictly increasing. (VALUE will be zero in any of these cases.) NOTE: The above errors are checked in the order listed, and following arguments have NOT been validated.
    • IERR = -4 in case of an error return from DPCHID (which should never occur).

Return Value real(kind=wp)

value of the requested integral.


Calls

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

Source Code

    function dpchia (n, x, f, d, incfd, skip, a, b, 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(inout)   :: ierr         !! error flag.
                                            !!
                                            !! * Normal return:
                                            !!     * IERR = 0  (no errors).
                                            !! * Warning errors:
                                            !!     * IERR = 1  if  A  is outside the interval [X(1),X(N)].
                                            !!     * IERR = 2  if  B  is outside the interval [X(1),X(N)].
                                            !!     * IERR = 3  if both of the above are true.  (Note that this
                                            !!               means that either [A,B] contains data interval
                                            !!               or the intervals do not intersect at all.)
                                            !! * "Recoverable" errors:
                                            !!     * IERR = -1  if N<2 .
                                            !!     * IERR = -2  if INCFD<1 .
                                            !!     * IERR = -3  if the X-array is not strictly increasing.
                                            !!       (VALUE will be zero in any of these cases.)
                                            !!      NOTE:  The above errors are checked in the order listed,
                                            !!          and following arguments have **NOT** been validated.
                                            !!     * IERR = -4  in case of an error return from DPCHID (which
                                            !!                should never occur).
    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).
    real(wp),intent(in)     :: a            !! the limits of integration.
                                            !! NOTE:  There is no requirement that [A,B] be contained in
                                            !!        [X(1),X(N)].  However, the resulting integral value
                                            !!        will be highly suspect, if not.
    real(wp),intent(in)     :: b            !! the limits of integration.
                                            !! NOTE:  There is no requirement that [A,B] be contained in
                                            !!        [X(1),X(N)].  However, the resulting integral value
                                            !!        will be highly suspect, if not.
    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 .
    real(wp) :: value !! value of the requested integral.

    integer :: i, ia, ib, ierd, il, ir
    real(wp) :: xa, xb

    value = zero

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

    !  function definition is ok, go on.
    skip = .true.
    ierr = 0
    if ( (a<x(1)) .or. (a>x(n)) )  ierr = ierr + 1
    if ( (b<x(1)) .or. (b>x(n)) )  ierr = ierr + 2

    !  compute integral value.

    if (a /= b) then
        xa = min (a, b)
        xb = max (a, b)
        if (xb <= x(2)) then
            ! interval is to left of x(2), so use first cubic.
            value = dchfie (x(1),x(2), f(1,1),f(1,2), &
                            d(1,1),d(1,2), a, b)
        else if (xa >= x(n-1)) then
            ! interval is to right of x(n-1), so use last cubic.
            value = dchfie(x(n-1),x(n), f(1,n-1),f(1,n), &
                        d(1,n-1),d(1,n), a, b)
        else
            ! 'normal' case -- xa<xb, xa<x(n-1), xb>x(2).
            ! locate ia and ib such that
            ! x(ia-1)<xa<=x(ia)<=x(ib)<=xb<=x(ib+1)
            ia = 1
            do i = 1, n-1
                if (xa > x(i))  ia = i + 1
            end do
            ! ia = 1 implies xa<x(1) .  otherwise,
            ! ia is largest index such that x(ia-1)<xa,.

            ib = n
            do i = n, ia, -1
                if (xb < x(i))  ib = i - 1
            end do
            ! ib = n implies xb>x(n) .  otherwise,
            ! ib is smallest index such that xb<x(ib+1) .

            ! compute the integral.
            if (ib < ia) then
                ! this means ib = ia-1 and
                ! (a,b) is a subset of (x(ib),x(ia)).
                value = dchfie (x(ib),x(ia), f(1,ib),f(1,ia), &
                                d(1,ib),d(1,ia), a, b)
            else

                ! first compute integral over (x(ia),x(ib)).
                ! (case (ib == ia) is taken care of by initialization
                ! of value to zero.)
                if (ib > ia) then
                    value = dpchid (n, x, f, d, incfd, skip, ia, ib, ierd)
                    if (ierd < 0) then
                        ! trouble in dpchid.  (should never occur.)
                        ierr = -4
                        call xermsg ('PCHIP', 'dpchia', 'trouble in dpchid', ierr, 1)
                        return
                    end if
                endif

                ! then add on integral over (xa,x(ia)).
                if (xa < x(ia)) then
                    il = max(1, ia-1)
                    ir = il + 1
                    value = value + dchfie (x(il),x(ir), f(1,il),f(1,ir), &
                                            d(1,il),d(1,ir), xa, x(ia))
                endif

                ! then add on integral over (x(ib),xb).
                if (xb > x(ib)) then
                    ir = min (ib+1, n)
                    il = ir - 1
                    value = value + dchfie (x(il),x(ir), f(1,il),f(1,ir), &
                                            d(1,il),d(1,ir), x(ib), xb)
                endif

                ! finally, adjust sign if necessary.
                if (a > b)  value = -value
            endif
        endif
    endif

    end function dpchia