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