[dcvresult,xval,ndata,nconst,nord,nbkpt,bkpt,w]=dcv(xval,ndata,nconst,nord,nbkpt,bkpt,w);

function [dcvresult,xval,ndata,nconst,nord,nbkpt,bkpt,w]=dcv(xval,ndata,nconst,nord,nbkpt,bkpt,w);
dcvresult=[];
persistent i ileft ip is last mdg mdw n v zero ; 

;
%***BEGIN PROLOGUE  DCV
%***PURPOSE  Evaluate the variance function of the curve obtained
%            by the constrained B-spline fitting subprogram DFC.
%***LIBRARY   SLATEC
%***CATEGORY  L7A3
%***TYPE      doubleprecision (CV-S, DCV-D)
%***KEYWORDS  ANALYSIS OF COVARIANCE, B-SPLINE,
%             CONSTRAINED LEAST SQUARES, CURVE FITTING
%***AUTHOR  Hanson, R. J., (SNLA)
%***DESCRIPTION
%
%     DCV( ) is a companion function subprogram for DFC( ).  The
%     documentation for DFC( ) has complete usage instructions.
%
%     DCV( ) is used to evaluate the variance function of the curve
%     obtained by the constrained B-spline fitting subprogram, DFC( ).
%     The variance function defines the square of the probable error
%     of the fitted curve at any point, XVAL.  One can use the square
%     root of this variance function to determine a probable error band
%     around the fitted curve.
%
%     DCV( ) is used after a call to DFC( ).  MODE, an input variable to
%     DFC( ), is used to indicate if the variance function is desired.
%     In order to use DCV( ), MODE must equal 2 or 4 on input to DFC( ).
%     MODE is also used as an output flag from DFC( ).  Check to make
%     sure that MODE = 0 after calling DFC( ), indicating a successful
%     constrained curve fit.  The array SDDATA, as input to DFC( ), must
%     also be defined with the standard deviation or uncertainty of the
%     Y values to use DCV( ).
%
%     To evaluate the variance function after calling DFC( ) as stated
%     above, use DCV( ) as shown here
%
%          VAR=DCV(XVAL,NDATA,NCONST,NORD,NBKPT,BKPT,W)
%
%     The variance function is given by
%
%      VAR=(transpose of B(XVAL))*C*B(XVAL)/DBLE(MAX(NDATA-N,1))
%
%     where N = NBKPT - NORD.
%
%     The vector B(XVAL) is the B-spline basis function values at
%     X=XVAL.  The covariance matrix, C, of the solution coefficients
%     accounts only for the least squares equations and the explicitly
%     stated equality constraints.  This fact must be considered when
%     interpreting the variance function from a data fitting problem
%     that has inequality constraints on the fitted curve.
%
%     All the variables in the calling sequence for DCV( ) are used in
%     DFC( ) except the variable XVAL.  Do not change the values of
%     these variables between the call to DFC( ) and the use of DCV( ).
%
%     The following is a brief description of the variables
%
%     XVAL    The point where the variance is desired, a double
%             precision variable.
%
%     NDATA   The number of discrete (X,Y) pairs for which DFC( )
%             calculated a piece-wise polynomial curve.
%
%     NCONST  The number of conditions that constrained the B-spline in
%             DFC( ).
%
%     NORD    The order of the B-spline used in DFC( ).
%             The value of NORD must satisfy 1 < NORD < 20 .
%
%             (The order of the spline is one more than the degree of
%             the piece-wise polynomial defined on each interval.  This
%             is consistent with the B-spline package convention.  For
%             example, NORD=4 when we are using piece-wise cubics.)
%
%     NBKPT   The number of knots in the array BKPT(*).
%             The value of NBKPT must satisfy NBKPT .GE. 2*NORD.
%
%     BKPT(*) The doubleprecision array of knots.  Normally the problem
%             data interval will be included between the limits
%             BKPT(NORD) and BKPT(NBKPT-NORD+1).  The additional end
%             knots BKPT(I),I=1,...,NORD-1 and I=NBKPT-NORD+2,...,NBKPT,
%             are required by DFC( ) to compute the functions used to
%             fit the data.
%
%     W(*)    doubleprecision work array as used in DFC( ).  See DFC( )
%             for the required length of W(*).  The contents of W(*)
%             must not be modified by the user if the variance function
%             is desired.
%
%***REFERENCES  R. J. Hanson, Constrained least squares curve fitting
%                 to discrete data using B-splines, a users guide,
%                 Report SAND78-1291, Sandia Laboratories, December
%                 1978.
%***ROUTINES CALLED  DDOT, DFSPVN
%***REVISION HISTORY  (YYMMDD)
%   780801  DATE WRITTEN
%   890531  Changed all specific intrinsics to generic.  (WRB)
%   890831  Modified array declarations.  (WRB)
%   890911  Removed unnecessary intrinsics.  (WRB)
%   891006  Cosmetic changes to prologue.  (WRB)
%   891006  REVISION DATE from Version 3.2
%   891214  Prologue converted to Version 4.0 format.  (BAB)
%   920501  Reformatted the REFERENCES section.  (WRB)
%***end PROLOGUE  DCV
if isempty(i), i=0; end;
if isempty(ileft), ileft=0; end;
if isempty(ip), ip=0; end;
if isempty(is), is=0; end;
if isempty(last), last=0; end;
if isempty(mdg), mdg=0; end;
if isempty(mdw), mdw=0; end;
if isempty(n), n=0; end;
if isempty(v), v=zeros(1,40); end;
if isempty(zero), zero=0; end;
bkpt_shape=size(bkpt);bkpt=reshape(bkpt,1,[]);
w_shape=size(w);w=reshape(w,1,[]);
%***FIRST EXECUTABLE STATEMENT  DCV
zero = 0.0d0;
mdg = fix(nbkpt - nord + 3);
mdw = fix(nbkpt - nord + 1 + nconst);
is = fix(mdg.*(nord+1) + 2.*max(ndata,nbkpt) + nbkpt + nord.^2);
last = fix(nbkpt - nord + 1);
ileft = fix(nord);
while( xval>=bkpt(ileft+1) & ileft<last-1 );
ileft = fix(ileft + 1);
end;
[bkpt,nord,dumvar3,xval,ileft,v(sub2ind(size(v),max(nord+1,1)):end)]=dfspvn(bkpt,nord,1,xval,ileft,v(sub2ind(size(v),max(nord+1,1)):end));
ileft = fix(ileft - nord + 1);
ip = fix(mdw.*(ileft-1) + ileft + is);
n = fix(nbkpt - nord);
for i = 1 : nord;
[v(i) ,nord,w(sub2ind(size(w),max(ip,1)):end),dumvar4,v(sub2ind(size(v),max(nord+1,1)):end)]=ddot(nord,w(sub2ind(size(w),max(ip,1)):end),1,v(sub2ind(size(v),max(nord+1,1)):end),1);
ip = fix(ip + mdw);
end; i = fix(nord+1);
dcvresult = max(ddot(nord,v,1,v(sub2ind(size(v),max(nord+1,1)):end),1),zero);
%
%     SCALE THE VARIANCE SO IT IS AN UNBIASED ESTIMATE.
dcvresult = dcvresult./max(ndata-n,1);
bkpt_shape=zeros(bkpt_shape);bkpt_shape(:)=bkpt(1:numel(bkpt_shape));bkpt=bkpt_shape;
w_shape=zeros(w_shape);w_shape(:)=w(1:numel(w_shape));w=w_shape;
csnil=dbstack(1); csnil=csnil(1).name(1)~='@';
if csnil&&~isempty(inputname(1)), assignin('caller','FUntemp',xval); evalin('caller',[inputname(1),'=FUntemp;']); end
if csnil&&~isempty(inputname(7)), assignin('caller','FUntemp',w); evalin('caller',[inputname(7),'=FUntemp;']); end
if csnil&&~isempty(inputname(4)), assignin('caller','FUntemp',nord); evalin('caller',[inputname(4),'=FUntemp;']); end
if csnil&&~isempty(inputname(2)), assignin('caller','FUntemp',ndata); evalin('caller',[inputname(2),'=FUntemp;']); end
if csnil&&~isempty(inputname(3)), assignin('caller','FUntemp',nconst); evalin('caller',[inputname(3),'=FUntemp;']); end
if csnil&&~isempty(inputname(5)), assignin('caller','FUntemp',nbkpt); evalin('caller',[inputname(5),'=FUntemp;']); end
if csnil&&~isempty(inputname(6)), assignin('caller','FUntemp',bkpt); evalin('caller',[inputname(6),'=FUntemp;']); end
end %function dcv
%DECK DDAINI