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