-
Notifications
You must be signed in to change notification settings - Fork 15
/
dchdd.f
202 lines (202 loc) · 6.71 KB
/
dchdd.f
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
*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