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.
| 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, | intent(in) | :: | n |
number of data points, n >= k |
||
| integer, | 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, | intent(out) | :: | iflag |
if 0: no errors. if 100: k does not satisfy k>=1. if 101: n does not satisfy n>=k. if 102: x(i) does not satisfy x(i)<x(i+1) for some i. if 103: some abscissa was not in the support of the. corresponding basis function and the system is singular. if 104: the system of solver detects a singular system. although the theoretical conditions for a solution were satisfied. |
subroutine dbintk(x,y,t,n,k,bcoef,q,work,iflag) implicit none integer,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,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,intent(out) :: iflag !! if 0: no errors. !! if 100: k does not satisfy k>=1. !! if 101: n does not satisfy n>=k. !! if 102: x(i) does not satisfy x(i)<x(i+1) for some i. !! if 103: some abscissa was not in the support of the. !! corresponding basis function and the system is singular. !! if 104: the system of solver detects a singular system. !! although the theoretical conditions for a solution were satisfied. integer :: iwork, i, ilp1mx, j, jj, km1, kpkm2, left,lenq, np1 real(wp) :: xi if (k<1) then write(error_unit,'(A)') 'dbintk - k does not satisfy k>=1' iflag = 100 return endif if (n<k) then write(error_unit,'(A)') 'dbintk - n does not satisfy n>=k' iflag = 101 return endif jj = n - 1 if (jj/=0) then do i=1,jj if (x(i)>=x(i+1)) then write(error_unit,'(A)') 'dbintk - x(i) does not satisfy x(i)<x(i+1) for some i' iflag = 102 return endif enddo endif np1 = n + 1 km1 = k - 1 kpkm2 = 2*km1 left = k ! zero out all entries of q lenq = n*(k+km1) do i=1,lenq q(i) = 0.0_wp enddo ! loop over i to construct the n interpolation equations do i=1,n xi = x(i) ilp1mx = min(i+k,np1) ! *** find left in the closed interval (i,i+k-1) such that ! t(left) <= x(i) < t(left+1) ! 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 return endif do if (xi<t(left+1)) go to 30 left = left + 1 if (left>=ilp1mx) exit enddo left = left - 1 if (xi>t(left+1)) 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 return endif ! *** 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 30 call dbspvn(t, k, k, 1, xi, left, bcoef, work, iwork, iflag) if (iflag/=0) 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 + (left-k)*(k+km1) do j=1,k jj = jj + kpkm2 q(jj) = bcoef(j) enddo enddo ! 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,n bcoef(i) = y(i) enddo call dbnslv(q, k+km1, n, km1, km1, bcoef) iflag = 0 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 endif endsubroutine dbintk