00001 SUBROUTINE ZHBGVX( JOBZ, RANGE, UPLO, N, KA, KB, AB, LDAB, BB, 00002 $ LDBB, Q, LDQ, VL, VU, IL, IU, ABSTOL, M, W, Z, 00003 $ LDZ, WORK, RWORK, IWORK, IFAIL, 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, KA, KB, LDAB, LDBB, LDQ, LDZ, M, 00013 $ N 00014 DOUBLE PRECISION ABSTOL, VL, VU 00015 * .. 00016 * .. Array Arguments .. 00017 INTEGER IFAIL( * ), IWORK( * ) 00018 DOUBLE PRECISION RWORK( * ), W( * ) 00019 COMPLEX*16 AB( LDAB, * ), BB( LDBB, * ), Q( LDQ, * ), 00020 $ WORK( * ), Z( LDZ, * ) 00021 * .. 00022 * 00023 * Purpose 00024 * ======= 00025 * 00026 * ZHBGVX computes all the eigenvalues, and optionally, the eigenvectors 00027 * of a complex generalized Hermitian-definite banded eigenproblem, of 00028 * the form A*x=(lambda)*B*x. Here A and B are assumed to be Hermitian 00029 * and banded, and B is also positive definite. Eigenvalues and 00030 * eigenvectors can be selected by specifying either all eigenvalues, 00031 * a range of values or a range of indices for the desired eigenvalues. 00032 * 00033 * Arguments 00034 * ========= 00035 * 00036 * JOBZ (input) CHARACTER*1 00037 * = 'N': Compute eigenvalues only; 00038 * = 'V': Compute eigenvalues and eigenvectors. 00039 * 00040 * RANGE (input) CHARACTER*1 00041 * = 'A': all eigenvalues will be found; 00042 * = 'V': all eigenvalues in the half-open interval (VL,VU] 00043 * will be found; 00044 * = 'I': the IL-th through IU-th eigenvalues will be found. 00045 * 00046 * UPLO (input) CHARACTER*1 00047 * = 'U': Upper triangles of A and B are stored; 00048 * = 'L': Lower triangles of A and B are stored. 00049 * 00050 * N (input) INTEGER 00051 * The order of the matrices A and B. N >= 0. 00052 * 00053 * KA (input) INTEGER 00054 * The number of superdiagonals of the matrix A if UPLO = 'U', 00055 * or the number of subdiagonals if UPLO = 'L'. KA >= 0. 00056 * 00057 * KB (input) INTEGER 00058 * The number of superdiagonals of the matrix B if UPLO = 'U', 00059 * or the number of subdiagonals if UPLO = 'L'. KB >= 0. 00060 * 00061 * AB (input/output) COMPLEX*16 array, dimension (LDAB, N) 00062 * On entry, the upper or lower triangle of the Hermitian band 00063 * matrix A, stored in the first ka+1 rows of the array. The 00064 * j-th column of A is stored in the j-th column of the array AB 00065 * as follows: 00066 * if UPLO = 'U', AB(ka+1+i-j,j) = A(i,j) for max(1,j-ka)<=i<=j; 00067 * if UPLO = 'L', AB(1+i-j,j) = A(i,j) for j<=i<=min(n,j+ka). 00068 * 00069 * On exit, the contents of AB are destroyed. 00070 * 00071 * LDAB (input) INTEGER 00072 * The leading dimension of the array AB. LDAB >= KA+1. 00073 * 00074 * BB (input/output) COMPLEX*16 array, dimension (LDBB, N) 00075 * On entry, the upper or lower triangle of the Hermitian band 00076 * matrix B, stored in the first kb+1 rows of the array. The 00077 * j-th column of B is stored in the j-th column of the array BB 00078 * as follows: 00079 * if UPLO = 'U', BB(kb+1+i-j,j) = B(i,j) for max(1,j-kb)<=i<=j; 00080 * if UPLO = 'L', BB(1+i-j,j) = B(i,j) for j<=i<=min(n,j+kb). 00081 * 00082 * On exit, the factor S from the split Cholesky factorization 00083 * B = S**H*S, as returned by ZPBSTF. 00084 * 00085 * LDBB (input) INTEGER 00086 * The leading dimension of the array BB. LDBB >= KB+1. 00087 * 00088 * Q (output) COMPLEX*16 array, dimension (LDQ, N) 00089 * If JOBZ = 'V', the n-by-n matrix used in the reduction of 00090 * A*x = (lambda)*B*x to standard form, i.e. C*x = (lambda)*x, 00091 * and consequently C to tridiagonal form. 00092 * If JOBZ = 'N', the array Q is not referenced. 00093 * 00094 * LDQ (input) INTEGER 00095 * The leading dimension of the array Q. If JOBZ = 'N', 00096 * LDQ >= 1. If JOBZ = 'V', LDQ >= max(1,N). 00097 * 00098 * VL (input) DOUBLE PRECISION 00099 * VU (input) DOUBLE PRECISION 00100 * If RANGE='V', the lower and upper bounds of the interval to 00101 * be searched for eigenvalues. VL < VU. 00102 * Not referenced if RANGE = 'A' or 'I'. 00103 * 00104 * IL (input) INTEGER 00105 * IU (input) INTEGER 00106 * If RANGE='I', the indices (in ascending order) of the 00107 * smallest and largest eigenvalues to be returned. 00108 * 1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0. 00109 * Not referenced if RANGE = 'A' or 'V'. 00110 * 00111 * ABSTOL (input) DOUBLE PRECISION 00112 * The absolute error tolerance for the eigenvalues. 00113 * An approximate eigenvalue is accepted as converged 00114 * when it is determined to lie in an interval [a,b] 00115 * of width less than or equal to 00116 * 00117 * ABSTOL + EPS * max( |a|,|b| ) , 00118 * 00119 * where EPS is the machine precision. If ABSTOL is less than 00120 * or equal to zero, then EPS*|T| will be used in its place, 00121 * where |T| is the 1-norm of the tridiagonal matrix obtained 00122 * by reducing AP to tridiagonal form. 00123 * 00124 * Eigenvalues will be computed most accurately when ABSTOL is 00125 * set to twice the underflow threshold 2*DLAMCH('S'), not zero. 00126 * If this routine returns with INFO>0, indicating that some 00127 * eigenvectors did not converge, try setting ABSTOL to 00128 * 2*DLAMCH('S'). 00129 * 00130 * M (output) INTEGER 00131 * The total number of eigenvalues found. 0 <= M <= N. 00132 * If RANGE = 'A', M = N, and if RANGE = 'I', M = IU-IL+1. 00133 * 00134 * W (output) DOUBLE PRECISION array, dimension (N) 00135 * If INFO = 0, the eigenvalues in ascending order. 00136 * 00137 * Z (output) COMPLEX*16 array, dimension (LDZ, N) 00138 * If JOBZ = 'V', then if INFO = 0, Z contains the matrix Z of 00139 * eigenvectors, with the i-th column of Z holding the 00140 * eigenvector associated with W(i). The eigenvectors are 00141 * normalized so that Z**H*B*Z = I. 00142 * If JOBZ = 'N', then Z is not referenced. 00143 * 00144 * LDZ (input) INTEGER 00145 * The leading dimension of the array Z. LDZ >= 1, and if 00146 * JOBZ = 'V', LDZ >= N. 00147 * 00148 * WORK (workspace) COMPLEX*16 array, dimension (N) 00149 * 00150 * RWORK (workspace) DOUBLE PRECISION array, dimension (7*N) 00151 * 00152 * IWORK (workspace) INTEGER array, dimension (5*N) 00153 * 00154 * IFAIL (output) INTEGER array, dimension (N) 00155 * If JOBZ = 'V', then if INFO = 0, the first M elements of 00156 * IFAIL are zero. If INFO > 0, then IFAIL contains the 00157 * indices of the eigenvectors that failed to converge. 00158 * If JOBZ = 'N', then IFAIL is not referenced. 00159 * 00160 * INFO (output) INTEGER 00161 * = 0: successful exit 00162 * < 0: if INFO = -i, the i-th argument had an illegal value 00163 * > 0: if INFO = i, and i is: 00164 * <= N: then i eigenvectors failed to converge. Their 00165 * indices are stored in array IFAIL. 00166 * > N: if INFO = N + i, for 1 <= i <= N, then ZPBSTF 00167 * returned INFO = i: B is not positive definite. 00168 * The factorization of B could not be completed and 00169 * no eigenvalues or eigenvectors were computed. 00170 * 00171 * Further Details 00172 * =============== 00173 * 00174 * Based on contributions by 00175 * Mark Fahey, Department of Mathematics, Univ. of Kentucky, USA 00176 * 00177 * ===================================================================== 00178 * 00179 * .. Parameters .. 00180 DOUBLE PRECISION ZERO 00181 PARAMETER ( ZERO = 0.0D+0 ) 00182 COMPLEX*16 CZERO, CONE 00183 PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ), 00184 $ CONE = ( 1.0D+0, 0.0D+0 ) ) 00185 * .. 00186 * .. Local Scalars .. 00187 LOGICAL ALLEIG, INDEIG, TEST, UPPER, VALEIG, WANTZ 00188 CHARACTER ORDER, VECT 00189 INTEGER I, IINFO, INDD, INDE, INDEE, INDIBL, INDISP, 00190 $ INDIWK, INDRWK, INDWRK, ITMP1, J, JJ, NSPLIT 00191 DOUBLE PRECISION TMP1 00192 * .. 00193 * .. External Functions .. 00194 LOGICAL LSAME 00195 EXTERNAL LSAME 00196 * .. 00197 * .. External Subroutines .. 00198 EXTERNAL DCOPY, DSTEBZ, DSTERF, XERBLA, ZCOPY, ZGEMV, 00199 $ ZHBGST, ZHBTRD, ZLACPY, ZPBSTF, ZSTEIN, ZSTEQR, 00200 $ ZSWAP 00201 * .. 00202 * .. Intrinsic Functions .. 00203 INTRINSIC MIN 00204 * .. 00205 * .. Executable Statements .. 00206 * 00207 * Test the input parameters. 00208 * 00209 WANTZ = LSAME( JOBZ, 'V' ) 00210 UPPER = LSAME( UPLO, 'U' ) 00211 ALLEIG = LSAME( RANGE, 'A' ) 00212 VALEIG = LSAME( RANGE, 'V' ) 00213 INDEIG = LSAME( RANGE, 'I' ) 00214 * 00215 INFO = 0 00216 IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN 00217 INFO = -1 00218 ELSE IF( .NOT.( ALLEIG .OR. VALEIG .OR. INDEIG ) ) THEN 00219 INFO = -2 00220 ELSE IF( .NOT.( UPPER .OR. LSAME( UPLO, 'L' ) ) ) THEN 00221 INFO = -3 00222 ELSE IF( N.LT.0 ) THEN 00223 INFO = -4 00224 ELSE IF( KA.LT.0 ) THEN 00225 INFO = -5 00226 ELSE IF( KB.LT.0 .OR. KB.GT.KA ) THEN 00227 INFO = -6 00228 ELSE IF( LDAB.LT.KA+1 ) THEN 00229 INFO = -8 00230 ELSE IF( LDBB.LT.KB+1 ) THEN 00231 INFO = -10 00232 ELSE IF( LDQ.LT.1 .OR. ( WANTZ .AND. LDQ.LT.N ) ) THEN 00233 INFO = -12 00234 ELSE 00235 IF( VALEIG ) THEN 00236 IF( N.GT.0 .AND. VU.LE.VL ) 00237 $ INFO = -14 00238 ELSE IF( INDEIG ) THEN 00239 IF( IL.LT.1 .OR. IL.GT.MAX( 1, N ) ) THEN 00240 INFO = -15 00241 ELSE IF ( IU.LT.MIN( N, IL ) .OR. IU.GT.N ) THEN 00242 INFO = -16 00243 END IF 00244 END IF 00245 END IF 00246 IF( INFO.EQ.0) THEN 00247 IF( LDZ.LT.1 .OR. ( WANTZ .AND. LDZ.LT.N ) ) THEN 00248 INFO = -21 00249 END IF 00250 END IF 00251 * 00252 IF( INFO.NE.0 ) THEN 00253 CALL XERBLA( 'ZHBGVX', -INFO ) 00254 RETURN 00255 END IF 00256 * 00257 * Quick return if possible 00258 * 00259 M = 0 00260 IF( N.EQ.0 ) 00261 $ RETURN 00262 * 00263 * Form a split Cholesky factorization of B. 00264 * 00265 CALL ZPBSTF( UPLO, N, KB, BB, LDBB, INFO ) 00266 IF( INFO.NE.0 ) THEN 00267 INFO = N + INFO 00268 RETURN 00269 END IF 00270 * 00271 * Transform problem to standard eigenvalue problem. 00272 * 00273 CALL ZHBGST( JOBZ, UPLO, N, KA, KB, AB, LDAB, BB, LDBB, Q, LDQ, 00274 $ WORK, RWORK, IINFO ) 00275 * 00276 * Solve the standard eigenvalue problem. 00277 * Reduce Hermitian band matrix to tridiagonal form. 00278 * 00279 INDD = 1 00280 INDE = INDD + N 00281 INDRWK = INDE + N 00282 INDWRK = 1 00283 IF( WANTZ ) THEN 00284 VECT = 'U' 00285 ELSE 00286 VECT = 'N' 00287 END IF 00288 CALL ZHBTRD( VECT, UPLO, N, KA, AB, LDAB, RWORK( INDD ), 00289 $ RWORK( INDE ), Q, LDQ, WORK( INDWRK ), IINFO ) 00290 * 00291 * If all eigenvalues are desired and ABSTOL is less than or equal 00292 * to zero, then call DSTERF or ZSTEQR. If this fails for some 00293 * eigenvalue, then try DSTEBZ. 00294 * 00295 TEST = .FALSE. 00296 IF( INDEIG ) THEN 00297 IF( IL.EQ.1 .AND. IU.EQ.N ) THEN 00298 TEST = .TRUE. 00299 END IF 00300 END IF 00301 IF( ( ALLEIG .OR. TEST ) .AND. ( ABSTOL.LE.ZERO ) ) THEN 00302 CALL DCOPY( N, RWORK( INDD ), 1, W, 1 ) 00303 INDEE = INDRWK + 2*N 00304 CALL DCOPY( N-1, RWORK( INDE ), 1, RWORK( INDEE ), 1 ) 00305 IF( .NOT.WANTZ ) THEN 00306 CALL DSTERF( N, W, RWORK( INDEE ), INFO ) 00307 ELSE 00308 CALL ZLACPY( 'A', N, N, Q, LDQ, Z, LDZ ) 00309 CALL ZSTEQR( JOBZ, N, W, RWORK( INDEE ), Z, LDZ, 00310 $ RWORK( INDRWK ), INFO ) 00311 IF( INFO.EQ.0 ) THEN 00312 DO 10 I = 1, N 00313 IFAIL( I ) = 0 00314 10 CONTINUE 00315 END IF 00316 END IF 00317 IF( INFO.EQ.0 ) THEN 00318 M = N 00319 GO TO 30 00320 END IF 00321 INFO = 0 00322 END IF 00323 * 00324 * Otherwise, call DSTEBZ and, if eigenvectors are desired, 00325 * call ZSTEIN. 00326 * 00327 IF( WANTZ ) THEN 00328 ORDER = 'B' 00329 ELSE 00330 ORDER = 'E' 00331 END IF 00332 INDIBL = 1 00333 INDISP = INDIBL + N 00334 INDIWK = INDISP + N 00335 CALL DSTEBZ( RANGE, ORDER, N, VL, VU, IL, IU, ABSTOL, 00336 $ RWORK( INDD ), RWORK( INDE ), M, NSPLIT, W, 00337 $ IWORK( INDIBL ), IWORK( INDISP ), RWORK( INDRWK ), 00338 $ IWORK( INDIWK ), INFO ) 00339 * 00340 IF( WANTZ ) THEN 00341 CALL ZSTEIN( N, RWORK( INDD ), RWORK( INDE ), M, W, 00342 $ IWORK( INDIBL ), IWORK( INDISP ), Z, LDZ, 00343 $ RWORK( INDRWK ), IWORK( INDIWK ), IFAIL, INFO ) 00344 * 00345 * Apply unitary matrix used in reduction to tridiagonal 00346 * form to eigenvectors returned by ZSTEIN. 00347 * 00348 DO 20 J = 1, M 00349 CALL ZCOPY( N, Z( 1, J ), 1, WORK( 1 ), 1 ) 00350 CALL ZGEMV( 'N', N, N, CONE, Q, LDQ, WORK, 1, CZERO, 00351 $ Z( 1, J ), 1 ) 00352 20 CONTINUE 00353 END IF 00354 * 00355 30 CONTINUE 00356 * 00357 * If eigenvalues are not in order, then sort them, along with 00358 * eigenvectors. 00359 * 00360 IF( WANTZ ) THEN 00361 DO 50 J = 1, M - 1 00362 I = 0 00363 TMP1 = W( J ) 00364 DO 40 JJ = J + 1, M 00365 IF( W( JJ ).LT.TMP1 ) THEN 00366 I = JJ 00367 TMP1 = W( JJ ) 00368 END IF 00369 40 CONTINUE 00370 * 00371 IF( I.NE.0 ) THEN 00372 ITMP1 = IWORK( INDIBL+I-1 ) 00373 W( I ) = W( J ) 00374 IWORK( INDIBL+I-1 ) = IWORK( INDIBL+J-1 ) 00375 W( J ) = TMP1 00376 IWORK( INDIBL+J-1 ) = ITMP1 00377 CALL ZSWAP( N, Z( 1, I ), 1, Z( 1, J ), 1 ) 00378 IF( INFO.NE.0 ) THEN 00379 ITMP1 = IFAIL( I ) 00380 IFAIL( I ) = IFAIL( J ) 00381 IFAIL( J ) = ITMP1 00382 END IF 00383 END IF 00384 50 CONTINUE 00385 END IF 00386 * 00387 RETURN 00388 * 00389 * End of ZHBGVX 00390 * 00391 END
1.7.2