Return to zhbev.f CVS log

Up to [local] / rpl / lapack / lapack

Cohérence

    1:       SUBROUTINE ZHBEV( JOBZ, UPLO, N, KD, AB, LDAB, W, Z, LDZ, WORK,
    2:      $                  RWORK, INFO )
    3: *
    4: *  -- LAPACK driver routine (version 3.2) --
    5: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
    6: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
    7: *     November 2006
    8: *
    9: *     .. Scalar Arguments ..
   10:       CHARACTER          JOBZ, UPLO
   11:       INTEGER            INFO, KD, LDAB, LDZ, N
   12: *     ..
   13: *     .. Array Arguments ..
   14:       DOUBLE PRECISION   RWORK( * ), W( * )
   15:       COMPLEX*16         AB( LDAB, * ), WORK( * ), Z( LDZ, * )
   16: *     ..
   17: *
   18: *  Purpose
   19: *  =======
   20: *
   21: *  ZHBEV computes all the eigenvalues and, optionally, eigenvectors of
   22: *  a complex Hermitian band matrix A.
   23: *
   24: *  Arguments
   25: *  =========
   26: *
   27: *  JOBZ    (input) CHARACTER*1
   28: *          = 'N':  Compute eigenvalues only;
   29: *          = 'V':  Compute eigenvalues and eigenvectors.
   30: *
   31: *  UPLO    (input) CHARACTER*1
   32: *          = 'U':  Upper triangle of A is stored;
   33: *          = 'L':  Lower triangle of A is stored.
   34: *
   35: *  N       (input) INTEGER
   36: *          The order of the matrix A.  N >= 0.
   37: *
   38: *  KD      (input) INTEGER
   39: *          The number of superdiagonals of the matrix A if UPLO = 'U',
   40: *          or the number of subdiagonals if UPLO = 'L'.  KD >= 0.
   41: *
   42: *  AB      (input/output) COMPLEX*16 array, dimension (LDAB, N)
   43: *          On entry, the upper or lower triangle of the Hermitian band
   44: *          matrix A, stored in the first KD+1 rows of the array.  The
   45: *          j-th column of A is stored in the j-th column of the array AB
   46: *          as follows:
   47: *          if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j;
   48: *          if UPLO = 'L', AB(1+i-j,j)    = A(i,j) for j<=i<=min(n,j+kd).
   49: *
   50: *          On exit, AB is overwritten by values generated during the
   51: *          reduction to tridiagonal form.  If UPLO = 'U', the first
   52: *          superdiagonal and the diagonal of the tridiagonal matrix T
   53: *          are returned in rows KD and KD+1 of AB, and if UPLO = 'L',
   54: *          the diagonal and first subdiagonal of T are returned in the
   55: *          first two rows of AB.
   56: *
   57: *  LDAB    (input) INTEGER
   58: *          The leading dimension of the array AB.  LDAB >= KD + 1.
   59: *
   60: *  W       (output) DOUBLE PRECISION array, dimension (N)
   61: *          If INFO = 0, the eigenvalues in ascending order.
   62: *
   63: *  Z       (output) COMPLEX*16 array, dimension (LDZ, N)
   64: *          If JOBZ = 'V', then if INFO = 0, Z contains the orthonormal
   65: *          eigenvectors of the matrix A, with the i-th column of Z
   66: *          holding the eigenvector associated with W(i).
   67: *          If JOBZ = 'N', then Z is not referenced.
   68: *
   69: *  LDZ     (input) INTEGER
   70: *          The leading dimension of the array Z.  LDZ >= 1, and if
   71: *          JOBZ = 'V', LDZ >= max(1,N).
   72: *
   73: *  WORK    (workspace) COMPLEX*16 array, dimension (N)
   74: *
   75: *  RWORK   (workspace) DOUBLE PRECISION array, dimension (max(1,3*N-2))
   76: *
   77: *  INFO    (output) INTEGER
   78: *          = 0:  successful exit.
   79: *          < 0:  if INFO = -i, the i-th argument had an illegal value.
   80: *          > 0:  if INFO = i, the algorithm failed to converge; i
   81: *                off-diagonal elements of an intermediate tridiagonal
   82: *                form did not converge to zero.
   83: *
   84: *  =====================================================================
   85: *
   86: *     .. Parameters ..
   87:       DOUBLE PRECISION   ZERO, ONE
   88:       PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0 )
   89: *     ..
   90: *     .. Local Scalars ..
   91:       LOGICAL            LOWER, WANTZ
   92:       INTEGER            IINFO, IMAX, INDE, INDRWK, ISCALE
   93:       DOUBLE PRECISION   ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA,
   94:      $                   SMLNUM
   95: *     ..
   96: *     .. External Functions ..
   97:       LOGICAL            LSAME
   98:       DOUBLE PRECISION   DLAMCH, ZLANHB
   99:       EXTERNAL           LSAME, DLAMCH, ZLANHB
  100: *     ..
  101: *     .. External Subroutines ..
  102:       EXTERNAL           DSCAL, DSTERF, XERBLA, ZHBTRD, ZLASCL, ZSTEQR
  103: *     ..
  104: *     .. Intrinsic Functions ..
  105:       INTRINSIC          SQRT
  106: *     ..
  107: *     .. Executable Statements ..
  108: *
  109: *     Test the input parameters.
  110: *
  111:       WANTZ = LSAME( JOBZ, 'V' )
  112:       LOWER = LSAME( UPLO, 'L' )
  113: *
  114:       INFO = 0
  115:       IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN
  116:          INFO = -1
  117:       ELSE IF( .NOT.( LOWER .OR. LSAME( UPLO, 'U' ) ) ) THEN
  118:          INFO = -2
  119:       ELSE IF( N.LT.0 ) THEN
  120:          INFO = -3
  121:       ELSE IF( KD.LT.0 ) THEN
  122:          INFO = -4
  123:       ELSE IF( LDAB.LT.KD+1 ) THEN
  124:          INFO = -6
  125:       ELSE IF( LDZ.LT.1 .OR. ( WANTZ .AND. LDZ.LT.N ) ) THEN
  126:          INFO = -9
  127:       END IF
  128: *
  129:       IF( INFO.NE.0 ) THEN
  130:          CALL XERBLA( 'ZHBEV ', -INFO )
  131:          RETURN
  132:       END IF
  133: *
  134: *     Quick return if possible
  135: *
  136:       IF( N.EQ.0 )
  137:      $   RETURN
  138: *
  139:       IF( N.EQ.1 ) THEN
  140:          IF( LOWER ) THEN
  141:             W( 1 ) = AB( 1, 1 )
  142:          ELSE
  143:             W( 1 ) = AB( KD+1, 1 )
  144:          END IF
  145:          IF( WANTZ )
  146:      $      Z( 1, 1 ) = ONE
  147:          RETURN
  148:       END IF
  149: *
  150: *     Get machine constants.
  151: *
  152:       SAFMIN = DLAMCH( 'Safe minimum' )
  153:       EPS = DLAMCH( 'Precision' )
  154:       SMLNUM = SAFMIN / EPS
  155:       BIGNUM = ONE / SMLNUM
  156:       RMIN = SQRT( SMLNUM )
  157:       RMAX = SQRT( BIGNUM )
  158: *
  159: *     Scale matrix to allowable range, if necessary.
  160: *
  161:       ANRM = ZLANHB( 'M', UPLO, N, KD, AB, LDAB, RWORK )
  162:       ISCALE = 0
  163:       IF( ANRM.GT.ZERO .AND. ANRM.LT.RMIN ) THEN
  164:          ISCALE = 1
  165:          SIGMA = RMIN / ANRM
  166:       ELSE IF( ANRM.GT.RMAX ) THEN
  167:          ISCALE = 1
  168:          SIGMA = RMAX / ANRM
  169:       END IF
  170:       IF( ISCALE.EQ.1 ) THEN
  171:          IF( LOWER ) THEN
  172:             CALL ZLASCL( 'B', KD, KD, ONE, SIGMA, N, N, AB, LDAB, INFO )
  173:          ELSE
  174:             CALL ZLASCL( 'Q', KD, KD, ONE, SIGMA, N, N, AB, LDAB, INFO )
  175:          END IF
  176:       END IF
  177: *
  178: *     Call ZHBTRD to reduce Hermitian band matrix to tridiagonal form.
  179: *
  180:       INDE = 1
  181:       CALL ZHBTRD( JOBZ, UPLO, N, KD, AB, LDAB, W, RWORK( INDE ), Z,
  182:      $             LDZ, WORK, IINFO )
  183: *
  184: *     For eigenvalues only, call DSTERF.  For eigenvectors, call ZSTEQR.
  185: *
  186:       IF( .NOT.WANTZ ) THEN
  187:          CALL DSTERF( N, W, RWORK( INDE ), INFO )
  188:       ELSE
  189:          INDRWK = INDE + N
  190:          CALL ZSTEQR( JOBZ, N, W, RWORK( INDE ), Z, LDZ,
  191:      $                RWORK( INDRWK ), INFO )
  192:       END IF
  193: *
  194: *     If matrix was scaled, then rescale eigenvalues appropriately.
  195: *
  196:       IF( ISCALE.EQ.1 ) THEN
  197:          IF( INFO.EQ.0 ) THEN
  198:             IMAX = N
  199:          ELSE
  200:             IMAX = INFO - 1
  201:          END IF
  202:          CALL DSCAL( IMAX, ONE / SIGMA, W, 1 )
  203:       END IF
  204: *
  205:       RETURN
  206: *
  207: *     End of ZHBEV
  208: *
  209:       END