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