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

Go to the source code of this file.

Functions/Subroutines

subroutine zheevx (JOBZ, RANGE, UPLO, N, A, LDA, VL, VU, IL, IU, ABSTOL, M, W, Z, LDZ, WORK, LWORK, RWORK, IWORK, IFAIL, INFO)
  ZHEEVX computes the eigenvalues and, optionally, the left and/or right eigenvectors for HE matrices More...
 

Function/Subroutine Documentation

subroutine zheevx ( character  JOBZ,
character  RANGE,
character  UPLO,
integer  N,
complex*16, dimension( lda, * )  A,
integer  LDA,
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,
integer  LWORK,
double precision, dimension( * )  RWORK,
integer, dimension( * )  IWORK,
integer, dimension( * )  IFAIL,
integer  INFO 
)

ZHEEVX computes the eigenvalues and, optionally, the left and/or right eigenvectors for HE matrices

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

Purpose:
 ZHEEVX computes selected eigenvalues and, optionally, eigenvectors
 of a complex Hermitian 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,out]A
          A is COMPLEX*16 array, dimension (LDA, N)
          On entry, the Hermitian matrix A.  If UPLO = 'U', the
          leading N-by-N upper triangular part of A contains the
          upper triangular part of the matrix A.  If UPLO = 'L',
          the leading N-by-N lower triangular part of A contains
          the lower triangular part of the matrix A.
          On exit, the lower triangle (if UPLO='L') or the upper
          triangle (if UPLO='U') of A, including the diagonal, is
          destroyed.
[in]LDA
          LDA is INTEGER
          The leading dimension of the array A.  LDA >= 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 A 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)
          On normal exit, 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 (MAX(1,LWORK))
          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
[in]LWORK
          LWORK is INTEGER
          The length of the array WORK.  LWORK >= 1, when N <= 1;
          otherwise 2*N.
          For optimal efficiency, LWORK >= (NB+1)*N,
          where NB is the max of the blocksize for ZHETRD and for
          ZUNMTR as returned by ILAENV.

          If LWORK = -1, then a workspace query is assumed; the routine
          only calculates the optimal size of the WORK array, returns
          this value as the first entry of the WORK array, and no error
          message related to LWORK is issued by XERBLA.
[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 254 of file zheevx.f.

254 *
255 * -- LAPACK driver routine (version 3.4.0) --
256 * -- LAPACK is a software package provided by Univ. of Tennessee, --
257 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
258 * November 2011
259 *
260 * .. Scalar Arguments ..
261  CHARACTER jobz, range, uplo
262  INTEGER il, info, iu, lda, ldz, lwork, m, n
263  DOUBLE PRECISION abstol, vl, vu
264 * ..
265 * .. Array Arguments ..
266  INTEGER ifail( * ), iwork( * )
267  DOUBLE PRECISION rwork( * ), w( * )
268  COMPLEX*16 a( lda, * ), work( * ), z( ldz, * )
269 * ..
270 *
271 * =====================================================================
272 *
273 * .. Parameters ..
274  DOUBLE PRECISION zero, one
275  parameter( zero = 0.0d+0, one = 1.0d+0 )
276  COMPLEX*16 cone
277  parameter( cone = ( 1.0d+0, 0.0d+0 ) )
278 * ..
279 * .. Local Scalars ..
280  LOGICAL alleig, indeig, lower, lquery, test, valeig,
281  $ wantz
282  CHARACTER order
283  INTEGER i, iinfo, imax, indd, inde, indee, indibl,
284  $ indisp, indiwk, indrwk, indtau, indwrk, iscale,
285  $ itmp1, j, jj, llwork, lwkmin, lwkopt, nb,
286  $ nsplit
287  DOUBLE PRECISION abstll, anrm, bignum, eps, rmax, rmin, safmin,
288  $ sigma, smlnum, tmp1, vll, vuu
289 * ..
290 * .. External Functions ..
291  LOGICAL lsame
292  INTEGER ilaenv
293  DOUBLE PRECISION dlamch, zlanhe
294  EXTERNAL lsame, ilaenv, dlamch, zlanhe
295 * ..
296 * .. External Subroutines ..
297  EXTERNAL dcopy, dscal, dstebz, dsterf, xerbla, zdscal,
299  $ zunmtr
300 * ..
301 * .. Intrinsic Functions ..
302  INTRINSIC dble, max, min, sqrt
303 * ..
304 * .. Executable Statements ..
305 *
306 * Test the input parameters.
307 *
308  lower = lsame( uplo, 'L' )
309  wantz = lsame( jobz, 'V' )
310  alleig = lsame( range, 'A' )
311  valeig = lsame( range, 'V' )
312  indeig = lsame( range, 'I' )
313  lquery = ( lwork.EQ.-1 )
314 *
315  info = 0
316  IF( .NOT.( wantz .OR. lsame( jobz, 'N' ) ) ) THEN
317  info = -1
318  ELSE IF( .NOT.( alleig .OR. valeig .OR. indeig ) ) THEN
319  info = -2
320  ELSE IF( .NOT.( lower .OR. lsame( uplo, 'U' ) ) ) THEN
321  info = -3
322  ELSE IF( n.LT.0 ) THEN
323  info = -4
324  ELSE IF( lda.LT.max( 1, n ) ) THEN
325  info = -6
326  ELSE
327  IF( valeig ) THEN
328  IF( n.GT.0 .AND. vu.LE.vl )
329  $ info = -8
330  ELSE IF( indeig ) THEN
331  IF( il.LT.1 .OR. il.GT.max( 1, n ) ) THEN
332  info = -9
333  ELSE IF( iu.LT.min( n, il ) .OR. iu.GT.n ) THEN
334  info = -10
335  END IF
336  END IF
337  END IF
338  IF( info.EQ.0 ) THEN
339  IF( ldz.LT.1 .OR. ( wantz .AND. ldz.LT.n ) ) THEN
340  info = -15
341  END IF
342  END IF
343 *
344  IF( info.EQ.0 ) THEN
345  IF( n.LE.1 ) THEN
346  lwkmin = 1
347  work( 1 ) = lwkmin
348  ELSE
349  lwkmin = 2*n
350  nb = ilaenv( 1, 'ZHETRD', uplo, n, -1, -1, -1 )
351  nb = max( nb, ilaenv( 1, 'ZUNMTR', uplo, n, -1, -1, -1 ) )
352  lwkopt = max( 1, ( nb + 1 )*n )
353  work( 1 ) = lwkopt
354  END IF
355 *
356  IF( lwork.LT.lwkmin .AND. .NOT.lquery )
357  $ info = -17
358  END IF
359 *
360  IF( info.NE.0 ) THEN
361  CALL xerbla( 'ZHEEVX', -info )
362  RETURN
363  ELSE IF( lquery ) THEN
364  RETURN
365  END IF
366 *
367 * Quick return if possible
368 *
369  m = 0
370  IF( n.EQ.0 ) THEN
371  RETURN
372  END IF
373 *
374  IF( n.EQ.1 ) THEN
375  IF( alleig .OR. indeig ) THEN
376  m = 1
377  w( 1 ) = a( 1, 1 )
378  ELSE IF( valeig ) THEN
379  IF( vl.LT.dble( a( 1, 1 ) ) .AND. vu.GE.dble( a( 1, 1 ) ) )
380  $ THEN
381  m = 1
382  w( 1 ) = a( 1, 1 )
383  END IF
384  END IF
385  IF( wantz )
386  $ z( 1, 1 ) = cone
387  RETURN
388  END IF
389 *
390 * Get machine constants.
391 *
392  safmin = dlamch( 'Safe minimum' )
393  eps = dlamch( 'Precision' )
394  smlnum = safmin / eps
395  bignum = one / smlnum
396  rmin = sqrt( smlnum )
397  rmax = min( sqrt( bignum ), one / sqrt( sqrt( safmin ) ) )
398 *
399 * Scale matrix to allowable range, if necessary.
400 *
401  iscale = 0
402  abstll = abstol
403  IF( valeig ) THEN
404  vll = vl
405  vuu = vu
406  END IF
407  anrm = zlanhe( 'M', uplo, n, a, lda, rwork )
408  IF( anrm.GT.zero .AND. anrm.LT.rmin ) THEN
409  iscale = 1
410  sigma = rmin / anrm
411  ELSE IF( anrm.GT.rmax ) THEN
412  iscale = 1
413  sigma = rmax / anrm
414  END IF
415  IF( iscale.EQ.1 ) THEN
416  IF( lower ) THEN
417  DO 10 j = 1, n
418  CALL zdscal( n-j+1, sigma, a( j, j ), 1 )
419  10 CONTINUE
420  ELSE
421  DO 20 j = 1, n
422  CALL zdscal( j, sigma, a( 1, j ), 1 )
423  20 CONTINUE
424  END IF
425  IF( abstol.GT.0 )
426  $ abstll = abstol*sigma
427  IF( valeig ) THEN
428  vll = vl*sigma
429  vuu = vu*sigma
430  END IF
431  END IF
432 *
433 * Call ZHETRD to reduce Hermitian matrix to tridiagonal form.
434 *
435  indd = 1
436  inde = indd + n
437  indrwk = inde + n
438  indtau = 1
439  indwrk = indtau + n
440  llwork = lwork - indwrk + 1
441  CALL zhetrd( uplo, n, a, lda, rwork( indd ), rwork( inde ),
442  $ work( indtau ), work( indwrk ), llwork, iinfo )
443 *
444 * If all eigenvalues are desired and ABSTOL is less than or equal to
445 * zero, then call DSTERF or ZUNGTR and ZSTEQR. If this fails for
446 * some eigenvalue, then try DSTEBZ.
447 *
448  test = .false.
449  IF( indeig ) THEN
450  IF( il.EQ.1 .AND. iu.EQ.n ) THEN
451  test = .true.
452  END IF
453  END IF
454  IF( ( alleig .OR. test ) .AND. ( abstol.LE.zero ) ) THEN
455  CALL dcopy( n, rwork( indd ), 1, w, 1 )
456  indee = indrwk + 2*n
457  IF( .NOT.wantz ) THEN
458  CALL dcopy( n-1, rwork( inde ), 1, rwork( indee ), 1 )
459  CALL dsterf( n, w, rwork( indee ), info )
460  ELSE
461  CALL zlacpy( 'A', n, n, a, lda, z, ldz )
462  CALL zungtr( uplo, n, z, ldz, work( indtau ),
463  $ work( indwrk ), llwork, iinfo )
464  CALL dcopy( n-1, rwork( inde ), 1, rwork( indee ), 1 )
465  CALL zsteqr( jobz, n, w, rwork( indee ), z, ldz,
466  $ rwork( indrwk ), info )
467  IF( info.EQ.0 ) THEN
468  DO 30 i = 1, n
469  ifail( i ) = 0
470  30 CONTINUE
471  END IF
472  END IF
473  IF( info.EQ.0 ) THEN
474  m = n
475  GO TO 40
476  END IF
477  info = 0
478  END IF
479 *
480 * Otherwise, call DSTEBZ and, if eigenvectors are desired, ZSTEIN.
481 *
482  IF( wantz ) THEN
483  order = 'B'
484  ELSE
485  order = 'E'
486  END IF
487  indibl = 1
488  indisp = indibl + n
489  indiwk = indisp + n
490  CALL dstebz( range, order, n, vll, vuu, il, iu, abstll,
491  $ rwork( indd ), rwork( inde ), m, nsplit, w,
492  $ iwork( indibl ), iwork( indisp ), rwork( indrwk ),
493  $ iwork( indiwk ), info )
494 *
495  IF( wantz ) THEN
496  CALL zstein( n, rwork( indd ), rwork( inde ), m, w,
497  $ iwork( indibl ), iwork( indisp ), z, ldz,
498  $ rwork( indrwk ), iwork( indiwk ), ifail, info )
499 *
500 * Apply unitary matrix used in reduction to tridiagonal
501 * form to eigenvectors returned by ZSTEIN.
502 *
503  CALL zunmtr( 'L', uplo, 'N', n, m, a, lda, work( indtau ), z,
504  $ ldz, work( indwrk ), llwork, iinfo )
505  END IF
506 *
507 * If matrix was scaled, then rescale eigenvalues appropriately.
508 *
509  40 CONTINUE
510  IF( iscale.EQ.1 ) THEN
511  IF( info.EQ.0 ) THEN
512  imax = m
513  ELSE
514  imax = info - 1
515  END IF
516  CALL dscal( imax, one / sigma, w, 1 )
517  END IF
518 *
519 * If eigenvalues are not in order, then sort them, along with
520 * eigenvectors.
521 *
522  IF( wantz ) THEN
523  DO 60 j = 1, m - 1
524  i = 0
525  tmp1 = w( j )
526  DO 50 jj = j + 1, m
527  IF( w( jj ).LT.tmp1 ) THEN
528  i = jj
529  tmp1 = w( jj )
530  END IF
531  50 CONTINUE
532 *
533  IF( i.NE.0 ) THEN
534  itmp1 = iwork( indibl+i-1 )
535  w( i ) = w( j )
536  iwork( indibl+i-1 ) = iwork( indibl+j-1 )
537  w( j ) = tmp1
538  iwork( indibl+j-1 ) = itmp1
539  CALL zswap( n, z( 1, i ), 1, z( 1, j ), 1 )
540  IF( info.NE.0 ) THEN
541  itmp1 = ifail( i )
542  ifail( i ) = ifail( j )
543  ifail( j ) = itmp1
544  END IF
545  END IF
546  60 CONTINUE
547  END IF
548 *
549 * Set WORK(1) to optimal complex workspace size.
550 *
551  work( 1 ) = lwkopt
552 *
553  RETURN
554 *
555 * End of ZHEEVX
556 *
subroutine zswap(N, ZX, INCX, ZY, INCY)
ZSWAP
Definition: zswap.f:52
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 zdscal(N, DA, ZX, INCX)
ZDSCAL
Definition: zdscal.f:54
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
double precision function zlanhe(NORM, UPLO, N, A, LDA, WORK)
ZLANHE returns the value of the 1-norm, or the Frobenius norm, or the infinity norm, or the element of largest absolute value of a complex Hermitian matrix.
Definition: zlanhe.f:126
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine zhetrd(UPLO, N, A, LDA, D, E, TAU, WORK, LWORK, INFO)
ZHETRD
Definition: zhetrd.f:194
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
subroutine dscal(N, DA, DX, INCX)
DSCAL
Definition: dscal.f:55
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
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
Definition: tstiee.f:83
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
subroutine zunmtr(SIDE, UPLO, TRANS, M, N, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
ZUNMTR
Definition: zunmtr.f:173
subroutine zungtr(UPLO, N, A, LDA, TAU, WORK, LWORK, INFO)
ZUNGTR
Definition: zungtr.f:125

Here is the call graph for this function:

Here is the caller graph for this function: