[ndata,xdata,ydata,sddata,nord,nbkpt,bkptin,mdein,mdeout,coeff,bf,xtemp,ptemp,bkpt,g,mdg,w,mdw,lw]=defcmn(ndata,xdata,ydata,sddata,nord,nbkpt,bkptin,mdein,mdeout,coeff,bf,xtemp,ptemp,bkpt,g,mdg,w,mdw,lw);

function [ndata,xdata,ydata,sddata,nord,nbkpt,bkptin,mdein,mdeout,coeff,bf,xtemp,ptemp,bkpt,g,mdg,w,mdw,lw]=defcmn(ndata,xdata,ydata,sddata,nord,nbkpt,bkptin,mdein,mdeout,coeff,bf,xtemp,ptemp,bkpt,g,mdg,w,mdw,lw);
%***BEGIN PROLOGUE  DEFCMN
%***SUBSIDIARY
%***PURPOSE  Subsidiary to DEFC
%***LIBRARY   SLATEC
%***TYPE      doubleprecision (EFCMN-S, DEFCMN-D)
%***AUTHOR  Hanson, R. J., (SNLA)
%***DESCRIPTION
%
%     This is a companion subprogram to DEFC( ).
%     This subprogram does weighted least squares fitting of data by
%     B-spline curves.
%     The documentation for DEFC( ) has complete usage instructions.
%
%***SEE ALSO  DEFC
%***ROUTINES CALLED  DBNDAC, DBNDSL, DCOPY, DFSPVN, DSCAL, DSORT, XERMSG
%***REVISION HISTORY  (YYMMDD)
%   800801  DATE WRITTEN
%   890531  Changed all specific intrinsics to generic.  (WRB)
%   890618  Completely restructured and extensively revised (WRBRWC)
%   891214  Prologue converted to Version 4.0 format.  (BAB)
%   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
%   900328  Added TYPE section.  (WRB)
%   900510  Convert XERRWV calls to XERMSG calls.  (RWC)
%   900604  DP version created from SP version.  (RWC)
%***end PROLOGUE  DEFCMN
persistent dummy i idata ileft intseq ip ir irow l mt n nb nordm1 nordp1 np1 rnorm xern1 xern2 xmax xmin xval ; 

bf_shape=size(bf);bf=reshape([bf(:).',zeros(1,ceil(numel(bf)./prod([nord])).*prod([nord])-numel(bf))],nord,[]);
bkpt_shape=size(bkpt);bkpt=reshape(bkpt,1,[]);
bkptin_shape=size(bkptin);bkptin=reshape(bkptin,1,[]);
coeff_shape=size(coeff);coeff=reshape(coeff,1,[]);
g_shape=size(g);g=reshape([g(:).',zeros(1,ceil(numel(g)./prod([mdg])).*prod([mdg])-numel(g))],mdg,[]);
ptemp_shape=size(ptemp);ptemp=reshape(ptemp,1,[]);
sddata_shape=size(sddata);sddata=reshape(sddata,1,[]);
w_shape=size(w);w=reshape([w(:).',zeros(1,ceil(numel(w)./prod([mdw])).*prod([mdw])-numel(w))],mdw,[]);
xdata_shape=size(xdata);xdata=reshape(xdata,1,[]);
xtemp_shape=size(xtemp);xtemp=reshape(xtemp,1,[]);
ydata_shape=size(ydata);ydata=reshape(ydata,1,[]);
%
%
if isempty(dummy), dummy=0; end;
if isempty(rnorm), rnorm=0; end;
if isempty(xmax), xmax=0; end;
if isempty(xmin), xmin=0; end;
if isempty(xval), xval=0; end;
if isempty(i), i=0; end;
if isempty(idata), idata=0; end;
if isempty(ileft), ileft=0; end;
if isempty(intseq), intseq=0; end;
if isempty(ip), ip=0; end;
if isempty(ir), ir=0; end;
if isempty(irow), irow=0; end;
if isempty(l), l=0; end;
if isempty(mt), mt=0; end;
if isempty(n), n=0; end;
if isempty(nb), nb=0; end;
if isempty(nordm1), nordm1=0; end;
if isempty(nordp1), nordp1=0; end;
if isempty(np1), np1=0; end;
if isempty(xern1), xern1=repmat(' ',1,8); end;
if isempty(xern2), xern2=repmat(' ',1,8); end;
%
%***FIRST EXECUTABLE STATEMENT  DEFCMN
%
%     Initialize variables and analyze input.
%
n = fix(nbkpt - nord);
np1 = fix(n + 1);
%
%     Initially set all output coefficients to zero.
%
[n,dumvar2,dumvar3,coeff]=dcopy(n,0.0d0,0,coeff,1);
mdeout = -1;
if( nord<1 || nord>20 )
xermsg('SLATEC','DEFCMN','IN DEFC, THE ORDER OF THE B-SPLINE MUST BE 1 THRU 20.',3,1);
bf_shape=zeros(bf_shape);bf_shape(:)=bf(1:numel(bf_shape));bf=bf_shape;
bkpt_shape=zeros(bkpt_shape);bkpt_shape(:)=bkpt(1:numel(bkpt_shape));bkpt=bkpt_shape;
bkptin_shape=zeros(bkptin_shape);bkptin_shape(:)=bkptin(1:numel(bkptin_shape));bkptin=bkptin_shape;
coeff_shape=zeros(coeff_shape);coeff_shape(:)=coeff(1:numel(coeff_shape));coeff=coeff_shape;
g_shape=zeros(g_shape);g_shape(:)=g(1:numel(g_shape));g=g_shape;
ptemp_shape=zeros(ptemp_shape);ptemp_shape(:)=ptemp(1:numel(ptemp_shape));ptemp=ptemp_shape;
sddata_shape=zeros(sddata_shape);sddata_shape(:)=sddata(1:numel(sddata_shape));sddata=sddata_shape;
w_shape=zeros(w_shape);w_shape(:)=w(1:numel(w_shape));w=w_shape;
xdata_shape=zeros(xdata_shape);xdata_shape(:)=xdata(1:numel(xdata_shape));xdata=xdata_shape;
xtemp_shape=zeros(xtemp_shape);xtemp_shape(:)=xtemp(1:numel(xtemp_shape));xtemp=xtemp_shape;
ydata_shape=zeros(ydata_shape);ydata_shape(:)=ydata(1:numel(ydata_shape));ydata=ydata_shape;
return;
end;
%
if( nbkpt<2.*nord )
xermsg('SLATEC','DEFCMN',['IN DEFC, THE NUMBER OF KNOTS MUST BE AT LEAST TWICE ','THE B-SPLINE ORDER.'],4,1);
bf_shape=zeros(bf_shape);bf_shape(:)=bf(1:numel(bf_shape));bf=bf_shape;
bkpt_shape=zeros(bkpt_shape);bkpt_shape(:)=bkpt(1:numel(bkpt_shape));bkpt=bkpt_shape;
bkptin_shape=zeros(bkptin_shape);bkptin_shape(:)=bkptin(1:numel(bkptin_shape));bkptin=bkptin_shape;
coeff_shape=zeros(coeff_shape);coeff_shape(:)=coeff(1:numel(coeff_shape));coeff=coeff_shape;
g_shape=zeros(g_shape);g_shape(:)=g(1:numel(g_shape));g=g_shape;
ptemp_shape=zeros(ptemp_shape);ptemp_shape(:)=ptemp(1:numel(ptemp_shape));ptemp=ptemp_shape;
sddata_shape=zeros(sddata_shape);sddata_shape(:)=sddata(1:numel(sddata_shape));sddata=sddata_shape;
w_shape=zeros(w_shape);w_shape(:)=w(1:numel(w_shape));w=w_shape;
xdata_shape=zeros(xdata_shape);xdata_shape(:)=xdata(1:numel(xdata_shape));xdata=xdata_shape;
xtemp_shape=zeros(xtemp_shape);xtemp_shape(:)=xtemp(1:numel(xtemp_shape));xtemp=xtemp_shape;
ydata_shape=zeros(ydata_shape);ydata_shape(:)=ydata(1:numel(ydata_shape));ydata=ydata_shape;
return;
end;
%
if( ndata<0 )
xermsg('SLATEC','DEFCMN','IN DEFC, THE NUMBER OF DATA POINTS MUST BE NONNEGATIVE.',5,1);
bf_shape=zeros(bf_shape);bf_shape(:)=bf(1:numel(bf_shape));bf=bf_shape;
bkpt_shape=zeros(bkpt_shape);bkpt_shape(:)=bkpt(1:numel(bkpt_shape));bkpt=bkpt_shape;
bkptin_shape=zeros(bkptin_shape);bkptin_shape(:)=bkptin(1:numel(bkptin_shape));bkptin=bkptin_shape;
coeff_shape=zeros(coeff_shape);coeff_shape(:)=coeff(1:numel(coeff_shape));coeff=coeff_shape;
g_shape=zeros(g_shape);g_shape(:)=g(1:numel(g_shape));g=g_shape;
ptemp_shape=zeros(ptemp_shape);ptemp_shape(:)=ptemp(1:numel(ptemp_shape));ptemp=ptemp_shape;
sddata_shape=zeros(sddata_shape);sddata_shape(:)=sddata(1:numel(sddata_shape));sddata=sddata_shape;
w_shape=zeros(w_shape);w_shape(:)=w(1:numel(w_shape));w=w_shape;
xdata_shape=zeros(xdata_shape);xdata_shape(:)=xdata(1:numel(xdata_shape));xdata=xdata_shape;
xtemp_shape=zeros(xtemp_shape);xtemp_shape(:)=xtemp(1:numel(xtemp_shape));xtemp=xtemp_shape;
ydata_shape=zeros(ydata_shape);ydata_shape(:)=ydata(1:numel(ydata_shape));ydata=ydata_shape;
return;
end;
%
nb =fix((nbkpt-nord+3).*(nord+1) +(nbkpt+1).*(nord+1)+ 2.*max(nbkpt,ndata) + nbkpt + nord.^2);
if( lw<nb )
xern1=sprintf(['%8i'], nb);
xern2=sprintf(['%8i'], lw);
xermsg('SLATEC','DEFCMN',['IN DEFC, INSUFFICIENT STORAGE FOR W(*).  CHECK FORMULA ',['THAT READS LW.GE. ... .  NEED = ',[xern1,[' GIVEN = ',xern2]]]],6,1);
mdeout = -1;
bf_shape=zeros(bf_shape);bf_shape(:)=bf(1:numel(bf_shape));bf=bf_shape;
bkpt_shape=zeros(bkpt_shape);bkpt_shape(:)=bkpt(1:numel(bkpt_shape));bkpt=bkpt_shape;
bkptin_shape=zeros(bkptin_shape);bkptin_shape(:)=bkptin(1:numel(bkptin_shape));bkptin=bkptin_shape;
coeff_shape=zeros(coeff_shape);coeff_shape(:)=coeff(1:numel(coeff_shape));coeff=coeff_shape;
g_shape=zeros(g_shape);g_shape(:)=g(1:numel(g_shape));g=g_shape;
ptemp_shape=zeros(ptemp_shape);ptemp_shape(:)=ptemp(1:numel(ptemp_shape));ptemp=ptemp_shape;
sddata_shape=zeros(sddata_shape);sddata_shape(:)=sddata(1:numel(sddata_shape));sddata=sddata_shape;
w_shape=zeros(w_shape);w_shape(:)=w(1:numel(w_shape));w=w_shape;
xdata_shape=zeros(xdata_shape);xdata_shape(:)=xdata(1:numel(xdata_shape));xdata=xdata_shape;
xtemp_shape=zeros(xtemp_shape);xtemp_shape(:)=xtemp(1:numel(xtemp_shape));xtemp=xtemp_shape;
ydata_shape=zeros(ydata_shape);ydata_shape(:)=ydata(1:numel(ydata_shape));ydata=ydata_shape;
return;
end;
%
if( mdein~=1 && mdein~=2 )
xermsg('SLATEC','DEFCMN','IN DEFC, INPUT VALUE OF MDEIN MUST BE 1-2.',7,1);
bf_shape=zeros(bf_shape);bf_shape(:)=bf(1:numel(bf_shape));bf=bf_shape;
bkpt_shape=zeros(bkpt_shape);bkpt_shape(:)=bkpt(1:numel(bkpt_shape));bkpt=bkpt_shape;
bkptin_shape=zeros(bkptin_shape);bkptin_shape(:)=bkptin(1:numel(bkptin_shape));bkptin=bkptin_shape;
coeff_shape=zeros(coeff_shape);coeff_shape(:)=coeff(1:numel(coeff_shape));coeff=coeff_shape;
g_shape=zeros(g_shape);g_shape(:)=g(1:numel(g_shape));g=g_shape;
ptemp_shape=zeros(ptemp_shape);ptemp_shape(:)=ptemp(1:numel(ptemp_shape));ptemp=ptemp_shape;
sddata_shape=zeros(sddata_shape);sddata_shape(:)=sddata(1:numel(sddata_shape));sddata=sddata_shape;
w_shape=zeros(w_shape);w_shape(:)=w(1:numel(w_shape));w=w_shape;
xdata_shape=zeros(xdata_shape);xdata_shape(:)=xdata(1:numel(xdata_shape));xdata=xdata_shape;
xtemp_shape=zeros(xtemp_shape);xtemp_shape(:)=xtemp(1:numel(xtemp_shape));xtemp=xtemp_shape;
ydata_shape=zeros(ydata_shape);ydata_shape(:)=ydata(1:numel(ydata_shape));ydata=ydata_shape;
return;
end;
%
%     Sort the breakpoints.
%
[nbkpt,bkptin,dumvar3,bkpt]=dcopy(nbkpt,bkptin,1,bkpt,1);
[bkpt,dummy,nbkpt]=dsort(bkpt,dummy,nbkpt,1);
%
%     Save interval containing knots.
%
xmin = bkpt(nord);
xmax = bkpt(np1);
nordm1 = fix(nord - 1);
nordp1 = fix(nord + 1);
%
%     Process least squares equations.
%
%     Sort data and an array of pointers.
%
[ndata,xdata,dumvar3,xtemp]=dcopy(ndata,xdata,1,xtemp,1);
for i = 1 : ndata;
ptemp(i) = i;
end; i = fix(ndata+1);
%
if( ndata>0 )
[xtemp,ptemp,ndata]=dsort(xtemp,ptemp,ndata,2);
xmin = min(xmin,xtemp(1));
xmax = max(xmax,xtemp(ndata));
end;
%
%     Fix breakpoint array if needed. This should only involve very
%     minor differences with the input array of breakpoints.
%
for i = 1 : nord;
bkpt(i) = min(bkpt(i),xmin);
end; i = fix(nord+1);
%
for i = np1 : nbkpt;
bkpt(i) = max(bkpt(i),xmax);
end; i = fix(nbkpt+1);
%
%     Initialize parameters of banded matrix processor, DBNDAC( ).
%
mt = 0;
ip = 1;
ir = 1;
ileft = fix(nord);
intseq = 1;
for idata = 1 : ndata;
%
%        Sorted indices are in PTEMP(*).
%
l = fix(ptemp(idata));
xval = xdata(l);
%
%        When interval changes, process equations in the last block.
%
if( xval>=bkpt(ileft+1) )
[g,mdg,nord,ip,ir,mt]=dbndac(g,mdg,nord,ip,ir,mt,ileft-nordm1);
mt = 0;
%
%           Move pointer up to have BKPT(ILEFT).LE.XVAL, ILEFT.LE.N.
%
for ileft = ileft : n;
if( xval<bkpt(ileft+1) )
break;
end;
if( mdein==2 )
%
%                 Data is being sequentially accumulated.
%                 Transfer previously accumulated rows from W(*,*) to
%                 G(*,*) and process them.
%
[nordp1,w(sub2ind(size(w),intseq,1):end),mdw,g(sub2ind(size(g),ir,1):end),mdg]=dcopy(nordp1,w(sub2ind(size(w),intseq,1):end),mdw,g(sub2ind(size(g),ir,1):end),mdg);
[g,mdg,nord,ip,ir,dumvar6,intseq]=dbndac(g,mdg,nord,ip,ir,1,intseq);
intseq = fix(intseq + 1);
end;
end;
end;
%
%        Obtain B-spline function value.
%
[bkpt,nord,dumvar3,xval,ileft,bf]=dfspvn(bkpt,nord,1,xval,ileft,bf);
%
%        Move row into place.
%
irow = fix(ir + mt);
mt = fix(mt + 1);
[nord,bf,dumvar3,g(sub2ind(size(g),irow,1):end),mdg]=dcopy(nord,bf,1,g(sub2ind(size(g),irow,1):end),mdg);
g(irow,nordp1) = ydata(l);
%
%        Scale data if uncertainty is nonzero.
%
if( sddata(l)~=0.0d0 )
[nordp1,dumvar2,g(sub2ind(size(g),irow,1):end),mdg]=dscal(nordp1,1.0d0./sddata(l),g(sub2ind(size(g),irow,1):end),mdg);
end;
%
%        When staging work area is exhausted, process rows.
%
if( irow==mdg-1 )
[g,mdg,nord,ip,ir,mt]=dbndac(g,mdg,nord,ip,ir,mt,ileft-nordm1);
mt = 0;
end;
end;
%
%     Process last block of equations.
%
[g,mdg,nord,ip,ir,mt]=dbndac(g,mdg,nord,ip,ir,mt,ileft-nordm1);
%
%     Finish processing any previously accumulated rows from W(*,*)
%     to G(*,*).
%
if( mdein==2 )
for i = intseq : np1;
[nordp1,w(sub2ind(size(w),i,1):end),mdw,g(sub2ind(size(g),ir,1):end),mdg]=dcopy(nordp1,w(sub2ind(size(w),i,1):end),mdw,g(sub2ind(size(g),ir,1):end),mdg);
[g,mdg,nord,ip,ir]=dbndac(g,mdg,nord,ip,ir,1,min(n,i));
end; i = fix(np1+1);
end;
%
%     Last call to adjust block positioning.
%
[nordp1,dumvar2,dumvar3,g(sub2ind(size(g),ir,1):end),mdg]=dcopy(nordp1,0.0d0,0,g(sub2ind(size(g),ir,1):end),mdg);
[g,mdg,nord,ip,ir,dumvar6,np1]=dbndac(g,mdg,nord,ip,ir,1,np1);
%
%     Transfer accumulated rows from G(*,*) to W(*,*) for
%     possible later sequential accumulation.
%
for i = 1 : np1;
[nordp1,g(sub2ind(size(g),i,1):end),mdg,w(sub2ind(size(w),i,1):end),mdw]=dcopy(nordp1,g(sub2ind(size(g),i,1):end),mdg,w(sub2ind(size(w),i,1):end),mdw);
end; i = fix(np1+1);
%
%     Solve for coefficients when possible.
%
for i = 1 : n;
if( g(i,1)==0.0d0 )
mdeout = 2;
bf_shape=zeros(bf_shape);bf_shape(:)=bf(1:numel(bf_shape));bf=bf_shape;
bkpt_shape=zeros(bkpt_shape);bkpt_shape(:)=bkpt(1:numel(bkpt_shape));bkpt=bkpt_shape;
bkptin_shape=zeros(bkptin_shape);bkptin_shape(:)=bkptin(1:numel(bkptin_shape));bkptin=bkptin_shape;
coeff_shape=zeros(coeff_shape);coeff_shape(:)=coeff(1:numel(coeff_shape));coeff=coeff_shape;
g_shape=zeros(g_shape);g_shape(:)=g(1:numel(g_shape));g=g_shape;
ptemp_shape=zeros(ptemp_shape);ptemp_shape(:)=ptemp(1:numel(ptemp_shape));ptemp=ptemp_shape;
sddata_shape=zeros(sddata_shape);sddata_shape(:)=sddata(1:numel(sddata_shape));sddata=sddata_shape;
w_shape=zeros(w_shape);w_shape(:)=w(1:numel(w_shape));w=w_shape;
xdata_shape=zeros(xdata_shape);xdata_shape(:)=xdata(1:numel(xdata_shape));xdata=xdata_shape;
xtemp_shape=zeros(xtemp_shape);xtemp_shape(:)=xtemp(1:numel(xtemp_shape));xtemp=xtemp_shape;
ydata_shape=zeros(ydata_shape);ydata_shape(:)=ydata(1:numel(ydata_shape));ydata=ydata_shape;
return;
end;
end; i = fix(n+1);
%
%     All the diagonal terms in the accumulated triangular
%     matrix are nonzero.  The solution can be computed but
%     it may be unsuitable for further use due to poor
%     conditioning or the lack of constraints.  No checking
%     for either of these is done here.
%
[dumvar1,g,mdg,nord,ip,ir,coeff,n,rnorm]=dbndsl(1,g,mdg,nord,ip,ir,coeff,n,rnorm);
mdeout = 1;
bf_shape=zeros(bf_shape);bf_shape(:)=bf(1:numel(bf_shape));bf=bf_shape;
bkpt_shape=zeros(bkpt_shape);bkpt_shape(:)=bkpt(1:numel(bkpt_shape));bkpt=bkpt_shape;
bkptin_shape=zeros(bkptin_shape);bkptin_shape(:)=bkptin(1:numel(bkptin_shape));bkptin=bkptin_shape;
coeff_shape=zeros(coeff_shape);coeff_shape(:)=coeff(1:numel(coeff_shape));coeff=coeff_shape;
g_shape=zeros(g_shape);g_shape(:)=g(1:numel(g_shape));g=g_shape;
ptemp_shape=zeros(ptemp_shape);ptemp_shape(:)=ptemp(1:numel(ptemp_shape));ptemp=ptemp_shape;
sddata_shape=zeros(sddata_shape);sddata_shape(:)=sddata(1:numel(sddata_shape));sddata=sddata_shape;
w_shape=zeros(w_shape);w_shape(:)=w(1:numel(w_shape));w=w_shape;
xdata_shape=zeros(xdata_shape);xdata_shape(:)=xdata(1:numel(xdata_shape));xdata=xdata_shape;
xtemp_shape=zeros(xtemp_shape);xtemp_shape(:)=xtemp(1:numel(xtemp_shape));xtemp=xtemp_shape;
ydata_shape=zeros(ydata_shape);ydata_shape(:)=ydata(1:numel(ydata_shape));ydata=ydata_shape;
end
%DECK DEFE4