[r,ldr,p,x,z,ldz,nz,y,rho,c,s]=cchud(r,ldr,p,x,z,ldz,nz,y,rho,c,s);

function [r,ldr,p,x,z,ldz,nz,y,rho,c,s]=cchud(r,ldr,p,x,z,ldz,nz,y,rho,c,s);
%***BEGIN PROLOGUE  CCHUD
%***PURPOSE  Update an augmented Cholesky decomposition of the
%            triangular part of an augmented QR decomposition.
%***LIBRARY   SLATEC (LINPACK)
%***CATEGORY  D7B
%***TYPE      COMPLEX (SCHUD-S, DCHUD-D, CCHUD-C)
%***KEYWORDS  CHOLESKY DECOMPOSITION, LINEAR ALGEBRA, LINPACK, MATRIX,
%             UPDATE
%***AUTHOR  Stewart, G. W., (U. of Maryland)
%***DESCRIPTION
%
%     CCHUD updates an augmented Cholesky decomposition of the
%     triangular part of an augmented QR decomposition.  Specifically,
%     given an upper triangular matrix R of order P, a row vector
%     X, a column vector Z, and a scalar Y, CCHUD determines a
%     unitary matrix U and a scalar ZETA such that
%
%
%                              (R  Z)     (RR   ZZ )
%                         U  * (    )  =  (        ) ,
%                              (X  Y)     ( 0  ZETA)
%
%     where RR is upper triangular.  If R and Z have been
%     obtained from the factorization of a least squares
%     problem, then RR and ZZ are the factors corresponding to
%     the problem with the observation (X,Y) appended.  In this
%     case, if RHO is the norm of the residual vector, then the
%     norm of the residual vector of the updated problem is
%     SQRT(RHO**2 + ZETA**2).  CCHUD will simultaneously update
%     several triplets (Z,Y,RHO).
%
%     For a less terse description of what CCHUD does and how
%     it may be applied see the LINPACK Guide.
%
%     The matrix U is determined as the product U(P)*...*U(1),
%     where U(I) is a rotation in the (I,P+1) plane of the
%     form
%
%                       (     (CI)      S(I) )
%                       (                    ) .
%                       ( -CONJG(S(I))  (CI) )
%
%     The rotations are chosen so that C(I) is real.
%
%     On Entry
%
%         R      COMPLEX(LDR,P), where LDR .GE. P.
%                R contains the upper triangular matrix
%                that is to be updated.  The part of R
%                below the diagonal is not referenced.
%
%         LDR    INTEGER.
%                LDR is the leading dimension of the array R.
%
%         P      INTEGER.
%                P is the order of the matrix R.
%
%         X      COMPLEX(P).
%                X contains the row to be added to R.  X is
%                not altered by CCHUD.
%
%         Z      COMPLEX(LDZ,NZ), where LDZ .GE. P.
%                Z is an array containing NZ P-vectors to
%                be updated with R.
%
%         LDZ    INTEGER.
%                LDZ is the leading dimension of the array Z.
%
%         NZ     INTEGER.
%                NZ is the number of vectors to be updated
%                NZ may be zero, in which case Z, Y, and RHO
%                are not referenced.
%
%         Y      COMPLEX(NZ).
%                Y contains the scalars for updating the vectors
%                Z.  Y is not altered by CCHUD.
%
%         RHO    REAL(NZ).
%                RHO contains the norms of the residual
%                vectors that are to be updated.  If RHO(J)
%                is negative, it is left unaltered.
%
%     On Return
%
%         RC
%         RHO    contain the updated quantities.
%         Z
%
%         C      REAL(P).
%                C contains the cosines of the transforming
%                rotations.
%
%         S      COMPLEX(P).
%                S contains the sines of the transforming
%                rotations.
%
%***REFERENCES  J. J. Dongarra, J. R. Bunch, C. B. Moler, and G. W.
%                 Stewart, LINPACK Users' Guide, SIAM, 1979.
%***ROUTINES CALLED  CROTG
%***REVISION HISTORY  (YYMMDD)
%   780814  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)
%   900326  Removed duplicate information from DESCRIPTION section.
%           (WRB)
%   920501  Reformatted the REFERENCES section.  (WRB)
%***end PROLOGUE  CCHUD
persistent azeta i j jm1 scalemlv t xj zeta ; 

rho_shape=size(rho);rho=reshape(rho,1,[]);
c_shape=size(c);c=reshape(c,1,[]);
r_shape=size(r);r=reshape([r(:).',zeros(1,ceil(numel(r)./prod([ldr])).*prod([ldr])-numel(r))],ldr,[]);
x_shape=size(x);x=reshape(x,1,[]);
z_shape=size(z);z=reshape([z(:).',zeros(1,ceil(numel(z)./prod([ldz])).*prod([ldz])-numel(z))],ldz,[]);
y_shape=size(y);y=reshape(y,1,[]);
s_shape=size(s);s=reshape(s,1,[]);
%
if isempty(i), i=0; end;
if isempty(j), j=0; end;
if isempty(jm1), jm1=0; end;
if isempty(azeta), azeta=0; end;
if isempty(scalemlv), scalemlv=0; end;
if isempty(t), t=0; end;
if isempty(xj), xj=0; end;
if isempty(zeta), zeta=0; end;
%
%     UPDATE R.
%
%***FIRST EXECUTABLE STATEMENT  CCHUD
for j = 1 : p;
xj = x(j);
%
%        APPLY THE PREVIOUS ROTATIONS.
%
jm1 = fix(j - 1);
if( jm1>=1 )
for i = 1 : jm1;
t = c(i).*r(i,j) + s(i).*xj;
xj = c(i).*xj - conj(s(i)).*r(i,j);
r(i,j) = t;
end; i = fix(jm1+1);
end;
%
%        COMPUTE THE NEXT ROTATION.
%
[r(j,j),xj,c(j),s(j)]=crotg(r(j,j),xj,c(j),s(j));
end; j = fix(p+1);
%
%     IF REQUIRED, UPDATE Z AND RHO.
%
if( nz>=1 )
for j = 1 : nz;
zeta = y(j);
for i = 1 : p;
t = c(i).*z(i,j) + s(i).*zeta;
zeta = c(i).*zeta - conj(s(i)).*z(i,j);
z(i,j) = t;
end; i = fix(p+1);
azeta = abs(zeta);
if( azeta~=0.0e0 && rho(j)>=0.0e0 )
scalemlv = azeta + rho(j);
rho(j) = scalemlv.*sqrt((azeta./scalemlv).^2+(rho(j)./scalemlv).^2);
end;
end; j = fix(nz+1);
end;
rho_shape=zeros(rho_shape);rho_shape(:)=rho(1:numel(rho_shape));rho=rho_shape;
c_shape=zeros(c_shape);c_shape(:)=c(1:numel(c_shape));c=c_shape;
r_shape=zeros(r_shape);r_shape(:)=r(1:numel(r_shape));r=r_shape;
x_shape=zeros(x_shape);x_shape(:)=x(1:numel(x_shape));x=x_shape;
z_shape=zeros(z_shape);z_shape(:)=z(1:numel(z_shape));z=z_shape;
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
s_shape=zeros(s_shape);s_shape(:)=s(1:numel(s_shape));s=s_shape;
end
%DECK CCMPB