3 * =========== DOCUMENTATION ===========
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
9 *> Download CHPTRD + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/chptrd.f">
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/chptrd.f">
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/chptrd.f">
21 * SUBROUTINE CHPTRD( UPLO, N, AP, D, E, TAU, INFO )
23 * .. Scalar Arguments ..
27 * .. Array Arguments ..
29 * COMPLEX AP( * ), TAU( * )
38 *> CHPTRD 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.
48 *> UPLO is CHARACTER*1
49 *> = 'U': Upper triangle of A is stored;
50 *> = 'L': Lower triangle of A is stored.
56 *> The order of the matrix A. N >= 0.
61 *> AP is COMPLEX 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.
81 *> D is REAL array, dimension (N)
82 *> The diagonal elements of the tridiagonal matrix T:
88 *> E is REAL 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'.
95 *> TAU is COMPLEX array, dimension (N-1)
96 *> The scalar factors of the elementary reflectors (see Further
103 *> = 0: successful exit
104 *> < 0: if INFO = -i, the i-th argument had an illegal value
110 *> \author Univ. of Tennessee
111 *> \author Univ. of California Berkeley
112 *> \author Univ. of Colorado Denver
115 *> \date November 2011
117 *> \ingroup complexOTHERcomputational
119 *> \par Further Details:
120 * =====================
124 *> If UPLO = 'U', the matrix Q is represented as a product of elementary
127 *> Q = H(n-1) . . . H(2) H(1).
129 *> Each H(i) has the form
131 *> H(i) = I - tau * v * v**H
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).
137 *> If UPLO = 'L', the matrix Q is represented as a product of elementary
140 *> Q = H(1) H(2) . . . H(n-1).
142 *> Each H(i) has the form
144 *> H(i) = I - tau * v * v**H
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).
151 * =====================================================================
152 SUBROUTINE CHPTRD( UPLO, N, AP, D, E, TAU, INFO )
154 * -- LAPACK computational routine (version 3.4.0) --
155 * -- LAPACK is a software package provided by Univ. of Tennessee, --
156 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
159 * .. Scalar Arguments ..
163 * .. Array Arguments ..
165 COMPLEX AP( * ), TAU( * )
168 * =====================================================================
171 COMPLEX ONE, ZERO, HALF
172 PARAMETER ( ONE = ( 1.0E+0, 0.0E+0 ),
173 $ ZERO = ( 0.0E+0, 0.0E+0 ),
174 $ HALF = ( 0.5E+0, 0.0E+0 ) )
176 * .. Local Scalars ..
178 INTEGER I, I1, I1I1, II
181 * .. External Subroutines ..
182 EXTERNAL CAXPY, CHPMV, CHPR2, CLARFG, XERBLA
184 * .. External Functions ..
187 EXTERNAL LSAME, CDOTC
189 * .. Intrinsic Functions ..
192 * .. Executable Statements ..
194 * Test the input parameters
197 UPPER = LSAME( UPLO, 'U' )
198 IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
200 ELSE IF( N.LT.0 ) THEN
204 CALL XERBLA( 'CHPTRD', -INFO )
208 * Quick return if possible
215 * Reduce the upper triangle of A.
216 * I1 is the index in AP of A(1,I+1).
218 I1 = N*( N-1 ) / 2 + 1
219 AP( I1+N-1 ) = REAL( AP( I1+N-1 ) )
220 DO 10 I = N - 1, 1, -1
222 * Generate elementary reflector H(i) = I - tau * v * v**H
223 * to annihilate A(1:i-1,i+1)
226 CALL CLARFG( I, ALPHA, AP( I1 ), 1, TAUI )
229 IF( TAUI.NE.ZERO ) THEN
231 * Apply H(i) from both sides to A(1:i,1:i)
235 * Compute y := tau * A * v storing y in TAU(1:i)
237 CALL CHPMV( UPLO, I, TAUI, AP, AP( I1 ), 1, ZERO, TAU,
240 * Compute w := y - 1/2 * tau * (y**H *v) * v
242 ALPHA = -HALF*TAUI*CDOTC( I, TAU, 1, AP( I1 ), 1 )
243 CALL CAXPY( I, ALPHA, AP( I1 ), 1, TAU, 1 )
245 * Apply the transformation as a rank-2 update:
246 * A := A - v * w**H - w * v**H
248 CALL CHPR2( UPLO, I, -ONE, AP( I1 ), 1, TAU, 1, AP )
251 AP( I1+I-1 ) = E( I )
252 D( I+1 ) = AP( I1+I )
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).
263 AP( 1 ) = REAL( AP( 1 ) )
265 I1I1 = II + N - I + 1
267 * Generate elementary reflector H(i) = I - tau * v * v**H
268 * to annihilate A(i+2:n,i)
271 CALL CLARFG( N-I, ALPHA, AP( II+2 ), 1, TAUI )
274 IF( TAUI.NE.ZERO ) THEN
276 * Apply H(i) from both sides to A(i+1:n,i+1:n)
280 * Compute y := tau * A * v storing y in TAU(i:n-1)
282 CALL CHPMV( UPLO, N-I, TAUI, AP( I1I1 ), AP( II+1 ), 1,
283 $ ZERO, TAU( I ), 1 )
285 * Compute w := y - 1/2 * tau * (y**H *v) * v
287 ALPHA = -HALF*TAUI*CDOTC( N-I, TAU( I ), 1, AP( II+1 ),
289 CALL CAXPY( N-I, ALPHA, AP( II+1 ), 1, TAU( I ), 1 )
291 * Apply the transformation as a rank-2 update:
292 * A := A - v * w**H - w * v**H
294 CALL CHPR2( UPLO, N-I, -ONE, AP( II+1 ), 1, TAU( I ), 1,