Return to dsyevd.f CVS log

Up to [local] / rpl / lapack / lapack

Initial revision

    1:       SUBROUTINE DSYEVD( JOBZ, UPLO, N, A, LDA, W, WORK, LWORK, IWORK,
    2:      $                   LIWORK, 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, LIWORK, LWORK, N
   12: *     ..
   13: *     .. Array Arguments ..
   14:       INTEGER            IWORK( * )
   15:       DOUBLE PRECISION   A( LDA, * ), W( * ), WORK( * )
   16: *     ..
   17: *
   18: *  Purpose
   19: *  =======
   20: *
   21: *  DSYEVD computes all eigenvalues and, optionally, eigenvectors of a
   22: *  real symmetric matrix A. If eigenvectors are desired, it uses a
   23: *  divide and conquer algorithm.
   24: *
   25: *  The divide and conquer algorithm makes very mild assumptions about
   26: *  floating point arithmetic. It will work on machines with a guard
   27: *  digit in add/subtract, or on those binary machines without guard
   28: *  digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or
   29: *  Cray-2. It could conceivably fail on hexadecimal or decimal machines
   30: *  without guard digits, but we know of none.
   31: *
   32: *  Because of large use of BLAS of level 3, DSYEVD needs N**2 more
   33: *  workspace than DSYEVX.
   34: *
   35: *  Arguments
   36: *  =========
   37: *
   38: *  JOBZ    (input) CHARACTER*1
   39: *          = 'N':  Compute eigenvalues only;
   40: *          = 'V':  Compute eigenvalues and eigenvectors.
   41: *
   42: *  UPLO    (input) CHARACTER*1
   43: *          = 'U':  Upper triangle of A is stored;
   44: *          = 'L':  Lower triangle of A is stored.
   45: *
   46: *  N       (input) INTEGER
   47: *          The order of the matrix A.  N >= 0.
   48: *
   49: *  A       (input/output) DOUBLE PRECISION array, dimension (LDA, N)
   50: *          On entry, the symmetric matrix A.  If UPLO = 'U', the
   51: *          leading N-by-N upper triangular part of A contains the
   52: *          upper triangular part of the matrix A.  If UPLO = 'L',
   53: *          the leading N-by-N lower triangular part of A contains
   54: *          the lower triangular part of the matrix A.
   55: *          On exit, if JOBZ = 'V', then if INFO = 0, A contains the
   56: *          orthonormal eigenvectors of the matrix A.
   57: *          If JOBZ = 'N', then on exit the lower triangle (if UPLO='L')
   58: *          or the upper triangle (if UPLO='U') of A, including the
   59: *          diagonal, is destroyed.
   60: *
   61: *  LDA     (input) INTEGER
   62: *          The leading dimension of the array A.  LDA >= max(1,N).
   63: *
   64: *  W       (output) DOUBLE PRECISION array, dimension (N)
   65: *          If INFO = 0, the eigenvalues in ascending order.
   66: *
   67: *  WORK    (workspace/output) DOUBLE PRECISION array,
   68: *                                         dimension (LWORK)
   69: *          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
   70: *
   71: *  LWORK   (input) INTEGER
   72: *          The dimension of the array WORK.
   73: *          If N <= 1,               LWORK must be at least 1.
   74: *          If JOBZ = 'N' and N > 1, LWORK must be at least 2*N+1.
   75: *          If JOBZ = 'V' and N > 1, LWORK must be at least
   76: *                                                1 + 6*N + 2*N**2.
   77: *
   78: *          If LWORK = -1, then a workspace query is assumed; the routine
   79: *          only calculates the optimal sizes of the WORK and IWORK
   80: *          arrays, returns these values as the first entries of the WORK
   81: *          and IWORK arrays, and no error message related to LWORK or
   82: *          LIWORK is issued by XERBLA.
   83: *
   84: *  IWORK   (workspace/output) INTEGER array, dimension (MAX(1,LIWORK))
   85: *          On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK.
   86: *
   87: *  LIWORK  (input) INTEGER
   88: *          The dimension of the array IWORK.
   89: *          If N <= 1,                LIWORK must be at least 1.
   90: *          If JOBZ  = 'N' and N > 1, LIWORK must be at least 1.
   91: *          If JOBZ  = 'V' and N > 1, LIWORK must be at least 3 + 5*N.
   92: *
   93: *          If LIWORK = -1, then a workspace query is assumed; the
   94: *          routine only calculates the optimal sizes of the WORK and
   95: *          IWORK arrays, returns these values as the first entries of
   96: *          the WORK and IWORK arrays, and no error message related to
   97: *          LWORK or LIWORK is issued by XERBLA.
   98: *
   99: *  INFO    (output) INTEGER
  100: *          = 0:  successful exit
  101: *          < 0:  if INFO = -i, the i-th argument had an illegal value
  102: *          > 0:  if INFO = i and JOBZ = 'N', then the algorithm failed
  103: *                to converge; i off-diagonal elements of an intermediate
  104: *                tridiagonal form did not converge to zero;
  105: *                if INFO = i and JOBZ = 'V', then the algorithm failed
  106: *                to compute an eigenvalue while working on the submatrix
  107: *                lying in rows and columns INFO/(N+1) through
  108: *                mod(INFO,N+1).
  109: *
  110: *  Further Details
  111: *  ===============
  112: *
  113: *  Based on contributions by
  114: *     Jeff Rutter, Computer Science Division, University of California
  115: *     at Berkeley, USA
  116: *  Modified by Francoise Tisseur, University of Tennessee.
  117: *
  118: *  Modified description of INFO. Sven, 16 Feb 05.
  119: *  =====================================================================
  120: *
  121: *     .. Parameters ..
  122:       DOUBLE PRECISION   ZERO, ONE
  123:       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
  124: *     ..
  125: *     .. Local Scalars ..
  126: *
  127:       LOGICAL            LOWER, LQUERY, WANTZ
  128:       INTEGER            IINFO, INDE, INDTAU, INDWK2, INDWRK, ISCALE,
  129:      $                   LIOPT, LIWMIN, LLWORK, LLWRK2, LOPT, LWMIN
  130:       DOUBLE PRECISION   ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA,
  131:      $                   SMLNUM
  132: *     ..
  133: *     .. External Functions ..
  134:       LOGICAL            LSAME
  135:       INTEGER            ILAENV
  136:       DOUBLE PRECISION   DLAMCH, DLANSY
  137:       EXTERNAL           LSAME, DLAMCH, DLANSY, ILAENV
  138: *     ..
  139: *     .. External Subroutines ..
  140:       EXTERNAL           DLACPY, DLASCL, DORMTR, DSCAL, DSTEDC, DSTERF,
  141:      $                   DSYTRD, XERBLA
  142: *     ..
  143: *     .. Intrinsic Functions ..
  144:       INTRINSIC          MAX, SQRT
  145: *     ..
  146: *     .. Executable Statements ..
  147: *
  148: *     Test the input parameters.
  149: *
  150:       WANTZ = LSAME( JOBZ, 'V' )
  151:       LOWER = LSAME( UPLO, 'L' )
  152:       LQUERY = ( LWORK.EQ.-1 .OR. LIWORK.EQ.-1 )
  153: *
  154:       INFO = 0
  155:       IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN
  156:          INFO = -1
  157:       ELSE IF( .NOT.( LOWER .OR. LSAME( UPLO, 'U' ) ) ) THEN
  158:          INFO = -2
  159:       ELSE IF( N.LT.0 ) THEN
  160:          INFO = -3
  161:       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
  162:          INFO = -5
  163:       END IF
  164: *
  165:       IF( INFO.EQ.0 ) THEN
  166:          IF( N.LE.1 ) THEN
  167:             LIWMIN = 1
  168:             LWMIN = 1
  169:             LOPT = LWMIN
  170:             LIOPT = LIWMIN
  171:          ELSE
  172:             IF( WANTZ ) THEN
  173:                LIWMIN = 3 + 5*N
  174:                LWMIN = 1 + 6*N + 2*N**2
  175:             ELSE
  176:                LIWMIN = 1
  177:                LWMIN = 2*N + 1
  178:             END IF
  179:             LOPT = MAX( LWMIN, 2*N +
  180:      $                  ILAENV( 1, 'DSYTRD', UPLO, N, -1, -1, -1 ) )
  181:             LIOPT = LIWMIN
  182:          END IF
  183:          WORK( 1 ) = LOPT
  184:          IWORK( 1 ) = LIOPT
  185: *
  186:          IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN
  187:             INFO = -8
  188:          ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT.LQUERY ) THEN
  189:             INFO = -10
  190:          END IF
  191:       END IF
  192: *
  193:       IF( INFO.NE.0 ) THEN
  194:          CALL XERBLA( 'DSYEVD', -INFO )
  195:          RETURN
  196:       ELSE IF( LQUERY ) THEN
  197:          RETURN
  198:       END IF
  199: *
  200: *     Quick return if possible
  201: *
  202:       IF( N.EQ.0 )
  203:      $   RETURN
  204: *
  205:       IF( N.EQ.1 ) THEN
  206:          W( 1 ) = A( 1, 1 )
  207:          IF( WANTZ )
  208:      $      A( 1, 1 ) = ONE
  209:          RETURN
  210:       END IF
  211: *
  212: *     Get machine constants.
  213: *
  214:       SAFMIN = DLAMCH( 'Safe minimum' )
  215:       EPS = DLAMCH( 'Precision' )
  216:       SMLNUM = SAFMIN / EPS
  217:       BIGNUM = ONE / SMLNUM
  218:       RMIN = SQRT( SMLNUM )
  219:       RMAX = SQRT( BIGNUM )
  220: *
  221: *     Scale matrix to allowable range, if necessary.
  222: *
  223:       ANRM = DLANSY( 'M', UPLO, N, A, LDA, WORK )
  224:       ISCALE = 0
  225:       IF( ANRM.GT.ZERO .AND. ANRM.LT.RMIN ) THEN
  226:          ISCALE = 1
  227:          SIGMA = RMIN / ANRM
  228:       ELSE IF( ANRM.GT.RMAX ) THEN
  229:          ISCALE = 1
  230:          SIGMA = RMAX / ANRM
  231:       END IF
  232:       IF( ISCALE.EQ.1 )
  233:      $   CALL DLASCL( UPLO, 0, 0, ONE, SIGMA, N, N, A, LDA, INFO )
  234: *
  235: *     Call DSYTRD to reduce symmetric matrix to tridiagonal form.
  236: *
  237:       INDE = 1
  238:       INDTAU = INDE + N
  239:       INDWRK = INDTAU + N
  240:       LLWORK = LWORK - INDWRK + 1
  241:       INDWK2 = INDWRK + N*N
  242:       LLWRK2 = LWORK - INDWK2 + 1
  243: *
  244:       CALL DSYTRD( UPLO, N, A, LDA, W, WORK( INDE ), WORK( INDTAU ),
  245:      $             WORK( INDWRK ), LLWORK, IINFO )
  246:       LOPT = 2*N + WORK( INDWRK )
  247: *
  248: *     For eigenvalues only, call DSTERF.  For eigenvectors, first call
  249: *     DSTEDC to generate the eigenvector matrix, WORK(INDWRK), of the
  250: *     tridiagonal matrix, then call DORMTR to multiply it by the
  251: *     Householder transformations stored in A.
  252: *
  253:       IF( .NOT.WANTZ ) THEN
  254:          CALL DSTERF( N, W, WORK( INDE ), INFO )
  255:       ELSE
  256:          CALL DSTEDC( 'I', N, W, WORK( INDE ), WORK( INDWRK ), N,
  257:      $                WORK( INDWK2 ), LLWRK2, IWORK, LIWORK, INFO )
  258:          CALL DORMTR( 'L', UPLO, 'N', N, N, A, LDA, WORK( INDTAU ),
  259:      $                WORK( INDWRK ), N, WORK( INDWK2 ), LLWRK2, IINFO )
  260:          CALL DLACPY( 'A', N, N, WORK( INDWRK ), N, A, LDA )
  261:          LOPT = MAX( LOPT, 1+6*N+2*N**2 )
  262:       END IF
  263: *
  264: *     If matrix was scaled, then rescale eigenvalues appropriately.
  265: *
  266:       IF( ISCALE.EQ.1 )
  267:      $   CALL DSCAL( N, ONE / SIGMA, W, 1 )
  268: *
  269:       WORK( 1 ) = LOPT
  270:       IWORK( 1 ) = LIOPT
  271: *
  272:       RETURN
  273: *
  274: *     End of DSYEVD
  275: *
  276:       END