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