LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

◆ cheevx()

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
          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 REAL
          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 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.

Definition at line 256 of file cheevx.f.

259 *
260 * -- LAPACK driver routine --
261 * -- LAPACK is a software package provided by Univ. of Tennessee, --
262 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
263 *
264 * .. Scalar Arguments ..
265  CHARACTER JOBZ, RANGE, UPLO
266  INTEGER IL, INFO, IU, LDA, LDZ, LWORK, M, N
267  REAL ABSTOL, VL, VU
268 * ..
269 * .. Array Arguments ..
270  INTEGER IFAIL( * ), IWORK( * )
271  REAL RWORK( * ), W( * )
272  COMPLEX A( LDA, * ), WORK( * ), Z( LDZ, * )
273 * ..
274 *
275 * =====================================================================
276 *
277 * .. Parameters ..
278  REAL ZERO, ONE
279  parameter( zero = 0.0e+0, one = 1.0e+0 )
280  COMPLEX CONE
281  parameter( cone = ( 1.0e+0, 0.0e+0 ) )
282 * ..
283 * .. Local Scalars ..
284  LOGICAL ALLEIG, INDEIG, LOWER, LQUERY, TEST, VALEIG,
285  $ WANTZ
286  CHARACTER ORDER
287  INTEGER I, IINFO, IMAX, INDD, INDE, INDEE, INDIBL,
288  $ INDISP, INDIWK, INDRWK, INDTAU, INDWRK, ISCALE,
289  $ ITMP1, J, JJ, LLWORK, LWKMIN, LWKOPT, NB,
290  $ NSPLIT
291  REAL ABSTLL, ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN,
292  $ SIGMA, SMLNUM, TMP1, VLL, VUU
293 * ..
294 * .. External Functions ..
295  LOGICAL LSAME
296  INTEGER ILAENV
297  REAL SLAMCH, CLANHE
298  EXTERNAL lsame, ilaenv, slamch, clanhe
299 * ..
300 * .. External Subroutines ..
301  EXTERNAL scopy, sscal, sstebz, ssterf, xerbla, csscal,
303  $ cunmtr
304 * ..
305 * .. Intrinsic Functions ..
306  INTRINSIC real, max, min, sqrt
307 * ..
308 * .. Executable Statements ..
309 *
310 * Test the input parameters.
311 *
312  lower = lsame( uplo, 'L' )
313  wantz = lsame( jobz, 'V' )
314  alleig = lsame( range, 'A' )
315  valeig = lsame( range, 'V' )
316  indeig = lsame( range, 'I' )
317  lquery = ( lwork.EQ.-1 )
318 *
319  info = 0
320  IF( .NOT.( wantz .OR. lsame( jobz, 'N' ) ) ) THEN
321  info = -1
322  ELSE IF( .NOT.( alleig .OR. valeig .OR. indeig ) ) THEN
323  info = -2
324  ELSE IF( .NOT.( lower .OR. lsame( uplo, 'U' ) ) ) THEN
325  info = -3
326  ELSE IF( n.LT.0 ) THEN
327  info = -4
328  ELSE IF( lda.LT.max( 1, n ) ) THEN
329  info = -6
330  ELSE
331  IF( valeig ) THEN
332  IF( n.GT.0 .AND. vu.LE.vl )
333  $ info = -8
334  ELSE IF( indeig ) THEN
335  IF( il.LT.1 .OR. il.GT.max( 1, n ) ) THEN
336  info = -9
337  ELSE IF( iu.LT.min( n, il ) .OR. iu.GT.n ) THEN
338  info = -10
339  END IF
340  END IF
341  END IF
342  IF( info.EQ.0 ) THEN
343  IF( ldz.LT.1 .OR. ( wantz .AND. ldz.LT.n ) ) THEN
344  info = -15
345  END IF
346  END IF
347 *
348  IF( info.EQ.0 ) THEN
349  IF( n.LE.1 ) THEN
350  lwkmin = 1
351  work( 1 ) = lwkmin
352  ELSE
353  lwkmin = 2*n
354  nb = ilaenv( 1, 'CHETRD', uplo, n, -1, -1, -1 )
355  nb = max( nb, ilaenv( 1, 'CUNMTR', uplo, n, -1, -1, -1 ) )
356  lwkopt = max( 1, ( nb + 1 )*n )
357  work( 1 ) = lwkopt
358  END IF
359 *
360  IF( lwork.LT.lwkmin .AND. .NOT.lquery )
361  $ info = -17
362  END IF
363 *
364  IF( info.NE.0 ) THEN
365  CALL xerbla( 'CHEEVX', -info )
366  RETURN
367  ELSE IF( lquery ) THEN
368  RETURN
369  END IF
370 *
371 * Quick return if possible
372 *
373  m = 0
374  IF( n.EQ.0 ) THEN
375  RETURN
376  END IF
377 *
378  IF( n.EQ.1 ) THEN
379  IF( alleig .OR. indeig ) THEN
380  m = 1
381  w( 1 ) = real( a( 1, 1 ) )
382  ELSE IF( valeig ) THEN
383  IF( vl.LT.real( a( 1, 1 ) ) .AND. vu.GE.real( a( 1, 1 ) ) )
384  $ THEN
385  m = 1
386  w( 1 ) = real( a( 1, 1 ) )
387  END IF
388  END IF
389  IF( wantz )
390  $ z( 1, 1 ) = cone
391  RETURN
392  END IF
393 *
394 * Get machine constants.
395 *
396  safmin = slamch( 'Safe minimum' )
397  eps = slamch( 'Precision' )
398  smlnum = safmin / eps
399  bignum = one / smlnum
400  rmin = sqrt( smlnum )
401  rmax = min( sqrt( bignum ), one / sqrt( sqrt( safmin ) ) )
402 *
403 * Scale matrix to allowable range, if necessary.
404 *
405  iscale = 0
406  abstll = abstol
407  IF( valeig ) THEN
408  vll = vl
409  vuu = vu
410  END IF
411  anrm = clanhe( 'M', uplo, n, a, lda, rwork )
412  IF( anrm.GT.zero .AND. anrm.LT.rmin ) THEN
413  iscale = 1
414  sigma = rmin / anrm
415  ELSE IF( anrm.GT.rmax ) THEN
416  iscale = 1
417  sigma = rmax / anrm
418  END IF
419  IF( iscale.EQ.1 ) THEN
420  IF( lower ) THEN
421  DO 10 j = 1, n
422  CALL csscal( n-j+1, sigma, a( j, j ), 1 )
423  10 CONTINUE
424  ELSE
425  DO 20 j = 1, n
426  CALL csscal( j, sigma, a( 1, j ), 1 )
427  20 CONTINUE
428  END IF
429  IF( abstol.GT.0 )
430  $ abstll = abstol*sigma
431  IF( valeig ) THEN
432  vll = vl*sigma
433  vuu = vu*sigma
434  END IF
435  END IF
436 *
437 * Call CHETRD to reduce Hermitian matrix to tridiagonal form.
438 *
439  indd = 1
440  inde = indd + n
441  indrwk = inde + n
442  indtau = 1
443  indwrk = indtau + n
444  llwork = lwork - indwrk + 1
445  CALL chetrd( uplo, n, a, lda, rwork( indd ), rwork( inde ),
446  $ work( indtau ), work( indwrk ), llwork, iinfo )
447 *
448 * If all eigenvalues are desired and ABSTOL is less than or equal to
449 * zero, then call SSTERF or CUNGTR and CSTEQR. If this fails for
450 * some eigenvalue, then try SSTEBZ.
451 *
452  test = .false.
453  IF( indeig ) THEN
454  IF( il.EQ.1 .AND. iu.EQ.n ) THEN
455  test = .true.
456  END IF
457  END IF
458  IF( ( alleig .OR. test ) .AND. ( abstol.LE.zero ) ) THEN
459  CALL scopy( n, rwork( indd ), 1, w, 1 )
460  indee = indrwk + 2*n
461  IF( .NOT.wantz ) THEN
462  CALL scopy( n-1, rwork( inde ), 1, rwork( indee ), 1 )
463  CALL ssterf( n, w, rwork( indee ), info )
464  ELSE
465  CALL clacpy( 'A', n, n, a, lda, z, ldz )
466  CALL cungtr( uplo, n, z, ldz, work( indtau ),
467  $ work( indwrk ), llwork, iinfo )
468  CALL scopy( n-1, rwork( inde ), 1, rwork( indee ), 1 )
469  CALL csteqr( jobz, n, w, rwork( indee ), z, ldz,
470  $ rwork( indrwk ), info )
471  IF( info.EQ.0 ) THEN
472  DO 30 i = 1, n
473  ifail( i ) = 0
474  30 CONTINUE
475  END IF
476  END IF
477  IF( info.EQ.0 ) THEN
478  m = n
479  GO TO 40
480  END IF
481  info = 0
482  END IF
483 *
484 * Otherwise, call SSTEBZ and, if eigenvectors are desired, CSTEIN.
485 *
486  IF( wantz ) THEN
487  order = 'B'
488  ELSE
489  order = 'E'
490  END IF
491  indibl = 1
492  indisp = indibl + n
493  indiwk = indisp + n
494  CALL sstebz( range, order, n, vll, vuu, il, iu, abstll,
495  $ rwork( indd ), rwork( inde ), m, nsplit, w,
496  $ iwork( indibl ), iwork( indisp ), rwork( indrwk ),
497  $ iwork( indiwk ), info )
498 *
499  IF( wantz ) THEN
500  CALL cstein( n, rwork( indd ), rwork( inde ), m, w,
501  $ iwork( indibl ), iwork( indisp ), z, ldz,
502  $ rwork( indrwk ), iwork( indiwk ), ifail, info )
503 *
504 * Apply unitary matrix used in reduction to tridiagonal
505 * form to eigenvectors returned by CSTEIN.
506 *
507  CALL cunmtr( 'L', uplo, 'N', n, m, a, lda, work( indtau ), z,
508  $ ldz, work( indwrk ), llwork, iinfo )
509  END IF
510 *
511 * If matrix was scaled, then rescale eigenvalues appropriately.
512 *
513  40 CONTINUE
514  IF( iscale.EQ.1 ) THEN
515  IF( info.EQ.0 ) THEN
516  imax = m
517  ELSE
518  imax = info - 1
519  END IF
520  CALL sscal( imax, one / sigma, w, 1 )
521  END IF
522 *
523 * If eigenvalues are not in order, then sort them, along with
524 * eigenvectors.
525 *
526  IF( wantz ) THEN
527  DO 60 j = 1, m - 1
528  i = 0
529  tmp1 = w( j )
530  DO 50 jj = j + 1, m
531  IF( w( jj ).LT.tmp1 ) THEN
532  i = jj
533  tmp1 = w( jj )
534  END IF
535  50 CONTINUE
536 *
537  IF( i.NE.0 ) THEN
538  itmp1 = iwork( indibl+i-1 )
539  w( i ) = w( j )
540  iwork( indibl+i-1 ) = iwork( indibl+j-1 )
541  w( j ) = tmp1
542  iwork( indibl+j-1 ) = itmp1
543  CALL cswap( n, z( 1, i ), 1, z( 1, j ), 1 )
544  IF( info.NE.0 ) THEN
545  itmp1 = ifail( i )
546  ifail( i ) = ifail( j )
547  ifail( j ) = itmp1
548  END IF
549  END IF
550  60 CONTINUE
551  END IF
552 *
553 * Set WORK(1) to optimal complex workspace size.
554 *
555  work( 1 ) = lwkopt
556 *
557  RETURN
558 *
559 * End of CHEEVX
560 *
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
ILAENV
Definition: ilaenv.f:162
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
subroutine ssterf(N, D, E, INFO)
SSTERF
Definition: ssterf.f:86
subroutine sstebz(RANGE, ORDER, N, VL, VU, IL, IU, ABSTOL, D, E, M, NSPLIT, W, IBLOCK, ISPLIT, WORK, IWORK, INFO)
SSTEBZ
Definition: sstebz.f:273
subroutine csscal(N, SA, CX, INCX)
CSSCAL
Definition: csscal.f:78
subroutine cswap(N, CX, INCX, CY, INCY)
CSWAP
Definition: cswap.f:81
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,...
Definition: clanhe.f:124
subroutine chetrd(UPLO, N, A, LDA, D, E, TAU, WORK, LWORK, INFO)
CHETRD
Definition: chetrd.f:192
subroutine clacpy(UPLO, M, N, A, LDA, B, LDB)
CLACPY copies all or part of one two-dimensional array to another.
Definition: clacpy.f:103
subroutine cunmtr(SIDE, UPLO, TRANS, M, N, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
CUNMTR
Definition: cunmtr.f:172
subroutine cstein(N, D, E, M, W, IBLOCK, ISPLIT, Z, LDZ, WORK, IWORK, IFAIL, INFO)
CSTEIN
Definition: cstein.f:182
subroutine csteqr(COMPZ, N, D, E, Z, LDZ, WORK, INFO)
CSTEQR
Definition: csteqr.f:132
subroutine cungtr(UPLO, N, A, LDA, TAU, WORK, LWORK, INFO)
CUNGTR
Definition: cungtr.f:123
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY
Definition: scopy.f:82
subroutine sscal(N, SA, SX, INCX)
SSCAL
Definition: sscal.f:79
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:68
Here is the call graph for this function:
Here is the caller graph for this function: