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