subroutine dchdd(r,ldr,p,x,z,ldz,nz,y,rho,c,s,info)
      integer ldr,p,ldz,nz,info
      double precision r(ldr,1),x(1),z(ldz,1),y(1),s(1)
      double precision rho(1),c(1)
c
c     dchdd downdates an augmented cholesky decomposition or the
c     triangular factor of an augmented qr decomposition.
c     specifically, given an upper triangular matrix r of order p,  a
c     row vector x, a column vector z, and a scalar y, dchdd
c     determineds a orthogonal matrix u and a scalar zeta such that
c
c                        (r   z )     (rr  zz)
c                    u * (      )  =  (      ) ,
c                        (0 zeta)     ( x   y)
c
c     where rr is upper triangular.  if r and z have been obtained
c     from the factorization of a least squares problem, then
c     rr and zz are the factors corresponding to the problem
c     with the observation (x,y) removed.  in this case, if rho
c     is the norm of the residual vector, then the norm of
c     the residual vector of the downdated problem is
c     dsqrt(rho**2 - zeta**2). dchdd will simultaneously downdate
c     several triplets (z,y,rho) along with r.
c     for a less terse description of what dchdd does and how
c     it may be applied, see the linpack guide.
c
c     the matrix u is determined as the product u(1)*...*u(p)
c     where u(i) is a rotation in the (p+1,i)-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     the user is warned that a given downdating problem may
c     be impossible to accomplish or may produce
c     inaccurate results.  for example, this can happen
c     if x is near a vector whose removal will reduce the
c     rank of r.  beware.
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 downdated.  the part of  r
c                below the diagonal is not referenced.
c
c         ldr    integer.
c                ldr is the leading dimension fo 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 vector that is to
c                be removed from r.  x is not altered by dchdd.
c
c         z      double precision(ldz,nz), where ldz .ge. p.
c                z is an array of nz p-vectors which
c                are to be downdated along 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 downdated
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 the downdating
c                of the vectors z.  y is not altered by dchdd.
c
c         rho    double precision(nz).
c                rho contains the norms of the residual
c                vectors that are to be downdated.
c
c     on return
c
c         r
c         z      contain the downdated quantities.
c         rho
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         info   integer.
c                info is set as follows.
c
c                   info = 0  if the entire downdating
c                             was successful.
c
c                   info =-1  if r could not be downdated.
c                             in this case, all quantities
c                             are left unaltered.
c
c                   info = 1  if some rho could not be
c                             downdated.  the offending rhos are
c                             set to -1.
c
c     linpack. this version dated 08/14/78 .
c     g.w. stewart, university of maryland, argonne national lab.
c
c     dchdd uses the following functions and subprograms.
c
c     fortran dabs
c     blas ddot, dnrm2
c
      integer i,ii,j
      double precision a,alpha,azeta,norm,dnrm2
      double precision ddot,t,zeta,b,xx
c
c     solve the system trans(r)*a = x, placing the result
c     in the array s.
c
      info = 0
      s(1) = x(1)/r(1,1)
      if (p .lt. 2) go to 20
      do 10 j = 2, p
         s(j) = x(j) - ddot(j-1,r(1,j),1,s,1)
         s(j) = s(j)/r(j,j)
   10 continue
   20 continue
      norm = dnrm2(p,s,1)
      if (norm .lt. 1.0d0) go to 30
         info = -1
      go to 120
   30 continue
         alpha = dsqrt(1.0d0-norm**2)
c
c        determine the transformations.
c
         do 40 ii = 1, p
            i = p - ii + 1
            scale = alpha + dabs(s(i))
            a = alpha/scale
            b = s(i)/scale
            norm = dsqrt(a**2+b**2+0.0d0**2)
            c(i) = a/norm
            s(i) = b/norm
            alpha = scale*norm
   40    continue
c
c        apply the transformations to r.
c
         do 60 j = 1, p
            xx = 0.0d0
            do 50 ii = 1, j
               i = j - ii + 1
               t = c(i)*xx + s(i)*r(i,j)
               r(i,j) = c(i)*r(i,j) - s(i)*xx
               xx = t
   50       continue
   60    continue
c
c        if required, downdate z and rho.
c
         if (nz .lt. 1) go to 110
         do 100 j = 1, nz
            zeta = y(j)
            do 70 i = 1, p
               z(i,j) = (z(i,j) - s(i)*zeta)/c(i)
               zeta = c(i)*zeta - s(i)*z(i,j)
   70       continue
            azeta = dabs(zeta)
            if (azeta .le. rho(j)) go to 80
               info = 1
               rho(j) = -1.0d0
            go to 90
   80       continue
               rho(j) = rho(j)*dsqrt(1.0d0-(azeta/rho(j))**2)
   90       continue
  100    continue
  110    continue
  120 continue
      return
      end