 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
00001 SUBROUTINE ZHPEVX( JOBZ, RANGE, UPLO, N, AP, VL, VU, IL, IU, 00002 $ ABSTOL, M, W, Z, LDZ, WORK, RWORK, IWORK, 00003 $ IFAIL, INFO ) 00004 * 00005 * -- LAPACK driver routine (version 3.2) -- 00006 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00007 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00008 * November 2006 00009 * 00010 * .. Scalar Arguments .. 00011 CHARACTER JOBZ, RANGE, UPLO 00012 INTEGER IL, INFO, IU, LDZ, M, N 00013 DOUBLE PRECISION ABSTOL, VL, VU 00014 * .. 00015 * .. Array Arguments .. 00016 INTEGER IFAIL( * ), IWORK( * ) 00017 DOUBLE PRECISION RWORK( * ), W( * ) 00018 COMPLEX*16 AP( * ), WORK( * ), Z( LDZ, * ) 00019 * .. 00020 * 00021 * Purpose 00022 * ======= 00023 * 00024 * ZHPEVX computes selected eigenvalues and, optionally, eigenvectors 00025 * of a complex Hermitian matrix A in packed storage. 00026 * Eigenvalues/vectors can be selected by specifying either a range of 00027 * values or a range of indices for the desired eigenvalues. 00028 * 00029 * Arguments 00030 * ========= 00031 * 00032 * JOBZ (input) CHARACTER*1 00033 * = 'N': Compute eigenvalues only; 00034 * = 'V': Compute eigenvalues and eigenvectors. 00035 * 00036 * RANGE (input) CHARACTER*1 00037 * = 'A': all eigenvalues will be found; 00038 * = 'V': all eigenvalues in the half-open interval (VL,VU] 00039 * will be found; 00040 * = 'I': the IL-th through IU-th eigenvalues will be found. 00041 * 00042 * UPLO (input) CHARACTER*1 00043 * = 'U': Upper triangle of A is stored; 00044 * = 'L': Lower triangle of A is stored. 00045 * 00046 * N (input) INTEGER 00047 * The order of the matrix A. N >= 0. 00048 * 00049 * AP (input/output) COMPLEX*16 array, dimension (N*(N+1)/2) 00050 * On entry, the upper or lower triangle of the Hermitian matrix 00051 * A, packed columnwise in a linear array. The j-th column of A 00052 * is stored in the array AP as follows: 00053 * if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j; 00054 * if UPLO = 'L', AP(i + (j-1)*(2*n-j)/2) = A(i,j) for j<=i<=n. 00055 * 00056 * On exit, AP is overwritten by values generated during the 00057 * reduction to tridiagonal form. If UPLO = 'U', the diagonal 00058 * and first superdiagonal of the tridiagonal matrix T overwrite 00059 * the corresponding elements of A, and if UPLO = 'L', the 00060 * diagonal and first subdiagonal of T overwrite the 00061 * corresponding elements of A. 00062 * 00063 * VL (input) DOUBLE PRECISION 00064 * VU (input) DOUBLE PRECISION 00065 * If RANGE='V', the lower and upper bounds of the interval to 00066 * be searched for eigenvalues. VL < VU. 00067 * Not referenced if RANGE = 'A' or 'I'. 00068 * 00069 * IL (input) INTEGER 00070 * IU (input) INTEGER 00071 * If RANGE='I', the indices (in ascending order) of the 00072 * smallest and largest eigenvalues to be returned. 00073 * 1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0. 00074 * Not referenced if RANGE = 'A' or 'V'. 00075 * 00076 * ABSTOL (input) DOUBLE PRECISION 00077 * The absolute error tolerance for the eigenvalues. 00078 * An approximate eigenvalue is accepted as converged 00079 * when it is determined to lie in an interval [a,b] 00080 * of width less than or equal to 00081 * 00082 * ABSTOL + EPS * max( |a|,|b| ) , 00083 * 00084 * where EPS is the machine precision. If ABSTOL is less than 00085 * or equal to zero, then EPS*|T| will be used in its place, 00086 * where |T| is the 1-norm of the tridiagonal matrix obtained 00087 * by reducing AP to tridiagonal form. 00088 * 00089 * Eigenvalues will be computed most accurately when ABSTOL is 00090 * set to twice the underflow threshold 2*DLAMCH('S'), not zero. 00091 * If this routine returns with INFO>0, indicating that some 00092 * eigenvectors did not converge, try setting ABSTOL to 00093 * 2*DLAMCH('S'). 00094 * 00095 * See "Computing Small Singular Values of Bidiagonal Matrices 00096 * with Guaranteed High Relative Accuracy," by Demmel and 00097 * Kahan, LAPACK Working Note #3. 00098 * 00099 * M (output) INTEGER 00100 * The total number of eigenvalues found. 0 <= M <= N. 00101 * If RANGE = 'A', M = N, and if RANGE = 'I', M = IU-IL+1. 00102 * 00103 * W (output) DOUBLE PRECISION array, dimension (N) 00104 * If INFO = 0, the selected eigenvalues in ascending order. 00105 * 00106 * Z (output) COMPLEX*16 array, dimension (LDZ, max(1,M)) 00107 * If JOBZ = 'V', then if INFO = 0, the first M columns of Z 00108 * contain the orthonormal eigenvectors of the matrix A 00109 * corresponding to the selected eigenvalues, with the i-th 00110 * column of Z holding the eigenvector associated with W(i). 00111 * If an eigenvector fails to converge, then that column of Z 00112 * contains the latest approximation to the eigenvector, and 00113 * the index of the eigenvector is returned in IFAIL. 00114 * If JOBZ = 'N', then Z is not referenced. 00115 * Note: the user must ensure that at least max(1,M) columns are 00116 * supplied in the array Z; if RANGE = 'V', the exact value of M 00117 * is not known in advance and an upper bound must be used. 00118 * 00119 * LDZ (input) INTEGER 00120 * The leading dimension of the array Z. LDZ >= 1, and if 00121 * JOBZ = 'V', LDZ >= max(1,N). 00122 * 00123 * WORK (workspace) COMPLEX*16 array, dimension (2*N) 00124 * 00125 * RWORK (workspace) DOUBLE PRECISION array, dimension (7*N) 00126 * 00127 * IWORK (workspace) INTEGER array, dimension (5*N) 00128 * 00129 * IFAIL (output) INTEGER array, dimension (N) 00130 * If JOBZ = 'V', then if INFO = 0, the first M elements of 00131 * IFAIL are zero. If INFO > 0, then IFAIL contains the 00132 * indices of the eigenvectors that failed to converge. 00133 * If JOBZ = 'N', then IFAIL is not referenced. 00134 * 00135 * INFO (output) INTEGER 00136 * = 0: successful exit 00137 * < 0: if INFO = -i, the i-th argument had an illegal value 00138 * > 0: if INFO = i, then i eigenvectors failed to converge. 00139 * Their indices are stored in array IFAIL. 00140 * 00141 * ===================================================================== 00142 * 00143 * .. Parameters .. 00144 DOUBLE PRECISION ZERO, ONE 00145 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 ) 00146 COMPLEX*16 CONE 00147 PARAMETER ( CONE = ( 1.0D0, 0.0D0 ) ) 00148 * .. 00149 * .. Local Scalars .. 00150 LOGICAL ALLEIG, INDEIG, TEST, VALEIG, WANTZ 00151 CHARACTER ORDER 00152 INTEGER I, IINFO, IMAX, INDD, INDE, INDEE, INDIBL, 00153 $ INDISP, INDIWK, INDRWK, INDTAU, INDWRK, ISCALE, 00154 $ ITMP1, J, JJ, NSPLIT 00155 DOUBLE PRECISION ABSTLL, ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, 00156 $ SIGMA, SMLNUM, TMP1, VLL, VUU 00157 * .. 00158 * .. External Functions .. 00159 LOGICAL LSAME 00160 DOUBLE PRECISION DLAMCH, ZLANHP 00161 EXTERNAL LSAME, DLAMCH, ZLANHP 00162 * .. 00163 * .. External Subroutines .. 00164 EXTERNAL DCOPY, DSCAL, DSTEBZ, DSTERF, XERBLA, ZDSCAL, 00165 $ ZHPTRD, ZSTEIN, ZSTEQR, ZSWAP, ZUPGTR, ZUPMTR 00166 * .. 00167 * .. Intrinsic Functions .. 00168 INTRINSIC DBLE, MAX, MIN, SQRT 00169 * .. 00170 * .. Executable Statements .. 00171 * 00172 * Test the input parameters. 00173 * 00174 WANTZ = LSAME( JOBZ, 'V' ) 00175 ALLEIG = LSAME( RANGE, 'A' ) 00176 VALEIG = LSAME( RANGE, 'V' ) 00177 INDEIG = LSAME( RANGE, 'I' ) 00178 * 00179 INFO = 0 00180 IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN 00181 INFO = -1 00182 ELSE IF( .NOT.( ALLEIG .OR. VALEIG .OR. INDEIG ) ) THEN 00183 INFO = -2 00184 ELSE IF( .NOT.( LSAME( UPLO, 'L' ) .OR. LSAME( UPLO, 'U' ) ) ) 00185 $ THEN 00186 INFO = -3 00187 ELSE IF( N.LT.0 ) THEN 00188 INFO = -4 00189 ELSE 00190 IF( VALEIG ) THEN 00191 IF( N.GT.0 .AND. VU.LE.VL ) 00192 $ INFO = -7 00193 ELSE IF( INDEIG ) THEN 00194 IF( IL.LT.1 .OR. IL.GT.MAX( 1, N ) ) THEN 00195 INFO = -8 00196 ELSE IF( IU.LT.MIN( N, IL ) .OR. IU.GT.N ) THEN 00197 INFO = -9 00198 END IF 00199 END IF 00200 END IF 00201 IF( INFO.EQ.0 ) THEN 00202 IF( LDZ.LT.1 .OR. ( WANTZ .AND. LDZ.LT.N ) ) 00203 $ INFO = -14 00204 END IF 00205 * 00206 IF( INFO.NE.0 ) THEN 00207 CALL XERBLA( 'ZHPEVX', -INFO ) 00208 RETURN 00209 END IF 00210 * 00211 * Quick return if possible 00212 * 00213 M = 0 00214 IF( N.EQ.0 ) 00215 $ RETURN 00216 * 00217 IF( N.EQ.1 ) THEN 00218 IF( ALLEIG .OR. INDEIG ) THEN 00219 M = 1 00220 W( 1 ) = AP( 1 ) 00221 ELSE 00222 IF( VL.LT.DBLE( AP( 1 ) ) .AND. VU.GE.DBLE( AP( 1 ) ) ) THEN 00223 M = 1 00224 W( 1 ) = AP( 1 ) 00225 END IF 00226 END IF 00227 IF( WANTZ ) 00228 $ Z( 1, 1 ) = CONE 00229 RETURN 00230 END IF 00231 * 00232 * Get machine constants. 00233 * 00234 SAFMIN = DLAMCH( 'Safe minimum' ) 00235 EPS = DLAMCH( 'Precision' ) 00236 SMLNUM = SAFMIN / EPS 00237 BIGNUM = ONE / SMLNUM 00238 RMIN = SQRT( SMLNUM ) 00239 RMAX = MIN( SQRT( BIGNUM ), ONE / SQRT( SQRT( SAFMIN ) ) ) 00240 * 00241 * Scale matrix to allowable range, if necessary. 00242 * 00243 ISCALE = 0 00244 ABSTLL = ABSTOL 00245 IF( VALEIG ) THEN 00246 VLL = VL 00247 VUU = VU 00248 ELSE 00249 VLL = ZERO 00250 VUU = ZERO 00251 END IF 00252 ANRM = ZLANHP( 'M', UPLO, N, AP, RWORK ) 00253 IF( ANRM.GT.ZERO .AND. ANRM.LT.RMIN ) THEN 00254 ISCALE = 1 00255 SIGMA = RMIN / ANRM 00256 ELSE IF( ANRM.GT.RMAX ) THEN 00257 ISCALE = 1 00258 SIGMA = RMAX / ANRM 00259 END IF 00260 IF( ISCALE.EQ.1 ) THEN 00261 CALL ZDSCAL( ( N*( N+1 ) ) / 2, SIGMA, AP, 1 ) 00262 IF( ABSTOL.GT.0 ) 00263 $ ABSTLL = ABSTOL*SIGMA 00264 IF( VALEIG ) THEN 00265 VLL = VL*SIGMA 00266 VUU = VU*SIGMA 00267 END IF 00268 END IF 00269 * 00270 * Call ZHPTRD to reduce Hermitian packed matrix to tridiagonal form. 00271 * 00272 INDD = 1 00273 INDE = INDD + N 00274 INDRWK = INDE + N 00275 INDTAU = 1 00276 INDWRK = INDTAU + N 00277 CALL ZHPTRD( UPLO, N, AP, RWORK( INDD ), RWORK( INDE ), 00278 $ WORK( INDTAU ), IINFO ) 00279 * 00280 * If all eigenvalues are desired and ABSTOL is less than or equal 00281 * to zero, then call DSTERF or ZUPGTR and ZSTEQR. If this fails 00282 * for some eigenvalue, then try DSTEBZ. 00283 * 00284 TEST = .FALSE. 00285 IF (INDEIG) THEN 00286 IF (IL.EQ.1 .AND. IU.EQ.N) THEN 00287 TEST = .TRUE. 00288 END IF 00289 END IF 00290 IF ((ALLEIG .OR. TEST) .AND. (ABSTOL.LE.ZERO)) THEN 00291 CALL DCOPY( N, RWORK( INDD ), 1, W, 1 ) 00292 INDEE = INDRWK + 2*N 00293 IF( .NOT.WANTZ ) THEN 00294 CALL DCOPY( N-1, RWORK( INDE ), 1, RWORK( INDEE ), 1 ) 00295 CALL DSTERF( N, W, RWORK( INDEE ), INFO ) 00296 ELSE 00297 CALL ZUPGTR( UPLO, N, AP, WORK( INDTAU ), Z, LDZ, 00298 $ WORK( INDWRK ), IINFO ) 00299 CALL DCOPY( N-1, RWORK( INDE ), 1, RWORK( INDEE ), 1 ) 00300 CALL ZSTEQR( JOBZ, N, W, RWORK( INDEE ), Z, LDZ, 00301 $ RWORK( INDRWK ), INFO ) 00302 IF( INFO.EQ.0 ) THEN 00303 DO 10 I = 1, N 00304 IFAIL( I ) = 0 00305 10 CONTINUE 00306 END IF 00307 END IF 00308 IF( INFO.EQ.0 ) THEN 00309 M = N 00310 GO TO 20 00311 END IF 00312 INFO = 0 00313 END IF 00314 * 00315 * Otherwise, call DSTEBZ and, if eigenvectors are desired, ZSTEIN. 00316 * 00317 IF( WANTZ ) THEN 00318 ORDER = 'B' 00319 ELSE 00320 ORDER = 'E' 00321 END IF 00322 INDIBL = 1 00323 INDISP = INDIBL + N 00324 INDIWK = INDISP + N 00325 CALL DSTEBZ( RANGE, ORDER, N, VLL, VUU, IL, IU, ABSTLL, 00326 $ RWORK( INDD ), RWORK( INDE ), M, NSPLIT, W, 00327 $ IWORK( INDIBL ), IWORK( INDISP ), RWORK( INDRWK ), 00328 $ IWORK( INDIWK ), INFO ) 00329 * 00330 IF( WANTZ ) THEN 00331 CALL ZSTEIN( N, RWORK( INDD ), RWORK( INDE ), M, W, 00332 $ IWORK( INDIBL ), IWORK( INDISP ), Z, LDZ, 00333 $ RWORK( INDRWK ), IWORK( INDIWK ), IFAIL, INFO ) 00334 * 00335 * Apply unitary matrix used in reduction to tridiagonal 00336 * form to eigenvectors returned by ZSTEIN. 00337 * 00338 INDWRK = INDTAU + N 00339 CALL ZUPMTR( 'L', UPLO, 'N', N, M, AP, WORK( INDTAU ), Z, LDZ, 00340 $ WORK( INDWRK ), IINFO ) 00341 END IF 00342 * 00343 * If matrix was scaled, then rescale eigenvalues appropriately. 00344 * 00345 20 CONTINUE 00346 IF( ISCALE.EQ.1 ) THEN 00347 IF( INFO.EQ.0 ) THEN 00348 IMAX = M 00349 ELSE 00350 IMAX = INFO - 1 00351 END IF 00352 CALL DSCAL( IMAX, ONE / SIGMA, W, 1 ) 00353 END IF 00354 * 00355 * If eigenvalues are not in order, then sort them, along with 00356 * eigenvectors. 00357 * 00358 IF( WANTZ ) THEN 00359 DO 40 J = 1, M - 1 00360 I = 0 00361 TMP1 = W( J ) 00362 DO 30 JJ = J + 1, M 00363 IF( W( JJ ).LT.TMP1 ) THEN 00364 I = JJ 00365 TMP1 = W( JJ ) 00366 END IF 00367 30 CONTINUE 00368 * 00369 IF( I.NE.0 ) THEN 00370 ITMP1 = IWORK( INDIBL+I-1 ) 00371 W( I ) = W( J ) 00372 IWORK( INDIBL+I-1 ) = IWORK( INDIBL+J-1 ) 00373 W( J ) = TMP1 00374 IWORK( INDIBL+J-1 ) = ITMP1 00375 CALL ZSWAP( N, Z( 1, I ), 1, Z( 1, J ), 1 ) 00376 IF( INFO.NE.0 ) THEN 00377 ITMP1 = IFAIL( I ) 00378 IFAIL( I ) = IFAIL( J ) 00379 IFAIL( J ) = ITMP1 00380 END IF 00381 END IF 00382 40 CONTINUE 00383 END IF 00384 * 00385 RETURN 00386 * 00387 * End of ZHPEVX 00388 * 00389 END
1.7.3 
