LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
zhbevd.f
Go to the documentation of this file.
00001 *> \brief <b> ZHBEVD computes the eigenvalues and, optionally, the left and/or right eigenvectors for OTHER matrices</b>
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *> \htmlonly
00009 *> Download ZHBEVD + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zhbevd.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zhbevd.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zhbevd.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       SUBROUTINE ZHBEVD( JOBZ, UPLO, N, KD, AB, LDAB, W, Z, LDZ, WORK,
00022 *                          LWORK, RWORK, LRWORK, IWORK, LIWORK, INFO )
00023 * 
00024 *       .. Scalar Arguments ..
00025 *       CHARACTER          JOBZ, UPLO
00026 *       INTEGER            INFO, KD, LDAB, LDZ, LIWORK, LRWORK, LWORK, N
00027 *       ..
00028 *       .. Array Arguments ..
00029 *       INTEGER            IWORK( * )
00030 *       DOUBLE PRECISION   RWORK( * ), W( * )
00031 *       COMPLEX*16         AB( LDAB, * ), WORK( * ), Z( LDZ, * )
00032 *       ..
00033 *  
00034 *
00035 *> \par Purpose:
00036 *  =============
00037 *>
00038 *> \verbatim
00039 *>
00040 *> ZHBEVD computes all the eigenvalues and, optionally, eigenvectors of
00041 *> a complex Hermitian band matrix A.  If eigenvectors are desired, it
00042 *> uses a divide and conquer algorithm.
00043 *>
00044 *> The divide and conquer algorithm makes very mild assumptions about
00045 *> floating point arithmetic. It will work on machines with a guard
00046 *> digit in add/subtract, or on those binary machines without guard
00047 *> digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or
00048 *> Cray-2. It could conceivably fail on hexadecimal or decimal machines
00049 *> without guard digits, but we know of none.
00050 *> \endverbatim
00051 *
00052 *  Arguments:
00053 *  ==========
00054 *
00055 *> \param[in] JOBZ
00056 *> \verbatim
00057 *>          JOBZ is CHARACTER*1
00058 *>          = 'N':  Compute eigenvalues only;
00059 *>          = 'V':  Compute eigenvalues and eigenvectors.
00060 *> \endverbatim
00061 *>
00062 *> \param[in] UPLO
00063 *> \verbatim
00064 *>          UPLO is CHARACTER*1
00065 *>          = 'U':  Upper triangle of A is stored;
00066 *>          = 'L':  Lower triangle of A is stored.
00067 *> \endverbatim
00068 *>
00069 *> \param[in] N
00070 *> \verbatim
00071 *>          N is INTEGER
00072 *>          The order of the matrix A.  N >= 0.
00073 *> \endverbatim
00074 *>
00075 *> \param[in] KD
00076 *> \verbatim
00077 *>          KD is INTEGER
00078 *>          The number of superdiagonals of the matrix A if UPLO = 'U',
00079 *>          or the number of subdiagonals if UPLO = 'L'.  KD >= 0.
00080 *> \endverbatim
00081 *>
00082 *> \param[in,out] AB
00083 *> \verbatim
00084 *>          AB is COMPLEX*16 array, dimension (LDAB, N)
00085 *>          On entry, the upper or lower triangle of the Hermitian band
00086 *>          matrix A, stored in the first KD+1 rows of the array.  The
00087 *>          j-th column of A is stored in the j-th column of the array AB
00088 *>          as follows:
00089 *>          if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j;
00090 *>          if UPLO = 'L', AB(1+i-j,j)    = A(i,j) for j<=i<=min(n,j+kd).
00091 *>
00092 *>          On exit, AB is overwritten by values generated during the
00093 *>          reduction to tridiagonal form.  If UPLO = 'U', the first
00094 *>          superdiagonal and the diagonal of the tridiagonal matrix T
00095 *>          are returned in rows KD and KD+1 of AB, and if UPLO = 'L',
00096 *>          the diagonal and first subdiagonal of T are returned in the
00097 *>          first two rows of AB.
00098 *> \endverbatim
00099 *>
00100 *> \param[in] LDAB
00101 *> \verbatim
00102 *>          LDAB is INTEGER
00103 *>          The leading dimension of the array AB.  LDAB >= KD + 1.
00104 *> \endverbatim
00105 *>
00106 *> \param[out] W
00107 *> \verbatim
00108 *>          W is DOUBLE PRECISION array, dimension (N)
00109 *>          If INFO = 0, the eigenvalues in ascending order.
00110 *> \endverbatim
00111 *>
00112 *> \param[out] Z
00113 *> \verbatim
00114 *>          Z is COMPLEX*16 array, dimension (LDZ, N)
00115 *>          If JOBZ = 'V', then if INFO = 0, Z contains the orthonormal
00116 *>          eigenvectors of the matrix A, with the i-th column of Z
00117 *>          holding the eigenvector associated with W(i).
00118 *>          If JOBZ = 'N', then Z is not referenced.
00119 *> \endverbatim
00120 *>
00121 *> \param[in] LDZ
00122 *> \verbatim
00123 *>          LDZ is INTEGER
00124 *>          The leading dimension of the array Z.  LDZ >= 1, and if
00125 *>          JOBZ = 'V', LDZ >= max(1,N).
00126 *> \endverbatim
00127 *>
00128 *> \param[out] WORK
00129 *> \verbatim
00130 *>          WORK is COMPLEX*16 array, dimension (MAX(1,LWORK))
00131 *>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
00132 *> \endverbatim
00133 *>
00134 *> \param[in] LWORK
00135 *> \verbatim
00136 *>          LWORK is INTEGER
00137 *>          The dimension of the array WORK.
00138 *>          If N <= 1,               LWORK must be at least 1.
00139 *>          If JOBZ = 'N' and N > 1, LWORK must be at least N.
00140 *>          If JOBZ = 'V' and N > 1, LWORK must be at least 2*N**2.
00141 *>
00142 *>          If LWORK = -1, then a workspace query is assumed; the routine
00143 *>          only calculates the optimal sizes of the WORK, RWORK and
00144 *>          IWORK arrays, returns these values as the first entries of
00145 *>          the WORK, RWORK and IWORK arrays, and no error message
00146 *>          related to LWORK or LRWORK or LIWORK is issued by XERBLA.
00147 *> \endverbatim
00148 *>
00149 *> \param[out] RWORK
00150 *> \verbatim
00151 *>          RWORK is DOUBLE PRECISION array,
00152 *>                                         dimension (LRWORK)
00153 *>          On exit, if INFO = 0, RWORK(1) returns the optimal LRWORK.
00154 *> \endverbatim
00155 *>
00156 *> \param[in] LRWORK
00157 *> \verbatim
00158 *>          LRWORK is INTEGER
00159 *>          The dimension of array RWORK.
00160 *>          If N <= 1,               LRWORK must be at least 1.
00161 *>          If JOBZ = 'N' and N > 1, LRWORK must be at least N.
00162 *>          If JOBZ = 'V' and N > 1, LRWORK must be at least
00163 *>                        1 + 5*N + 2*N**2.
00164 *>
00165 *>          If LRWORK = -1, then a workspace query is assumed; the
00166 *>          routine only calculates the optimal sizes of the WORK, RWORK
00167 *>          and IWORK arrays, returns these values as the first entries
00168 *>          of the WORK, RWORK and IWORK arrays, and no error message
00169 *>          related to LWORK or LRWORK or LIWORK is issued by XERBLA.
00170 *> \endverbatim
00171 *>
00172 *> \param[out] IWORK
00173 *> \verbatim
00174 *>          IWORK is INTEGER array, dimension (MAX(1,LIWORK))
00175 *>          On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK.
00176 *> \endverbatim
00177 *>
00178 *> \param[in] LIWORK
00179 *> \verbatim
00180 *>          LIWORK is INTEGER
00181 *>          The dimension of array IWORK.
00182 *>          If JOBZ = 'N' or N <= 1, LIWORK must be at least 1.
00183 *>          If JOBZ = 'V' and N > 1, LIWORK must be at least 3 + 5*N .
00184 *>
00185 *>          If LIWORK = -1, then a workspace query is assumed; the
00186 *>          routine only calculates the optimal sizes of the WORK, RWORK
00187 *>          and IWORK arrays, returns these values as the first entries
00188 *>          of the WORK, RWORK and IWORK arrays, and no error message
00189 *>          related to LWORK or LRWORK or LIWORK is issued by XERBLA.
00190 *> \endverbatim
00191 *>
00192 *> \param[out] INFO
00193 *> \verbatim
00194 *>          INFO is INTEGER
00195 *>          = 0:  successful exit.
00196 *>          < 0:  if INFO = -i, the i-th argument had an illegal value.
00197 *>          > 0:  if INFO = i, the algorithm failed to converge; i
00198 *>                off-diagonal elements of an intermediate tridiagonal
00199 *>                form did not converge to zero.
00200 *> \endverbatim
00201 *
00202 *  Authors:
00203 *  ========
00204 *
00205 *> \author Univ. of Tennessee 
00206 *> \author Univ. of California Berkeley 
00207 *> \author Univ. of Colorado Denver 
00208 *> \author NAG Ltd. 
00209 *
00210 *> \date November 2011
00211 *
00212 *> \ingroup complex16OTHEReigen
00213 *
00214 *  =====================================================================
00215       SUBROUTINE ZHBEVD( JOBZ, UPLO, N, KD, AB, LDAB, W, Z, LDZ, WORK,
00216      $                   LWORK, RWORK, LRWORK, IWORK, LIWORK, INFO )
00217 *
00218 *  -- LAPACK driver routine (version 3.4.0) --
00219 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00220 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00221 *     November 2011
00222 *
00223 *     .. Scalar Arguments ..
00224       CHARACTER          JOBZ, UPLO
00225       INTEGER            INFO, KD, LDAB, LDZ, LIWORK, LRWORK, LWORK, N
00226 *     ..
00227 *     .. Array Arguments ..
00228       INTEGER            IWORK( * )
00229       DOUBLE PRECISION   RWORK( * ), W( * )
00230       COMPLEX*16         AB( LDAB, * ), WORK( * ), Z( LDZ, * )
00231 *     ..
00232 *
00233 *  =====================================================================
00234 *
00235 *     .. Parameters ..
00236       DOUBLE PRECISION   ZERO, ONE
00237       PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0 )
00238       COMPLEX*16         CZERO, CONE
00239       PARAMETER          ( CZERO = ( 0.0D0, 0.0D0 ),
00240      $                   CONE = ( 1.0D0, 0.0D0 ) )
00241 *     ..
00242 *     .. Local Scalars ..
00243       LOGICAL            LOWER, LQUERY, WANTZ
00244       INTEGER            IINFO, IMAX, INDE, INDWK2, INDWRK, ISCALE,
00245      $                   LIWMIN, LLRWK, LLWK2, LRWMIN, LWMIN
00246       DOUBLE PRECISION   ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA,
00247      $                   SMLNUM
00248 *     ..
00249 *     .. External Functions ..
00250       LOGICAL            LSAME
00251       DOUBLE PRECISION   DLAMCH, ZLANHB
00252       EXTERNAL           LSAME, DLAMCH, ZLANHB
00253 *     ..
00254 *     .. External Subroutines ..
00255       EXTERNAL           DSCAL, DSTERF, XERBLA, ZGEMM, ZHBTRD, ZLACPY,
00256      $                   ZLASCL, ZSTEDC
00257 *     ..
00258 *     .. Intrinsic Functions ..
00259       INTRINSIC          SQRT
00260 *     ..
00261 *     .. Executable Statements ..
00262 *
00263 *     Test the input parameters.
00264 *
00265       WANTZ = LSAME( JOBZ, 'V' )
00266       LOWER = LSAME( UPLO, 'L' )
00267       LQUERY = ( LWORK.EQ.-1 .OR. LIWORK.EQ.-1 .OR. LRWORK.EQ.-1 )
00268 *
00269       INFO = 0
00270       IF( N.LE.1 ) THEN
00271          LWMIN = 1
00272          LRWMIN = 1
00273          LIWMIN = 1
00274       ELSE
00275          IF( WANTZ ) THEN
00276             LWMIN = 2*N**2
00277             LRWMIN = 1 + 5*N + 2*N**2
00278             LIWMIN = 3 + 5*N
00279          ELSE
00280             LWMIN = N
00281             LRWMIN = N
00282             LIWMIN = 1
00283          END IF
00284       END IF
00285       IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN
00286          INFO = -1
00287       ELSE IF( .NOT.( LOWER .OR. LSAME( UPLO, 'U' ) ) ) THEN
00288          INFO = -2
00289       ELSE IF( N.LT.0 ) THEN
00290          INFO = -3
00291       ELSE IF( KD.LT.0 ) THEN
00292          INFO = -4
00293       ELSE IF( LDAB.LT.KD+1 ) THEN
00294          INFO = -6
00295       ELSE IF( LDZ.LT.1 .OR. ( WANTZ .AND. LDZ.LT.N ) ) THEN
00296          INFO = -9
00297       END IF
00298 *
00299       IF( INFO.EQ.0 ) THEN
00300          WORK( 1 ) = LWMIN
00301          RWORK( 1 ) = LRWMIN
00302          IWORK( 1 ) = LIWMIN
00303 *
00304          IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN
00305             INFO = -11
00306          ELSE IF( LRWORK.LT.LRWMIN .AND. .NOT.LQUERY ) THEN
00307             INFO = -13
00308          ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT.LQUERY ) THEN
00309             INFO = -15
00310          END IF
00311       END IF
00312 *
00313       IF( INFO.NE.0 ) THEN
00314          CALL XERBLA( 'ZHBEVD', -INFO )
00315          RETURN
00316       ELSE IF( LQUERY ) THEN
00317          RETURN
00318       END IF
00319 *
00320 *     Quick return if possible
00321 *
00322       IF( N.EQ.0 )
00323      $   RETURN
00324 *
00325       IF( N.EQ.1 ) THEN
00326          W( 1 ) = AB( 1, 1 )
00327          IF( WANTZ )
00328      $      Z( 1, 1 ) = CONE
00329          RETURN
00330       END IF
00331 *
00332 *     Get machine constants.
00333 *
00334       SAFMIN = DLAMCH( 'Safe minimum' )
00335       EPS = DLAMCH( 'Precision' )
00336       SMLNUM = SAFMIN / EPS
00337       BIGNUM = ONE / SMLNUM
00338       RMIN = SQRT( SMLNUM )
00339       RMAX = SQRT( BIGNUM )
00340 *
00341 *     Scale matrix to allowable range, if necessary.
00342 *
00343       ANRM = ZLANHB( 'M', UPLO, N, KD, AB, LDAB, RWORK )
00344       ISCALE = 0
00345       IF( ANRM.GT.ZERO .AND. ANRM.LT.RMIN ) THEN
00346          ISCALE = 1
00347          SIGMA = RMIN / ANRM
00348       ELSE IF( ANRM.GT.RMAX ) THEN
00349          ISCALE = 1
00350          SIGMA = RMAX / ANRM
00351       END IF
00352       IF( ISCALE.EQ.1 ) THEN
00353          IF( LOWER ) THEN
00354             CALL ZLASCL( 'B', KD, KD, ONE, SIGMA, N, N, AB, LDAB, INFO )
00355          ELSE
00356             CALL ZLASCL( 'Q', KD, KD, ONE, SIGMA, N, N, AB, LDAB, INFO )
00357          END IF
00358       END IF
00359 *
00360 *     Call ZHBTRD to reduce Hermitian band matrix to tridiagonal form.
00361 *
00362       INDE = 1
00363       INDWRK = INDE + N
00364       INDWK2 = 1 + N*N
00365       LLWK2 = LWORK - INDWK2 + 1
00366       LLRWK = LRWORK - INDWRK + 1
00367       CALL ZHBTRD( JOBZ, UPLO, N, KD, AB, LDAB, W, RWORK( INDE ), Z,
00368      $             LDZ, WORK, IINFO )
00369 *
00370 *     For eigenvalues only, call DSTERF.  For eigenvectors, call ZSTEDC.
00371 *
00372       IF( .NOT.WANTZ ) THEN
00373          CALL DSTERF( N, W, RWORK( INDE ), INFO )
00374       ELSE
00375          CALL ZSTEDC( 'I', N, W, RWORK( INDE ), WORK, N, WORK( INDWK2 ),
00376      $                LLWK2, RWORK( INDWRK ), LLRWK, IWORK, LIWORK,
00377      $                INFO )
00378          CALL ZGEMM( 'N', 'N', N, N, N, CONE, Z, LDZ, WORK, N, CZERO,
00379      $               WORK( INDWK2 ), N )
00380          CALL ZLACPY( 'A', N, N, WORK( INDWK2 ), N, Z, LDZ )
00381       END IF
00382 *
00383 *     If matrix was scaled, then rescale eigenvalues appropriately.
00384 *
00385       IF( ISCALE.EQ.1 ) THEN
00386          IF( INFO.EQ.0 ) THEN
00387             IMAX = N
00388          ELSE
00389             IMAX = INFO - 1
00390          END IF
00391          CALL DSCAL( IMAX, ONE / SIGMA, W, 1 )
00392       END IF
00393 *
00394       WORK( 1 ) = LWMIN
00395       RWORK( 1 ) = LRWMIN
00396       IWORK( 1 ) = LIWMIN
00397       RETURN
00398 *
00399 *     End of ZHBEVD
00400 *
00401       END
 All Files Functions