LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

◆ zhpevx()

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.

Definition at line 237 of file zhpevx.f.

240 *
241 * -- LAPACK driver routine --
242 * -- LAPACK is a software package provided by Univ. of Tennessee, --
243 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
244 *
245 * .. Scalar Arguments ..
246  CHARACTER JOBZ, RANGE, UPLO
247  INTEGER IL, INFO, IU, LDZ, M, N
248  DOUBLE PRECISION ABSTOL, VL, VU
249 * ..
250 * .. Array Arguments ..
251  INTEGER IFAIL( * ), IWORK( * )
252  DOUBLE PRECISION RWORK( * ), W( * )
253  COMPLEX*16 AP( * ), WORK( * ), Z( LDZ, * )
254 * ..
255 *
256 * =====================================================================
257 *
258 * .. Parameters ..
259  DOUBLE PRECISION ZERO, ONE
260  parameter( zero = 0.0d0, one = 1.0d0 )
261  COMPLEX*16 CONE
262  parameter( cone = ( 1.0d0, 0.0d0 ) )
263 * ..
264 * .. Local Scalars ..
265  LOGICAL ALLEIG, INDEIG, TEST, VALEIG, WANTZ
266  CHARACTER ORDER
267  INTEGER I, IINFO, IMAX, INDD, INDE, INDEE, INDIBL,
268  $ INDISP, INDIWK, INDRWK, INDTAU, INDWRK, ISCALE,
269  $ ITMP1, J, JJ, NSPLIT
270  DOUBLE PRECISION ABSTLL, ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN,
271  $ SIGMA, SMLNUM, TMP1, VLL, VUU
272 * ..
273 * .. External Functions ..
274  LOGICAL LSAME
275  DOUBLE PRECISION DLAMCH, ZLANHP
276  EXTERNAL lsame, dlamch, zlanhp
277 * ..
278 * .. External Subroutines ..
279  EXTERNAL dcopy, dscal, dstebz, dsterf, xerbla, zdscal,
281 * ..
282 * .. Intrinsic Functions ..
283  INTRINSIC dble, max, min, sqrt
284 * ..
285 * .. Executable Statements ..
286 *
287 * Test the input parameters.
288 *
289  wantz = lsame( jobz, 'V' )
290  alleig = lsame( range, 'A' )
291  valeig = lsame( range, 'V' )
292  indeig = lsame( range, 'I' )
293 *
294  info = 0
295  IF( .NOT.( wantz .OR. lsame( jobz, 'N' ) ) ) THEN
296  info = -1
297  ELSE IF( .NOT.( alleig .OR. valeig .OR. indeig ) ) THEN
298  info = -2
299  ELSE IF( .NOT.( lsame( uplo, 'L' ) .OR. lsame( uplo, 'U' ) ) )
300  $ THEN
301  info = -3
302  ELSE IF( n.LT.0 ) THEN
303  info = -4
304  ELSE
305  IF( valeig ) THEN
306  IF( n.GT.0 .AND. vu.LE.vl )
307  $ info = -7
308  ELSE IF( indeig ) THEN
309  IF( il.LT.1 .OR. il.GT.max( 1, n ) ) THEN
310  info = -8
311  ELSE IF( iu.LT.min( n, il ) .OR. iu.GT.n ) THEN
312  info = -9
313  END IF
314  END IF
315  END IF
316  IF( info.EQ.0 ) THEN
317  IF( ldz.LT.1 .OR. ( wantz .AND. ldz.LT.n ) )
318  $ info = -14
319  END IF
320 *
321  IF( info.NE.0 ) THEN
322  CALL xerbla( 'ZHPEVX', -info )
323  RETURN
324  END IF
325 *
326 * Quick return if possible
327 *
328  m = 0
329  IF( n.EQ.0 )
330  $ RETURN
331 *
332  IF( n.EQ.1 ) THEN
333  IF( alleig .OR. indeig ) THEN
334  m = 1
335  w( 1 ) = dble( ap( 1 ) )
336  ELSE
337  IF( vl.LT.dble( ap( 1 ) ) .AND. vu.GE.dble( ap( 1 ) ) ) THEN
338  m = 1
339  w( 1 ) = dble( ap( 1 ) )
340  END IF
341  END IF
342  IF( wantz )
343  $ z( 1, 1 ) = cone
344  RETURN
345  END IF
346 *
347 * Get machine constants.
348 *
349  safmin = dlamch( 'Safe minimum' )
350  eps = dlamch( 'Precision' )
351  smlnum = safmin / eps
352  bignum = one / smlnum
353  rmin = sqrt( smlnum )
354  rmax = min( sqrt( bignum ), one / sqrt( sqrt( safmin ) ) )
355 *
356 * Scale matrix to allowable range, if necessary.
357 *
358  iscale = 0
359  abstll = abstol
360  IF( valeig ) THEN
361  vll = vl
362  vuu = vu
363  ELSE
364  vll = zero
365  vuu = zero
366  END IF
367  anrm = zlanhp( 'M', uplo, n, ap, rwork )
368  IF( anrm.GT.zero .AND. anrm.LT.rmin ) THEN
369  iscale = 1
370  sigma = rmin / anrm
371  ELSE IF( anrm.GT.rmax ) THEN
372  iscale = 1
373  sigma = rmax / anrm
374  END IF
375  IF( iscale.EQ.1 ) THEN
376  CALL zdscal( ( n*( n+1 ) ) / 2, sigma, ap, 1 )
377  IF( abstol.GT.0 )
378  $ abstll = abstol*sigma
379  IF( valeig ) THEN
380  vll = vl*sigma
381  vuu = vu*sigma
382  END IF
383  END IF
384 *
385 * Call ZHPTRD to reduce Hermitian packed matrix to tridiagonal form.
386 *
387  indd = 1
388  inde = indd + n
389  indrwk = inde + n
390  indtau = 1
391  indwrk = indtau + n
392  CALL zhptrd( uplo, n, ap, rwork( indd ), rwork( inde ),
393  $ work( indtau ), iinfo )
394 *
395 * If all eigenvalues are desired and ABSTOL is less than or equal
396 * to zero, then call DSTERF or ZUPGTR and ZSTEQR. If this fails
397 * for some eigenvalue, then try DSTEBZ.
398 *
399  test = .false.
400  IF (indeig) THEN
401  IF (il.EQ.1 .AND. iu.EQ.n) THEN
402  test = .true.
403  END IF
404  END IF
405  IF ((alleig .OR. test) .AND. (abstol.LE.zero)) THEN
406  CALL dcopy( n, rwork( indd ), 1, w, 1 )
407  indee = indrwk + 2*n
408  IF( .NOT.wantz ) THEN
409  CALL dcopy( n-1, rwork( inde ), 1, rwork( indee ), 1 )
410  CALL dsterf( n, w, rwork( indee ), info )
411  ELSE
412  CALL zupgtr( uplo, n, ap, work( indtau ), z, ldz,
413  $ work( indwrk ), iinfo )
414  CALL dcopy( n-1, rwork( inde ), 1, rwork( indee ), 1 )
415  CALL zsteqr( jobz, n, w, rwork( indee ), z, ldz,
416  $ rwork( indrwk ), info )
417  IF( info.EQ.0 ) THEN
418  DO 10 i = 1, n
419  ifail( i ) = 0
420  10 CONTINUE
421  END IF
422  END IF
423  IF( info.EQ.0 ) THEN
424  m = n
425  GO TO 20
426  END IF
427  info = 0
428  END IF
429 *
430 * Otherwise, call DSTEBZ and, if eigenvectors are desired, ZSTEIN.
431 *
432  IF( wantz ) THEN
433  order = 'B'
434  ELSE
435  order = 'E'
436  END IF
437  indibl = 1
438  indisp = indibl + n
439  indiwk = indisp + n
440  CALL dstebz( range, order, n, vll, vuu, il, iu, abstll,
441  $ rwork( indd ), rwork( inde ), m, nsplit, w,
442  $ iwork( indibl ), iwork( indisp ), rwork( indrwk ),
443  $ iwork( indiwk ), info )
444 *
445  IF( wantz ) THEN
446  CALL zstein( n, rwork( indd ), rwork( inde ), m, w,
447  $ iwork( indibl ), iwork( indisp ), z, ldz,
448  $ rwork( indrwk ), iwork( indiwk ), ifail, info )
449 *
450 * Apply unitary matrix used in reduction to tridiagonal
451 * form to eigenvectors returned by ZSTEIN.
452 *
453  indwrk = indtau + n
454  CALL zupmtr( 'L', uplo, 'N', n, m, ap, work( indtau ), z, ldz,
455  $ work( indwrk ), iinfo )
456  END IF
457 *
458 * If matrix was scaled, then rescale eigenvalues appropriately.
459 *
460  20 CONTINUE
461  IF( iscale.EQ.1 ) THEN
462  IF( info.EQ.0 ) THEN
463  imax = m
464  ELSE
465  imax = info - 1
466  END IF
467  CALL dscal( imax, one / sigma, w, 1 )
468  END IF
469 *
470 * If eigenvalues are not in order, then sort them, along with
471 * eigenvectors.
472 *
473  IF( wantz ) THEN
474  DO 40 j = 1, m - 1
475  i = 0
476  tmp1 = w( j )
477  DO 30 jj = j + 1, m
478  IF( w( jj ).LT.tmp1 ) THEN
479  i = jj
480  tmp1 = w( jj )
481  END IF
482  30 CONTINUE
483 *
484  IF( i.NE.0 ) THEN
485  itmp1 = iwork( indibl+i-1 )
486  w( i ) = w( j )
487  iwork( indibl+i-1 ) = iwork( indibl+j-1 )
488  w( j ) = tmp1
489  iwork( indibl+j-1 ) = itmp1
490  CALL zswap( n, z( 1, i ), 1, z( 1, j ), 1 )
491  IF( info.NE.0 ) THEN
492  itmp1 = ifail( i )
493  ifail( i ) = ifail( j )
494  ifail( j ) = itmp1
495  END IF
496  END IF
497  40 CONTINUE
498  END IF
499 *
500  RETURN
501 *
502 * End of ZHPEVX
503 *
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:69
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
subroutine dstebz(RANGE, ORDER, N, VL, VU, IL, IU, ABSTOL, D, E, M, NSPLIT, W, IBLOCK, ISPLIT, WORK, IWORK, INFO)
DSTEBZ
Definition: dstebz.f:273
subroutine dsterf(N, D, E, INFO)
DSTERF
Definition: dsterf.f:86
subroutine zswap(N, ZX, INCX, ZY, INCY)
ZSWAP
Definition: zswap.f:81
subroutine zdscal(N, DA, ZX, INCX)
ZDSCAL
Definition: zdscal.f:78
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,...
Definition: zlanhp.f:117
subroutine zsteqr(COMPZ, N, D, E, Z, LDZ, WORK, INFO)
ZSTEQR
Definition: zsteqr.f:132
subroutine zupmtr(SIDE, UPLO, TRANS, M, N, AP, TAU, C, LDC, WORK, INFO)
ZUPMTR
Definition: zupmtr.f:150
subroutine zhptrd(UPLO, N, AP, D, E, TAU, INFO)
ZHPTRD
Definition: zhptrd.f:151
subroutine zstein(N, D, E, M, W, IBLOCK, ISPLIT, Z, LDZ, WORK, IWORK, IFAIL, INFO)
ZSTEIN
Definition: zstein.f:182
subroutine zupgtr(UPLO, N, AP, TAU, Q, LDQ, WORK, INFO)
ZUPGTR
Definition: zupgtr.f:114
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
Definition: dcopy.f:82
subroutine dscal(N, DA, DX, INCX)
DSCAL
Definition: dscal.f:79
Here is the call graph for this function:
Here is the caller graph for this function: