LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
ssytd2.f
Go to the documentation of this file.
1*> \brief \b SSYTD2 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 SSYTD2 + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ssytd2.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ssytd2.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ssytd2.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE SSYTD2( 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* REAL A( LDA, * ), D( * ), E( * ), TAU( * )
29* ..
30*
31*
32*> \par Purpose:
33* =============
34*>
35*> \verbatim
36*>
37*> SSYTD2 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 REAL 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 REAL 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 REAL 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 REAL 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 ssytd2( 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 REAL A( LDA, * ), D( * ), E( * ), TAU( * )
184* ..
185*
186* =====================================================================
187*
188* .. Parameters ..
189 REAL ONE, ZERO, HALF
190 parameter( one = 1.0, zero = 0.0, half = 1.0 / 2.0 )
191* ..
192* .. Local Scalars ..
193 LOGICAL UPPER
194 INTEGER I
195 REAL ALPHA, TAUI
196* ..
197* .. External Subroutines ..
198 EXTERNAL saxpy, slarfg, ssymv, ssyr2, xerbla
199* ..
200* .. External Functions ..
201 LOGICAL LSAME
202 REAL SDOT
203 EXTERNAL lsame, sdot
204* ..
205* .. Intrinsic Functions ..
206 INTRINSIC max, min
207* ..
208* .. Executable Statements ..
209*
210* Test the input parameters
211*
212 info = 0
213 upper = lsame( uplo, 'U' )
214 IF( .NOT.upper .AND. .NOT.lsame( uplo, 'L' ) ) THEN
215 info = -1
216 ELSE IF( n.LT.0 ) THEN
217 info = -2
218 ELSE IF( lda.LT.max( 1, n ) ) THEN
219 info = -4
220 END IF
221 IF( info.NE.0 ) THEN
222 CALL xerbla( 'SSYTD2', -info )
223 RETURN
224 END IF
225*
226* Quick return if possible
227*
228 IF( n.LE.0 )
229 $ RETURN
230*
231 IF( upper ) THEN
232*
233* Reduce the upper triangle of A
234*
235 DO 10 i = n - 1, 1, -1
236*
237* Generate elementary reflector H(i) = I - tau * v * v**T
238* to annihilate A(1:i-1,i+1)
239*
240 CALL slarfg( i, a( i, i+1 ), a( 1, i+1 ), 1, taui )
241 e( i ) = a( i, i+1 )
242*
243 IF( taui.NE.zero ) THEN
244*
245* Apply H(i) from both sides to A(1:i,1:i)
246*
247 a( i, i+1 ) = one
248*
249* Compute x := tau * A * v storing x in TAU(1:i)
250*
251 CALL ssymv( uplo, i, taui, a, lda, a( 1, i+1 ), 1, zero,
252 $ tau, 1 )
253*
254* Compute w := x - 1/2 * tau * (x**T * v) * v
255*
256 alpha = -half*taui*sdot( i, tau, 1, a( 1, i+1 ), 1 )
257 CALL saxpy( i, alpha, a( 1, i+1 ), 1, tau, 1 )
258*
259* Apply the transformation as a rank-2 update:
260* A := A - v * w**T - w * v**T
261*
262 CALL ssyr2( uplo, i, -one, a( 1, i+1 ), 1, tau, 1, a,
263 $ lda )
264*
265 a( i, i+1 ) = e( i )
266 END IF
267 d( i+1 ) = a( i+1, i+1 )
268 tau( i ) = taui
269 10 CONTINUE
270 d( 1 ) = a( 1, 1 )
271 ELSE
272*
273* Reduce the lower triangle of A
274*
275 DO 20 i = 1, n - 1
276*
277* Generate elementary reflector H(i) = I - tau * v * v**T
278* to annihilate A(i+2:n,i)
279*
280 CALL slarfg( n-i, a( i+1, i ), a( min( i+2, n ), i ), 1,
281 $ taui )
282 e( i ) = a( i+1, i )
283*
284 IF( taui.NE.zero ) THEN
285*
286* Apply H(i) from both sides to A(i+1:n,i+1:n)
287*
288 a( i+1, i ) = one
289*
290* Compute x := tau * A * v storing y in TAU(i:n-1)
291*
292 CALL ssymv( uplo, n-i, taui, a( i+1, i+1 ), lda,
293 $ a( i+1, i ), 1, zero, tau( i ), 1 )
294*
295* Compute w := x - 1/2 * tau * (x**T * v) * v
296*
297 alpha = -half*taui*sdot( n-i, tau( i ), 1, a( i+1, i ),
298 $ 1 )
299 CALL saxpy( n-i, alpha, a( i+1, i ), 1, tau( i ), 1 )
300*
301* Apply the transformation as a rank-2 update:
302* A := A - v * w**T - w * v**T
303*
304 CALL ssyr2( uplo, n-i, -one, a( i+1, i ), 1, tau( i ), 1,
305 $ a( i+1, i+1 ), lda )
306*
307 a( i+1, i ) = e( i )
308 END IF
309 d( i ) = a( i, i )
310 tau( i ) = taui
311 20 CONTINUE
312 d( n ) = a( n, n )
313 END IF
314*
315 RETURN
316*
317* End of SSYTD2
318*
319 END
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine saxpy(n, sa, sx, incx, sy, incy)
SAXPY
Definition saxpy.f:89
subroutine ssymv(uplo, n, alpha, a, lda, x, incx, beta, y, incy)
SSYMV
Definition ssymv.f:152
subroutine ssyr2(uplo, n, alpha, x, incx, y, incy, a, lda)
SSYR2
Definition ssyr2.f:147
subroutine ssytd2(uplo, n, a, lda, d, e, tau, info)
SSYTD2 reduces a symmetric matrix to real symmetric tridiagonal form by an orthogonal similarity tran...
Definition ssytd2.f:173
subroutine slarfg(n, alpha, x, incx, tau)
SLARFG generates an elementary reflector (Householder matrix).
Definition slarfg.f:106