subroutine dchud(r,ldr,p,x,z,ldz,nz,y,rho,c,s)
integer ldr,p,ldz,nz
double precision rho(1),c(1)
double precision r(ldr,1),x(1),z(ldz,1),y(1),s(1)
c
c dchud updates an augmented cholesky decomposition of the
c triangular part of an augmented qr decomposition. specifically,
c given an upper triangular matrix r of order p, a row vector
c x, a column vector z, and a scalar y, dchud determines a
c untiary matrix u and a scalar zeta such that
c
c
c (r z) (rr zz )
c u * ( ) = ( ) ,
c (x y) ( 0 zeta)
c
c where rr is upper triangular. if r and z have been
c obtained from the factorization of a least squares
c problem, then rr and zz are the factors corresponding to
c the problem with the observation (x,y) appended. in this
c case, if rho is the norm of the residual vector, then the
c norm of the residual vector of the updated problem is
c dsqrt(rho**2 + zeta**2). dchud will simultaneously update
c several triplets (z,y,rho).
c for a less terse description of what dchud does and how
c it may be applied, see the linpack guide.
c
c the matrix u is determined as the product u(p)*...*u(1),
c where u(i) is a rotation in the (i,p+1) plane of the
c form
c
c ( c(i) s(i) )
c ( ) .
c ( -s(i) c(i) )
c
c the rotations are chosen so that c(i) is double precision.
c
c on entry
c
c r double precision(ldr,p), where ldr .ge. p.
c r contains the upper triangular matrix
c that is to be updated. the part of r
c below the diagonal is not referenced.
c
c ldr integer.
c ldr is the leading dimension of the array r.
c
c p integer.
c p is the order of the matrix r.
c
c x double precision(p).
c x contains the row to be added to r. x is
c not altered by dchud.
c
c z double precision(ldz,nz), where ldz .ge. p.
c z is an array containing nz p-vectors to
c be updated with r.
c
c ldz integer.
c ldz is the leading dimension of the array z.
c
c nz integer.
c nz is the number of vectors to be updated
c nz may be zero, in which case z, y, and rho
c are not referenced.
c
c y double precision(nz).
c y contains the scalars for updating the vectors
c z. y is not altered by dchud.
c
c rho double precision(nz).
c rho contains the norms of the residual
c vectors that are to be updated. if rho(j)
c is negative, it is left unaltered.
c
c on return
c
c rc
c rho contain the updated quantities.
c z
c
c c double precision(p).
c c contains the cosines of the transforming
c rotations.
c
c s double precision(p).
c s contains the sines of the transforming
c rotations.
c
c linpack. this version dated 08/14/78 .
c g.w. stewart, university of maryland, argonne national lab.
c
c dchud uses the following functions and subroutines.
c
c extended blas drotg
c fortran dsqrt
c
integer i,j,jm1
double precision azeta,scale
double precision t,xj,zeta
c
c update r.
c
do 30 j = 1, p
xj = x(j)
c
c apply the previous rotations.
c
jm1 = j - 1
if (jm1 .lt. 1) go to 20
do 10 i = 1, jm1
t = c(i)*r(i,j) + s(i)*xj
xj = c(i)*xj - s(i)*r(i,j)
r(i,j) = t
10 continue
20 continue
c
c compute the next rotation.
c
call drotg(r(j,j),xj,c(j),s(j))
30 continue
c
c if required, update z and rho.
c
if (nz .lt. 1) go to 70
do 60 j = 1, nz
zeta = y(j)
do 40 i = 1, p
t = c(i)*z(i,j) + s(i)*zeta
zeta = c(i)*zeta - s(i)*z(i,j)
z(i,j) = t
40 continue
azeta = dabs(zeta)
if (azeta .eq. 0.0d0 .or. rho(j) .lt. 0.0d0) go to 50
scale = azeta + rho(j)
rho(j) = scale*dsqrt((azeta/scale)**2+(rho(j)/scale)**2)
50 continue
60 continue
70 continue
return
end