LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

◆ zhbevx()

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
          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 DOUBLE PRECISION
          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 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.

Definition at line 264 of file zhbevx.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  DOUBLE PRECISION ABSTOL, VL, VU
276 * ..
277 * .. Array Arguments ..
278  INTEGER IFAIL( * ), IWORK( * )
279  DOUBLE PRECISION RWORK( * ), W( * )
280  COMPLEX*16 AB( LDAB, * ), Q( LDQ, * ), WORK( * ),
281  $ Z( LDZ, * )
282 * ..
283 *
284 * =====================================================================
285 *
286 * .. Parameters ..
287  DOUBLE PRECISION ZERO, ONE
288  parameter( zero = 0.0d0, one = 1.0d0 )
289  COMPLEX*16 CZERO, CONE
290  parameter( czero = ( 0.0d0, 0.0d0 ),
291  $ cone = ( 1.0d0, 0.0d0 ) )
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  DOUBLE PRECISION ABSTLL, ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN,
300  $ SIGMA, SMLNUM, TMP1, VLL, VUU
301  COMPLEX*16 CTMP1
302 * ..
303 * .. External Functions ..
304  LOGICAL LSAME
305  DOUBLE PRECISION DLAMCH, ZLANHB
306  EXTERNAL lsame, dlamch, zlanhb
307 * ..
308 * .. External Subroutines ..
309  EXTERNAL dcopy, dscal, dstebz, dsterf, xerbla, zcopy,
311  $ zswap
312 * ..
313 * .. Intrinsic Functions ..
314  INTRINSIC dble, max, min, 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( 'ZHBEVX', -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 = dble( 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 ) = dble( 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 = dlamch( 'Safe minimum' )
392  eps = dlamch( '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  END IF
409  anrm = zlanhb( '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 zlascl( 'B', kd, kd, one, sigma, n, n, ab, ldab, info )
420  ELSE
421  CALL zlascl( '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 ZHBTRD 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 zhbtrd( 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 DSTERF or ZSTEQR. If this fails for some
442 * eigenvalue, then try DSTEBZ.
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 dcopy( n, rwork( indd ), 1, w, 1 )
452  indee = indrwk + 2*n
453  IF( .NOT.wantz ) THEN
454  CALL dcopy( n-1, rwork( inde ), 1, rwork( indee ), 1 )
455  CALL dsterf( n, w, rwork( indee ), info )
456  ELSE
457  CALL zlacpy( 'A', n, n, q, ldq, z, ldz )
458  CALL dcopy( n-1, rwork( inde ), 1, rwork( indee ), 1 )
459  CALL zsteqr( 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 DSTEBZ and, if eigenvectors are desired, ZSTEIN.
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 dstebz( 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 zstein( 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 ZSTEIN.
496 *
497  DO 20 j = 1, m
498  CALL zcopy( n, z( 1, j ), 1, work( 1 ), 1 )
499  CALL zgemv( '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 dscal( 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 zswap( 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 ZHBEVX
549 *
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:69
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
subroutine dstebz(RANGE, ORDER, N, VL, VU, IL, IU, ABSTOL, D, E, M, NSPLIT, W, IBLOCK, ISPLIT, WORK, IWORK, INFO)
DSTEBZ
Definition: dstebz.f:273
subroutine dsterf(N, D, E, INFO)
DSTERF
Definition: dsterf.f:86
subroutine zswap(N, ZX, INCX, ZY, INCY)
ZSWAP
Definition: zswap.f:81
subroutine zcopy(N, ZX, INCX, ZY, INCY)
ZCOPY
Definition: zcopy.f:81
subroutine zgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
ZGEMV
Definition: zgemv.f:158
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:143
subroutine zlacpy(UPLO, M, N, A, LDA, B, LDB)
ZLACPY copies all or part of one two-dimensional array to another.
Definition: zlacpy.f:103
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,...
Definition: zlanhb.f:132
subroutine zsteqr(COMPZ, N, D, E, Z, LDZ, WORK, INFO)
ZSTEQR
Definition: zsteqr.f:132
subroutine zstein(N, D, E, M, W, IBLOCK, ISPLIT, Z, LDZ, WORK, IWORK, IFAIL, INFO)
ZSTEIN
Definition: zstein.f:182
subroutine zhbtrd(VECT, UPLO, N, KD, AB, LDAB, D, E, Q, LDQ, WORK, INFO)
ZHBTRD
Definition: zhbtrd.f:163
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
Definition: dcopy.f:82
subroutine dscal(N, DA, DX, INCX)
DSCAL
Definition: dscal.f:79
Here is the call graph for this function:
Here is the caller graph for this function: