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

Go to the source code of this file.

Functions/Subroutines

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

Function/Subroutine Documentation

subroutine zhbevx ( character  JOBZ,
character  RANGE,
character  UPLO,
integer  N,
integer  KD,
complex*16, dimension( ldab, * )  AB,
integer  LDAB,
complex*16, dimension( ldq, * )  Q,
integer  LDQ,
double precision  VL,
double precision  VU,
integer  IL,
integer  IU,
double precision  ABSTOL,
integer  M,
double precision, dimension( * )  W,
complex*16, dimension( ldz, * )  Z,
integer  LDZ,
complex*16, dimension( * )  WORK,
double precision, dimension( * )  RWORK,
integer, dimension( * )  IWORK,
integer, dimension( * )  IFAIL,
integer  INFO 
)

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

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

Purpose:
 ZHBEVX 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*16 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*16 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 DOUBLE PRECISION
[in]VU
          VU is DOUBLE PRECISION
          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 DOUBLE PRECISION
          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*DLAMCH('S'), not zero.
          If this routine returns with INFO>0, indicating that some
          eigenvectors did not converge, try setting ABSTOL to
          2*DLAMCH('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 DOUBLE PRECISION array, dimension (N)
          The first M elements contain the selected eigenvalues in
          ascending order.
[out]Z
          Z is COMPLEX*16 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*16 array, dimension (N)
[out]RWORK
          RWORK is DOUBLE PRECISION 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 zhbevx.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  DOUBLE PRECISION abstol, vl, vu
272 * ..
273 * .. Array Arguments ..
274  INTEGER ifail( * ), iwork( * )
275  DOUBLE PRECISION rwork( * ), w( * )
276  COMPLEX*16 ab( ldab, * ), q( ldq, * ), work( * ),
277  $ z( ldz, * )
278 * ..
279 *
280 * =====================================================================
281 *
282 * .. Parameters ..
283  DOUBLE PRECISION zero, one
284  parameter( zero = 0.0d0, one = 1.0d0 )
285  COMPLEX*16 czero, cone
286  parameter( czero = ( 0.0d0, 0.0d0 ),
287  $ cone = ( 1.0d0, 0.0d0 ) )
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  DOUBLE PRECISION abstll, anrm, bignum, eps, rmax, rmin, safmin,
296  $ sigma, smlnum, tmp1, vll, vuu
297  COMPLEX*16 ctmp1
298 * ..
299 * .. External Functions ..
300  LOGICAL lsame
301  DOUBLE PRECISION dlamch, zlanhb
302  EXTERNAL lsame, dlamch, zlanhb
303 * ..
304 * .. External Subroutines ..
305  EXTERNAL dcopy, dscal, dstebz, dsterf, xerbla, zcopy,
307  $ zswap
308 * ..
309 * .. Intrinsic Functions ..
310  INTRINSIC dble, max, min, 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( 'ZHBEVX', -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 = dble( 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 = dlamch( 'Safe minimum' )
388  eps = dlamch( '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  END IF
405  anrm = zlanhb( '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 zlascl( 'B', kd, kd, one, sigma, n, n, ab, ldab, info )
416  ELSE
417  CALL zlascl( '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 ZHBTRD 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 zhbtrd( 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 DSTERF or ZSTEQR. If this fails for some
438 * eigenvalue, then try DSTEBZ.
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 dcopy( n, rwork( indd ), 1, w, 1 )
448  indee = indrwk + 2*n
449  IF( .NOT.wantz ) THEN
450  CALL dcopy( n-1, rwork( inde ), 1, rwork( indee ), 1 )
451  CALL dsterf( n, w, rwork( indee ), info )
452  ELSE
453  CALL zlacpy( 'A', n, n, q, ldq, z, ldz )
454  CALL dcopy( n-1, rwork( inde ), 1, rwork( indee ), 1 )
455  CALL zsteqr( 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 DSTEBZ and, if eigenvectors are desired, ZSTEIN.
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 dstebz( 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 zstein( 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 ZSTEIN.
492 *
493  DO 20 j = 1, m
494  CALL zcopy( n, z( 1, j ), 1, work( 1 ), 1 )
495  CALL zgemv( '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 dscal( 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 zswap( 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 ZHBEVX
545 *
subroutine zswap(N, ZX, INCX, ZY, INCY)
ZSWAP
Definition: zswap.f:52
subroutine zgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
ZGEMV
Definition: zgemv.f:160
subroutine zlacpy(UPLO, M, N, A, LDA, B, LDB)
ZLACPY copies all or part of one two-dimensional array to another.
Definition: zlacpy.f:105
subroutine zstein(N, D, E, M, W, IBLOCK, ISPLIT, Z, LDZ, WORK, IWORK, IFAIL, INFO)
ZSTEIN
Definition: zstein.f:184
subroutine dsterf(N, D, E, INFO)
DSTERF
Definition: dsterf.f:88
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
Definition: dcopy.f:53
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine zhbtrd(VECT, UPLO, N, KD, AB, LDAB, D, E, Q, LDQ, WORK, INFO)
ZHBTRD
Definition: zhbtrd.f:165
subroutine zlascl(TYPE, KL, KU, CFROM, CTO, M, N, A, LDA, INFO)
ZLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
Definition: zlascl.f:141
subroutine zcopy(N, ZX, INCX, ZY, INCY)
ZCOPY
Definition: zcopy.f:52
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
subroutine dscal(N, DA, DX, INCX)
DSCAL
Definition: dscal.f:55
double precision function zlanhb(NORM, UPLO, N, K, AB, LDAB, WORK)
ZLANHB 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: zlanhb.f:134
subroutine zsteqr(COMPZ, N, D, E, Z, LDZ, WORK, INFO)
ZSTEQR
Definition: zsteqr.f:134
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:65
subroutine dstebz(RANGE, ORDER, N, VL, VU, IL, IU, ABSTOL, D, E, M, NSPLIT, W, IBLOCK, ISPLIT, WORK, IWORK, INFO)
DSTEBZ
Definition: dstebz.f:265

Here is the call graph for this function:

Here is the caller graph for this function: