![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief <b> CHBEVX computes the eigenvalues and, optionally, the left and/or right eigenvectors for OTHER matrices</b> 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download CHBEVX + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/chbevx.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/chbevx.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/chbevx.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE CHBEVX( JOBZ, RANGE, UPLO, N, KD, AB, LDAB, Q, LDQ, VL, 00022 * VU, IL, IU, ABSTOL, M, W, Z, LDZ, WORK, RWORK, 00023 * IWORK, IFAIL, INFO ) 00024 * 00025 * .. Scalar Arguments .. 00026 * CHARACTER JOBZ, RANGE, UPLO 00027 * INTEGER IL, INFO, IU, KD, LDAB, LDQ, LDZ, M, N 00028 * REAL ABSTOL, VL, VU 00029 * .. 00030 * .. Array Arguments .. 00031 * INTEGER IFAIL( * ), IWORK( * ) 00032 * REAL RWORK( * ), W( * ) 00033 * COMPLEX AB( LDAB, * ), Q( LDQ, * ), WORK( * ), 00034 * $ Z( LDZ, * ) 00035 * .. 00036 * 00037 * 00038 *> \par Purpose: 00039 * ============= 00040 *> 00041 *> \verbatim 00042 *> 00043 *> CHBEVX computes selected eigenvalues and, optionally, eigenvectors 00044 *> of a complex Hermitian band matrix A. Eigenvalues and eigenvectors 00045 *> can be selected by specifying either a range of values or a range of 00046 *> indices for the desired eigenvalues. 00047 *> \endverbatim 00048 * 00049 * Arguments: 00050 * ========== 00051 * 00052 *> \param[in] JOBZ 00053 *> \verbatim 00054 *> JOBZ is CHARACTER*1 00055 *> = 'N': Compute eigenvalues only; 00056 *> = 'V': Compute eigenvalues and eigenvectors. 00057 *> \endverbatim 00058 *> 00059 *> \param[in] RANGE 00060 *> \verbatim 00061 *> RANGE is CHARACTER*1 00062 *> = 'A': all eigenvalues will be found; 00063 *> = 'V': all eigenvalues in the half-open interval (VL,VU] 00064 *> will be found; 00065 *> = 'I': the IL-th through IU-th eigenvalues will be found. 00066 *> \endverbatim 00067 *> 00068 *> \param[in] UPLO 00069 *> \verbatim 00070 *> UPLO is CHARACTER*1 00071 *> = 'U': Upper triangle of A is stored; 00072 *> = 'L': Lower triangle of A is stored. 00073 *> \endverbatim 00074 *> 00075 *> \param[in] N 00076 *> \verbatim 00077 *> N is INTEGER 00078 *> The order of the matrix A. N >= 0. 00079 *> \endverbatim 00080 *> 00081 *> \param[in] KD 00082 *> \verbatim 00083 *> KD is INTEGER 00084 *> The number of superdiagonals of the matrix A if UPLO = 'U', 00085 *> or the number of subdiagonals if UPLO = 'L'. KD >= 0. 00086 *> \endverbatim 00087 *> 00088 *> \param[in,out] AB 00089 *> \verbatim 00090 *> AB is COMPLEX array, dimension (LDAB, N) 00091 *> On entry, the upper or lower triangle of the Hermitian band 00092 *> matrix A, stored in the first KD+1 rows of the array. The 00093 *> j-th column of A is stored in the j-th column of the array AB 00094 *> as follows: 00095 *> if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j; 00096 *> if UPLO = 'L', AB(1+i-j,j) = A(i,j) for j<=i<=min(n,j+kd). 00097 *> 00098 *> On exit, AB is overwritten by values generated during the 00099 *> reduction to tridiagonal form. 00100 *> \endverbatim 00101 *> 00102 *> \param[in] LDAB 00103 *> \verbatim 00104 *> LDAB is INTEGER 00105 *> The leading dimension of the array AB. LDAB >= KD + 1. 00106 *> \endverbatim 00107 *> 00108 *> \param[out] Q 00109 *> \verbatim 00110 *> Q is COMPLEX array, dimension (LDQ, N) 00111 *> If JOBZ = 'V', the N-by-N unitary matrix used in the 00112 *> reduction to tridiagonal form. 00113 *> If JOBZ = 'N', the array Q is not referenced. 00114 *> \endverbatim 00115 *> 00116 *> \param[in] LDQ 00117 *> \verbatim 00118 *> LDQ is INTEGER 00119 *> The leading dimension of the array Q. If JOBZ = 'V', then 00120 *> LDQ >= max(1,N). 00121 *> \endverbatim 00122 *> 00123 *> \param[in] VL 00124 *> \verbatim 00125 *> VL is REAL 00126 *> \endverbatim 00127 *> 00128 *> \param[in] VU 00129 *> \verbatim 00130 *> VU is REAL 00131 *> If RANGE='V', the lower and upper bounds of the interval to 00132 *> be searched for eigenvalues. VL < VU. 00133 *> Not referenced if RANGE = 'A' or 'I'. 00134 *> \endverbatim 00135 *> 00136 *> \param[in] IL 00137 *> \verbatim 00138 *> IL is INTEGER 00139 *> \endverbatim 00140 *> 00141 *> \param[in] IU 00142 *> \verbatim 00143 *> IU is INTEGER 00144 *> If RANGE='I', the indices (in ascending order) of the 00145 *> smallest and largest eigenvalues to be returned. 00146 *> 1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0. 00147 *> Not referenced if RANGE = 'A' or 'V'. 00148 *> \endverbatim 00149 *> 00150 *> \param[in] ABSTOL 00151 *> \verbatim 00152 *> ABSTOL is REAL 00153 *> The absolute error tolerance for the eigenvalues. 00154 *> An approximate eigenvalue is accepted as converged 00155 *> when it is determined to lie in an interval [a,b] 00156 *> of width less than or equal to 00157 *> 00158 *> ABSTOL + EPS * max( |a|,|b| ) , 00159 *> 00160 *> where EPS is the machine precision. If ABSTOL is less than 00161 *> or equal to zero, then EPS*|T| will be used in its place, 00162 *> where |T| is the 1-norm of the tridiagonal matrix obtained 00163 *> by reducing AB to tridiagonal form. 00164 *> 00165 *> Eigenvalues will be computed most accurately when ABSTOL is 00166 *> set to twice the underflow threshold 2*SLAMCH('S'), not zero. 00167 *> If this routine returns with INFO>0, indicating that some 00168 *> eigenvectors did not converge, try setting ABSTOL to 00169 *> 2*SLAMCH('S'). 00170 *> 00171 *> See "Computing Small Singular Values of Bidiagonal Matrices 00172 *> with Guaranteed High Relative Accuracy," by Demmel and 00173 *> Kahan, LAPACK Working Note #3. 00174 *> \endverbatim 00175 *> 00176 *> \param[out] M 00177 *> \verbatim 00178 *> M is INTEGER 00179 *> The total number of eigenvalues found. 0 <= M <= N. 00180 *> If RANGE = 'A', M = N, and if RANGE = 'I', M = IU-IL+1. 00181 *> \endverbatim 00182 *> 00183 *> \param[out] W 00184 *> \verbatim 00185 *> W is REAL array, dimension (N) 00186 *> The first M elements contain the selected eigenvalues in 00187 *> ascending order. 00188 *> \endverbatim 00189 *> 00190 *> \param[out] Z 00191 *> \verbatim 00192 *> Z is COMPLEX array, dimension (LDZ, max(1,M)) 00193 *> If JOBZ = 'V', then if INFO = 0, the first M columns of Z 00194 *> contain the orthonormal eigenvectors of the matrix A 00195 *> corresponding to the selected eigenvalues, with the i-th 00196 *> column of Z holding the eigenvector associated with W(i). 00197 *> If an eigenvector fails to converge, then that column of Z 00198 *> contains the latest approximation to the eigenvector, and the 00199 *> index of the eigenvector is returned in IFAIL. 00200 *> If JOBZ = 'N', then Z is not referenced. 00201 *> Note: the user must ensure that at least max(1,M) columns are 00202 *> supplied in the array Z; if RANGE = 'V', the exact value of M 00203 *> is not known in advance and an upper bound must be used. 00204 *> \endverbatim 00205 *> 00206 *> \param[in] LDZ 00207 *> \verbatim 00208 *> LDZ is INTEGER 00209 *> The leading dimension of the array Z. LDZ >= 1, and if 00210 *> JOBZ = 'V', LDZ >= max(1,N). 00211 *> \endverbatim 00212 *> 00213 *> \param[out] WORK 00214 *> \verbatim 00215 *> WORK is COMPLEX array, dimension (N) 00216 *> \endverbatim 00217 *> 00218 *> \param[out] RWORK 00219 *> \verbatim 00220 *> RWORK is REAL array, dimension (7*N) 00221 *> \endverbatim 00222 *> 00223 *> \param[out] IWORK 00224 *> \verbatim 00225 *> IWORK is INTEGER array, dimension (5*N) 00226 *> \endverbatim 00227 *> 00228 *> \param[out] IFAIL 00229 *> \verbatim 00230 *> IFAIL is INTEGER array, dimension (N) 00231 *> If JOBZ = 'V', then if INFO = 0, the first M elements of 00232 *> IFAIL are zero. If INFO > 0, then IFAIL contains the 00233 *> indices of the eigenvectors that failed to converge. 00234 *> If JOBZ = 'N', then IFAIL is not referenced. 00235 *> \endverbatim 00236 *> 00237 *> \param[out] INFO 00238 *> \verbatim 00239 *> INFO is INTEGER 00240 *> = 0: successful exit 00241 *> < 0: if INFO = -i, the i-th argument had an illegal value 00242 *> > 0: if INFO = i, then i eigenvectors failed to converge. 00243 *> Their indices are stored in array IFAIL. 00244 *> \endverbatim 00245 * 00246 * Authors: 00247 * ======== 00248 * 00249 *> \author Univ. of Tennessee 00250 *> \author Univ. of California Berkeley 00251 *> \author Univ. of Colorado Denver 00252 *> \author NAG Ltd. 00253 * 00254 *> \date November 2011 00255 * 00256 *> \ingroup complexOTHEReigen 00257 * 00258 * ===================================================================== 00259 SUBROUTINE CHBEVX( JOBZ, RANGE, UPLO, N, KD, AB, LDAB, Q, LDQ, VL, 00260 $ VU, IL, IU, ABSTOL, M, W, Z, LDZ, WORK, RWORK, 00261 $ IWORK, IFAIL, INFO ) 00262 * 00263 * -- LAPACK driver routine (version 3.4.0) -- 00264 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00265 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00266 * November 2011 00267 * 00268 * .. Scalar Arguments .. 00269 CHARACTER JOBZ, RANGE, UPLO 00270 INTEGER IL, INFO, IU, KD, LDAB, LDQ, LDZ, M, N 00271 REAL ABSTOL, VL, VU 00272 * .. 00273 * .. Array Arguments .. 00274 INTEGER IFAIL( * ), IWORK( * ) 00275 REAL RWORK( * ), W( * ) 00276 COMPLEX AB( LDAB, * ), Q( LDQ, * ), WORK( * ), 00277 $ Z( LDZ, * ) 00278 * .. 00279 * 00280 * ===================================================================== 00281 * 00282 * .. Parameters .. 00283 REAL ZERO, ONE 00284 PARAMETER ( ZERO = 0.0E0, ONE = 1.0E0 ) 00285 COMPLEX CZERO, CONE 00286 PARAMETER ( CZERO = ( 0.0E0, 0.0E0 ), 00287 $ CONE = ( 1.0E0, 0.0E0 ) ) 00288 * .. 00289 * .. Local Scalars .. 00290 LOGICAL ALLEIG, INDEIG, LOWER, TEST, VALEIG, WANTZ 00291 CHARACTER ORDER 00292 INTEGER I, IINFO, IMAX, INDD, INDE, INDEE, INDIBL, 00293 $ INDISP, INDIWK, INDRWK, INDWRK, ISCALE, ITMP1, 00294 $ J, JJ, NSPLIT 00295 REAL ABSTLL, ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, 00296 $ SIGMA, SMLNUM, TMP1, VLL, VUU 00297 COMPLEX CTMP1 00298 * .. 00299 * .. External Functions .. 00300 LOGICAL LSAME 00301 REAL CLANHB, SLAMCH 00302 EXTERNAL LSAME, CLANHB, SLAMCH 00303 * .. 00304 * .. External Subroutines .. 00305 EXTERNAL CCOPY, CGEMV, CHBTRD, CLACPY, CLASCL, CSTEIN, 00306 $ CSTEQR, CSWAP, SCOPY, SSCAL, SSTEBZ, SSTERF, 00307 $ XERBLA 00308 * .. 00309 * .. Intrinsic Functions .. 00310 INTRINSIC MAX, MIN, REAL, SQRT 00311 * .. 00312 * .. Executable Statements .. 00313 * 00314 * Test the input parameters. 00315 * 00316 WANTZ = LSAME( JOBZ, 'V' ) 00317 ALLEIG = LSAME( RANGE, 'A' ) 00318 VALEIG = LSAME( RANGE, 'V' ) 00319 INDEIG = LSAME( RANGE, 'I' ) 00320 LOWER = LSAME( UPLO, 'L' ) 00321 * 00322 INFO = 0 00323 IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN 00324 INFO = -1 00325 ELSE IF( .NOT.( ALLEIG .OR. VALEIG .OR. INDEIG ) ) THEN 00326 INFO = -2 00327 ELSE IF( .NOT.( LOWER .OR. LSAME( UPLO, 'U' ) ) ) THEN 00328 INFO = -3 00329 ELSE IF( N.LT.0 ) THEN 00330 INFO = -4 00331 ELSE IF( KD.LT.0 ) THEN 00332 INFO = -5 00333 ELSE IF( LDAB.LT.KD+1 ) THEN 00334 INFO = -7 00335 ELSE IF( WANTZ .AND. LDQ.LT.MAX( 1, N ) ) THEN 00336 INFO = -9 00337 ELSE 00338 IF( VALEIG ) THEN 00339 IF( N.GT.0 .AND. VU.LE.VL ) 00340 $ INFO = -11 00341 ELSE IF( INDEIG ) THEN 00342 IF( IL.LT.1 .OR. IL.GT.MAX( 1, N ) ) THEN 00343 INFO = -12 00344 ELSE IF( IU.LT.MIN( N, IL ) .OR. IU.GT.N ) THEN 00345 INFO = -13 00346 END IF 00347 END IF 00348 END IF 00349 IF( INFO.EQ.0 ) THEN 00350 IF( LDZ.LT.1 .OR. ( WANTZ .AND. LDZ.LT.N ) ) 00351 $ INFO = -18 00352 END IF 00353 * 00354 IF( INFO.NE.0 ) THEN 00355 CALL XERBLA( 'CHBEVX', -INFO ) 00356 RETURN 00357 END IF 00358 * 00359 * Quick return if possible 00360 * 00361 M = 0 00362 IF( N.EQ.0 ) 00363 $ RETURN 00364 * 00365 IF( N.EQ.1 ) THEN 00366 M = 1 00367 IF( LOWER ) THEN 00368 CTMP1 = AB( 1, 1 ) 00369 ELSE 00370 CTMP1 = AB( KD+1, 1 ) 00371 END IF 00372 TMP1 = REAL( CTMP1 ) 00373 IF( VALEIG ) THEN 00374 IF( .NOT.( VL.LT.TMP1 .AND. VU.GE.TMP1 ) ) 00375 $ M = 0 00376 END IF 00377 IF( M.EQ.1 ) THEN 00378 W( 1 ) = CTMP1 00379 IF( WANTZ ) 00380 $ Z( 1, 1 ) = CONE 00381 END IF 00382 RETURN 00383 END IF 00384 * 00385 * Get machine constants. 00386 * 00387 SAFMIN = SLAMCH( 'Safe minimum' ) 00388 EPS = SLAMCH( 'Precision' ) 00389 SMLNUM = SAFMIN / EPS 00390 BIGNUM = ONE / SMLNUM 00391 RMIN = SQRT( SMLNUM ) 00392 RMAX = MIN( SQRT( BIGNUM ), ONE / SQRT( SQRT( SAFMIN ) ) ) 00393 * 00394 * Scale matrix to allowable range, if necessary. 00395 * 00396 ISCALE = 0 00397 ABSTLL = ABSTOL 00398 IF ( VALEIG ) THEN 00399 VLL = VL 00400 VUU = VU 00401 ELSE 00402 VLL = ZERO 00403 VUU = ZERO 00404 ENDIF 00405 ANRM = CLANHB( 'M', UPLO, N, KD, AB, LDAB, RWORK ) 00406 IF( ANRM.GT.ZERO .AND. ANRM.LT.RMIN ) THEN 00407 ISCALE = 1 00408 SIGMA = RMIN / ANRM 00409 ELSE IF( ANRM.GT.RMAX ) THEN 00410 ISCALE = 1 00411 SIGMA = RMAX / ANRM 00412 END IF 00413 IF( ISCALE.EQ.1 ) THEN 00414 IF( LOWER ) THEN 00415 CALL CLASCL( 'B', KD, KD, ONE, SIGMA, N, N, AB, LDAB, INFO ) 00416 ELSE 00417 CALL CLASCL( 'Q', KD, KD, ONE, SIGMA, N, N, AB, LDAB, INFO ) 00418 END IF 00419 IF( ABSTOL.GT.0 ) 00420 $ ABSTLL = ABSTOL*SIGMA 00421 IF( VALEIG ) THEN 00422 VLL = VL*SIGMA 00423 VUU = VU*SIGMA 00424 END IF 00425 END IF 00426 * 00427 * Call CHBTRD to reduce Hermitian band matrix to tridiagonal form. 00428 * 00429 INDD = 1 00430 INDE = INDD + N 00431 INDRWK = INDE + N 00432 INDWRK = 1 00433 CALL CHBTRD( JOBZ, UPLO, N, KD, AB, LDAB, RWORK( INDD ), 00434 $ RWORK( INDE ), Q, LDQ, WORK( INDWRK ), IINFO ) 00435 * 00436 * If all eigenvalues are desired and ABSTOL is less than or equal 00437 * to zero, then call SSTERF or CSTEQR. If this fails for some 00438 * eigenvalue, then try SSTEBZ. 00439 * 00440 TEST = .FALSE. 00441 IF (INDEIG) THEN 00442 IF (IL.EQ.1 .AND. IU.EQ.N) THEN 00443 TEST = .TRUE. 00444 END IF 00445 END IF 00446 IF ((ALLEIG .OR. TEST) .AND. (ABSTOL.LE.ZERO)) THEN 00447 CALL SCOPY( N, RWORK( INDD ), 1, W, 1 ) 00448 INDEE = INDRWK + 2*N 00449 IF( .NOT.WANTZ ) THEN 00450 CALL SCOPY( N-1, RWORK( INDE ), 1, RWORK( INDEE ), 1 ) 00451 CALL SSTERF( N, W, RWORK( INDEE ), INFO ) 00452 ELSE 00453 CALL CLACPY( 'A', N, N, Q, LDQ, Z, LDZ ) 00454 CALL SCOPY( N-1, RWORK( INDE ), 1, RWORK( INDEE ), 1 ) 00455 CALL CSTEQR( JOBZ, N, W, RWORK( INDEE ), Z, LDZ, 00456 $ RWORK( INDRWK ), INFO ) 00457 IF( INFO.EQ.0 ) THEN 00458 DO 10 I = 1, N 00459 IFAIL( I ) = 0 00460 10 CONTINUE 00461 END IF 00462 END IF 00463 IF( INFO.EQ.0 ) THEN 00464 M = N 00465 GO TO 30 00466 END IF 00467 INFO = 0 00468 END IF 00469 * 00470 * Otherwise, call SSTEBZ and, if eigenvectors are desired, CSTEIN. 00471 * 00472 IF( WANTZ ) THEN 00473 ORDER = 'B' 00474 ELSE 00475 ORDER = 'E' 00476 END IF 00477 INDIBL = 1 00478 INDISP = INDIBL + N 00479 INDIWK = INDISP + N 00480 CALL SSTEBZ( RANGE, ORDER, N, VLL, VUU, IL, IU, ABSTLL, 00481 $ RWORK( INDD ), RWORK( INDE ), M, NSPLIT, W, 00482 $ IWORK( INDIBL ), IWORK( INDISP ), RWORK( INDRWK ), 00483 $ IWORK( INDIWK ), INFO ) 00484 * 00485 IF( WANTZ ) THEN 00486 CALL CSTEIN( N, RWORK( INDD ), RWORK( INDE ), M, W, 00487 $ IWORK( INDIBL ), IWORK( INDISP ), Z, LDZ, 00488 $ RWORK( INDRWK ), IWORK( INDIWK ), IFAIL, INFO ) 00489 * 00490 * Apply unitary matrix used in reduction to tridiagonal 00491 * form to eigenvectors returned by CSTEIN. 00492 * 00493 DO 20 J = 1, M 00494 CALL CCOPY( N, Z( 1, J ), 1, WORK( 1 ), 1 ) 00495 CALL CGEMV( 'N', N, N, CONE, Q, LDQ, WORK, 1, CZERO, 00496 $ Z( 1, J ), 1 ) 00497 20 CONTINUE 00498 END IF 00499 * 00500 * If matrix was scaled, then rescale eigenvalues appropriately. 00501 * 00502 30 CONTINUE 00503 IF( ISCALE.EQ.1 ) THEN 00504 IF( INFO.EQ.0 ) THEN 00505 IMAX = M 00506 ELSE 00507 IMAX = INFO - 1 00508 END IF 00509 CALL SSCAL( IMAX, ONE / SIGMA, W, 1 ) 00510 END IF 00511 * 00512 * If eigenvalues are not in order, then sort them, along with 00513 * eigenvectors. 00514 * 00515 IF( WANTZ ) THEN 00516 DO 50 J = 1, M - 1 00517 I = 0 00518 TMP1 = W( J ) 00519 DO 40 JJ = J + 1, M 00520 IF( W( JJ ).LT.TMP1 ) THEN 00521 I = JJ 00522 TMP1 = W( JJ ) 00523 END IF 00524 40 CONTINUE 00525 * 00526 IF( I.NE.0 ) THEN 00527 ITMP1 = IWORK( INDIBL+I-1 ) 00528 W( I ) = W( J ) 00529 IWORK( INDIBL+I-1 ) = IWORK( INDIBL+J-1 ) 00530 W( J ) = TMP1 00531 IWORK( INDIBL+J-1 ) = ITMP1 00532 CALL CSWAP( N, Z( 1, I ), 1, Z( 1, J ), 1 ) 00533 IF( INFO.NE.0 ) THEN 00534 ITMP1 = IFAIL( I ) 00535 IFAIL( I ) = IFAIL( J ) 00536 IFAIL( J ) = ITMP1 00537 END IF 00538 END IF 00539 50 CONTINUE 00540 END IF 00541 * 00542 RETURN 00543 * 00544 * End of CHBEVX 00545 * 00546 END