LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
dsytd2.f
Go to the documentation of this file.
1*> \brief \b DSYTD2 reduces a symmetric matrix to real symmetric tridiagonal form by an orthogonal similarity transformation (unblocked algorithm).
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download DSYTD2 + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dsytd2.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dsytd2.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dsytd2.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE DSYTD2( UPLO, N, A, LDA, D, E, TAU, INFO )
22*
23* .. Scalar Arguments ..
24* CHARACTER UPLO
25* INTEGER INFO, LDA, N
26* ..
27* .. Array Arguments ..
28* DOUBLE PRECISION A( LDA, * ), D( * ), E( * ), TAU( * )
29* ..
30*
31*
32*> \par Purpose:
33* =============
34*>
35*> \verbatim
36*>
37*> DSYTD2 reduces a real symmetric matrix A to symmetric tridiagonal
38*> form T by an orthogonal similarity transformation: Q**T * A * Q = T.
39*> \endverbatim
40*
41* Arguments:
42* ==========
43*
44*> \param[in] UPLO
45*> \verbatim
46*> UPLO is CHARACTER*1
47*> Specifies whether the upper or lower triangular part of the
48*> symmetric matrix A is stored:
49*> = 'U': Upper triangular
50*> = 'L': Lower triangular
51*> \endverbatim
52*>
53*> \param[in] N
54*> \verbatim
55*> N is INTEGER
56*> The order of the matrix A. N >= 0.
57*> \endverbatim
58*>
59*> \param[in,out] A
60*> \verbatim
61*> A is DOUBLE PRECISION array, dimension (LDA,N)
62*> On entry, the symmetric matrix A. If UPLO = 'U', the leading
63*> n-by-n upper triangular part of A contains the upper
64*> triangular part of the matrix A, and the strictly lower
65*> triangular part of A is not referenced. If UPLO = 'L', the
66*> leading n-by-n lower triangular part of A contains the lower
67*> triangular part of the matrix A, and the strictly upper
68*> triangular part of A is not referenced.
69*> On exit, if UPLO = 'U', the diagonal and first superdiagonal
70*> of A are overwritten by the corresponding elements of the
71*> tridiagonal matrix T, and the elements above the first
72*> superdiagonal, with the array TAU, represent the orthogonal
73*> matrix Q as a product of elementary reflectors; if UPLO
74*> = 'L', the diagonal and first subdiagonal of A are over-
75*> written by the corresponding elements of the tridiagonal
76*> matrix T, and the elements below the first subdiagonal, with
77*> the array TAU, represent the orthogonal matrix Q as a product
78*> of elementary reflectors. See Further Details.
79*> \endverbatim
80*>
81*> \param[in] LDA
82*> \verbatim
83*> LDA is INTEGER
84*> The leading dimension of the array A. LDA >= max(1,N).
85*> \endverbatim
86*>
87*> \param[out] D
88*> \verbatim
89*> D is DOUBLE PRECISION array, dimension (N)
90*> The diagonal elements of the tridiagonal matrix T:
91*> D(i) = A(i,i).
92*> \endverbatim
93*>
94*> \param[out] E
95*> \verbatim
96*> E is DOUBLE PRECISION array, dimension (N-1)
97*> The off-diagonal elements of the tridiagonal matrix T:
98*> E(i) = A(i,i+1) if UPLO = 'U', E(i) = A(i+1,i) if UPLO = 'L'.
99*> \endverbatim
100*>
101*> \param[out] TAU
102*> \verbatim
103*> TAU is DOUBLE PRECISION array, dimension (N-1)
104*> The scalar factors of the elementary reflectors (see Further
105*> Details).
106*> \endverbatim
107*>
108*> \param[out] INFO
109*> \verbatim
110*> INFO is INTEGER
111*> = 0: successful exit
112*> < 0: if INFO = -i, the i-th argument had an illegal value.
113*> \endverbatim
114*
115* Authors:
116* ========
117*
118*> \author Univ. of Tennessee
119*> \author Univ. of California Berkeley
120*> \author Univ. of Colorado Denver
121*> \author NAG Ltd.
122*
123*> \ingroup hetd2
124*
125*> \par Further Details:
126* =====================
127*>
128*> \verbatim
129*>
130*> If UPLO = 'U', the matrix Q is represented as a product of elementary
131*> reflectors
132*>
133*> Q = H(n-1) . . . H(2) H(1).
134*>
135*> Each H(i) has the form
136*>
137*> H(i) = I - tau * v * v**T
138*>
139*> where tau is a real scalar, and v is a real vector with
140*> v(i+1:n) = 0 and v(i) = 1; v(1:i-1) is stored on exit in
141*> A(1:i-1,i+1), and tau in TAU(i).
142*>
143*> If UPLO = 'L', the matrix Q is represented as a product of elementary
144*> reflectors
145*>
146*> Q = H(1) H(2) . . . H(n-1).
147*>
148*> Each H(i) has the form
149*>
150*> H(i) = I - tau * v * v**T
151*>
152*> where tau is a real scalar, and v is a real vector with
153*> v(1:i) = 0 and v(i+1) = 1; v(i+2:n) is stored on exit in A(i+2:n,i),
154*> and tau in TAU(i).
155*>
156*> The contents of A on exit are illustrated by the following examples
157*> with n = 5:
158*>
159*> if UPLO = 'U': if UPLO = 'L':
160*>
161*> ( d e v2 v3 v4 ) ( d )
162*> ( d e v3 v4 ) ( e d )
163*> ( d e v4 ) ( v1 e d )
164*> ( d e ) ( v1 v2 e d )
165*> ( d ) ( v1 v2 v3 e d )
166*>
167*> where d and e denote diagonal and off-diagonal elements of T, and vi
168*> denotes an element of the vector defining H(i).
169*> \endverbatim
170*>
171* =====================================================================
172 SUBROUTINE dsytd2( UPLO, N, A, LDA, D, E, TAU, INFO )
173*
174* -- LAPACK computational routine --
175* -- LAPACK is a software package provided by Univ. of Tennessee, --
176* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
177*
178* .. Scalar Arguments ..
179 CHARACTER UPLO
180 INTEGER INFO, LDA, N
181* ..
182* .. Array Arguments ..
183 DOUBLE PRECISION A( LDA, * ), D( * ), E( * ), TAU( * )
184* ..
185*
186* =====================================================================
187*
188* .. Parameters ..
189 DOUBLE PRECISION ONE, ZERO, HALF
190 parameter( one = 1.0d0, zero = 0.0d0,
191 $ half = 1.0d0 / 2.0d0 )
192* ..
193* .. Local Scalars ..
194 LOGICAL UPPER
195 INTEGER I
196 DOUBLE PRECISION ALPHA, TAUI
197* ..
198* .. External Subroutines ..
199 EXTERNAL daxpy, dlarfg, dsymv, dsyr2, xerbla
200* ..
201* .. External Functions ..
202 LOGICAL LSAME
203 DOUBLE PRECISION DDOT
204 EXTERNAL lsame, ddot
205* ..
206* .. Intrinsic Functions ..
207 INTRINSIC max, min
208* ..
209* .. Executable Statements ..
210*
211* Test the input parameters
212*
213 info = 0
214 upper = lsame( uplo, 'U' )
215 IF( .NOT.upper .AND. .NOT.lsame( uplo, 'L' ) ) THEN
216 info = -1
217 ELSE IF( n.LT.0 ) THEN
218 info = -2
219 ELSE IF( lda.LT.max( 1, n ) ) THEN
220 info = -4
221 END IF
222 IF( info.NE.0 ) THEN
223 CALL xerbla( 'DSYTD2', -info )
224 RETURN
225 END IF
226*
227* Quick return if possible
228*
229 IF( n.LE.0 )
230 $ RETURN
231*
232 IF( upper ) THEN
233*
234* Reduce the upper triangle of A
235*
236 DO 10 i = n - 1, 1, -1
237*
238* Generate elementary reflector H(i) = I - tau * v * v**T
239* to annihilate A(1:i-1,i+1)
240*
241 CALL dlarfg( i, a( i, i+1 ), a( 1, i+1 ), 1, taui )
242 e( i ) = a( i, i+1 )
243*
244 IF( taui.NE.zero ) THEN
245*
246* Apply H(i) from both sides to A(1:i,1:i)
247*
248 a( i, i+1 ) = one
249*
250* Compute x := tau * A * v storing x in TAU(1:i)
251*
252 CALL dsymv( uplo, i, taui, a, lda, a( 1, i+1 ), 1, zero,
253 $ tau, 1 )
254*
255* Compute w := x - 1/2 * tau * (x**T * v) * v
256*
257 alpha = -half*taui*ddot( i, tau, 1, a( 1, i+1 ), 1 )
258 CALL daxpy( i, alpha, a( 1, i+1 ), 1, tau, 1 )
259*
260* Apply the transformation as a rank-2 update:
261* A := A - v * w**T - w * v**T
262*
263 CALL dsyr2( uplo, i, -one, a( 1, i+1 ), 1, tau, 1, a,
264 $ lda )
265*
266 a( i, i+1 ) = e( i )
267 END IF
268 d( i+1 ) = a( i+1, i+1 )
269 tau( i ) = taui
270 10 CONTINUE
271 d( 1 ) = a( 1, 1 )
272 ELSE
273*
274* Reduce the lower triangle of A
275*
276 DO 20 i = 1, n - 1
277*
278* Generate elementary reflector H(i) = I - tau * v * v**T
279* to annihilate A(i+2:n,i)
280*
281 CALL dlarfg( n-i, a( i+1, i ), a( min( i+2, n ), i ), 1,
282 $ taui )
283 e( i ) = a( i+1, i )
284*
285 IF( taui.NE.zero ) THEN
286*
287* Apply H(i) from both sides to A(i+1:n,i+1:n)
288*
289 a( i+1, i ) = one
290*
291* Compute x := tau * A * v storing y in TAU(i:n-1)
292*
293 CALL dsymv( uplo, n-i, taui, a( i+1, i+1 ), lda,
294 $ a( i+1, i ), 1, zero, tau( i ), 1 )
295*
296* Compute w := x - 1/2 * tau * (x**T * v) * v
297*
298 alpha = -half*taui*ddot( n-i, tau( i ), 1, a( i+1, i ),
299 $ 1 )
300 CALL daxpy( n-i, alpha, a( i+1, i ), 1, tau( i ), 1 )
301*
302* Apply the transformation as a rank-2 update:
303* A := A - v * w**T - w * v**T
304*
305 CALL dsyr2( uplo, n-i, -one, a( i+1, i ), 1, tau( i ), 1,
306 $ a( i+1, i+1 ), lda )
307*
308 a( i+1, i ) = e( i )
309 END IF
310 d( i ) = a( i, i )
311 tau( i ) = taui
312 20 CONTINUE
313 d( n ) = a( n, n )
314 END IF
315*
316 RETURN
317*
318* End of DSYTD2
319*
320 END
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine daxpy(n, da, dx, incx, dy, incy)
DAXPY
Definition daxpy.f:89
subroutine dsymv(uplo, n, alpha, a, lda, x, incx, beta, y, incy)
DSYMV
Definition dsymv.f:152
subroutine dsyr2(uplo, n, alpha, x, incx, y, incy, a, lda)
DSYR2
Definition dsyr2.f:147
subroutine dsytd2(uplo, n, a, lda, d, e, tau, info)
DSYTD2 reduces a symmetric matrix to real symmetric tridiagonal form by an orthogonal similarity tran...
Definition dsytd2.f:173
subroutine dlarfg(n, alpha, x, incx, tau)
DLARFG generates an elementary reflector (Householder matrix).
Definition dlarfg.f:106