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

Go to the source code of this file.

Functions/Subroutines

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

Function/Subroutine Documentation

subroutine cheevx ( character  JOBZ,
character  RANGE,
character  UPLO,
integer  N,
complex, dimension( lda, * )  A,
integer  LDA,
real  VL,
real  VU,
integer  IL,
integer  IU,
real  ABSTOL,
integer  M,
real, dimension( * )  W,
complex, dimension( ldz, * )  Z,
integer  LDZ,
complex, dimension( * )  WORK,
integer  LWORK,
real, dimension( * )  RWORK,
integer, dimension( * )  IWORK,
integer, dimension( * )  IFAIL,
integer  INFO 
)

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

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

Purpose:
 CHEEVX 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 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 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 A 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)
          On normal exit, 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 (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 CHETRD and for
          CUNMTR 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 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 254 of file cheevx.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  REAL abstol, vl, vu
264 * ..
265 * .. Array Arguments ..
266  INTEGER ifail( * ), iwork( * )
267  REAL rwork( * ), w( * )
268  COMPLEX a( lda, * ), work( * ), z( ldz, * )
269 * ..
270 *
271 * =====================================================================
272 *
273 * .. Parameters ..
274  REAL zero, one
275  parameter( zero = 0.0e+0, one = 1.0e+0 )
276  COMPLEX cone
277  parameter( cone = ( 1.0e+0, 0.0e+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  REAL abstll, anrm, bignum, eps, rmax, rmin, safmin,
288  $ sigma, smlnum, tmp1, vll, vuu
289 * ..
290 * .. External Functions ..
291  LOGICAL lsame
292  INTEGER ilaenv
293  REAL slamch, clanhe
294  EXTERNAL lsame, ilaenv, slamch, clanhe
295 * ..
296 * .. External Subroutines ..
297  EXTERNAL scopy, sscal, sstebz, ssterf, xerbla, csscal,
299  $ cunmtr
300 * ..
301 * .. Intrinsic Functions ..
302  INTRINSIC REAL, 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, 'CHETRD', uplo, n, -1, -1, -1 )
351  nb = max( nb, ilaenv( 1, 'CUNMTR', 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( 'CHEEVX', -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.REAL( A( 1, 1 ) ) .AND. vu.GE.REAL( 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 = slamch( 'Safe minimum' )
393  eps = slamch( '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 = clanhe( '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 csscal( n-j+1, sigma, a( j, j ), 1 )
419  10 CONTINUE
420  ELSE
421  DO 20 j = 1, n
422  CALL csscal( 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 CHETRD 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 chetrd( 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 SSTERF or CUNGTR and CSTEQR. If this fails for
446 * some eigenvalue, then try SSTEBZ.
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 scopy( n, rwork( indd ), 1, w, 1 )
456  indee = indrwk + 2*n
457  IF( .NOT.wantz ) THEN
458  CALL scopy( n-1, rwork( inde ), 1, rwork( indee ), 1 )
459  CALL ssterf( n, w, rwork( indee ), info )
460  ELSE
461  CALL clacpy( 'A', n, n, a, lda, z, ldz )
462  CALL cungtr( uplo, n, z, ldz, work( indtau ),
463  $ work( indwrk ), llwork, iinfo )
464  CALL scopy( n-1, rwork( inde ), 1, rwork( indee ), 1 )
465  CALL csteqr( 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 SSTEBZ and, if eigenvectors are desired, CSTEIN.
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 sstebz( 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 cstein( 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 CSTEIN.
502 *
503  CALL cunmtr( '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 sscal( 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 cswap( 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 CHEEVX
556 *
real function clanhe(NORM, UPLO, N, A, LDA, WORK)
CLANHE 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: clanhe.f:126
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine cunmtr(SIDE, UPLO, TRANS, M, N, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
CUNMTR
Definition: cunmtr.f:174
subroutine cungtr(UPLO, N, A, LDA, TAU, WORK, LWORK, INFO)
CUNGTR
Definition: cungtr.f:125
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 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
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
Definition: tstiee.f:83
subroutine csscal(N, SA, CX, INCX)
CSSCAL
Definition: csscal.f:54
subroutine ssterf(N, D, E, INFO)
SSTERF
Definition: ssterf.f:88
subroutine chetrd(UPLO, N, A, LDA, D, E, TAU, WORK, LWORK, INFO)
CHETRD
Definition: chetrd.f:194
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: