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