|
LAPACK 3.3.0
|
|
LAPACK 3.3.0
|
00001 SUBROUTINE CHEEVX( JOBZ, RANGE, UPLO, N, A, LDA, VL, VU, IL, IU, 00002 $ ABSTOL, M, W, Z, LDZ, WORK, LWORK, RWORK, 00003 $ IWORK, 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, LDA, LDZ, LWORK, M, N 00013 REAL ABSTOL, VL, VU 00014 * .. 00015 * .. Array Arguments .. 00016 INTEGER IFAIL( * ), IWORK( * ) 00017 REAL RWORK( * ), W( * ) 00018 COMPLEX A( LDA, * ), WORK( * ), Z( LDZ, * ) 00019 * .. 00020 * 00021 * Purpose 00022 * ======= 00023 * 00024 * CHEEVX computes selected eigenvalues and, optionally, eigenvectors 00025 * of a complex Hermitian matrix A. Eigenvalues and eigenvectors can 00026 * be selected by specifying either a range of values or a range of 00027 * 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 * A (input/output) COMPLEX array, dimension (LDA, N) 00050 * On entry, the Hermitian matrix A. If UPLO = 'U', the 00051 * leading N-by-N upper triangular part of A contains the 00052 * upper triangular part of the matrix A. If UPLO = 'L', 00053 * the leading N-by-N lower triangular part of A contains 00054 * the lower triangular part of the matrix A. 00055 * On exit, the lower triangle (if UPLO='L') or the upper 00056 * triangle (if UPLO='U') of A, including the diagonal, is 00057 * destroyed. 00058 * 00059 * LDA (input) INTEGER 00060 * The leading dimension of the array A. LDA >= max(1,N). 00061 * 00062 * VL (input) REAL 00063 * VU (input) REAL 00064 * If RANGE='V', the lower and upper bounds of the interval to 00065 * be searched for eigenvalues. VL < VU. 00066 * Not referenced if RANGE = 'A' or 'I'. 00067 * 00068 * IL (input) INTEGER 00069 * IU (input) INTEGER 00070 * If RANGE='I', the indices (in ascending order) of the 00071 * smallest and largest eigenvalues to be returned. 00072 * 1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0. 00073 * Not referenced if RANGE = 'A' or 'V'. 00074 * 00075 * ABSTOL (input) REAL 00076 * The absolute error tolerance for the eigenvalues. 00077 * An approximate eigenvalue is accepted as converged 00078 * when it is determined to lie in an interval [a,b] 00079 * of width less than or equal to 00080 * 00081 * ABSTOL + EPS * max( |a|,|b| ) , 00082 * 00083 * where EPS is the machine precision. If ABSTOL is less than 00084 * or equal to zero, then EPS*|T| will be used in its place, 00085 * where |T| is the 1-norm of the tridiagonal matrix obtained 00086 * by reducing A to tridiagonal form. 00087 * 00088 * Eigenvalues will be computed most accurately when ABSTOL is 00089 * set to twice the underflow threshold 2*SLAMCH('S'), not zero. 00090 * If this routine returns with INFO>0, indicating that some 00091 * eigenvectors did not converge, try setting ABSTOL to 00092 * 2*SLAMCH('S'). 00093 * 00094 * See "Computing Small Singular Values of Bidiagonal Matrices 00095 * with Guaranteed High Relative Accuracy," by Demmel and 00096 * Kahan, LAPACK Working Note #3. 00097 * 00098 * M (output) INTEGER 00099 * The total number of eigenvalues found. 0 <= M <= N. 00100 * If RANGE = 'A', M = N, and if RANGE = 'I', M = IU-IL+1. 00101 * 00102 * W (output) REAL array, dimension (N) 00103 * On normal exit, the first M elements contain the selected 00104 * eigenvalues in ascending order. 00105 * 00106 * Z (output) COMPLEX 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 the 00113 * 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/output) COMPLEX array, dimension (MAX(1,LWORK)) 00124 * On exit, if INFO = 0, WORK(1) returns the optimal LWORK. 00125 * 00126 * LWORK (input) INTEGER 00127 * The length of the array WORK. LWORK >= 1, when N <= 1; 00128 * otherwise 2*N. 00129 * For optimal efficiency, LWORK >= (NB+1)*N, 00130 * where NB is the max of the blocksize for CHETRD and for 00131 * CUNMTR as returned by ILAENV. 00132 * 00133 * If LWORK = -1, then a workspace query is assumed; the routine 00134 * only calculates the optimal size of the WORK array, returns 00135 * this value as the first entry of the WORK array, and no error 00136 * message related to LWORK is issued by XERBLA. 00137 * 00138 * RWORK (workspace) REAL array, dimension (7*N) 00139 * 00140 * IWORK (workspace) INTEGER array, dimension (5*N) 00141 * 00142 * IFAIL (output) INTEGER array, dimension (N) 00143 * If JOBZ = 'V', then if INFO = 0, the first M elements of 00144 * IFAIL are zero. If INFO > 0, then IFAIL contains the 00145 * indices of the eigenvectors that failed to converge. 00146 * If JOBZ = 'N', then IFAIL is not referenced. 00147 * 00148 * INFO (output) INTEGER 00149 * = 0: successful exit 00150 * < 0: if INFO = -i, the i-th argument had an illegal value 00151 * > 0: if INFO = i, then i eigenvectors failed to converge. 00152 * Their indices are stored in array IFAIL. 00153 * 00154 * ===================================================================== 00155 * 00156 * .. Parameters .. 00157 REAL ZERO, ONE 00158 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) 00159 COMPLEX CONE 00160 PARAMETER ( CONE = ( 1.0E+0, 0.0E+0 ) ) 00161 * .. 00162 * .. Local Scalars .. 00163 LOGICAL ALLEIG, INDEIG, LOWER, LQUERY, TEST, VALEIG, 00164 $ WANTZ 00165 CHARACTER ORDER 00166 INTEGER I, IINFO, IMAX, INDD, INDE, INDEE, INDIBL, 00167 $ INDISP, INDIWK, INDRWK, INDTAU, INDWRK, ISCALE, 00168 $ ITMP1, J, JJ, LLWORK, LWKMIN, LWKOPT, NB, 00169 $ NSPLIT 00170 REAL ABSTLL, ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, 00171 $ SIGMA, SMLNUM, TMP1, VLL, VUU 00172 * .. 00173 * .. External Functions .. 00174 LOGICAL LSAME 00175 INTEGER ILAENV 00176 REAL CLANHE, SLAMCH 00177 EXTERNAL LSAME, ILAENV, CLANHE, SLAMCH 00178 * .. 00179 * .. External Subroutines .. 00180 EXTERNAL CHETRD, CLACPY, CSSCAL, CSTEIN, CSTEQR, CSWAP, 00181 $ CUNGTR, CUNMTR, SCOPY, SSCAL, SSTEBZ, SSTERF, 00182 $ XERBLA 00183 * .. 00184 * .. Intrinsic Functions .. 00185 INTRINSIC MAX, MIN, REAL, SQRT 00186 * .. 00187 * .. Executable Statements .. 00188 * 00189 * Test the input parameters. 00190 * 00191 LOWER = LSAME( UPLO, 'L' ) 00192 WANTZ = LSAME( JOBZ, 'V' ) 00193 ALLEIG = LSAME( RANGE, 'A' ) 00194 VALEIG = LSAME( RANGE, 'V' ) 00195 INDEIG = LSAME( RANGE, 'I' ) 00196 LQUERY = ( LWORK.EQ.-1 ) 00197 * 00198 INFO = 0 00199 IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN 00200 INFO = -1 00201 ELSE IF( .NOT.( ALLEIG .OR. VALEIG .OR. INDEIG ) ) THEN 00202 INFO = -2 00203 ELSE IF( .NOT.( LOWER .OR. LSAME( UPLO, 'U' ) ) ) THEN 00204 INFO = -3 00205 ELSE IF( N.LT.0 ) THEN 00206 INFO = -4 00207 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 00208 INFO = -6 00209 ELSE 00210 IF( VALEIG ) THEN 00211 IF( N.GT.0 .AND. VU.LE.VL ) 00212 $ INFO = -8 00213 ELSE IF( INDEIG ) THEN 00214 IF( IL.LT.1 .OR. IL.GT.MAX( 1, N ) ) THEN 00215 INFO = -9 00216 ELSE IF( IU.LT.MIN( N, IL ) .OR. IU.GT.N ) THEN 00217 INFO = -10 00218 END IF 00219 END IF 00220 END IF 00221 IF( INFO.EQ.0 ) THEN 00222 IF( LDZ.LT.1 .OR. ( WANTZ .AND. LDZ.LT.N ) ) THEN 00223 INFO = -15 00224 END IF 00225 END IF 00226 * 00227 IF( INFO.EQ.0 ) THEN 00228 IF( N.LE.1 ) THEN 00229 LWKMIN = 1 00230 WORK( 1 ) = LWKMIN 00231 ELSE 00232 LWKMIN = 2*N 00233 NB = ILAENV( 1, 'CHETRD', UPLO, N, -1, -1, -1 ) 00234 NB = MAX( NB, ILAENV( 1, 'CUNMTR', UPLO, N, -1, -1, -1 ) ) 00235 LWKOPT = MAX( 1, ( NB + 1 )*N ) 00236 WORK( 1 ) = LWKOPT 00237 END IF 00238 * 00239 IF( LWORK.LT.LWKMIN .AND. .NOT.LQUERY ) 00240 $ INFO = -17 00241 END IF 00242 * 00243 IF( INFO.NE.0 ) THEN 00244 CALL XERBLA( 'CHEEVX', -INFO ) 00245 RETURN 00246 ELSE IF( LQUERY ) THEN 00247 RETURN 00248 END IF 00249 * 00250 * Quick return if possible 00251 * 00252 M = 0 00253 IF( N.EQ.0 ) THEN 00254 RETURN 00255 END IF 00256 * 00257 IF( N.EQ.1 ) THEN 00258 IF( ALLEIG .OR. INDEIG ) THEN 00259 M = 1 00260 W( 1 ) = A( 1, 1 ) 00261 ELSE IF( VALEIG ) THEN 00262 IF( VL.LT.REAL( A( 1, 1 ) ) .AND. VU.GE.REAL( A( 1, 1 ) ) ) 00263 $ THEN 00264 M = 1 00265 W( 1 ) = A( 1, 1 ) 00266 END IF 00267 END IF 00268 IF( WANTZ ) 00269 $ Z( 1, 1 ) = CONE 00270 RETURN 00271 END IF 00272 * 00273 * Get machine constants. 00274 * 00275 SAFMIN = SLAMCH( 'Safe minimum' ) 00276 EPS = SLAMCH( 'Precision' ) 00277 SMLNUM = SAFMIN / EPS 00278 BIGNUM = ONE / SMLNUM 00279 RMIN = SQRT( SMLNUM ) 00280 RMAX = MIN( SQRT( BIGNUM ), ONE / SQRT( SQRT( SAFMIN ) ) ) 00281 * 00282 * Scale matrix to allowable range, if necessary. 00283 * 00284 ISCALE = 0 00285 ABSTLL = ABSTOL 00286 IF( VALEIG ) THEN 00287 VLL = VL 00288 VUU = VU 00289 END IF 00290 ANRM = CLANHE( 'M', UPLO, N, A, LDA, RWORK ) 00291 IF( ANRM.GT.ZERO .AND. ANRM.LT.RMIN ) THEN 00292 ISCALE = 1 00293 SIGMA = RMIN / ANRM 00294 ELSE IF( ANRM.GT.RMAX ) THEN 00295 ISCALE = 1 00296 SIGMA = RMAX / ANRM 00297 END IF 00298 IF( ISCALE.EQ.1 ) THEN 00299 IF( LOWER ) THEN 00300 DO 10 J = 1, N 00301 CALL CSSCAL( N-J+1, SIGMA, A( J, J ), 1 ) 00302 10 CONTINUE 00303 ELSE 00304 DO 20 J = 1, N 00305 CALL CSSCAL( J, SIGMA, A( 1, J ), 1 ) 00306 20 CONTINUE 00307 END IF 00308 IF( ABSTOL.GT.0 ) 00309 $ ABSTLL = ABSTOL*SIGMA 00310 IF( VALEIG ) THEN 00311 VLL = VL*SIGMA 00312 VUU = VU*SIGMA 00313 END IF 00314 END IF 00315 * 00316 * Call CHETRD to reduce Hermitian matrix to tridiagonal form. 00317 * 00318 INDD = 1 00319 INDE = INDD + N 00320 INDRWK = INDE + N 00321 INDTAU = 1 00322 INDWRK = INDTAU + N 00323 LLWORK = LWORK - INDWRK + 1 00324 CALL CHETRD( UPLO, N, A, LDA, RWORK( INDD ), RWORK( INDE ), 00325 $ WORK( INDTAU ), WORK( INDWRK ), LLWORK, IINFO ) 00326 * 00327 * If all eigenvalues are desired and ABSTOL is less than or equal to 00328 * zero, then call SSTERF or CUNGTR and CSTEQR. If this fails for 00329 * some eigenvalue, then try SSTEBZ. 00330 * 00331 TEST = .FALSE. 00332 IF( INDEIG ) THEN 00333 IF( IL.EQ.1 .AND. IU.EQ.N ) THEN 00334 TEST = .TRUE. 00335 END IF 00336 END IF 00337 IF( ( ALLEIG .OR. TEST ) .AND. ( ABSTOL.LE.ZERO ) ) THEN 00338 CALL SCOPY( N, RWORK( INDD ), 1, W, 1 ) 00339 INDEE = INDRWK + 2*N 00340 IF( .NOT.WANTZ ) THEN 00341 CALL SCOPY( N-1, RWORK( INDE ), 1, RWORK( INDEE ), 1 ) 00342 CALL SSTERF( N, W, RWORK( INDEE ), INFO ) 00343 ELSE 00344 CALL CLACPY( 'A', N, N, A, LDA, Z, LDZ ) 00345 CALL CUNGTR( UPLO, N, Z, LDZ, WORK( INDTAU ), 00346 $ WORK( INDWRK ), LLWORK, IINFO ) 00347 CALL SCOPY( N-1, RWORK( INDE ), 1, RWORK( INDEE ), 1 ) 00348 CALL CSTEQR( JOBZ, N, W, RWORK( INDEE ), Z, LDZ, 00349 $ RWORK( INDRWK ), INFO ) 00350 IF( INFO.EQ.0 ) THEN 00351 DO 30 I = 1, N 00352 IFAIL( I ) = 0 00353 30 CONTINUE 00354 END IF 00355 END IF 00356 IF( INFO.EQ.0 ) THEN 00357 M = N 00358 GO TO 40 00359 END IF 00360 INFO = 0 00361 END IF 00362 * 00363 * Otherwise, call SSTEBZ and, if eigenvectors are desired, CSTEIN. 00364 * 00365 IF( WANTZ ) THEN 00366 ORDER = 'B' 00367 ELSE 00368 ORDER = 'E' 00369 END IF 00370 INDIBL = 1 00371 INDISP = INDIBL + N 00372 INDIWK = INDISP + N 00373 CALL SSTEBZ( RANGE, ORDER, N, VLL, VUU, IL, IU, ABSTLL, 00374 $ RWORK( INDD ), RWORK( INDE ), M, NSPLIT, W, 00375 $ IWORK( INDIBL ), IWORK( INDISP ), RWORK( INDRWK ), 00376 $ IWORK( INDIWK ), INFO ) 00377 * 00378 IF( WANTZ ) THEN 00379 CALL CSTEIN( N, RWORK( INDD ), RWORK( INDE ), M, W, 00380 $ IWORK( INDIBL ), IWORK( INDISP ), Z, LDZ, 00381 $ RWORK( INDRWK ), IWORK( INDIWK ), IFAIL, INFO ) 00382 * 00383 * Apply unitary matrix used in reduction to tridiagonal 00384 * form to eigenvectors returned by CSTEIN. 00385 * 00386 CALL CUNMTR( 'L', UPLO, 'N', N, M, A, LDA, WORK( INDTAU ), Z, 00387 $ LDZ, WORK( INDWRK ), LLWORK, IINFO ) 00388 END IF 00389 * 00390 * If matrix was scaled, then rescale eigenvalues appropriately. 00391 * 00392 40 CONTINUE 00393 IF( ISCALE.EQ.1 ) THEN 00394 IF( INFO.EQ.0 ) THEN 00395 IMAX = M 00396 ELSE 00397 IMAX = INFO - 1 00398 END IF 00399 CALL SSCAL( IMAX, ONE / SIGMA, W, 1 ) 00400 END IF 00401 * 00402 * If eigenvalues are not in order, then sort them, along with 00403 * eigenvectors. 00404 * 00405 IF( WANTZ ) THEN 00406 DO 60 J = 1, M - 1 00407 I = 0 00408 TMP1 = W( J ) 00409 DO 50 JJ = J + 1, M 00410 IF( W( JJ ).LT.TMP1 ) THEN 00411 I = JJ 00412 TMP1 = W( JJ ) 00413 END IF 00414 50 CONTINUE 00415 * 00416 IF( I.NE.0 ) THEN 00417 ITMP1 = IWORK( INDIBL+I-1 ) 00418 W( I ) = W( J ) 00419 IWORK( INDIBL+I-1 ) = IWORK( INDIBL+J-1 ) 00420 W( J ) = TMP1 00421 IWORK( INDIBL+J-1 ) = ITMP1 00422 CALL CSWAP( N, Z( 1, I ), 1, Z( 1, J ), 1 ) 00423 IF( INFO.NE.0 ) THEN 00424 ITMP1 = IFAIL( I ) 00425 IFAIL( I ) = IFAIL( J ) 00426 IFAIL( J ) = ITMP1 00427 END IF 00428 END IF 00429 60 CONTINUE 00430 END IF 00431 * 00432 * Set WORK(1) to optimal complex workspace size. 00433 * 00434 WORK( 1 ) = LWKOPT 00435 * 00436 RETURN 00437 * 00438 * End of CHEEVX 00439 * 00440 END
1.7.3 
