1 SUBROUTINE DSBGVX( JOBZ, RANGE, UPLO, N, KA, KB, AB, LDAB, BB, 2 $ LDBB, Q, LDQ, VL, VU, IL, IU, ABSTOL, M, W, Z, 3 $ LDZ, WORK, IWORK, IFAIL, INFO ) 4* 5* -- LAPACK driver routine (version 3.0) -- 6* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., 7* Courant Institute, Argonne National Lab, and Rice University 8* June 30, 1999 9* 10* .. Scalar Arguments .. 11 CHARACTER JOBZ, RANGE, UPLO 12 INTEGER IL, INFO, IU, KA, KB, LDAB, LDBB, LDQ, LDZ, M, 13 $ N 14 DOUBLE PRECISION ABSTOL, VL, VU 15* .. 16* .. Array Arguments .. 17 INTEGER IFAIL( * ), IWORK( * ) 18 DOUBLE PRECISION AB( LDAB, * ), BB( LDBB, * ), Q( LDQ, * ), 19 $ W( * ), WORK( * ), Z( LDZ, * ) 20* .. 21* 22* Purpose 23* ======= 24* 25* DSBGVX computes selected eigenvalues, and optionally, eigenvectors 26* of a real generalized symmetric-definite banded eigenproblem, of 27* the form A*x=(lambda)*B*x. Here A and B are assumed to be symmetric 28* and banded, and B is also positive definite. Eigenvalues and 29* eigenvectors can be selected by specifying either all eigenvalues, 30* a range of values or a range of indices for the desired eigenvalues. 31* 32* Arguments 33* ========= 34* 35* JOBZ (input) CHARACTER*1 36* = 'N': Compute eigenvalues only; 37* = 'V': Compute eigenvalues and eigenvectors. 38* 39* RANGE (input) CHARACTER*1 40* = 'A': all eigenvalues will be found. 41* = 'V': all eigenvalues in the half-open interval (VL,VU] 42* will be found. 43* = 'I': the IL-th through IU-th eigenvalues will be found. 44* 45* UPLO (input) CHARACTER*1 46* = 'U': Upper triangles of A and B are stored; 47* = 'L': Lower triangles of A and B are stored. 48* 49* N (input) INTEGER 50* The order of the matrices A and B. N >= 0. 51* 52* KA (input) INTEGER 53* The number of superdiagonals of the matrix A if UPLO = 'U', 54* or the number of subdiagonals if UPLO = 'L'. KA >= 0. 55* 56* KB (input) INTEGER 57* The number of superdiagonals of the matrix B if UPLO = 'U', 58* or the number of subdiagonals if UPLO = 'L'. KB >= 0. 59* 60* AB (input/output) DOUBLE PRECISION array, dimension (LDAB, N) 61* On entry, the upper or lower triangle of the symmetric band 62* matrix A, stored in the first ka+1 rows of the array. The 63* j-th column of A is stored in the j-th column of the array AB 64* as follows: 65* if UPLO = 'U', AB(ka+1+i-j,j) = A(i,j) for max(1,j-ka)<=i<=j; 66* if UPLO = 'L', AB(1+i-j,j) = A(i,j) for j<=i<=min(n,j+ka). 67* 68* On exit, the contents of AB are destroyed. 69* 70* LDAB (input) INTEGER 71* The leading dimension of the array AB. LDAB >= KA+1. 72* 73* BB (input/output) DOUBLE PRECISION array, dimension (LDBB, N) 74* On entry, the upper or lower triangle of the symmetric band 75* matrix B, stored in the first kb+1 rows of the array. The 76* j-th column of B is stored in the j-th column of the array BB 77* as follows: 78* if UPLO = 'U', BB(ka+1+i-j,j) = B(i,j) for max(1,j-kb)<=i<=j; 79* if UPLO = 'L', BB(1+i-j,j) = B(i,j) for j<=i<=min(n,j+kb). 80* 81* On exit, the factor S from the split Cholesky factorization 82* B = S**T*S, as returned by DPBSTF. 83* 84* LDBB (input) INTEGER 85* The leading dimension of the array BB. LDBB >= KB+1. 86* 87* Q (output) DOUBLE PRECISION array, dimension (LDQ, N) 88* If JOBZ = 'V', the n-by-n matrix used in the reduction of 89* A*x = (lambda)*B*x to standard form, i.e. C*x = (lambda)*x, 90* and consequently C to tridiagonal form. 91* If JOBZ = 'N', the array Q is not referenced. 92* 93* LDQ (input) INTEGER 94* The leading dimension of the array Q. If JOBZ = 'N', 95* LDQ >= 1. If JOBZ = 'V', LDQ >= max(1,N). 96* 97* VL (input) DOUBLE PRECISION 98* VU (input) DOUBLE PRECISION 99* If RANGE='V', the lower and upper bounds of the interval to 100* be searched for eigenvalues. VL < VU. 101* Not referenced if RANGE = 'A' or 'I'. 102* 103* IL (input) INTEGER 104* IU (input) INTEGER 105* If RANGE='I', the indices (in ascending order) of the 106* smallest and largest eigenvalues to be returned. 107* 1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0. 108* Not referenced if RANGE = 'A' or 'V'. 109* 110* ABSTOL (input) DOUBLE PRECISION 111* The absolute error tolerance for the eigenvalues. 112* An approximate eigenvalue is accepted as converged 113* when it is determined to lie in an interval [a,b] 114* of width less than or equal to 115* 116* ABSTOL + EPS * max( |a|,|b| ) , 117* 118* where EPS is the machine precision. If ABSTOL is less than 119* or equal to zero, then EPS*|T| will be used in its place, 120* where |T| is the 1-norm of the tridiagonal matrix obtained 121* by reducing A to tridiagonal form. 122* 123* Eigenvalues will be computed most accurately when ABSTOL is 124* set to twice the underflow threshold 2*DLAMCH('S'), not zero. 125* If this routine returns with INFO>0, indicating that some 126* eigenvectors did not converge, try setting ABSTOL to 127* 2*DLAMCH('S'). 128* 129* M (output) INTEGER 130* The total number of eigenvalues found. 0 <= M <= N. 131* If RANGE = 'A', M = N, and if RANGE = 'I', M = IU-IL+1. 132* 133* W (output) DOUBLE PRECISION array, dimension (N) 134* If INFO = 0, the eigenvalues in ascending order. 135* 136* Z (output) DOUBLE PRECISION array, dimension (LDZ, N) 137* If JOBZ = 'V', then if INFO = 0, Z contains the matrix Z of 138* eigenvectors, with the i-th column of Z holding the 139* eigenvector associated with W(i). The eigenvectors are 140* normalized so Z**T*B*Z = I. 141* If JOBZ = 'N', then Z is not referenced. 142* 143* LDZ (input) INTEGER 144* The leading dimension of the array Z. LDZ >= 1, and if 145* JOBZ = 'V', LDZ >= max(1,N). 146* 147* WORK (workspace/output) DOUBLE PRECISION array, dimension (7N) 148* 149* IWORK (workspace/output) INTEGER array, dimension (5N) 150* 151* IFAIL (input) INTEGER array, dimension (M) 152* If JOBZ = 'V', then if INFO = 0, the first M elements of 153* IFAIL are zero. If INFO > 0, then IFAIL contains the 154* indices of the eigenvalues that failed to converge. 155* If JOBZ = 'N', then IFAIL is not referenced. 156* 157* INFO (output) INTEGER 158* = 0 : successful exit 159* < 0 : if INFO = -i, the i-th argument had an illegal value 160* <= N: if INFO = i, then i eigenvectors failed to converge. 161* Their indices are stored in IFAIL. 162* > N : DPBSTF returned an error code; i.e., 163* if INFO = N + i, for 1 <= i <= N, then the leading 164* minor of order i of B is not positive definite. 165* The factorization of B could not be completed and 166* no eigenvalues or eigenvectors were computed. 167* 168* Further Details 169* =============== 170* 171* Based on contributions by 172* Mark Fahey, Department of Mathematics, Univ. of Kentucky, USA 173* 174* ===================================================================== 175* 176* .. Parameters .. 177 DOUBLE PRECISION ZERO, ONE 178 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) 179* .. 180* .. Local Scalars .. 181 LOGICAL ALLEIG, INDEIG, UPPER, VALEIG, WANTZ 182 CHARACTER ORDER, VECT 183 INTEGER I, IINFO, INDD, INDE, INDEE, INDIBL, INDISP, 184 $ INDIWO, INDWRK, ITMP1, J, JJ, NSPLIT 185 DOUBLE PRECISION TMP1 186* .. 187* .. External Functions .. 188 LOGICAL LSAME 189 EXTERNAL LSAME 190* .. 191* .. External Subroutines .. 192 EXTERNAL DCOPY, DGEMV, DLACPY, DPBSTF, DSBGST, DSBTRD, 193 $ DSTEBZ, DSTEIN, DSTEQR, DSTERF, DSWAP, XERBLA 194* .. 195* .. Intrinsic Functions .. 196 INTRINSIC MIN 197* .. 198* .. Executable Statements .. 199* 200* Test the input parameters. 201* 202 WANTZ = LSAME( JOBZ, 'V' ) 203 UPPER = LSAME( UPLO, 'U' ) 204 ALLEIG = LSAME( RANGE, 'A' ) 205 VALEIG = LSAME( RANGE, 'V' ) 206 INDEIG = LSAME( RANGE, 'I' ) 207* 208 INFO = 0 209 IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN 210 INFO = -1 211 ELSE IF( .NOT.( ALLEIG .OR. VALEIG .OR. INDEIG ) ) THEN 212 INFO = -2 213 ELSE IF( .NOT.( UPPER .OR. LSAME( UPLO, 'L' ) ) ) THEN 214 INFO = -3 215 ELSE IF( N.LT.0 ) THEN 216 INFO = -4 217 ELSE IF( KA.LT.0 ) THEN 218 INFO = -5 219 ELSE IF( KB.LT.0 .OR. KB.GT.KA ) THEN 220 INFO = -6 221 ELSE IF( LDAB.LT.KA+1 ) THEN 222 INFO = -8 223 ELSE IF( LDBB.LT.KB+1 ) THEN 224 INFO = -10 225 ELSE IF( LDQ.LT.1 ) THEN 226 INFO = -12 227 ELSE IF( VALEIG .AND. N.GT.0 .AND. VU.LE.VL ) THEN 228 INFO = -14 229 ELSE IF( INDEIG .AND. IL.LT.1 ) THEN 230 INFO = -15 231 ELSE IF( INDEIG .AND. ( IU.LT.MIN( N, IL ) .OR. IU.GT.N ) ) THEN 232 INFO = -16 233 ELSE IF( LDZ.LT.1 .OR. ( WANTZ .AND. LDZ.LT.N ) ) THEN 234 INFO = -21 235 END IF 236* 237 IF( INFO.NE.0 ) THEN 238 CALL XERBLA( 'DSBGVX', -INFO ) 239 RETURN 240 END IF 241* 242* Quick return if possible 243* 244 M = 0 245 IF( N.EQ.0 ) THEN 246 WORK( 1 ) = 1 247 RETURN 248 END IF 249* 250* Form a split Cholesky factorization of B. 251* 252 CALL DPBSTF( UPLO, N, KB, BB, LDBB, INFO ) 253 IF( INFO.NE.0 ) THEN 254 INFO = N + INFO 255 RETURN 256 END IF 257* 258* Transform problem to standard eigenvalue problem. 259* 260 CALL DSBGST( JOBZ, UPLO, N, KA, KB, AB, LDAB, BB, LDBB, Q, LDQ, 261 $ WORK, IINFO ) 262* 263* Reduce symmetric band matrix to tridiagonal form. 264* 265 INDD = 1 266 INDE = INDD + N 267 INDWRK = INDE + N 268 IF( WANTZ ) THEN 269 VECT = 'U' 270 ELSE 271 VECT = 'N' 272 END IF 273 CALL DSBTRD( VECT, UPLO, N, KA, AB, LDAB, WORK( INDD ), 274 $ WORK( INDE ), Q, LDQ, WORK( INDWRK ), IINFO ) 275* 276* If all eigenvalues are desired and ABSTOL is less than or equal 277* to zero, then call DSTERF or SSTEQR. If this fails for some 278* eigenvalue, then try DSTEBZ. 279* 280 IF( ( ALLEIG .OR. ( INDEIG .AND. IL.EQ.1 .AND. IU.EQ.N ) ) .AND. 281 $ ( ABSTOL.LE.ZERO ) ) THEN 282 CALL DCOPY( N, WORK( INDD ), 1, W, 1 ) 283 INDEE = INDWRK + 2*N 284 CALL DCOPY( N-1, WORK( INDE ), 1, WORK( INDEE ), 1 ) 285 IF( .NOT.WANTZ ) THEN 286 CALL DSTERF( N, W, WORK( INDEE ), INFO ) 287 ELSE 288 CALL DLACPY( 'A', N, N, Q, LDQ, Z, LDZ ) 289 CALL DSTEQR( JOBZ, N, W, WORK( INDEE ), Z, LDZ, 290 $ WORK( INDWRK ), INFO ) 291 IF( INFO.EQ.0 ) THEN 292 DO 10 I = 1, N 293 IFAIL( I ) = 0 294 10 CONTINUE 295 END IF 296 END IF 297 IF( INFO.EQ.0 ) THEN 298 M = N 299 GO TO 30 300 END IF 301 INFO = 0 302 END IF 303* 304* Otherwise, call DSTEBZ and, if eigenvectors are desired, 305* call DSTEIN. 306* 307 IF( WANTZ ) THEN 308 ORDER = 'B' 309 ELSE 310 ORDER = 'E' 311 END IF 312 INDIBL = 1 313 INDISP = INDIBL + N 314 INDIWO = INDISP + N 315 CALL DSTEBZ( RANGE, ORDER, N, VL, VU, IL, IU, ABSTOL, 316 $ WORK( INDD ), WORK( INDE ), M, NSPLIT, W, 317 $ IWORK( INDIBL ), IWORK( INDISP ), WORK( INDWRK ), 318 $ IWORK( INDIWO ), INFO ) 319* 320 IF( WANTZ ) THEN 321 CALL DSTEIN( N, WORK( INDD ), WORK( INDE ), M, W, 322 $ IWORK( INDIBL ), IWORK( INDISP ), Z, LDZ, 323 $ WORK( INDWRK ), IWORK( INDIWO ), IFAIL, INFO ) 324* 325* Apply transformation matrix used in reduction to tridiagonal 326* form to eigenvectors returned by DSTEIN. 327* 328 DO 20 J = 1, M 329 CALL DCOPY( N, Z( 1, J ), 1, WORK( 1 ), 1 ) 330 CALL DGEMV( 'N', N, N, ONE, Q, LDQ, WORK, 1, ZERO, 331 $ Z( 1, J ), 1 ) 332 20 CONTINUE 333 END IF 334* 335 30 CONTINUE 336* 337* If eigenvalues are not in order, then sort them, along with 338* eigenvectors. 339* 340 IF( WANTZ ) THEN 341 DO 50 J = 1, M - 1 342 I = 0 343 TMP1 = W( J ) 344 DO 40 JJ = J + 1, M 345 IF( W( JJ ).LT.TMP1 ) THEN 346 I = JJ 347 TMP1 = W( JJ ) 348 END IF 349 40 CONTINUE 350* 351 IF( I.NE.0 ) THEN 352 ITMP1 = IWORK( INDIBL+I-1 ) 353 W( I ) = W( J ) 354 IWORK( INDIBL+I-1 ) = IWORK( INDIBL+J-1 ) 355 W( J ) = TMP1 356 IWORK( INDIBL+J-1 ) = ITMP1 357 CALL DSWAP( N, Z( 1, I ), 1, Z( 1, J ), 1 ) 358 IF( INFO.NE.0 ) THEN 359 ITMP1 = IFAIL( I ) 360 IFAIL( I ) = IFAIL( J ) 361 IFAIL( J ) = ITMP1 362 END IF 363 END IF 364 50 CONTINUE 365 END IF 366* 367 RETURN 368* 369* End of DSBGVX 370* 371 END 372