| [x,y,t,n,k,bcoef,q,work]=bintk(x,y,t,n,k,bcoef,q,work); |
function [x,y,t,n,k,bcoef,q,work]=bintk(x,y,t,n,k,bcoef,q,work);
%***BEGIN PROLOGUE BINTK
%***PURPOSE Compute the B-representation of a spline which interpolates
% given data.
%***LIBRARY SLATEC
%***CATEGORY E1A
%***TYPE SINGLE PRECISION (BINTK-S, DBINTK-D)
%***KEYWORDS B-SPLINE, DATA FITTING, INTERPOLATION
%***AUTHOR Amos, D. E., (SNLA)
%***DESCRIPTION
%
% Written by Carl de Boor and modified by D. E. Amos
%
% Abstract
%
% BINTK is the SPLINT routine of the reference.
%
% BINTK 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 BVALU.
% 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), 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 BNFAC (which
% constructs the triangular factorization for A and stores it
% again in Q), followed by a call to BNSLV (which then
% obtains the solution BCOEF by substitution). BNFAC 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) .LT. X(I)) .LT. T(I+K), 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.
%
% Description of Arguments
% Input
% X - vector of length N containing data point abscissa
% in strictly increasing order.
% Y - corresponding vector of length N containing data
% point ordinates.
% T - knot vector of length N+K
% since T(1),..,T(K) .LE. X(1) and T(N+1),..,T(N+K)
% .GE. X(N), this leaves only N-K knots (not nec-
% essarily X(I)) values) interior to (X(1),X(N))
% N - number of data points, N .GE. K
% K - order of the spline, K .GE. 1
%
% Output
% BCOEF - a vector of length N containing the B-spline
% coefficients
% 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 BNSLV (Q,2K-1,N,K-1,K-1,BCOEF)
% WORK - work vector of length 2*K
%
% Error Conditions
% Improper input is a fatal error
% Singular system of equations is a fatal error
%
%***REFERENCES D. E. Amos, Computation with splines and B-splines,
% Report SAND78-1968, Sandia Laboratories, March 1979.
% Carl de Boor, Package for calculating with B-splines,
% SIAM Journal on Numerical Analysis 14, 3 (June 1977),
% pp. 441-472.
% Carl de Boor, A Practical Guide to Splines, Applied
% Mathematics Series 27, Springer-Verlag, New York,
% 1978.
%***ROUTINES CALLED BNFAC, BNSLV, BSPVN, XERMSG
%***REVISION HISTORY (YYMMDD)
% 800901 DATE WRITTEN
% 890531 Changed all specific intrinsics to generic. (WRB)
% 890831 Modified array declarations. (WRB)
% 890831 REVISION DATE from Version 3.2
% 891214 Prologue converted to Version 4.0 format. (BAB)
% 900315 CALLs to XERROR changed to CALLs to XERMSG. (THJ)
% 900326 Removed duplicate information from DESCRIPTION section.
% (WRB)
% 920501 Reformatted the REFERENCES section. (WRB)
%***end PROLOGUE BINTK
%
persistent i iflag igo ilp1mx iwork j jj km1 kpkm2 left lenq np1 xi ;
if isempty(iflag), iflag=0; end;
if isempty(iwork), iwork=0; end;
if isempty(i), i=0; end;
if isempty(ilp1mx), ilp1mx=0; end;
if isempty(j), j=0; end;
if isempty(jj), jj=0; end;
if isempty(km1), km1=0; end;
if isempty(kpkm2), kpkm2=0; end;
if isempty(left), left=0; end;
if isempty(lenq), lenq=0; end;
if isempty(np1), np1=0; end;
if isempty(igo), igo=0; end;
bcoef_shape=size(bcoef);bcoef=reshape(bcoef,1,[]);
y_shape=size(y);y=reshape(y,1,[]);
q_shape=size(q);q=reshape(q,1,[]);
t_shape=size(t);t=reshape(t,1,[]);
x_shape=size(x);x=reshape(x,1,[]);
if isempty(xi), xi=0; end;
work_shape=size(work);work=reshape(work,1,[]);
% DIMENSION Q(2*K-1,N), T(N+K)
%***FIRST EXECUTABLE STATEMENT BINTK
if( k<1 )
xermsg('SLATEC','BINTK','K DOES NOT SATISFY K.GE.1',2,1);
bcoef_shape=zeros(bcoef_shape);bcoef_shape(:)=bcoef(1:numel(bcoef_shape));bcoef=bcoef_shape;
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
q_shape=zeros(q_shape);q_shape(:)=q(1:numel(q_shape));q=q_shape;
t_shape=zeros(t_shape);t_shape(:)=t(1:numel(t_shape));t=t_shape;
x_shape=zeros(x_shape);x_shape(:)=x(1:numel(x_shape));x=x_shape;
work_shape=zeros(work_shape);work_shape(:)=work(1:numel(work_shape));work=work_shape;
return;
elseif( n<k ) ;
xermsg('SLATEC','BINTK','N DOES NOT SATISFY N.GE.K',2,1);
bcoef_shape=zeros(bcoef_shape);bcoef_shape(:)=bcoef(1:numel(bcoef_shape));bcoef=bcoef_shape;
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
q_shape=zeros(q_shape);q_shape(:)=q(1:numel(q_shape));q=q_shape;
t_shape=zeros(t_shape);t_shape(:)=t(1:numel(t_shape));t=t_shape;
x_shape=zeros(x_shape);x_shape(:)=x(1:numel(x_shape));x=x_shape;
work_shape=zeros(work_shape);work_shape(:)=work(1:numel(work_shape));work=work_shape;
return;
else;
jj = fix(n - 1);
if( jj~=0 )
for i = 1 : jj;
if( x(i)>=x(i+1) )
xermsg('SLATEC','BINTK','X(I) DOES NOT SATISFY X(I).LT.X(I+1) FOR SOME I',2,1);
bcoef_shape=zeros(bcoef_shape);bcoef_shape(:)=bcoef(1:numel(bcoef_shape));bcoef=bcoef_shape;
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
q_shape=zeros(q_shape);q_shape(:)=q(1:numel(q_shape));q=q_shape;
t_shape=zeros(t_shape);t_shape(:)=t(1:numel(t_shape));t=t_shape;
x_shape=zeros(x_shape);x_shape(:)=x(1:numel(x_shape));x=x_shape;
work_shape=zeros(work_shape);work_shape(:)=work(1:numel(work_shape));work=work_shape;
return;
end;
end; i = fix(jj+1);
end;
np1 = fix(n + 1);
km1 = fix(k - 1);
kpkm2 = fix(2.*km1);
left = fix(k);
% ZERO OUT ALL ENTRIES OF Q
lenq = fix(n.*(k+km1));
for i = 1 : lenq;
q(i) = 0.0e0;
end; i = fix(lenq+1);
%
% *** LOOP OVER I TO CONSTRUCT THE N INTERPOLATION EQUATIONS
for i = 1 : n;
xi = x(i);
ilp1mx = fix(min(i+k,np1));
% *** FIND LEFT IN THE CLOSED INTERVAL (I,I+K-1) SUCH THAT
% T(LEFT) .LE. X(I) .LT. T(LEFT+1)
% MATRIX IS SINGULAR IF THIS IS NOT POSSIBLE
left = fix(max(left,i));
if( xi<t(left) )
xermsg('SLATEC','BINTK',['SOME ABSCISSA WAS NOT IN THE SUPPORT OF THE CORRESPONDING ','BASIS FUNCTION AND THE SYSTEM IS SINGULAR.'],2,1);
bcoef_shape=zeros(bcoef_shape);bcoef_shape(:)=bcoef(1:numel(bcoef_shape));bcoef=bcoef_shape;
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
q_shape=zeros(q_shape);q_shape(:)=q(1:numel(q_shape));q=q_shape;
t_shape=zeros(t_shape);t_shape(:)=t(1:numel(t_shape));t=t_shape;
x_shape=zeros(x_shape);x_shape(:)=x(1:numel(x_shape));x=x_shape;
work_shape=zeros(work_shape);work_shape(:)=work(1:numel(work_shape));work=work_shape;
return;
end;
igo=1;
while(igo==1);
igo=0;
if( xi>=t(left+1) )
left = fix(left + 1);
if( left<ilp1mx )
igo=1;
continue;
end;
left = fix(left - 1);
if( xi>t(left+1) )
xermsg('SLATEC','BINTK',['SOME ABSCISSA WAS NOT IN THE SUPPORT OF THE CORRESPONDING ','BASIS FUNCTION AND THE SYSTEM IS SINGULAR.'],2,1);
bcoef_shape=zeros(bcoef_shape);bcoef_shape(:)=bcoef(1:numel(bcoef_shape));bcoef=bcoef_shape;
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
q_shape=zeros(q_shape);q_shape(:)=q(1:numel(q_shape));q=q_shape;
t_shape=zeros(t_shape);t_shape(:)=t(1:numel(t_shape));t=t_shape;
x_shape=zeros(x_shape);x_shape(:)=x(1:numel(x_shape));x=x_shape;
work_shape=zeros(work_shape);work_shape(:)=work(1:numel(work_shape));work=work_shape;
return;
end;
end;
end;
% *** 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
k_orig=k; [t,k,dumvar3,dumvar4,xi,left,bcoef,work,iwork]=bspvn(t,k,k,1,xi,left,bcoef,work,iwork); k(dumvar3~=k_orig)=dumvar3(dumvar3~=k_orig);
% 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
% BNFAC). 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 = fix(i - left + 1 +(left-k).*(k+km1));
for j = 1 : k;
jj = fix(jj + kpkm2);
q(jj) = bcoef(j);
end; j = fix(k+1);
end; i = fix(n+1);
%
% ***OBTAIN FACTORIZATION OF A , STORED AGAIN IN Q.
km1_orig=km1; [q,dumvar2,n,km1,dumvar5,iflag]=bnfac(q,k+km1,n,km1,km1,iflag); km1(dumvar5~=km1_orig)=dumvar5(dumvar5~=km1_orig);
if( iflag==2 )
xermsg('SLATEC','BINTK',['THE SYSTEM OF SOLVER DETECTS A SINGULAR SYSTEM ALTHOUGH ','THE THEORETICAL CONDITIONS FOR A SOLUTION WERE SATISFIED.'],8,1);
bcoef_shape=zeros(bcoef_shape);bcoef_shape(:)=bcoef(1:numel(bcoef_shape));bcoef=bcoef_shape;
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
q_shape=zeros(q_shape);q_shape(:)=q(1:numel(q_shape));q=q_shape;
t_shape=zeros(t_shape);t_shape(:)=t(1:numel(t_shape));t=t_shape;
x_shape=zeros(x_shape);x_shape(:)=x(1:numel(x_shape));x=x_shape;
work_shape=zeros(work_shape);work_shape(:)=work(1:numel(work_shape));work=work_shape;
return;
else;
% *** SOLVE A*BCOEF = Y BY BACKSUBSTITUTION
for i = 1 : n;
bcoef(i) = y(i);
end; i = fix(n+1);
km1_orig=km1; [q,dumvar2,n,km1,dumvar5,bcoef]=bnslv(q,k+km1,n,km1,km1,bcoef); km1(dumvar5~=km1_orig)=dumvar5(dumvar5~=km1_orig);
bcoef_shape=zeros(bcoef_shape);bcoef_shape(:)=bcoef(1:numel(bcoef_shape));bcoef=bcoef_shape;
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
q_shape=zeros(q_shape);q_shape(:)=q(1:numel(q_shape));q=q_shape;
t_shape=zeros(t_shape);t_shape(:)=t(1:numel(t_shape));t=t_shape;
x_shape=zeros(x_shape);x_shape(:)=x(1:numel(x_shape));x=x_shape;
work_shape=zeros(work_shape);work_shape(:)=work(1:numel(work_shape));work=work_shape;
return;
end;
%
%
xermsg('SLATEC','BINTK',['SOME ABSCISSA WAS NOT IN THE SUPPORT OF THE CORRESPONDING ','BASIS FUNCTION AND THE SYSTEM IS SINGULAR.'],2,1);
bcoef_shape=zeros(bcoef_shape);bcoef_shape(:)=bcoef(1:numel(bcoef_shape));bcoef=bcoef_shape;
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
q_shape=zeros(q_shape);q_shape(:)=q(1:numel(q_shape));q=q_shape;
t_shape=zeros(t_shape);t_shape(:)=t(1:numel(t_shape));t=t_shape;
x_shape=zeros(x_shape);x_shape(:)=x(1:numel(x_shape));x=x_shape;
work_shape=zeros(work_shape);work_shape(:)=work(1:numel(work_shape));work=work_shape;
return;
end;
bcoef_shape=zeros(bcoef_shape);bcoef_shape(:)=bcoef(1:numel(bcoef_shape));bcoef=bcoef_shape;
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
q_shape=zeros(q_shape);q_shape(:)=q(1:numel(q_shape));q=q_shape;
t_shape=zeros(t_shape);t_shape(:)=t(1:numel(t_shape));t=t_shape;
x_shape=zeros(x_shape);x_shape(:)=x(1:numel(x_shape));x=x_shape;
work_shape=zeros(work_shape);work_shape(:)=work(1:numel(work_shape));work=work_shape;
end
%DECK BISECT
|
|
