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