![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief <b> CHBEVD computes the eigenvalues and, optionally, the left and/or right eigenvectors for OTHER matrices</b> 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download CHBEVD + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/chbevd.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/chbevd.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/chbevd.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE CHBEVD( JOBZ, UPLO, N, KD, AB, LDAB, W, Z, LDZ, WORK, 00022 * LWORK, RWORK, LRWORK, IWORK, LIWORK, INFO ) 00023 * 00024 * .. Scalar Arguments .. 00025 * CHARACTER JOBZ, UPLO 00026 * INTEGER INFO, KD, LDAB, LDZ, LIWORK, LRWORK, LWORK, N 00027 * .. 00028 * .. Array Arguments .. 00029 * INTEGER IWORK( * ) 00030 * REAL RWORK( * ), W( * ) 00031 * COMPLEX AB( LDAB, * ), WORK( * ), Z( LDZ, * ) 00032 * .. 00033 * 00034 * 00035 *> \par Purpose: 00036 * ============= 00037 *> 00038 *> \verbatim 00039 *> 00040 *> CHBEVD computes all the eigenvalues and, optionally, eigenvectors of 00041 *> a complex Hermitian band matrix A. If eigenvectors are desired, it 00042 *> uses a divide and conquer algorithm. 00043 *> 00044 *> The divide and conquer algorithm makes very mild assumptions about 00045 *> floating point arithmetic. It will work on machines with a guard 00046 *> digit in add/subtract, or on those binary machines without guard 00047 *> digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or 00048 *> Cray-2. It could conceivably fail on hexadecimal or decimal machines 00049 *> without guard digits, but we know of none. 00050 *> \endverbatim 00051 * 00052 * Arguments: 00053 * ========== 00054 * 00055 *> \param[in] JOBZ 00056 *> \verbatim 00057 *> JOBZ is CHARACTER*1 00058 *> = 'N': Compute eigenvalues only; 00059 *> = 'V': Compute eigenvalues and eigenvectors. 00060 *> \endverbatim 00061 *> 00062 *> \param[in] UPLO 00063 *> \verbatim 00064 *> UPLO is CHARACTER*1 00065 *> = 'U': Upper triangle of A is stored; 00066 *> = 'L': Lower triangle of A is stored. 00067 *> \endverbatim 00068 *> 00069 *> \param[in] N 00070 *> \verbatim 00071 *> N is INTEGER 00072 *> The order of the matrix A. N >= 0. 00073 *> \endverbatim 00074 *> 00075 *> \param[in] KD 00076 *> \verbatim 00077 *> KD is INTEGER 00078 *> The number of superdiagonals of the matrix A if UPLO = 'U', 00079 *> or the number of subdiagonals if UPLO = 'L'. KD >= 0. 00080 *> \endverbatim 00081 *> 00082 *> \param[in,out] AB 00083 *> \verbatim 00084 *> AB is COMPLEX array, dimension (LDAB, N) 00085 *> On entry, the upper or lower triangle of the Hermitian band 00086 *> matrix A, stored in the first KD+1 rows of the array. The 00087 *> j-th column of A is stored in the j-th column of the array AB 00088 *> as follows: 00089 *> if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j; 00090 *> if UPLO = 'L', AB(1+i-j,j) = A(i,j) for j<=i<=min(n,j+kd). 00091 *> 00092 *> On exit, AB is overwritten by values generated during the 00093 *> reduction to tridiagonal form. If UPLO = 'U', the first 00094 *> superdiagonal and the diagonal of the tridiagonal matrix T 00095 *> are returned in rows KD and KD+1 of AB, and if UPLO = 'L', 00096 *> the diagonal and first subdiagonal of T are returned in the 00097 *> first two rows of AB. 00098 *> \endverbatim 00099 *> 00100 *> \param[in] LDAB 00101 *> \verbatim 00102 *> LDAB is INTEGER 00103 *> The leading dimension of the array AB. LDAB >= KD + 1. 00104 *> \endverbatim 00105 *> 00106 *> \param[out] W 00107 *> \verbatim 00108 *> W is REAL array, dimension (N) 00109 *> If INFO = 0, the eigenvalues in ascending order. 00110 *> \endverbatim 00111 *> 00112 *> \param[out] Z 00113 *> \verbatim 00114 *> Z is COMPLEX array, dimension (LDZ, N) 00115 *> If JOBZ = 'V', then if INFO = 0, Z contains the orthonormal 00116 *> eigenvectors of the matrix A, with the i-th column of Z 00117 *> holding the eigenvector associated with W(i). 00118 *> If JOBZ = 'N', then Z is not referenced. 00119 *> \endverbatim 00120 *> 00121 *> \param[in] LDZ 00122 *> \verbatim 00123 *> LDZ is INTEGER 00124 *> The leading dimension of the array Z. LDZ >= 1, and if 00125 *> JOBZ = 'V', LDZ >= max(1,N). 00126 *> \endverbatim 00127 *> 00128 *> \param[out] WORK 00129 *> \verbatim 00130 *> WORK is COMPLEX array, dimension (MAX(1,LWORK)) 00131 *> On exit, if INFO = 0, WORK(1) returns the optimal LWORK. 00132 *> \endverbatim 00133 *> 00134 *> \param[in] LWORK 00135 *> \verbatim 00136 *> LWORK is INTEGER 00137 *> The dimension of the array WORK. 00138 *> If N <= 1, LWORK must be at least 1. 00139 *> If JOBZ = 'N' and N > 1, LWORK must be at least N. 00140 *> If JOBZ = 'V' and N > 1, LWORK must be at least 2*N**2. 00141 *> 00142 *> If LWORK = -1, then a workspace query is assumed; the routine 00143 *> only calculates the optimal sizes of the WORK, RWORK and 00144 *> IWORK arrays, returns these values as the first entries of 00145 *> the WORK, RWORK and IWORK arrays, and no error message 00146 *> related to LWORK or LRWORK or LIWORK is issued by XERBLA. 00147 *> \endverbatim 00148 *> 00149 *> \param[out] RWORK 00150 *> \verbatim 00151 *> RWORK is REAL array, 00152 *> dimension (LRWORK) 00153 *> On exit, if INFO = 0, RWORK(1) returns the optimal LRWORK. 00154 *> \endverbatim 00155 *> 00156 *> \param[in] LRWORK 00157 *> \verbatim 00158 *> LRWORK is INTEGER 00159 *> The dimension of array RWORK. 00160 *> If N <= 1, LRWORK must be at least 1. 00161 *> If JOBZ = 'N' and N > 1, LRWORK must be at least N. 00162 *> If JOBZ = 'V' and N > 1, LRWORK must be at least 00163 *> 1 + 5*N + 2*N**2. 00164 *> 00165 *> If LRWORK = -1, then a workspace query is assumed; the 00166 *> routine only calculates the optimal sizes of the WORK, RWORK 00167 *> and IWORK arrays, returns these values as the first entries 00168 *> of the WORK, RWORK and IWORK arrays, and no error message 00169 *> related to LWORK or LRWORK or LIWORK is issued by XERBLA. 00170 *> \endverbatim 00171 *> 00172 *> \param[out] IWORK 00173 *> \verbatim 00174 *> IWORK is INTEGER array, dimension (MAX(1,LIWORK)) 00175 *> On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK. 00176 *> \endverbatim 00177 *> 00178 *> \param[in] LIWORK 00179 *> \verbatim 00180 *> LIWORK is INTEGER 00181 *> The dimension of array IWORK. 00182 *> If JOBZ = 'N' or N <= 1, LIWORK must be at least 1. 00183 *> If JOBZ = 'V' and N > 1, LIWORK must be at least 3 + 5*N . 00184 *> 00185 *> If LIWORK = -1, then a workspace query is assumed; the 00186 *> routine only calculates the optimal sizes of the WORK, RWORK 00187 *> and IWORK arrays, returns these values as the first entries 00188 *> of the WORK, RWORK and IWORK arrays, and no error message 00189 *> related to LWORK or LRWORK or LIWORK is issued by XERBLA. 00190 *> \endverbatim 00191 *> 00192 *> \param[out] INFO 00193 *> \verbatim 00194 *> INFO is INTEGER 00195 *> = 0: successful exit. 00196 *> < 0: if INFO = -i, the i-th argument had an illegal value. 00197 *> > 0: if INFO = i, the algorithm failed to converge; i 00198 *> off-diagonal elements of an intermediate tridiagonal 00199 *> form did not converge to zero. 00200 *> \endverbatim 00201 * 00202 * Authors: 00203 * ======== 00204 * 00205 *> \author Univ. of Tennessee 00206 *> \author Univ. of California Berkeley 00207 *> \author Univ. of Colorado Denver 00208 *> \author NAG Ltd. 00209 * 00210 *> \date November 2011 00211 * 00212 *> \ingroup complexOTHEReigen 00213 * 00214 * ===================================================================== 00215 SUBROUTINE CHBEVD( JOBZ, UPLO, N, KD, AB, LDAB, W, Z, LDZ, WORK, 00216 $ LWORK, RWORK, LRWORK, IWORK, LIWORK, INFO ) 00217 * 00218 * -- LAPACK driver routine (version 3.4.0) -- 00219 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00220 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00221 * November 2011 00222 * 00223 * .. Scalar Arguments .. 00224 CHARACTER JOBZ, UPLO 00225 INTEGER INFO, KD, LDAB, LDZ, LIWORK, LRWORK, LWORK, N 00226 * .. 00227 * .. Array Arguments .. 00228 INTEGER IWORK( * ) 00229 REAL RWORK( * ), W( * ) 00230 COMPLEX AB( LDAB, * ), WORK( * ), Z( LDZ, * ) 00231 * .. 00232 * 00233 * ===================================================================== 00234 * 00235 * .. Parameters .. 00236 REAL ZERO, ONE 00237 PARAMETER ( ZERO = 0.0E0, ONE = 1.0E0 ) 00238 COMPLEX CZERO, CONE 00239 PARAMETER ( CZERO = ( 0.0E0, 0.0E0 ), 00240 $ CONE = ( 1.0E0, 0.0E0 ) ) 00241 * .. 00242 * .. Local Scalars .. 00243 LOGICAL LOWER, LQUERY, WANTZ 00244 INTEGER IINFO, IMAX, INDE, INDWK2, INDWRK, ISCALE, 00245 $ LIWMIN, LLRWK, LLWK2, LRWMIN, LWMIN 00246 REAL ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA, 00247 $ SMLNUM 00248 * .. 00249 * .. External Functions .. 00250 LOGICAL LSAME 00251 REAL CLANHB, SLAMCH 00252 EXTERNAL LSAME, CLANHB, SLAMCH 00253 * .. 00254 * .. External Subroutines .. 00255 EXTERNAL CGEMM, CHBTRD, CLACPY, CLASCL, CSTEDC, SSCAL, 00256 $ SSTERF, XERBLA 00257 * .. 00258 * .. Intrinsic Functions .. 00259 INTRINSIC SQRT 00260 * .. 00261 * .. Executable Statements .. 00262 * 00263 * Test the input parameters. 00264 * 00265 WANTZ = LSAME( JOBZ, 'V' ) 00266 LOWER = LSAME( UPLO, 'L' ) 00267 LQUERY = ( LWORK.EQ.-1 .OR. LIWORK.EQ.-1 .OR. LRWORK.EQ.-1 ) 00268 * 00269 INFO = 0 00270 IF( N.LE.1 ) THEN 00271 LWMIN = 1 00272 LRWMIN = 1 00273 LIWMIN = 1 00274 ELSE 00275 IF( WANTZ ) THEN 00276 LWMIN = 2*N**2 00277 LRWMIN = 1 + 5*N + 2*N**2 00278 LIWMIN = 3 + 5*N 00279 ELSE 00280 LWMIN = N 00281 LRWMIN = N 00282 LIWMIN = 1 00283 END IF 00284 END IF 00285 IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN 00286 INFO = -1 00287 ELSE IF( .NOT.( LOWER .OR. LSAME( UPLO, 'U' ) ) ) THEN 00288 INFO = -2 00289 ELSE IF( N.LT.0 ) THEN 00290 INFO = -3 00291 ELSE IF( KD.LT.0 ) THEN 00292 INFO = -4 00293 ELSE IF( LDAB.LT.KD+1 ) THEN 00294 INFO = -6 00295 ELSE IF( LDZ.LT.1 .OR. ( WANTZ .AND. LDZ.LT.N ) ) THEN 00296 INFO = -9 00297 END IF 00298 * 00299 IF( INFO.EQ.0 ) THEN 00300 WORK( 1 ) = LWMIN 00301 RWORK( 1 ) = LRWMIN 00302 IWORK( 1 ) = LIWMIN 00303 * 00304 IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN 00305 INFO = -11 00306 ELSE IF( LRWORK.LT.LRWMIN .AND. .NOT.LQUERY ) THEN 00307 INFO = -13 00308 ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT.LQUERY ) THEN 00309 INFO = -15 00310 END IF 00311 END IF 00312 * 00313 IF( INFO.NE.0 ) THEN 00314 CALL XERBLA( 'CHBEVD', -INFO ) 00315 RETURN 00316 ELSE IF( LQUERY ) THEN 00317 RETURN 00318 END IF 00319 * 00320 * Quick return if possible 00321 * 00322 IF( N.EQ.0 ) 00323 $ RETURN 00324 * 00325 IF( N.EQ.1 ) THEN 00326 W( 1 ) = AB( 1, 1 ) 00327 IF( WANTZ ) 00328 $ Z( 1, 1 ) = CONE 00329 RETURN 00330 END IF 00331 * 00332 * Get machine constants. 00333 * 00334 SAFMIN = SLAMCH( 'Safe minimum' ) 00335 EPS = SLAMCH( 'Precision' ) 00336 SMLNUM = SAFMIN / EPS 00337 BIGNUM = ONE / SMLNUM 00338 RMIN = SQRT( SMLNUM ) 00339 RMAX = SQRT( BIGNUM ) 00340 * 00341 * Scale matrix to allowable range, if necessary. 00342 * 00343 ANRM = CLANHB( 'M', UPLO, N, KD, AB, LDAB, RWORK ) 00344 ISCALE = 0 00345 IF( ANRM.GT.ZERO .AND. ANRM.LT.RMIN ) THEN 00346 ISCALE = 1 00347 SIGMA = RMIN / ANRM 00348 ELSE IF( ANRM.GT.RMAX ) THEN 00349 ISCALE = 1 00350 SIGMA = RMAX / ANRM 00351 END IF 00352 IF( ISCALE.EQ.1 ) THEN 00353 IF( LOWER ) THEN 00354 CALL CLASCL( 'B', KD, KD, ONE, SIGMA, N, N, AB, LDAB, INFO ) 00355 ELSE 00356 CALL CLASCL( 'Q', KD, KD, ONE, SIGMA, N, N, AB, LDAB, INFO ) 00357 END IF 00358 END IF 00359 * 00360 * Call CHBTRD to reduce Hermitian band matrix to tridiagonal form. 00361 * 00362 INDE = 1 00363 INDWRK = INDE + N 00364 INDWK2 = 1 + N*N 00365 LLWK2 = LWORK - INDWK2 + 1 00366 LLRWK = LRWORK - INDWRK + 1 00367 CALL CHBTRD( JOBZ, UPLO, N, KD, AB, LDAB, W, RWORK( INDE ), Z, 00368 $ LDZ, WORK, IINFO ) 00369 * 00370 * For eigenvalues only, call SSTERF. For eigenvectors, call CSTEDC. 00371 * 00372 IF( .NOT.WANTZ ) THEN 00373 CALL SSTERF( N, W, RWORK( INDE ), INFO ) 00374 ELSE 00375 CALL CSTEDC( 'I', N, W, RWORK( INDE ), WORK, N, WORK( INDWK2 ), 00376 $ LLWK2, RWORK( INDWRK ), LLRWK, IWORK, LIWORK, 00377 $ INFO ) 00378 CALL CGEMM( 'N', 'N', N, N, N, CONE, Z, LDZ, WORK, N, CZERO, 00379 $ WORK( INDWK2 ), N ) 00380 CALL CLACPY( 'A', N, N, WORK( INDWK2 ), N, Z, LDZ ) 00381 END IF 00382 * 00383 * If matrix was scaled, then rescale eigenvalues appropriately. 00384 * 00385 IF( ISCALE.EQ.1 ) THEN 00386 IF( INFO.EQ.0 ) THEN 00387 IMAX = N 00388 ELSE 00389 IMAX = INFO - 1 00390 END IF 00391 CALL SSCAL( IMAX, ONE / SIGMA, W, 1 ) 00392 END IF 00393 * 00394 WORK( 1 ) = LWMIN 00395 RWORK( 1 ) = LRWMIN 00396 IWORK( 1 ) = LIWMIN 00397 RETURN 00398 * 00399 * End of CHBEVD 00400 * 00401 END