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