00001 SUBROUTINE CHBGVD( JOBZ, UPLO, N, KA, KB, AB, LDAB, BB, LDBB, W, 00002 $ Z, LDZ, WORK, LWORK, RWORK, LRWORK, IWORK, 00003 $ LIWORK, 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, UPLO 00012 INTEGER INFO, KA, KB, LDAB, LDBB, LDZ, LIWORK, LRWORK, 00013 $ LWORK, N 00014 * .. 00015 * .. Array Arguments .. 00016 INTEGER IWORK( * ) 00017 REAL RWORK( * ), W( * ) 00018 COMPLEX AB( LDAB, * ), BB( LDBB, * ), WORK( * ), 00019 $ Z( LDZ, * ) 00020 * .. 00021 * 00022 * Purpose 00023 * ======= 00024 * 00025 * CHBGVD computes all the eigenvalues, and optionally, the eigenvectors 00026 * of a complex generalized Hermitian-definite banded eigenproblem, of 00027 * the form A*x=(lambda)*B*x. Here A and B are assumed to be Hermitian 00028 * and banded, and B is also positive definite. If eigenvectors are 00029 * desired, it uses a divide and conquer algorithm. 00030 * 00031 * The divide and conquer algorithm makes very mild assumptions about 00032 * floating point arithmetic. It will work on machines with a guard 00033 * digit in add/subtract, or on those binary machines without guard 00034 * digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or 00035 * Cray-2. It could conceivably fail on hexadecimal or decimal machines 00036 * without guard digits, but we know of none. 00037 * 00038 * Arguments 00039 * ========= 00040 * 00041 * JOBZ (input) CHARACTER*1 00042 * = 'N': Compute eigenvalues only; 00043 * = 'V': Compute eigenvalues and eigenvectors. 00044 * 00045 * UPLO (input) CHARACTER*1 00046 * = 'U': Upper triangles of A and B are stored; 00047 * = 'L': Lower triangles of A and B are stored. 00048 * 00049 * N (input) INTEGER 00050 * The order of the matrices A and B. N >= 0. 00051 * 00052 * KA (input) INTEGER 00053 * The number of superdiagonals of the matrix A if UPLO = 'U', 00054 * or the number of subdiagonals if UPLO = 'L'. KA >= 0. 00055 * 00056 * KB (input) INTEGER 00057 * The number of superdiagonals of the matrix B if UPLO = 'U', 00058 * or the number of subdiagonals if UPLO = 'L'. KB >= 0. 00059 * 00060 * AB (input/output) COMPLEX array, dimension (LDAB, N) 00061 * On entry, the upper or lower triangle of the Hermitian band 00062 * matrix A, stored in the first ka+1 rows of the array. The 00063 * j-th column of A is stored in the j-th column of the array AB 00064 * as follows: 00065 * if UPLO = 'U', AB(ka+1+i-j,j) = A(i,j) for max(1,j-ka)<=i<=j; 00066 * if UPLO = 'L', AB(1+i-j,j) = A(i,j) for j<=i<=min(n,j+ka). 00067 * 00068 * On exit, the contents of AB are destroyed. 00069 * 00070 * LDAB (input) INTEGER 00071 * The leading dimension of the array AB. LDAB >= KA+1. 00072 * 00073 * BB (input/output) COMPLEX array, dimension (LDBB, N) 00074 * On entry, the upper or lower triangle of the Hermitian band 00075 * matrix B, stored in the first kb+1 rows of the array. The 00076 * j-th column of B is stored in the j-th column of the array BB 00077 * as follows: 00078 * if UPLO = 'U', BB(kb+1+i-j,j) = B(i,j) for max(1,j-kb)<=i<=j; 00079 * if UPLO = 'L', BB(1+i-j,j) = B(i,j) for j<=i<=min(n,j+kb). 00080 * 00081 * On exit, the factor S from the split Cholesky factorization 00082 * B = S**H*S, as returned by CPBSTF. 00083 * 00084 * LDBB (input) INTEGER 00085 * The leading dimension of the array BB. LDBB >= KB+1. 00086 * 00087 * W (output) REAL array, dimension (N) 00088 * If INFO = 0, the eigenvalues in ascending order. 00089 * 00090 * Z (output) COMPLEX array, dimension (LDZ, N) 00091 * If JOBZ = 'V', then if INFO = 0, Z contains the matrix Z of 00092 * eigenvectors, with the i-th column of Z holding the 00093 * eigenvector associated with W(i). The eigenvectors are 00094 * normalized so that Z**H*B*Z = I. 00095 * If JOBZ = 'N', then Z is not referenced. 00096 * 00097 * LDZ (input) INTEGER 00098 * The leading dimension of the array Z. LDZ >= 1, and if 00099 * JOBZ = 'V', LDZ >= N. 00100 * 00101 * WORK (workspace/output) COMPLEX array, dimension (MAX(1,LWORK)) 00102 * On exit, if INFO=0, WORK(1) returns the optimal LWORK. 00103 * 00104 * LWORK (input) INTEGER 00105 * The dimension of the array WORK. 00106 * If N <= 1, LWORK >= 1. 00107 * If JOBZ = 'N' and N > 1, LWORK >= N. 00108 * If JOBZ = 'V' and N > 1, LWORK >= 2*N**2. 00109 * 00110 * If LWORK = -1, then a workspace query is assumed; the routine 00111 * only calculates the optimal sizes of the WORK, RWORK and 00112 * IWORK arrays, returns these values as the first entries of 00113 * the WORK, RWORK and IWORK arrays, and no error message 00114 * related to LWORK or LRWORK or LIWORK is issued by XERBLA. 00115 * 00116 * RWORK (workspace/output) REAL array, dimension (MAX(1,LRWORK)) 00117 * On exit, if INFO=0, RWORK(1) returns the optimal LRWORK. 00118 * 00119 * LRWORK (input) INTEGER 00120 * The dimension of array RWORK. 00121 * If N <= 1, LRWORK >= 1. 00122 * If JOBZ = 'N' and N > 1, LRWORK >= N. 00123 * If JOBZ = 'V' and N > 1, LRWORK >= 1 + 5*N + 2*N**2. 00124 * 00125 * If LRWORK = -1, then a workspace query is assumed; the 00126 * routine only calculates the optimal sizes of the WORK, RWORK 00127 * and IWORK arrays, returns these values as the first entries 00128 * of the WORK, RWORK and IWORK arrays, and no error message 00129 * related to LWORK or LRWORK or LIWORK is issued by XERBLA. 00130 * 00131 * IWORK (workspace/output) INTEGER array, dimension (MAX(1,LIWORK)) 00132 * On exit, if INFO=0, IWORK(1) returns the optimal LIWORK. 00133 * 00134 * LIWORK (input) INTEGER 00135 * The dimension of array IWORK. 00136 * If JOBZ = 'N' or N <= 1, LIWORK >= 1. 00137 * If JOBZ = 'V' and N > 1, LIWORK >= 3 + 5*N. 00138 * 00139 * If LIWORK = -1, then a workspace query is assumed; the 00140 * routine only calculates the optimal sizes of the WORK, RWORK 00141 * and IWORK arrays, returns these values as the first entries 00142 * of the WORK, RWORK and IWORK arrays, and no error message 00143 * related to LWORK or LRWORK or LIWORK is issued by XERBLA. 00144 * 00145 * INFO (output) INTEGER 00146 * = 0: successful exit 00147 * < 0: if INFO = -i, the i-th argument had an illegal value 00148 * > 0: if INFO = i, and i is: 00149 * <= N: the algorithm failed to converge: 00150 * i off-diagonal elements of an intermediate 00151 * tridiagonal form did not converge to zero; 00152 * > N: if INFO = N + i, for 1 <= i <= N, then CPBSTF 00153 * returned INFO = i: B is not positive definite. 00154 * The factorization of B could not be completed and 00155 * no eigenvalues or eigenvectors were computed. 00156 * 00157 * Further Details 00158 * =============== 00159 * 00160 * Based on contributions by 00161 * Mark Fahey, Department of Mathematics, Univ. of Kentucky, USA 00162 * 00163 * ===================================================================== 00164 * 00165 * .. Parameters .. 00166 COMPLEX CONE, CZERO 00167 PARAMETER ( CONE = ( 1.0E+0, 0.0E+0 ), 00168 $ CZERO = ( 0.0E+0, 0.0E+0 ) ) 00169 * .. 00170 * .. Local Scalars .. 00171 LOGICAL LQUERY, UPPER, WANTZ 00172 CHARACTER VECT 00173 INTEGER IINFO, INDE, INDWK2, INDWRK, LIWMIN, LLRWK, 00174 $ LLWK2, LRWMIN, LWMIN 00175 * .. 00176 * .. External Functions .. 00177 LOGICAL LSAME 00178 EXTERNAL LSAME 00179 * .. 00180 * .. External Subroutines .. 00181 EXTERNAL CGEMM, CHBGST, CHBTRD, CLACPY, CPBSTF, CSTEDC, 00182 $ SSTERF, XERBLA 00183 * .. 00184 * .. Executable Statements .. 00185 * 00186 * Test the input parameters. 00187 * 00188 WANTZ = LSAME( JOBZ, 'V' ) 00189 UPPER = LSAME( UPLO, 'U' ) 00190 LQUERY = ( LWORK.EQ.-1 .OR. LRWORK.EQ.-1 .OR. LIWORK.EQ.-1 ) 00191 * 00192 INFO = 0 00193 IF( N.LE.1 ) THEN 00194 LWMIN = 1 00195 LRWMIN = 1 00196 LIWMIN = 1 00197 ELSE IF( WANTZ ) THEN 00198 LWMIN = 2*N**2 00199 LRWMIN = 1 + 5*N + 2*N**2 00200 LIWMIN = 3 + 5*N 00201 ELSE 00202 LWMIN = N 00203 LRWMIN = N 00204 LIWMIN = 1 00205 END IF 00206 IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN 00207 INFO = -1 00208 ELSE IF( .NOT.( UPPER .OR. LSAME( UPLO, 'L' ) ) ) THEN 00209 INFO = -2 00210 ELSE IF( N.LT.0 ) THEN 00211 INFO = -3 00212 ELSE IF( KA.LT.0 ) THEN 00213 INFO = -4 00214 ELSE IF( KB.LT.0 .OR. KB.GT.KA ) THEN 00215 INFO = -5 00216 ELSE IF( LDAB.LT.KA+1 ) THEN 00217 INFO = -7 00218 ELSE IF( LDBB.LT.KB+1 ) THEN 00219 INFO = -9 00220 ELSE IF( LDZ.LT.1 .OR. ( WANTZ .AND. LDZ.LT.N ) ) THEN 00221 INFO = -12 00222 END IF 00223 * 00224 IF( INFO.EQ.0 ) THEN 00225 WORK( 1 ) = LWMIN 00226 RWORK( 1 ) = LRWMIN 00227 IWORK( 1 ) = LIWMIN 00228 * 00229 IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN 00230 INFO = -14 00231 ELSE IF( LRWORK.LT.LRWMIN .AND. .NOT.LQUERY ) THEN 00232 INFO = -16 00233 ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT.LQUERY ) THEN 00234 INFO = -18 00235 END IF 00236 END IF 00237 * 00238 IF( INFO.NE.0 ) THEN 00239 CALL XERBLA( 'CHBGVD', -INFO ) 00240 RETURN 00241 ELSE IF( LQUERY ) THEN 00242 RETURN 00243 END IF 00244 * 00245 * Quick return if possible 00246 * 00247 IF( N.EQ.0 ) 00248 $ RETURN 00249 * 00250 * Form a split Cholesky factorization of B. 00251 * 00252 CALL CPBSTF( UPLO, N, KB, BB, LDBB, INFO ) 00253 IF( INFO.NE.0 ) THEN 00254 INFO = N + INFO 00255 RETURN 00256 END IF 00257 * 00258 * Transform problem to standard eigenvalue problem. 00259 * 00260 INDE = 1 00261 INDWRK = INDE + N 00262 INDWK2 = 1 + N*N 00263 LLWK2 = LWORK - INDWK2 + 2 00264 LLRWK = LRWORK - INDWRK + 2 00265 CALL CHBGST( JOBZ, UPLO, N, KA, KB, AB, LDAB, BB, LDBB, Z, LDZ, 00266 $ WORK, RWORK( INDWRK ), IINFO ) 00267 * 00268 * Reduce Hermitian band matrix to tridiagonal form. 00269 * 00270 IF( WANTZ ) THEN 00271 VECT = 'U' 00272 ELSE 00273 VECT = 'N' 00274 END IF 00275 CALL CHBTRD( VECT, UPLO, N, KA, AB, LDAB, W, RWORK( INDE ), Z, 00276 $ LDZ, WORK, IINFO ) 00277 * 00278 * For eigenvalues only, call SSTERF. For eigenvectors, call CSTEDC. 00279 * 00280 IF( .NOT.WANTZ ) THEN 00281 CALL SSTERF( N, W, RWORK( INDE ), INFO ) 00282 ELSE 00283 CALL CSTEDC( 'I', N, W, RWORK( INDE ), WORK, N, WORK( INDWK2 ), 00284 $ LLWK2, RWORK( INDWRK ), LLRWK, IWORK, LIWORK, 00285 $ INFO ) 00286 CALL CGEMM( 'N', 'N', N, N, N, CONE, Z, LDZ, WORK, N, CZERO, 00287 $ WORK( INDWK2 ), N ) 00288 CALL CLACPY( 'A', N, N, WORK( INDWK2 ), N, Z, LDZ ) 00289 END IF 00290 * 00291 WORK( 1 ) = LWMIN 00292 RWORK( 1 ) = LRWMIN 00293 IWORK( 1 ) = LIWMIN 00294 RETURN 00295 * 00296 * End of CHBGVD 00297 * 00298 END
1.7.2