Return to dsyev.f CVS log

Up to [local] / rpl / lapack / lapack

Patches pour OS/2

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