*DECK DCHDD SUBROUTINE DCHDD (R, LDR, P, X, Z, LDZ, NZ, Y, RHO, C, S, INFO) C***BEGIN PROLOGUE DCHDD C***PURPOSE Downdate an augmented Cholesky decomposition or the C triangular factor of an augmented QR decomposition. C***LIBRARY SLATEC (LINPACK) C***CATEGORY D7B C***TYPE DOUBLE PRECISION (SCHDD-S, DCHDD-D, CCHDD-C) C***KEYWORDS CHOLESKY DECOMPOSITION, DOWNDATE, LINEAR ALGEBRA, LINPACK, C MATRIX C***AUTHOR Stewart, G. W., (U. of Maryland) C***DESCRIPTION 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 determines an 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 SQRT(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 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 vector that is to C be removed from R. X is not altered by DCHDD. C C Z DOUBLE PRECISION(LDZ,N)Z), 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 RHO's are C set to -1. C C***REFERENCES J. J. Dongarra, J. R. Bunch, C. B. Moler, and G. W. C Stewart, LINPACK Users' Guide, SIAM, 1979. C***ROUTINES CALLED DDOT, DNRM2 C***REVISION HISTORY (YYMMDD) C 780814 DATE WRITTEN C 890531 Changed all specific intrinsics to generic. (WRB) C 890831 Modified array declarations. (WRB) C 890831 REVISION DATE from Version 3.2 C 891214 Prologue converted to Version 4.0 format. (BAB) C 900326 Removed duplicate information from DESCRIPTION section. C (WRB) C 920501 Reformatted the REFERENCES section. (WRB) C***END PROLOGUE DCHDD INTEGER LDR,P,LDZ,NZ,INFO DOUBLE PRECISION R(LDR,*),X(*),Z(LDZ,*),Y(*),S(*) DOUBLE PRECISION RHO(*),C(*) C INTEGER I,II,J DOUBLE PRECISION A,ALPHA,AZETA,NORM,DNRM2 DOUBLE PRECISION DDOT,T,ZETA,B,XX,SCALE C C SOLVE THE SYSTEM TRANS(R)*A = X, PLACING THE RESULT C IN THE ARRAY S. C C***FIRST EXECUTABLE STATEMENT DCHDD 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 = SQRT(1.0D0-NORM**2) C C DETERMINE THE TRANSFORMATIONS. C DO 40 II = 1, P I = P - II + 1 SCALE = ALPHA + ABS(S(I)) A = ALPHA/SCALE B = S(I)/SCALE NORM = SQRT(A**2+B**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 = ABS(ZETA) IF (AZETA .LE. RHO(J)) GO TO 80 INFO = 1 RHO(J) = -1.0D0 GO TO 90 80 CONTINUE RHO(J) = RHO(J)*SQRT(1.0D0-(AZETA/RHO(J))**2) 90 CONTINUE 100 CONTINUE 110 CONTINUE 120 CONTINUE RETURN END