1 SUBROUTINE CHBGVD( JOBZ, UPLO, N, KA, KB, AB, LDAB, BB, LDBB, W, 2 $ Z, LDZ, WORK, LWORK, RWORK, LRWORK, IWORK, 3 $ LIWORK, 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* October 31, 1999 9* 10* .. Scalar Arguments .. 11 CHARACTER JOBZ, UPLO 12 INTEGER INFO, KA, KB, LDAB, LDBB, LDZ, LIWORK, LRWORK, 13 $ LWORK, N 14* .. 15* .. Array Arguments .. 16 INTEGER IWORK( * ) 17 REAL RWORK( * ), W( * ) 18 COMPLEX AB( LDAB, * ), BB( LDBB, * ), WORK( * ), 19 $ Z( LDZ, * ) 20* .. 21* 22* Purpose 23* ======= 24* 25* CHBGVD computes all the eigenvalues, and optionally, the eigenvectors 26* of a complex generalized Hermitian-definite banded eigenproblem, of 27* the form A*x=(lambda)*B*x. Here A and B are assumed to be Hermitian 28* and banded, and B is also positive definite. If eigenvectors are 29* desired, it uses a divide and conquer algorithm. 30* 31* The divide and conquer algorithm makes very mild assumptions about 32* floating point arithmetic. It will work on machines with a guard 33* digit in add/subtract, or on those binary machines without guard 34* digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or 35* Cray-2. It could conceivably fail on hexadecimal or decimal machines 36* without guard digits, but we know of none. 37* 38* Arguments 39* ========= 40* 41* JOBZ (input) CHARACTER*1 42* = 'N': Compute eigenvalues only; 43* = 'V': Compute eigenvalues and eigenvectors. 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) COMPLEX array, dimension (LDAB, N) 61* On entry, the upper or lower triangle of the Hermitian 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) COMPLEX array, dimension (LDBB, N) 74* On entry, the upper or lower triangle of the Hermitian 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(kb+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**H*S, as returned by CPBSTF. 83* 84* LDBB (input) INTEGER 85* The leading dimension of the array BB. LDBB >= KB+1. 86* 87* W (output) REAL array, dimension (N) 88* If INFO = 0, the eigenvalues in ascending order. 89* 90* Z (output) COMPLEX array, dimension (LDZ, N) 91* If JOBZ = 'V', then if INFO = 0, Z contains the matrix Z of 92* eigenvectors, with the i-th column of Z holding the 93* eigenvector associated with W(i). The eigenvectors are 94* normalized so that Z**H*B*Z = I. 95* If JOBZ = 'N', then Z is not referenced. 96* 97* LDZ (input) INTEGER 98* The leading dimension of the array Z. LDZ >= 1, and if 99* JOBZ = 'V', LDZ >= N. 100* 101* WORK (workspace/output) COMPLEX array, dimension (LWORK) 102* On exit, if INFO=0, WORK(1) returns the optimal LWORK. 103* 104* LWORK (input) INTEGER 105* The dimension of the array WORK. 106* If N <= 1, LWORK >= 1. 107* If JOBZ = 'N' and N > 1, LWORK >= N. 108* If JOBZ = 'V' and N > 1, LWORK >= 2*N**2. 109* 110* If LWORK = -1, then a workspace query is assumed; the routine 111* only calculates the optimal size of the WORK array, returns 112* this value as the first entry of the WORK array, and no error 113* message related to LWORK is issued by XERBLA. 114* 115* RWORK (workspace/output) REAL array, dimension (LRWORK) 116* On exit, if INFO=0, RWORK(1) returns the optimal LRWORK. 117* 118* LRWORK (input) INTEGER 119* The dimension of array RWORK. 120* If N <= 1, LRWORK >= 1. 121* If JOBZ = 'N' and N > 1, LRWORK >= N. 122* If JOBZ = 'V' and N > 1, LRWORK >= 1 + 5*N + 2*N**2. 123* 124* If LRWORK = -1, then a workspace query is assumed; the 125* routine only calculates the optimal size of the RWORK array, 126* returns this value as the first entry of the RWORK array, and 127* no error message related to LRWORK is issued by XERBLA. 128* 129* IWORK (workspace/output) INTEGER array, dimension (LIWORK) 130* On exit, if INFO=0, IWORK(1) returns the optimal LIWORK. 131* 132* LIWORK (input) INTEGER 133* The dimension of array IWORK. 134* If JOBZ = 'N' or N <= 1, LIWORK >= 1. 135* If JOBZ = 'V' and N > 1, LIWORK >= 3 + 5*N. 136* 137* If LIWORK = -1, then a workspace query is assumed; the 138* routine only calculates the optimal size of the IWORK array, 139* returns this value as the first entry of the IWORK array, and 140* no error message related to LIWORK is issued by XERBLA. 141* 142* INFO (output) INTEGER 143* = 0: successful exit 144* < 0: if INFO = -i, the i-th argument had an illegal value 145* > 0: if INFO = i, and i is: 146* <= N: the algorithm failed to converge: 147* i off-diagonal elements of an intermediate 148* tridiagonal form did not converge to zero; 149* > N: if INFO = N + i, for 1 <= i <= N, then CPBSTF 150* returned INFO = i: B is not positive definite. 151* The factorization of B could not be completed and 152* no eigenvalues or eigenvectors were computed. 153* 154* Further Details 155* =============== 156* 157* Based on contributions by 158* Mark Fahey, Department of Mathematics, Univ. of Kentucky, USA 159* 160* ===================================================================== 161* 162* .. Parameters .. 163 COMPLEX CONE, CZERO 164 PARAMETER ( CONE = ( 1.0E+0, 0.0E+0 ), 165 $ CZERO = ( 0.0E+0, 0.0E+0 ) ) 166* .. 167* .. Local Scalars .. 168 LOGICAL LQUERY, UPPER, WANTZ 169 CHARACTER VECT 170 INTEGER IINFO, INDE, INDWK2, INDWRK, LIWMIN, LLRWK, 171 $ LLWK2, LRWMIN, LWMIN 172* .. 173* .. External Functions .. 174 LOGICAL LSAME 175 EXTERNAL LSAME 176* .. 177* .. External Subroutines .. 178 EXTERNAL CGEMM, CHBGST, CHBTRD, CLACPY, CPBSTF, CSTEDC, 179 $ SSTERF, XERBLA 180* .. 181* .. Executable Statements .. 182* 183* Test the input parameters. 184* 185 WANTZ = LSAME( JOBZ, 'V' ) 186 UPPER = LSAME( UPLO, 'U' ) 187 LQUERY = ( LWORK.EQ.-1 .OR. LRWORK.EQ.-1 .OR. LIWORK.EQ.-1 ) 188* 189 INFO = 0 190 IF( N.LE.1 ) THEN 191 LWMIN = 1 192 LRWMIN = 1 193 LIWMIN = 1 194 ELSE 195 IF( WANTZ ) THEN 196 LWMIN = 2*N**2 197 LRWMIN = 1 + 5*N + 2*N**2 198 LIWMIN = 3 + 5*N 199 ELSE 200 LWMIN = N 201 LRWMIN = N 202 LIWMIN = 1 203 END IF 204 END IF 205 IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN 206 INFO = -1 207 ELSE IF( .NOT.( UPPER .OR. LSAME( UPLO, 'L' ) ) ) THEN 208 INFO = -2 209 ELSE IF( N.LT.0 ) THEN 210 INFO = -3 211 ELSE IF( KA.LT.0 ) THEN 212 INFO = -4 213 ELSE IF( KB.LT.0 .OR. KB.GT.KA ) THEN 214 INFO = -5 215 ELSE IF( LDAB.LT.KA+1 ) THEN 216 INFO = -7 217 ELSE IF( LDBB.LT.KB+1 ) THEN 218 INFO = -9 219 ELSE IF( LDZ.LT.1 .OR. ( WANTZ .AND. LDZ.LT.N ) ) THEN 220 INFO = -12 221 ELSE IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN 222 INFO = -14 223 ELSE IF( LRWORK.LT.LRWMIN .AND. .NOT.LQUERY ) THEN 224 INFO = -16 225 ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT.LQUERY ) THEN 226 INFO = -18 227 END IF 228* 229 IF( INFO.EQ.0 ) THEN 230 WORK( 1 ) = LWMIN 231 RWORK( 1 ) = LRWMIN 232 IWORK( 1 ) = LIWMIN 233 END IF 234* 235 IF( INFO.NE.0 ) THEN 236 CALL XERBLA( 'CHBGVD', -INFO ) 237 RETURN 238 ELSE IF( LQUERY ) THEN 239 RETURN 240 END IF 241* 242* Quick return if possible 243* 244 IF( N.EQ.0 ) 245 $ RETURN 246* 247* Form a split Cholesky factorization of B. 248* 249 CALL CPBSTF( UPLO, N, KB, BB, LDBB, INFO ) 250 IF( INFO.NE.0 ) THEN 251 INFO = N + INFO 252 RETURN 253 END IF 254* 255* Transform problem to standard eigenvalue problem. 256* 257 INDE = 1 258 INDWRK = INDE + N 259 INDWK2 = 1 + N*N 260 LLWK2 = LWORK - INDWK2 + 2 261 LLRWK = LRWORK - INDWRK + 2 262 CALL CHBGST( JOBZ, UPLO, N, KA, KB, AB, LDAB, BB, LDBB, Z, LDZ, 263 $ WORK, RWORK( INDWRK ), IINFO ) 264* 265* Reduce Hermitian band matrix to tridiagonal form. 266* 267 IF( WANTZ ) THEN 268 VECT = 'U' 269 ELSE 270 VECT = 'N' 271 END IF 272 CALL CHBTRD( VECT, UPLO, N, KA, AB, LDAB, W, RWORK( INDE ), Z, 273 $ LDZ, WORK, IINFO ) 274* 275* For eigenvalues only, call SSTERF. For eigenvectors, call CSTEDC. 276* 277 IF( .NOT.WANTZ ) THEN 278 CALL SSTERF( N, W, RWORK( INDE ), INFO ) 279 ELSE 280 CALL CSTEDC( 'I', N, W, RWORK( INDE ), WORK, N, WORK( INDWK2 ), 281 $ LLWK2, RWORK( INDWRK ), LLRWK, IWORK, LIWORK, 282 $ INFO ) 283 CALL CGEMM( 'N', 'N', N, N, N, CONE, Z, LDZ, WORK, N, CZERO, 284 $ WORK( INDWK2 ), N ) 285 CALL CLACPY( 'A', N, N, WORK( INDWK2 ), N, Z, LDZ ) 286 END IF 287* 288 WORK( 1 ) = LWMIN 289 RWORK( 1 ) = LRWMIN 290 IWORK( 1 ) = LIWMIN 291 RETURN 292* 293* End of CHBGVD 294* 295 END 296