LAPACK  3.5.0
LAPACK: Linear Algebra PACKage
chbevx.f File Reference

Go to the source code of this file.

Functions/Subroutines

subroutine chbevx (JOBZ, RANGE, UPLO, N, KD, AB, LDAB, Q, LDQ, VL, VU, IL, IU, ABSTOL, M, W, Z, LDZ, WORK, RWORK, IWORK, IFAIL, INFO)
  CHBEVX computes the eigenvalues and, optionally, the left and/or right eigenvectors for OTHER matrices More...
 

Function/Subroutine Documentation

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
[in]VU
          VU is REAL
          If RANGE='V', the lower and upper bounds of the interval to
          be searched for eigenvalues. VL < VU.
          Not referenced if RANGE = 'A' or 'I'.
[in]IL
          IL is INTEGER
[in]IU
          IU is INTEGER
          If RANGE='I', the indices (in ascending order) of the
          smallest and largest eigenvalues 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.
Date
November 2011

Definition at line 262 of file chbevx.f.

262 *
263 * -- LAPACK driver routine (version 3.4.0) --
264 * -- LAPACK is a software package provided by Univ. of Tennessee, --
265 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
266 * November 2011
267 *
268 * .. Scalar Arguments ..
269  CHARACTER jobz, range, uplo
270  INTEGER il, info, iu, kd, ldab, ldq, ldz, m, n
271  REAL abstol, vl, vu
272 * ..
273 * .. Array Arguments ..
274  INTEGER ifail( * ), iwork( * )
275  REAL rwork( * ), w( * )
276  COMPLEX ab( ldab, * ), q( ldq, * ), work( * ),
277  $ z( ldz, * )
278 * ..
279 *
280 * =====================================================================
281 *
282 * .. Parameters ..
283  REAL zero, one
284  parameter( zero = 0.0e0, one = 1.0e0 )
285  COMPLEX czero, cone
286  parameter( czero = ( 0.0e0, 0.0e0 ),
287  $ cone = ( 1.0e0, 0.0e0 ) )
288 * ..
289 * .. Local Scalars ..
290  LOGICAL alleig, indeig, lower, test, valeig, wantz
291  CHARACTER order
292  INTEGER i, iinfo, imax, indd, inde, indee, indibl,
293  $ indisp, indiwk, indrwk, indwrk, iscale, itmp1,
294  $ j, jj, nsplit
295  REAL abstll, anrm, bignum, eps, rmax, rmin, safmin,
296  $ sigma, smlnum, tmp1, vll, vuu
297  COMPLEX ctmp1
298 * ..
299 * .. External Functions ..
300  LOGICAL lsame
301  REAL clanhb, slamch
302  EXTERNAL lsame, clanhb, slamch
303 * ..
304 * .. External Subroutines ..
305  EXTERNAL ccopy, cgemv, chbtrd, clacpy, clascl, cstein,
307  $ xerbla
308 * ..
309 * .. Intrinsic Functions ..
310  INTRINSIC max, min, REAL, sqrt
311 * ..
312 * .. Executable Statements ..
313 *
314 * Test the input parameters.
315 *
316  wantz = lsame( jobz, 'V' )
317  alleig = lsame( range, 'A' )
318  valeig = lsame( range, 'V' )
319  indeig = lsame( range, 'I' )
320  lower = lsame( uplo, 'L' )
321 *
322  info = 0
323  IF( .NOT.( wantz .OR. lsame( jobz, 'N' ) ) ) THEN
324  info = -1
325  ELSE IF( .NOT.( alleig .OR. valeig .OR. indeig ) ) THEN
326  info = -2
327  ELSE IF( .NOT.( lower .OR. lsame( uplo, 'U' ) ) ) THEN
328  info = -3
329  ELSE IF( n.LT.0 ) THEN
330  info = -4
331  ELSE IF( kd.LT.0 ) THEN
332  info = -5
333  ELSE IF( ldab.LT.kd+1 ) THEN
334  info = -7
335  ELSE IF( wantz .AND. ldq.LT.max( 1, n ) ) THEN
336  info = -9
337  ELSE
338  IF( valeig ) THEN
339  IF( n.GT.0 .AND. vu.LE.vl )
340  $ info = -11
341  ELSE IF( indeig ) THEN
342  IF( il.LT.1 .OR. il.GT.max( 1, n ) ) THEN
343  info = -12
344  ELSE IF( iu.LT.min( n, il ) .OR. iu.GT.n ) THEN
345  info = -13
346  END IF
347  END IF
348  END IF
349  IF( info.EQ.0 ) THEN
350  IF( ldz.LT.1 .OR. ( wantz .AND. ldz.LT.n ) )
351  $ info = -18
352  END IF
353 *
354  IF( info.NE.0 ) THEN
355  CALL xerbla( 'CHBEVX', -info )
356  RETURN
357  END IF
358 *
359 * Quick return if possible
360 *
361  m = 0
362  IF( n.EQ.0 )
363  $ RETURN
364 *
365  IF( n.EQ.1 ) THEN
366  m = 1
367  IF( lower ) THEN
368  ctmp1 = ab( 1, 1 )
369  ELSE
370  ctmp1 = ab( kd+1, 1 )
371  END IF
372  tmp1 = REAL( ctmp1 )
373  IF( valeig ) THEN
374  IF( .NOT.( vl.LT.tmp1 .AND. vu.GE.tmp1 ) )
375  $ m = 0
376  END IF
377  IF( m.EQ.1 ) THEN
378  w( 1 ) = ctmp1
379  IF( wantz )
380  $ z( 1, 1 ) = cone
381  END IF
382  RETURN
383  END IF
384 *
385 * Get machine constants.
386 *
387  safmin = slamch( 'Safe minimum' )
388  eps = slamch( 'Precision' )
389  smlnum = safmin / eps
390  bignum = one / smlnum
391  rmin = sqrt( smlnum )
392  rmax = min( sqrt( bignum ), one / sqrt( sqrt( safmin ) ) )
393 *
394 * Scale matrix to allowable range, if necessary.
395 *
396  iscale = 0
397  abstll = abstol
398  IF ( valeig ) THEN
399  vll = vl
400  vuu = vu
401  ELSE
402  vll = zero
403  vuu = zero
404  ENDIF
405  anrm = clanhb( 'M', uplo, n, kd, ab, ldab, rwork )
406  IF( anrm.GT.zero .AND. anrm.LT.rmin ) THEN
407  iscale = 1
408  sigma = rmin / anrm
409  ELSE IF( anrm.GT.rmax ) THEN
410  iscale = 1
411  sigma = rmax / anrm
412  END IF
413  IF( iscale.EQ.1 ) THEN
414  IF( lower ) THEN
415  CALL clascl( 'B', kd, kd, one, sigma, n, n, ab, ldab, info )
416  ELSE
417  CALL clascl( 'Q', kd, kd, one, sigma, n, n, ab, ldab, info )
418  END IF
419  IF( abstol.GT.0 )
420  $ abstll = abstol*sigma
421  IF( valeig ) THEN
422  vll = vl*sigma
423  vuu = vu*sigma
424  END IF
425  END IF
426 *
427 * Call CHBTRD to reduce Hermitian band matrix to tridiagonal form.
428 *
429  indd = 1
430  inde = indd + n
431  indrwk = inde + n
432  indwrk = 1
433  CALL chbtrd( jobz, uplo, n, kd, ab, ldab, rwork( indd ),
434  $ rwork( inde ), q, ldq, work( indwrk ), iinfo )
435 *
436 * If all eigenvalues are desired and ABSTOL is less than or equal
437 * to zero, then call SSTERF or CSTEQR. If this fails for some
438 * eigenvalue, then try SSTEBZ.
439 *
440  test = .false.
441  IF (indeig) THEN
442  IF (il.EQ.1 .AND. iu.EQ.n) THEN
443  test = .true.
444  END IF
445  END IF
446  IF ((alleig .OR. test) .AND. (abstol.LE.zero)) THEN
447  CALL scopy( n, rwork( indd ), 1, w, 1 )
448  indee = indrwk + 2*n
449  IF( .NOT.wantz ) THEN
450  CALL scopy( n-1, rwork( inde ), 1, rwork( indee ), 1 )
451  CALL ssterf( n, w, rwork( indee ), info )
452  ELSE
453  CALL clacpy( 'A', n, n, q, ldq, z, ldz )
454  CALL scopy( n-1, rwork( inde ), 1, rwork( indee ), 1 )
455  CALL csteqr( jobz, n, w, rwork( indee ), z, ldz,
456  $ rwork( indrwk ), info )
457  IF( info.EQ.0 ) THEN
458  DO 10 i = 1, n
459  ifail( i ) = 0
460  10 CONTINUE
461  END IF
462  END IF
463  IF( info.EQ.0 ) THEN
464  m = n
465  GO TO 30
466  END IF
467  info = 0
468  END IF
469 *
470 * Otherwise, call SSTEBZ and, if eigenvectors are desired, CSTEIN.
471 *
472  IF( wantz ) THEN
473  order = 'B'
474  ELSE
475  order = 'E'
476  END IF
477  indibl = 1
478  indisp = indibl + n
479  indiwk = indisp + n
480  CALL sstebz( range, order, n, vll, vuu, il, iu, abstll,
481  $ rwork( indd ), rwork( inde ), m, nsplit, w,
482  $ iwork( indibl ), iwork( indisp ), rwork( indrwk ),
483  $ iwork( indiwk ), info )
484 *
485  IF( wantz ) THEN
486  CALL cstein( n, rwork( indd ), rwork( inde ), m, w,
487  $ iwork( indibl ), iwork( indisp ), z, ldz,
488  $ rwork( indrwk ), iwork( indiwk ), ifail, info )
489 *
490 * Apply unitary matrix used in reduction to tridiagonal
491 * form to eigenvectors returned by CSTEIN.
492 *
493  DO 20 j = 1, m
494  CALL ccopy( n, z( 1, j ), 1, work( 1 ), 1 )
495  CALL cgemv( 'N', n, n, cone, q, ldq, work, 1, czero,
496  $ z( 1, j ), 1 )
497  20 CONTINUE
498  END IF
499 *
500 * If matrix was scaled, then rescale eigenvalues appropriately.
501 *
502  30 CONTINUE
503  IF( iscale.EQ.1 ) THEN
504  IF( info.EQ.0 ) THEN
505  imax = m
506  ELSE
507  imax = info - 1
508  END IF
509  CALL sscal( imax, one / sigma, w, 1 )
510  END IF
511 *
512 * If eigenvalues are not in order, then sort them, along with
513 * eigenvectors.
514 *
515  IF( wantz ) THEN
516  DO 50 j = 1, m - 1
517  i = 0
518  tmp1 = w( j )
519  DO 40 jj = j + 1, m
520  IF( w( jj ).LT.tmp1 ) THEN
521  i = jj
522  tmp1 = w( jj )
523  END IF
524  40 CONTINUE
525 *
526  IF( i.NE.0 ) THEN
527  itmp1 = iwork( indibl+i-1 )
528  w( i ) = w( j )
529  iwork( indibl+i-1 ) = iwork( indibl+j-1 )
530  w( j ) = tmp1
531  iwork( indibl+j-1 ) = itmp1
532  CALL cswap( n, z( 1, i ), 1, z( 1, j ), 1 )
533  IF( info.NE.0 ) THEN
534  itmp1 = ifail( i )
535  ifail( i ) = ifail( j )
536  ifail( j ) = itmp1
537  END IF
538  END IF
539  50 CONTINUE
540  END IF
541 *
542  RETURN
543 *
544 * End of CHBEVX
545 *
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine chbtrd(VECT, UPLO, N, KD, AB, LDAB, D, E, Q, LDQ, WORK, INFO)
CHBTRD
Definition: chbtrd.f:165
subroutine cswap(N, CX, INCX, CY, INCY)
CSWAP
Definition: cswap.f:52
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY
Definition: scopy.f:53
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
subroutine cgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
CGEMV
Definition: cgemv.f:160
subroutine clacpy(UPLO, M, N, A, LDA, B, LDB)
CLACPY copies all or part of one two-dimensional array to another.
Definition: clacpy.f:105
subroutine csteqr(COMPZ, N, D, E, Z, LDZ, WORK, INFO)
CSTEQR
Definition: csteqr.f:134
subroutine cstein(N, D, E, M, W, IBLOCK, ISPLIT, Z, LDZ, WORK, IWORK, IFAIL, INFO)
CSTEIN
Definition: cstein.f:184
subroutine sstebz(RANGE, ORDER, N, VL, VU, IL, IU, ABSTOL, D, E, M, NSPLIT, W, IBLOCK, ISPLIT, WORK, IWORK, INFO)
SSTEBZ
Definition: sstebz.f:265
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:69
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, or the element of largest absolute value of a Hermitian band matrix.
Definition: clanhb.f:134
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:141
subroutine ccopy(N, CX, INCX, CY, INCY)
CCOPY
Definition: ccopy.f:52
subroutine ssterf(N, D, E, INFO)
SSTERF
Definition: ssterf.f:88
subroutine sscal(N, SA, SX, INCX)
SSCAL
Definition: sscal.f:55

Here is the call graph for this function:

Here is the caller graph for this function: