| blob | f7e43f257504263ace1a7c5233d33aa41592e7ef |
1 SUBROUTINE DSYTD2( UPLO, N, A, LDA, D, E, TAU, INFO )
2 !
3 ! -- LAPACK routine (version 3.1) --
4 ! Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
5 ! November 2006
6 !
7 ! .. Scalar Arguments ..
8 CHARACTER UPLO
9 INTEGER INFO, LDA, N
10 ! ..
11 ! .. Array Arguments ..
12 DOUBLE PRECISION A( LDA, * ), D( * ), E( * ), TAU( * )
13 ! ..
14 !
15 ! Purpose
16 ! =======
17 !
18 ! DSYTD2 reduces a real symmetric matrix A to symmetric tridiagonal
19 ! form T by an orthogonal similarity transformation: Q' * A * Q = T.
20 !
21 ! Arguments
22 ! =========
23 !
24 ! UPLO (input) CHARACTER*1
25 ! Specifies whether the upper or lower triangular part of the
26 ! symmetric matrix A is stored:
27 ! = 'U': Upper triangular
28 ! = 'L': Lower triangular
29 !
30 ! N (input) INTEGER
31 ! The order of the matrix A. N >= 0.
32 !
33 ! A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
34 ! On entry, the symmetric matrix A. If UPLO = 'U', the leading
35 ! n-by-n upper triangular part of A contains the upper
36 ! triangular part of the matrix A, and the strictly lower
37 ! triangular part of A is not referenced. If UPLO = 'L', the
38 ! leading n-by-n lower triangular part of A contains the lower
39 ! triangular part of the matrix A, and the strictly upper
40 ! triangular part of A is not referenced.
41 ! On exit, if UPLO = 'U', the diagonal and first superdiagonal
42 ! of A are overwritten by the corresponding elements of the
43 ! tridiagonal matrix T, and the elements above the first
44 ! superdiagonal, with the array TAU, represent the orthogonal
45 ! matrix Q as a product of elementary reflectors; if UPLO
46 ! = 'L', the diagonal and first subdiagonal of A are over-
47 ! written by the corresponding elements of the tridiagonal
48 ! matrix T, and the elements below the first subdiagonal, with
49 ! the array TAU, represent the orthogonal matrix Q as a product
50 ! of elementary reflectors. See Further Details.
51 !
52 ! LDA (input) INTEGER
53 ! The leading dimension of the array A. LDA >= max(1,N).
54 !
55 ! D (output) DOUBLE PRECISION array, dimension (N)
56 ! The diagonal elements of the tridiagonal matrix T:
57 ! D(i) = A(i,i).
58 !
59 ! E (output) DOUBLE PRECISION array, dimension (N-1)
60 ! The off-diagonal elements of the tridiagonal matrix T:
61 ! E(i) = A(i,i+1) if UPLO = 'U', E(i) = A(i+1,i) if UPLO = 'L'.
62 !
63 ! TAU (output) DOUBLE PRECISION array, dimension (N-1)
64 ! The scalar factors of the elementary reflectors (see Further
65 ! Details).
66 !
67 ! INFO (output) INTEGER
68 ! = 0: successful exit
69 ! < 0: if INFO = -i, the i-th argument had an illegal value.
70 !
71 ! Further Details
72 ! ===============
73 !
74 ! If UPLO = 'U', the matrix Q is represented as a product of elementary
75 ! reflectors
76 !
77 ! Q = H(n-1) . . . H(2) H(1).
78 !
79 ! Each H(i) has the form
80 !
81 ! H(i) = I - tau * v * v'
82 !
83 ! where tau is a real scalar, and v is a real vector with
84 ! v(i+1:n) = 0 and v(i) = 1; v(1:i-1) is stored on exit in
85 ! A(1:i-1,i+1), and tau in TAU(i).
86 !
87 ! If UPLO = 'L', the matrix Q is represented as a product of elementary
88 ! reflectors
89 !
90 ! Q = H(1) H(2) . . . H(n-1).
91 !
92 ! Each H(i) has the form
93 !
94 ! H(i) = I - tau * v * v'
95 !
96 ! where tau is a real scalar, and v is a real vector with
97 ! v(1:i) = 0 and v(i+1) = 1; v(i+2:n) is stored on exit in A(i+2:n,i),
98 ! and tau in TAU(i).
99 !
100 ! The contents of A on exit are illustrated by the following examples
101 ! with n = 5:
102 !
103 ! if UPLO = 'U': if UPLO = 'L':
104 !
105 ! ( d e v2 v3 v4 ) ( d )
106 ! ( d e v3 v4 ) ( e d )
107 ! ( d e v4 ) ( v1 e d )
108 ! ( d e ) ( v1 v2 e d )
109 ! ( d ) ( v1 v2 v3 e d )
110 !
111 ! where d and e denote diagonal and off-diagonal elements of T, and vi
112 ! denotes an element of the vector defining H(i).
113 !
114 ! =====================================================================
115 !
116 ! .. Parameters ..
117 DOUBLE PRECISION ONE, ZERO, HALF
118 PARAMETER ( ONE = 1.0D0, ZERO = 0.0D0, &
119 HALF = 1.0D0 / 2.0D0 )
120 ! ..
121 ! .. Local Scalars ..
122 LOGICAL UPPER
123 INTEGER I
124 DOUBLE PRECISION ALPHA, TAUI
125 ! ..
126 ! .. External Subroutines ..
127 ! EXTERNAL DAXPY, DLARFG, DSYMV, DSYR2, XERBLA
128 ! ..
129 ! .. External Functions ..
130 ! LOGICAL LSAME
131 ! DOUBLE PRECISION DDOT
132 ! EXTERNAL LSAME, DDOT
133 ! ..
134 ! .. Intrinsic Functions ..
135 INTRINSIC MAX, MIN
136 ! ..
137 ! .. Executable Statements ..
138 !
139 ! Test the input parameters
140 !
141 INFO = 0
142 UPPER = LSAME( UPLO, 'U' )
143 IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
144 INFO = -1
145 ELSE IF( N.LT.0 ) THEN
146 INFO = -2
147 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
148 INFO = -4
149 END IF
150 IF( INFO.NE.0 ) THEN
151 CALL XERBLA( 'DSYTD2', -INFO )
152 RETURN
153 END IF
154 !
155 ! Quick return if possible
156 !
157 IF( N.LE.0 ) &
158 RETURN
159 !
160 IF( UPPER ) THEN
161 !
162 ! Reduce the upper triangle of A
163 !
164 DO 10 I = N - 1, 1, -1
165 !
166 ! Generate elementary reflector H(i) = I - tau * v * v'
167 ! to annihilate A(1:i-1,i+1)
168 !
169 CALL DLARFG( I, A( I, I+1 ), A( 1, I+1 ), 1, TAUI )
170 E( I ) = A( I, I+1 )
171 !
172 IF( TAUI.NE.ZERO ) THEN
173 !
174 ! Apply H(i) from both sides to A(1:i,1:i)
175 !
176 A( I, I+1 ) = ONE
177 !
178 ! Compute x := tau * A * v storing x in TAU(1:i)
179 !
180 CALL DSYMV( UPLO, I, TAUI, A, LDA, A( 1, I+1 ), 1, ZERO, &
181 TAU, 1 )
182 !
183 ! Compute w := x - 1/2 * tau * (x'*v) * v
184 !
185 ALPHA = -HALF*TAUI*DDOT( I, TAU, 1, A( 1, I+1 ), 1 )
186 CALL DAXPY( I, ALPHA, A( 1, I+1 ), 1, TAU, 1 )
187 !
188 ! Apply the transformation as a rank-2 update:
189 ! A := A - v * w' - w * v'
190 !
191 CALL DSYR2( UPLO, I, -ONE, A( 1, I+1 ), 1, TAU, 1, A, &
192 LDA )
193 !
194 A( I, I+1 ) = E( I )
195 END IF
196 D( I+1 ) = A( I+1, I+1 )
197 TAU( I ) = TAUI
198 10 CONTINUE
199 D( 1 ) = A( 1, 1 )
200 ELSE
201 !
202 ! Reduce the lower triangle of A
203 !
204 DO 20 I = 1, N - 1
205 !
206 ! Generate elementary reflector H(i) = I - tau * v * v'
207 ! to annihilate A(i+2:n,i)
208 !
209 CALL DLARFG( N-I, A( I+1, I ), A( MIN( I+2, N ), I ), 1, &
210 TAUI )
211 E( I ) = A( I+1, I )
212 !
213 IF( TAUI.NE.ZERO ) THEN
214 !
215 ! Apply H(i) from both sides to A(i+1:n,i+1:n)
216 !
217 A( I+1, I ) = ONE
218 !
219 ! Compute x := tau * A * v storing y in TAU(i:n-1)
220 !
221 CALL DSYMV( UPLO, N-I, TAUI, A( I+1, I+1 ), LDA, &
222 A( I+1, I ), 1, ZERO, TAU( I ), 1 )
223 !
224 ! Compute w := x - 1/2 * tau * (x'*v) * v
225 !
226 ALPHA = -HALF*TAUI*DDOT( N-I, TAU( I ), 1, A( I+1, I ), &
227 1 )
228 CALL DAXPY( N-I, ALPHA, A( I+1, I ), 1, TAU( I ), 1 )
229 !
230 ! Apply the transformation as a rank-2 update:
231 ! A := A - v * w' - w * v'
232 !
233 CALL DSYR2( UPLO, N-I, -ONE, A( I+1, I ), 1, TAU( I ), 1, &
234 A( I+1, I+1 ), LDA )
235 !
236 A( I+1, I ) = E( I )
237 END IF
238 D( I ) = A( I, I )
239 TAU( I ) = TAUI
240 20 CONTINUE
241 D( N ) = A( N, N )
242 END IF
243 !
244 RETURN
245 !
246 ! End of DSYTD2
247 !
248 END SUBROUTINE DSYTD2