LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ chbevx()

subroutine chbevx ( character  jobz,
character  range,
character  uplo,
integer  n,
integer  kd,
complex, dimension( ldab, * )  ab,
integer  ldab,
complex, dimension( ldq, * )  q,
integer  ldq,
real  vl,
real  vu,
integer  il,
integer  iu,
real  abstol,
integer  m,
real, dimension( * )  w,
complex, dimension( ldz, * )  z,
integer  ldz,
complex, dimension( * )  work,
real, dimension( * )  rwork,
integer, dimension( * )  iwork,
integer, dimension( * )  ifail,
integer  info 
)

CHBEVX computes the eigenvalues and, optionally, the left and/or right eigenvectors for OTHER matrices

Download CHBEVX + dependencies [TGZ] [ZIP] [TXT]

Purpose:
 CHBEVX computes selected eigenvalues and, optionally, eigenvectors
 of a complex Hermitian band matrix A.  Eigenvalues and eigenvectors
 can be selected by specifying either a range of values or a range of
 indices for the desired eigenvalues.
Parameters
[in]JOBZ
          JOBZ is CHARACTER*1
          = 'N':  Compute eigenvalues only;
          = 'V':  Compute eigenvalues and eigenvectors.
[in]RANGE
          RANGE is CHARACTER*1
          = 'A': all eigenvalues will be found;
          = 'V': all eigenvalues in the half-open interval (VL,VU]
                 will be found;
          = 'I': the IL-th through IU-th eigenvalues will be found.
[in]UPLO
          UPLO is CHARACTER*1
          = 'U':  Upper triangle of A is stored;
          = 'L':  Lower triangle of A is stored.
[in]N
          N is INTEGER
          The order of the matrix A.  N >= 0.
[in]KD
          KD is INTEGER
          The number of superdiagonals of the matrix A if UPLO = 'U',
          or the number of subdiagonals if UPLO = 'L'.  KD >= 0.
[in,out]AB
          AB is COMPLEX array, dimension (LDAB, N)
          On entry, the upper or lower triangle of the Hermitian band
          matrix A, stored in the first KD+1 rows of the array.  The
          j-th column of A is stored in the j-th column of the array AB
          as follows:
          if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j;
          if UPLO = 'L', AB(1+i-j,j)    = A(i,j) for j<=i<=min(n,j+kd).

          On exit, AB is overwritten by values generated during the
          reduction to tridiagonal form.
[in]LDAB
          LDAB is INTEGER
          The leading dimension of the array AB.  LDAB >= KD + 1.
[out]Q
          Q is COMPLEX array, dimension (LDQ, N)
          If JOBZ = 'V', the N-by-N unitary matrix used in the
                          reduction to tridiagonal form.
          If JOBZ = 'N', the array Q is not referenced.
[in]LDQ
          LDQ is INTEGER
          The leading dimension of the array Q.  If JOBZ = 'V', then
          LDQ >= max(1,N).
[in]VL
          VL is REAL
          If RANGE='V', the lower bound of the interval to
          be searched for eigenvalues. VL < VU.
          Not referenced if RANGE = 'A' or 'I'.
[in]VU
          VU is REAL
          If RANGE='V', the upper bound of the interval to
          be searched for eigenvalues. VL < VU.
          Not referenced if RANGE = 'A' or 'I'.
[in]IL
          IL is INTEGER
          If RANGE='I', the index of the
          smallest eigenvalue to be returned.
          1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
          Not referenced if RANGE = 'A' or 'V'.
[in]IU
          IU is INTEGER
          If RANGE='I', the index of the
          largest eigenvalue to be returned.
          1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
          Not referenced if RANGE = 'A' or 'V'.
[in]ABSTOL
          ABSTOL is REAL
          The absolute error tolerance for the eigenvalues.
          An approximate eigenvalue is accepted as converged
          when it is determined to lie in an interval [a,b]
          of width less than or equal to

                  ABSTOL + EPS *   max( |a|,|b| ) ,

          where EPS is the machine precision.  If ABSTOL is less than
          or equal to zero, then  EPS*|T|  will be used in its place,
          where |T| is the 1-norm of the tridiagonal matrix obtained
          by reducing AB to tridiagonal form.

          Eigenvalues will be computed most accurately when ABSTOL is
          set to twice the underflow threshold 2*SLAMCH('S'), not zero.
          If this routine returns with INFO>0, indicating that some
          eigenvectors did not converge, try setting ABSTOL to
          2*SLAMCH('S').

          See "Computing Small Singular Values of Bidiagonal Matrices
          with Guaranteed High Relative Accuracy," by Demmel and
          Kahan, LAPACK Working Note #3.
[out]M
          M is INTEGER
          The total number of eigenvalues found.  0 <= M <= N.
          If RANGE = 'A', M = N, and if RANGE = 'I', M = IU-IL+1.
[out]W
          W is REAL array, dimension (N)
          The first M elements contain the selected eigenvalues in
          ascending order.
[out]Z
          Z is COMPLEX array, dimension (LDZ, max(1,M))
          If JOBZ = 'V', then if INFO = 0, the first M columns of Z
          contain the orthonormal eigenvectors of the matrix A
          corresponding to the selected eigenvalues, with the i-th
          column of Z holding the eigenvector associated with W(i).
          If an eigenvector fails to converge, then that column of Z
          contains the latest approximation to the eigenvector, and the
          index of the eigenvector is returned in IFAIL.
          If JOBZ = 'N', then Z is not referenced.
          Note: the user must ensure that at least max(1,M) columns are
          supplied in the array Z; if RANGE = 'V', the exact value of M
          is not known in advance and an upper bound must be used.
[in]LDZ
          LDZ is INTEGER
          The leading dimension of the array Z.  LDZ >= 1, and if
          JOBZ = 'V', LDZ >= max(1,N).
[out]WORK
          WORK is COMPLEX array, dimension (N)
[out]RWORK
          RWORK is REAL array, dimension (7*N)
[out]IWORK
          IWORK is INTEGER array, dimension (5*N)
[out]IFAIL
          IFAIL is INTEGER array, dimension (N)
          If JOBZ = 'V', then if INFO = 0, the first M elements of
          IFAIL are zero.  If INFO > 0, then IFAIL contains the
          indices of the eigenvectors that failed to converge.
          If JOBZ = 'N', then IFAIL is not referenced.
[out]INFO
          INFO is INTEGER
          = 0:  successful exit
          < 0:  if INFO = -i, the i-th argument had an illegal value
          > 0:  if INFO = i, then i eigenvectors failed to converge.
                Their indices are stored in array IFAIL.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 264 of file chbevx.f.

267*
268* -- LAPACK driver routine --
269* -- LAPACK is a software package provided by Univ. of Tennessee, --
270* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
271*
272* .. Scalar Arguments ..
273 CHARACTER JOBZ, RANGE, UPLO
274 INTEGER IL, INFO, IU, KD, LDAB, LDQ, LDZ, M, N
275 REAL ABSTOL, VL, VU
276* ..
277* .. Array Arguments ..
278 INTEGER IFAIL( * ), IWORK( * )
279 REAL RWORK( * ), W( * )
280 COMPLEX AB( LDAB, * ), Q( LDQ, * ), WORK( * ),
281 $ Z( LDZ, * )
282* ..
283*
284* =====================================================================
285*
286* .. Parameters ..
287 REAL ZERO, ONE
288 parameter( zero = 0.0e0, one = 1.0e0 )
289 COMPLEX CZERO, CONE
290 parameter( czero = ( 0.0e0, 0.0e0 ),
291 $ cone = ( 1.0e0, 0.0e0 ) )
292* ..
293* .. Local Scalars ..
294 LOGICAL ALLEIG, INDEIG, LOWER, TEST, VALEIG, WANTZ
295 CHARACTER ORDER
296 INTEGER I, IINFO, IMAX, INDD, INDE, INDEE, INDIBL,
297 $ INDISP, INDIWK, INDRWK, INDWRK, ISCALE, ITMP1,
298 $ J, JJ, NSPLIT
299 REAL ABSTLL, ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN,
300 $ SIGMA, SMLNUM, TMP1, VLL, VUU
301 COMPLEX CTMP1
302* ..
303* .. External Functions ..
304 LOGICAL LSAME
305 REAL CLANHB, SLAMCH
306 EXTERNAL lsame, clanhb, slamch
307* ..
308* .. External Subroutines ..
309 EXTERNAL ccopy, cgemv, chbtrd, clacpy, clascl, cstein,
311 $ xerbla
312* ..
313* .. Intrinsic Functions ..
314 INTRINSIC max, min, real, sqrt
315* ..
316* .. Executable Statements ..
317*
318* Test the input parameters.
319*
320 wantz = lsame( jobz, 'V' )
321 alleig = lsame( range, 'A' )
322 valeig = lsame( range, 'V' )
323 indeig = lsame( range, 'I' )
324 lower = lsame( uplo, 'L' )
325*
326 info = 0
327 IF( .NOT.( wantz .OR. lsame( jobz, 'N' ) ) ) THEN
328 info = -1
329 ELSE IF( .NOT.( alleig .OR. valeig .OR. indeig ) ) THEN
330 info = -2
331 ELSE IF( .NOT.( lower .OR. lsame( uplo, 'U' ) ) ) THEN
332 info = -3
333 ELSE IF( n.LT.0 ) THEN
334 info = -4
335 ELSE IF( kd.LT.0 ) THEN
336 info = -5
337 ELSE IF( ldab.LT.kd+1 ) THEN
338 info = -7
339 ELSE IF( wantz .AND. ldq.LT.max( 1, n ) ) THEN
340 info = -9
341 ELSE
342 IF( valeig ) THEN
343 IF( n.GT.0 .AND. vu.LE.vl )
344 $ info = -11
345 ELSE IF( indeig ) THEN
346 IF( il.LT.1 .OR. il.GT.max( 1, n ) ) THEN
347 info = -12
348 ELSE IF( iu.LT.min( n, il ) .OR. iu.GT.n ) THEN
349 info = -13
350 END IF
351 END IF
352 END IF
353 IF( info.EQ.0 ) THEN
354 IF( ldz.LT.1 .OR. ( wantz .AND. ldz.LT.n ) )
355 $ info = -18
356 END IF
357*
358 IF( info.NE.0 ) THEN
359 CALL xerbla( 'CHBEVX', -info )
360 RETURN
361 END IF
362*
363* Quick return if possible
364*
365 m = 0
366 IF( n.EQ.0 )
367 $ RETURN
368*
369 IF( n.EQ.1 ) THEN
370 m = 1
371 IF( lower ) THEN
372 ctmp1 = ab( 1, 1 )
373 ELSE
374 ctmp1 = ab( kd+1, 1 )
375 END IF
376 tmp1 = real( ctmp1 )
377 IF( valeig ) THEN
378 IF( .NOT.( vl.LT.tmp1 .AND. vu.GE.tmp1 ) )
379 $ m = 0
380 END IF
381 IF( m.EQ.1 ) THEN
382 w( 1 ) = real( ctmp1 )
383 IF( wantz )
384 $ z( 1, 1 ) = cone
385 END IF
386 RETURN
387 END IF
388*
389* Get machine constants.
390*
391 safmin = slamch( 'Safe minimum' )
392 eps = slamch( 'Precision' )
393 smlnum = safmin / eps
394 bignum = one / smlnum
395 rmin = sqrt( smlnum )
396 rmax = min( sqrt( bignum ), one / sqrt( sqrt( safmin ) ) )
397*
398* Scale matrix to allowable range, if necessary.
399*
400 iscale = 0
401 abstll = abstol
402 IF ( valeig ) THEN
403 vll = vl
404 vuu = vu
405 ELSE
406 vll = zero
407 vuu = zero
408 ENDIF
409 anrm = clanhb( 'M', uplo, n, kd, ab, ldab, rwork )
410 IF( anrm.GT.zero .AND. anrm.LT.rmin ) THEN
411 iscale = 1
412 sigma = rmin / anrm
413 ELSE IF( anrm.GT.rmax ) THEN
414 iscale = 1
415 sigma = rmax / anrm
416 END IF
417 IF( iscale.EQ.1 ) THEN
418 IF( lower ) THEN
419 CALL clascl( 'B', kd, kd, one, sigma, n, n, ab, ldab, info )
420 ELSE
421 CALL clascl( 'Q', kd, kd, one, sigma, n, n, ab, ldab, info )
422 END IF
423 IF( abstol.GT.0 )
424 $ abstll = abstol*sigma
425 IF( valeig ) THEN
426 vll = vl*sigma
427 vuu = vu*sigma
428 END IF
429 END IF
430*
431* Call CHBTRD to reduce Hermitian band matrix to tridiagonal form.
432*
433 indd = 1
434 inde = indd + n
435 indrwk = inde + n
436 indwrk = 1
437 CALL chbtrd( jobz, uplo, n, kd, ab, ldab, rwork( indd ),
438 $ rwork( inde ), q, ldq, work( indwrk ), iinfo )
439*
440* If all eigenvalues are desired and ABSTOL is less than or equal
441* to zero, then call SSTERF or CSTEQR. If this fails for some
442* eigenvalue, then try SSTEBZ.
443*
444 test = .false.
445 IF (indeig) THEN
446 IF (il.EQ.1 .AND. iu.EQ.n) THEN
447 test = .true.
448 END IF
449 END IF
450 IF ((alleig .OR. test) .AND. (abstol.LE.zero)) THEN
451 CALL scopy( n, rwork( indd ), 1, w, 1 )
452 indee = indrwk + 2*n
453 IF( .NOT.wantz ) THEN
454 CALL scopy( n-1, rwork( inde ), 1, rwork( indee ), 1 )
455 CALL ssterf( n, w, rwork( indee ), info )
456 ELSE
457 CALL clacpy( 'A', n, n, q, ldq, z, ldz )
458 CALL scopy( n-1, rwork( inde ), 1, rwork( indee ), 1 )
459 CALL csteqr( jobz, n, w, rwork( indee ), z, ldz,
460 $ rwork( indrwk ), info )
461 IF( info.EQ.0 ) THEN
462 DO 10 i = 1, n
463 ifail( i ) = 0
464 10 CONTINUE
465 END IF
466 END IF
467 IF( info.EQ.0 ) THEN
468 m = n
469 GO TO 30
470 END IF
471 info = 0
472 END IF
473*
474* Otherwise, call SSTEBZ and, if eigenvectors are desired, CSTEIN.
475*
476 IF( wantz ) THEN
477 order = 'B'
478 ELSE
479 order = 'E'
480 END IF
481 indibl = 1
482 indisp = indibl + n
483 indiwk = indisp + n
484 CALL sstebz( range, order, n, vll, vuu, il, iu, abstll,
485 $ rwork( indd ), rwork( inde ), m, nsplit, w,
486 $ iwork( indibl ), iwork( indisp ), rwork( indrwk ),
487 $ iwork( indiwk ), info )
488*
489 IF( wantz ) THEN
490 CALL cstein( n, rwork( indd ), rwork( inde ), m, w,
491 $ iwork( indibl ), iwork( indisp ), z, ldz,
492 $ rwork( indrwk ), iwork( indiwk ), ifail, info )
493*
494* Apply unitary matrix used in reduction to tridiagonal
495* form to eigenvectors returned by CSTEIN.
496*
497 DO 20 j = 1, m
498 CALL ccopy( n, z( 1, j ), 1, work( 1 ), 1 )
499 CALL cgemv( 'N', n, n, cone, q, ldq, work, 1, czero,
500 $ z( 1, j ), 1 )
501 20 CONTINUE
502 END IF
503*
504* If matrix was scaled, then rescale eigenvalues appropriately.
505*
506 30 CONTINUE
507 IF( iscale.EQ.1 ) THEN
508 IF( info.EQ.0 ) THEN
509 imax = m
510 ELSE
511 imax = info - 1
512 END IF
513 CALL sscal( imax, one / sigma, w, 1 )
514 END IF
515*
516* If eigenvalues are not in order, then sort them, along with
517* eigenvectors.
518*
519 IF( wantz ) THEN
520 DO 50 j = 1, m - 1
521 i = 0
522 tmp1 = w( j )
523 DO 40 jj = j + 1, m
524 IF( w( jj ).LT.tmp1 ) THEN
525 i = jj
526 tmp1 = w( jj )
527 END IF
528 40 CONTINUE
529*
530 IF( i.NE.0 ) THEN
531 itmp1 = iwork( indibl+i-1 )
532 w( i ) = w( j )
533 iwork( indibl+i-1 ) = iwork( indibl+j-1 )
534 w( j ) = tmp1
535 iwork( indibl+j-1 ) = itmp1
536 CALL cswap( n, z( 1, i ), 1, z( 1, j ), 1 )
537 IF( info.NE.0 ) THEN
538 itmp1 = ifail( i )
539 ifail( i ) = ifail( j )
540 ifail( j ) = itmp1
541 END IF
542 END IF
543 50 CONTINUE
544 END IF
545*
546 RETURN
547*
548* End of CHBEVX
549*
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine scopy(n, sx, incx, sy, incy)
SCOPY
Definition scopy.f:82
subroutine ccopy(n, cx, incx, cy, incy)
CCOPY
Definition ccopy.f:81
subroutine cgemv(trans, m, n, alpha, a, lda, x, incx, beta, y, incy)
CGEMV
Definition cgemv.f:160
subroutine chbtrd(vect, uplo, n, kd, ab, ldab, d, e, q, ldq, work, info)
CHBTRD
Definition chbtrd.f:163
subroutine clacpy(uplo, m, n, a, lda, b, ldb)
CLACPY copies all or part of one two-dimensional array to another.
Definition clacpy.f:103
real function slamch(cmach)
SLAMCH
Definition slamch.f:68
real function clanhb(norm, uplo, n, k, ab, ldab, work)
CLANHB returns the value of the 1-norm, or the Frobenius norm, or the infinity norm,...
Definition clanhb.f:132
subroutine clascl(type, kl, ku, cfrom, cto, m, n, a, lda, info)
CLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
Definition clascl.f:143
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
subroutine sscal(n, sa, sx, incx)
SSCAL
Definition sscal.f:79
subroutine sstebz(range, order, n, vl, vu, il, iu, abstol, d, e, m, nsplit, w, iblock, isplit, work, iwork, info)
SSTEBZ
Definition sstebz.f:273
subroutine cstein(n, d, e, m, w, iblock, isplit, z, ldz, work, iwork, ifail, info)
CSTEIN
Definition cstein.f:182
subroutine csteqr(compz, n, d, e, z, ldz, work, info)
CSTEQR
Definition csteqr.f:132
subroutine ssterf(n, d, e, info)
SSTERF
Definition ssterf.f:86
subroutine cswap(n, cx, incx, cy, incy)
CSWAP
Definition cswap.f:81
Here is the call graph for this function:
Here is the caller graph for this function: