SRC/chbev.f

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