LAPACK  3.10.1 LAPACK: Linear Algebra PACKage

## ◆ zheevr_2stage()

 subroutine zheevr_2stage ( 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, integer, dimension( * ) ISUPPZ, complex*16, dimension( * ) WORK, integer LWORK, double precision, dimension( * ) RWORK, integer LRWORK, integer, dimension( * ) IWORK, integer LIWORK, integer INFO )

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

Purpose:
``` ZHEEVR_2STAGE computes selected eigenvalues and, optionally, eigenvectors
of a complex Hermitian matrix A using the 2stage technique for
the reduction to tridiagonal.  Eigenvalues and eigenvectors can
be selected by specifying either a range of values or a range of
indices for the desired eigenvalues.

ZHEEVR_2STAGE first reduces the matrix A to tridiagonal form T with a call
to ZHETRD.  Then, whenever possible, ZHEEVR_2STAGE calls ZSTEMR to compute
eigenspectrum using Relatively Robust Representations.  ZSTEMR
computes eigenvalues by the dqds algorithm, while orthogonal
eigenvectors are computed from various "good" L D L^T representations
(also known as Relatively Robust Representations). Gram-Schmidt
orthogonalization is avoided as far as possible. More specifically,
the various steps of the algorithm are as follows.

For each unreduced block (submatrix) of T,
(a) Compute T - sigma I  = L D L^T, so that L and D
define all the wanted eigenvalues to high relative accuracy.
This means that small relative changes in the entries of D and L
cause only small relative changes in the eigenvalues and
eigenvectors. The standard (unfactored) representation of the
tridiagonal matrix T does not have this property in general.
(b) Compute the eigenvalues to suitable accuracy.
If the eigenvectors are desired, the algorithm attains full
accuracy of the computed eigenvalues only right before
the corresponding vectors have to be computed, see steps c) and d).
(c) For each cluster of close eigenvalues, select a new
shift close to the cluster, find a new factorization, and refine
the shifted eigenvalues to suitable accuracy.
(d) For each eigenvalue with a large enough relative separation compute
the corresponding eigenvector by forming a rank revealing twisted
factorization. Go back to (c) for any clusters that remain.

The desired accuracy of the output can be specified by the input
parameter ABSTOL.

For more details, see ZSTEMR's documentation and:
- Inderjit S. Dhillon and Beresford N. Parlett: "Multiple representations
to compute orthogonal eigenvectors of symmetric tridiagonal matrices,"
Linear Algebra and its Applications, 387(1), pp. 1-28, August 2004.
- Inderjit Dhillon and Beresford Parlett: "Orthogonal Eigenvectors and
Relative Gaps," SIAM Journal on Matrix Analysis and Applications, Vol. 25,
2004.  Also LAPACK Working Note 154.
- Inderjit Dhillon: "A new O(n^2) algorithm for the symmetric
tridiagonal eigenvalue/eigenvector problem",
Computer Science Division Technical Report No. UCB/CSD-97-971,
UC Berkeley, May 1997.

Note 1 : ZHEEVR_2STAGE calls ZSTEMR when the full spectrum is requested
on machines which conform to the ieee-754 floating point standard.
ZHEEVR_2STAGE calls DSTEBZ and ZSTEIN on non-ieee machines and
when partial spectrum requests are made.

Normal execution of ZSTEMR may create NaNs and infinities and
hence may abort due to a floating point exception in environments
which do not handle NaNs and infinities in the ieee standard default
manner.```
Parameters
 [in] JOBZ ``` JOBZ is CHARACTER*1 = 'N': Compute eigenvalues only; = 'V': Compute eigenvalues and eigenvectors. Not available in this release.``` [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. For RANGE = 'V' or 'I' and IU - IL < N - 1, DSTEBZ and ZSTEIN are called``` [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. See "Computing Small Singular Values of Bidiagonal Matrices with Guaranteed High Relative Accuracy," by Demmel and Kahan, LAPACK Working Note #3. If high relative accuracy is important, set ABSTOL to DLAMCH( 'Safe minimum' ). Doing so will guarantee that eigenvalues are computed to high relative accuracy when possible in future releases. The current code does not make any guarantees about high relative accuracy, but future releases will. See J. Barlow and J. Demmel, "Computing Accurate Eigensystems of Scaled Diagonally Dominant Matrices", LAPACK Working Note #7, for a discussion of which matrices define their eigenvalues to high relative accuracy.``` [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) 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 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] ISUPPZ ``` ISUPPZ is INTEGER array, dimension ( 2*max(1,M) ) The support of the eigenvectors in Z, i.e., the indices indicating the nonzero elements in Z. The i-th eigenvector is nonzero only in elements ISUPPZ( 2*i-1 ) through ISUPPZ( 2*i ). This is an output of ZSTEMR (tridiagonal matrix). The support of the eigenvectors of A is typically 1:N because of the unitary transformations applied by ZUNMTR. Implemented only for RANGE = 'A' or 'I' and IU - IL = N - 1``` [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 dimension of the array WORK. If JOBZ = 'N' and N > 1, LWORK must be queried. LWORK = MAX(1, 26*N, dimension) where dimension = max(stage1,stage2) + (KD+1)*N + N = N*KD + N*max(KD+1,FACTOPTNB) + max(2*KD*KD, KD*NTHREADS) + (KD+1)*N + N where KD is the blocking size of the reduction, FACTOPTNB is the blocking used by the QR or LQ algorithm, usually FACTOPTNB=128 is a good choice NTHREADS is the number of threads used when openMP compilation is enabled, otherwise =1. If JOBZ = 'V' and N > 1, LWORK must be queried. Not yet available If LWORK = -1, then a workspace query is assumed; the routine only calculates the optimal sizes of the WORK, RWORK and IWORK arrays, returns these values as the first entries of the WORK, RWORK and IWORK arrays, and no error message related to LWORK or LRWORK or LIWORK is issued by XERBLA.``` [out] RWORK ``` RWORK is DOUBLE PRECISION array, dimension (MAX(1,LRWORK)) On exit, if INFO = 0, RWORK(1) returns the optimal (and minimal) LRWORK.``` [in] LRWORK ``` LRWORK is INTEGER The length of the array RWORK. LRWORK >= max(1,24*N). If LRWORK = -1, then a workspace query is assumed; the routine only calculates the optimal sizes of the WORK, RWORK and IWORK arrays, returns these values as the first entries of the WORK, RWORK and IWORK arrays, and no error message related to LWORK or LRWORK or LIWORK is issued by XERBLA.``` [out] IWORK ``` IWORK is INTEGER array, dimension (MAX(1,LIWORK)) On exit, if INFO = 0, IWORK(1) returns the optimal (and minimal) LIWORK.``` [in] LIWORK ``` LIWORK is INTEGER The dimension of the array IWORK. LIWORK >= max(1,10*N). If LIWORK = -1, then a workspace query is assumed; the routine only calculates the optimal sizes of the WORK, RWORK and IWORK arrays, returns these values as the first entries of the WORK, RWORK and IWORK arrays, and no error message related to LWORK or LRWORK or LIWORK is issued by XERBLA.``` [out] INFO ``` INFO is INTEGER = 0: successful exit < 0: if INFO = -i, the i-th argument had an illegal value > 0: Internal error```
Contributors:
```Inderjit Dhillon, IBM Almaden, USA \n
Osni Marques, LBNL/NERSC, USA \n
Ken Stanley, Computer Science Division, University of
California at Berkeley, USA \n
Jason Riedy, Computer Science Division, University of
California at Berkeley, USA \n
```
Further Details:
```  All details about the 2stage techniques are available in:

Azzam Haidar, Hatem Ltaief, and Jack Dongarra.
Parallel reduction to condensed forms for symmetric eigenvalue problems
using aggregated fine-grained and memory-aware kernels. In Proceedings
of 2011 International Conference for High Performance Computing,
Networking, Storage and Analysis (SC '11), New York, NY, USA,
Article 8 , 11 pages.
http://doi.acm.org/10.1145/2063384.2063394

A. Haidar, J. Kurzak, P. Luszczek, 2013.
An improved parallel singular value algorithm and its implementation
for multicore hardware, In Proceedings of 2013 International Conference
for High Performance Computing, Networking, Storage and Analysis (SC '13).
Article 90, 12 pages.
http://doi.acm.org/10.1145/2503210.2503292

A. Haidar, R. Solca, S. Tomov, T. Schulthess and J. Dongarra.
A novel hybrid CPU-GPU generalized eigensolver for electronic structure
calculations based on fine-grained memory aware tasks.
International Journal of High Performance Computing Applications.
Volume 28 Issue 2, Pages 196-209, May 2014.
http://hpc.sagepub.com/content/28/2/196 ```

Definition at line 402 of file zheevr_2stage.f.

406 *
407  IMPLICIT NONE
408 *
409 * -- LAPACK driver routine --
410 * -- LAPACK is a software package provided by Univ. of Tennessee, --
411 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
412 *
413 * .. Scalar Arguments ..
414  CHARACTER JOBZ, RANGE, UPLO
415  INTEGER IL, INFO, IU, LDA, LDZ, LIWORK, LRWORK, LWORK,
416  \$ M, N
417  DOUBLE PRECISION ABSTOL, VL, VU
418 * ..
419 * .. Array Arguments ..
420  INTEGER ISUPPZ( * ), IWORK( * )
421  DOUBLE PRECISION RWORK( * ), W( * )
422  COMPLEX*16 A( LDA, * ), WORK( * ), Z( LDZ, * )
423 * ..
424 *
425 * =====================================================================
426 *
427 * .. Parameters ..
428  DOUBLE PRECISION ZERO, ONE, TWO
429  parameter( zero = 0.0d+0, one = 1.0d+0, two = 2.0d+0 )
430 * ..
431 * .. Local Scalars ..
432  LOGICAL ALLEIG, INDEIG, LOWER, LQUERY, TEST, VALEIG,
433  \$ WANTZ, TRYRAC
434  CHARACTER ORDER
435  INTEGER I, IEEEOK, IINFO, IMAX, INDIBL, INDIFL, INDISP,
436  \$ INDIWO, INDRD, INDRDD, INDRE, INDREE, INDRWK,
437  \$ INDTAU, INDWK, INDWKN, ISCALE, ITMP1, J, JJ,
438  \$ LIWMIN, LLWORK, LLRWORK, LLWRKN, LRWMIN,
439  \$ LWMIN, NSPLIT, LHTRD, LWTRD, KD, IB, INDHOUS
440  DOUBLE PRECISION ABSTLL, ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN,
441  \$ SIGMA, SMLNUM, TMP1, VLL, VUU
442 * ..
443 * .. External Functions ..
444  LOGICAL LSAME
445  INTEGER ILAENV, ILAENV2STAGE
446  DOUBLE PRECISION DLAMCH, ZLANSY
447  EXTERNAL lsame, dlamch, zlansy, ilaenv, ilaenv2stage
448 * ..
449 * .. External Subroutines ..
450  EXTERNAL dcopy, dscal, dstebz, dsterf, xerbla, zdscal,
452 * ..
453 * .. Intrinsic Functions ..
454  INTRINSIC dble, max, min, sqrt
455 * ..
456 * .. Executable Statements ..
457 *
458 * Test the input parameters.
459 *
460  ieeeok = ilaenv( 10, 'ZHEEVR', 'N', 1, 2, 3, 4 )
461 *
462  lower = lsame( uplo, 'L' )
463  wantz = lsame( jobz, 'V' )
464  alleig = lsame( range, 'A' )
465  valeig = lsame( range, 'V' )
466  indeig = lsame( range, 'I' )
467 *
468  lquery = ( ( lwork.EQ.-1 ) .OR. ( lrwork.EQ.-1 ) .OR.
469  \$ ( liwork.EQ.-1 ) )
470 *
471  kd = ilaenv2stage( 1, 'ZHETRD_2STAGE', jobz, n, -1, -1, -1 )
472  ib = ilaenv2stage( 2, 'ZHETRD_2STAGE', jobz, n, kd, -1, -1 )
473  lhtrd = ilaenv2stage( 3, 'ZHETRD_2STAGE', jobz, n, kd, ib, -1 )
474  lwtrd = ilaenv2stage( 4, 'ZHETRD_2STAGE', jobz, n, kd, ib, -1 )
475  lwmin = n + lhtrd + lwtrd
476  lrwmin = max( 1, 24*n )
477  liwmin = max( 1, 10*n )
478 *
479  info = 0
480  IF( .NOT.( lsame( jobz, 'N' ) ) ) THEN
481  info = -1
482  ELSE IF( .NOT.( alleig .OR. valeig .OR. indeig ) ) THEN
483  info = -2
484  ELSE IF( .NOT.( lower .OR. lsame( uplo, 'U' ) ) ) THEN
485  info = -3
486  ELSE IF( n.LT.0 ) THEN
487  info = -4
488  ELSE IF( lda.LT.max( 1, n ) ) THEN
489  info = -6
490  ELSE
491  IF( valeig ) THEN
492  IF( n.GT.0 .AND. vu.LE.vl )
493  \$ info = -8
494  ELSE IF( indeig ) THEN
495  IF( il.LT.1 .OR. il.GT.max( 1, n ) ) THEN
496  info = -9
497  ELSE IF( iu.LT.min( n, il ) .OR. iu.GT.n ) THEN
498  info = -10
499  END IF
500  END IF
501  END IF
502  IF( info.EQ.0 ) THEN
503  IF( ldz.LT.1 .OR. ( wantz .AND. ldz.LT.n ) ) THEN
504  info = -15
505  END IF
506  END IF
507 *
508  IF( info.EQ.0 ) THEN
509  work( 1 ) = lwmin
510  rwork( 1 ) = lrwmin
511  iwork( 1 ) = liwmin
512 *
513  IF( lwork.LT.lwmin .AND. .NOT.lquery ) THEN
514  info = -18
515  ELSE IF( lrwork.LT.lrwmin .AND. .NOT.lquery ) THEN
516  info = -20
517  ELSE IF( liwork.LT.liwmin .AND. .NOT.lquery ) THEN
518  info = -22
519  END IF
520  END IF
521 *
522  IF( info.NE.0 ) THEN
523  CALL xerbla( 'ZHEEVR_2STAGE', -info )
524  RETURN
525  ELSE IF( lquery ) THEN
526  RETURN
527  END IF
528 *
529 * Quick return if possible
530 *
531  m = 0
532  IF( n.EQ.0 ) THEN
533  work( 1 ) = 1
534  RETURN
535  END IF
536 *
537  IF( n.EQ.1 ) THEN
538  work( 1 ) = 2
539  IF( alleig .OR. indeig ) THEN
540  m = 1
541  w( 1 ) = dble( a( 1, 1 ) )
542  ELSE
543  IF( vl.LT.dble( a( 1, 1 ) ) .AND. vu.GE.dble( a( 1, 1 ) ) )
544  \$ THEN
545  m = 1
546  w( 1 ) = dble( a( 1, 1 ) )
547  END IF
548  END IF
549  IF( wantz ) THEN
550  z( 1, 1 ) = one
551  isuppz( 1 ) = 1
552  isuppz( 2 ) = 1
553  END IF
554  RETURN
555  END IF
556 *
557 * Get machine constants.
558 *
559  safmin = dlamch( 'Safe minimum' )
560  eps = dlamch( 'Precision' )
561  smlnum = safmin / eps
562  bignum = one / smlnum
563  rmin = sqrt( smlnum )
564  rmax = min( sqrt( bignum ), one / sqrt( sqrt( safmin ) ) )
565 *
566 * Scale matrix to allowable range, if necessary.
567 *
568  iscale = 0
569  abstll = abstol
570  IF (valeig) THEN
571  vll = vl
572  vuu = vu
573  END IF
574  anrm = zlansy( 'M', uplo, n, a, lda, rwork )
575  IF( anrm.GT.zero .AND. anrm.LT.rmin ) THEN
576  iscale = 1
577  sigma = rmin / anrm
578  ELSE IF( anrm.GT.rmax ) THEN
579  iscale = 1
580  sigma = rmax / anrm
581  END IF
582  IF( iscale.EQ.1 ) THEN
583  IF( lower ) THEN
584  DO 10 j = 1, n
585  CALL zdscal( n-j+1, sigma, a( j, j ), 1 )
586  10 CONTINUE
587  ELSE
588  DO 20 j = 1, n
589  CALL zdscal( j, sigma, a( 1, j ), 1 )
590  20 CONTINUE
591  END IF
592  IF( abstol.GT.0 )
593  \$ abstll = abstol*sigma
594  IF( valeig ) THEN
595  vll = vl*sigma
596  vuu = vu*sigma
597  END IF
598  END IF
599
600 * Initialize indices into workspaces. Note: The IWORK indices are
601 * used only if DSTERF or ZSTEMR fail.
602
603 * WORK(INDTAU:INDTAU+N-1) stores the complex scalar factors of the
604 * elementary reflectors used in ZHETRD.
605  indtau = 1
606 * INDWK is the starting offset of the remaining complex workspace,
607 * and LLWORK is the remaining complex workspace size.
608  indhous = indtau + n
609  indwk = indhous + lhtrd
610  llwork = lwork - indwk + 1
611
612 * RWORK(INDRD:INDRD+N-1) stores the real tridiagonal's diagonal
613 * entries.
614  indrd = 1
615 * RWORK(INDRE:INDRE+N-1) stores the off-diagonal entries of the
616 * tridiagonal matrix from ZHETRD.
617  indre = indrd + n
618 * RWORK(INDRDD:INDRDD+N-1) is a copy of the diagonal entries over
619 * -written by ZSTEMR (the DSTERF path copies the diagonal to W).
620  indrdd = indre + n
621 * RWORK(INDREE:INDREE+N-1) is a copy of the off-diagonal entries over
622 * -written while computing the eigenvalues in DSTERF and ZSTEMR.
623  indree = indrdd + n
624 * INDRWK is the starting offset of the left-over real workspace, and
625 * LLRWORK is the remaining workspace size.
626  indrwk = indree + n
627  llrwork = lrwork - indrwk + 1
628
629 * IWORK(INDIBL:INDIBL+M-1) corresponds to IBLOCK in DSTEBZ and
630 * stores the block indices of each of the M<=N eigenvalues.
631  indibl = 1
632 * IWORK(INDISP:INDISP+NSPLIT-1) corresponds to ISPLIT in DSTEBZ and
633 * stores the starting and finishing indices of each block.
634  indisp = indibl + n
635 * IWORK(INDIFL:INDIFL+N-1) stores the indices of eigenvectors
636 * that corresponding to eigenvectors that fail to converge in
637 * ZSTEIN. This information is discarded; if any fail, the driver
638 * returns INFO > 0.
639  indifl = indisp + n
640 * INDIWO is the offset of the remaining integer workspace.
641  indiwo = indifl + n
642
643 *
644 * Call ZHETRD_2STAGE to reduce Hermitian matrix to tridiagonal form.
645 *
646  CALL zhetrd_2stage( jobz, uplo, n, a, lda, rwork( indrd ),
647  \$ rwork( indre ), work( indtau ),
648  \$ work( indhous ), lhtrd,
649  \$ work( indwk ), llwork, iinfo )
650 *
651 * If all eigenvalues are desired
652 * then call DSTERF or ZSTEMR and ZUNMTR.
653 *
654  test = .false.
655  IF( indeig ) THEN
656  IF( il.EQ.1 .AND. iu.EQ.n ) THEN
657  test = .true.
658  END IF
659  END IF
660  IF( ( alleig.OR.test ) .AND. ( ieeeok.EQ.1 ) ) THEN
661  IF( .NOT.wantz ) THEN
662  CALL dcopy( n, rwork( indrd ), 1, w, 1 )
663  CALL dcopy( n-1, rwork( indre ), 1, rwork( indree ), 1 )
664  CALL dsterf( n, w, rwork( indree ), info )
665  ELSE
666  CALL dcopy( n-1, rwork( indre ), 1, rwork( indree ), 1 )
667  CALL dcopy( n, rwork( indrd ), 1, rwork( indrdd ), 1 )
668 *
669  IF (abstol .LE. two*n*eps) THEN
670  tryrac = .true.
671  ELSE
672  tryrac = .false.
673  END IF
674  CALL zstemr( jobz, 'A', n, rwork( indrdd ),
675  \$ rwork( indree ), vl, vu, il, iu, m, w,
676  \$ z, ldz, n, isuppz, tryrac,
677  \$ rwork( indrwk ), llrwork,
678  \$ iwork, liwork, info )
679 *
680 * Apply unitary matrix used in reduction to tridiagonal
681 * form to eigenvectors returned by ZSTEMR.
682 *
683  IF( wantz .AND. info.EQ.0 ) THEN
684  indwkn = indwk
685  llwrkn = lwork - indwkn + 1
686  CALL zunmtr( 'L', uplo, 'N', n, m, a, lda,
687  \$ work( indtau ), z, ldz, work( indwkn ),
688  \$ llwrkn, iinfo )
689  END IF
690  END IF
691 *
692 *
693  IF( info.EQ.0 ) THEN
694  m = n
695  GO TO 30
696  END IF
697  info = 0
698  END IF
699 *
700 * Otherwise, call DSTEBZ and, if eigenvectors are desired, ZSTEIN.
701 * Also call DSTEBZ and ZSTEIN if ZSTEMR fails.
702 *
703  IF( wantz ) THEN
704  order = 'B'
705  ELSE
706  order = 'E'
707  END IF
708
709  CALL dstebz( range, order, n, vll, vuu, il, iu, abstll,
710  \$ rwork( indrd ), rwork( indre ), m, nsplit, w,
711  \$ iwork( indibl ), iwork( indisp ), rwork( indrwk ),
712  \$ iwork( indiwo ), info )
713 *
714  IF( wantz ) THEN
715  CALL zstein( n, rwork( indrd ), rwork( indre ), m, w,
716  \$ iwork( indibl ), iwork( indisp ), z, ldz,
717  \$ rwork( indrwk ), iwork( indiwo ), iwork( indifl ),
718  \$ info )
719 *
720 * Apply unitary matrix used in reduction to tridiagonal
721 * form to eigenvectors returned by ZSTEIN.
722 *
723  indwkn = indwk
724  llwrkn = lwork - indwkn + 1
725  CALL zunmtr( 'L', uplo, 'N', n, m, a, lda, work( indtau ), z,
726  \$ ldz, work( indwkn ), llwrkn, iinfo )
727  END IF
728 *
729 * If matrix was scaled, then rescale eigenvalues appropriately.
730 *
731  30 CONTINUE
732  IF( iscale.EQ.1 ) THEN
733  IF( info.EQ.0 ) THEN
734  imax = m
735  ELSE
736  imax = info - 1
737  END IF
738  CALL dscal( imax, one / sigma, w, 1 )
739  END IF
740 *
741 * If eigenvalues are not in order, then sort them, along with
742 * eigenvectors.
743 *
744  IF( wantz ) THEN
745  DO 50 j = 1, m - 1
746  i = 0
747  tmp1 = w( j )
748  DO 40 jj = j + 1, m
749  IF( w( jj ).LT.tmp1 ) THEN
750  i = jj
751  tmp1 = w( jj )
752  END IF
753  40 CONTINUE
754 *
755  IF( i.NE.0 ) THEN
756  itmp1 = iwork( indibl+i-1 )
757  w( i ) = w( j )
758  iwork( indibl+i-1 ) = iwork( indibl+j-1 )
759  w( j ) = tmp1
760  iwork( indibl+j-1 ) = itmp1
761  CALL zswap( n, z( 1, i ), 1, z( 1, j ), 1 )
762  END IF
763  50 CONTINUE
764  END IF
765 *
766 * Set WORK(1) to optimal workspace size.
767 *
768  work( 1 ) = lwmin
769  rwork( 1 ) = lrwmin
770  iwork( 1 ) = liwmin
771 *
772  RETURN
773 *
774 * End of ZHEEVR_2STAGE
775 *
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:69
integer function ilaenv2stage(ISPEC, NAME, OPTS, N1, N2, N3, N4)
ILAENV2STAGE
Definition: ilaenv2stage.f:149
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
ILAENV
Definition: ilaenv.f:162
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
subroutine zhetrd_2stage(VECT, UPLO, N, A, LDA, D, E, TAU, HOUS2, LHOUS2, WORK, LWORK, INFO)
ZHETRD_2STAGE
subroutine zstemr(JOBZ, RANGE, N, D, E, VL, VU, IL, IU, M, W, Z, LDZ, NZC, ISUPPZ, TRYRAC, WORK, LWORK, IWORK, LIWORK, INFO)
ZSTEMR
Definition: zstemr.f:338
subroutine zunmtr(SIDE, UPLO, TRANS, M, N, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
ZUNMTR
Definition: zunmtr.f:171
subroutine zstein(N, D, E, M, W, IBLOCK, ISPLIT, Z, LDZ, WORK, IWORK, IFAIL, INFO)
ZSTEIN
Definition: zstein.f:182
double precision function zlansy(NORM, UPLO, N, A, LDA, WORK)
ZLANSY returns the value of the 1-norm, or the Frobenius norm, or the infinity norm,...
Definition: zlansy.f:123
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: