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

Go to the source code of this file.

Functions/Subroutines

subroutine zhpevx (JOBZ, RANGE, UPLO, N, AP, VL, VU, IL, IU, ABSTOL, M, W, Z, LDZ, WORK, RWORK, IWORK, IFAIL, INFO)
  ZHPEVX computes the eigenvalues and, optionally, the left and/or right eigenvectors for OTHER matrices More...
 

Function/Subroutine Documentation

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
[in]VU
          VU is DOUBLE PRECISION
          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 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
November 2011

Definition at line 235 of file zhpevx.f.

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

Here is the call graph for this function:

Here is the caller graph for this function: