Annotation of rpl/lapack/lapack/zhptrd.f, revision 1.9
1.9 ! bertrand 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: *> \date November 2011
! 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.9 ! bertrand 154: * -- LAPACK computational routine (version 3.4.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.9 ! bertrand 157: * November 2011
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