1* 2* $Id$ 3* 4#include "blas_lapackf.h" 5*====================================================================== 6* 7* DISCLAIMER 8* 9* This material was prepared as an account of work sponsored by an 10* agency of the United States Government. Neither the United States 11* Government nor the United States Department of Energy, nor Battelle, 12* nor any of their employees, MAKES ANY WARRANTY, EXPRESS OR IMPLIED, OR 13* ASSUMES ANY LEGAL LIABILITY OR RESPONSIBILITY FOR THE ACCURACY, 14* COMPLETENESS, OR USEFULNESS OF ANY INFORMATION, APPARATUS, PRODUCT, 15* SOFTWARE, OR PROCESS DISCLOSED, OR REPRESENTS THAT ITS USE WOULD NOT 16* INFRINGE PRIVATELY OWNED RIGHTS. 17* 18* ACKNOWLEDGMENT 19* 20* This software and its documentation were produced with Government 21* support under Contract Number DE-AC06-76RLO-1830 awarded by the United 22* States Department of Energy. The Government retains a paid-up 23* non-exclusive, irrevocable worldwide license to reproduce, prepare 24* derivative works, perform publicly and display publicly by or for the 25* Government, including the right to distribute to other Government 26* contractors. 27* 28*====================================================================== 29* 30* -- PEIGS routine (version 2.1) -- 31* Pacific Northwest Laboratory 32* July 28, 1995 33* 34*====================================================================== 35 SUBROUTINE DSPEVX2( JOBZ, RANGE, UPLO, N, AP, VL, VU, IL, IU, 36 $ ABSTOL, M, W, Z, LDZ, WORK, IWORK, IFAIL, 37 $ INFO ) 38* 39* -- LAPACK driver routine (version 2.0) -- 40* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., 41* Courant Institute, Argonne National Lab, and Rice University 42* September 30, 1994 43* 44* .. Scalar Arguments .. 45 CHARACTER JOBZ, RANGE, UPLO 46 INTEGER IL, INFO, IU, LDZ, M, N 47 DOUBLE PRECISION ABSTOL, VL, VU 48* .. 49* .. Array Arguments .. 50 INTEGER IFAIL( * ), IWORK( * ) 51 DOUBLE PRECISION AP( * ), W( * ), WORK( * ), Z( LDZ, * ) 52* .. 53* 54* Purpose 55* ======= 56* 57* DSPEVX computes selected eigenvalues and, optionally, eigenvectors 58* of a real symmetric matrix A in packed storage. Eigenvalues/vectors 59* can be selected by specifying either a range of values or a range of 60* indices for the desired eigenvalues. 61* 62* Arguments 63* ========= 64* 65* JOBZ (input) CHARACTER*1 66* = 'N': Compute eigenvalues only; 67* = 'V': Compute eigenvalues and eigenvectors. 68* 69* RANGE (input) CHARACTER*1 70* = 'A': all eigenvalues will be found; 71* = 'V': all eigenvalues in the half-open interval (VL,VU] 72* will be found; 73* = 'I': the IL-th through IU-th eigenvalues will be found. 74* 75* UPLO (input) CHARACTER*1 76* = 'U': Upper triangle of A is stored; 77* = 'L': Lower triangle of A is stored. 78* 79* N (input) INTEGER 80* The order of the matrix A. N >= 0. 81* 82* AP (input/output) DOUBLE PRECISION array, dimension (N*(N+1)/2) 83* On entry, the upper or lower triangle of the symmetric matrix 84* A, packed columnwise in a linear array. The j-th column of A 85* is stored in the array AP as follows: 86* if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j; 87* if UPLO = 'L', AP(i + (j-1)*(2*n-j)/2) = A(i,j) for j<=i<=n. 88* 89* On exit, AP is overwritten by values generated during the 90* reduction to tridiagonal form. If UPLO = 'U', the diagonal 91* and first superdiagonal of the tridiagonal matrix T overwrite 92* the corresponding elements of A, and if UPLO = 'L', the 93* diagonal and first subdiagonal of T overwrite the 94* corresponding elements of A. 95* 96* VL (input) DOUBLE PRECISION 97* VU (input) DOUBLE PRECISION 98* If RANGE='V', the lower and upper bounds of the interval to 99* be searched for eigenvalues. VL < VU. 100* Not referenced if RANGE = 'A' or 'I'. 101* 102* IL (input) INTEGER 103* IU (input) INTEGER 104* If RANGE='I', the indices (in ascending order) of the 105* smallest and largest eigenvalues to be returned. 106* 1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0. 107* Not referenced if RANGE = 'A' or 'V'. 108* 109* ABSTOL (input) DOUBLE PRECISION 110* The absolute error tolerance for the eigenvalues. 111* An approximate eigenvalue is accepted as converged 112* when it is determined to lie in an interval [a,b] 113* of width less than or equal to 114* 115* ABSTOL + EPS * max( |a|,|b| ) , 116* 117* where EPS is the machine precision. If ABSTOL is less than 118* or equal to zero, then EPS*|T| will be used in its place, 119* where |T| is the 1-norm of the tridiagonal matrix obtained 120* by reducing AP to tridiagonal form. 121* 122* Eigenvalues will be computed most accurately when ABSTOL is 123* set to twice the underflow threshold 2*DLAMCH('S'), not zero. 124* If this routine returns with INFO>0, indicating that some 125* eigenvectors did not converge, try setting ABSTOL to 126* 2*DLAMCH('S'). 127* 128* See "Computing Small Singular Values of Bidiagonal Matrices 129* with Guaranteed High Relative Accuracy," by Demmel and 130* Kahan, LAPACK Working Note #3. 131* 132* M (output) INTEGER 133* The total number of eigenvalues found. 0 <= M <= N. 134* If RANGE = 'A', M = N, and if RANGE = 'I', M = IU-IL+1. 135* 136* W (output) DOUBLE PRECISION array, dimension (N) 137* If INFO = 0, the selected eigenvalues in ascending order. 138* 139* Z (output) DOUBLE PRECISION array, dimension (LDZ, max(1,M)) 140* If JOBZ = 'V', then if INFO = 0, the first M columns of Z 141* contain the orthonormal eigenvectors of the matrix A 142* corresponding to the selected eigenvalues, with the i-th 143* column of Z holding the eigenvector associated with W(i). 144* If an eigenvector fails to converge, then that column of Z 145* contains the latest approximation to the eigenvector, and the 146* index of the eigenvector is returned in IFAIL. 147* If JOBZ = 'N', then Z is not referenced. 148* Note: the user must ensure that at least max(1,M) columns are 149* supplied in the array Z; if RANGE = 'V', the exact value of M 150* is not known in advance and an upper bound must be used. 151* 152* LDZ (input) INTEGER 153* The leading dimension of the array Z. LDZ >= 1, and if 154* JOBZ = 'V', LDZ >= max(1,N). 155* 156* WORK (workspace) DOUBLE PRECISION array, dimension (8*N) 157* 158* IWORK (workspace) INTEGER array, dimension (5*N) 159* 160* IFAIL (output) INTEGER array, dimension (N) 161* If JOBZ = 'V', then if INFO = 0, the first M elements of 162* IFAIL are zero. If INFO > 0, then IFAIL contains the 163* indices of the eigenvectors that failed to converge. 164* If JOBZ = 'N', then IFAIL is not referenced. 165* 166* INFO (output) INTEGER 167* = 0: successful exit 168* < 0: if INFO = -i, the i-th argument had an illegal value 169* > 0: if INFO = i, then i eigenvectors failed to converge. 170* Their indices are stored in array IFAIL. 171* 172* ===================================================================== 173* 174* .. Parameters .. 175 DOUBLE PRECISION ZERO, ONE 176 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 ) 177* .. 178* .. Local Scalars .. 179 LOGICAL ALLEIG, INDEIG, VALEIG, WANTZ 180 CHARACTER ORDER 181 INTEGER I, IINFO, IMAX, INDD, INDE, INDIBL, 182 $ INDISP, INDIWO, INDTAU, INDWRK, ISCALE, ITMP1, 183 $ J, JJ, NSPLIT 184 DOUBLE PRECISION ABSTLL, ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, 185 $ SIGMA, SMLNUM, TMP1, VLL, VUU 186 DOUBLE PRECISION MXCLOCK, T1, T2 187* .. 188* .. External Functions .. 189 LOGICAL LSAME 190 DOUBLE PRECISION dlamch, dlansp 191 EXTERNAL LSAME, dlamch, dlansp 192* .. 193* .. External Subroutines .. 194 EXTERNAL dcopy, DOPGTR, DOPMTR, dscal, DSPTRD, DSTEBZ, 195 $ DSTEIN, dsteqr, dsterf, dswap, XERBLA 196* .. 197* .. Intrinsic Functions .. 198 INTRINSIC MIN, SQRT 199#include "blas_lapack.data" 200* .. 201* .. Executable Statements .. 202* 203* Test the input parameters. 204* 205 WANTZ = LSAME( JOBZ, 'V' ) 206 ALLEIG = LSAME( RANGE, 'A' ) 207 VALEIG = LSAME( RANGE, 'V' ) 208 INDEIG = LSAME( RANGE, 'I' ) 209* 210 INFO = 0 211 IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN 212 INFO = -1 213 ELSE IF( .NOT.( ALLEIG .OR. VALEIG .OR. INDEIG ) ) THEN 214 INFO = -2 215 ELSE IF( .NOT.( LSAME( UPLO, 'L' ) .OR. LSAME( UPLO, 'U' ) ) ) 216 $ THEN 217 INFO = -3 218 ELSE IF( N.LT.0 ) THEN 219 INFO = -4 220 ELSE IF( VALEIG .AND. N.GT.0 .AND. VU.LE.VL ) THEN 221 INFO = -7 222 ELSE IF( INDEIG .AND. IL.LT.1 ) THEN 223 INFO = -8 224 ELSE IF( INDEIG .AND. ( IU.LT.MIN( N, IL ) .OR. IU.GT.N ) ) THEN 225 INFO = -9 226 ELSE IF( LDZ.LT.1 .OR. ( WANTZ .AND. LDZ.LT.N ) ) THEN 227 INFO = -14 228 END IF 229* 230 IF( INFO.NE.0 ) THEN 231 CALL XERBLA( 'DSPEVX', -INFO ) 232 RETURN 233 END IF 234* 235* Quick return if possible 236* 237 M = 0 238 IF( N.EQ.0 ) 239 $ RETURN 240* 241 IF( N.EQ.1 ) THEN 242 IF( ALLEIG .OR. INDEIG ) THEN 243 M = 1 244 W( 1 ) = AP( 1 ) 245 ELSE 246 IF( VL.LT.AP( 1 ) .AND. VU.GE.AP( 1 ) ) THEN 247 M = 1 248 W( 1 ) = AP( 1 ) 249 END IF 250 END IF 251 IF( WANTZ ) 252 $ Z( 1, 1 ) = ONE 253 RETURN 254 END IF 255* 256* Get machine constants. 257* 258c SAFMIN = dlamch( 'Safe minimum' ) 259c EPS = dlamch( 'Precision' ) 260c 261 SAFMIN = DLAMCHS 262 EPS = DLAMCHP 263c 264 SMLNUM = SAFMIN / EPS 265 BIGNUM = ONE / SMLNUM 266 RMIN = SQRT( SMLNUM ) 267 RMAX = MIN( SQRT( BIGNUM ), ONE / SQRT( SQRT( SAFMIN ) ) ) 268* 269* Scale matrix to allowable range, if necessary. 270* 271 ISCALE = 0 272 ABSTLL = ABSTOL 273 IF( VALEIG ) THEN 274 VLL = VL 275 VUU = VU 276 END IF 277 ANRM = dlansp( 'M', UPLO, N, AP, WORK ) 278 IF( ANRM.GT.ZERO .AND. ANRM.LT.RMIN ) THEN 279 ISCALE = 1 280 SIGMA = RMIN / ANRM 281 ELSE IF( ANRM.GT.RMAX ) THEN 282 ISCALE = 1 283 SIGMA = RMAX / ANRM 284 END IF 285 IF( ISCALE.EQ.1 ) THEN 286 CALL dscal( ( N*( N+1 ) ) / 2, SIGMA, AP, 1 ) 287 IF( ABSTOL.GT.0 ) 288 $ ABSTLL = ABSTOL*SIGMA 289 IF( VALEIG ) THEN 290 VLL = VL*SIGMA 291 VUU = VU*SIGMA 292 END IF 293 END IF 294 295 T1 = MXCLOCK() 296 T1 = MXCLOCK() 297* 298* Call DSPTRD to reduce symmetric packed matrix to tridiagonal form. 299* 300 INDTAU = 1 301 INDE = INDTAU + N 302 INDD = INDE + N 303 INDWRK = INDD + N 304 CALL DSPTRD( UPLO, N, AP, WORK( INDD ), WORK( INDE ), 305 $ WORK( INDTAU ), IINFO ) 306* 307 T2 = MXCLOCK() 308 WRITE(*,*) ' reduce to tridiag time = ', T2-T1 309 T1 = MXCLOCK() 310 311* 312* Otherwise, call DSTEBZ and, if eigenvectors are desired, SSTEIN. 313* 314 IF( WANTZ ) THEN 315 ORDER = 'B' 316 ELSE 317 ORDER = 'E' 318 END IF 319 INDIBL = 1 320 INDISP = INDIBL + N 321 INDIWO = INDISP + N 322 CALL DSTEBZ( RANGE, ORDER, N, VLL, VUU, IL, IU, ABSTLL, 323 $ WORK( INDD ), WORK( INDE ), M, NSPLIT, W, 324 $ IWORK( INDIBL ), IWORK( INDISP ), WORK( INDWRK ), 325 $ IWORK( INDIWO ), INFO ) 326* 327 T2 = MXCLOCK() 328 WRITE(*,*) ' eigenvalue time = ', T2-T1 329 T1 = MXCLOCK() 330 331 IF( WANTZ ) THEN 332 CALL DSTEIN( N, WORK( INDD ), WORK( INDE ), M, W, 333 $ IWORK( INDIBL ), IWORK( INDISP ), Z, LDZ, 334 $ WORK( INDWRK ), IWORK( INDIWO ), IFAIL, INFO ) 335* 336 T2 = MXCLOCK() 337 WRITE(*,*) ' eigenvector time = ', T2-T1 338 T1 = MXCLOCK() 339 340* Apply orthogonal matrix used in reduction to tridiagonal 341* form to eigenvectors returned by DSTEIN. 342* 343 CALL DOPMTR( 'L', UPLO, 'N', N, M, AP, WORK( INDTAU ), Z, LDZ, 344 $ WORK( INDWRK ), INFO ) 345 T2 = MXCLOCK() 346 WRITE(*,*) ' std. back. trans. time = ', T2-T1 347 T1 = MXCLOCK() 348 349 END IF 350* 351* If matrix was scaled, then rescale eigenvalues appropriately. 352* 353 20 CONTINUE 354 IF( ISCALE.EQ.1 ) THEN 355 IF( INFO.EQ.0 ) THEN 356 IMAX = M 357 ELSE 358 IMAX = INFO - 1 359 END IF 360 CALL dscal( IMAX, ONE / SIGMA, W, 1 ) 361 END IF 362* 363* If eigenvalues are not in order, then sort them, along with 364* eigenvectors. 365* 366 IF( WANTZ ) THEN 367 DO 40 J = 1, M - 1 368 I = 0 369 TMP1 = W( J ) 370 DO 30 JJ = J + 1, M 371 IF( W( JJ ).LT.TMP1 ) THEN 372 I = JJ 373 TMP1 = W( JJ ) 374 END IF 375 30 CONTINUE 376* 377 IF( I.NE.0 ) THEN 378 ITMP1 = IWORK( INDIBL+I-1 ) 379 W( I ) = W( J ) 380 IWORK( INDIBL+I-1 ) = IWORK( INDIBL+J-1 ) 381 W( J ) = TMP1 382 IWORK( INDIBL+J-1 ) = ITMP1 383 CALL dswap( N, Z( 1, I ), 1, Z( 1, J ), 1 ) 384 IF( INFO.NE.0 ) THEN 385 ITMP1 = IFAIL( I ) 386 IFAIL( I ) = IFAIL( J ) 387 IFAIL( J ) = ITMP1 388 END IF 389 END IF 390 40 CONTINUE 391 END IF 392* 393 RETURN 394* 395* End of DSPEVX 396* 397 END 398