1 *> \brief \b CHETD2 reduces a Hermitian matrix to real symmetric tridiagonal form by an unitary similarity transformation (unblocked algorithm).
3 * =========== DOCUMENTATION ===========
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
9 *> Download CHETD2 + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/chetd2.f">
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/chetd2.f">
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/chetd2.f">
21 * SUBROUTINE CHETD2( UPLO, N, A, LDA, D, E, TAU, INFO )
23 * .. Scalar Arguments ..
25 * INTEGER INFO, LDA, N
27 * .. Array Arguments ..
29 * COMPLEX A( LDA, * ), TAU( * )
38 *> CHETD2 reduces a complex Hermitian matrix A to real symmetric
39 *> tridiagonal form T by a unitary similarity transformation:
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
58 *> The order of the matrix A. N >= 0.
63 *> A is COMPLEX 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.
86 *> The leading dimension of the array A. LDA >= max(1,N).
91 *> D is REAL array, dimension (N)
92 *> The diagonal elements of the tridiagonal matrix T:
98 *> E is REAL 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'.
105 *> TAU is COMPLEX array, dimension (N-1)
106 *> The scalar factors of the elementary reflectors (see Further
113 *> = 0: successful exit
114 *> < 0: if INFO = -i, the i-th argument had an illegal value.
120 *> \author Univ. of Tennessee
121 *> \author Univ. of California Berkeley
122 *> \author Univ. of Colorado Denver
125 *> \date September 2012
127 *> \ingroup complexHEcomputational
129 *> \par Further Details:
130 * =====================
134 *> If UPLO = 'U', the matrix Q is represented as a product of elementary
137 *> Q = H(n-1) . . . H(2) H(1).
139 *> Each H(i) has the form
141 *> H(i) = I - tau * v * v**H
143 *> where tau is a complex scalar, and v is a complex vector with
144 *> v(i+1:n) = 0 and v(i) = 1; v(1:i-1) is stored on exit in
145 *> A(1:i-1,i+1), and tau in TAU(i).
147 *> If UPLO = 'L', the matrix Q is represented as a product of elementary
150 *> Q = H(1) H(2) . . . H(n-1).
152 *> Each H(i) has the form
154 *> H(i) = I - tau * v * v**H
156 *> where tau is a complex scalar, and v is a complex vector with
157 *> v(1:i) = 0 and v(i+1) = 1; v(i+2:n) is stored on exit in A(i+2:n,i),
158 *> and tau in TAU(i).
160 *> The contents of A on exit are illustrated by the following examples
163 *> if UPLO = 'U': if UPLO = 'L':
165 *> ( d e v2 v3 v4 ) ( d )
166 *> ( d e v3 v4 ) ( e d )
167 *> ( d e v4 ) ( v1 e d )
168 *> ( d e ) ( v1 v2 e d )
169 *> ( d ) ( v1 v2 v3 e d )
171 *> where d and e denote diagonal and off-diagonal elements of T, and vi
172 *> denotes an element of the vector defining H(i).
175 * =====================================================================
176 SUBROUTINE CHETD2( UPLO, N, A, LDA, D, E, TAU, INFO )
178 * -- LAPACK computational routine (version 3.4.2) --
179 * -- LAPACK is a software package provided by Univ. of Tennessee, --
180 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
183 * .. Scalar Arguments ..
187 * .. Array Arguments ..
189 COMPLEX A( LDA, * ), TAU( * )
192 * =====================================================================
195 COMPLEX ONE, ZERO, HALF
196 PARAMETER ( ONE = ( 1.0E+0, 0.0E+0 ),
197 $ ZERO = ( 0.0E+0, 0.0E+0 ),
198 $ HALF = ( 0.5E+0, 0.0E+0 ) )
200 * .. Local Scalars ..
205 * .. External Subroutines ..
206 EXTERNAL CAXPY, CHEMV, CHER2, CLARFG, XERBLA
208 * .. External Functions ..
211 EXTERNAL LSAME, CDOTC
213 * .. Intrinsic Functions ..
214 INTRINSIC MAX, MIN, REAL
216 * .. Executable Statements ..
218 * Test the input parameters
221 UPPER = LSAME( UPLO, 'U' )
222 IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
224 ELSE IF( N.LT.0 ) THEN
226 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
230 CALL XERBLA( 'CHETD2', -INFO )
234 * Quick return if possible
241 * Reduce the upper triangle of A
243 A( N, N ) = REAL( A( N, N ) )
244 DO 10 I = N - 1, 1, -1
246 * Generate elementary reflector H(i) = I - tau * v * v**H
247 * to annihilate A(1:i-1,i+1)
250 CALL CLARFG( I, ALPHA, A( 1, I+1 ), 1, TAUI )
253 IF( TAUI.NE.ZERO ) THEN
255 * Apply H(i) from both sides to A(1:i,1:i)
259 * Compute x := tau * A * v storing x in TAU(1:i)
261 CALL CHEMV( UPLO, I, TAUI, A, LDA, A( 1, I+1 ), 1, ZERO,
264 * Compute w := x - 1/2 * tau * (x**H * v) * v
266 ALPHA = -HALF*TAUI*CDOTC( I, TAU, 1, A( 1, I+1 ), 1 )
267 CALL CAXPY( I, ALPHA, A( 1, I+1 ), 1, TAU, 1 )
269 * Apply the transformation as a rank-2 update:
270 * A := A - v * w**H - w * v**H
272 CALL CHER2( UPLO, I, -ONE, A( 1, I+1 ), 1, TAU, 1, A,
276 A( I, I ) = REAL( A( I, I ) )
279 D( I+1 ) = A( I+1, I+1 )
285 * Reduce the lower triangle of A
287 A( 1, 1 ) = REAL( A( 1, 1 ) )
290 * Generate elementary reflector H(i) = I - tau * v * v**H
291 * to annihilate A(i+2:n,i)
294 CALL CLARFG( N-I, ALPHA, A( MIN( I+2, N ), I ), 1, TAUI )
297 IF( TAUI.NE.ZERO ) THEN
299 * Apply H(i) from both sides to A(i+1:n,i+1:n)
303 * Compute x := tau * A * v storing y in TAU(i:n-1)
305 CALL CHEMV( UPLO, N-I, TAUI, A( I+1, I+1 ), LDA,
306 $ A( I+1, I ), 1, ZERO, TAU( I ), 1 )
308 * Compute w := x - 1/2 * tau * (x**H * v) * v
310 ALPHA = -HALF*TAUI*CDOTC( N-I, TAU( I ), 1, A( I+1, I ),
312 CALL CAXPY( N-I, ALPHA, A( I+1, I ), 1, TAU( I ), 1 )
314 * Apply the transformation as a rank-2 update:
315 * A := A - v * w**H - w * v**H
317 CALL CHER2( UPLO, N-I, -ONE, A( I+1, I ), 1, TAU( I ), 1,
318 $ A( I+1, I+1 ), LDA )
321 A( I+1, I+1 ) = REAL( A( I+1, I+1 ) )