Return to zheev.f CVS log

Up to [local] / rpl / lapack / lapack

Cohérence

    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