dbintk Subroutine

private pure subroutine dbintk(x, y, t, n, k, bcoef, q, work, iflag)

dbintk produces the b-spline coefficients, bcoef, of the b-spline of order k with knots t(i), i=1,...,n+k, which takes on the value y(i) at x(i), i=1,...,n. the spline or any of its derivatives can be evaluated by calls to dbvalu.

the i-th equation of the linear system a*bcoef = b for the coefficients of the interpolant enforces interpolation at x(i), i=1,...,n. hence, b(i) = y(i), for all i, and a is a band matrix with 2k-1 bands if a is invertible. the matrix a is generated row by row and stored, diagonal by diagonal, in the rows of q, with the main diagonal going into row k. the banded system is then solved by a call to dbnfac (which constructs the triangular factorization for a and stores it again in q), followed by a call to dbnslv (which then obtains the solution bcoef by substitution). dbnfac does no pivoting, since the total positivity of the matrix a makes this unnecessary. the linear system to be solved is (theoretically) invertible if and only if t(i) < x(i) < t(i+k), for all i. equality is permitted on the left for i=1 and on the right for i=n when k knots are used at x(1) or x(n). otherwise, violation of this condition is certain to lead to an error.

Error conditions

  • improper input
  • singular system of equations

History

  • splint written by carl de boor [5]
  • dbintk author: amos, d. e., (snla) : date written 800901
  • revision date 820801
  • 000330 modified array declarations. (jec)
  • Jacob Williams, 5/10/2015 : converted to free-form Fortran.

Arguments

Type IntentOptional Attributes Name
real(kind=wp), intent(in), dimension(n) :: x

vector of length n containing data point abscissa in strictly increasing order.
real(kind=wp), intent(in), dimension(n) :: y

corresponding vector of length n containing data point ordinates.
real(kind=wp), intent(in), dimension(*) :: t

knot vector of length n+k since t(1),..,t(k) <= x(1) and t(n+1),..,t(n+k)

= x(n), this leaves only n-k knots (not necessarily x(i) values) interior to (x(1),x(n))
integer(kind=ip), intent(in) :: n

number of data points, n >= k
integer(kind=ip), intent(in) :: k

order of the spline, k >= 1
real(kind=wp), intent(out), dimension(n) :: bcoef

a vector of length n containing the b-spline coefficients
real(kind=wp), intent(out), dimension(*) :: q

a work vector of length (2k-1)n, containing the triangular factorization of the coefficient matrix of the linear system being solved. the coefficients for the interpolant of an additional data set (x(i),yy(i)), i=1,...,n with the same abscissa can be obtained by loading yy into bcoef and then executing call dbnslv(q,2k-1,n,k-1,k-1,bcoef)
real(kind=wp), intent(out), dimension(*) :: work

work vector of length 2*k
integer(kind=ip), intent(out) :: iflag
  • 0: no errors.
  • 100: k does not satisfy k>=1.
  • 101: n does not satisfy n>=k.
  • 102: x(i) does not satisfy x(i)<x(i+1) for some i.
  • 103: some abscissa was not in the support of the corresponding basis function and the system is singular.
  • 104: the system of solver detects a singular system. although the theoretical conditions for a solution were satisfied.

Calls

proc~~dbintk~~CallsGraph proc~dbintk bspline_sub_module::dbintk proc~dbnfac bspline_sub_module::dbnfac proc~dbintk->proc~dbnfac proc~dbnslv bspline_sub_module::dbnslv proc~dbintk->proc~dbnslv proc~dbspvn bspline_sub_module::dbspvn proc~dbintk->proc~dbspvn

Called by

proc~~dbintk~~CalledByGraph proc~dbintk bspline_sub_module::dbintk proc~dbtpcf bspline_sub_module::dbtpcf proc~dbtpcf->proc~dbintk proc~db1ink_default bspline_sub_module::db1ink_default proc~db1ink_default->proc~dbtpcf proc~db2ink bspline_sub_module::db2ink proc~db2ink->proc~dbtpcf proc~db3ink bspline_sub_module::db3ink proc~db3ink->proc~dbtpcf proc~db4ink bspline_sub_module::db4ink proc~db4ink->proc~dbtpcf proc~db5ink bspline_sub_module::db5ink proc~db5ink->proc~dbtpcf proc~db6ink bspline_sub_module::db6ink proc~db6ink->proc~dbtpcf interface~db1ink bspline_sub_module::db1ink interface~db1ink->proc~db1ink_default proc~initialize_2d_auto_knots bspline_oo_module::bspline_2d%initialize_2d_auto_knots proc~initialize_2d_auto_knots->proc~db2ink proc~initialize_2d_specify_knots bspline_oo_module::bspline_2d%initialize_2d_specify_knots proc~initialize_2d_specify_knots->proc~db2ink proc~initialize_3d_auto_knots bspline_oo_module::bspline_3d%initialize_3d_auto_knots proc~initialize_3d_auto_knots->proc~db3ink proc~initialize_3d_specify_knots bspline_oo_module::bspline_3d%initialize_3d_specify_knots proc~initialize_3d_specify_knots->proc~db3ink proc~initialize_4d_auto_knots bspline_oo_module::bspline_4d%initialize_4d_auto_knots proc~initialize_4d_auto_knots->proc~db4ink proc~initialize_4d_specify_knots bspline_oo_module::bspline_4d%initialize_4d_specify_knots proc~initialize_4d_specify_knots->proc~db4ink proc~initialize_5d_auto_knots bspline_oo_module::bspline_5d%initialize_5d_auto_knots proc~initialize_5d_auto_knots->proc~db5ink proc~initialize_5d_specify_knots bspline_oo_module::bspline_5d%initialize_5d_specify_knots proc~initialize_5d_specify_knots->proc~db5ink proc~initialize_6d_auto_knots bspline_oo_module::bspline_6d%initialize_6d_auto_knots proc~initialize_6d_auto_knots->proc~db6ink proc~initialize_6d_specify_knots bspline_oo_module::bspline_6d%initialize_6d_specify_knots proc~initialize_6d_specify_knots->proc~db6ink proc~bspline_2d_constructor_auto_knots bspline_oo_module::bspline_2d_constructor_auto_knots proc~bspline_2d_constructor_auto_knots->proc~initialize_2d_auto_knots proc~bspline_2d_constructor_specify_knots bspline_oo_module::bspline_2d_constructor_specify_knots proc~bspline_2d_constructor_specify_knots->proc~initialize_2d_specify_knots proc~bspline_3d_constructor_auto_knots bspline_oo_module::bspline_3d_constructor_auto_knots proc~bspline_3d_constructor_auto_knots->proc~initialize_3d_auto_knots proc~bspline_3d_constructor_specify_knots bspline_oo_module::bspline_3d_constructor_specify_knots proc~bspline_3d_constructor_specify_knots->proc~initialize_3d_specify_knots proc~bspline_4d_constructor_auto_knots bspline_oo_module::bspline_4d_constructor_auto_knots proc~bspline_4d_constructor_auto_knots->proc~initialize_4d_auto_knots proc~bspline_4d_constructor_specify_knots bspline_oo_module::bspline_4d_constructor_specify_knots proc~bspline_4d_constructor_specify_knots->proc~initialize_4d_specify_knots proc~bspline_5d_constructor_auto_knots bspline_oo_module::bspline_5d_constructor_auto_knots proc~bspline_5d_constructor_auto_knots->proc~initialize_5d_auto_knots proc~bspline_5d_constructor_specify_knots bspline_oo_module::bspline_5d_constructor_specify_knots proc~bspline_5d_constructor_specify_knots->proc~initialize_5d_specify_knots proc~bspline_6d_constructor_auto_knots bspline_oo_module::bspline_6d_constructor_auto_knots proc~bspline_6d_constructor_auto_knots->proc~initialize_6d_auto_knots proc~bspline_6d_constructor_specify_knots bspline_oo_module::bspline_6d_constructor_specify_knots proc~bspline_6d_constructor_specify_knots->proc~initialize_6d_specify_knots proc~initialize_1d_auto_knots bspline_oo_module::bspline_1d%initialize_1d_auto_knots proc~initialize_1d_auto_knots->interface~db1ink proc~initialize_1d_specify_knots bspline_oo_module::bspline_1d%initialize_1d_specify_knots proc~initialize_1d_specify_knots->interface~db1ink interface~bspline_2d bspline_oo_module::bspline_2d interface~bspline_2d->proc~bspline_2d_constructor_auto_knots interface~bspline_2d->proc~bspline_2d_constructor_specify_knots interface~bspline_3d bspline_oo_module::bspline_3d interface~bspline_3d->proc~bspline_3d_constructor_auto_knots interface~bspline_3d->proc~bspline_3d_constructor_specify_knots interface~bspline_4d bspline_oo_module::bspline_4d interface~bspline_4d->proc~bspline_4d_constructor_auto_knots interface~bspline_4d->proc~bspline_4d_constructor_specify_knots interface~bspline_5d bspline_oo_module::bspline_5d interface~bspline_5d->proc~bspline_5d_constructor_auto_knots interface~bspline_5d->proc~bspline_5d_constructor_specify_knots interface~bspline_6d bspline_oo_module::bspline_6d interface~bspline_6d->proc~bspline_6d_constructor_auto_knots interface~bspline_6d->proc~bspline_6d_constructor_specify_knots proc~bspline_1d_constructor_auto_knots bspline_oo_module::bspline_1d_constructor_auto_knots proc~bspline_1d_constructor_auto_knots->proc~initialize_1d_auto_knots proc~bspline_1d_constructor_specify_knots bspline_oo_module::bspline_1d_constructor_specify_knots proc~bspline_1d_constructor_specify_knots->proc~initialize_1d_specify_knots interface~bspline_1d bspline_oo_module::bspline_1d interface~bspline_1d->proc~bspline_1d_constructor_auto_knots interface~bspline_1d->proc~bspline_1d_constructor_specify_knots

Source Code

    pure subroutine dbintk(x,y,t,n,k,bcoef,q,work,iflag)

    implicit none

    integer(ip),intent(in)            :: n      !! number of data points, n >= k
    real(wp),dimension(n),intent(in)  :: x      !! vector of length n containing data point abscissa
                                                !! in strictly increasing order.
    real(wp),dimension(n),intent(in)  :: y      !! corresponding vector of length n containing data
                                                !! point ordinates.
    real(wp),dimension(*),intent(in)  :: t      !! knot vector of length n+k
                                                !! since t(1),..,t(k) <= x(1) and t(n+1),..,t(n+k)
                                                !! >= x(n), this leaves only n-k knots (not
                                                !! necessarily x(i) values) interior to (x(1),x(n))
    integer(ip),intent(in)            :: k      !! order of the spline, k >= 1
    real(wp),dimension(n),intent(out) :: bcoef  !! a vector of length n containing the b-spline coefficients
    real(wp),dimension(*),intent(out) :: q      !! a work vector of length (2*k-1)*n, containing
                                                !! the triangular factorization of the coefficient
                                                !! matrix of the linear system being solved.  the
                                                !! coefficients for the interpolant of an
                                                !! additional data set (x(i),yy(i)), i=1,...,n
                                                !! with the same abscissa can be obtained by loading
                                                !! yy into bcoef and then executing
                                                !! call dbnslv(q,2k-1,n,k-1,k-1,bcoef)
    real(wp),dimension(*),intent(out) :: work   !! work vector of length 2*k
    integer(ip),intent(out)           :: iflag  !! *   0: no errors.
                                                !! * 100: k does not satisfy k>=1.
                                                !! * 101: n does not satisfy n>=k.
                                                !! * 102: x(i) does not satisfy x(i)<x(i+1) for some i.
                                                !! * 103: some abscissa was not in the support of the
                                                !! corresponding basis function and the system is singular.
                                                !! * 104: the system of solver detects a singular system.
                                                !! although the theoretical conditions for a solution were satisfied.

    integer(ip) :: iwork, i, ilp1mx, j, jj, km1, kpkm2, left,lenq, np1
    real(wp) :: xi
    logical :: found

    if (k<1_ip) then
        !write(error_unit,'(A)') 'dbintk - k does not satisfy k>=1'
        iflag = 100_ip
        return
    end if

    if (n<k) then
        !write(error_unit,'(A)') 'dbintk - n does not satisfy n>=k'
        iflag = 101_ip
        return
    end if

    jj = n - 1_ip
    if (jj/=0_ip) then
        do i=1_ip,jj
            if (x(i)>=x(i+1_ip)) then
                !write(error_unit,'(A)') 'dbintk - x(i) does not satisfy x(i)<x(i+1) for some i'
                iflag = 102_ip
                return
            end if
        end do
    end if

    np1 = n + 1_ip
    km1 = k - 1_ip
    kpkm2 = 2_ip*km1
    left = k
    ! zero out all entries of q
    lenq = n*(k+km1)
    do i=1_ip,lenq
        q(i) = 0.0_wp
    end do

    ! loop over i to construct the n interpolation equations
    do i=1_ip,n

        xi = x(i)
        ilp1mx = min(i+k,np1)
        ! find left in the closed interval (i,i+k-1_ip) such that
        !         t(left) <= x(i) < t(left+1_ip)
        ! matrix is singular if this is not possible
        left = max(left,i)
        if (xi<t(left)) then
            !write(error_unit,'(A)') 'dbintk - some abscissa was not in the support of the'//&
            !             ' corresponding basis function and the system is singular'
            iflag = 103_ip
            return
        end if
        found = .false.
        do
            found = (xi<t(left+1_ip))
            if (found) exit
            left = left + 1_ip
            if (left>=ilp1mx) exit
        end do
        if (.not. found) then
            left = left - 1_ip
            if (xi>t(left+1_ip)) then
                !write(error_unit,'(A)') 'dbintk - some abscissa was not in the support of the'//&
                !             ' corresponding basis function and the system is singular'
                iflag = 103_ip
                return
            end if
        end if
        ! the i-th equation enforces interpolation at xi, hence
        ! a(i,j) = b(j,k,t)(xi), all j. only the  k  entries with  j =
        ! left-k+1,...,left actually might be nonzero. these  k  numbers
        ! are returned, in  bcoef (used for temp.storage here), by the
        ! following
        call dbspvn(t, k, k, 1_ip, xi, left, bcoef, work, iwork, iflag)
        if (iflag/=0_ip) return

        ! we therefore want  bcoef(j) = b(left-k+j)(xi) to go into
        ! a(i,left-k+j), i.e., into  q(i-(left+j)+2*k,(left+j)-k) since
        ! a(i+j,j)  is to go into  q(i+k,j), all i,j,  if we consider  q
        ! as a two-dim. array , with  2*k-1  rows (see comments in
        ! dbnfac). in the present program, we treat  q  as an equivalent
        ! one-dimensional array (because of fortran restrictions on
        ! dimension statements) . we therefore want  bcoef(j) to go into
        ! entry
        !     i -(left+j) + 2*k + ((left+j) - k-1)*(2*k-1)
        !            = i-left+1 + (left -k)*(2*k-1) + (2*k-2)*j
        ! of q.
        jj = i - left + 1_ip + (left-k)*(k+km1)
        do j=1_ip,k
            jj = jj + kpkm2
            q(jj) = bcoef(j)
        end do

    end do

    ! obtain factorization of a, stored again in q.
    call dbnfac(q, k+km1, n, km1, km1, iflag)

    if (iflag==1) then !success
        ! solve  a*bcoef = y  by backsubstitution
        do i=1_ip,n
            bcoef(i) = y(i)
        end do
        call dbnslv(q, k+km1, n, km1, km1, bcoef)
        iflag = 0_ip
    else  !failure
        !write(error_unit,'(A)') 'dbintk - the system of solver detects a singular system'//&
        !             ' although the theoretical conditions for a solution were satisfied'
        iflag = 104_ip
    end if

    end subroutine dbintk