Annotation of rpl/lapack/lapack/zhptrd.f, revision 1.16
1.9 bertrand 1: *> \brief \b ZHPTRD
2: *
3: * =========== DOCUMENTATION ===========
4: *
1.15 bertrand 5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
1.9 bertrand 7: *
8: *> \htmlonly
1.15 bertrand 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">
1.9 bertrand 15: *> [TXT]</a>
1.15 bertrand 16: *> \endhtmlonly
1.9 bertrand 17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE ZHPTRD( UPLO, N, AP, D, E, TAU, INFO )
1.15 bertrand 22: *
1.9 bertrand 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: * ..
1.15 bertrand 31: *
1.9 bertrand 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: *
1.15 bertrand 110: *> \author Univ. of Tennessee
111: *> \author Univ. of California Berkeley
112: *> \author Univ. of Colorado Denver
113: *> \author NAG Ltd.
1.9 bertrand 114: *
1.15 bertrand 115: *> \date December 2016
1.9 bertrand 116: *
117: *> \ingroup complex16OTHERcomputational
118: *
119: *> \par Further Details:
120: * =====================
121: *>
122: *> \verbatim
123: *>
124: *> If UPLO = 'U', the matrix Q is represented as a product of elementary
125: *> reflectors
126: *>
127: *> Q = H(n-1) . . . H(2) H(1).
128: *>
129: *> Each H(i) has the form
130: *>
131: *> H(i) = I - tau * v * v**H
132: *>
133: *> where tau is a complex scalar, and v is a complex vector with
134: *> v(i+1:n) = 0 and v(i) = 1; v(1:i-1) is stored on exit in AP,
135: *> overwriting A(1:i-1,i+1), and tau is stored in TAU(i).
136: *>
137: *> If UPLO = 'L', the matrix Q is represented as a product of elementary
138: *> reflectors
139: *>
140: *> Q = H(1) H(2) . . . H(n-1).
141: *>
142: *> Each H(i) has the form
143: *>
144: *> H(i) = I - tau * v * v**H
145: *>
146: *> where tau is a complex scalar, and v is a complex vector with
147: *> v(1:i) = 0 and v(i+1) = 1; v(i+2:n) is stored on exit in AP,
148: *> overwriting A(i+2:n,i), and tau is stored in TAU(i).
149: *> \endverbatim
150: *>
151: * =====================================================================
1.1 bertrand 152: SUBROUTINE ZHPTRD( UPLO, N, AP, D, E, TAU, INFO )
153: *
1.15 bertrand 154: * -- LAPACK computational routine (version 3.7.0) --
1.1 bertrand 155: * -- LAPACK is a software package provided by Univ. of Tennessee, --
156: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.15 bertrand 157: * December 2016
1.1 bertrand 158: *
159: * .. Scalar Arguments ..
160: CHARACTER UPLO
161: INTEGER INFO, N
162: * ..
163: * .. Array Arguments ..
164: DOUBLE PRECISION D( * ), E( * )
165: COMPLEX*16 AP( * ), TAU( * )
166: * ..
167: *
168: * =====================================================================
169: *
170: * .. Parameters ..
171: COMPLEX*16 ONE, ZERO, HALF
172: PARAMETER ( ONE = ( 1.0D+0, 0.0D+0 ),
173: $ ZERO = ( 0.0D+0, 0.0D+0 ),
174: $ HALF = ( 0.5D+0, 0.0D+0 ) )
175: * ..
176: * .. Local Scalars ..
177: LOGICAL UPPER
178: INTEGER I, I1, I1I1, II
179: COMPLEX*16 ALPHA, TAUI
180: * ..
181: * .. External Subroutines ..
182: EXTERNAL XERBLA, ZAXPY, ZHPMV, ZHPR2, ZLARFG
183: * ..
184: * .. External Functions ..
185: LOGICAL LSAME
186: COMPLEX*16 ZDOTC
187: EXTERNAL LSAME, ZDOTC
188: * ..
189: * .. Intrinsic Functions ..
190: INTRINSIC DBLE
191: * ..
192: * .. Executable Statements ..
193: *
194: * Test the input parameters
195: *
196: INFO = 0
197: UPPER = LSAME( UPLO, 'U' )
198: IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
199: INFO = -1
200: ELSE IF( N.LT.0 ) THEN
201: INFO = -2
202: END IF
203: IF( INFO.NE.0 ) THEN
204: CALL XERBLA( 'ZHPTRD', -INFO )
205: RETURN
206: END IF
207: *
208: * Quick return if possible
209: *
210: IF( N.LE.0 )
211: $ RETURN
212: *
213: IF( UPPER ) THEN
214: *
215: * Reduce the upper triangle of A.
216: * I1 is the index in AP of A(1,I+1).
217: *
218: I1 = N*( N-1 ) / 2 + 1
219: AP( I1+N-1 ) = DBLE( AP( I1+N-1 ) )
220: DO 10 I = N - 1, 1, -1
221: *
1.8 bertrand 222: * Generate elementary reflector H(i) = I - tau * v * v**H
1.1 bertrand 223: * to annihilate A(1:i-1,i+1)
224: *
225: ALPHA = AP( I1+I-1 )
226: CALL ZLARFG( I, ALPHA, AP( I1 ), 1, TAUI )
227: E( I ) = ALPHA
228: *
229: IF( TAUI.NE.ZERO ) THEN
230: *
231: * Apply H(i) from both sides to A(1:i,1:i)
232: *
233: AP( I1+I-1 ) = ONE
234: *
235: * Compute y := tau * A * v storing y in TAU(1:i)
236: *
237: CALL ZHPMV( UPLO, I, TAUI, AP, AP( I1 ), 1, ZERO, TAU,
238: $ 1 )
239: *
1.8 bertrand 240: * Compute w := y - 1/2 * tau * (y**H *v) * v
1.1 bertrand 241: *
242: ALPHA = -HALF*TAUI*ZDOTC( I, TAU, 1, AP( I1 ), 1 )
243: CALL ZAXPY( I, ALPHA, AP( I1 ), 1, TAU, 1 )
244: *
245: * Apply the transformation as a rank-2 update:
1.8 bertrand 246: * A := A - v * w**H - w * v**H
1.1 bertrand 247: *
248: CALL ZHPR2( UPLO, I, -ONE, AP( I1 ), 1, TAU, 1, AP )
249: *
250: END IF
251: AP( I1+I-1 ) = E( I )
252: D( I+1 ) = AP( I1+I )
253: TAU( I ) = TAUI
254: I1 = I1 - I
255: 10 CONTINUE
256: D( 1 ) = AP( 1 )
257: ELSE
258: *
259: * Reduce the lower triangle of A. II is the index in AP of
260: * A(i,i) and I1I1 is the index of A(i+1,i+1).
261: *
262: II = 1
263: AP( 1 ) = DBLE( AP( 1 ) )
264: DO 20 I = 1, N - 1
265: I1I1 = II + N - I + 1
266: *
1.8 bertrand 267: * Generate elementary reflector H(i) = I - tau * v * v**H
1.1 bertrand 268: * to annihilate A(i+2:n,i)
269: *
270: ALPHA = AP( II+1 )
271: CALL ZLARFG( N-I, ALPHA, AP( II+2 ), 1, TAUI )
272: E( I ) = ALPHA
273: *
274: IF( TAUI.NE.ZERO ) THEN
275: *
276: * Apply H(i) from both sides to A(i+1:n,i+1:n)
277: *
278: AP( II+1 ) = ONE
279: *
280: * Compute y := tau * A * v storing y in TAU(i:n-1)
281: *
282: CALL ZHPMV( UPLO, N-I, TAUI, AP( I1I1 ), AP( II+1 ), 1,
283: $ ZERO, TAU( I ), 1 )
284: *
1.8 bertrand 285: * Compute w := y - 1/2 * tau * (y**H *v) * v
1.1 bertrand 286: *
287: ALPHA = -HALF*TAUI*ZDOTC( N-I, TAU( I ), 1, AP( II+1 ),
288: $ 1 )
289: CALL ZAXPY( N-I, ALPHA, AP( II+1 ), 1, TAU( I ), 1 )
290: *
291: * Apply the transformation as a rank-2 update:
1.8 bertrand 292: * A := A - v * w**H - w * v**H
1.1 bertrand 293: *
294: CALL ZHPR2( UPLO, N-I, -ONE, AP( II+1 ), 1, TAU( I ), 1,
295: $ AP( I1I1 ) )
296: *
297: END IF
298: AP( II+1 ) = E( I )
299: D( I ) = AP( II )
300: TAU( I ) = TAUI
301: II = I1I1
302: 20 CONTINUE
303: D( N ) = AP( II )
304: END IF
305: *
306: RETURN
307: *
308: * End of ZHPTRD
309: *
310: 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