1: *> \brief \b ZHPTRD
2: *
3: * =========== DOCUMENTATION ===========
4: *
5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
7: *
8: *> \htmlonly
9: *> Download ZHPTRD + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zhptrd.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zhptrd.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zhptrd.f">
15: *> [TXT]</a>
16: *> \endhtmlonly
17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE ZHPTRD( UPLO, N, AP, D, E, TAU, INFO )
22: *
23: * .. Scalar Arguments ..
24: * CHARACTER UPLO
25: * INTEGER INFO, N
26: * ..
27: * .. Array Arguments ..
28: * DOUBLE PRECISION D( * ), E( * )
29: * COMPLEX*16 AP( * ), TAU( * )
30: * ..
31: *
32: *
33: *> \par Purpose:
34: * =============
35: *>
36: *> \verbatim
37: *>
38: *> ZHPTRD reduces a complex Hermitian matrix A stored in packed form to
39: *> real symmetric tridiagonal form T by a unitary similarity
40: *> transformation: Q**H * A * Q = T.
41: *> \endverbatim
42: *
43: * Arguments:
44: * ==========
45: *
46: *> \param[in] UPLO
47: *> \verbatim
48: *> UPLO is CHARACTER*1
49: *> = 'U': Upper triangle of A is stored;
50: *> = 'L': Lower triangle of A is stored.
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] AP
60: *> \verbatim
61: *> AP is COMPLEX*16 array, dimension (N*(N+1)/2)
62: *> On entry, the upper or lower triangle of the Hermitian matrix
63: *> A, packed columnwise in a linear array. The j-th column of A
64: *> is stored in the array AP as follows:
65: *> if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
66: *> if UPLO = 'L', AP(i + (j-1)*(2*n-j)/2) = A(i,j) for j<=i<=n.
67: *> On exit, if UPLO = 'U', the diagonal and first superdiagonal
68: *> of A are overwritten by the corresponding elements of the
69: *> tridiagonal matrix T, and the elements above the first
70: *> superdiagonal, with the array TAU, represent the unitary
71: *> matrix Q as a product of elementary reflectors; if UPLO
72: *> = 'L', the diagonal and first subdiagonal of A are over-
73: *> written by the corresponding elements of the tridiagonal
74: *> matrix T, and the elements below the first subdiagonal, with
75: *> the array TAU, represent the unitary matrix Q as a product
76: *> of elementary reflectors. See Further Details.
77: *> \endverbatim
78: *>
79: *> \param[out] D
80: *> \verbatim
81: *> D is DOUBLE PRECISION array, dimension (N)
82: *> The diagonal elements of the tridiagonal matrix T:
83: *> D(i) = A(i,i).
84: *> \endverbatim
85: *>
86: *> \param[out] E
87: *> \verbatim
88: *> E is DOUBLE PRECISION array, dimension (N-1)
89: *> The off-diagonal elements of the tridiagonal matrix T:
90: *> E(i) = A(i,i+1) if UPLO = 'U', E(i) = A(i+1,i) if UPLO = 'L'.
91: *> \endverbatim
92: *>
93: *> \param[out] TAU
94: *> \verbatim
95: *> TAU is COMPLEX*16 array, dimension (N-1)
96: *> The scalar factors of the elementary reflectors (see Further
97: *> Details).
98: *> \endverbatim
99: *>
100: *> \param[out] INFO
101: *> \verbatim
102: *> INFO is INTEGER
103: *> = 0: successful exit
104: *> < 0: if INFO = -i, the i-th argument had an illegal value
105: *> \endverbatim
106: *
107: * Authors:
108: * ========
109: *
110: *> \author Univ. of Tennessee
111: *> \author Univ. of California Berkeley
112: *> \author Univ. of Colorado Denver
113: *> \author NAG Ltd.
114: *
115: *> \ingroup complex16OTHERcomputational
116: *
117: *> \par Further Details:
118: * =====================
119: *>
120: *> \verbatim
121: *>
122: *> If UPLO = 'U', the matrix Q is represented as a product of elementary
123: *> reflectors
124: *>
125: *> Q = H(n-1) . . . H(2) H(1).
126: *>
127: *> Each H(i) has the form
128: *>
129: *> H(i) = I - tau * v * v**H
130: *>
131: *> where tau is a complex scalar, and v is a complex vector with
132: *> v(i+1:n) = 0 and v(i) = 1; v(1:i-1) is stored on exit in AP,
133: *> overwriting A(1:i-1,i+1), and tau is stored in TAU(i).
134: *>
135: *> If UPLO = 'L', the matrix Q is represented as a product of elementary
136: *> reflectors
137: *>
138: *> Q = H(1) H(2) . . . H(n-1).
139: *>
140: *> Each H(i) has the form
141: *>
142: *> H(i) = I - tau * v * v**H
143: *>
144: *> where tau is a complex scalar, and v is a complex vector with
145: *> v(1:i) = 0 and v(i+1) = 1; v(i+2:n) is stored on exit in AP,
146: *> overwriting A(i+2:n,i), and tau is stored in TAU(i).
147: *> \endverbatim
148: *>
149: * =====================================================================
150: SUBROUTINE ZHPTRD( UPLO, N, AP, D, E, TAU, INFO )
151: *
152: * -- LAPACK computational routine --
153: * -- LAPACK is a software package provided by Univ. of Tennessee, --
154: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
155: *
156: * .. Scalar Arguments ..
157: CHARACTER UPLO
158: INTEGER INFO, N
159: * ..
160: * .. Array Arguments ..
161: DOUBLE PRECISION D( * ), E( * )
162: COMPLEX*16 AP( * ), TAU( * )
163: * ..
164: *
165: * =====================================================================
166: *
167: * .. Parameters ..
168: COMPLEX*16 ONE, ZERO, HALF
169: PARAMETER ( ONE = ( 1.0D+0, 0.0D+0 ),
170: $ ZERO = ( 0.0D+0, 0.0D+0 ),
171: $ HALF = ( 0.5D+0, 0.0D+0 ) )
172: * ..
173: * .. Local Scalars ..
174: LOGICAL UPPER
175: INTEGER I, I1, I1I1, II
176: COMPLEX*16 ALPHA, TAUI
177: * ..
178: * .. External Subroutines ..
179: EXTERNAL XERBLA, ZAXPY, ZHPMV, ZHPR2, ZLARFG
180: * ..
181: * .. External Functions ..
182: LOGICAL LSAME
183: COMPLEX*16 ZDOTC
184: EXTERNAL LSAME, ZDOTC
185: * ..
186: * .. Intrinsic Functions ..
187: INTRINSIC DBLE
188: * ..
189: * .. Executable Statements ..
190: *
191: * Test the input parameters
192: *
193: INFO = 0
194: UPPER = LSAME( UPLO, 'U' )
195: IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
196: INFO = -1
197: ELSE IF( N.LT.0 ) THEN
198: INFO = -2
199: END IF
200: IF( INFO.NE.0 ) THEN
201: CALL XERBLA( 'ZHPTRD', -INFO )
202: RETURN
203: END IF
204: *
205: * Quick return if possible
206: *
207: IF( N.LE.0 )
208: $ RETURN
209: *
210: IF( UPPER ) THEN
211: *
212: * Reduce the upper triangle of A.
213: * I1 is the index in AP of A(1,I+1).
214: *
215: I1 = N*( N-1 ) / 2 + 1
216: AP( I1+N-1 ) = DBLE( AP( I1+N-1 ) )
217: DO 10 I = N - 1, 1, -1
218: *
219: * Generate elementary reflector H(i) = I - tau * v * v**H
220: * to annihilate A(1:i-1,i+1)
221: *
222: ALPHA = AP( I1+I-1 )
223: CALL ZLARFG( I, ALPHA, AP( I1 ), 1, TAUI )
224: E( I ) = DBLE( ALPHA )
225: *
226: IF( TAUI.NE.ZERO ) THEN
227: *
228: * Apply H(i) from both sides to A(1:i,1:i)
229: *
230: AP( I1+I-1 ) = ONE
231: *
232: * Compute y := tau * A * v storing y in TAU(1:i)
233: *
234: CALL ZHPMV( UPLO, I, TAUI, AP, AP( I1 ), 1, ZERO, TAU,
235: $ 1 )
236: *
237: * Compute w := y - 1/2 * tau * (y**H *v) * v
238: *
239: ALPHA = -HALF*TAUI*ZDOTC( I, TAU, 1, AP( I1 ), 1 )
240: CALL ZAXPY( I, ALPHA, AP( I1 ), 1, TAU, 1 )
241: *
242: * Apply the transformation as a rank-2 update:
243: * A := A - v * w**H - w * v**H
244: *
245: CALL ZHPR2( UPLO, I, -ONE, AP( I1 ), 1, TAU, 1, AP )
246: *
247: END IF
248: AP( I1+I-1 ) = E( I )
249: D( I+1 ) = DBLE( AP( I1+I ) )
250: TAU( I ) = TAUI
251: I1 = I1 - I
252: 10 CONTINUE
253: D( 1 ) = DBLE( AP( 1 ) )
254: ELSE
255: *
256: * Reduce the lower triangle of A. II is the index in AP of
257: * A(i,i) and I1I1 is the index of A(i+1,i+1).
258: *
259: II = 1
260: AP( 1 ) = DBLE( AP( 1 ) )
261: DO 20 I = 1, N - 1
262: I1I1 = II + N - I + 1
263: *
264: * Generate elementary reflector H(i) = I - tau * v * v**H
265: * to annihilate A(i+2:n,i)
266: *
267: ALPHA = AP( II+1 )
268: CALL ZLARFG( N-I, ALPHA, AP( II+2 ), 1, TAUI )
269: E( I ) = DBLE( ALPHA )
270: *
271: IF( TAUI.NE.ZERO ) THEN
272: *
273: * Apply H(i) from both sides to A(i+1:n,i+1:n)
274: *
275: AP( II+1 ) = ONE
276: *
277: * Compute y := tau * A * v storing y in TAU(i:n-1)
278: *
279: CALL ZHPMV( UPLO, N-I, TAUI, AP( I1I1 ), AP( II+1 ), 1,
280: $ ZERO, TAU( I ), 1 )
281: *
282: * Compute w := y - 1/2 * tau * (y**H *v) * v
283: *
284: ALPHA = -HALF*TAUI*ZDOTC( N-I, TAU( I ), 1, AP( II+1 ),
285: $ 1 )
286: CALL ZAXPY( N-I, ALPHA, AP( II+1 ), 1, TAU( I ), 1 )
287: *
288: * Apply the transformation as a rank-2 update:
289: * A := A - v * w**H - w * v**H
290: *
291: CALL ZHPR2( UPLO, N-I, -ONE, AP( II+1 ), 1, TAU( I ), 1,
292: $ AP( I1I1 ) )
293: *
294: END IF
295: AP( II+1 ) = E( I )
296: D( I ) = DBLE( AP( II ) )
297: TAU( I ) = TAUI
298: II = I1I1
299: 20 CONTINUE
300: D( N ) = DBLE( AP( II ) )
301: END IF
302: *
303: RETURN
304: *
305: * End of ZHPTRD
306: *
307: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zhptrd.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack