LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

◆ cheevr_2stage()

subroutine cheevr_2stage ( character  JOBZ,
character  RANGE,
character  UPLO,
integer  N,
complex, dimension( lda, * )  A,
integer  LDA,
real  VL,
real  VU,
integer  IL,
integer  IU,
real  ABSTOL,
integer  M,
real, dimension( * )  W,
complex, dimension( ldz, * )  Z,
integer  LDZ,
integer, dimension( * )  ISUPPZ,
complex, dimension( * )  WORK,
integer  LWORK,
real, dimension( * )  RWORK,
integer  LRWORK,
integer, dimension( * )  IWORK,
integer  LIWORK,
integer  INFO 
)

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

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

Purpose:
 CHEEVR_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.

 CHEEVR_2STAGE first reduces the matrix A to tridiagonal form T with a call
 to CHETRD.  Then, whenever possible, CHEEVR_2STAGE calls CSTEMR to compute
 eigenspectrum using Relatively Robust Representations.  CSTEMR
 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 CSTEMR'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 : CHEEVR_2STAGE calls CSTEMR when the full spectrum is requested
 on machines which conform to the ieee-754 floating point standard.
 CHEEVR_2STAGE calls SSTEBZ and CSTEIN on non-ieee machines and
 when partial spectrum requests are made.

 Normal execution of CSTEMR 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, SSTEBZ and
          CSTEIN 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 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 REAL
          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 REAL
          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 REAL
          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
          SLAMCH( '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 REAL array, dimension (N)
          The first M elements contain the selected eigenvalues in
          ascending order.
[out]Z
          Z is COMPLEX 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 CSTEMR (tridiagonal
          matrix). The support of the eigenvectors of A is typically 
          1:N because of the unitary transformations applied by CUNMTR.
          Implemented only for RANGE = 'A' or 'I' and IU - IL = N - 1
[out]WORK
          WORK is COMPLEX 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 REAL 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
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
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).
  Denver, Colorado, USA, 2013.
  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 cheevr_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  REAL ABSTOL, VL, VU
418 * ..
419 * .. Array Arguments ..
420  INTEGER ISUPPZ( * ), IWORK( * )
421  REAL RWORK( * ), W( * )
422  COMPLEX A( LDA, * ), WORK( * ), Z( LDZ, * )
423 * ..
424 *
425 * =====================================================================
426 *
427 * .. Parameters ..
428  REAL ZERO, ONE, TWO
429  parameter( zero = 0.0e+0, one = 1.0e+0, two = 2.0e+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  REAL 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  REAL SLAMCH, CLANSY
447  EXTERNAL lsame, slamch, clansy, ilaenv, ilaenv2stage
448 * ..
449 * .. External Subroutines ..
450  EXTERNAL scopy, sscal, sstebz, ssterf, xerbla, csscal,
452 * ..
453 * .. Intrinsic Functions ..
454  INTRINSIC real, max, min, sqrt
455 * ..
456 * .. Executable Statements ..
457 *
458 * Test the input parameters.
459 *
460  ieeeok = ilaenv( 10, 'CHEEVR', '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, 'CHETRD_2STAGE', jobz, n, -1, -1, -1 )
472  ib = ilaenv2stage( 2, 'CHETRD_2STAGE', jobz, n, kd, -1, -1 )
473  lhtrd = ilaenv2stage( 3, 'CHETRD_2STAGE', jobz, n, kd, ib, -1 )
474  lwtrd = ilaenv2stage( 4, 'CHETRD_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( 'CHEEVR_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 ) = real( a( 1, 1 ) )
542  ELSE
543  IF( vl.LT.real( a( 1, 1 ) ) .AND. vu.GE.real( a( 1, 1 ) ) )
544  $ THEN
545  m = 1
546  w( 1 ) = real( 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 = slamch( 'Safe minimum' )
560  eps = slamch( '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 = clansy( '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 csscal( n-j+1, sigma, a( j, j ), 1 )
586  10 CONTINUE
587  ELSE
588  DO 20 j = 1, n
589  CALL csscal( 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 SSTERF or CSTEMR fail.
602 
603 * WORK(INDTAU:INDTAU+N-1) stores the complex scalar factors of the
604 * elementary reflectors used in CHETRD.
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 CHETRD.
617  indre = indrd + n
618 * RWORK(INDRDD:INDRDD+N-1) is a copy of the diagonal entries over
619 * -written by CSTEMR (the SSTERF 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 SSTERF and CSTEMR.
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 SSTEBZ 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 SSTEBZ 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 * CSTEIN. 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 CHETRD_2STAGE to reduce Hermitian matrix to tridiagonal form.
645 *
646  CALL chetrd_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 SSTERF or CSTEMR and CUNMTR.
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 scopy( n, rwork( indrd ), 1, w, 1 )
663  CALL scopy( n-1, rwork( indre ), 1, rwork( indree ), 1 )
664  CALL ssterf( n, w, rwork( indree ), info )
665  ELSE
666  CALL scopy( n-1, rwork( indre ), 1, rwork( indree ), 1 )
667  CALL scopy( 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 cstemr( 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 CSTEMR.
682 *
683  IF( wantz .AND. info.EQ.0 ) THEN
684  indwkn = indwk
685  llwrkn = lwork - indwkn + 1
686  CALL cunmtr( '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 SSTEBZ and, if eigenvectors are desired, CSTEIN.
701 * Also call SSTEBZ and CSTEIN if CSTEMR fails.
702 *
703  IF( wantz ) THEN
704  order = 'B'
705  ELSE
706  order = 'E'
707  END IF
708 
709  CALL sstebz( 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 cstein( 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 CSTEIN.
722 *
723  indwkn = indwk
724  llwrkn = lwork - indwkn + 1
725  CALL cunmtr( '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 sscal( 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 cswap( 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 CHEEVR_2STAGE
775 *
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 ssterf(N, D, E, INFO)
SSTERF
Definition: ssterf.f:86
subroutine sstebz(RANGE, ORDER, N, VL, VU, IL, IU, ABSTOL, D, E, M, NSPLIT, W, IBLOCK, ISPLIT, WORK, IWORK, INFO)
SSTEBZ
Definition: sstebz.f:273
subroutine csscal(N, SA, CX, INCX)
CSSCAL
Definition: csscal.f:78
subroutine cswap(N, CX, INCX, CY, INCY)
CSWAP
Definition: cswap.f:81
subroutine chetrd_2stage(VECT, UPLO, N, A, LDA, D, E, TAU, HOUS2, LHOUS2, WORK, LWORK, INFO)
CHETRD_2STAGE
subroutine cunmtr(SIDE, UPLO, TRANS, M, N, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
CUNMTR
Definition: cunmtr.f:172
subroutine cstein(N, D, E, M, W, IBLOCK, ISPLIT, Z, LDZ, WORK, IWORK, IFAIL, INFO)
CSTEIN
Definition: cstein.f:182
subroutine cstemr(JOBZ, RANGE, N, D, E, VL, VU, IL, IU, M, W, Z, LDZ, NZC, ISUPPZ, TRYRAC, WORK, LWORK, IWORK, LIWORK, INFO)
CSTEMR
Definition: cstemr.f:338
real function clansy(NORM, UPLO, N, A, LDA, WORK)
CLANSY returns the value of the 1-norm, or the Frobenius norm, or the infinity norm,...
Definition: clansy.f:123
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY
Definition: scopy.f:82
subroutine sscal(N, SA, SX, INCX)
SSCAL
Definition: sscal.f:79
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:68
Here is the call graph for this function:
Here is the caller graph for this function: