LAPACK: SRC/zhbev.f Source File

Go to the documentation of this file.
00001       SUBROUTINE ZHBEV( JOBZ, UPLO, N, KD, AB, LDAB, W, Z, LDZ, WORK,
00002      $                  RWORK, INFO )
00003 *
00004 *  -- LAPACK driver routine (version 3.2) --
00005 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00006 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00007 *     November 2006
00008 *
00009 *     .. Scalar Arguments ..
00010       CHARACTER          JOBZ, UPLO
00011       INTEGER            INFO, KD, LDAB, LDZ, N
00012 *     ..
00013 *     .. Array Arguments ..
00014       DOUBLE PRECISION   RWORK( * ), W( * )
00015       COMPLEX*16         AB( LDAB, * ), WORK( * ), Z( LDZ, * )
00016 *     ..
00017 *
00018 *  Purpose
00019 *  =======
00020 *
00021 *  ZHBEV computes all the eigenvalues and, optionally, eigenvectors of
00022 *  a complex Hermitian band matrix A.
00023 *
00024 *  Arguments
00025 *  =========
00026 *
00027 *  JOBZ    (input) CHARACTER*1
00028 *          = 'N':  Compute eigenvalues only;
00029 *          = 'V':  Compute eigenvalues and eigenvectors.
00030 *
00031 *  UPLO    (input) CHARACTER*1
00032 *          = 'U':  Upper triangle of A is stored;
00033 *          = 'L':  Lower triangle of A is stored.
00034 *
00035 *  N       (input) INTEGER
00036 *          The order of the matrix A.  N >= 0.
00037 *
00038 *  KD      (input) INTEGER
00039 *          The number of superdiagonals of the matrix A if UPLO = 'U',
00040 *          or the number of subdiagonals if UPLO = 'L'.  KD >= 0.
00041 *
00042 *  AB      (input/output) COMPLEX*16 array, dimension (LDAB, N)
00043 *          On entry, the upper or lower triangle of the Hermitian band
00044 *          matrix A, stored in the first KD+1 rows of the array.  The
00045 *          j-th column of A is stored in the j-th column of the array AB
00046 *          as follows:
00047 *          if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j;
00048 *          if UPLO = 'L', AB(1+i-j,j)    = A(i,j) for j<=i<=min(n,j+kd).
00049 *
00050 *          On exit, AB is overwritten by values generated during the
00051 *          reduction to tridiagonal form.  If UPLO = 'U', the first
00052 *          superdiagonal and the diagonal of the tridiagonal matrix T
00053 *          are returned in rows KD and KD+1 of AB, and if UPLO = 'L',
00054 *          the diagonal and first subdiagonal of T are returned in the
00055 *          first two rows of AB.
00056 *
00057 *  LDAB    (input) INTEGER
00058 *          The leading dimension of the array AB.  LDAB >= KD + 1.
00059 *
00060 *  W       (output) DOUBLE PRECISION array, dimension (N)
00061 *          If INFO = 0, the eigenvalues in ascending order.
00062 *
00063 *  Z       (output) COMPLEX*16 array, dimension (LDZ, N)
00064 *          If JOBZ = 'V', then if INFO = 0, Z contains the orthonormal
00065 *          eigenvectors of the matrix A, with the i-th column of Z
00066 *          holding the eigenvector associated with W(i).
00067 *          If JOBZ = 'N', then Z is not referenced.
00068 *
00069 *  LDZ     (input) INTEGER
00070 *          The leading dimension of the array Z.  LDZ >= 1, and if
00071 *          JOBZ = 'V', LDZ >= max(1,N).
00072 *
00073 *  WORK    (workspace) COMPLEX*16 array, dimension (N)
00074 *
00075 *  RWORK   (workspace) DOUBLE PRECISION array, dimension (max(1,3*N-2))
00076 *
00077 *  INFO    (output) INTEGER
00078 *          = 0:  successful exit.
00079 *          < 0:  if INFO = -i, the i-th argument had an illegal value.
00080 *          > 0:  if INFO = i, the algorithm failed to converge; i
00081 *                off-diagonal elements of an intermediate tridiagonal
00082 *                form did not converge to zero.
00083 *
00084 *  =====================================================================
00085 *
00086 *     .. Parameters ..
00087       DOUBLE PRECISION   ZERO, ONE
00088       PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0 )
00089 *     ..
00090 *     .. Local Scalars ..
00091       LOGICAL            LOWER, WANTZ
00092       INTEGER            IINFO, IMAX, INDE, INDRWK, ISCALE
00093       DOUBLE PRECISION   ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA,
00094      $                   SMLNUM
00095 *     ..
00096 *     .. External Functions ..
00097       LOGICAL            LSAME
00098       DOUBLE PRECISION   DLAMCH, ZLANHB
00099       EXTERNAL           LSAME, DLAMCH, ZLANHB
00100 *     ..
00101 *     .. External Subroutines ..
00102       EXTERNAL           DSCAL, DSTERF, XERBLA, ZHBTRD, ZLASCL, ZSTEQR
00103 *     ..
00104 *     .. Intrinsic Functions ..
00105       INTRINSIC          SQRT
00106 *     ..
00107 *     .. Executable Statements ..
00108 *
00109 *     Test the input parameters.
00110 *
00111       WANTZ = LSAME( JOBZ, 'V' )
00112       LOWER = LSAME( UPLO, 'L' )
00113 *
00114       INFO = 0
00115       IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN
00116          INFO = -1
00117       ELSE IF( .NOT.( LOWER .OR. LSAME( UPLO, 'U' ) ) ) THEN
00118          INFO = -2
00119       ELSE IF( N.LT.0 ) THEN
00120          INFO = -3
00121       ELSE IF( KD.LT.0 ) THEN
00122          INFO = -4
00123       ELSE IF( LDAB.LT.KD+1 ) THEN
00124          INFO = -6
00125       ELSE IF( LDZ.LT.1 .OR. ( WANTZ .AND. LDZ.LT.N ) ) THEN
00126          INFO = -9
00127       END IF
00128 *
00129       IF( INFO.NE.0 ) THEN
00130          CALL XERBLA( 'ZHBEV ', -INFO )
00131          RETURN
00132       END IF
00133 *
00134 *     Quick return if possible
00135 *
00136       IF( N.EQ.0 )
00137      $   RETURN
00138 *
00139       IF( N.EQ.1 ) THEN
00140          IF( LOWER ) THEN
00141             W( 1 ) = AB( 1, 1 )
00142          ELSE
00143             W( 1 ) = AB( KD+1, 1 )
00144          END IF
00145          IF( WANTZ )
00146      $      Z( 1, 1 ) = ONE
00147          RETURN
00148       END IF
00149 *
00150 *     Get machine constants.
00151 *
00152       SAFMIN = DLAMCH( 'Safe minimum' )
00153       EPS = DLAMCH( 'Precision' )
00154       SMLNUM = SAFMIN / EPS
00155       BIGNUM = ONE / SMLNUM
00156       RMIN = SQRT( SMLNUM )
00157       RMAX = SQRT( BIGNUM )
00158 *
00159 *     Scale matrix to allowable range, if necessary.
00160 *
00161       ANRM = ZLANHB( 'M', UPLO, N, KD, AB, LDAB, RWORK )
00162       ISCALE = 0
00163       IF( ANRM.GT.ZERO .AND. ANRM.LT.RMIN ) THEN
00164          ISCALE = 1
00165          SIGMA = RMIN / ANRM
00166       ELSE IF( ANRM.GT.RMAX ) THEN
00167          ISCALE = 1
00168          SIGMA = RMAX / ANRM
00169       END IF
00170       IF( ISCALE.EQ.1 ) THEN
00171          IF( LOWER ) THEN
00172             CALL ZLASCL( 'B', KD, KD, ONE, SIGMA, N, N, AB, LDAB, INFO )
00173          ELSE
00174             CALL ZLASCL( 'Q', KD, KD, ONE, SIGMA, N, N, AB, LDAB, INFO )
00175          END IF
00176       END IF
00177 *
00178 *     Call ZHBTRD to reduce Hermitian band matrix to tridiagonal form.
00179 *
00180       INDE = 1
00181       CALL ZHBTRD( JOBZ, UPLO, N, KD, AB, LDAB, W, RWORK( INDE ), Z,
00182      $             LDZ, WORK, IINFO )
00183 *
00184 *     For eigenvalues only, call DSTERF.  For eigenvectors, call ZSTEQR.
00185 *
00186       IF( .NOT.WANTZ ) THEN
00187          CALL DSTERF( N, W, RWORK( INDE ), INFO )
00188       ELSE
00189          INDRWK = INDE + N
00190          CALL ZSTEQR( JOBZ, N, W, RWORK( INDE ), Z, LDZ,
00191      $                RWORK( INDRWK ), INFO )
00192       END IF
00193 *
00194 *     If matrix was scaled, then rescale eigenvalues appropriately.
00195 *
00196       IF( ISCALE.EQ.1 ) THEN
00197          IF( INFO.EQ.0 ) THEN
00198             IMAX = N
00199          ELSE
00200             IMAX = INFO - 1
00201          END IF
00202          CALL DSCAL( IMAX, ONE / SIGMA, W, 1 )
00203       END IF
00204 *
00205       RETURN
00206 *
00207 *     End of ZHBEV
00208 *
00209       END