LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ zheevx()

subroutine zheevx ( character  jobz,
character  range,
character  uplo,
integer  n,
complex*16, dimension( lda, * )  a,
integer  lda,
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,
integer  lwork,
double precision, dimension( * )  rwork,
integer, dimension( * )  iwork,
integer, dimension( * )  ifail,
integer  info 
)

ZHEEVX computes the eigenvalues and, optionally, the left and/or right eigenvectors for HE matrices

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

Purpose:
 ZHEEVX 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*16 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 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 A 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)
          On normal exit, the first M elements contain 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 (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 ZHETRD and for
          ZUNMTR 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 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 256 of file zheevx.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 DOUBLE PRECISION ABSTOL, VL, VU
268* ..
269* .. Array Arguments ..
270 INTEGER IFAIL( * ), IWORK( * )
271 DOUBLE PRECISION RWORK( * ), W( * )
272 COMPLEX*16 A( LDA, * ), WORK( * ), Z( LDZ, * )
273* ..
274*
275* =====================================================================
276*
277* .. Parameters ..
278 DOUBLE PRECISION ZERO, ONE
279 parameter( zero = 0.0d+0, one = 1.0d+0 )
280 COMPLEX*16 CONE
281 parameter( cone = ( 1.0d+0, 0.0d+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 DOUBLE PRECISION ABSTLL, ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN,
292 $ SIGMA, SMLNUM, TMP1, VLL, VUU
293* ..
294* .. External Functions ..
295 LOGICAL LSAME
296 INTEGER ILAENV
297 DOUBLE PRECISION DLAMCH, ZLANHE
298 EXTERNAL lsame, ilaenv, dlamch, zlanhe
299* ..
300* .. External Subroutines ..
301 EXTERNAL dcopy, dscal, dstebz, dsterf, xerbla, zdscal,
303 $ zunmtr
304* ..
305* .. Intrinsic Functions ..
306 INTRINSIC dble, 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, 'ZHETRD', uplo, n, -1, -1, -1 )
355 nb = max( nb, ilaenv( 1, 'ZUNMTR', 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( 'ZHEEVX', -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 ) = dble( a( 1, 1 ) )
382 ELSE IF( valeig ) THEN
383 IF( vl.LT.dble( a( 1, 1 ) ) .AND. vu.GE.dble( a( 1, 1 ) ) )
384 $ THEN
385 m = 1
386 w( 1 ) = dble( 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 = dlamch( 'Safe minimum' )
397 eps = dlamch( '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 = zlanhe( '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 zdscal( n-j+1, sigma, a( j, j ), 1 )
423 10 CONTINUE
424 ELSE
425 DO 20 j = 1, n
426 CALL zdscal( 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 ZHETRD 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 zhetrd( 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 DSTERF or ZUNGTR and ZSTEQR. If this fails for
450* some eigenvalue, then try DSTEBZ.
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 dcopy( n, rwork( indd ), 1, w, 1 )
460 indee = indrwk + 2*n
461 IF( .NOT.wantz ) THEN
462 CALL dcopy( n-1, rwork( inde ), 1, rwork( indee ), 1 )
463 CALL dsterf( n, w, rwork( indee ), info )
464 ELSE
465 CALL zlacpy( 'A', n, n, a, lda, z, ldz )
466 CALL zungtr( uplo, n, z, ldz, work( indtau ),
467 $ work( indwrk ), llwork, iinfo )
468 CALL dcopy( n-1, rwork( inde ), 1, rwork( indee ), 1 )
469 CALL zsteqr( 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 DSTEBZ and, if eigenvectors are desired, ZSTEIN.
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 dstebz( 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 zstein( 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 ZSTEIN.
506*
507 CALL zunmtr( '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 dscal( 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 zswap( 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 ZHEEVX
560*
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine dcopy(n, dx, incx, dy, incy)
DCOPY
Definition dcopy.f:82
subroutine zhetrd(uplo, n, a, lda, d, e, tau, work, lwork, info)
ZHETRD
Definition zhetrd.f:192
integer function ilaenv(ispec, name, opts, n1, n2, n3, n4)
ILAENV
Definition ilaenv.f:162
subroutine zlacpy(uplo, m, n, a, lda, b, ldb)
ZLACPY copies all or part of one two-dimensional array to another.
Definition zlacpy.f:103
double precision function dlamch(cmach)
DLAMCH
Definition dlamch.f:69
double precision function zlanhe(norm, uplo, n, a, lda, work)
ZLANHE returns the value of the 1-norm, or the Frobenius norm, or the infinity norm,...
Definition zlanhe.f:124
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
subroutine dscal(n, da, dx, incx)
DSCAL
Definition dscal.f:79
subroutine zdscal(n, da, zx, incx)
ZDSCAL
Definition zdscal.f:78
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 zstein(n, d, e, m, w, iblock, isplit, z, ldz, work, iwork, ifail, info)
ZSTEIN
Definition zstein.f:182
subroutine zsteqr(compz, n, d, e, z, ldz, work, info)
ZSTEQR
Definition zsteqr.f:132
subroutine dsterf(n, d, e, info)
DSTERF
Definition dsterf.f:86
subroutine zswap(n, zx, incx, zy, incy)
ZSWAP
Definition zswap.f:81
subroutine zungtr(uplo, n, a, lda, tau, work, lwork, info)
ZUNGTR
Definition zungtr.f:123
subroutine zunmtr(side, uplo, trans, m, n, a, lda, tau, c, ldc, work, lwork, info)
ZUNMTR
Definition zunmtr.f:171
Here is the call graph for this function:
Here is the caller graph for this function: