LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
subroutine chpevx ( character  JOBZ,
character  RANGE,
character  UPLO,
integer  N,
complex, dimension( * )  AP,
real  VL,
real  VU,
integer  IL,
integer  IU,
real  ABSTOL,
integer  M,
real, dimension( * )  W,
complex, dimension( ldz, * )  Z,
integer  LDZ,
complex, dimension( * )  WORK,
real, dimension( * )  RWORK,
integer, dimension( * )  IWORK,
integer, dimension( * )  IFAIL,
integer  INFO 
)

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

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

Purpose:
 CHPEVX computes selected eigenvalues and, optionally, eigenvectors
 of a complex Hermitian matrix A in packed storage.
 Eigenvalues/vectors 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]AP
          AP is COMPLEX array, dimension (N*(N+1)/2)
          On entry, the upper or lower triangle of the Hermitian matrix
          A, packed columnwise in a linear array.  The j-th column of A
          is stored in the array AP as follows:
          if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
          if UPLO = 'L', AP(i + (j-1)*(2*n-j)/2) = A(i,j) for j<=i<=n.

          On exit, AP is overwritten by values generated during the
          reduction to tridiagonal form.  If UPLO = 'U', the diagonal
          and first superdiagonal of the tridiagonal matrix T overwrite
          the corresponding elements of A, and if UPLO = 'L', the
          diagonal and first subdiagonal of T overwrite the
          corresponding elements of A.
[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 AP 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)
          If INFO = 0, 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 (2*N)
[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
June 2016

Definition at line 242 of file chpevx.f.

242 *
243 * -- LAPACK driver routine (version 3.6.1) --
244 * -- LAPACK is a software package provided by Univ. of Tennessee, --
245 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
246 * June 2016
247 *
248 * .. Scalar Arguments ..
249  CHARACTER jobz, range, uplo
250  INTEGER il, info, iu, ldz, m, n
251  REAL abstol, vl, vu
252 * ..
253 * .. Array Arguments ..
254  INTEGER ifail( * ), iwork( * )
255  REAL rwork( * ), w( * )
256  COMPLEX ap( * ), work( * ), z( ldz, * )
257 * ..
258 *
259 * =====================================================================
260 *
261 * .. Parameters ..
262  REAL zero, one
263  parameter ( zero = 0.0e0, one = 1.0e0 )
264  COMPLEX cone
265  parameter ( cone = ( 1.0e0, 0.0e0 ) )
266 * ..
267 * .. Local Scalars ..
268  LOGICAL alleig, indeig, test, valeig, wantz
269  CHARACTER order
270  INTEGER i, iinfo, imax, indd, inde, indee, indibl,
271  $ indisp, indiwk, indrwk, indtau, indwrk, iscale,
272  $ itmp1, j, jj, nsplit
273  REAL abstll, anrm, bignum, eps, rmax, rmin, safmin,
274  $ sigma, smlnum, tmp1, vll, vuu
275 * ..
276 * .. External Functions ..
277  LOGICAL lsame
278  REAL clanhp, slamch
279  EXTERNAL lsame, clanhp, slamch
280 * ..
281 * .. External Subroutines ..
282  EXTERNAL chptrd, csscal, cstein, csteqr, cswap, cupgtr,
284 * ..
285 * .. Intrinsic Functions ..
286  INTRINSIC max, min, REAL, sqrt
287 * ..
288 * .. Executable Statements ..
289 *
290 * Test the input parameters.
291 *
292  wantz = lsame( jobz, 'V' )
293  alleig = lsame( range, 'A' )
294  valeig = lsame( range, 'V' )
295  indeig = lsame( range, 'I' )
296 *
297  info = 0
298  IF( .NOT.( wantz .OR. lsame( jobz, 'N' ) ) ) THEN
299  info = -1
300  ELSE IF( .NOT.( alleig .OR. valeig .OR. indeig ) ) THEN
301  info = -2
302  ELSE IF( .NOT.( lsame( uplo, 'L' ) .OR. lsame( uplo, 'U' ) ) )
303  $ THEN
304  info = -3
305  ELSE IF( n.LT.0 ) THEN
306  info = -4
307  ELSE
308  IF( valeig ) THEN
309  IF( n.GT.0 .AND. vu.LE.vl )
310  $ info = -7
311  ELSE IF( indeig ) THEN
312  IF( il.LT.1 .OR. il.GT.max( 1, n ) ) THEN
313  info = -8
314  ELSE IF( iu.LT.min( n, il ) .OR. iu.GT.n ) THEN
315  info = -9
316  END IF
317  END IF
318  END IF
319  IF( info.EQ.0 ) THEN
320  IF( ldz.LT.1 .OR. ( wantz .AND. ldz.LT.n ) )
321  $ info = -14
322  END IF
323 *
324  IF( info.NE.0 ) THEN
325  CALL xerbla( 'CHPEVX', -info )
326  RETURN
327  END IF
328 *
329 * Quick return if possible
330 *
331  m = 0
332  IF( n.EQ.0 )
333  $ RETURN
334 *
335  IF( n.EQ.1 ) THEN
336  IF( alleig .OR. indeig ) THEN
337  m = 1
338  w( 1 ) = ap( 1 )
339  ELSE
340  IF( vl.LT.REAL( AP( 1 ) ) .AND. vu.GE.REAL( AP( 1 ) ) ) then
341  m = 1
342  w( 1 ) = ap( 1 )
343  END IF
344  END IF
345  IF( wantz )
346  $ z( 1, 1 ) = cone
347  RETURN
348  END IF
349 *
350 * Get machine constants.
351 *
352  safmin = slamch( 'Safe minimum' )
353  eps = slamch( 'Precision' )
354  smlnum = safmin / eps
355  bignum = one / smlnum
356  rmin = sqrt( smlnum )
357  rmax = min( sqrt( bignum ), one / sqrt( sqrt( safmin ) ) )
358 *
359 * Scale matrix to allowable range, if necessary.
360 *
361  iscale = 0
362  abstll = abstol
363  IF ( valeig ) THEN
364  vll = vl
365  vuu = vu
366  ELSE
367  vll = zero
368  vuu = zero
369  ENDIF
370  anrm = clanhp( 'M', uplo, n, ap, rwork )
371  IF( anrm.GT.zero .AND. anrm.LT.rmin ) THEN
372  iscale = 1
373  sigma = rmin / anrm
374  ELSE IF( anrm.GT.rmax ) THEN
375  iscale = 1
376  sigma = rmax / anrm
377  END IF
378  IF( iscale.EQ.1 ) THEN
379  CALL csscal( ( n*( n+1 ) ) / 2, sigma, ap, 1 )
380  IF( abstol.GT.0 )
381  $ abstll = abstol*sigma
382  IF( valeig ) THEN
383  vll = vl*sigma
384  vuu = vu*sigma
385  END IF
386  END IF
387 *
388 * Call CHPTRD to reduce Hermitian packed matrix to tridiagonal form.
389 *
390  indd = 1
391  inde = indd + n
392  indrwk = inde + n
393  indtau = 1
394  indwrk = indtau + n
395  CALL chptrd( uplo, n, ap, rwork( indd ), rwork( inde ),
396  $ work( indtau ), iinfo )
397 *
398 * If all eigenvalues are desired and ABSTOL is less than or equal
399 * to zero, then call SSTERF or CUPGTR and CSTEQR. If this fails
400 * for some eigenvalue, then try SSTEBZ.
401 *
402  test = .false.
403  IF (indeig) THEN
404  IF (il.EQ.1 .AND. iu.EQ.n) THEN
405  test = .true.
406  END IF
407  END IF
408  IF ((alleig .OR. test) .AND. (abstol.LE.zero)) THEN
409  CALL scopy( n, rwork( indd ), 1, w, 1 )
410  indee = indrwk + 2*n
411  IF( .NOT.wantz ) THEN
412  CALL scopy( n-1, rwork( inde ), 1, rwork( indee ), 1 )
413  CALL ssterf( n, w, rwork( indee ), info )
414  ELSE
415  CALL cupgtr( uplo, n, ap, work( indtau ), z, ldz,
416  $ work( indwrk ), iinfo )
417  CALL scopy( n-1, rwork( inde ), 1, rwork( indee ), 1 )
418  CALL csteqr( jobz, n, w, rwork( indee ), z, ldz,
419  $ rwork( indrwk ), info )
420  IF( info.EQ.0 ) THEN
421  DO 10 i = 1, n
422  ifail( i ) = 0
423  10 CONTINUE
424  END IF
425  END IF
426  IF( info.EQ.0 ) THEN
427  m = n
428  GO TO 20
429  END IF
430  info = 0
431  END IF
432 *
433 * Otherwise, call SSTEBZ and, if eigenvectors are desired, CSTEIN.
434 *
435  IF( wantz ) THEN
436  order = 'B'
437  ELSE
438  order = 'E'
439  END IF
440  indibl = 1
441  indisp = indibl + n
442  indiwk = indisp + n
443  CALL sstebz( range, order, n, vll, vuu, il, iu, abstll,
444  $ rwork( indd ), rwork( inde ), m, nsplit, w,
445  $ iwork( indibl ), iwork( indisp ), rwork( indrwk ),
446  $ iwork( indiwk ), info )
447 *
448  IF( wantz ) THEN
449  CALL cstein( n, rwork( indd ), rwork( inde ), m, w,
450  $ iwork( indibl ), iwork( indisp ), z, ldz,
451  $ rwork( indrwk ), iwork( indiwk ), ifail, info )
452 *
453 * Apply unitary matrix used in reduction to tridiagonal
454 * form to eigenvectors returned by CSTEIN.
455 *
456  indwrk = indtau + n
457  CALL cupmtr( 'L', uplo, 'N', n, m, ap, work( indtau ), z, ldz,
458  $ work( indwrk ), iinfo )
459  END IF
460 *
461 * If matrix was scaled, then rescale eigenvalues appropriately.
462 *
463  20 CONTINUE
464  IF( iscale.EQ.1 ) THEN
465  IF( info.EQ.0 ) THEN
466  imax = m
467  ELSE
468  imax = info - 1
469  END IF
470  CALL sscal( imax, one / sigma, w, 1 )
471  END IF
472 *
473 * If eigenvalues are not in order, then sort them, along with
474 * eigenvectors.
475 *
476  IF( wantz ) THEN
477  DO 40 j = 1, m - 1
478  i = 0
479  tmp1 = w( j )
480  DO 30 jj = j + 1, m
481  IF( w( jj ).LT.tmp1 ) THEN
482  i = jj
483  tmp1 = w( jj )
484  END IF
485  30 CONTINUE
486 *
487  IF( i.NE.0 ) THEN
488  itmp1 = iwork( indibl+i-1 )
489  w( i ) = w( j )
490  iwork( indibl+i-1 ) = iwork( indibl+j-1 )
491  w( j ) = tmp1
492  iwork( indibl+j-1 ) = itmp1
493  CALL cswap( n, z( 1, i ), 1, z( 1, j ), 1 )
494  IF( info.NE.0 ) THEN
495  itmp1 = ifail( i )
496  ifail( i ) = ifail( j )
497  ifail( j ) = itmp1
498  END IF
499  END IF
500  40 CONTINUE
501  END IF
502 *
503  RETURN
504 *
505 * End of CHPEVX
506 *
subroutine sstebz(RANGE, ORDER, N, VL, VU, IL, IU, ABSTOL, D, E, M, NSPLIT, W, IBLOCK, ISPLIT, WORK, IWORK, INFO)
SSTEBZ
Definition: sstebz.f:275
subroutine cupgtr(UPLO, N, AP, TAU, Q, LDQ, WORK, INFO)
CUPGTR
Definition: cupgtr.f:116
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
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 chptrd(UPLO, N, AP, D, E, TAU, INFO)
CHPTRD
Definition: chptrd.f:153
subroutine sscal(N, SA, SX, INCX)
SSCAL
Definition: sscal.f:55
real function clanhp(NORM, UPLO, N, AP, WORK)
CLANHP 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 supplied in packed form.
Definition: clanhp.f:119
subroutine cupmtr(SIDE, UPLO, TRANS, M, N, AP, TAU, C, LDC, WORK, INFO)
CUPMTR
Definition: cupmtr.f:152
subroutine cswap(N, CX, INCX, CY, INCY)
CSWAP
Definition: cswap.f:52
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:69
subroutine ssterf(N, D, E, INFO)
SSTERF
Definition: ssterf.f:88
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY
Definition: scopy.f:53
subroutine csscal(N, SA, CX, INCX)
CSSCAL
Definition: csscal.f:54

Here is the call graph for this function:

Here is the caller graph for this function: