Annotation of rpl/lapack/lapack/zheev.f, revision 1.4

1.1       bertrand    1:       SUBROUTINE ZHEEV( JOBZ, UPLO, N, A, LDA, W, WORK, LWORK, RWORK,
                      2:      $                  INFO )
                      3: *
                      4: *  -- LAPACK driver routine (version 3.2) --
                      5: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
                      6: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
                      7: *     November 2006
                      8: *
                      9: *     .. Scalar Arguments ..
                     10:       CHARACTER          JOBZ, UPLO
                     11:       INTEGER            INFO, LDA, LWORK, N
                     12: *     ..
                     13: *     .. Array Arguments ..
                     14:       DOUBLE PRECISION   RWORK( * ), W( * )
                     15:       COMPLEX*16         A( LDA, * ), WORK( * )
                     16: *     ..
                     17: *
                     18: *  Purpose
                     19: *  =======
                     20: *
                     21: *  ZHEEV computes all eigenvalues and, optionally, eigenvectors of a
                     22: *  complex Hermitian matrix A.
                     23: *
                     24: *  Arguments
                     25: *  =========
                     26: *
                     27: *  JOBZ    (input) CHARACTER*1
                     28: *          = 'N':  Compute eigenvalues only;
                     29: *          = 'V':  Compute eigenvalues and eigenvectors.
                     30: *
                     31: *  UPLO    (input) CHARACTER*1
                     32: *          = 'U':  Upper triangle of A is stored;
                     33: *          = 'L':  Lower triangle of A is stored.
                     34: *
                     35: *  N       (input) INTEGER
                     36: *          The order of the matrix A.  N >= 0.
                     37: *
                     38: *  A       (input/output) COMPLEX*16 array, dimension (LDA, N)
                     39: *          On entry, the Hermitian matrix A.  If UPLO = 'U', the
                     40: *          leading N-by-N upper triangular part of A contains the
                     41: *          upper triangular part of the matrix A.  If UPLO = 'L',
                     42: *          the leading N-by-N lower triangular part of A contains
                     43: *          the lower triangular part of the matrix A.
                     44: *          On exit, if JOBZ = 'V', then if INFO = 0, A contains the
                     45: *          orthonormal eigenvectors of the matrix A.
                     46: *          If JOBZ = 'N', then on exit the lower triangle (if UPLO='L')
                     47: *          or the upper triangle (if UPLO='U') of A, including the
                     48: *          diagonal, is destroyed.
                     49: *
                     50: *  LDA     (input) INTEGER
                     51: *          The leading dimension of the array A.  LDA >= max(1,N).
                     52: *
                     53: *  W       (output) DOUBLE PRECISION array, dimension (N)
                     54: *          If INFO = 0, the eigenvalues in ascending order.
                     55: *
                     56: *  WORK    (workspace/output) COMPLEX*16 array, dimension (MAX(1,LWORK))
                     57: *          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
                     58: *
                     59: *  LWORK   (input) INTEGER
                     60: *          The length of the array WORK.  LWORK >= max(1,2*N-1).
                     61: *          For optimal efficiency, LWORK >= (NB+1)*N,
                     62: *          where NB is the blocksize for ZHETRD returned by ILAENV.
                     63: *
                     64: *          If LWORK = -1, then a workspace query is assumed; the routine
                     65: *          only calculates the optimal size of the WORK array, returns
                     66: *          this value as the first entry of the WORK array, and no error
                     67: *          message related to LWORK is issued by XERBLA.
                     68: *
                     69: *  RWORK   (workspace) DOUBLE PRECISION array, dimension (max(1, 3*N-2))
                     70: *
                     71: *  INFO    (output) INTEGER
                     72: *          = 0:  successful exit
                     73: *          < 0:  if INFO = -i, the i-th argument had an illegal value
                     74: *          > 0:  if INFO = i, the algorithm failed to converge; i
                     75: *                off-diagonal elements of an intermediate tridiagonal
                     76: *                form did not converge to zero.
                     77: *
                     78: *  =====================================================================
                     79: *
                     80: *     .. Parameters ..
                     81:       DOUBLE PRECISION   ZERO, ONE
                     82:       PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0 )
                     83:       COMPLEX*16         CONE
                     84:       PARAMETER          ( CONE = ( 1.0D0, 0.0D0 ) )
                     85: *     ..
                     86: *     .. Local Scalars ..
                     87:       LOGICAL            LOWER, LQUERY, WANTZ
                     88:       INTEGER            IINFO, IMAX, INDE, INDTAU, INDWRK, ISCALE,
                     89:      $                   LLWORK, LWKOPT, NB
                     90:       DOUBLE PRECISION   ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA,
                     91:      $                   SMLNUM
                     92: *     ..
                     93: *     .. External Functions ..
                     94:       LOGICAL            LSAME
                     95:       INTEGER            ILAENV
                     96:       DOUBLE PRECISION   DLAMCH, ZLANHE
                     97:       EXTERNAL           LSAME, ILAENV, DLAMCH, ZLANHE
                     98: *     ..
                     99: *     .. External Subroutines ..
                    100:       EXTERNAL           DSCAL, DSTERF, XERBLA, ZHETRD, ZLASCL, ZSTEQR,
                    101:      $                   ZUNGTR
                    102: *     ..
                    103: *     .. Intrinsic Functions ..
                    104:       INTRINSIC          MAX, SQRT
                    105: *     ..
                    106: *     .. Executable Statements ..
                    107: *
                    108: *     Test the input parameters.
                    109: *
                    110:       WANTZ = LSAME( JOBZ, 'V' )
                    111:       LOWER = LSAME( UPLO, 'L' )
                    112:       LQUERY = ( LWORK.EQ.-1 )
                    113: *
                    114:       INFO = 0
                    115:       IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN
                    116:          INFO = -1
                    117:       ELSE IF( .NOT.( LOWER .OR. LSAME( UPLO, 'U' ) ) ) THEN
                    118:          INFO = -2
                    119:       ELSE IF( N.LT.0 ) THEN
                    120:          INFO = -3
                    121:       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
                    122:          INFO = -5
                    123:       END IF
                    124: *
                    125:       IF( INFO.EQ.0 ) THEN
                    126:          NB = ILAENV( 1, 'ZHETRD', UPLO, N, -1, -1, -1 )
                    127:          LWKOPT = MAX( 1, ( NB+1 )*N )
                    128:          WORK( 1 ) = LWKOPT
                    129: *
                    130:          IF( LWORK.LT.MAX( 1, 2*N-1 ) .AND. .NOT.LQUERY )
                    131:      $      INFO = -8
                    132:       END IF
                    133: *
                    134:       IF( INFO.NE.0 ) THEN
                    135:          CALL XERBLA( 'ZHEEV ', -INFO )
                    136:          RETURN
                    137:       ELSE IF( LQUERY ) THEN
                    138:          RETURN
                    139:       END IF
                    140: *
                    141: *     Quick return if possible
                    142: *
                    143:       IF( N.EQ.0 ) THEN
                    144:          RETURN
                    145:       END IF
                    146: *
                    147:       IF( N.EQ.1 ) THEN
                    148:          W( 1 ) = A( 1, 1 )
                    149:          WORK( 1 ) = 1
                    150:          IF( WANTZ )
                    151:      $      A( 1, 1 ) = CONE
                    152:          RETURN
                    153:       END IF
                    154: *
                    155: *     Get machine constants.
                    156: *
                    157:       SAFMIN = DLAMCH( 'Safe minimum' )
                    158:       EPS = DLAMCH( 'Precision' )
                    159:       SMLNUM = SAFMIN / EPS
                    160:       BIGNUM = ONE / SMLNUM
                    161:       RMIN = SQRT( SMLNUM )
                    162:       RMAX = SQRT( BIGNUM )
                    163: *
                    164: *     Scale matrix to allowable range, if necessary.
                    165: *
                    166:       ANRM = ZLANHE( 'M', UPLO, N, A, LDA, RWORK )
                    167:       ISCALE = 0
                    168:       IF( ANRM.GT.ZERO .AND. ANRM.LT.RMIN ) THEN
                    169:          ISCALE = 1
                    170:          SIGMA = RMIN / ANRM
                    171:       ELSE IF( ANRM.GT.RMAX ) THEN
                    172:          ISCALE = 1
                    173:          SIGMA = RMAX / ANRM
                    174:       END IF
                    175:       IF( ISCALE.EQ.1 )
                    176:      $   CALL ZLASCL( UPLO, 0, 0, ONE, SIGMA, N, N, A, LDA, INFO )
                    177: *
                    178: *     Call ZHETRD to reduce Hermitian matrix to tridiagonal form.
                    179: *
                    180:       INDE = 1
                    181:       INDTAU = 1
                    182:       INDWRK = INDTAU + N
                    183:       LLWORK = LWORK - INDWRK + 1
                    184:       CALL ZHETRD( UPLO, N, A, LDA, W, RWORK( INDE ), WORK( INDTAU ),
                    185:      $             WORK( INDWRK ), LLWORK, IINFO )
                    186: *
                    187: *     For eigenvalues only, call DSTERF.  For eigenvectors, first call
                    188: *     ZUNGTR to generate the unitary matrix, then call ZSTEQR.
                    189: *
                    190:       IF( .NOT.WANTZ ) THEN
                    191:          CALL DSTERF( N, W, RWORK( INDE ), INFO )
                    192:       ELSE
                    193:          CALL ZUNGTR( UPLO, N, A, LDA, WORK( INDTAU ), WORK( INDWRK ),
                    194:      $                LLWORK, IINFO )
                    195:          INDWRK = INDE + N
                    196:          CALL ZSTEQR( JOBZ, N, W, RWORK( INDE ), A, LDA,
                    197:      $                RWORK( INDWRK ), INFO )
                    198:       END IF
                    199: *
                    200: *     If matrix was scaled, then rescale eigenvalues appropriately.
                    201: *
                    202:       IF( ISCALE.EQ.1 ) THEN
                    203:          IF( INFO.EQ.0 ) THEN
                    204:             IMAX = N
                    205:          ELSE
                    206:             IMAX = INFO - 1
                    207:          END IF
                    208:          CALL DSCAL( IMAX, ONE / SIGMA, W, 1 )
                    209:       END IF
                    210: *
                    211: *     Set WORK(1) to optimal complex workspace size.
                    212: *
                    213:       WORK( 1 ) = LWKOPT
                    214: *
                    215:       RETURN
                    216: *
                    217: *     End of ZHEEV
                    218: *
                    219:       END

CVSweb interface <joel.bertrand@systella.fr>