subroutine cchdd(r,ldr,p,x,z,ldz,nz,y,rho,c,s,info) integer ldr,p,ldz,nz,info complex r(ldr,1),x(1),z(ldz,1),y(1),s(1) real rho(1),c(1) c c cchdd 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, cchdd c determineds a unitary 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 sqrt(rho**2 - zeta**2). cchdd will simultaneously downdate c several triplets (z,y,rho) along with r. c for a less terse description of what cchdd 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) -conjg(s(i)) ) c ( ) . c ( s(i) c(i) ) c c the rotations are chosen so that c(i) is real. 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 complex(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 complex(p). c x contains the row vector that is to c be removed from r. x is not altered by cchdd. c c z complex(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 complex(nz). c y contains the scalars for the downdating c of the vectors z. y is not altered by cchdd. c c rho real(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 real(p). c c contains the cosines of the transforming c rotations. c c s complex(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 cchdd uses the following functions and subprograms. c c fortran aimag,cabs,conjg c blas cdotc, scnrm2 c integer i,ii,j real a,alpha,azeta,norm,scnrm2 complex cdotc,t,zeta,b,xx c c solve the system ctrans(r)*a = x, placing the result c in the array s. c info = 0 s(1) = conjg(x(1))/conjg(r(1,1)) if (p .lt. 2) go to 20 do 10 j = 2, p s(j) = conjg(x(j)) - cdotc(j-1,r(1,j),1,s,1) s(j) = s(j)/conjg(r(j,j)) 10 continue 20 continue norm = scnrm2(p,s,1) if (norm .lt. 1.0e0) go to 30 info = -1 go to 120 30 continue alpha = sqrt(1.0e0-norm**2) c c determine the transformations. c do 40 ii = 1, p i = p - ii + 1 scale = alpha + cabs(s(i)) a = alpha/scale b = s(i)/scale norm = sqrt(a**2+real(b)**2+aimag(b)**2) c(i) = a/norm s(i) = conjg(b)/norm alpha = scale*norm 40 continue c c apply the transformations to r. c do 60 j = 1, p xx = (0.0e0,0.0e0) 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) - conjg(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) - conjg(s(i))*zeta)/c(i) zeta = c(i)*zeta - s(i)*z(i,j) 70 continue azeta = cabs(zeta) if (azeta .le. rho(j)) go to 80 info = 1 rho(j) = -1.0e0 go to 90 80 continue rho(j) = rho(j)*sqrt(1.0e0-(azeta/rho(j))**2) 90 continue 100 continue 110 continue 120 continue return end