00001 SUBROUTINE CHBEVD( JOBZ, UPLO, N, KD, AB, LDAB, W, Z, LDZ, WORK, 00002 $ LWORK, RWORK, LRWORK, IWORK, LIWORK, INFO ) 00003 * 00004 * -- LAPACK driver routine (version 3.2) -- 00005 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00006 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00007 * November 2006 00008 * 00009 * .. Scalar Arguments .. 00010 CHARACTER JOBZ, UPLO 00011 INTEGER INFO, KD, LDAB, LDZ, LIWORK, LRWORK, LWORK, N 00012 * .. 00013 * .. Array Arguments .. 00014 INTEGER IWORK( * ) 00015 REAL RWORK( * ), W( * ) 00016 COMPLEX AB( LDAB, * ), WORK( * ), Z( LDZ, * ) 00017 * .. 00018 * 00019 * Purpose 00020 * ======= 00021 * 00022 * CHBEVD computes all the eigenvalues and, optionally, eigenvectors of 00023 * a complex Hermitian band matrix A. If eigenvectors are desired, it 00024 * uses a divide and conquer algorithm. 00025 * 00026 * The divide and conquer algorithm makes very mild assumptions about 00027 * floating point arithmetic. It will work on machines with a guard 00028 * digit in add/subtract, or on those binary machines without guard 00029 * digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or 00030 * Cray-2. It could conceivably fail on hexadecimal or decimal machines 00031 * without guard digits, but we know of none. 00032 * 00033 * Arguments 00034 * ========= 00035 * 00036 * JOBZ (input) CHARACTER*1 00037 * = 'N': Compute eigenvalues only; 00038 * = 'V': Compute eigenvalues and eigenvectors. 00039 * 00040 * UPLO (input) CHARACTER*1 00041 * = 'U': Upper triangle of A is stored; 00042 * = 'L': Lower triangle of A is stored. 00043 * 00044 * N (input) INTEGER 00045 * The order of the matrix A. N >= 0. 00046 * 00047 * KD (input) INTEGER 00048 * The number of superdiagonals of the matrix A if UPLO = 'U', 00049 * or the number of subdiagonals if UPLO = 'L'. KD >= 0. 00050 * 00051 * AB (input/output) COMPLEX array, dimension (LDAB, N) 00052 * On entry, the upper or lower triangle of the Hermitian band 00053 * matrix A, stored in the first KD+1 rows of the array. The 00054 * j-th column of A is stored in the j-th column of the array AB 00055 * as follows: 00056 * if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j; 00057 * if UPLO = 'L', AB(1+i-j,j) = A(i,j) for j<=i<=min(n,j+kd). 00058 * 00059 * On exit, AB is overwritten by values generated during the 00060 * reduction to tridiagonal form. If UPLO = 'U', the first 00061 * superdiagonal and the diagonal of the tridiagonal matrix T 00062 * are returned in rows KD and KD+1 of AB, and if UPLO = 'L', 00063 * the diagonal and first subdiagonal of T are returned in the 00064 * first two rows of AB. 00065 * 00066 * LDAB (input) INTEGER 00067 * The leading dimension of the array AB. LDAB >= KD + 1. 00068 * 00069 * W (output) REAL array, dimension (N) 00070 * If INFO = 0, the eigenvalues in ascending order. 00071 * 00072 * Z (output) COMPLEX array, dimension (LDZ, N) 00073 * If JOBZ = 'V', then if INFO = 0, Z contains the orthonormal 00074 * eigenvectors of the matrix A, with the i-th column of Z 00075 * holding the eigenvector associated with W(i). 00076 * If JOBZ = 'N', then Z is not referenced. 00077 * 00078 * LDZ (input) INTEGER 00079 * The leading dimension of the array Z. LDZ >= 1, and if 00080 * JOBZ = 'V', LDZ >= max(1,N). 00081 * 00082 * WORK (workspace/output) COMPLEX array, dimension (MAX(1,LWORK)) 00083 * On exit, if INFO = 0, WORK(1) returns the optimal LWORK. 00084 * 00085 * LWORK (input) INTEGER 00086 * The dimension of the array WORK. 00087 * If N <= 1, LWORK must be at least 1. 00088 * If JOBZ = 'N' and N > 1, LWORK must be at least N. 00089 * If JOBZ = 'V' and N > 1, LWORK must be at least 2*N**2. 00090 * 00091 * If LWORK = -1, then a workspace query is assumed; the routine 00092 * only calculates the optimal sizes of the WORK, RWORK and 00093 * IWORK arrays, returns these values as the first entries of 00094 * the WORK, RWORK and IWORK arrays, and no error message 00095 * related to LWORK or LRWORK or LIWORK is issued by XERBLA. 00096 * 00097 * RWORK (workspace/output) REAL array, 00098 * dimension (LRWORK) 00099 * On exit, if INFO = 0, RWORK(1) returns the optimal LRWORK. 00100 * 00101 * LRWORK (input) INTEGER 00102 * The dimension of array RWORK. 00103 * If N <= 1, LRWORK must be at least 1. 00104 * If JOBZ = 'N' and N > 1, LRWORK must be at least N. 00105 * If JOBZ = 'V' and N > 1, LRWORK must be at least 00106 * 1 + 5*N + 2*N**2. 00107 * 00108 * If LRWORK = -1, then a workspace query is assumed; the 00109 * routine only calculates the optimal sizes of the WORK, RWORK 00110 * and IWORK arrays, returns these values as the first entries 00111 * of the WORK, RWORK and IWORK arrays, and no error message 00112 * related to LWORK or LRWORK or LIWORK is issued by XERBLA. 00113 * 00114 * IWORK (workspace/output) INTEGER array, dimension (MAX(1,LIWORK)) 00115 * On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK. 00116 * 00117 * LIWORK (input) INTEGER 00118 * The dimension of array IWORK. 00119 * If JOBZ = 'N' or N <= 1, LIWORK must be at least 1. 00120 * If JOBZ = 'V' and N > 1, LIWORK must be at least 3 + 5*N . 00121 * 00122 * If LIWORK = -1, then a workspace query is assumed; the 00123 * routine only calculates the optimal sizes of the WORK, RWORK 00124 * and IWORK arrays, returns these values as the first entries 00125 * of the WORK, RWORK and IWORK arrays, and no error message 00126 * related to LWORK or LRWORK or LIWORK is issued by XERBLA. 00127 * 00128 * INFO (output) INTEGER 00129 * = 0: successful exit. 00130 * < 0: if INFO = -i, the i-th argument had an illegal value. 00131 * > 0: if INFO = i, the algorithm failed to converge; i 00132 * off-diagonal elements of an intermediate tridiagonal 00133 * form did not converge to zero. 00134 * 00135 * ===================================================================== 00136 * 00137 * .. Parameters .. 00138 REAL ZERO, ONE 00139 PARAMETER ( ZERO = 0.0E0, ONE = 1.0E0 ) 00140 COMPLEX CZERO, CONE 00141 PARAMETER ( CZERO = ( 0.0E0, 0.0E0 ), 00142 $ CONE = ( 1.0E0, 0.0E0 ) ) 00143 * .. 00144 * .. Local Scalars .. 00145 LOGICAL LOWER, LQUERY, WANTZ 00146 INTEGER IINFO, IMAX, INDE, INDWK2, INDWRK, ISCALE, 00147 $ LIWMIN, LLRWK, LLWK2, LRWMIN, LWMIN 00148 REAL ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA, 00149 $ SMLNUM 00150 * .. 00151 * .. External Functions .. 00152 LOGICAL LSAME 00153 REAL CLANHB, SLAMCH 00154 EXTERNAL LSAME, CLANHB, SLAMCH 00155 * .. 00156 * .. External Subroutines .. 00157 EXTERNAL CGEMM, CHBTRD, CLACPY, CLASCL, CSTEDC, SSCAL, 00158 $ SSTERF, XERBLA 00159 * .. 00160 * .. Intrinsic Functions .. 00161 INTRINSIC SQRT 00162 * .. 00163 * .. Executable Statements .. 00164 * 00165 * Test the input parameters. 00166 * 00167 WANTZ = LSAME( JOBZ, 'V' ) 00168 LOWER = LSAME( UPLO, 'L' ) 00169 LQUERY = ( LWORK.EQ.-1 .OR. LIWORK.EQ.-1 .OR. LRWORK.EQ.-1 ) 00170 * 00171 INFO = 0 00172 IF( N.LE.1 ) THEN 00173 LWMIN = 1 00174 LRWMIN = 1 00175 LIWMIN = 1 00176 ELSE 00177 IF( WANTZ ) THEN 00178 LWMIN = 2*N**2 00179 LRWMIN = 1 + 5*N + 2*N**2 00180 LIWMIN = 3 + 5*N 00181 ELSE 00182 LWMIN = N 00183 LRWMIN = N 00184 LIWMIN = 1 00185 END IF 00186 END IF 00187 IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN 00188 INFO = -1 00189 ELSE IF( .NOT.( LOWER .OR. LSAME( UPLO, 'U' ) ) ) THEN 00190 INFO = -2 00191 ELSE IF( N.LT.0 ) THEN 00192 INFO = -3 00193 ELSE IF( KD.LT.0 ) THEN 00194 INFO = -4 00195 ELSE IF( LDAB.LT.KD+1 ) THEN 00196 INFO = -6 00197 ELSE IF( LDZ.LT.1 .OR. ( WANTZ .AND. LDZ.LT.N ) ) THEN 00198 INFO = -9 00199 END IF 00200 * 00201 IF( INFO.EQ.0 ) THEN 00202 WORK( 1 ) = LWMIN 00203 RWORK( 1 ) = LRWMIN 00204 IWORK( 1 ) = LIWMIN 00205 * 00206 IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN 00207 INFO = -11 00208 ELSE IF( LRWORK.LT.LRWMIN .AND. .NOT.LQUERY ) THEN 00209 INFO = -13 00210 ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT.LQUERY ) THEN 00211 INFO = -15 00212 END IF 00213 END IF 00214 * 00215 IF( INFO.NE.0 ) THEN 00216 CALL XERBLA( 'CHBEVD', -INFO ) 00217 RETURN 00218 ELSE IF( LQUERY ) THEN 00219 RETURN 00220 END IF 00221 * 00222 * Quick return if possible 00223 * 00224 IF( N.EQ.0 ) 00225 $ RETURN 00226 * 00227 IF( N.EQ.1 ) THEN 00228 W( 1 ) = AB( 1, 1 ) 00229 IF( WANTZ ) 00230 $ Z( 1, 1 ) = CONE 00231 RETURN 00232 END IF 00233 * 00234 * Get machine constants. 00235 * 00236 SAFMIN = SLAMCH( 'Safe minimum' ) 00237 EPS = SLAMCH( 'Precision' ) 00238 SMLNUM = SAFMIN / EPS 00239 BIGNUM = ONE / SMLNUM 00240 RMIN = SQRT( SMLNUM ) 00241 RMAX = SQRT( BIGNUM ) 00242 * 00243 * Scale matrix to allowable range, if necessary. 00244 * 00245 ANRM = CLANHB( 'M', UPLO, N, KD, AB, LDAB, RWORK ) 00246 ISCALE = 0 00247 IF( ANRM.GT.ZERO .AND. ANRM.LT.RMIN ) THEN 00248 ISCALE = 1 00249 SIGMA = RMIN / ANRM 00250 ELSE IF( ANRM.GT.RMAX ) THEN 00251 ISCALE = 1 00252 SIGMA = RMAX / ANRM 00253 END IF 00254 IF( ISCALE.EQ.1 ) THEN 00255 IF( LOWER ) THEN 00256 CALL CLASCL( 'B', KD, KD, ONE, SIGMA, N, N, AB, LDAB, INFO ) 00257 ELSE 00258 CALL CLASCL( 'Q', KD, KD, ONE, SIGMA, N, N, AB, LDAB, INFO ) 00259 END IF 00260 END IF 00261 * 00262 * Call CHBTRD to reduce Hermitian band matrix to tridiagonal form. 00263 * 00264 INDE = 1 00265 INDWRK = INDE + N 00266 INDWK2 = 1 + N*N 00267 LLWK2 = LWORK - INDWK2 + 1 00268 LLRWK = LRWORK - INDWRK + 1 00269 CALL CHBTRD( JOBZ, UPLO, N, KD, AB, LDAB, W, RWORK( INDE ), Z, 00270 $ LDZ, WORK, IINFO ) 00271 * 00272 * For eigenvalues only, call SSTERF. For eigenvectors, call CSTEDC. 00273 * 00274 IF( .NOT.WANTZ ) THEN 00275 CALL SSTERF( N, W, RWORK( INDE ), INFO ) 00276 ELSE 00277 CALL CSTEDC( 'I', N, W, RWORK( INDE ), WORK, N, WORK( INDWK2 ), 00278 $ LLWK2, RWORK( INDWRK ), LLRWK, IWORK, LIWORK, 00279 $ INFO ) 00280 CALL CGEMM( 'N', 'N', N, N, N, CONE, Z, LDZ, WORK, N, CZERO, 00281 $ WORK( INDWK2 ), N ) 00282 CALL CLACPY( 'A', N, N, WORK( INDWK2 ), N, Z, LDZ ) 00283 END IF 00284 * 00285 * If matrix was scaled, then rescale eigenvalues appropriately. 00286 * 00287 IF( ISCALE.EQ.1 ) THEN 00288 IF( INFO.EQ.0 ) THEN 00289 IMAX = N 00290 ELSE 00291 IMAX = INFO - 1 00292 END IF 00293 CALL SSCAL( IMAX, ONE / SIGMA, W, 1 ) 00294 END IF 00295 * 00296 WORK( 1 ) = LWMIN 00297 RWORK( 1 ) = LRWMIN 00298 IWORK( 1 ) = LIWMIN 00299 RETURN 00300 * 00301 * End of CHBEVD 00302 * 00303 END
1.7.2