SRC/dsbev.f

   1 *> \brief <b> DSBEV 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 DSBEV + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dsbev.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dsbev.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dsbev.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE DSBEV( JOBZ, UPLO, N, KD, AB, LDAB, W, Z, LDZ, WORK,
  22 *                         INFO )
  23 *
  24 *       .. Scalar Arguments ..
  25 *       CHARACTER          JOBZ, UPLO
  26 *       INTEGER            INFO, KD, LDAB, LDZ, N
  27 *       ..
  28 *       .. Array Arguments ..
  29 *       DOUBLE PRECISION   AB( LDAB, * ), W( * ), WORK( * ), Z( LDZ, * )
  30 *       ..
  31 *
  32 *
  33 *> \par Purpose:
  34 *  =============
  35 *>
  36 *> \verbatim
  37 *>
  38 *> DSBEV computes all the eigenvalues and, optionally, eigenvectors of
  39 *> a real symmetric band matrix A.
  40 *> \endverbatim
  41 *
  42 *  Arguments:
  43 *  ==========
  44 *
  45 *> \param[in] JOBZ
  46 *> \verbatim
  47 *>          JOBZ is CHARACTER*1
  48 *>          = 'N':  Compute eigenvalues only;
  49 *>          = 'V':  Compute eigenvalues and eigenvectors.
  50 *> \endverbatim
  51 *>
  52 *> \param[in] UPLO
  53 *> \verbatim
  54 *>          UPLO is CHARACTER*1
  55 *>          = 'U':  Upper triangle of A is stored;
  56 *>          = 'L':  Lower triangle of A is stored.
  57 *> \endverbatim
  58 *>
  59 *> \param[in] N
  60 *> \verbatim
  61 *>          N is INTEGER
  62 *>          The order of the matrix A.  N >= 0.
  63 *> \endverbatim
  64 *>
  65 *> \param[in] KD
  66 *> \verbatim
  67 *>          KD is INTEGER
  68 *>          The number of superdiagonals of the matrix A if UPLO = 'U',
  69 *>          or the number of subdiagonals if UPLO = 'L'.  KD >= 0.
  70 *> \endverbatim
  71 *>
  72 *> \param[in,out] AB
  73 *> \verbatim
  74 *>          AB is DOUBLE PRECISION array, dimension (LDAB, N)
  75 *>          On entry, the upper or lower triangle of the symmetric band
  76 *>          matrix A, stored in the first KD+1 rows of the array.  The
  77 *>          j-th column of A is stored in the j-th column of the array AB
  78 *>          as follows:
  79 *>          if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j;
  80 *>          if UPLO = 'L', AB(1+i-j,j)    = A(i,j) for j<=i<=min(n,j+kd).
  81 *>
  82 *>          On exit, AB is overwritten by values generated during the
  83 *>          reduction to tridiagonal form.  If UPLO = 'U', the first
  84 *>          superdiagonal and the diagonal of the tridiagonal matrix T
  85 *>          are returned in rows KD and KD+1 of AB, and if UPLO = 'L',
  86 *>          the diagonal and first subdiagonal of T are returned in the
  87 *>          first two rows of AB.
  88 *> \endverbatim
  89 *>
  90 *> \param[in] LDAB
  91 *> \verbatim
  92 *>          LDAB is INTEGER
  93 *>          The leading dimension of the array AB.  LDAB >= KD + 1.
  94 *> \endverbatim
  95 *>
  96 *> \param[out] W
  97 *> \verbatim
  98 *>          W is DOUBLE PRECISION array, dimension (N)
  99 *>          If INFO = 0, the eigenvalues in ascending order.
 100 *> \endverbatim
 101 *>
 102 *> \param[out] Z
 103 *> \verbatim
 104 *>          Z is DOUBLE PRECISION array, dimension (LDZ, N)
 105 *>          If JOBZ = 'V', then if INFO = 0, Z contains the orthonormal
 106 *>          eigenvectors of the matrix A, with the i-th column of Z
 107 *>          holding the eigenvector associated with W(i).
 108 *>          If JOBZ = 'N', then Z is not referenced.
 109 *> \endverbatim
 110 *>
 111 *> \param[in] LDZ
 112 *> \verbatim
 113 *>          LDZ is INTEGER
 114 *>          The leading dimension of the array Z.  LDZ >= 1, and if
 115 *>          JOBZ = 'V', LDZ >= max(1,N).
 116 *> \endverbatim
 117 *>
 118 *> \param[out] WORK
 119 *> \verbatim
 120 *>          WORK is DOUBLE PRECISION array, dimension (max(1,3*N-2))
 121 *> \endverbatim
 122 *>
 123 *> \param[out] INFO
 124 *> \verbatim
 125 *>          INFO is INTEGER
 126 *>          = 0:  successful exit
 127 *>          < 0:  if INFO = -i, the i-th argument had an illegal value
 128 *>          > 0:  if INFO = i, the algorithm failed to converge; i
 129 *>                off-diagonal elements of an intermediate tridiagonal
 130 *>                form did not converge to zero.
 131 *> \endverbatim
 132 *
 133 *  Authors:
 134 *  ========
 135 *
 136 *> \author Univ. of Tennessee
 137 *> \author Univ. of California Berkeley
 138 *> \author Univ. of Colorado Denver
 139 *> \author NAG Ltd.
 140 *
 141 *> \date November 2011
 142 *
 143 *> \ingroup doubleOTHEReigen
 144 *
 145 *  =====================================================================
 146       SUBROUTINE DSBEV( JOBZ, UPLO, N, KD, AB, LDAB, W, Z, LDZ, WORK,
 147      $                  INFO )
 148 *
 149 *  -- LAPACK driver routine (version 3.4.0) --
 150 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 151 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 152 *     November 2011
 153 *
 154 *     .. Scalar Arguments ..
 155       CHARACTER          JOBZ, UPLO
 156       INTEGER            INFO, KD, LDAB, LDZ, N
 157 *     ..
 158 *     .. Array Arguments ..
 159       DOUBLE PRECISION   AB( LDAB, * ), W( * ), WORK( * ), Z( LDZ, * )
 160 *     ..
 161 *
 162 *  =====================================================================
 163 *
 164 *     .. Parameters ..
 165       DOUBLE PRECISION   ZERO, ONE
 166       PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0 )
 167 *     ..
 168 *     .. Local Scalars ..
 169       LOGICAL            LOWER, WANTZ
 170       INTEGER            IINFO, IMAX, INDE, INDWRK, ISCALE
 171       DOUBLE PRECISION   ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA,
 172      $                   SMLNUM
 173 *     ..
 174 *     .. External Functions ..
 175       LOGICAL            LSAME
 176       DOUBLE PRECISION   DLAMCH, DLANSB
 177       EXTERNAL           LSAME, DLAMCH, DLANSB
 178 *     ..
 179 *     .. External Subroutines ..
 180       EXTERNAL           DLASCL, DSBTRD, DSCAL, DSTEQR, DSTERF, XERBLA
 181 *     ..
 182 *     .. Intrinsic Functions ..
 183       INTRINSIC          SQRT
 184 *     ..
 185 *     .. Executable Statements ..
 186 *
 187 *     Test the input parameters.
 188 *
 189       WANTZ = LSAME( JOBZ, 'V' )
 190       LOWER = LSAME( UPLO, 'L' )
 191 *
 192       INFO = 0
 193       IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN
 194          INFO = -1
 195       ELSE IF( .NOT.( LOWER .OR. LSAME( UPLO, 'U' ) ) ) THEN
 196          INFO = -2
 197       ELSE IF( N.LT.0 ) THEN
 198          INFO = -3
 199       ELSE IF( KD.LT.0 ) THEN
 200          INFO = -4
 201       ELSE IF( LDAB.LT.KD+1 ) THEN
 202          INFO = -6
 203       ELSE IF( LDZ.LT.1 .OR. ( WANTZ .AND. LDZ.LT.N ) ) THEN
 204          INFO = -9
 205       END IF
 206 *
 207       IF( INFO.NE.0 ) THEN
 208          CALL XERBLA( 'DSBEV ', -INFO )
 209          RETURN
 210       END IF
 211 *
 212 *     Quick return if possible
 213 *
 214       IF( N.EQ.0 )
 215      $   RETURN
 216 *
 217       IF( N.EQ.1 ) THEN
 218          IF( LOWER ) THEN
 219             W( 1 ) = AB( 1, 1 )
 220          ELSE
 221             W( 1 ) = AB( KD+1, 1 )
 222          END IF
 223          IF( WANTZ )
 224      $      Z( 1, 1 ) = ONE
 225          RETURN
 226       END IF
 227 *
 228 *     Get machine constants.
 229 *
 230       SAFMIN = DLAMCH( 'Safe minimum' )
 231       EPS = DLAMCH( 'Precision' )
 232       SMLNUM = SAFMIN / EPS
 233       BIGNUM = ONE / SMLNUM
 234       RMIN = SQRT( SMLNUM )
 235       RMAX = SQRT( BIGNUM )
 236 *
 237 *     Scale matrix to allowable range, if necessary.
 238 *
 239       ANRM = DLANSB( 'M', UPLO, N, KD, AB, LDAB, WORK )
 240       ISCALE = 0
 241       IF( ANRM.GT.ZERO .AND. ANRM.LT.RMIN ) THEN
 242          ISCALE = 1
 243          SIGMA = RMIN / ANRM
 244       ELSE IF( ANRM.GT.RMAX ) THEN
 245          ISCALE = 1
 246          SIGMA = RMAX / ANRM
 247       END IF
 248       IF( ISCALE.EQ.1 ) THEN
 249          IF( LOWER ) THEN
 250             CALL DLASCL( 'B', KD, KD, ONE, SIGMA, N, N, AB, LDAB, INFO )
 251          ELSE
 252             CALL DLASCL( 'Q', KD, KD, ONE, SIGMA, N, N, AB, LDAB, INFO )
 253          END IF
 254       END IF
 255 *
 256 *     Call DSBTRD to reduce symmetric band matrix to tridiagonal form.
 257 *
 258       INDE = 1
 259       INDWRK = INDE + N
 260       CALL DSBTRD( JOBZ, UPLO, N, KD, AB, LDAB, W, WORK( INDE ), Z, LDZ,
 261      $             WORK( INDWRK ), IINFO )
 262 *
 263 *     For eigenvalues only, call DSTERF.  For eigenvectors, call SSTEQR.
 264 *
 265       IF( .NOT.WANTZ ) THEN
 266          CALL DSTERF( N, W, WORK( INDE ), INFO )
 267       ELSE
 268          CALL DSTEQR( JOBZ, N, W, WORK( INDE ), Z, LDZ, WORK( INDWRK ),
 269      $                INFO )
 270       END IF
 271 *
 272 *     If matrix was scaled, then rescale eigenvalues appropriately.
 273 *
 274       IF( ISCALE.EQ.1 ) THEN
 275          IF( INFO.EQ.0 ) THEN
 276             IMAX = N
 277          ELSE
 278             IMAX = INFO - 1
 279          END IF
 280          CALL DSCAL( IMAX, ONE / SIGMA, W, 1 )
 281       END IF
 282 *
 283       RETURN
 284 *
 285 *     End of DSBEV
 286 *
 287       END