1*> \brief <b> CHBEVD computes the eigenvalues and, optionally, the left and/or right eigenvectors for OTHER matrices</b> 2* 3* =========== DOCUMENTATION =========== 4* 5* Online html documentation available at 6* http://www.netlib.org/lapack/explore-html/ 7* 8*> \htmlonly 9*> Download CHBEVD + dependencies 10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/chbevd.f"> 11*> [TGZ]</a> 12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/chbevd.f"> 13*> [ZIP]</a> 14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/chbevd.f"> 15*> [TXT]</a> 16*> \endhtmlonly 17* 18* Definition: 19* =========== 20* 21* SUBROUTINE CHBEVD( JOBZ, UPLO, N, KD, AB, LDAB, W, Z, LDZ, WORK, 22* LWORK, RWORK, LRWORK, IWORK, LIWORK, INFO ) 23* 24* .. Scalar Arguments .. 25* CHARACTER JOBZ, UPLO 26* INTEGER INFO, KD, LDAB, LDZ, LIWORK, LRWORK, LWORK, N 27* .. 28* .. Array Arguments .. 29* INTEGER IWORK( * ) 30* REAL RWORK( * ), W( * ) 31* COMPLEX AB( LDAB, * ), WORK( * ), Z( LDZ, * ) 32* .. 33* 34* 35*> \par Purpose: 36* ============= 37*> 38*> \verbatim 39*> 40*> CHBEVD computes all the eigenvalues and, optionally, eigenvectors of 41*> a complex Hermitian band matrix A. If eigenvectors are desired, it 42*> uses a divide and conquer algorithm. 43*> 44*> The divide and conquer algorithm makes very mild assumptions about 45*> floating point arithmetic. It will work on machines with a guard 46*> digit in add/subtract, or on those binary machines without guard 47*> digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or 48*> Cray-2. It could conceivably fail on hexadecimal or decimal machines 49*> without guard digits, but we know of none. 50*> \endverbatim 51* 52* Arguments: 53* ========== 54* 55*> \param[in] JOBZ 56*> \verbatim 57*> JOBZ is CHARACTER*1 58*> = 'N': Compute eigenvalues only; 59*> = 'V': Compute eigenvalues and eigenvectors. 60*> \endverbatim 61*> 62*> \param[in] UPLO 63*> \verbatim 64*> UPLO is CHARACTER*1 65*> = 'U': Upper triangle of A is stored; 66*> = 'L': Lower triangle of A is stored. 67*> \endverbatim 68*> 69*> \param[in] N 70*> \verbatim 71*> N is INTEGER 72*> The order of the matrix A. N >= 0. 73*> \endverbatim 74*> 75*> \param[in] KD 76*> \verbatim 77*> KD is INTEGER 78*> The number of superdiagonals of the matrix A if UPLO = 'U', 79*> or the number of subdiagonals if UPLO = 'L'. KD >= 0. 80*> \endverbatim 81*> 82*> \param[in,out] AB 83*> \verbatim 84*> AB is COMPLEX array, dimension (LDAB, N) 85*> On entry, the upper or lower triangle of the Hermitian band 86*> matrix A, stored in the first KD+1 rows of the array. The 87*> j-th column of A is stored in the j-th column of the array AB 88*> as follows: 89*> if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j; 90*> if UPLO = 'L', AB(1+i-j,j) = A(i,j) for j<=i<=min(n,j+kd). 91*> 92*> On exit, AB is overwritten by values generated during the 93*> reduction to tridiagonal form. If UPLO = 'U', the first 94*> superdiagonal and the diagonal of the tridiagonal matrix T 95*> are returned in rows KD and KD+1 of AB, and if UPLO = 'L', 96*> the diagonal and first subdiagonal of T are returned in the 97*> first two rows of AB. 98*> \endverbatim 99*> 100*> \param[in] LDAB 101*> \verbatim 102*> LDAB is INTEGER 103*> The leading dimension of the array AB. LDAB >= KD + 1. 104*> \endverbatim 105*> 106*> \param[out] W 107*> \verbatim 108*> W is REAL array, dimension (N) 109*> If INFO = 0, the eigenvalues in ascending order. 110*> \endverbatim 111*> 112*> \param[out] Z 113*> \verbatim 114*> Z is COMPLEX array, dimension (LDZ, N) 115*> If JOBZ = 'V', then if INFO = 0, Z contains the orthonormal 116*> eigenvectors of the matrix A, with the i-th column of Z 117*> holding the eigenvector associated with W(i). 118*> If JOBZ = 'N', then Z is not referenced. 119*> \endverbatim 120*> 121*> \param[in] LDZ 122*> \verbatim 123*> LDZ is INTEGER 124*> The leading dimension of the array Z. LDZ >= 1, and if 125*> JOBZ = 'V', LDZ >= max(1,N). 126*> \endverbatim 127*> 128*> \param[out] WORK 129*> \verbatim 130*> WORK is COMPLEX array, dimension (MAX(1,LWORK)) 131*> On exit, if INFO = 0, WORK(1) returns the optimal LWORK. 132*> \endverbatim 133*> 134*> \param[in] LWORK 135*> \verbatim 136*> LWORK is INTEGER 137*> The dimension of the array WORK. 138*> If N <= 1, LWORK must be at least 1. 139*> If JOBZ = 'N' and N > 1, LWORK must be at least N. 140*> If JOBZ = 'V' and N > 1, LWORK must be at least 2*N**2. 141*> 142*> If LWORK = -1, then a workspace query is assumed; the routine 143*> only calculates the optimal sizes of the WORK, RWORK and 144*> IWORK arrays, returns these values as the first entries of 145*> the WORK, RWORK and IWORK arrays, and no error message 146*> related to LWORK or LRWORK or LIWORK is issued by XERBLA. 147*> \endverbatim 148*> 149*> \param[out] RWORK 150*> \verbatim 151*> RWORK is REAL array, 152*> dimension (LRWORK) 153*> On exit, if INFO = 0, RWORK(1) returns the optimal LRWORK. 154*> \endverbatim 155*> 156*> \param[in] LRWORK 157*> \verbatim 158*> LRWORK is INTEGER 159*> The dimension of array RWORK. 160*> If N <= 1, LRWORK must be at least 1. 161*> If JOBZ = 'N' and N > 1, LRWORK must be at least N. 162*> If JOBZ = 'V' and N > 1, LRWORK must be at least 163*> 1 + 5*N + 2*N**2. 164*> 165*> If LRWORK = -1, then a workspace query is assumed; the 166*> routine only calculates the optimal sizes of the WORK, RWORK 167*> and IWORK arrays, returns these values as the first entries 168*> of the WORK, RWORK and IWORK arrays, and no error message 169*> related to LWORK or LRWORK or LIWORK is issued by XERBLA. 170*> \endverbatim 171*> 172*> \param[out] IWORK 173*> \verbatim 174*> IWORK is INTEGER array, dimension (MAX(1,LIWORK)) 175*> On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK. 176*> \endverbatim 177*> 178*> \param[in] LIWORK 179*> \verbatim 180*> LIWORK is INTEGER 181*> The dimension of array IWORK. 182*> If JOBZ = 'N' or N <= 1, LIWORK must be at least 1. 183*> If JOBZ = 'V' and N > 1, LIWORK must be at least 3 + 5*N . 184*> 185*> If LIWORK = -1, then a workspace query is assumed; the 186*> routine only calculates the optimal sizes of the WORK, RWORK 187*> and IWORK arrays, returns these values as the first entries 188*> of the WORK, RWORK and IWORK arrays, and no error message 189*> related to LWORK or LRWORK or LIWORK is issued by XERBLA. 190*> \endverbatim 191*> 192*> \param[out] INFO 193*> \verbatim 194*> INFO is INTEGER 195*> = 0: successful exit. 196*> < 0: if INFO = -i, the i-th argument had an illegal value. 197*> > 0: if INFO = i, the algorithm failed to converge; i 198*> off-diagonal elements of an intermediate tridiagonal 199*> form did not converge to zero. 200*> \endverbatim 201* 202* Authors: 203* ======== 204* 205*> \author Univ. of Tennessee 206*> \author Univ. of California Berkeley 207*> \author Univ. of Colorado Denver 208*> \author NAG Ltd. 209* 210*> \date December 2016 211* 212*> \ingroup complexOTHEReigen 213* 214* ===================================================================== 215 SUBROUTINE CHBEVD( JOBZ, UPLO, N, KD, AB, LDAB, W, Z, LDZ, WORK, 216 $ LWORK, RWORK, LRWORK, IWORK, LIWORK, INFO ) 217* 218* -- LAPACK driver routine (version 3.7.0) -- 219* -- LAPACK is a software package provided by Univ. of Tennessee, -- 220* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 221* December 2016 222* 223* .. Scalar Arguments .. 224 CHARACTER JOBZ, UPLO 225 INTEGER INFO, KD, LDAB, LDZ, LIWORK, LRWORK, LWORK, N 226* .. 227* .. Array Arguments .. 228 INTEGER IWORK( * ) 229 REAL RWORK( * ), W( * ) 230 COMPLEX AB( LDAB, * ), WORK( * ), Z( LDZ, * ) 231* .. 232* 233* ===================================================================== 234* 235* .. Parameters .. 236 REAL ZERO, ONE 237 PARAMETER ( ZERO = 0.0E0, ONE = 1.0E0 ) 238 COMPLEX CZERO, CONE 239 PARAMETER ( CZERO = ( 0.0E0, 0.0E0 ), 240 $ CONE = ( 1.0E0, 0.0E0 ) ) 241* .. 242* .. Local Scalars .. 243 LOGICAL LOWER, LQUERY, WANTZ 244 INTEGER IINFO, IMAX, INDE, INDWK2, INDWRK, ISCALE, 245 $ LIWMIN, LLRWK, LLWK2, LRWMIN, LWMIN 246 REAL ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA, 247 $ SMLNUM 248* .. 249* .. External Functions .. 250 LOGICAL LSAME 251 REAL CLANHB, SLAMCH 252 EXTERNAL LSAME, CLANHB, SLAMCH 253* .. 254* .. External Subroutines .. 255 EXTERNAL CGEMM, CHBTRD, CLACPY, CLASCL, CSTEDC, SSCAL, 256 $ SSTERF, XERBLA 257* .. 258* .. Intrinsic Functions .. 259 INTRINSIC SQRT 260* .. 261* .. Executable Statements .. 262* 263* Test the input parameters. 264* 265 WANTZ = LSAME( JOBZ, 'V' ) 266 LOWER = LSAME( UPLO, 'L' ) 267 LQUERY = ( LWORK.EQ.-1 .OR. LIWORK.EQ.-1 .OR. LRWORK.EQ.-1 ) 268* 269 INFO = 0 270 IF( N.LE.1 ) THEN 271 LWMIN = 1 272 LRWMIN = 1 273 LIWMIN = 1 274 ELSE 275 IF( WANTZ ) THEN 276 LWMIN = 2*N**2 277 LRWMIN = 1 + 5*N + 2*N**2 278 LIWMIN = 3 + 5*N 279 ELSE 280 LWMIN = N 281 LRWMIN = N 282 LIWMIN = 1 283 END IF 284 END IF 285 IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN 286 INFO = -1 287 ELSE IF( .NOT.( LOWER .OR. LSAME( UPLO, 'U' ) ) ) THEN 288 INFO = -2 289 ELSE IF( N.LT.0 ) THEN 290 INFO = -3 291 ELSE IF( KD.LT.0 ) THEN 292 INFO = -4 293 ELSE IF( LDAB.LT.KD+1 ) THEN 294 INFO = -6 295 ELSE IF( LDZ.LT.1 .OR. ( WANTZ .AND. LDZ.LT.N ) ) THEN 296 INFO = -9 297 END IF 298* 299 IF( INFO.EQ.0 ) THEN 300 WORK( 1 ) = LWMIN 301 RWORK( 1 ) = LRWMIN 302 IWORK( 1 ) = LIWMIN 303* 304 IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN 305 INFO = -11 306 ELSE IF( LRWORK.LT.LRWMIN .AND. .NOT.LQUERY ) THEN 307 INFO = -13 308 ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT.LQUERY ) THEN 309 INFO = -15 310 END IF 311 END IF 312* 313 IF( INFO.NE.0 ) THEN 314 CALL XERBLA( 'CHBEVD', -INFO ) 315 RETURN 316 ELSE IF( LQUERY ) THEN 317 RETURN 318 END IF 319* 320* Quick return if possible 321* 322 IF( N.EQ.0 ) 323 $ RETURN 324* 325 IF( N.EQ.1 ) THEN 326 W( 1 ) = AB( 1, 1 ) 327 IF( WANTZ ) 328 $ Z( 1, 1 ) = CONE 329 RETURN 330 END IF 331* 332* Get machine constants. 333* 334 SAFMIN = SLAMCH( 'Safe minimum' ) 335 EPS = SLAMCH( 'Precision' ) 336 SMLNUM = SAFMIN / EPS 337 BIGNUM = ONE / SMLNUM 338 RMIN = SQRT( SMLNUM ) 339 RMAX = SQRT( BIGNUM ) 340* 341* Scale matrix to allowable range, if necessary. 342* 343 ANRM = CLANHB( 'M', UPLO, N, KD, AB, LDAB, RWORK ) 344 ISCALE = 0 345 IF( ANRM.GT.ZERO .AND. ANRM.LT.RMIN ) THEN 346 ISCALE = 1 347 SIGMA = RMIN / ANRM 348 ELSE IF( ANRM.GT.RMAX ) THEN 349 ISCALE = 1 350 SIGMA = RMAX / ANRM 351 END IF 352 IF( ISCALE.EQ.1 ) THEN 353 IF( LOWER ) THEN 354 CALL CLASCL( 'B', KD, KD, ONE, SIGMA, N, N, AB, LDAB, INFO ) 355 ELSE 356 CALL CLASCL( 'Q', KD, KD, ONE, SIGMA, N, N, AB, LDAB, INFO ) 357 END IF 358 END IF 359* 360* Call CHBTRD to reduce Hermitian band matrix to tridiagonal form. 361* 362 INDE = 1 363 INDWRK = INDE + N 364 INDWK2 = 1 + N*N 365 LLWK2 = LWORK - INDWK2 + 1 366 LLRWK = LRWORK - INDWRK + 1 367 CALL CHBTRD( JOBZ, UPLO, N, KD, AB, LDAB, W, RWORK( INDE ), Z, 368 $ LDZ, WORK, IINFO ) 369* 370* For eigenvalues only, call SSTERF. For eigenvectors, call CSTEDC. 371* 372 IF( .NOT.WANTZ ) THEN 373 CALL SSTERF( N, W, RWORK( INDE ), INFO ) 374 ELSE 375 CALL CSTEDC( 'I', N, W, RWORK( INDE ), WORK, N, WORK( INDWK2 ), 376 $ LLWK2, RWORK( INDWRK ), LLRWK, IWORK, LIWORK, 377 $ INFO ) 378 CALL CGEMM( 'N', 'N', N, N, N, CONE, Z, LDZ, WORK, N, CZERO, 379 $ WORK( INDWK2 ), N ) 380 CALL CLACPY( 'A', N, N, WORK( INDWK2 ), N, Z, LDZ ) 381 END IF 382* 383* If matrix was scaled, then rescale eigenvalues appropriately. 384* 385 IF( ISCALE.EQ.1 ) THEN 386 IF( INFO.EQ.0 ) THEN 387 IMAX = N 388 ELSE 389 IMAX = INFO - 1 390 END IF 391 CALL SSCAL( IMAX, ONE / SIGMA, W, 1 ) 392 END IF 393* 394 WORK( 1 ) = LWMIN 395 RWORK( 1 ) = LRWMIN 396 IWORK( 1 ) = LIWMIN 397 RETURN 398* 399* End of CHBEVD 400* 401 END 402