![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b ZHETD2 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download ZHETD2 + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zhetd2.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zhetd2.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zhetd2.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE ZHETD2( UPLO, N, A, LDA, D, E, TAU, INFO ) 00022 * 00023 * .. Scalar Arguments .. 00024 * CHARACTER UPLO 00025 * INTEGER INFO, LDA, N 00026 * .. 00027 * .. Array Arguments .. 00028 * DOUBLE PRECISION D( * ), E( * ) 00029 * COMPLEX*16 A( LDA, * ), TAU( * ) 00030 * .. 00031 * 00032 * 00033 *> \par Purpose: 00034 * ============= 00035 *> 00036 *> \verbatim 00037 *> 00038 *> ZHETD2 reduces a complex Hermitian matrix A to real symmetric 00039 *> tridiagonal form T by a unitary similarity transformation: 00040 *> Q**H * A * Q = T. 00041 *> \endverbatim 00042 * 00043 * Arguments: 00044 * ========== 00045 * 00046 *> \param[in] UPLO 00047 *> \verbatim 00048 *> UPLO is CHARACTER*1 00049 *> Specifies whether the upper or lower triangular part of the 00050 *> Hermitian matrix A is stored: 00051 *> = 'U': Upper triangular 00052 *> = 'L': Lower triangular 00053 *> \endverbatim 00054 *> 00055 *> \param[in] N 00056 *> \verbatim 00057 *> N is INTEGER 00058 *> The order of the matrix A. N >= 0. 00059 *> \endverbatim 00060 *> 00061 *> \param[in,out] A 00062 *> \verbatim 00063 *> A is COMPLEX*16 array, dimension (LDA,N) 00064 *> On entry, the Hermitian matrix A. If UPLO = 'U', the leading 00065 *> n-by-n upper triangular part of A contains the upper 00066 *> triangular part of the matrix A, and the strictly lower 00067 *> triangular part of A is not referenced. If UPLO = 'L', the 00068 *> leading n-by-n lower triangular part of A contains the lower 00069 *> triangular part of the matrix A, and the strictly upper 00070 *> triangular part of A is not referenced. 00071 *> On exit, if UPLO = 'U', the diagonal and first superdiagonal 00072 *> of A are overwritten by the corresponding elements of the 00073 *> tridiagonal matrix T, and the elements above the first 00074 *> superdiagonal, with the array TAU, represent the unitary 00075 *> matrix Q as a product of elementary reflectors; if UPLO 00076 *> = 'L', the diagonal and first subdiagonal of A are over- 00077 *> written by the corresponding elements of the tridiagonal 00078 *> matrix T, and the elements below the first subdiagonal, with 00079 *> the array TAU, represent the unitary matrix Q as a product 00080 *> of elementary reflectors. See Further Details. 00081 *> \endverbatim 00082 *> 00083 *> \param[in] LDA 00084 *> \verbatim 00085 *> LDA is INTEGER 00086 *> The leading dimension of the array A. LDA >= max(1,N). 00087 *> \endverbatim 00088 *> 00089 *> \param[out] D 00090 *> \verbatim 00091 *> D is DOUBLE PRECISION array, dimension (N) 00092 *> The diagonal elements of the tridiagonal matrix T: 00093 *> D(i) = A(i,i). 00094 *> \endverbatim 00095 *> 00096 *> \param[out] E 00097 *> \verbatim 00098 *> E is DOUBLE PRECISION array, dimension (N-1) 00099 *> The off-diagonal elements of the tridiagonal matrix T: 00100 *> E(i) = A(i,i+1) if UPLO = 'U', E(i) = A(i+1,i) if UPLO = 'L'. 00101 *> \endverbatim 00102 *> 00103 *> \param[out] TAU 00104 *> \verbatim 00105 *> TAU is COMPLEX*16 array, dimension (N-1) 00106 *> The scalar factors of the elementary reflectors (see Further 00107 *> Details). 00108 *> \endverbatim 00109 *> 00110 *> \param[out] INFO 00111 *> \verbatim 00112 *> INFO is INTEGER 00113 *> = 0: successful exit 00114 *> < 0: if INFO = -i, the i-th argument had an illegal value. 00115 *> \endverbatim 00116 * 00117 * Authors: 00118 * ======== 00119 * 00120 *> \author Univ. of Tennessee 00121 *> \author Univ. of California Berkeley 00122 *> \author Univ. of Colorado Denver 00123 *> \author NAG Ltd. 00124 * 00125 *> \date November 2011 00126 * 00127 *> \ingroup complex16HEcomputational 00128 * 00129 *> \par Further Details: 00130 * ===================== 00131 *> 00132 *> \verbatim 00133 *> 00134 *> If UPLO = 'U', the matrix Q is represented as a product of elementary 00135 *> reflectors 00136 *> 00137 *> Q = H(n-1) . . . H(2) H(1). 00138 *> 00139 *> Each H(i) has the form 00140 *> 00141 *> H(i) = I - tau * v * v**H 00142 *> 00143 *> where tau is a complex scalar, and v is a complex vector with 00144 *> v(i+1:n) = 0 and v(i) = 1; v(1:i-1) is stored on exit in 00145 *> A(1:i-1,i+1), and tau in TAU(i). 00146 *> 00147 *> If UPLO = 'L', the matrix Q is represented as a product of elementary 00148 *> reflectors 00149 *> 00150 *> Q = H(1) H(2) . . . H(n-1). 00151 *> 00152 *> Each H(i) has the form 00153 *> 00154 *> H(i) = I - tau * v * v**H 00155 *> 00156 *> where tau is a complex scalar, and v is a complex vector with 00157 *> v(1:i) = 0 and v(i+1) = 1; v(i+2:n) is stored on exit in A(i+2:n,i), 00158 *> and tau in TAU(i). 00159 *> 00160 *> The contents of A on exit are illustrated by the following examples 00161 *> with n = 5: 00162 *> 00163 *> if UPLO = 'U': if UPLO = 'L': 00164 *> 00165 *> ( d e v2 v3 v4 ) ( d ) 00166 *> ( d e v3 v4 ) ( e d ) 00167 *> ( d e v4 ) ( v1 e d ) 00168 *> ( d e ) ( v1 v2 e d ) 00169 *> ( d ) ( v1 v2 v3 e d ) 00170 *> 00171 *> where d and e denote diagonal and off-diagonal elements of T, and vi 00172 *> denotes an element of the vector defining H(i). 00173 *> \endverbatim 00174 *> 00175 * ===================================================================== 00176 SUBROUTINE ZHETD2( UPLO, N, A, LDA, D, E, TAU, INFO ) 00177 * 00178 * -- LAPACK computational routine (version 3.4.0) -- 00179 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00180 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00181 * November 2011 00182 * 00183 * .. Scalar Arguments .. 00184 CHARACTER UPLO 00185 INTEGER INFO, LDA, N 00186 * .. 00187 * .. Array Arguments .. 00188 DOUBLE PRECISION D( * ), E( * ) 00189 COMPLEX*16 A( LDA, * ), TAU( * ) 00190 * .. 00191 * 00192 * ===================================================================== 00193 * 00194 * .. Parameters .. 00195 COMPLEX*16 ONE, ZERO, HALF 00196 PARAMETER ( ONE = ( 1.0D+0, 0.0D+0 ), 00197 $ ZERO = ( 0.0D+0, 0.0D+0 ), 00198 $ HALF = ( 0.5D+0, 0.0D+0 ) ) 00199 * .. 00200 * .. Local Scalars .. 00201 LOGICAL UPPER 00202 INTEGER I 00203 COMPLEX*16 ALPHA, TAUI 00204 * .. 00205 * .. External Subroutines .. 00206 EXTERNAL XERBLA, ZAXPY, ZHEMV, ZHER2, ZLARFG 00207 * .. 00208 * .. External Functions .. 00209 LOGICAL LSAME 00210 COMPLEX*16 ZDOTC 00211 EXTERNAL LSAME, ZDOTC 00212 * .. 00213 * .. Intrinsic Functions .. 00214 INTRINSIC DBLE, MAX, MIN 00215 * .. 00216 * .. Executable Statements .. 00217 * 00218 * Test the input parameters 00219 * 00220 INFO = 0 00221 UPPER = LSAME( UPLO, 'U') 00222 IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN 00223 INFO = -1 00224 ELSE IF( N.LT.0 ) THEN 00225 INFO = -2 00226 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 00227 INFO = -4 00228 END IF 00229 IF( INFO.NE.0 ) THEN 00230 CALL XERBLA( 'ZHETD2', -INFO ) 00231 RETURN 00232 END IF 00233 * 00234 * Quick return if possible 00235 * 00236 IF( N.LE.0 ) 00237 $ RETURN 00238 * 00239 IF( UPPER ) THEN 00240 * 00241 * Reduce the upper triangle of A 00242 * 00243 A( N, N ) = DBLE( A( N, N ) ) 00244 DO 10 I = N - 1, 1, -1 00245 * 00246 * Generate elementary reflector H(i) = I - tau * v * v**H 00247 * to annihilate A(1:i-1,i+1) 00248 * 00249 ALPHA = A( I, I+1 ) 00250 CALL ZLARFG( I, ALPHA, A( 1, I+1 ), 1, TAUI ) 00251 E( I ) = ALPHA 00252 * 00253 IF( TAUI.NE.ZERO ) THEN 00254 * 00255 * Apply H(i) from both sides to A(1:i,1:i) 00256 * 00257 A( I, I+1 ) = ONE 00258 * 00259 * Compute x := tau * A * v storing x in TAU(1:i) 00260 * 00261 CALL ZHEMV( UPLO, I, TAUI, A, LDA, A( 1, I+1 ), 1, ZERO, 00262 $ TAU, 1 ) 00263 * 00264 * Compute w := x - 1/2 * tau * (x**H * v) * v 00265 * 00266 ALPHA = -HALF*TAUI*ZDOTC( I, TAU, 1, A( 1, I+1 ), 1 ) 00267 CALL ZAXPY( I, ALPHA, A( 1, I+1 ), 1, TAU, 1 ) 00268 * 00269 * Apply the transformation as a rank-2 update: 00270 * A := A - v * w**H - w * v**H 00271 * 00272 CALL ZHER2( UPLO, I, -ONE, A( 1, I+1 ), 1, TAU, 1, A, 00273 $ LDA ) 00274 * 00275 ELSE 00276 A( I, I ) = DBLE( A( I, I ) ) 00277 END IF 00278 A( I, I+1 ) = E( I ) 00279 D( I+1 ) = A( I+1, I+1 ) 00280 TAU( I ) = TAUI 00281 10 CONTINUE 00282 D( 1 ) = A( 1, 1 ) 00283 ELSE 00284 * 00285 * Reduce the lower triangle of A 00286 * 00287 A( 1, 1 ) = DBLE( A( 1, 1 ) ) 00288 DO 20 I = 1, N - 1 00289 * 00290 * Generate elementary reflector H(i) = I - tau * v * v**H 00291 * to annihilate A(i+2:n,i) 00292 * 00293 ALPHA = A( I+1, I ) 00294 CALL ZLARFG( N-I, ALPHA, A( MIN( I+2, N ), I ), 1, TAUI ) 00295 E( I ) = ALPHA 00296 * 00297 IF( TAUI.NE.ZERO ) THEN 00298 * 00299 * Apply H(i) from both sides to A(i+1:n,i+1:n) 00300 * 00301 A( I+1, I ) = ONE 00302 * 00303 * Compute x := tau * A * v storing y in TAU(i:n-1) 00304 * 00305 CALL ZHEMV( UPLO, N-I, TAUI, A( I+1, I+1 ), LDA, 00306 $ A( I+1, I ), 1, ZERO, TAU( I ), 1 ) 00307 * 00308 * Compute w := x - 1/2 * tau * (x**H * v) * v 00309 * 00310 ALPHA = -HALF*TAUI*ZDOTC( N-I, TAU( I ), 1, A( I+1, I ), 00311 $ 1 ) 00312 CALL ZAXPY( N-I, ALPHA, A( I+1, I ), 1, TAU( I ), 1 ) 00313 * 00314 * Apply the transformation as a rank-2 update: 00315 * A := A - v * w**H - w * v**H 00316 * 00317 CALL ZHER2( UPLO, N-I, -ONE, A( I+1, I ), 1, TAU( I ), 1, 00318 $ A( I+1, I+1 ), LDA ) 00319 * 00320 ELSE 00321 A( I+1, I+1 ) = DBLE( A( I+1, I+1 ) ) 00322 END IF 00323 A( I+1, I ) = E( I ) 00324 D( I ) = A( I, I ) 00325 TAU( I ) = TAUI 00326 20 CONTINUE 00327 D( N ) = A( N, N ) 00328 END IF 00329 * 00330 RETURN 00331 * 00332 * End of ZHETD2 00333 * 00334 END