SRC/dstevd.f

   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 November 2011
 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.4.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 *     November 2011
 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