1*> \brief <b> ZHPEVD computes the eigenvalues and, optionally, the left and/or right eigenvectors for OTHER matrices</b> 2* 3* =========== DOCUMENTATION =========== 4* 5* Online html documentation available at 6* http://www.netlib.org/lapack/explore-html/ 7* 8*> \htmlonly 9*> Download ZHPEVD + dependencies 10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zhpevd.f"> 11*> [TGZ]</a> 12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zhpevd.f"> 13*> [ZIP]</a> 14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zhpevd.f"> 15*> [TXT]</a> 16*> \endhtmlonly 17* 18* Definition: 19* =========== 20* 21* SUBROUTINE ZHPEVD( JOBZ, UPLO, N, AP, W, Z, LDZ, WORK, LWORK, 22* RWORK, LRWORK, IWORK, LIWORK, INFO ) 23* 24* .. Scalar Arguments .. 25* CHARACTER JOBZ, UPLO 26* INTEGER INFO, LDZ, LIWORK, LRWORK, LWORK, N 27* .. 28* .. Array Arguments .. 29* INTEGER IWORK( * ) 30* DOUBLE PRECISION RWORK( * ), W( * ) 31* COMPLEX*16 AP( * ), WORK( * ), Z( LDZ, * ) 32* .. 33* 34* 35*> \par Purpose: 36* ============= 37*> 38*> \verbatim 39*> 40*> ZHPEVD computes all the eigenvalues and, optionally, eigenvectors of 41*> a complex Hermitian matrix A in packed storage. If eigenvectors are 42*> desired, it uses a divide and conquer algorithm. 43*> 44*> The divide and conquer algorithm makes very mild assumptions about 45*> floating point arithmetic. It will work on machines with a guard 46*> digit in add/subtract, or on those binary machines without guard 47*> digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or 48*> Cray-2. It could conceivably fail on hexadecimal or decimal machines 49*> without guard digits, but we know of none. 50*> \endverbatim 51* 52* Arguments: 53* ========== 54* 55*> \param[in] JOBZ 56*> \verbatim 57*> JOBZ is CHARACTER*1 58*> = 'N': Compute eigenvalues only; 59*> = 'V': Compute eigenvalues and eigenvectors. 60*> \endverbatim 61*> 62*> \param[in] UPLO 63*> \verbatim 64*> UPLO is CHARACTER*1 65*> = 'U': Upper triangle of A is stored; 66*> = 'L': Lower triangle of A is stored. 67*> \endverbatim 68*> 69*> \param[in] N 70*> \verbatim 71*> N is INTEGER 72*> The order of the matrix A. N >= 0. 73*> \endverbatim 74*> 75*> \param[in,out] AP 76*> \verbatim 77*> AP is COMPLEX*16 array, dimension (N*(N+1)/2) 78*> On entry, the upper or lower triangle of the Hermitian matrix 79*> A, packed columnwise in a linear array. The j-th column of A 80*> is stored in the array AP as follows: 81*> if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j; 82*> if UPLO = 'L', AP(i + (j-1)*(2*n-j)/2) = A(i,j) for j<=i<=n. 83*> 84*> On exit, AP is overwritten by values generated during the 85*> reduction to tridiagonal form. If UPLO = 'U', the diagonal 86*> and first superdiagonal of the tridiagonal matrix T overwrite 87*> the corresponding elements of A, and if UPLO = 'L', the 88*> diagonal and first subdiagonal of T overwrite the 89*> corresponding elements of A. 90*> \endverbatim 91*> 92*> \param[out] W 93*> \verbatim 94*> W is DOUBLE PRECISION array, dimension (N) 95*> If INFO = 0, the eigenvalues in ascending order. 96*> \endverbatim 97*> 98*> \param[out] Z 99*> \verbatim 100*> Z is COMPLEX*16 array, dimension (LDZ, N) 101*> If JOBZ = 'V', then if INFO = 0, Z contains the orthonormal 102*> eigenvectors of the matrix A, with the i-th column of Z 103*> holding the eigenvector associated with W(i). 104*> If JOBZ = 'N', then Z is not referenced. 105*> \endverbatim 106*> 107*> \param[in] LDZ 108*> \verbatim 109*> LDZ is INTEGER 110*> The leading dimension of the array Z. LDZ >= 1, and if 111*> JOBZ = 'V', LDZ >= max(1,N). 112*> \endverbatim 113*> 114*> \param[out] WORK 115*> \verbatim 116*> WORK is COMPLEX*16 array, dimension (MAX(1,LWORK)) 117*> On exit, if INFO = 0, WORK(1) returns the required LWORK. 118*> \endverbatim 119*> 120*> \param[in] LWORK 121*> \verbatim 122*> LWORK is INTEGER 123*> The dimension of array WORK. 124*> If N <= 1, LWORK must be at least 1. 125*> If JOBZ = 'N' and N > 1, LWORK must be at least N. 126*> If JOBZ = 'V' and N > 1, LWORK must be at least 2*N. 127*> 128*> If LWORK = -1, then a workspace query is assumed; the routine 129*> only calculates the required sizes of the WORK, RWORK and 130*> IWORK arrays, returns these values as the first entries of 131*> the WORK, RWORK and IWORK arrays, and no error message 132*> related to LWORK or LRWORK or LIWORK is issued by XERBLA. 133*> \endverbatim 134*> 135*> \param[out] RWORK 136*> \verbatim 137*> RWORK is DOUBLE PRECISION array, dimension (MAX(1,LRWORK)) 138*> On exit, if INFO = 0, RWORK(1) returns the required LRWORK. 139*> \endverbatim 140*> 141*> \param[in] LRWORK 142*> \verbatim 143*> LRWORK is INTEGER 144*> The dimension of array RWORK. 145*> If N <= 1, LRWORK must be at least 1. 146*> If JOBZ = 'N' and N > 1, LRWORK must be at least N. 147*> If JOBZ = 'V' and N > 1, LRWORK must be at least 148*> 1 + 5*N + 2*N**2. 149*> 150*> If LRWORK = -1, then a workspace query is assumed; the 151*> routine only calculates the required sizes of the WORK, RWORK 152*> and IWORK arrays, returns these values as the first entries 153*> of the WORK, RWORK and IWORK arrays, and no error message 154*> related to LWORK or LRWORK or LIWORK is issued by XERBLA. 155*> \endverbatim 156*> 157*> \param[out] IWORK 158*> \verbatim 159*> IWORK is INTEGER array, dimension (MAX(1,LIWORK)) 160*> On exit, if INFO = 0, IWORK(1) returns the required LIWORK. 161*> \endverbatim 162*> 163*> \param[in] LIWORK 164*> \verbatim 165*> LIWORK is INTEGER 166*> The dimension of array IWORK. 167*> If JOBZ = 'N' or N <= 1, LIWORK must be at least 1. 168*> If JOBZ = 'V' and N > 1, LIWORK must be at least 3 + 5*N. 169*> 170*> If LIWORK = -1, then a workspace query is assumed; the 171*> routine only calculates the required sizes of the WORK, RWORK 172*> and IWORK arrays, returns these values as the first entries 173*> of the WORK, RWORK and IWORK arrays, and no error message 174*> related to LWORK or LRWORK or LIWORK is issued by XERBLA. 175*> \endverbatim 176*> 177*> \param[out] INFO 178*> \verbatim 179*> INFO is INTEGER 180*> = 0: successful exit 181*> < 0: if INFO = -i, the i-th argument had an illegal value. 182*> > 0: if INFO = i, the algorithm failed to converge; i 183*> off-diagonal elements of an intermediate tridiagonal 184*> form did not converge to zero. 185*> \endverbatim 186* 187* Authors: 188* ======== 189* 190*> \author Univ. of Tennessee 191*> \author Univ. of California Berkeley 192*> \author Univ. of Colorado Denver 193*> \author NAG Ltd. 194* 195*> \ingroup complex16OTHEReigen 196* 197* ===================================================================== 198 SUBROUTINE ZHPEVD( JOBZ, UPLO, N, AP, W, Z, LDZ, WORK, LWORK, 199 $ RWORK, LRWORK, IWORK, LIWORK, INFO ) 200* 201* -- LAPACK driver routine -- 202* -- LAPACK is a software package provided by Univ. of Tennessee, -- 203* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 204* 205* .. Scalar Arguments .. 206 CHARACTER JOBZ, UPLO 207 INTEGER INFO, LDZ, LIWORK, LRWORK, LWORK, N 208* .. 209* .. Array Arguments .. 210 INTEGER IWORK( * ) 211 DOUBLE PRECISION RWORK( * ), W( * ) 212 COMPLEX*16 AP( * ), WORK( * ), Z( LDZ, * ) 213* .. 214* 215* ===================================================================== 216* 217* .. Parameters .. 218 DOUBLE PRECISION ZERO, ONE 219 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) 220 COMPLEX*16 CONE 221 PARAMETER ( CONE = ( 1.0D+0, 0.0D+0 ) ) 222* .. 223* .. Local Scalars .. 224 LOGICAL LQUERY, WANTZ 225 INTEGER IINFO, IMAX, INDE, INDRWK, INDTAU, INDWRK, 226 $ ISCALE, LIWMIN, LLRWK, LLWRK, LRWMIN, LWMIN 227 DOUBLE PRECISION ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA, 228 $ SMLNUM 229* .. 230* .. External Functions .. 231 LOGICAL LSAME 232 DOUBLE PRECISION DLAMCH, ZLANHP 233 EXTERNAL LSAME, DLAMCH, ZLANHP 234* .. 235* .. External Subroutines .. 236 EXTERNAL DSCAL, DSTERF, XERBLA, ZDSCAL, ZHPTRD, ZSTEDC, 237 $ ZUPMTR 238* .. 239* .. Intrinsic Functions .. 240 INTRINSIC SQRT 241* .. 242* .. Executable Statements .. 243* 244* Test the input parameters. 245* 246 WANTZ = LSAME( JOBZ, 'V' ) 247 LQUERY = ( LWORK.EQ.-1 .OR. LRWORK.EQ.-1 .OR. LIWORK.EQ.-1 ) 248* 249 INFO = 0 250 IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN 251 INFO = -1 252 ELSE IF( .NOT.( LSAME( UPLO, 'L' ) .OR. LSAME( UPLO, 'U' ) ) ) 253 $ THEN 254 INFO = -2 255 ELSE IF( N.LT.0 ) THEN 256 INFO = -3 257 ELSE IF( LDZ.LT.1 .OR. ( WANTZ .AND. LDZ.LT.N ) ) THEN 258 INFO = -7 259 END IF 260* 261 IF( INFO.EQ.0 ) THEN 262 IF( N.LE.1 ) THEN 263 LWMIN = 1 264 LIWMIN = 1 265 LRWMIN = 1 266 ELSE 267 IF( WANTZ ) THEN 268 LWMIN = 2*N 269 LRWMIN = 1 + 5*N + 2*N**2 270 LIWMIN = 3 + 5*N 271 ELSE 272 LWMIN = N 273 LRWMIN = N 274 LIWMIN = 1 275 END IF 276 END IF 277 WORK( 1 ) = LWMIN 278 RWORK( 1 ) = LRWMIN 279 IWORK( 1 ) = LIWMIN 280* 281 IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN 282 INFO = -9 283 ELSE IF( LRWORK.LT.LRWMIN .AND. .NOT.LQUERY ) THEN 284 INFO = -11 285 ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT.LQUERY ) THEN 286 INFO = -13 287 END IF 288 END IF 289* 290 IF( INFO.NE.0 ) THEN 291 CALL XERBLA( 'ZHPEVD', -INFO ) 292 RETURN 293 ELSE IF( LQUERY ) THEN 294 RETURN 295 END IF 296* 297* Quick return if possible 298* 299 IF( N.EQ.0 ) 300 $ RETURN 301* 302 IF( N.EQ.1 ) THEN 303 W( 1 ) = DBLE( AP( 1 ) ) 304 IF( WANTZ ) 305 $ Z( 1, 1 ) = CONE 306 RETURN 307 END IF 308* 309* Get machine constants. 310* 311 SAFMIN = DLAMCH( 'Safe minimum' ) 312 EPS = DLAMCH( 'Precision' ) 313 SMLNUM = SAFMIN / EPS 314 BIGNUM = ONE / SMLNUM 315 RMIN = SQRT( SMLNUM ) 316 RMAX = SQRT( BIGNUM ) 317* 318* Scale matrix to allowable range, if necessary. 319* 320 ANRM = ZLANHP( 'M', UPLO, N, AP, RWORK ) 321 ISCALE = 0 322 IF( ANRM.GT.ZERO .AND. ANRM.LT.RMIN ) THEN 323 ISCALE = 1 324 SIGMA = RMIN / ANRM 325 ELSE IF( ANRM.GT.RMAX ) THEN 326 ISCALE = 1 327 SIGMA = RMAX / ANRM 328 END IF 329 IF( ISCALE.EQ.1 ) THEN 330 CALL ZDSCAL( ( N*( N+1 ) ) / 2, SIGMA, AP, 1 ) 331 END IF 332* 333* Call ZHPTRD to reduce Hermitian packed matrix to tridiagonal form. 334* 335 INDE = 1 336 INDTAU = 1 337 INDRWK = INDE + N 338 INDWRK = INDTAU + N 339 LLWRK = LWORK - INDWRK + 1 340 LLRWK = LRWORK - INDRWK + 1 341 CALL ZHPTRD( UPLO, N, AP, W, RWORK( INDE ), WORK( INDTAU ), 342 $ IINFO ) 343* 344* For eigenvalues only, call DSTERF. For eigenvectors, first call 345* ZUPGTR to generate the orthogonal matrix, then call ZSTEDC. 346* 347 IF( .NOT.WANTZ ) THEN 348 CALL DSTERF( N, W, RWORK( INDE ), INFO ) 349 ELSE 350 CALL ZSTEDC( 'I', N, W, RWORK( INDE ), Z, LDZ, WORK( INDWRK ), 351 $ LLWRK, RWORK( INDRWK ), LLRWK, IWORK, LIWORK, 352 $ INFO ) 353 CALL ZUPMTR( 'L', UPLO, 'N', N, N, AP, WORK( INDTAU ), Z, LDZ, 354 $ WORK( INDWRK ), IINFO ) 355 END IF 356* 357* If matrix was scaled, then rescale eigenvalues appropriately. 358* 359 IF( ISCALE.EQ.1 ) THEN 360 IF( INFO.EQ.0 ) THEN 361 IMAX = N 362 ELSE 363 IMAX = INFO - 1 364 END IF 365 CALL DSCAL( IMAX, ONE / SIGMA, W, 1 ) 366 END IF 367* 368 WORK( 1 ) = LWMIN 369 RWORK( 1 ) = LRWMIN 370 IWORK( 1 ) = LIWMIN 371 RETURN 372* 373* End of ZHPEVD 374* 375 END 376