1: SUBROUTINE ZHETD2( UPLO, N, A, LDA, D, E, TAU, INFO )
2: *
3: * -- LAPACK routine (version 3.2) --
4: * -- LAPACK is a software package provided by Univ. of Tennessee, --
5: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
6: * November 2006
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' * 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'
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'
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'
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'*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' - w * v'
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'
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'*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' - w * v'
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
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zhetd2.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack