--- rpl/lapack/lapack/zhptrd.f 2010/08/06 15:28:55 1.3 +++ rpl/lapack/lapack/zhptrd.f 2017/06/17 10:54:16 1.15 @@ -1,9 +1,160 @@ +*> \brief \b ZHPTRD +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +*> \htmlonly +*> Download ZHPTRD + dependencies +*> +*> [TGZ] +*> +*> [ZIP] +*> +*> [TXT] +*> \endhtmlonly +* +* Definition: +* =========== +* +* SUBROUTINE ZHPTRD( UPLO, N, AP, D, E, TAU, INFO ) +* +* .. Scalar Arguments .. +* CHARACTER UPLO +* INTEGER INFO, N +* .. +* .. Array Arguments .. +* DOUBLE PRECISION D( * ), E( * ) +* COMPLEX*16 AP( * ), TAU( * ) +* .. +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> ZHPTRD reduces a complex Hermitian matrix A stored in packed form to +*> real symmetric tridiagonal form T by a unitary similarity +*> transformation: Q**H * A * Q = T. +*> \endverbatim +* +* Arguments: +* ========== +* +*> \param[in] UPLO +*> \verbatim +*> UPLO is CHARACTER*1 +*> = 'U': Upper triangle of A is stored; +*> = 'L': Lower triangle of A is stored. +*> \endverbatim +*> +*> \param[in] N +*> \verbatim +*> N is INTEGER +*> The order of the matrix A. N >= 0. +*> \endverbatim +*> +*> \param[in,out] AP +*> \verbatim +*> AP is COMPLEX*16 array, dimension (N*(N+1)/2) +*> On entry, the upper or lower triangle of the Hermitian matrix +*> A, packed columnwise in a linear array. The j-th column of A +*> is stored in the array AP as follows: +*> if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j; +*> if UPLO = 'L', AP(i + (j-1)*(2*n-j)/2) = A(i,j) for j<=i<=n. +*> On exit, if UPLO = 'U', the diagonal and first superdiagonal +*> of A are overwritten by the corresponding elements of the +*> tridiagonal matrix T, and the elements above the first +*> superdiagonal, with the array TAU, represent the unitary +*> matrix Q as a product of elementary reflectors; if UPLO +*> = 'L', the diagonal and first subdiagonal of A are over- +*> written by the corresponding elements of the tridiagonal +*> matrix T, and the elements below the first subdiagonal, with +*> the array TAU, represent the unitary matrix Q as a product +*> of elementary reflectors. See Further Details. +*> \endverbatim +*> +*> \param[out] D +*> \verbatim +*> D is DOUBLE PRECISION array, dimension (N) +*> The diagonal elements of the tridiagonal matrix T: +*> D(i) = A(i,i). +*> \endverbatim +*> +*> \param[out] E +*> \verbatim +*> E is DOUBLE PRECISION array, dimension (N-1) +*> The off-diagonal elements of the tridiagonal matrix T: +*> E(i) = A(i,i+1) if UPLO = 'U', E(i) = A(i+1,i) if UPLO = 'L'. +*> \endverbatim +*> +*> \param[out] TAU +*> \verbatim +*> TAU is COMPLEX*16 array, dimension (N-1) +*> The scalar factors of the elementary reflectors (see Further +*> Details). +*> \endverbatim +*> +*> \param[out] INFO +*> \verbatim +*> INFO is INTEGER +*> = 0: successful exit +*> < 0: if INFO = -i, the i-th argument had an illegal value +*> \endverbatim +* +* Authors: +* ======== +* +*> \author Univ. of Tennessee +*> \author Univ. of California Berkeley +*> \author Univ. of Colorado Denver +*> \author NAG Ltd. +* +*> \date December 2016 +* +*> \ingroup complex16OTHERcomputational +* +*> \par Further Details: +* ===================== +*> +*> \verbatim +*> +*> If UPLO = 'U', the matrix Q is represented as a product of elementary +*> reflectors +*> +*> Q = H(n-1) . . . H(2) H(1). +*> +*> Each H(i) has the form +*> +*> H(i) = I - tau * v * v**H +*> +*> where tau is a complex scalar, and v is a complex vector with +*> v(i+1:n) = 0 and v(i) = 1; v(1:i-1) is stored on exit in AP, +*> overwriting A(1:i-1,i+1), and tau is stored in TAU(i). +*> +*> If UPLO = 'L', the matrix Q is represented as a product of elementary +*> reflectors +*> +*> Q = H(1) H(2) . . . H(n-1). +*> +*> Each H(i) has the form +*> +*> H(i) = I - tau * v * v**H +*> +*> where tau is a complex scalar, and v is a complex vector with +*> v(1:i) = 0 and v(i+1) = 1; v(i+2:n) is stored on exit in AP, +*> overwriting A(i+2:n,i), and tau is stored in TAU(i). +*> \endverbatim +*> +* ===================================================================== SUBROUTINE ZHPTRD( UPLO, N, AP, D, E, TAU, INFO ) * -* -- LAPACK routine (version 3.2) -- +* -- LAPACK computational routine (version 3.7.0) -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- -* November 2006 +* December 2016 * * .. Scalar Arguments .. CHARACTER UPLO @@ -14,85 +165,6 @@ COMPLEX*16 AP( * ), TAU( * ) * .. * -* Purpose -* ======= -* -* ZHPTRD reduces a complex Hermitian matrix A stored in packed form to -* real symmetric tridiagonal form T by a unitary similarity -* transformation: Q**H * A * Q = T. -* -* Arguments -* ========= -* -* UPLO (input) CHARACTER*1 -* = 'U': Upper triangle of A is stored; -* = 'L': Lower triangle of A is stored. -* -* N (input) INTEGER -* The order of the matrix A. N >= 0. -* -* AP (input/output) COMPLEX*16 array, dimension (N*(N+1)/2) -* On entry, the upper or lower triangle of the Hermitian matrix -* A, packed columnwise in a linear array. The j-th column of A -* is stored in the array AP as follows: -* if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j; -* if UPLO = 'L', AP(i + (j-1)*(2*n-j)/2) = A(i,j) for j<=i<=n. -* On exit, if UPLO = 'U', the diagonal and first superdiagonal -* of A are overwritten by the corresponding elements of the -* tridiagonal matrix T, and the elements above the first -* superdiagonal, with the array TAU, represent the unitary -* matrix Q as a product of elementary reflectors; if UPLO -* = 'L', the diagonal and first subdiagonal of A are over- -* written by the corresponding elements of the tridiagonal -* matrix T, and the elements below the first subdiagonal, with -* the array TAU, represent the unitary matrix Q as a product -* of elementary reflectors. See Further Details. -* -* D (output) DOUBLE PRECISION array, dimension (N) -* The diagonal elements of the tridiagonal matrix T: -* D(i) = A(i,i). -* -* E (output) DOUBLE PRECISION array, dimension (N-1) -* The off-diagonal elements of the tridiagonal matrix T: -* E(i) = A(i,i+1) if UPLO = 'U', E(i) = A(i+1,i) if UPLO = 'L'. -* -* TAU (output) COMPLEX*16 array, dimension (N-1) -* The scalar factors of the elementary reflectors (see Further -* Details). -* -* INFO (output) INTEGER -* = 0: successful exit -* < 0: if INFO = -i, the i-th argument had an illegal value -* -* Further Details -* =============== -* -* If UPLO = 'U', the matrix Q is represented as a product of elementary -* reflectors -* -* Q = H(n-1) . . . H(2) H(1). -* -* Each H(i) has the form -* -* H(i) = I - tau * v * v' -* -* where tau is a complex scalar, and v is a complex vector with -* v(i+1:n) = 0 and v(i) = 1; v(1:i-1) is stored on exit in AP, -* overwriting A(1:i-1,i+1), and tau is stored in TAU(i). -* -* If UPLO = 'L', the matrix Q is represented as a product of elementary -* reflectors -* -* Q = H(1) H(2) . . . H(n-1). -* -* Each H(i) has the form -* -* H(i) = I - tau * v * v' -* -* where tau is a complex scalar, and v is a complex vector with -* v(1:i) = 0 and v(i+1) = 1; v(i+2:n) is stored on exit in AP, -* overwriting A(i+2:n,i), and tau is stored in TAU(i). -* * ===================================================================== * * .. Parameters .. @@ -147,7 +219,7 @@ AP( I1+N-1 ) = DBLE( AP( I1+N-1 ) ) DO 10 I = N - 1, 1, -1 * -* Generate elementary reflector H(i) = I - tau * v * v' +* Generate elementary reflector H(i) = I - tau * v * v**H * to annihilate A(1:i-1,i+1) * ALPHA = AP( I1+I-1 ) @@ -165,13 +237,13 @@ CALL ZHPMV( UPLO, I, TAUI, AP, AP( I1 ), 1, ZERO, TAU, $ 1 ) * -* Compute w := y - 1/2 * tau * (y'*v) * v +* Compute w := y - 1/2 * tau * (y**H *v) * v * ALPHA = -HALF*TAUI*ZDOTC( I, TAU, 1, AP( I1 ), 1 ) CALL ZAXPY( I, ALPHA, AP( I1 ), 1, TAU, 1 ) * * Apply the transformation as a rank-2 update: -* A := A - v * w' - w * v' +* A := A - v * w**H - w * v**H * CALL ZHPR2( UPLO, I, -ONE, AP( I1 ), 1, TAU, 1, AP ) * @@ -192,7 +264,7 @@ DO 20 I = 1, N - 1 I1I1 = II + N - I + 1 * -* Generate elementary reflector H(i) = I - tau * v * v' +* Generate elementary reflector H(i) = I - tau * v * v**H * to annihilate A(i+2:n,i) * ALPHA = AP( II+1 ) @@ -210,14 +282,14 @@ CALL ZHPMV( UPLO, N-I, TAUI, AP( I1I1 ), AP( II+1 ), 1, $ ZERO, TAU( I ), 1 ) * -* Compute w := y - 1/2 * tau * (y'*v) * v +* Compute w := y - 1/2 * tau * (y**H *v) * v * ALPHA = -HALF*TAUI*ZDOTC( N-I, TAU( I ), 1, AP( II+1 ), $ 1 ) CALL ZAXPY( N-I, ALPHA, AP( II+1 ), 1, TAU( I ), 1 ) * * Apply the transformation as a rank-2 update: -* A := A - v * w' - w * v' +* A := A - v * w**H - w * v**H * CALL ZHPR2( UPLO, N-I, -ONE, AP( II+1 ), 1, TAU( I ), 1, $ AP( I1I1 ) )