dbintk produces the bspline coefficients, bcoef, of the bspline 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 ith 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 2k1 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.
Type  Intent  Optional  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)


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 bspline coefficients 

real(kind=wp),  intent(out),  dimension(*)  ::  q 
a work vector of length (2k1)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,2k1,n,k1,k1,bcoef) 

real(kind=wp),  intent(out),  dimension(*)  ::  work 
work vector of length 2*k 

integer(kind=ip),  intent(out)  ::  iflag 

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 nk 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 bspline coefficients real(wp),dimension(*),intent(out) :: q !! a work vector of length (2*k1)*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,2k1,n,k1,k1,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+k1_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 ith equation enforces interpolation at xi, hence ! a(i,j) = b(j,k,t)(xi), all j. only the k entries with j = ! leftk+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(leftk+j)(xi) to go into ! a(i,leftk+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 twodim. array , with 2*k1 rows (see comments in ! dbnfac). in the present program, we treat q as an equivalent ! onedimensional array (because of fortran restrictions on ! dimension statements) . we therefore want bcoef(j) to go into ! entry ! i (left+j) + 2*k + ((left+j)  k1)*(2*k1) ! = ileft+1 + (left k)*(2*k1) + (2*k2)*j ! of q. jj = i  left + 1_ip + (leftk)*(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