 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
00001 SUBROUTINE DSTEVX( JOBZ, RANGE, N, D, E, VL, VU, IL, IU, ABSTOL, 00002 $ M, W, Z, LDZ, WORK, IWORK, IFAIL, INFO ) 00003 * 00004 * -- LAPACK driver routine (version 3.2) -- 00005 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00006 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00007 * November 2006 00008 * 00009 * .. Scalar Arguments .. 00010 CHARACTER JOBZ, RANGE 00011 INTEGER IL, INFO, IU, LDZ, M, N 00012 DOUBLE PRECISION ABSTOL, VL, VU 00013 * .. 00014 * .. Array Arguments .. 00015 INTEGER IFAIL( * ), IWORK( * ) 00016 DOUBLE PRECISION D( * ), E( * ), W( * ), WORK( * ), Z( LDZ, * ) 00017 * .. 00018 * 00019 * Purpose 00020 * ======= 00021 * 00022 * DSTEVX computes selected eigenvalues and, optionally, eigenvectors 00023 * of a real symmetric tridiagonal matrix A. Eigenvalues and 00024 * eigenvectors can be selected by specifying either a range of values 00025 * or a range of indices for the desired eigenvalues. 00026 * 00027 * Arguments 00028 * ========= 00029 * 00030 * JOBZ (input) CHARACTER*1 00031 * = 'N': Compute eigenvalues only; 00032 * = 'V': Compute eigenvalues and eigenvectors. 00033 * 00034 * RANGE (input) CHARACTER*1 00035 * = 'A': all eigenvalues will be found. 00036 * = 'V': all eigenvalues in the half-open interval (VL,VU] 00037 * will be found. 00038 * = 'I': the IL-th through IU-th eigenvalues will be found. 00039 * 00040 * N (input) INTEGER 00041 * The order of the matrix. N >= 0. 00042 * 00043 * D (input/output) DOUBLE PRECISION array, dimension (N) 00044 * On entry, the n diagonal elements of the tridiagonal matrix 00045 * A. 00046 * On exit, D may be multiplied by a constant factor chosen 00047 * to avoid over/underflow in computing the eigenvalues. 00048 * 00049 * E (input/output) DOUBLE PRECISION array, dimension (max(1,N-1)) 00050 * On entry, the (n-1) subdiagonal elements of the tridiagonal 00051 * matrix A in elements 1 to N-1 of E. 00052 * On exit, E may be multiplied by a constant factor chosen 00053 * to avoid over/underflow in computing the eigenvalues. 00054 * 00055 * VL (input) DOUBLE PRECISION 00056 * VU (input) DOUBLE PRECISION 00057 * If RANGE='V', the lower and upper bounds of the interval to 00058 * be searched for eigenvalues. VL < VU. 00059 * Not referenced if RANGE = 'A' or 'I'. 00060 * 00061 * IL (input) INTEGER 00062 * IU (input) INTEGER 00063 * If RANGE='I', the indices (in ascending order) of the 00064 * smallest and largest eigenvalues to be returned. 00065 * 1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0. 00066 * Not referenced if RANGE = 'A' or 'V'. 00067 * 00068 * ABSTOL (input) DOUBLE PRECISION 00069 * The absolute error tolerance for the eigenvalues. 00070 * An approximate eigenvalue is accepted as converged 00071 * when it is determined to lie in an interval [a,b] 00072 * of width less than or equal to 00073 * 00074 * ABSTOL + EPS * max( |a|,|b| ) , 00075 * 00076 * where EPS is the machine precision. If ABSTOL is less 00077 * than or equal to zero, then EPS*|T| will be used in 00078 * its place, where |T| is the 1-norm of the tridiagonal 00079 * matrix. 00080 * 00081 * Eigenvalues will be computed most accurately when ABSTOL is 00082 * set to twice the underflow threshold 2*DLAMCH('S'), not zero. 00083 * If this routine returns with INFO>0, indicating that some 00084 * eigenvectors did not converge, try setting ABSTOL to 00085 * 2*DLAMCH('S'). 00086 * 00087 * See "Computing Small Singular Values of Bidiagonal Matrices 00088 * with Guaranteed High Relative Accuracy," by Demmel and 00089 * Kahan, LAPACK Working Note #3. 00090 * 00091 * M (output) INTEGER 00092 * The total number of eigenvalues found. 0 <= M <= N. 00093 * If RANGE = 'A', M = N, and if RANGE = 'I', M = IU-IL+1. 00094 * 00095 * W (output) DOUBLE PRECISION array, dimension (N) 00096 * The first M elements contain the selected eigenvalues in 00097 * ascending order. 00098 * 00099 * Z (output) DOUBLE PRECISION array, dimension (LDZ, max(1,M) ) 00100 * If JOBZ = 'V', then if INFO = 0, the first M columns of Z 00101 * contain the orthonormal eigenvectors of the matrix A 00102 * corresponding to the selected eigenvalues, with the i-th 00103 * column of Z holding the eigenvector associated with W(i). 00104 * If an eigenvector fails to converge (INFO > 0), then that 00105 * column of Z contains the latest approximation to the 00106 * eigenvector, and the index of the eigenvector is returned 00107 * in IFAIL. If JOBZ = 'N', then Z is not referenced. 00108 * Note: the user must ensure that at least max(1,M) columns are 00109 * supplied in the array Z; if RANGE = 'V', the exact value of M 00110 * is not known in advance and an upper bound must be used. 00111 * 00112 * LDZ (input) INTEGER 00113 * The leading dimension of the array Z. LDZ >= 1, and if 00114 * JOBZ = 'V', LDZ >= max(1,N). 00115 * 00116 * WORK (workspace) DOUBLE PRECISION array, dimension (5*N) 00117 * 00118 * IWORK (workspace) INTEGER array, dimension (5*N) 00119 * 00120 * IFAIL (output) INTEGER array, dimension (N) 00121 * If JOBZ = 'V', then if INFO = 0, the first M elements of 00122 * IFAIL are zero. If INFO > 0, then IFAIL contains the 00123 * indices of the eigenvectors that failed to converge. 00124 * If JOBZ = 'N', then IFAIL is not referenced. 00125 * 00126 * INFO (output) INTEGER 00127 * = 0: successful exit 00128 * < 0: if INFO = -i, the i-th argument had an illegal value 00129 * > 0: if INFO = i, then i eigenvectors failed to converge. 00130 * Their indices are stored in array IFAIL. 00131 * 00132 * ===================================================================== 00133 * 00134 * .. Parameters .. 00135 DOUBLE PRECISION ZERO, ONE 00136 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 ) 00137 * .. 00138 * .. Local Scalars .. 00139 LOGICAL ALLEIG, INDEIG, TEST, VALEIG, WANTZ 00140 CHARACTER ORDER 00141 INTEGER I, IMAX, INDIBL, INDISP, INDIWO, INDWRK, 00142 $ ISCALE, ITMP1, J, JJ, NSPLIT 00143 DOUBLE PRECISION BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA, SMLNUM, 00144 $ TMP1, TNRM, VLL, VUU 00145 * .. 00146 * .. External Functions .. 00147 LOGICAL LSAME 00148 DOUBLE PRECISION DLAMCH, DLANST 00149 EXTERNAL LSAME, DLAMCH, DLANST 00150 * .. 00151 * .. External Subroutines .. 00152 EXTERNAL DCOPY, DSCAL, DSTEBZ, DSTEIN, DSTEQR, DSTERF, 00153 $ DSWAP, XERBLA 00154 * .. 00155 * .. Intrinsic Functions .. 00156 INTRINSIC MAX, MIN, SQRT 00157 * .. 00158 * .. Executable Statements .. 00159 * 00160 * Test the input parameters. 00161 * 00162 WANTZ = LSAME( JOBZ, 'V' ) 00163 ALLEIG = LSAME( RANGE, 'A' ) 00164 VALEIG = LSAME( RANGE, 'V' ) 00165 INDEIG = LSAME( RANGE, 'I' ) 00166 * 00167 INFO = 0 00168 IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN 00169 INFO = -1 00170 ELSE IF( .NOT.( ALLEIG .OR. VALEIG .OR. INDEIG ) ) THEN 00171 INFO = -2 00172 ELSE IF( N.LT.0 ) THEN 00173 INFO = -3 00174 ELSE 00175 IF( VALEIG ) THEN 00176 IF( N.GT.0 .AND. VU.LE.VL ) 00177 $ INFO = -7 00178 ELSE IF( INDEIG ) THEN 00179 IF( IL.LT.1 .OR. IL.GT.MAX( 1, N ) ) THEN 00180 INFO = -8 00181 ELSE IF( IU.LT.MIN( N, IL ) .OR. IU.GT.N ) THEN 00182 INFO = -9 00183 END IF 00184 END IF 00185 END IF 00186 IF( INFO.EQ.0 ) THEN 00187 IF( LDZ.LT.1 .OR. ( WANTZ .AND. LDZ.LT.N ) ) 00188 $ INFO = -14 00189 END IF 00190 * 00191 IF( INFO.NE.0 ) THEN 00192 CALL XERBLA( 'DSTEVX', -INFO ) 00193 RETURN 00194 END IF 00195 * 00196 * Quick return if possible 00197 * 00198 M = 0 00199 IF( N.EQ.0 ) 00200 $ RETURN 00201 * 00202 IF( N.EQ.1 ) THEN 00203 IF( ALLEIG .OR. INDEIG ) THEN 00204 M = 1 00205 W( 1 ) = D( 1 ) 00206 ELSE 00207 IF( VL.LT.D( 1 ) .AND. VU.GE.D( 1 ) ) THEN 00208 M = 1 00209 W( 1 ) = D( 1 ) 00210 END IF 00211 END IF 00212 IF( WANTZ ) 00213 $ Z( 1, 1 ) = ONE 00214 RETURN 00215 END IF 00216 * 00217 * Get machine constants. 00218 * 00219 SAFMIN = DLAMCH( 'Safe minimum' ) 00220 EPS = DLAMCH( 'Precision' ) 00221 SMLNUM = SAFMIN / EPS 00222 BIGNUM = ONE / SMLNUM 00223 RMIN = SQRT( SMLNUM ) 00224 RMAX = MIN( SQRT( BIGNUM ), ONE / SQRT( SQRT( SAFMIN ) ) ) 00225 * 00226 * Scale matrix to allowable range, if necessary. 00227 * 00228 ISCALE = 0 00229 IF( VALEIG ) THEN 00230 VLL = VL 00231 VUU = VU 00232 ELSE 00233 VLL = ZERO 00234 VUU = ZERO 00235 END IF 00236 TNRM = DLANST( 'M', N, D, E ) 00237 IF( TNRM.GT.ZERO .AND. TNRM.LT.RMIN ) THEN 00238 ISCALE = 1 00239 SIGMA = RMIN / TNRM 00240 ELSE IF( TNRM.GT.RMAX ) THEN 00241 ISCALE = 1 00242 SIGMA = RMAX / TNRM 00243 END IF 00244 IF( ISCALE.EQ.1 ) THEN 00245 CALL DSCAL( N, SIGMA, D, 1 ) 00246 CALL DSCAL( N-1, SIGMA, E( 1 ), 1 ) 00247 IF( VALEIG ) THEN 00248 VLL = VL*SIGMA 00249 VUU = VU*SIGMA 00250 END IF 00251 END IF 00252 * 00253 * If all eigenvalues are desired and ABSTOL is less than zero, then 00254 * call DSTERF or SSTEQR. If this fails for some eigenvalue, then 00255 * try DSTEBZ. 00256 * 00257 TEST = .FALSE. 00258 IF( INDEIG ) THEN 00259 IF( IL.EQ.1 .AND. IU.EQ.N ) THEN 00260 TEST = .TRUE. 00261 END IF 00262 END IF 00263 IF( ( ALLEIG .OR. TEST ) .AND. ( ABSTOL.LE.ZERO ) ) THEN 00264 CALL DCOPY( N, D, 1, W, 1 ) 00265 CALL DCOPY( N-1, E( 1 ), 1, WORK( 1 ), 1 ) 00266 INDWRK = N + 1 00267 IF( .NOT.WANTZ ) THEN 00268 CALL DSTERF( N, W, WORK, INFO ) 00269 ELSE 00270 CALL DSTEQR( 'I', N, W, WORK, Z, LDZ, WORK( INDWRK ), INFO ) 00271 IF( INFO.EQ.0 ) THEN 00272 DO 10 I = 1, N 00273 IFAIL( I ) = 0 00274 10 CONTINUE 00275 END IF 00276 END IF 00277 IF( INFO.EQ.0 ) THEN 00278 M = N 00279 GO TO 20 00280 END IF 00281 INFO = 0 00282 END IF 00283 * 00284 * Otherwise, call DSTEBZ and, if eigenvectors are desired, SSTEIN. 00285 * 00286 IF( WANTZ ) THEN 00287 ORDER = 'B' 00288 ELSE 00289 ORDER = 'E' 00290 END IF 00291 INDWRK = 1 00292 INDIBL = 1 00293 INDISP = INDIBL + N 00294 INDIWO = INDISP + N 00295 CALL DSTEBZ( RANGE, ORDER, N, VLL, VUU, IL, IU, ABSTOL, D, E, M, 00296 $ NSPLIT, W, IWORK( INDIBL ), IWORK( INDISP ), 00297 $ WORK( INDWRK ), IWORK( INDIWO ), INFO ) 00298 * 00299 IF( WANTZ ) THEN 00300 CALL DSTEIN( N, D, E, M, W, IWORK( INDIBL ), IWORK( INDISP ), 00301 $ Z, LDZ, WORK( INDWRK ), IWORK( INDIWO ), IFAIL, 00302 $ INFO ) 00303 END IF 00304 * 00305 * If matrix was scaled, then rescale eigenvalues appropriately. 00306 * 00307 20 CONTINUE 00308 IF( ISCALE.EQ.1 ) THEN 00309 IF( INFO.EQ.0 ) THEN 00310 IMAX = M 00311 ELSE 00312 IMAX = INFO - 1 00313 END IF 00314 CALL DSCAL( IMAX, ONE / SIGMA, W, 1 ) 00315 END IF 00316 * 00317 * If eigenvalues are not in order, then sort them, along with 00318 * eigenvectors. 00319 * 00320 IF( WANTZ ) THEN 00321 DO 40 J = 1, M - 1 00322 I = 0 00323 TMP1 = W( J ) 00324 DO 30 JJ = J + 1, M 00325 IF( W( JJ ).LT.TMP1 ) THEN 00326 I = JJ 00327 TMP1 = W( JJ ) 00328 END IF 00329 30 CONTINUE 00330 * 00331 IF( I.NE.0 ) THEN 00332 ITMP1 = IWORK( INDIBL+I-1 ) 00333 W( I ) = W( J ) 00334 IWORK( INDIBL+I-1 ) = IWORK( INDIBL+J-1 ) 00335 W( J ) = TMP1 00336 IWORK( INDIBL+J-1 ) = ITMP1 00337 CALL DSWAP( N, Z( 1, I ), 1, Z( 1, J ), 1 ) 00338 IF( INFO.NE.0 ) THEN 00339 ITMP1 = IFAIL( I ) 00340 IFAIL( I ) = IFAIL( J ) 00341 IFAIL( J ) = ITMP1 00342 END IF 00343 END IF 00344 40 CONTINUE 00345 END IF 00346 * 00347 RETURN 00348 * 00349 * End of DSTEVX 00350 * 00351 END
1.7.3 
