File:  [local] / rpl / lapack / lapack / dstevd.f
Revision 1.17: download - view: text, annotated - select for diffs - revision graph
Mon Aug 7 08:39:07 2023 UTC (8 months, 3 weeks ago) by bertrand
Branches: MAIN
CVS tags: rpl-4_1_35, rpl-4_1_34, HEAD
Première mise à jour de lapack et blas.

    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: *> \ingroup doubleOTHEReigen
  159: *
  160: *  =====================================================================
  161:       SUBROUTINE DSTEVD( JOBZ, N, D, E, Z, LDZ, WORK, LWORK, IWORK,
  162:      $                   LIWORK, INFO )
  163: *
  164: *  -- LAPACK driver routine --
  165: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
  166: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  167: *
  168: *     .. Scalar Arguments ..
  169:       CHARACTER          JOBZ
  170:       INTEGER            INFO, LDZ, LIWORK, LWORK, N
  171: *     ..
  172: *     .. Array Arguments ..
  173:       INTEGER            IWORK( * )
  174:       DOUBLE PRECISION   D( * ), E( * ), WORK( * ), Z( LDZ, * )
  175: *     ..
  176: *
  177: *  =====================================================================
  178: *
  179: *     .. Parameters ..
  180:       DOUBLE PRECISION   ZERO, ONE
  181:       PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0 )
  182: *     ..
  183: *     .. Local Scalars ..
  184:       LOGICAL            LQUERY, WANTZ
  185:       INTEGER            ISCALE, LIWMIN, LWMIN
  186:       DOUBLE PRECISION   BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA, SMLNUM,
  187:      $                   TNRM
  188: *     ..
  189: *     .. External Functions ..
  190:       LOGICAL            LSAME
  191:       DOUBLE PRECISION   DLAMCH, DLANST
  192:       EXTERNAL           LSAME, DLAMCH, DLANST
  193: *     ..
  194: *     .. External Subroutines ..
  195:       EXTERNAL           DSCAL, DSTEDC, DSTERF, XERBLA
  196: *     ..
  197: *     .. Intrinsic Functions ..
  198:       INTRINSIC          SQRT
  199: *     ..
  200: *     .. Executable Statements ..
  201: *
  202: *     Test the input parameters.
  203: *
  204:       WANTZ = LSAME( JOBZ, 'V' )
  205:       LQUERY = ( LWORK.EQ.-1 .OR. LIWORK.EQ.-1 )
  206: *
  207:       INFO = 0
  208:       LIWMIN = 1
  209:       LWMIN = 1
  210:       IF( N.GT.1 .AND. WANTZ ) THEN
  211:          LWMIN = 1 + 4*N + N**2
  212:          LIWMIN = 3 + 5*N
  213:       END IF
  214: *
  215:       IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN
  216:          INFO = -1
  217:       ELSE IF( N.LT.0 ) THEN
  218:          INFO = -2
  219:       ELSE IF( LDZ.LT.1 .OR. ( WANTZ .AND. LDZ.LT.N ) ) THEN
  220:          INFO = -6
  221:       END IF
  222: *
  223:       IF( INFO.EQ.0 ) THEN
  224:          WORK( 1 ) = LWMIN
  225:          IWORK( 1 ) = LIWMIN
  226: *
  227:          IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN
  228:             INFO = -8
  229:          ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT.LQUERY ) THEN
  230:             INFO = -10
  231:          END IF
  232:       END IF
  233: *
  234:       IF( INFO.NE.0 ) THEN
  235:          CALL XERBLA( 'DSTEVD', -INFO )
  236:          RETURN
  237:       ELSE IF( LQUERY ) THEN
  238:          RETURN
  239:       END IF
  240: *
  241: *     Quick return if possible
  242: *
  243:       IF( N.EQ.0 )
  244:      $   RETURN
  245: *
  246:       IF( N.EQ.1 ) THEN
  247:          IF( WANTZ )
  248:      $      Z( 1, 1 ) = ONE
  249:          RETURN
  250:       END IF
  251: *
  252: *     Get machine constants.
  253: *
  254:       SAFMIN = DLAMCH( 'Safe minimum' )
  255:       EPS = DLAMCH( 'Precision' )
  256:       SMLNUM = SAFMIN / EPS
  257:       BIGNUM = ONE / SMLNUM
  258:       RMIN = SQRT( SMLNUM )
  259:       RMAX = SQRT( BIGNUM )
  260: *
  261: *     Scale matrix to allowable range, if necessary.
  262: *
  263:       ISCALE = 0
  264:       TNRM = DLANST( 'M', N, D, E )
  265:       IF( TNRM.GT.ZERO .AND. TNRM.LT.RMIN ) THEN
  266:          ISCALE = 1
  267:          SIGMA = RMIN / TNRM
  268:       ELSE IF( TNRM.GT.RMAX ) THEN
  269:          ISCALE = 1
  270:          SIGMA = RMAX / TNRM
  271:       END IF
  272:       IF( ISCALE.EQ.1 ) THEN
  273:          CALL DSCAL( N, SIGMA, D, 1 )
  274:          CALL DSCAL( N-1, SIGMA, E( 1 ), 1 )
  275:       END IF
  276: *
  277: *     For eigenvalues only, call DSTERF.  For eigenvalues and
  278: *     eigenvectors, call DSTEDC.
  279: *
  280:       IF( .NOT.WANTZ ) THEN
  281:          CALL DSTERF( N, D, E, INFO )
  282:       ELSE
  283:          CALL DSTEDC( 'I', N, D, E, Z, LDZ, WORK, LWORK, IWORK, LIWORK,
  284:      $                INFO )
  285:       END IF
  286: *
  287: *     If matrix was scaled, then rescale eigenvalues appropriately.
  288: *
  289:       IF( ISCALE.EQ.1 )
  290:      $   CALL DSCAL( N, ONE / SIGMA, D, 1 )
  291: *
  292:       WORK( 1 ) = LWMIN
  293:       IWORK( 1 ) = LIWMIN
  294: *
  295:       RETURN
  296: *
  297: *     End of DSTEVD
  298: *
  299:       END

CVSweb interface <joel.bertrand@systella.fr>