 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
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.3.1) -- 00006 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00007 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00008 * -- April 2011 -- 00009 * @generated c 00010 * 00011 * .. Scalar Arguments .. 00012 CHARACTER JOBZ, RANGE, UPLO 00013 INTEGER IL, INFO, IU, LDA, LDZ, LWORK, M, N 00014 REAL ABSTOL, VL, VU 00015 * .. 00016 * .. Array Arguments .. 00017 INTEGER IFAIL( * ), IWORK( * ) 00018 REAL RWORK( * ), W( * ) 00019 COMPLEX A( LDA, * ), WORK( * ), Z( LDZ, * ) 00020 * .. 00021 * 00022 * Purpose 00023 * ======= 00024 * 00025 * CHEEVX computes selected eigenvalues and, optionally, eigenvectors 00026 * of a complex Hermitian matrix A. Eigenvalues and eigenvectors can 00027 * be selected by specifying either a range of values or a range of 00028 * indices for the desired eigenvalues. 00029 * 00030 * Arguments 00031 * ========= 00032 * 00033 * JOBZ (input) CHARACTER*1 00034 * = 'N': Compute eigenvalues only; 00035 * = 'V': Compute eigenvalues and eigenvectors. 00036 * 00037 * RANGE (input) CHARACTER*1 00038 * = 'A': all eigenvalues will be found. 00039 * = 'V': all eigenvalues in the half-open interval (VL,VU] 00040 * will be found. 00041 * = 'I': the IL-th through IU-th eigenvalues will be found. 00042 * 00043 * UPLO (input) CHARACTER*1 00044 * = 'U': Upper triangle of A is stored; 00045 * = 'L': Lower triangle of A is stored. 00046 * 00047 * N (input) INTEGER 00048 * The order of the matrix A. N >= 0. 00049 * 00050 * A (input/output) COMPLEX array, dimension (LDA, N) 00051 * On entry, the Hermitian matrix A. If UPLO = 'U', the 00052 * leading N-by-N upper triangular part of A contains the 00053 * upper triangular part of the matrix A. If UPLO = 'L', 00054 * the leading N-by-N lower triangular part of A contains 00055 * the lower triangular part of the matrix A. 00056 * On exit, the lower triangle (if UPLO='L') or the upper 00057 * triangle (if UPLO='U') of A, including the diagonal, is 00058 * destroyed. 00059 * 00060 * LDA (input) INTEGER 00061 * The leading dimension of the array A. LDA >= max(1,N). 00062 * 00063 * VL (input) REAL 00064 * VU (input) REAL 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) REAL 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 A to tridiagonal form. 00088 * 00089 * Eigenvalues will be computed most accurately when ABSTOL is 00090 * set to twice the underflow threshold 2*SLAMCH('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*SLAMCH('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) REAL array, dimension (N) 00104 * On normal exit, the first M elements contain the selected 00105 * eigenvalues in ascending order. 00106 * 00107 * Z (output) COMPLEX array, dimension (LDZ, max(1,M)) 00108 * If JOBZ = 'V', then if INFO = 0, the first M columns of Z 00109 * contain the orthonormal eigenvectors of the matrix A 00110 * corresponding to the selected eigenvalues, with the i-th 00111 * column of Z holding the eigenvector associated with W(i). 00112 * If an eigenvector fails to converge, then that column of Z 00113 * contains the latest approximation to the eigenvector, and the 00114 * index of the eigenvector is returned in IFAIL. 00115 * If JOBZ = 'N', then Z is not referenced. 00116 * Note: the user must ensure that at least max(1,M) columns are 00117 * supplied in the array Z; if RANGE = 'V', the exact value of M 00118 * is not known in advance and an upper bound must be used. 00119 * 00120 * LDZ (input) INTEGER 00121 * The leading dimension of the array Z. LDZ >= 1, and if 00122 * JOBZ = 'V', LDZ >= max(1,N). 00123 * 00124 * WORK (workspace/output) COMPLEX array, dimension (MAX(1,LWORK)) 00125 * On exit, if INFO = 0, WORK(1) returns the optimal LWORK. 00126 * 00127 * LWORK (input) INTEGER 00128 * The length of the array WORK. LWORK >= 1, when N <= 1; 00129 * otherwise 2*N. 00130 * For optimal efficiency, LWORK >= (NB+1)*N, 00131 * where NB is the max of the blocksize for CHETRD and for 00132 * CUNMTR as returned by ILAENV. 00133 * 00134 * If LWORK = -1, then a workspace query is assumed; the routine 00135 * only calculates the optimal size of the WORK array, returns 00136 * this value as the first entry of the WORK array, and no error 00137 * message related to LWORK is issued by XERBLA. 00138 * 00139 * RWORK (workspace) REAL array, dimension (7*N) 00140 * 00141 * IWORK (workspace) INTEGER array, dimension (5*N) 00142 * 00143 * IFAIL (output) INTEGER array, dimension (N) 00144 * If JOBZ = 'V', then if INFO = 0, the first M elements of 00145 * IFAIL are zero. If INFO > 0, then IFAIL contains the 00146 * indices of the eigenvectors that failed to converge. 00147 * If JOBZ = 'N', then IFAIL is not referenced. 00148 * 00149 * INFO (output) INTEGER 00150 * = 0: successful exit 00151 * < 0: if INFO = -i, the i-th argument had an illegal value 00152 * > 0: if INFO = i, then i eigenvectors failed to converge. 00153 * Their indices are stored in array IFAIL. 00154 * 00155 * ===================================================================== 00156 * 00157 * .. Parameters .. 00158 REAL ZERO, ONE 00159 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) 00160 COMPLEX CONE 00161 PARAMETER ( CONE = ( 1.0E+0, 0.0E+0 ) ) 00162 * .. 00163 * .. Local Scalars .. 00164 LOGICAL ALLEIG, INDEIG, LOWER, LQUERY, TEST, VALEIG, 00165 $ WANTZ 00166 CHARACTER ORDER 00167 INTEGER I, IINFO, IMAX, INDD, INDE, INDEE, INDIBL, 00168 $ INDISP, INDIWK, INDRWK, INDTAU, INDWRK, ISCALE, 00169 $ ITMP1, J, JJ, LLWORK, LWKMIN, LWKOPT, NB, 00170 $ NSPLIT 00171 REAL ABSTLL, ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, 00172 $ SIGMA, SMLNUM, TMP1, VLL, VUU 00173 * .. 00174 * .. External Functions .. 00175 LOGICAL LSAME 00176 INTEGER ILAENV 00177 REAL SLAMCH, CLANHE 00178 EXTERNAL LSAME, ILAENV, SLAMCH, CLANHE 00179 * .. 00180 * .. External Subroutines .. 00181 EXTERNAL SCOPY, SSCAL, SSTEBZ, SSTERF, XERBLA, CSSCAL, 00182 $ CHETRD, CLACPY, CSTEIN, CSTEQR, CSWAP, CUNGTR, 00183 $ CUNMTR 00184 * .. 00185 * .. Intrinsic Functions .. 00186 INTRINSIC REAL, MAX, MIN, SQRT 00187 * .. 00188 * .. Executable Statements .. 00189 * 00190 * Test the input parameters. 00191 * 00192 LOWER = LSAME( UPLO, 'L' ) 00193 WANTZ = LSAME( JOBZ, 'V' ) 00194 ALLEIG = LSAME( RANGE, 'A' ) 00195 VALEIG = LSAME( RANGE, 'V' ) 00196 INDEIG = LSAME( RANGE, 'I' ) 00197 LQUERY = ( LWORK.EQ.-1 ) 00198 * 00199 INFO = 0 00200 IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN 00201 INFO = -1 00202 ELSE IF( .NOT.( ALLEIG .OR. VALEIG .OR. INDEIG ) ) THEN 00203 INFO = -2 00204 ELSE IF( .NOT.( LOWER .OR. LSAME( UPLO, 'U' ) ) ) THEN 00205 INFO = -3 00206 ELSE IF( N.LT.0 ) THEN 00207 INFO = -4 00208 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 00209 INFO = -6 00210 ELSE 00211 IF( VALEIG ) THEN 00212 IF( N.GT.0 .AND. VU.LE.VL ) 00213 $ INFO = -8 00214 ELSE IF( INDEIG ) THEN 00215 IF( IL.LT.1 .OR. IL.GT.MAX( 1, N ) ) THEN 00216 INFO = -9 00217 ELSE IF( IU.LT.MIN( N, IL ) .OR. IU.GT.N ) THEN 00218 INFO = -10 00219 END IF 00220 END IF 00221 END IF 00222 IF( INFO.EQ.0 ) THEN 00223 IF( LDZ.LT.1 .OR. ( WANTZ .AND. LDZ.LT.N ) ) THEN 00224 INFO = -15 00225 END IF 00226 END IF 00227 * 00228 IF( INFO.EQ.0 ) THEN 00229 IF( N.LE.1 ) THEN 00230 LWKMIN = 1 00231 WORK( 1 ) = LWKMIN 00232 ELSE 00233 LWKMIN = 2*N 00234 NB = ILAENV( 1, 'CHETRD', UPLO, N, -1, -1, -1 ) 00235 NB = MAX( NB, ILAENV( 1, 'CUNMTR', UPLO, N, -1, -1, -1 ) ) 00236 LWKOPT = MAX( 1, ( NB + 1 )*N ) 00237 WORK( 1 ) = LWKOPT 00238 END IF 00239 * 00240 IF( LWORK.LT.LWKMIN .AND. .NOT.LQUERY ) 00241 $ INFO = -17 00242 END IF 00243 * 00244 IF( INFO.NE.0 ) THEN 00245 CALL XERBLA( 'CHEEVX', -INFO ) 00246 RETURN 00247 ELSE IF( LQUERY ) THEN 00248 RETURN 00249 END IF 00250 * 00251 * Quick return if possible 00252 * 00253 M = 0 00254 IF( N.EQ.0 ) THEN 00255 RETURN 00256 END IF 00257 * 00258 IF( N.EQ.1 ) THEN 00259 IF( ALLEIG .OR. INDEIG ) THEN 00260 M = 1 00261 W( 1 ) = A( 1, 1 ) 00262 ELSE IF( VALEIG ) THEN 00263 IF( VL.LT.REAL( A( 1, 1 ) ) .AND. VU.GE.REAL( A( 1, 1 ) ) ) 00264 $ THEN 00265 M = 1 00266 W( 1 ) = A( 1, 1 ) 00267 END IF 00268 END IF 00269 IF( WANTZ ) 00270 $ Z( 1, 1 ) = CONE 00271 RETURN 00272 END IF 00273 * 00274 * Get machine constants. 00275 * 00276 SAFMIN = SLAMCH( 'Safe minimum' ) 00277 EPS = SLAMCH( 'Precision' ) 00278 SMLNUM = SAFMIN / EPS 00279 BIGNUM = ONE / SMLNUM 00280 RMIN = SQRT( SMLNUM ) 00281 RMAX = MIN( SQRT( BIGNUM ), ONE / SQRT( SQRT( SAFMIN ) ) ) 00282 * 00283 * Scale matrix to allowable range, if necessary. 00284 * 00285 ISCALE = 0 00286 ABSTLL = ABSTOL 00287 IF( VALEIG ) THEN 00288 VLL = VL 00289 VUU = VU 00290 END IF 00291 ANRM = CLANHE( 'M', UPLO, N, A, LDA, RWORK ) 00292 IF( ANRM.GT.ZERO .AND. ANRM.LT.RMIN ) THEN 00293 ISCALE = 1 00294 SIGMA = RMIN / ANRM 00295 ELSE IF( ANRM.GT.RMAX ) THEN 00296 ISCALE = 1 00297 SIGMA = RMAX / ANRM 00298 END IF 00299 IF( ISCALE.EQ.1 ) THEN 00300 IF( LOWER ) THEN 00301 DO 10 J = 1, N 00302 CALL CSSCAL( N-J+1, SIGMA, A( J, J ), 1 ) 00303 10 CONTINUE 00304 ELSE 00305 DO 20 J = 1, N 00306 CALL CSSCAL( J, SIGMA, A( 1, J ), 1 ) 00307 20 CONTINUE 00308 END IF 00309 IF( ABSTOL.GT.0 ) 00310 $ ABSTLL = ABSTOL*SIGMA 00311 IF( VALEIG ) THEN 00312 VLL = VL*SIGMA 00313 VUU = VU*SIGMA 00314 END IF 00315 END IF 00316 * 00317 * Call CHETRD to reduce Hermitian matrix to tridiagonal form. 00318 * 00319 INDD = 1 00320 INDE = INDD + N 00321 INDRWK = INDE + N 00322 INDTAU = 1 00323 INDWRK = INDTAU + N 00324 LLWORK = LWORK - INDWRK + 1 00325 CALL CHETRD( UPLO, N, A, LDA, RWORK( INDD ), RWORK( INDE ), 00326 $ WORK( INDTAU ), WORK( INDWRK ), LLWORK, IINFO ) 00327 * 00328 * If all eigenvalues are desired and ABSTOL is less than or equal to 00329 * zero, then call SSTERF or CUNGTR and CSTEQR. If this fails for 00330 * some eigenvalue, then try SSTEBZ. 00331 * 00332 TEST = .FALSE. 00333 IF( INDEIG ) THEN 00334 IF( IL.EQ.1 .AND. IU.EQ.N ) THEN 00335 TEST = .TRUE. 00336 END IF 00337 END IF 00338 IF( ( ALLEIG .OR. TEST ) .AND. ( ABSTOL.LE.ZERO ) ) THEN 00339 CALL SCOPY( N, RWORK( INDD ), 1, W, 1 ) 00340 INDEE = INDRWK + 2*N 00341 IF( .NOT.WANTZ ) THEN 00342 CALL SCOPY( N-1, RWORK( INDE ), 1, RWORK( INDEE ), 1 ) 00343 CALL SSTERF( N, W, RWORK( INDEE ), INFO ) 00344 ELSE 00345 CALL CLACPY( 'A', N, N, A, LDA, Z, LDZ ) 00346 CALL CUNGTR( UPLO, N, Z, LDZ, WORK( INDTAU ), 00347 $ WORK( INDWRK ), LLWORK, IINFO ) 00348 CALL SCOPY( N-1, RWORK( INDE ), 1, RWORK( INDEE ), 1 ) 00349 CALL CSTEQR( JOBZ, N, W, RWORK( INDEE ), Z, LDZ, 00350 $ RWORK( INDRWK ), INFO ) 00351 IF( INFO.EQ.0 ) THEN 00352 DO 30 I = 1, N 00353 IFAIL( I ) = 0 00354 30 CONTINUE 00355 END IF 00356 END IF 00357 IF( INFO.EQ.0 ) THEN 00358 M = N 00359 GO TO 40 00360 END IF 00361 INFO = 0 00362 END IF 00363 * 00364 * Otherwise, call SSTEBZ and, if eigenvectors are desired, CSTEIN. 00365 * 00366 IF( WANTZ ) THEN 00367 ORDER = 'B' 00368 ELSE 00369 ORDER = 'E' 00370 END IF 00371 INDIBL = 1 00372 INDISP = INDIBL + N 00373 INDIWK = INDISP + N 00374 CALL SSTEBZ( RANGE, ORDER, N, VLL, VUU, IL, IU, ABSTLL, 00375 $ RWORK( INDD ), RWORK( INDE ), M, NSPLIT, W, 00376 $ IWORK( INDIBL ), IWORK( INDISP ), RWORK( INDRWK ), 00377 $ IWORK( INDIWK ), INFO ) 00378 * 00379 IF( WANTZ ) THEN 00380 CALL CSTEIN( N, RWORK( INDD ), RWORK( INDE ), M, W, 00381 $ IWORK( INDIBL ), IWORK( INDISP ), Z, LDZ, 00382 $ RWORK( INDRWK ), IWORK( INDIWK ), IFAIL, INFO ) 00383 * 00384 * Apply unitary matrix used in reduction to tridiagonal 00385 * form to eigenvectors returned by CSTEIN. 00386 * 00387 CALL CUNMTR( 'L', UPLO, 'N', N, M, A, LDA, WORK( INDTAU ), Z, 00388 $ LDZ, WORK( INDWRK ), LLWORK, IINFO ) 00389 END IF 00390 * 00391 * If matrix was scaled, then rescale eigenvalues appropriately. 00392 * 00393 40 CONTINUE 00394 IF( ISCALE.EQ.1 ) THEN 00395 IF( INFO.EQ.0 ) THEN 00396 IMAX = M 00397 ELSE 00398 IMAX = INFO - 1 00399 END IF 00400 CALL SSCAL( IMAX, ONE / SIGMA, W, 1 ) 00401 END IF 00402 * 00403 * If eigenvalues are not in order, then sort them, along with 00404 * eigenvectors. 00405 * 00406 IF( WANTZ ) THEN 00407 DO 60 J = 1, M - 1 00408 I = 0 00409 TMP1 = W( J ) 00410 DO 50 JJ = J + 1, M 00411 IF( W( JJ ).LT.TMP1 ) THEN 00412 I = JJ 00413 TMP1 = W( JJ ) 00414 END IF 00415 50 CONTINUE 00416 * 00417 IF( I.NE.0 ) THEN 00418 ITMP1 = IWORK( INDIBL+I-1 ) 00419 W( I ) = W( J ) 00420 IWORK( INDIBL+I-1 ) = IWORK( INDIBL+J-1 ) 00421 W( J ) = TMP1 00422 IWORK( INDIBL+J-1 ) = ITMP1 00423 CALL CSWAP( N, Z( 1, I ), 1, Z( 1, J ), 1 ) 00424 IF( INFO.NE.0 ) THEN 00425 ITMP1 = IFAIL( I ) 00426 IFAIL( I ) = IFAIL( J ) 00427 IFAIL( J ) = ITMP1 00428 END IF 00429 END IF 00430 60 CONTINUE 00431 END IF 00432 * 00433 * Set WORK(1) to optimal complex workspace size. 00434 * 00435 WORK( 1 ) = LWKOPT 00436 * 00437 RETURN 00438 * 00439 * End of CHEEVX 00440 * 00441 END
1.7.3 
