Annotation of rpl/lapack/lapack/zheevr.f, revision 1.22

1.10      bertrand    1: *> \brief <b> ZHEEVR computes the eigenvalues and, optionally, the left and/or right eigenvectors for HE matrices</b>
                      2: *
                      3: *  =========== DOCUMENTATION ===========
                      4: *
1.18      bertrand    5: * Online html documentation available at
                      6: *            http://www.netlib.org/lapack/explore-html/
1.10      bertrand    7: *
                      8: *> \htmlonly
1.18      bertrand    9: *> Download ZHEEVR + dependencies
                     10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zheevr.f">
                     11: *> [TGZ]</a>
                     12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zheevr.f">
                     13: *> [ZIP]</a>
                     14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zheevr.f">
1.10      bertrand   15: *> [TXT]</a>
1.18      bertrand   16: *> \endhtmlonly
1.10      bertrand   17: *
                     18: *  Definition:
                     19: *  ===========
                     20: *
                     21: *       SUBROUTINE ZHEEVR( JOBZ, RANGE, UPLO, N, A, LDA, VL, VU, IL, IU,
                     22: *                          ABSTOL, M, W, Z, LDZ, ISUPPZ, WORK, LWORK,
                     23: *                          RWORK, LRWORK, IWORK, LIWORK, INFO )
1.18      bertrand   24: *
1.10      bertrand   25: *       .. Scalar Arguments ..
                     26: *       CHARACTER          JOBZ, RANGE, UPLO
                     27: *       INTEGER            IL, INFO, IU, LDA, LDZ, LIWORK, LRWORK, LWORK,
                     28: *      $                   M, N
                     29: *       DOUBLE PRECISION   ABSTOL, VL, VU
                     30: *       ..
                     31: *       .. Array Arguments ..
                     32: *       INTEGER            ISUPPZ( * ), IWORK( * )
                     33: *       DOUBLE PRECISION   RWORK( * ), W( * )
                     34: *       COMPLEX*16         A( LDA, * ), WORK( * ), Z( LDZ, * )
                     35: *       ..
1.18      bertrand   36: *
1.10      bertrand   37: *
                     38: *> \par Purpose:
                     39: *  =============
                     40: *>
                     41: *> \verbatim
                     42: *>
                     43: *> ZHEEVR computes selected eigenvalues and, optionally, eigenvectors
                     44: *> of a complex Hermitian matrix A.  Eigenvalues and eigenvectors can
                     45: *> be selected by specifying either a range of values or a range of
                     46: *> indices for the desired eigenvalues.
                     47: *>
                     48: *> ZHEEVR first reduces the matrix A to tridiagonal form T with a call
                     49: *> to ZHETRD.  Then, whenever possible, ZHEEVR calls ZSTEMR to compute
                     50: *> eigenspectrum using Relatively Robust Representations.  ZSTEMR
                     51: *> computes eigenvalues by the dqds algorithm, while orthogonal
                     52: *> eigenvectors are computed from various "good" L D L^T representations
                     53: *> (also known as Relatively Robust Representations). Gram-Schmidt
                     54: *> orthogonalization is avoided as far as possible. More specifically,
                     55: *> the various steps of the algorithm are as follows.
                     56: *>
                     57: *> For each unreduced block (submatrix) of T,
                     58: *>    (a) Compute T - sigma I  = L D L^T, so that L and D
                     59: *>        define all the wanted eigenvalues to high relative accuracy.
                     60: *>        This means that small relative changes in the entries of D and L
                     61: *>        cause only small relative changes in the eigenvalues and
                     62: *>        eigenvectors. The standard (unfactored) representation of the
                     63: *>        tridiagonal matrix T does not have this property in general.
                     64: *>    (b) Compute the eigenvalues to suitable accuracy.
                     65: *>        If the eigenvectors are desired, the algorithm attains full
                     66: *>        accuracy of the computed eigenvalues only right before
                     67: *>        the corresponding vectors have to be computed, see steps c) and d).
                     68: *>    (c) For each cluster of close eigenvalues, select a new
                     69: *>        shift close to the cluster, find a new factorization, and refine
                     70: *>        the shifted eigenvalues to suitable accuracy.
                     71: *>    (d) For each eigenvalue with a large enough relative separation compute
                     72: *>        the corresponding eigenvector by forming a rank revealing twisted
                     73: *>        factorization. Go back to (c) for any clusters that remain.
                     74: *>
                     75: *> The desired accuracy of the output can be specified by the input
                     76: *> parameter ABSTOL.
                     77: *>
1.22    ! bertrand   78: *> For more details, see ZSTEMR's documentation and:
1.10      bertrand   79: *> - Inderjit S. Dhillon and Beresford N. Parlett: "Multiple representations
                     80: *>   to compute orthogonal eigenvectors of symmetric tridiagonal matrices,"
                     81: *>   Linear Algebra and its Applications, 387(1), pp. 1-28, August 2004.
                     82: *> - Inderjit Dhillon and Beresford Parlett: "Orthogonal Eigenvectors and
                     83: *>   Relative Gaps," SIAM Journal on Matrix Analysis and Applications, Vol. 25,
                     84: *>   2004.  Also LAPACK Working Note 154.
                     85: *> - Inderjit Dhillon: "A new O(n^2) algorithm for the symmetric
                     86: *>   tridiagonal eigenvalue/eigenvector problem",
                     87: *>   Computer Science Division Technical Report No. UCB/CSD-97-971,
                     88: *>   UC Berkeley, May 1997.
                     89: *>
                     90: *>
                     91: *> Note 1 : ZHEEVR calls ZSTEMR when the full spectrum is requested
                     92: *> on machines which conform to the ieee-754 floating point standard.
                     93: *> ZHEEVR calls DSTEBZ and ZSTEIN on non-ieee machines and
                     94: *> when partial spectrum requests are made.
                     95: *>
                     96: *> Normal execution of ZSTEMR may create NaNs and infinities and
                     97: *> hence may abort due to a floating point exception in environments
                     98: *> which do not handle NaNs and infinities in the ieee standard default
                     99: *> manner.
                    100: *> \endverbatim
                    101: *
                    102: *  Arguments:
                    103: *  ==========
                    104: *
                    105: *> \param[in] JOBZ
                    106: *> \verbatim
                    107: *>          JOBZ is CHARACTER*1
                    108: *>          = 'N':  Compute eigenvalues only;
                    109: *>          = 'V':  Compute eigenvalues and eigenvectors.
                    110: *> \endverbatim
                    111: *>
                    112: *> \param[in] RANGE
                    113: *> \verbatim
                    114: *>          RANGE is CHARACTER*1
                    115: *>          = 'A': all eigenvalues will be found.
                    116: *>          = 'V': all eigenvalues in the half-open interval (VL,VU]
                    117: *>                 will be found.
                    118: *>          = 'I': the IL-th through IU-th eigenvalues will be found.
                    119: *>          For RANGE = 'V' or 'I' and IU - IL < N - 1, DSTEBZ and
                    120: *>          ZSTEIN are called
                    121: *> \endverbatim
                    122: *>
                    123: *> \param[in] UPLO
                    124: *> \verbatim
                    125: *>          UPLO is CHARACTER*1
                    126: *>          = 'U':  Upper triangle of A is stored;
                    127: *>          = 'L':  Lower triangle of A is stored.
                    128: *> \endverbatim
                    129: *>
                    130: *> \param[in] N
                    131: *> \verbatim
                    132: *>          N is INTEGER
                    133: *>          The order of the matrix A.  N >= 0.
                    134: *> \endverbatim
                    135: *>
                    136: *> \param[in,out] A
                    137: *> \verbatim
                    138: *>          A is COMPLEX*16 array, dimension (LDA, N)
                    139: *>          On entry, the Hermitian matrix A.  If UPLO = 'U', the
                    140: *>          leading N-by-N upper triangular part of A contains the
                    141: *>          upper triangular part of the matrix A.  If UPLO = 'L',
                    142: *>          the leading N-by-N lower triangular part of A contains
                    143: *>          the lower triangular part of the matrix A.
                    144: *>          On exit, the lower triangle (if UPLO='L') or the upper
                    145: *>          triangle (if UPLO='U') of A, including the diagonal, is
                    146: *>          destroyed.
                    147: *> \endverbatim
                    148: *>
                    149: *> \param[in] LDA
                    150: *> \verbatim
                    151: *>          LDA is INTEGER
                    152: *>          The leading dimension of the array A.  LDA >= max(1,N).
                    153: *> \endverbatim
                    154: *>
                    155: *> \param[in] VL
                    156: *> \verbatim
                    157: *>          VL is DOUBLE PRECISION
1.16      bertrand  158: *>          If RANGE='V', the lower bound of the interval to
                    159: *>          be searched for eigenvalues. VL < VU.
                    160: *>          Not referenced if RANGE = 'A' or 'I'.
1.10      bertrand  161: *> \endverbatim
                    162: *>
                    163: *> \param[in] VU
                    164: *> \verbatim
                    165: *>          VU is DOUBLE PRECISION
1.16      bertrand  166: *>          If RANGE='V', the upper bound of the interval to
1.10      bertrand  167: *>          be searched for eigenvalues. VL < VU.
                    168: *>          Not referenced if RANGE = 'A' or 'I'.
                    169: *> \endverbatim
                    170: *>
                    171: *> \param[in] IL
                    172: *> \verbatim
                    173: *>          IL is INTEGER
1.16      bertrand  174: *>          If RANGE='I', the index of the
                    175: *>          smallest eigenvalue to be returned.
                    176: *>          1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
                    177: *>          Not referenced if RANGE = 'A' or 'V'.
1.10      bertrand  178: *> \endverbatim
                    179: *>
                    180: *> \param[in] IU
                    181: *> \verbatim
                    182: *>          IU is INTEGER
1.16      bertrand  183: *>          If RANGE='I', the index of the
                    184: *>          largest eigenvalue to be returned.
1.10      bertrand  185: *>          1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
                    186: *>          Not referenced if RANGE = 'A' or 'V'.
                    187: *> \endverbatim
                    188: *>
                    189: *> \param[in] ABSTOL
                    190: *> \verbatim
                    191: *>          ABSTOL is DOUBLE PRECISION
                    192: *>          The absolute error tolerance for the eigenvalues.
                    193: *>          An approximate eigenvalue is accepted as converged
                    194: *>          when it is determined to lie in an interval [a,b]
                    195: *>          of width less than or equal to
                    196: *>
                    197: *>                  ABSTOL + EPS *   max( |a|,|b| ) ,
                    198: *>
                    199: *>          where EPS is the machine precision.  If ABSTOL is less than
                    200: *>          or equal to zero, then  EPS*|T|  will be used in its place,
                    201: *>          where |T| is the 1-norm of the tridiagonal matrix obtained
                    202: *>          by reducing A to tridiagonal form.
                    203: *>
                    204: *>          See "Computing Small Singular Values of Bidiagonal Matrices
                    205: *>          with Guaranteed High Relative Accuracy," by Demmel and
                    206: *>          Kahan, LAPACK Working Note #3.
                    207: *>
                    208: *>          If high relative accuracy is important, set ABSTOL to
                    209: *>          DLAMCH( 'Safe minimum' ).  Doing so will guarantee that
                    210: *>          eigenvalues are computed to high relative accuracy when
                    211: *>          possible in future releases.  The current code does not
                    212: *>          make any guarantees about high relative accuracy, but
1.21      bertrand  213: *>          future releases will. See J. Barlow and J. Demmel,
1.10      bertrand  214: *>          "Computing Accurate Eigensystems of Scaled Diagonally
                    215: *>          Dominant Matrices", LAPACK Working Note #7, for a discussion
                    216: *>          of which matrices define their eigenvalues to high relative
                    217: *>          accuracy.
                    218: *> \endverbatim
                    219: *>
                    220: *> \param[out] M
                    221: *> \verbatim
                    222: *>          M is INTEGER
                    223: *>          The total number of eigenvalues found.  0 <= M <= N.
                    224: *>          If RANGE = 'A', M = N, and if RANGE = 'I', M = IU-IL+1.
                    225: *> \endverbatim
                    226: *>
                    227: *> \param[out] W
                    228: *> \verbatim
                    229: *>          W is DOUBLE PRECISION array, dimension (N)
                    230: *>          The first M elements contain the selected eigenvalues in
                    231: *>          ascending order.
                    232: *> \endverbatim
                    233: *>
                    234: *> \param[out] Z
                    235: *> \verbatim
                    236: *>          Z is COMPLEX*16 array, dimension (LDZ, max(1,M))
                    237: *>          If JOBZ = 'V', then if INFO = 0, the first M columns of Z
                    238: *>          contain the orthonormal eigenvectors of the matrix A
                    239: *>          corresponding to the selected eigenvalues, with the i-th
                    240: *>          column of Z holding the eigenvector associated with W(i).
                    241: *>          If JOBZ = 'N', then Z is not referenced.
                    242: *>          Note: the user must ensure that at least max(1,M) columns are
                    243: *>          supplied in the array Z; if RANGE = 'V', the exact value of M
                    244: *>          is not known in advance and an upper bound must be used.
                    245: *> \endverbatim
                    246: *>
                    247: *> \param[in] LDZ
                    248: *> \verbatim
                    249: *>          LDZ is INTEGER
                    250: *>          The leading dimension of the array Z.  LDZ >= 1, and if
                    251: *>          JOBZ = 'V', LDZ >= max(1,N).
                    252: *> \endverbatim
                    253: *>
                    254: *> \param[out] ISUPPZ
                    255: *> \verbatim
                    256: *>          ISUPPZ is INTEGER array, dimension ( 2*max(1,M) )
                    257: *>          The support of the eigenvectors in Z, i.e., the indices
                    258: *>          indicating the nonzero elements in Z. The i-th eigenvector
                    259: *>          is nonzero only in elements ISUPPZ( 2*i-1 ) through
1.18      bertrand  260: *>          ISUPPZ( 2*i ). This is an output of ZSTEMR (tridiagonal
                    261: *>          matrix). The support of the eigenvectors of A is typically
                    262: *>          1:N because of the unitary transformations applied by ZUNMTR.
1.10      bertrand  263: *>          Implemented only for RANGE = 'A' or 'I' and IU - IL = N - 1
                    264: *> \endverbatim
                    265: *>
                    266: *> \param[out] WORK
                    267: *> \verbatim
                    268: *>          WORK is COMPLEX*16 array, dimension (MAX(1,LWORK))
                    269: *>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
                    270: *> \endverbatim
                    271: *>
                    272: *> \param[in] LWORK
                    273: *> \verbatim
                    274: *>          LWORK is INTEGER
                    275: *>          The length of the array WORK.  LWORK >= max(1,2*N).
                    276: *>          For optimal efficiency, LWORK >= (NB+1)*N,
                    277: *>          where NB is the max of the blocksize for ZHETRD and for
                    278: *>          ZUNMTR as returned by ILAENV.
                    279: *>
                    280: *>          If LWORK = -1, then a workspace query is assumed; the routine
                    281: *>          only calculates the optimal sizes of the WORK, RWORK and
                    282: *>          IWORK arrays, returns these values as the first entries of
                    283: *>          the WORK, RWORK and IWORK arrays, and no error message
                    284: *>          related to LWORK or LRWORK or LIWORK is issued by XERBLA.
                    285: *> \endverbatim
                    286: *>
                    287: *> \param[out] RWORK
                    288: *> \verbatim
                    289: *>          RWORK is DOUBLE PRECISION array, dimension (MAX(1,LRWORK))
                    290: *>          On exit, if INFO = 0, RWORK(1) returns the optimal
                    291: *>          (and minimal) LRWORK.
                    292: *> \endverbatim
                    293: *>
                    294: *> \param[in] LRWORK
                    295: *> \verbatim
                    296: *>          LRWORK is INTEGER
                    297: *>          The length of the array RWORK.  LRWORK >= max(1,24*N).
                    298: *>
                    299: *>          If LRWORK = -1, then a workspace query is assumed; the
                    300: *>          routine only calculates the optimal sizes of the WORK, RWORK
                    301: *>          and IWORK arrays, returns these values as the first entries
                    302: *>          of the WORK, RWORK and IWORK arrays, and no error message
                    303: *>          related to LWORK or LRWORK or LIWORK is issued by XERBLA.
                    304: *> \endverbatim
                    305: *>
                    306: *> \param[out] IWORK
                    307: *> \verbatim
                    308: *>          IWORK is INTEGER array, dimension (MAX(1,LIWORK))
                    309: *>          On exit, if INFO = 0, IWORK(1) returns the optimal
                    310: *>          (and minimal) LIWORK.
                    311: *> \endverbatim
                    312: *>
                    313: *> \param[in] LIWORK
                    314: *> \verbatim
                    315: *>          LIWORK is INTEGER
                    316: *>          The dimension of the array IWORK.  LIWORK >= max(1,10*N).
                    317: *>
                    318: *>          If LIWORK = -1, then a workspace query is assumed; the
                    319: *>          routine only calculates the optimal sizes of the WORK, RWORK
                    320: *>          and IWORK arrays, returns these values as the first entries
                    321: *>          of the WORK, RWORK and IWORK arrays, and no error message
                    322: *>          related to LWORK or LRWORK or LIWORK is issued by XERBLA.
                    323: *> \endverbatim
                    324: *>
                    325: *> \param[out] INFO
                    326: *> \verbatim
                    327: *>          INFO is INTEGER
                    328: *>          = 0:  successful exit
                    329: *>          < 0:  if INFO = -i, the i-th argument had an illegal value
                    330: *>          > 0:  Internal error
                    331: *> \endverbatim
                    332: *
                    333: *  Authors:
                    334: *  ========
                    335: *
1.18      bertrand  336: *> \author Univ. of Tennessee
                    337: *> \author Univ. of California Berkeley
                    338: *> \author Univ. of Colorado Denver
                    339: *> \author NAG Ltd.
1.10      bertrand  340: *
                    341: *> \ingroup complex16HEeigen
                    342: *
                    343: *> \par Contributors:
                    344: *  ==================
                    345: *>
                    346: *>     Inderjit Dhillon, IBM Almaden, USA \n
                    347: *>     Osni Marques, LBNL/NERSC, USA \n
                    348: *>     Ken Stanley, Computer Science Division, University of
                    349: *>       California at Berkeley, USA \n
                    350: *>     Jason Riedy, Computer Science Division, University of
                    351: *>       California at Berkeley, USA \n
                    352: *>
                    353: *  =====================================================================
1.1       bertrand  354:       SUBROUTINE ZHEEVR( JOBZ, RANGE, UPLO, N, A, LDA, VL, VU, IL, IU,
                    355:      $                   ABSTOL, M, W, Z, LDZ, ISUPPZ, WORK, LWORK,
                    356:      $                   RWORK, LRWORK, IWORK, LIWORK, INFO )
                    357: *
1.22    ! bertrand  358: *  -- LAPACK driver routine --
1.1       bertrand  359: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
                    360: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
                    361: *
                    362: *     .. Scalar Arguments ..
                    363:       CHARACTER          JOBZ, RANGE, UPLO
                    364:       INTEGER            IL, INFO, IU, LDA, LDZ, LIWORK, LRWORK, LWORK,
                    365:      $                   M, N
                    366:       DOUBLE PRECISION   ABSTOL, VL, VU
                    367: *     ..
                    368: *     .. Array Arguments ..
                    369:       INTEGER            ISUPPZ( * ), IWORK( * )
                    370:       DOUBLE PRECISION   RWORK( * ), W( * )
                    371:       COMPLEX*16         A( LDA, * ), WORK( * ), Z( LDZ, * )
                    372: *     ..
                    373: *
1.9       bertrand  374: *  =====================================================================
1.1       bertrand  375: *
                    376: *     .. Parameters ..
                    377:       DOUBLE PRECISION   ZERO, ONE, TWO
                    378:       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0, TWO = 2.0D+0 )
                    379: *     ..
                    380: *     .. Local Scalars ..
                    381:       LOGICAL            ALLEIG, INDEIG, LOWER, LQUERY, TEST, VALEIG,
                    382:      $                   WANTZ, TRYRAC
                    383:       CHARACTER          ORDER
                    384:       INTEGER            I, IEEEOK, IINFO, IMAX, INDIBL, INDIFL, INDISP,
                    385:      $                   INDIWO, INDRD, INDRDD, INDRE, INDREE, INDRWK,
                    386:      $                   INDTAU, INDWK, INDWKN, ISCALE, ITMP1, J, JJ,
                    387:      $                   LIWMIN, LLWORK, LLRWORK, LLWRKN, LRWMIN,
                    388:      $                   LWKOPT, LWMIN, NB, NSPLIT
                    389:       DOUBLE PRECISION   ABSTLL, ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN,
                    390:      $                   SIGMA, SMLNUM, TMP1, VLL, VUU
                    391: *     ..
                    392: *     .. External Functions ..
                    393:       LOGICAL            LSAME
                    394:       INTEGER            ILAENV
                    395:       DOUBLE PRECISION   DLAMCH, ZLANSY
                    396:       EXTERNAL           LSAME, ILAENV, DLAMCH, ZLANSY
                    397: *     ..
                    398: *     .. External Subroutines ..
                    399:       EXTERNAL           DCOPY, DSCAL, DSTEBZ, DSTERF, XERBLA, ZDSCAL,
                    400:      $                   ZHETRD, ZSTEMR, ZSTEIN, ZSWAP, ZUNMTR
                    401: *     ..
                    402: *     .. Intrinsic Functions ..
                    403:       INTRINSIC          DBLE, MAX, MIN, SQRT
                    404: *     ..
                    405: *     .. Executable Statements ..
                    406: *
                    407: *     Test the input parameters.
                    408: *
                    409:       IEEEOK = ILAENV( 10, 'ZHEEVR', 'N', 1, 2, 3, 4 )
                    410: *
                    411:       LOWER = LSAME( UPLO, 'L' )
                    412:       WANTZ = LSAME( JOBZ, 'V' )
                    413:       ALLEIG = LSAME( RANGE, 'A' )
                    414:       VALEIG = LSAME( RANGE, 'V' )
                    415:       INDEIG = LSAME( RANGE, 'I' )
                    416: *
                    417:       LQUERY = ( ( LWORK.EQ.-1 ) .OR. ( LRWORK.EQ.-1 ) .OR.
                    418:      $         ( LIWORK.EQ.-1 ) )
                    419: *
                    420:       LRWMIN = MAX( 1, 24*N )
                    421:       LIWMIN = MAX( 1, 10*N )
                    422:       LWMIN = MAX( 1, 2*N )
                    423: *
                    424:       INFO = 0
                    425:       IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN
                    426:          INFO = -1
                    427:       ELSE IF( .NOT.( ALLEIG .OR. VALEIG .OR. INDEIG ) ) THEN
                    428:          INFO = -2
                    429:       ELSE IF( .NOT.( LOWER .OR. LSAME( UPLO, 'U' ) ) ) THEN
                    430:          INFO = -3
                    431:       ELSE IF( N.LT.0 ) THEN
                    432:          INFO = -4
                    433:       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
                    434:          INFO = -6
                    435:       ELSE
                    436:          IF( VALEIG ) THEN
                    437:             IF( N.GT.0 .AND. VU.LE.VL )
                    438:      $         INFO = -8
                    439:          ELSE IF( INDEIG ) THEN
                    440:             IF( IL.LT.1 .OR. IL.GT.MAX( 1, N ) ) THEN
                    441:                INFO = -9
                    442:             ELSE IF( IU.LT.MIN( N, IL ) .OR. IU.GT.N ) THEN
                    443:                INFO = -10
                    444:             END IF
                    445:          END IF
                    446:       END IF
                    447:       IF( INFO.EQ.0 ) THEN
                    448:          IF( LDZ.LT.1 .OR. ( WANTZ .AND. LDZ.LT.N ) ) THEN
                    449:             INFO = -15
                    450:          END IF
                    451:       END IF
                    452: *
                    453:       IF( INFO.EQ.0 ) THEN
                    454:          NB = ILAENV( 1, 'ZHETRD', UPLO, N, -1, -1, -1 )
                    455:          NB = MAX( NB, ILAENV( 1, 'ZUNMTR', UPLO, N, -1, -1, -1 ) )
                    456:          LWKOPT = MAX( ( NB+1 )*N, LWMIN )
                    457:          WORK( 1 ) = LWKOPT
                    458:          RWORK( 1 ) = LRWMIN
                    459:          IWORK( 1 ) = LIWMIN
                    460: *
                    461:          IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN
                    462:             INFO = -18
                    463:          ELSE IF( LRWORK.LT.LRWMIN .AND. .NOT.LQUERY ) THEN
                    464:             INFO = -20
                    465:          ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT.LQUERY ) THEN
                    466:             INFO = -22
                    467:          END IF
                    468:       END IF
                    469: *
                    470:       IF( INFO.NE.0 ) THEN
                    471:          CALL XERBLA( 'ZHEEVR', -INFO )
                    472:          RETURN
                    473:       ELSE IF( LQUERY ) THEN
                    474:          RETURN
                    475:       END IF
                    476: *
                    477: *     Quick return if possible
                    478: *
                    479:       M = 0
                    480:       IF( N.EQ.0 ) THEN
                    481:          WORK( 1 ) = 1
                    482:          RETURN
                    483:       END IF
                    484: *
                    485:       IF( N.EQ.1 ) THEN
                    486:          WORK( 1 ) = 2
                    487:          IF( ALLEIG .OR. INDEIG ) THEN
                    488:             M = 1
                    489:             W( 1 ) = DBLE( A( 1, 1 ) )
                    490:          ELSE
                    491:             IF( VL.LT.DBLE( A( 1, 1 ) ) .AND. VU.GE.DBLE( A( 1, 1 ) ) )
                    492:      $           THEN
                    493:                M = 1
                    494:                W( 1 ) = DBLE( A( 1, 1 ) )
                    495:             END IF
                    496:          END IF
1.5       bertrand  497:          IF( WANTZ ) THEN
                    498:             Z( 1, 1 ) = ONE
                    499:             ISUPPZ( 1 ) = 1
                    500:             ISUPPZ( 2 ) = 1
                    501:          END IF
1.1       bertrand  502:          RETURN
                    503:       END IF
                    504: *
                    505: *     Get machine constants.
                    506: *
                    507:       SAFMIN = DLAMCH( 'Safe minimum' )
                    508:       EPS = DLAMCH( 'Precision' )
                    509:       SMLNUM = SAFMIN / EPS
                    510:       BIGNUM = ONE / SMLNUM
                    511:       RMIN = SQRT( SMLNUM )
                    512:       RMAX = MIN( SQRT( BIGNUM ), ONE / SQRT( SQRT( SAFMIN ) ) )
                    513: *
                    514: *     Scale matrix to allowable range, if necessary.
                    515: *
                    516:       ISCALE = 0
                    517:       ABSTLL = ABSTOL
                    518:       IF (VALEIG) THEN
                    519:          VLL = VL
                    520:          VUU = VU
                    521:       END IF
                    522:       ANRM = ZLANSY( 'M', UPLO, N, A, LDA, RWORK )
                    523:       IF( ANRM.GT.ZERO .AND. ANRM.LT.RMIN ) THEN
                    524:          ISCALE = 1
                    525:          SIGMA = RMIN / ANRM
                    526:       ELSE IF( ANRM.GT.RMAX ) THEN
                    527:          ISCALE = 1
                    528:          SIGMA = RMAX / ANRM
                    529:       END IF
                    530:       IF( ISCALE.EQ.1 ) THEN
                    531:          IF( LOWER ) THEN
                    532:             DO 10 J = 1, N
                    533:                CALL ZDSCAL( N-J+1, SIGMA, A( J, J ), 1 )
                    534:    10       CONTINUE
                    535:          ELSE
                    536:             DO 20 J = 1, N
                    537:                CALL ZDSCAL( J, SIGMA, A( 1, J ), 1 )
                    538:    20       CONTINUE
                    539:          END IF
                    540:          IF( ABSTOL.GT.0 )
                    541:      $      ABSTLL = ABSTOL*SIGMA
                    542:          IF( VALEIG ) THEN
                    543:             VLL = VL*SIGMA
                    544:             VUU = VU*SIGMA
                    545:          END IF
                    546:       END IF
                    547: 
                    548: *     Initialize indices into workspaces.  Note: The IWORK indices are
                    549: *     used only if DSTERF or ZSTEMR fail.
                    550: 
                    551: *     WORK(INDTAU:INDTAU+N-1) stores the complex scalar factors of the
                    552: *     elementary reflectors used in ZHETRD.
                    553:       INDTAU = 1
                    554: *     INDWK is the starting offset of the remaining complex workspace,
                    555: *     and LLWORK is the remaining complex workspace size.
                    556:       INDWK = INDTAU + N
                    557:       LLWORK = LWORK - INDWK + 1
                    558: 
                    559: *     RWORK(INDRD:INDRD+N-1) stores the real tridiagonal's diagonal
                    560: *     entries.
                    561:       INDRD = 1
                    562: *     RWORK(INDRE:INDRE+N-1) stores the off-diagonal entries of the
                    563: *     tridiagonal matrix from ZHETRD.
                    564:       INDRE = INDRD + N
                    565: *     RWORK(INDRDD:INDRDD+N-1) is a copy of the diagonal entries over
                    566: *     -written by ZSTEMR (the DSTERF path copies the diagonal to W).
                    567:       INDRDD = INDRE + N
                    568: *     RWORK(INDREE:INDREE+N-1) is a copy of the off-diagonal entries over
                    569: *     -written while computing the eigenvalues in DSTERF and ZSTEMR.
                    570:       INDREE = INDRDD + N
                    571: *     INDRWK is the starting offset of the left-over real workspace, and
                    572: *     LLRWORK is the remaining workspace size.
                    573:       INDRWK = INDREE + N
                    574:       LLRWORK = LRWORK - INDRWK + 1
                    575: 
                    576: *     IWORK(INDIBL:INDIBL+M-1) corresponds to IBLOCK in DSTEBZ and
                    577: *     stores the block indices of each of the M<=N eigenvalues.
                    578:       INDIBL = 1
                    579: *     IWORK(INDISP:INDISP+NSPLIT-1) corresponds to ISPLIT in DSTEBZ and
                    580: *     stores the starting and finishing indices of each block.
                    581:       INDISP = INDIBL + N
                    582: *     IWORK(INDIFL:INDIFL+N-1) stores the indices of eigenvectors
                    583: *     that corresponding to eigenvectors that fail to converge in
                    584: *     DSTEIN.  This information is discarded; if any fail, the driver
                    585: *     returns INFO > 0.
                    586:       INDIFL = INDISP + N
                    587: *     INDIWO is the offset of the remaining integer workspace.
1.13      bertrand  588:       INDIWO = INDIFL + N
1.1       bertrand  589: 
                    590: *
                    591: *     Call ZHETRD to reduce Hermitian matrix to tridiagonal form.
                    592: *
                    593:       CALL ZHETRD( UPLO, N, A, LDA, RWORK( INDRD ), RWORK( INDRE ),
                    594:      $             WORK( INDTAU ), WORK( INDWK ), LLWORK, IINFO )
                    595: *
                    596: *     If all eigenvalues are desired
                    597: *     then call DSTERF or ZSTEMR and ZUNMTR.
                    598: *
                    599:       TEST = .FALSE.
                    600:       IF( INDEIG ) THEN
                    601:          IF( IL.EQ.1 .AND. IU.EQ.N ) THEN
                    602:             TEST = .TRUE.
                    603:          END IF
                    604:       END IF
                    605:       IF( ( ALLEIG.OR.TEST ) .AND. ( IEEEOK.EQ.1 ) ) THEN
                    606:          IF( .NOT.WANTZ ) THEN
                    607:             CALL DCOPY( N, RWORK( INDRD ), 1, W, 1 )
                    608:             CALL DCOPY( N-1, RWORK( INDRE ), 1, RWORK( INDREE ), 1 )
                    609:             CALL DSTERF( N, W, RWORK( INDREE ), INFO )
                    610:          ELSE
                    611:             CALL DCOPY( N-1, RWORK( INDRE ), 1, RWORK( INDREE ), 1 )
                    612:             CALL DCOPY( N, RWORK( INDRD ), 1, RWORK( INDRDD ), 1 )
                    613: *
                    614:             IF (ABSTOL .LE. TWO*N*EPS) THEN
                    615:                TRYRAC = .TRUE.
                    616:             ELSE
                    617:                TRYRAC = .FALSE.
                    618:             END IF
                    619:             CALL ZSTEMR( JOBZ, 'A', N, RWORK( INDRDD ),
                    620:      $                   RWORK( INDREE ), VL, VU, IL, IU, M, W,
                    621:      $                   Z, LDZ, N, ISUPPZ, TRYRAC,
                    622:      $                   RWORK( INDRWK ), LLRWORK,
                    623:      $                   IWORK, LIWORK, INFO )
                    624: *
                    625: *           Apply unitary matrix used in reduction to tridiagonal
1.18      bertrand  626: *           form to eigenvectors returned by ZSTEMR.
1.1       bertrand  627: *
                    628:             IF( WANTZ .AND. INFO.EQ.0 ) THEN
                    629:                INDWKN = INDWK
                    630:                LLWRKN = LWORK - INDWKN + 1
                    631:                CALL ZUNMTR( 'L', UPLO, 'N', N, M, A, LDA,
                    632:      $                      WORK( INDTAU ), Z, LDZ, WORK( INDWKN ),
                    633:      $                      LLWRKN, IINFO )
                    634:             END IF
                    635:          END IF
                    636: *
                    637: *
                    638:          IF( INFO.EQ.0 ) THEN
                    639:             M = N
                    640:             GO TO 30
                    641:          END IF
                    642:          INFO = 0
                    643:       END IF
                    644: *
                    645: *     Otherwise, call DSTEBZ and, if eigenvectors are desired, ZSTEIN.
                    646: *     Also call DSTEBZ and ZSTEIN if ZSTEMR fails.
                    647: *
                    648:       IF( WANTZ ) THEN
                    649:          ORDER = 'B'
                    650:       ELSE
                    651:          ORDER = 'E'
                    652:       END IF
                    653: 
                    654:       CALL DSTEBZ( RANGE, ORDER, N, VLL, VUU, IL, IU, ABSTLL,
                    655:      $             RWORK( INDRD ), RWORK( INDRE ), M, NSPLIT, W,
                    656:      $             IWORK( INDIBL ), IWORK( INDISP ), RWORK( INDRWK ),
                    657:      $             IWORK( INDIWO ), INFO )
                    658: *
                    659:       IF( WANTZ ) THEN
                    660:          CALL ZSTEIN( N, RWORK( INDRD ), RWORK( INDRE ), M, W,
                    661:      $                IWORK( INDIBL ), IWORK( INDISP ), Z, LDZ,
                    662:      $                RWORK( INDRWK ), IWORK( INDIWO ), IWORK( INDIFL ),
                    663:      $                INFO )
                    664: *
                    665: *        Apply unitary matrix used in reduction to tridiagonal
                    666: *        form to eigenvectors returned by ZSTEIN.
                    667: *
                    668:          INDWKN = INDWK
                    669:          LLWRKN = LWORK - INDWKN + 1
                    670:          CALL ZUNMTR( 'L', UPLO, 'N', N, M, A, LDA, WORK( INDTAU ), Z,
                    671:      $                LDZ, WORK( INDWKN ), LLWRKN, IINFO )
                    672:       END IF
                    673: *
                    674: *     If matrix was scaled, then rescale eigenvalues appropriately.
                    675: *
                    676:    30 CONTINUE
                    677:       IF( ISCALE.EQ.1 ) THEN
                    678:          IF( INFO.EQ.0 ) THEN
                    679:             IMAX = M
                    680:          ELSE
                    681:             IMAX = INFO - 1
                    682:          END IF
                    683:          CALL DSCAL( IMAX, ONE / SIGMA, W, 1 )
                    684:       END IF
                    685: *
                    686: *     If eigenvalues are not in order, then sort them, along with
                    687: *     eigenvectors.
                    688: *
                    689:       IF( WANTZ ) THEN
                    690:          DO 50 J = 1, M - 1
                    691:             I = 0
                    692:             TMP1 = W( J )
                    693:             DO 40 JJ = J + 1, M
                    694:                IF( W( JJ ).LT.TMP1 ) THEN
                    695:                   I = JJ
                    696:                   TMP1 = W( JJ )
                    697:                END IF
                    698:    40       CONTINUE
                    699: *
                    700:             IF( I.NE.0 ) THEN
                    701:                ITMP1 = IWORK( INDIBL+I-1 )
                    702:                W( I ) = W( J )
                    703:                IWORK( INDIBL+I-1 ) = IWORK( INDIBL+J-1 )
                    704:                W( J ) = TMP1
                    705:                IWORK( INDIBL+J-1 ) = ITMP1
                    706:                CALL ZSWAP( N, Z( 1, I ), 1, Z( 1, J ), 1 )
                    707:             END IF
                    708:    50    CONTINUE
                    709:       END IF
                    710: *
                    711: *     Set WORK(1) to optimal workspace size.
                    712: *
                    713:       WORK( 1 ) = LWKOPT
                    714:       RWORK( 1 ) = LRWMIN
                    715:       IWORK( 1 ) = LIWMIN
                    716: *
                    717:       RETURN
                    718: *
                    719: *     End of ZHEEVR
                    720: *
                    721:       END

CVSweb interface <joel.bertrand@systella.fr>