dbint4 – bspline-fortran

private pure subroutine dbint4(x, y, ndata, ibcl, ibcr, fbcl, fbcr, kntopt, tleft, tright, t, bcoef, n, k, w, iflag)

DBINT4 computes the B representation (t,bcoef,n,k) of a cubic spline (k=4) which interpolates data (x(i),y(i)),i=1,ndata.

Parameters ibcl, ibcr, fbcl, fbcr allow the specification of the spline first or second derivative at both x(1) and x(ndata). When this data is not specified by the problem, it is common practice to use a natural spline by setting second derivatives at x(1) and x(ndata) to zero (ibcl=ibcr=2,fbcl=fbcr=0.0).

The spline is defined on t(4) <= x <= t(n+1) with (ordered) interior knots at x(i) values where n=ndata+2. The knots t(1),t(2),t(3) lie to the left of t(4)=x(1) and the knots t(n+2), t(n+3), t(n+4) lie to the right of t(n+1)=x(ndata) in increasing order.

If no extrapolation outside (x(1),x(ndata)) is anticipated, the knots t(1)=t(2)=t(3)=t(4)=x(1) and t(n+2)=t(n+3)=t(n+4)=t(n+1)=x(ndata) can be specified by kntopt=1.
kntopt=2 selects a knot placement for t(1), t(2), t(3) to make the first 7 knots symmetric about t(4)=x(1) and similarly for t(n+2), t(n+3), t(n+4) about t(n+1)=x(ndata).
kntopt=3 allows the user to make his own selection, in increasing order, for t(1), t(2), t(3) to the left of x(1) and t(n+2), t(n+3), t(n+4) to the right of x(ndata).

In any case, the interpolation on t(4) <= x <= t(n+1) by using function dbvalu is unique for given boundary conditions.

Error conditions

improper input
singular system of equations

History

Written by D. E. Amos (SNLA), August, 1979.
date written 800901
revision date 820801
000330 Modified array declarations. (JEC)
Jacob Williams, 8/30/2018 : refactored to modern Fortran.

Arguments

Type	Intent	Attributes		Name
real(kind=wp),	intent(in),	dimension(:)	::	x	x vector of abscissae of length `ndata`, distinct and in increasing order
real(kind=wp),	intent(in),	dimension(:)	::	y	y vector of ordinates of length ndata
integer(kind=ip),	intent(in)		::	ndata	number of data points, `ndata >= 2`
integer(kind=ip),	intent(in)		::	ibcl	selection parameter for left boundary condition: `ibcl = 1` constrain the first derivative at `x(1)` to `fbcl` `ibcl = 2` constrain the second derivative at `x(1)` to `fbcl`
integer(kind=ip),	intent(in)		::	ibcr	selection parameter for right boundary condition: `ibcr = 1` constrain first derivative at `x(ndata)` to `fbcr` `ibcr = 2` constrain second derivative at `x(ndata)` to `fbcr`
real(kind=wp),	intent(in)		::	fbcl	left boundary values governed by `ibcl`
real(kind=wp),	intent(in)		::	fbcr	right boundary values governed by `ibcr`
integer(kind=ip),	intent(in)		::	kntopt	knot selection parameter: `kntopt = 1` sets knot multiplicity at `t(4)` and `t(n+1)` to 4 `kntopt = 2` sets a symmetric placement of knots about `t(4)` and `t(n+1)` `kntopt = 3` sets `t(i)=tleft(i)` and `t(n+1+i)=tright(i)`,`i=1,3`
real(kind=wp),	intent(in),	dimension(3)	::	tleft	when `kntopt = 3`: `t(1:3)` in increasing order to be supplied by the user.
real(kind=wp),	intent(in),	dimension(3)	::	tright	when `kntopt = 3`: `t(n+2:n+4)` in increasing order to be supplied by the user.
real(kind=wp),	intent(out),	dimension(:)	::	t	knot array of length `n+4`
real(kind=wp),	intent(out),	dimension(:)	::	bcoef	b spline coefficient array of length `n`
integer(kind=ip),	intent(out)		::	n	number of coefficients, `n=ndata+2`
integer(kind=ip),	intent(out)		::	k	order of spline, `k=4`
real(kind=wp),	intent(inout),	dimension(5,ndata+2)	::	w	work array
integer(kind=ip),	intent(out)		::	iflag	status flag: 0: no errors 2001: `ndata` is less than 2 2002: `x` values are not distinct or not ordered 2003: `ibcl` is not 1 or 2 2004: `ibcr` is not 1 or 2 2005: `kntopt` is not 1, 2, or 3 2006: knot input through `tleft`, `tright` is not ordered properly 2007: the system of equations is singular

Calls

Help

Called by

Help

Source Code

    pure subroutine dbint4(x,y,ndata,ibcl,ibcr,fbcl,fbcr,kntopt,tleft,tright,t,bcoef,n,k,w,iflag)

    implicit none

    real(wp),dimension(:),intent(in)   :: x       !! x vector of abscissae of length `ndata`, distinct
                                                  !! and in increasing order
    real(wp),dimension(:),intent(in)   :: y       !! y vector of ordinates of length ndata
    integer(ip),intent(in)             :: ndata   !! number of data points, `ndata >= 2`
    integer(ip),intent(in)             :: ibcl    !! selection parameter for left boundary condition:
                                                  !!
                                                  !! * `ibcl = 1` constrain the first derivative at `x(1)` to `fbcl`
                                                  !! * `ibcl = 2` constrain the second derivative at `x(1)` to `fbcl`
    integer(ip),intent(in)             :: ibcr    !! selection parameter for right boundary condition:
                                                  !!
                                                  !! * `ibcr = 1` constrain first derivative at `x(ndata)` to `fbcr`
                                                  !! * `ibcr = 2` constrain second derivative at `x(ndata)` to `fbcr`
    real(wp),intent(in)                :: fbcl    !! left boundary values governed by `ibcl`
    real(wp),intent(in)                :: fbcr    !! right boundary values governed by `ibcr`
    integer(ip),intent(in)             :: kntopt  !! knot selection parameter:
                                                  !!
                                                  !! * `kntopt = 1` sets knot multiplicity at `t(4)` and
                                                  !!   `t(n+1)` to 4
                                                  !! * `kntopt = 2` sets a symmetric placement of knots
                                                  !!   about `t(4)` and `t(n+1)`
                                                  !! * `kntopt = 3` sets `t(i)=tleft(i)` and
                                                  !!   `t(n+1+i)=tright(i)`,`i=1,3`
    real(wp),dimension(3),intent(in)  :: tleft    !! when `kntopt = 3`: `t(1:3)` in increasing
                                                  !! order to be supplied by the user.
    real(wp),dimension(3),intent(in)  :: tright   !! when `kntopt = 3`: `t(n+2:n+4)` in increasing
                                                  !! order to be supplied by the user.
    real(wp),dimension(:),intent(out)  :: t       !! knot array of length `n+4`
    real(wp),dimension(:),intent(out)  :: bcoef   !! b spline coefficient array of length `n`
    integer(ip),intent(out)            :: n       !! number of coefficients, `n=ndata+2`
    integer(ip),intent(out)            :: k       !! order of spline, `k=4`
    real(wp),dimension(5,ndata+2),intent(inout) :: w  !! work array
    integer(ip),intent(out)            :: iflag   !! status flag:
                                                  !!
                                                  !! * 0: no errors
                                                  !! * 2001: `ndata` is less than 2
                                                  !! * 2002: `x` values are not distinct or not ordered
                                                  !! * 2003: `ibcl` is not 1 or 2
                                                  !! * 2004: `ibcr` is not 1 or 2
                                                  !! * 2005: `kntopt` is not 1, 2, or 3
                                                  !! * 2006: knot input through `tleft`, `tright` is
                                                  !!   not ordered properly
                                                  !! * 2007: the system of equations is singular

    integer(ip)  :: i, ilb, ileft, it, iub, iw, iwp, j, jw, ndm, np, nwrow
    real(wp) :: txn, tx1, xl
    real(wp),dimension(4,4) :: vnikx
    real(wp),dimension(15) :: work  !! work array for [[dbspvd]] -- length `(k+1)*(k+2)/2`

    real(wp),parameter :: wdtol = epsilon(1.0_wp) !! d1mach(4)
    real(wp),parameter :: tol = sqrt(wdtol)

    if (ndata<2_ip) then
        iflag = 2001_ip ! ndata is less than 2
        return
    end if

    ndm = ndata - 1_ip
    do i=1_ip,ndm
        if (x(i)>=x(i+1_ip)) then
            iflag = 2002_ip ! x values are not distinct or not ordered
            return
        end if
    end do

    if (ibcl<1_ip .or. ibcl>2_ip) then
        iflag = 2003_ip ! ibcl is not 1 or 2
        return
    end if

    if (ibcr<1_ip .or. ibcr>2_ip) then
        iflag = 2004_ip ! ibcr is not 1 or 2
        return
    end if

    if (kntopt<1_ip .or. kntopt>3_ip) then
        iflag = 2005_ip ! kntopt is not 1, 2, or 3
        return
    end if

    iflag = 0_ip
    k = 4_ip
    n = ndata + 2_ip
    np = n + 1_ip
    do i=1_ip,ndata
        t(i+3) = x(i)
    end do

    select case (kntopt)
    case(1_ip)
        ! set up knot array with multiplicity 4 at x(1) and x(ndata)
        do i=1,3_ip
            t(4-i) = x(1)
            t(np+i) = x(ndata)
        end do
    case(2_ip)
        !set up knot array with symmetric placement about end points
        if (ndata>3) then
            tx1 = x(1) + x(1)
            txn = x(ndata) + x(ndata)
            do i=1,3
                t(4-i) = tx1 - x(i+1)
                t(np+i) = txn - x(ndata-i)
            end do
        else
            xl = (x(ndata)-x(1))/3.0_wp
            do i=1,3
                t(4-i) = t(5-i) - xl
                t(np+i) = t(np+i-1) + xl
            end do
        end if
    case(3)
        ! set up knot array less than x(1) and greater than x(ndata) to be
        ! supplied by user in tleft & tright when kntopt=3
        t(1:3)             = tleft
        t(ndata+4:ndata+6) = tright
        do i=1,3
            if ((t(4-i)>t(5-i)) .or. (t(np+i)<t(np+i-1))) then
                iflag = 2006_ip ! knot input through tleft, tright is not ordered properly
                return
            end if
        end do
    end select

    w = 0.0_wp

    ! set up left interpolation point and left boundary condition for
    ! right limits
    it = ibcl + 1
    call dbspvd(t, k, it, x(1), k, 4_ip, vnikx, work, iflag)
    if (iflag/=0_ip) return ! error check
    iw = 0_ip
    if (abs(vnikx(3,1))<tol) iw = 1_ip
    do j=1,3
        w(j+1,4-j) = vnikx(4-j,it)
        w(j,4-j) = vnikx(4-j,1)
    end do
    bcoef(1) = y(1)
    bcoef(2) = fbcl
    ! set up interpolation equations for points i=2 to i=ndata-1
    ileft = 4_ip
    if (ndm>=2) then
        do i=2,ndm
            ileft = ileft + 1_ip
            call dbspvd(t, k, 1_ip, x(i), ileft, 4_ip, vnikx, work, iflag)
            if (iflag/=0_ip) return ! error check
            do j=1,3
                w(j+1,3+i-j) = vnikx(4-j,1)
            end do
            bcoef(i+1) = y(i)
        end do
    end if

    ! set up right interpolation point and right boundary condition for
    ! left limits(ileft is associated with t(n)=x(ndata-1))
    it = ibcr + 1_ip
    call dbspvd(t, k, it, x(ndata), ileft, 4_ip, vnikx, work, iflag)
    if (iflag/=0_ip) return ! error check
    jw = 0_ip
    if (abs(vnikx(2,1))<tol) jw = 1_ip
    do j=1,3
        w(j+1,3+ndata-j) = vnikx(5-j,it)
        w(j+2,3+ndata-j) = vnikx(5-j,1)
    end do
    bcoef(n-1) = fbcr
    bcoef(n) = y(ndata)
    ! solve system of equations
    ilb = 2_ip - jw
    iub = 2_ip - iw
    nwrow = 5_ip
    iwp = iw + 1_ip
    call dbnfac(w(iwp,1), nwrow, n, ilb, iub, iflag)
    if (iflag==2_ip) then
        iflag = 2007_ip  ! the system of equations is singular
    else
        iflag = 0_ip  ! success
        call dbnslv(w(iwp,1), nwrow, n, ilb, iub, bcoef)
    end if

    end subroutine dbint4

dbint4 Subroutine

Contents

Source Code

private pure subroutine dbint4(x, y, ndata, ibcl, ibcr, fbcl, fbcr, kntopt, tleft, tright, t, bcoef, n, k, w, iflag)

Error conditions

See also

History

Arguments

Calls

Called by

Source Code