SRC/ssbevd.f

   1 *> \brief <b> SSBEVD 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 SSBEVD + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ssbevd.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ssbevd.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ssbevd.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE SSBEVD( JOBZ, UPLO, N, KD, AB, LDAB, W, Z, LDZ, WORK,
  22 *                          LWORK, IWORK, LIWORK, INFO )
  23 *
  24 *       .. Scalar Arguments ..
  25 *       CHARACTER          JOBZ, UPLO
  26 *       INTEGER            INFO, KD, LDAB, LDZ, LIWORK, LWORK, N
  27 *       ..
  28 *       .. Array Arguments ..
  29 *       INTEGER            IWORK( * )
  30 *       REAL               AB( LDAB, * ), W( * ), WORK( * ), Z( LDZ, * )
  31 *       ..
  32 *
  33 *
  34 *> \par Purpose:
  35 *  =============
  36 *>
  37 *> \verbatim
  38 *>
  39 *> SSBEVD computes all the eigenvalues and, optionally, eigenvectors of
  40 *> a real symmetric band matrix A. If eigenvectors are desired, it uses
  41 *> 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] UPLO
  62 *> \verbatim
  63 *>          UPLO is CHARACTER*1
  64 *>          = 'U':  Upper triangle of A is stored;
  65 *>          = 'L':  Lower triangle of A is stored.
  66 *> \endverbatim
  67 *>
  68 *> \param[in] N
  69 *> \verbatim
  70 *>          N is INTEGER
  71 *>          The order of the matrix A.  N >= 0.
  72 *> \endverbatim
  73 *>
  74 *> \param[in] KD
  75 *> \verbatim
  76 *>          KD is INTEGER
  77 *>          The number of superdiagonals of the matrix A if UPLO = 'U',
  78 *>          or the number of subdiagonals if UPLO = 'L'.  KD >= 0.
  79 *> \endverbatim
  80 *>
  81 *> \param[in,out] AB
  82 *> \verbatim
  83 *>          AB is REAL array, dimension (LDAB, N)
  84 *>          On entry, the upper or lower triangle of the symmetric band
  85 *>          matrix A, stored in the first KD+1 rows of the array.  The
  86 *>          j-th column of A is stored in the j-th column of the array AB
  87 *>          as follows:
  88 *>          if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j;
  89 *>          if UPLO = 'L', AB(1+i-j,j)    = A(i,j) for j<=i<=min(n,j+kd).
  90 *>
  91 *>          On exit, AB is overwritten by values generated during the
  92 *>          reduction to tridiagonal form.  If UPLO = 'U', the first
  93 *>          superdiagonal and the diagonal of the tridiagonal matrix T
  94 *>          are returned in rows KD and KD+1 of AB, and if UPLO = 'L',
  95 *>          the diagonal and first subdiagonal of T are returned in the
  96 *>          first two rows of AB.
  97 *> \endverbatim
  98 *>
  99 *> \param[in] LDAB
 100 *> \verbatim
 101 *>          LDAB is INTEGER
 102 *>          The leading dimension of the array AB.  LDAB >= KD + 1.
 103 *> \endverbatim
 104 *>
 105 *> \param[out] W
 106 *> \verbatim
 107 *>          W is REAL array, dimension (N)
 108 *>          If INFO = 0, the eigenvalues in ascending order.
 109 *> \endverbatim
 110 *>
 111 *> \param[out] Z
 112 *> \verbatim
 113 *>          Z is REAL array, dimension (LDZ, N)
 114 *>          If JOBZ = 'V', then if INFO = 0, Z contains the orthonormal
 115 *>          eigenvectors of the matrix A, with the i-th column of Z
 116 *>          holding the eigenvector associated with W(i).
 117 *>          If JOBZ = 'N', then Z is not referenced.
 118 *> \endverbatim
 119 *>
 120 *> \param[in] LDZ
 121 *> \verbatim
 122 *>          LDZ is INTEGER
 123 *>          The leading dimension of the array Z.  LDZ >= 1, and if
 124 *>          JOBZ = 'V', LDZ >= max(1,N).
 125 *> \endverbatim
 126 *>
 127 *> \param[out] WORK
 128 *> \verbatim
 129 *>          WORK is REAL array,
 130 *>                                         dimension (LWORK)
 131 *>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
 132 *> \endverbatim
 133 *>
 134 *> \param[in] LWORK
 135 *> \verbatim
 136 *>          LWORK is INTEGER
 137 *>          The dimension of the array WORK.
 138 *>          IF N <= 1,                LWORK must be at least 1.
 139 *>          If JOBZ  = 'N' and N > 2, LWORK must be at least 2*N.
 140 *>          If JOBZ  = 'V' and N > 2, LWORK must be at least
 141 *>                         ( 1 + 5*N + 2*N**2 ).
 142 *>
 143 *>          If LWORK = -1, then a workspace query is assumed; the routine
 144 *>          only calculates the optimal sizes of the WORK and IWORK
 145 *>          arrays, returns these values as the first entries of the WORK
 146 *>          and IWORK arrays, and no error message related to LWORK or
 147 *>          LIWORK is issued by XERBLA.
 148 *> \endverbatim
 149 *>
 150 *> \param[out] IWORK
 151 *> \verbatim
 152 *>          IWORK is INTEGER array, dimension (MAX(1,LIWORK))
 153 *>          On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK.
 154 *> \endverbatim
 155 *>
 156 *> \param[in] LIWORK
 157 *> \verbatim
 158 *>          LIWORK is INTEGER
 159 *>          The dimension of the array IWORK.
 160 *>          If JOBZ  = 'N' or N <= 1, LIWORK must be at least 1.
 161 *>          If JOBZ  = 'V' and N > 2, LIWORK must be at least 3 + 5*N.
 162 *>
 163 *>          If LIWORK = -1, then a workspace query is assumed; the
 164 *>          routine only calculates the optimal sizes of the WORK and
 165 *>          IWORK arrays, returns these values as the first entries of
 166 *>          the WORK and IWORK arrays, and no error message related to
 167 *>          LWORK or LIWORK is issued by XERBLA.
 168 *> \endverbatim
 169 *>
 170 *> \param[out] INFO
 171 *> \verbatim
 172 *>          INFO is INTEGER
 173 *>          = 0:  successful exit
 174 *>          < 0:  if INFO = -i, the i-th argument had an illegal value
 175 *>          > 0:  if INFO = i, the algorithm failed to converge; i
 176 *>                off-diagonal elements of an intermediate tridiagonal
 177 *>                form did not converge to zero.
 178 *> \endverbatim
 179 *
 180 *  Authors:
 181 *  ========
 182 *
 183 *> \author Univ. of Tennessee
 184 *> \author Univ. of California Berkeley
 185 *> \author Univ. of Colorado Denver
 186 *> \author NAG Ltd.
 187 *
 188 *> \date November 2011
 189 *
 190 *> \ingroup realOTHEReigen
 191 *
 192 *  =====================================================================
 193       SUBROUTINE SSBEVD( JOBZ, UPLO, N, KD, AB, LDAB, W, Z, LDZ, WORK,
 194      $                   LWORK, IWORK, LIWORK, INFO )
 195 *
 196 *  -- LAPACK driver routine (version 3.4.0) --
 197 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 198 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 199 *     November 2011
 200 *
 201 *     .. Scalar Arguments ..
 202       CHARACTER          JOBZ, UPLO
 203       INTEGER            INFO, KD, LDAB, LDZ, LIWORK, LWORK, N
 204 *     ..
 205 *     .. Array Arguments ..
 206       INTEGER            IWORK( * )
 207       REAL               AB( LDAB, * ), W( * ), WORK( * ), Z( LDZ, * )
 208 *     ..
 209 *
 210 *  =====================================================================
 211 *
 212 *     .. Parameters ..
 213       REAL               ZERO, ONE
 214       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
 215 *     ..
 216 *     .. Local Scalars ..
 217       LOGICAL            LOWER, LQUERY, WANTZ
 218       INTEGER            IINFO, INDE, INDWK2, INDWRK, ISCALE, LIWMIN,
 219      $                   LLWRK2, LWMIN
 220       REAL               ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA,
 221      $                   SMLNUM
 222 *     ..
 223 *     .. External Functions ..
 224       LOGICAL            LSAME
 225       REAL               SLAMCH, SLANSB
 226       EXTERNAL           LSAME, SLAMCH, SLANSB
 227 *     ..
 228 *     .. External Subroutines ..
 229       EXTERNAL           SGEMM, SLACPY, SLASCL, SSBTRD, SSCAL, SSTEDC,
 230      $                   SSTERF, XERBLA
 231 *     ..
 232 *     .. Intrinsic Functions ..
 233       INTRINSIC          SQRT
 234 *     ..
 235 *     .. Executable Statements ..
 236 *
 237 *     Test the input parameters.
 238 *
 239       WANTZ = LSAME( JOBZ, 'V' )
 240       LOWER = LSAME( UPLO, 'L' )
 241       LQUERY = ( LWORK.EQ.-1 .OR. LIWORK.EQ.-1 )
 242 *
 243       INFO = 0
 244       IF( N.LE.1 ) THEN
 245          LIWMIN = 1
 246          LWMIN = 1
 247       ELSE
 248          IF( WANTZ ) THEN
 249             LIWMIN = 3 + 5*N
 250             LWMIN = 1 + 5*N + 2*N**2
 251          ELSE
 252             LIWMIN = 1
 253             LWMIN = 2*N
 254          END IF
 255       END IF
 256       IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN
 257          INFO = -1
 258       ELSE IF( .NOT.( LOWER .OR. LSAME( UPLO, 'U' ) ) ) THEN
 259          INFO = -2
 260       ELSE IF( N.LT.0 ) THEN
 261          INFO = -3
 262       ELSE IF( KD.LT.0 ) THEN
 263          INFO = -4
 264       ELSE IF( LDAB.LT.KD+1 ) THEN
 265          INFO = -6
 266       ELSE IF( LDZ.LT.1 .OR. ( WANTZ .AND. LDZ.LT.N ) ) THEN
 267          INFO = -9
 268       END IF
 269 *
 270       IF( INFO.EQ.0 ) THEN
 271          WORK( 1 ) = LWMIN
 272          IWORK( 1 ) = LIWMIN
 273 *
 274          IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN
 275             INFO = -11
 276          ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT.LQUERY ) THEN
 277             INFO = -13
 278          END IF
 279       END IF
 280 *
 281       IF( INFO.NE.0 ) THEN
 282          CALL XERBLA( 'SSBEVD', -INFO )
 283          RETURN
 284       ELSE IF( LQUERY ) THEN
 285          RETURN
 286       END IF
 287 *
 288 *     Quick return if possible
 289 *
 290       IF( N.EQ.0 )
 291      $   RETURN
 292 *
 293       IF( N.EQ.1 ) THEN
 294          W( 1 ) = AB( 1, 1 )
 295          IF( WANTZ )
 296      $      Z( 1, 1 ) = ONE
 297          RETURN
 298       END IF
 299 *
 300 *     Get machine constants.
 301 *
 302       SAFMIN = SLAMCH( 'Safe minimum' )
 303       EPS = SLAMCH( 'Precision' )
 304       SMLNUM = SAFMIN / EPS
 305       BIGNUM = ONE / SMLNUM
 306       RMIN = SQRT( SMLNUM )
 307       RMAX = SQRT( BIGNUM )
 308 *
 309 *     Scale matrix to allowable range, if necessary.
 310 *
 311       ANRM = SLANSB( 'M', UPLO, N, KD, AB, LDAB, WORK )
 312       ISCALE = 0
 313       IF( ANRM.GT.ZERO .AND. ANRM.LT.RMIN ) THEN
 314          ISCALE = 1
 315          SIGMA = RMIN / ANRM
 316       ELSE IF( ANRM.GT.RMAX ) THEN
 317          ISCALE = 1
 318          SIGMA = RMAX / ANRM
 319       END IF
 320       IF( ISCALE.EQ.1 ) THEN
 321          IF( LOWER ) THEN
 322             CALL SLASCL( 'B', KD, KD, ONE, SIGMA, N, N, AB, LDAB, INFO )
 323          ELSE
 324             CALL SLASCL( 'Q', KD, KD, ONE, SIGMA, N, N, AB, LDAB, INFO )
 325          END IF
 326       END IF
 327 *
 328 *     Call SSBTRD to reduce symmetric band matrix to tridiagonal form.
 329 *
 330       INDE = 1
 331       INDWRK = INDE + N
 332       INDWK2 = INDWRK + N*N
 333       LLWRK2 = LWORK - INDWK2 + 1
 334       CALL SSBTRD( JOBZ, UPLO, N, KD, AB, LDAB, W, WORK( INDE ), Z, LDZ,
 335      $             WORK( INDWRK ), IINFO )
 336 *
 337 *     For eigenvalues only, call SSTERF.  For eigenvectors, call SSTEDC.
 338 *
 339       IF( .NOT.WANTZ ) THEN
 340          CALL SSTERF( N, W, WORK( INDE ), INFO )
 341       ELSE
 342          CALL SSTEDC( 'I', N, W, WORK( INDE ), WORK( INDWRK ), N,
 343      $                WORK( INDWK2 ), LLWRK2, IWORK, LIWORK, INFO )
 344          CALL SGEMM( 'N', 'N', N, N, N, ONE, Z, LDZ, WORK( INDWRK ), N,
 345      $               ZERO, WORK( INDWK2 ), N )
 346          CALL SLACPY( 'A', N, N, WORK( INDWK2 ), N, Z, LDZ )
 347       END IF
 348 *
 349 *     If matrix was scaled, then rescale eigenvalues appropriately.
 350 *
 351       IF( ISCALE.EQ.1 )
 352      $   CALL SSCAL( N, ONE / SIGMA, W, 1 )
 353 *
 354       WORK( 1 ) = LWMIN
 355       IWORK( 1 ) = LIWMIN
 356       RETURN
 357 *
 358 *     End of SSBEVD
 359 *
 360       END