dbint4 – bspline-fortran

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

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

Parameters ibcl, ibcr, fbcl, fbcr allow the specification of one derivative of the spline (of order 1 to k-1) 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 has n=ndata+2 coefficients and is defined on t(k) <= x <= t(n+1), where t(k)=x(1) and t(n+1)=x(ndata). The interior knots are placed at the interior x(i) values for even k, and at midpoints between them for odd k (for k=4 the interior knots are x(2),...,x(ndata-1)). The knots t(1),...,t(k-1) lie to the left of t(k)=x(1) and the knots t(n+2),...,t(n+k) 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(k)=x(1) and t(n+1)=...=t(n+k)=x(ndata) can be specified by kntopt=1.
kntopt=2 selects a knot placement for t(1),...,t(k-1) to make the first 2k-1 knots symmetric about t(k)=x(1) and similarly for t(n+2),...,t(n+k) about t(n+1)=x(ndata).
kntopt=3 allows the user to make their own selection, in increasing order, for t(1),...,t(k-1) to the left of x(1) and t(n+2),...,t(n+k) to the right of x(ndata).

In any case, the interpolation on t(k) <= 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.
Generalized from cubic (k=4) to any order k >= 3.

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 >= max(2,k-2)`
integer(kind=ip),	intent(in)		::	ibcl	selection parameter for left boundary condition: constrain the `ibcl`-th derivative at `x(1)` to `fbcl` (`1 <= ibcl <= k-1`)
integer(kind=ip),	intent(in)		::	ibcr	selection parameter for right boundary condition: constrain the `ibcr`-th derivative at `x(ndata)` to `fbcr` (`1 <= ibcr <= k-1`)
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(k)` and `t(n+1)` to `k` `kntopt = 2` sets a symmetric placement of knots about `t(k)` and `t(n+1)` `kntopt = 3` sets `t(i)=tleft(i)` and `t(n+1+i)=tright(i)`,`i=1,k-1`
real(kind=wp),	intent(in),	dimension(:)	::	tleft	when `kntopt = 3`: `t(1:k-1)` in increasing order to be supplied by the user. (not used for other `kntopt` values)
real(kind=wp),	intent(in),	dimension(:)	::	tright	when `kntopt = 3`: `t(n+2:n+k)` in increasing order to be supplied by the user. (not used for other `kntopt` values)
real(kind=wp),	intent(out),	dimension(:)	::	t	knot array of length `n+k`
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(in)		::	k	order of spline, `k >= 3`
integer(kind=ip),	intent(out)		::	iflag	status flag: 0: no errors 806: `k` is less than 3 2001: `ndata` is less than `max(2,k-2)` 2002: `x` values are not distinct or not ordered 2003: `ibcl` is not in the range `1` to `k-1` 2004: `ibcr` is not in the range `1` to `k-1` 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 2008: `tleft` or `tright` does not have size `k-1`

Calls

Help

Called by

Help

Source Code

    pure subroutine dbint4(x,y,ndata,ibcl,ibcr,fbcl,fbcr,kntopt,tleft,tright,t,bcoef,n,k,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 >= max(2,k-2)`
    integer(ip),intent(in)             :: ibcl    !! selection parameter for left boundary condition:
                                                  !! constrain the `ibcl`-th derivative at `x(1)` to `fbcl`
                                                  !! (`1 <= ibcl <= k-1`)
    integer(ip),intent(in)             :: ibcr    !! selection parameter for right boundary condition:
                                                  !! constrain the `ibcr`-th derivative at `x(ndata)` to `fbcr`
                                                  !! (`1 <= ibcr <= k-1`)
    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(k)` and
                                                  !!   `t(n+1)` to `k`
                                                  !! * `kntopt = 2` sets a symmetric placement of knots
                                                  !!   about `t(k)` and `t(n+1)`
                                                  !! * `kntopt = 3` sets `t(i)=tleft(i)` and
                                                  !!   `t(n+1+i)=tright(i)`,`i=1,k-1`
    real(wp),dimension(:),intent(in)  :: tleft    !! when `kntopt = 3`: `t(1:k-1)` in increasing
                                                  !! order to be supplied by the user.
                                                  !! (not used for other `kntopt` values)
    real(wp),dimension(:),intent(in)  :: tright   !! when `kntopt = 3`: `t(n+2:n+k)` in increasing
                                                  !! order to be supplied by the user.
                                                  !! (not used for other `kntopt` values)
    real(wp),dimension(:),intent(out)  :: t       !! knot array of length `n+k`
    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(in)             :: k       !! order of spline, `k >= 3`
    integer(ip),intent(out)            :: iflag   !! status flag:
                                                  !!
                                                  !! * 0: no errors
                                                  !! * 806: `k` is less than 3
                                                  !! * 2001: `ndata` is less than `max(2,k-2)`
                                                  !! * 2002: `x` values are not distinct or not ordered
                                                  !! * 2003: `ibcl` is not in the range `1` to `k-1`
                                                  !! * 2004: `ibcr` is not in the range `1` to `k-1`
                                                  !! * 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
                                                  !! * 2008: `tleft` or `tright` does not have size `k-1`

    integer(ip) :: i, j, ilo, mflag, np, ml, mu, m, nderiv, icol
    real(wp) :: xl, rowmax
    real(wp),dimension(:),allocatable :: tau       !! abscissa of each equation
    integer(ip),dimension(:),allocatable :: ider   !! derivative order of each equation
    integer(ip),dimension(:),allocatable :: ileft  !! knot interval of each equation
    real(wp),dimension(:,:),allocatable :: abd     !! banded matrix of the system
    integer(ip),dimension(:),allocatable :: ipvt   !! pivot indices
    real(wp),dimension(:,:),allocatable :: vnikx   !! basis function values from [[dbspvd]]
    real(wp),dimension(:),allocatable :: work      !! work array for [[dbspvd]]
    real(wp),dimension(:),allocatable :: colmax    !! scale factor of each column of the system

    n = ndata + 2_ip

    if (k<3_ip) then
        iflag = 806_ip ! k is less than 3
        return
    end if

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

    do i=1_ip,ndata-1_ip
        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>k-1_ip) then
        iflag = 2003_ip ! ibcl is not in the range 1 to k-1
        return
    end if

    if (ibcr<1_ip .or. ibcr>k-1_ip) then
        iflag = 2004_ip ! ibcr is not in the range 1 to k-1
        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

    if (kntopt==3_ip) then
        if (size(tleft,kind=ip)/=k-1_ip .or. size(tright,kind=ip)/=k-1_ip) then
            iflag = 2008_ip ! tleft or tright does not have size k-1
            return
        end if
    end if

    iflag = 0_ip
    np = n + 1_ip

    ! the equations, in order: f(x(1)), the left boundary condition,
    ! f(x(2)),...,f(x(ndata-1)), the right boundary condition, f(x(ndata))
    allocate(tau(n), ider(n), ileft(n))
    tau(1:2) = x(1)
    tau(3:ndata) = x(2:ndata-1)
    tau(n-1:n) = x(ndata)
    ider = 0_ip
    ider(2) = ibcl
    ider(n-1) = ibcr

    ! interior knots, centered on the equation abscissae
    t(k) = x(1)
    t(np) = x(ndata)
    do i=1_ip,n-k
        if (mod(k,2_ip)==0_ip) then
            t(k+i) = tau(i+k/2_ip)
        else
            t(k+i) = 0.5_wp*(tau(i+(k-1_ip)/2_ip)+tau(i+(k+1_ip)/2_ip))
        end if
    end do

    select case (kntopt)
    case(1_ip)
        ! set up knot array with multiplicity k at x(1) and x(ndata)
        do i=1_ip,k-1_ip
            t(k-i) = x(1)
            t(np+i) = x(ndata)
        end do
    case(2_ip)
        ! set up knot array with symmetric placement about end points
        if (n>=2_ip*k-2_ip) then
            do i=1_ip,k-1_ip
                t(k-i) = 2.0_wp*t(k) - t(k+i)
                t(np+i) = 2.0_wp*t(np) - t(np-i)
            end do
        else
            xl = (x(ndata)-x(1))/real(k-1_ip,wp)
            do i=1_ip,k-1_ip
                t(k-i) = t(k-i+1_ip) - xl
                t(np+i) = t(np+i-1_ip) + xl
            end do
        end if
    case(3_ip)
        ! 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:k-1) = tleft
        t(n+2:n+k) = tright
        do i=1_ip,k-1_ip
            if ((t(k-i)>t(k-i+1_ip)) .or. (t(np+i)<t(np+i-1_ip))) then
                iflag = 2006_ip ! knot input through tleft, tright is not ordered properly
                return
            end if
        end do
    end select

    ! knot interval for each equation. right limits at x(1) and
    ! left limits at x(ndata), so the end intervals are k and n.
    ilo = 1_ip
    do i=1_ip,n
        if (i<=2_ip) then
            ileft(i) = k
        else if (i>=n-1_ip) then
            ileft(i) = n
        else
            call dintrv(t,np,tau(i),ilo,ileft(i),mflag)
        end if
    end do

    ! bandwidth of the system (row i has nonzeros in columns ileft-k+1:ileft)
    ml = 0_ip
    mu = 0_ip
    do i=1_ip,n
        ml = max(ml, i-(ileft(i)-k+1_ip))
        mu = max(mu, ileft(i)-i)
    end do
    m = ml + mu + 1_ip

    ! set up the banded system: A(i,j) is stored in abd(i-j+m,j)
    allocate(abd(2_ip*ml+mu+1_ip,n), ipvt(n))
    allocate(vnikx(k,k), work((k+1_ip)*(k+2_ip)/2_ip))
    ! right hand side
    bcoef(1) = y(1)
    bcoef(2) = fbcl
    bcoef(3:ndata) = y(2:ndata-1)
    bcoef(n-1) = fbcr
    bcoef(n) = y(ndata)

    abd = 0.0_wp
    do i=1_ip,n
        nderiv = ider(i) + 1_ip
        call dbspvd(t,k,nderiv,tau(i),ileft(i),k,vnikx,work,iflag)
        if (iflag/=0_ip) return
        ! scale each equation by its largest coefficient, since the derivative
        ! rows can be orders of magnitude larger than the interpolation rows
        rowmax = maxval(abs(vnikx(:,nderiv)))
        if (rowmax==0.0_wp) rowmax = 1.0_wp
        do j=1_ip,k
            icol = ileft(i) - k + j
            abd(i-icol+m,icol) = vnikx(j,nderiv) / rowmax
        end do
        bcoef(i) = bcoef(i) / rowmax
    end do

    ! also scale each column (unknown), since a basis function can be small
    ! at all of the data points (e.g., when the data spacing near an end is
    ! very uneven), and its pivot would then look singular. The scale factors
    ! are the largest powers of 2 not greater than the largest entry of each
    ! column, so the scaling is exact and does not change the solution (partial
    ! pivoting is not affected by column scaling), and does not enlarge the
    ! scaled solution (which could overflow). The solution is unscaled below.
    allocate(colmax(n))
    do j=1_ip,n
        colmax(j) = maxval(abs(abd(:,j)))
        if (colmax(j)==0.0_wp) then
            colmax(j) = 1.0_wp
        else
            colmax(j) = scale(1.0_wp, exponent(colmax(j))-1)
        end if
        abd(:,j) = abd(:,j) / colmax(j)
    end do

    ! solve system of equations
    ! (the rows and columns are scaled so the largest entry of each is about 1,
    !  so an absolute tolerance on the pivots detects a numerically singular system)
    call dgbfa(abd,n,ml,mu,ipvt,1000.0_wp*epsilon(1.0_wp),iflag)
    if (iflag/=0_ip) then
        iflag = 2007_ip  ! the system of equations is singular
    else
        call dgbsl(abd,n,ml,mu,ipvt,bcoef(1:n))
        bcoef(1:n) = bcoef(1:n) / colmax
    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, iflag)

Error conditions

See also

History

Arguments

Calls

Called by

Source Code