Return to dstevd.f CVS log

Up to [local] / rpl / lapack / lapack

Mise à jour de Lapack.

    1: *> \brief <b> DSTEVD computes the eigenvalues and, optionally, the left and/or right eigenvectors for OTHER matrices</b>
    2: *
    3: *  =========== DOCUMENTATION ===========
    4: *
    5: * Online html documentation available at
    6: *            http://www.netlib.org/lapack/explore-html/
    7: *
    8: *> \htmlonly
    9: *> Download DSTEVD + dependencies
   10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dstevd.f">
   11: *> [TGZ]</a>
   12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dstevd.f">
   13: *> [ZIP]</a>
   14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dstevd.f">
   15: *> [TXT]</a>
   16: *> \endhtmlonly
   17: *
   18: *  Definition:
   19: *  ===========
   20: *
   21: *       SUBROUTINE DSTEVD( JOBZ, N, D, E, Z, LDZ, WORK, LWORK, IWORK,
   22: *                          LIWORK, INFO )
   23: *
   24: *       .. Scalar Arguments ..
   25: *       CHARACTER          JOBZ
   26: *       INTEGER            INFO, LDZ, LIWORK, LWORK, N
   27: *       ..
   28: *       .. Array Arguments ..
   29: *       INTEGER            IWORK( * )
   30: *       DOUBLE PRECISION   D( * ), E( * ), WORK( * ), Z( LDZ, * )
   31: *       ..
   32: *
   33: *
   34: *> \par Purpose:
   35: *  =============
   36: *>
   37: *> \verbatim
   38: *>
   39: *> DSTEVD computes all eigenvalues and, optionally, eigenvectors of a
   40: *> real symmetric tridiagonal matrix. If eigenvectors are desired, it
   41: *> uses a divide and conquer algorithm.
   42: *>
   43: *> The divide and conquer algorithm makes very mild assumptions about
   44: *> floating point arithmetic. It will work on machines with a guard
   45: *> digit in add/subtract, or on those binary machines without guard
   46: *> digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or
   47: *> Cray-2. It could conceivably fail on hexadecimal or decimal machines
   48: *> without guard digits, but we know of none.
   49: *> \endverbatim
   50: *
   51: *  Arguments:
   52: *  ==========
   53: *
   54: *> \param[in] JOBZ
   55: *> \verbatim
   56: *>          JOBZ is CHARACTER*1
   57: *>          = 'N':  Compute eigenvalues only;
   58: *>          = 'V':  Compute eigenvalues and eigenvectors.
   59: *> \endverbatim
   60: *>
   61: *> \param[in] N
   62: *> \verbatim
   63: *>          N is INTEGER
   64: *>          The order of the matrix.  N >= 0.
   65: *> \endverbatim
   66: *>
   67: *> \param[in,out] D
   68: *> \verbatim
   69: *>          D is DOUBLE PRECISION array, dimension (N)
   70: *>          On entry, the n diagonal elements of the tridiagonal matrix
   71: *>          A.
   72: *>          On exit, if INFO = 0, the eigenvalues in ascending order.
   73: *> \endverbatim
   74: *>
   75: *> \param[in,out] E
   76: *> \verbatim
   77: *>          E is DOUBLE PRECISION array, dimension (N-1)
   78: *>          On entry, the (n-1) subdiagonal elements of the tridiagonal
   79: *>          matrix A, stored in elements 1 to N-1 of E.
   80: *>          On exit, the contents of E are destroyed.
   81: *> \endverbatim
   82: *>
   83: *> \param[out] Z
   84: *> \verbatim
   85: *>          Z is DOUBLE PRECISION array, dimension (LDZ, N)
   86: *>          If JOBZ = 'V', then if INFO = 0, Z contains the orthonormal
   87: *>          eigenvectors of the matrix A, with the i-th column of Z
   88: *>          holding the eigenvector associated with D(i).
   89: *>          If JOBZ = 'N', then Z is not referenced.
   90: *> \endverbatim
   91: *>
   92: *> \param[in] LDZ
   93: *> \verbatim
   94: *>          LDZ is INTEGER
   95: *>          The leading dimension of the array Z.  LDZ >= 1, and if
   96: *>          JOBZ = 'V', LDZ >= max(1,N).
   97: *> \endverbatim
   98: *>
   99: *> \param[out] WORK
  100: *> \verbatim
  101: *>          WORK is DOUBLE PRECISION array,
  102: *>                                         dimension (LWORK)
  103: *>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
  104: *> \endverbatim
  105: *>
  106: *> \param[in] LWORK
  107: *> \verbatim
  108: *>          LWORK is INTEGER
  109: *>          The dimension of the array WORK.
  110: *>          If JOBZ  = 'N' or N <= 1 then LWORK must be at least 1.
  111: *>          If JOBZ  = 'V' and N > 1 then LWORK must be at least
  112: *>                         ( 1 + 4*N + N**2 ).
  113: *>
  114: *>          If LWORK = -1, then a workspace query is assumed; the routine
  115: *>          only calculates the optimal sizes of the WORK and IWORK
  116: *>          arrays, returns these values as the first entries of the WORK
  117: *>          and IWORK arrays, and no error message related to LWORK or
  118: *>          LIWORK is issued by XERBLA.
  119: *> \endverbatim
  120: *>
  121: *> \param[out] IWORK
  122: *> \verbatim
  123: *>          IWORK is INTEGER array, dimension (MAX(1,LIWORK))
  124: *>          On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK.
  125: *> \endverbatim
  126: *>
  127: *> \param[in] LIWORK
  128: *> \verbatim
  129: *>          LIWORK is INTEGER
  130: *>          The dimension of the array IWORK.
  131: *>          If JOBZ  = 'N' or N <= 1 then LIWORK must be at least 1.
  132: *>          If JOBZ  = 'V' and N > 1 then LIWORK must be at least 3+5*N.
  133: *>
  134: *>          If LIWORK = -1, then a workspace query is assumed; the
  135: *>          routine only calculates the optimal sizes of the WORK and
  136: *>          IWORK arrays, returns these values as the first entries of
  137: *>          the WORK and IWORK arrays, and no error message related to
  138: *>          LWORK or LIWORK is issued by XERBLA.
  139: *> \endverbatim
  140: *>
  141: *> \param[out] INFO
  142: *> \verbatim
  143: *>          INFO is INTEGER
  144: *>          = 0:  successful exit
  145: *>          < 0:  if INFO = -i, the i-th argument had an illegal value
  146: *>          > 0:  if INFO = i, the algorithm failed to converge; i
  147: *>                off-diagonal elements of E did not converge to zero.
  148: *> \endverbatim
  149: *
  150: *  Authors:
  151: *  ========
  152: *
  153: *> \author Univ. of Tennessee
  154: *> \author Univ. of California Berkeley
  155: *> \author Univ. of Colorado Denver
  156: *> \author NAG Ltd.
  157: *
  158: *> \date December 2016
  159: *
  160: *> \ingroup doubleOTHEReigen
  161: *
  162: *  =====================================================================
  163:       SUBROUTINE DSTEVD( JOBZ, N, D, E, Z, LDZ, WORK, LWORK, IWORK,
  164:      $                   LIWORK, INFO )
  165: *
  166: *  -- LAPACK driver routine (version 3.7.0) --
  167: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
  168: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  169: *     December 2016
  170: *
  171: *     .. Scalar Arguments ..
  172:       CHARACTER          JOBZ
  173:       INTEGER            INFO, LDZ, LIWORK, LWORK, N
  174: *     ..
  175: *     .. Array Arguments ..
  176:       INTEGER            IWORK( * )
  177:       DOUBLE PRECISION   D( * ), E( * ), WORK( * ), Z( LDZ, * )
  178: *     ..
  179: *
  180: *  =====================================================================
  181: *
  182: *     .. Parameters ..
  183:       DOUBLE PRECISION   ZERO, ONE
  184:       PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0 )
  185: *     ..
  186: *     .. Local Scalars ..
  187:       LOGICAL            LQUERY, WANTZ
  188:       INTEGER            ISCALE, LIWMIN, LWMIN
  189:       DOUBLE PRECISION   BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA, SMLNUM,
  190:      $                   TNRM
  191: *     ..
  192: *     .. External Functions ..
  193:       LOGICAL            LSAME
  194:       DOUBLE PRECISION   DLAMCH, DLANST
  195:       EXTERNAL           LSAME, DLAMCH, DLANST
  196: *     ..
  197: *     .. External Subroutines ..
  198:       EXTERNAL           DSCAL, DSTEDC, DSTERF, XERBLA
  199: *     ..
  200: *     .. Intrinsic Functions ..
  201:       INTRINSIC          SQRT
  202: *     ..
  203: *     .. Executable Statements ..
  204: *
  205: *     Test the input parameters.
  206: *
  207:       WANTZ = LSAME( JOBZ, 'V' )
  208:       LQUERY = ( LWORK.EQ.-1 .OR. LIWORK.EQ.-1 )
  209: *
  210:       INFO = 0
  211:       LIWMIN = 1
  212:       LWMIN = 1
  213:       IF( N.GT.1 .AND. WANTZ ) THEN
  214:          LWMIN = 1 + 4*N + N**2
  215:          LIWMIN = 3 + 5*N
  216:       END IF
  217: *
  218:       IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN
  219:          INFO = -1
  220:       ELSE IF( N.LT.0 ) THEN
  221:          INFO = -2
  222:       ELSE IF( LDZ.LT.1 .OR. ( WANTZ .AND. LDZ.LT.N ) ) THEN
  223:          INFO = -6
  224:       END IF
  225: *
  226:       IF( INFO.EQ.0 ) THEN
  227:          WORK( 1 ) = LWMIN
  228:          IWORK( 1 ) = LIWMIN
  229: *
  230:          IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN
  231:             INFO = -8
  232:          ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT.LQUERY ) THEN
  233:             INFO = -10
  234:          END IF
  235:       END IF
  236: *
  237:       IF( INFO.NE.0 ) THEN
  238:          CALL XERBLA( 'DSTEVD', -INFO )
  239:          RETURN
  240:       ELSE IF( LQUERY ) THEN
  241:          RETURN
  242:       END IF
  243: *
  244: *     Quick return if possible
  245: *
  246:       IF( N.EQ.0 )
  247:      $   RETURN
  248: *
  249:       IF( N.EQ.1 ) THEN
  250:          IF( WANTZ )
  251:      $      Z( 1, 1 ) = ONE
  252:          RETURN
  253:       END IF
  254: *
  255: *     Get machine constants.
  256: *
  257:       SAFMIN = DLAMCH( 'Safe minimum' )
  258:       EPS = DLAMCH( 'Precision' )
  259:       SMLNUM = SAFMIN / EPS
  260:       BIGNUM = ONE / SMLNUM
  261:       RMIN = SQRT( SMLNUM )
  262:       RMAX = SQRT( BIGNUM )
  263: *
  264: *     Scale matrix to allowable range, if necessary.
  265: *
  266:       ISCALE = 0
  267:       TNRM = DLANST( 'M', N, D, E )
  268:       IF( TNRM.GT.ZERO .AND. TNRM.LT.RMIN ) THEN
  269:          ISCALE = 1
  270:          SIGMA = RMIN / TNRM
  271:       ELSE IF( TNRM.GT.RMAX ) THEN
  272:          ISCALE = 1
  273:          SIGMA = RMAX / TNRM
  274:       END IF
  275:       IF( ISCALE.EQ.1 ) THEN
  276:          CALL DSCAL( N, SIGMA, D, 1 )
  277:          CALL DSCAL( N-1, SIGMA, E( 1 ), 1 )
  278:       END IF
  279: *
  280: *     For eigenvalues only, call DSTERF.  For eigenvalues and
  281: *     eigenvectors, call DSTEDC.
  282: *
  283:       IF( .NOT.WANTZ ) THEN
  284:          CALL DSTERF( N, D, E, INFO )
  285:       ELSE
  286:          CALL DSTEDC( 'I', N, D, E, Z, LDZ, WORK, LWORK, IWORK, LIWORK,
  287:      $                INFO )
  288:       END IF
  289: *
  290: *     If matrix was scaled, then rescale eigenvalues appropriately.
  291: *
  292:       IF( ISCALE.EQ.1 )
  293:      $   CALL DSCAL( N, ONE / SIGMA, D, 1 )
  294: *
  295:       WORK( 1 ) = LWMIN
  296:       IWORK( 1 ) = LIWMIN
  297: *
  298:       RETURN
  299: *
  300: *     End of DSTEVD
  301: *
  302:       END