LAPACK  3.10.1
LAPACK: Linear Algebra PACKage

◆ 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)
XERBLA
Definition: xerbla.f:60
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
subroutine ssterf(N, D, E, INFO)
SSTERF
Definition: ssterf.f:86
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 ccopy(N, CX, INCX, CY, INCY)
CCOPY
Definition: ccopy.f:81
subroutine cswap(N, CX, INCX, CY, INCY)
CSWAP
Definition: cswap.f:81
subroutine cgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
CGEMV
Definition: cgemv.f:158
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
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 clacpy(UPLO, M, N, A, LDA, B, LDB)
CLACPY copies all or part of one two-dimensional array to another.
Definition: clacpy.f:103
subroutine cstein(N, D, E, M, W, IBLOCK, ISPLIT, Z, LDZ, WORK, IWORK, IFAIL, INFO)
CSTEIN
Definition: cstein.f:182
subroutine chbtrd(VECT, UPLO, N, KD, AB, LDAB, D, E, Q, LDQ, WORK, INFO)
CHBTRD
Definition: chbtrd.f:163
subroutine csteqr(COMPZ, N, D, E, Z, LDZ, WORK, INFO)
CSTEQR
Definition: csteqr.f:132
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY
Definition: scopy.f:82
subroutine sscal(N, SA, SX, INCX)
SSCAL
Definition: sscal.f:79
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:68
Here is the call graph for this function:
Here is the caller graph for this function: