LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
subroutine zhpevx ( character  JOBZ,
character  RANGE,
character  UPLO,
integer  N,
complex*16, dimension( * )  AP,
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 
)

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

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

Purpose:
 ZHPEVX 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*16 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 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 AP 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)
          If INFO = 0, 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 (2*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.
Date
June 2016

Definition at line 242 of file zhpevx.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  DOUBLE PRECISION abstol, vl, vu
252 * ..
253 * .. Array Arguments ..
254  INTEGER ifail( * ), iwork( * )
255  DOUBLE PRECISION rwork( * ), w( * )
256  COMPLEX*16 ap( * ), work( * ), z( ldz, * )
257 * ..
258 *
259 * =====================================================================
260 *
261 * .. Parameters ..
262  DOUBLE PRECISION zero, one
263  parameter ( zero = 0.0d0, one = 1.0d0 )
264  COMPLEX*16 cone
265  parameter ( cone = ( 1.0d0, 0.0d0 ) )
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  DOUBLE PRECISION abstll, anrm, bignum, eps, rmax, rmin, safmin,
274  $ sigma, smlnum, tmp1, vll, vuu
275 * ..
276 * .. External Functions ..
277  LOGICAL lsame
278  DOUBLE PRECISION dlamch, zlanhp
279  EXTERNAL lsame, dlamch, zlanhp
280 * ..
281 * .. External Subroutines ..
282  EXTERNAL dcopy, dscal, dstebz, dsterf, xerbla, zdscal,
284 * ..
285 * .. Intrinsic Functions ..
286  INTRINSIC dble, max, min, 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( 'ZHPEVX', -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.dble( ap( 1 ) ) .AND. vu.GE.dble( 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 = dlamch( 'Safe minimum' )
353  eps = dlamch( '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  END IF
370  anrm = zlanhp( '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 zdscal( ( 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 ZHPTRD 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 zhptrd( 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 DSTERF or ZUPGTR and ZSTEQR. If this fails
400 * for some eigenvalue, then try DSTEBZ.
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 dcopy( n, rwork( indd ), 1, w, 1 )
410  indee = indrwk + 2*n
411  IF( .NOT.wantz ) THEN
412  CALL dcopy( n-1, rwork( inde ), 1, rwork( indee ), 1 )
413  CALL dsterf( n, w, rwork( indee ), info )
414  ELSE
415  CALL zupgtr( uplo, n, ap, work( indtau ), z, ldz,
416  $ work( indwrk ), iinfo )
417  CALL dcopy( n-1, rwork( inde ), 1, rwork( indee ), 1 )
418  CALL zsteqr( 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 DSTEBZ and, if eigenvectors are desired, ZSTEIN.
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 dstebz( 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 zstein( 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 ZSTEIN.
455 *
456  indwrk = indtau + n
457  CALL zupmtr( '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 dscal( 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 zswap( 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 ZHPEVX
506 *
subroutine dsterf(N, D, E, INFO)
DSTERF
Definition: dsterf.f:88
subroutine dstebz(RANGE, ORDER, N, VL, VU, IL, IU, ABSTOL, D, E, M, NSPLIT, W, IBLOCK, ISPLIT, WORK, IWORK, INFO)
DSTEBZ
Definition: dstebz.f:275
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
Definition: dcopy.f:53
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:65
double precision function zlanhp(NORM, UPLO, N, AP, WORK)
ZLANHP 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: zlanhp.f:119
subroutine zswap(N, ZX, INCX, ZY, INCY)
ZSWAP
Definition: zswap.f:52
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine zsteqr(COMPZ, N, D, E, Z, LDZ, WORK, INFO)
ZSTEQR
Definition: zsteqr.f:134
subroutine zstein(N, D, E, M, W, IBLOCK, ISPLIT, Z, LDZ, WORK, IWORK, IFAIL, INFO)
ZSTEIN
Definition: zstein.f:184
subroutine zhptrd(UPLO, N, AP, D, E, TAU, INFO)
ZHPTRD
Definition: zhptrd.f:153
subroutine dscal(N, DA, DX, INCX)
DSCAL
Definition: dscal.f:55
subroutine zdscal(N, DA, ZX, INCX)
ZDSCAL
Definition: zdscal.f:54
subroutine zupmtr(SIDE, UPLO, TRANS, M, N, AP, TAU, C, LDC, WORK, INFO)
ZUPMTR
Definition: zupmtr.f:152
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
subroutine zupgtr(UPLO, N, AP, TAU, Q, LDQ, WORK, INFO)
ZUPGTR
Definition: zupgtr.f:116

Here is the call graph for this function:

Here is the caller graph for this function: