SRC/zhegvd.f

   1 *> \brief \b ZHEGVD
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download ZHEGVD + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zhegvd.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zhegvd.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zhegvd.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE ZHEGVD( ITYPE, JOBZ, UPLO, N, A, LDA, B, LDB, W, WORK,
  22 *                          LWORK, RWORK, LRWORK, IWORK, LIWORK, INFO )
  23 *
  24 *       .. Scalar Arguments ..
  25 *       CHARACTER          JOBZ, UPLO
  26 *       INTEGER            INFO, ITYPE, LDA, LDB, LIWORK, LRWORK, LWORK, N
  27 *       ..
  28 *       .. Array Arguments ..
  29 *       INTEGER            IWORK( * )
  30 *       DOUBLE PRECISION   RWORK( * ), W( * )
  31 *       COMPLEX*16         A( LDA, * ), B( LDB, * ), WORK( * )
  32 *       ..
  33 *
  34 *
  35 *> \par Purpose:
  36 *  =============
  37 *>
  38 *> \verbatim
  39 *>
  40 *> ZHEGVD computes all the eigenvalues, and optionally, the eigenvectors
  41 *> of a complex generalized Hermitian-definite eigenproblem, of the form
  42 *> A*x=(lambda)*B*x,  A*Bx=(lambda)*x,  or B*A*x=(lambda)*x.  Here A and
  43 *> B are assumed to be Hermitian and B is also positive definite.
  44 *> If eigenvectors are desired, it uses a divide and conquer algorithm.
  45 *>
  46 *> The divide and conquer algorithm makes very mild assumptions about
  47 *> floating point arithmetic. It will work on machines with a guard
  48 *> digit in add/subtract, or on those binary machines without guard
  49 *> digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or
  50 *> Cray-2. It could conceivably fail on hexadecimal or decimal machines
  51 *> without guard digits, but we know of none.
  52 *> \endverbatim
  53 *
  54 *  Arguments:
  55 *  ==========
  56 *
  57 *> \param[in] ITYPE
  58 *> \verbatim
  59 *>          ITYPE is INTEGER
  60 *>          Specifies the problem type to be solved:
  61 *>          = 1:  A*x = (lambda)*B*x
  62 *>          = 2:  A*B*x = (lambda)*x
  63 *>          = 3:  B*A*x = (lambda)*x
  64 *> \endverbatim
  65 *>
  66 *> \param[in] JOBZ
  67 *> \verbatim
  68 *>          JOBZ is CHARACTER*1
  69 *>          = 'N':  Compute eigenvalues only;
  70 *>          = 'V':  Compute eigenvalues and eigenvectors.
  71 *> \endverbatim
  72 *>
  73 *> \param[in] UPLO
  74 *> \verbatim
  75 *>          UPLO is CHARACTER*1
  76 *>          = 'U':  Upper triangles of A and B are stored;
  77 *>          = 'L':  Lower triangles of A and B are stored.
  78 *> \endverbatim
  79 *>
  80 *> \param[in] N
  81 *> \verbatim
  82 *>          N is INTEGER
  83 *>          The order of the matrices A and B.  N >= 0.
  84 *> \endverbatim
  85 *>
  86 *> \param[in,out] A
  87 *> \verbatim
  88 *>          A is COMPLEX*16 array, dimension (LDA, N)
  89 *>          On entry, the Hermitian matrix A.  If UPLO = 'U', the
  90 *>          leading N-by-N upper triangular part of A contains the
  91 *>          upper triangular part of the matrix A.  If UPLO = 'L',
  92 *>          the leading N-by-N lower triangular part of A contains
  93 *>          the lower triangular part of the matrix A.
  94 *>
  95 *>          On exit, if JOBZ = 'V', then if INFO = 0, A contains the
  96 *>          matrix Z of eigenvectors.  The eigenvectors are normalized
  97 *>          as follows:
  98 *>          if ITYPE = 1 or 2, Z**H*B*Z = I;
  99 *>          if ITYPE = 3, Z**H*inv(B)*Z = I.
 100 *>          If JOBZ = 'N', then on exit the upper triangle (if UPLO='U')
 101 *>          or the lower triangle (if UPLO='L') of A, including the
 102 *>          diagonal, is destroyed.
 103 *> \endverbatim
 104 *>
 105 *> \param[in] LDA
 106 *> \verbatim
 107 *>          LDA is INTEGER
 108 *>          The leading dimension of the array A.  LDA >= max(1,N).
 109 *> \endverbatim
 110 *>
 111 *> \param[in,out] B
 112 *> \verbatim
 113 *>          B is COMPLEX*16 array, dimension (LDB, N)
 114 *>          On entry, the Hermitian matrix B.  If UPLO = 'U', the
 115 *>          leading N-by-N upper triangular part of B contains the
 116 *>          upper triangular part of the matrix B.  If UPLO = 'L',
 117 *>          the leading N-by-N lower triangular part of B contains
 118 *>          the lower triangular part of the matrix B.
 119 *>
 120 *>          On exit, if INFO <= N, the part of B containing the matrix is
 121 *>          overwritten by the triangular factor U or L from the Cholesky
 122 *>          factorization B = U**H*U or B = L*L**H.
 123 *> \endverbatim
 124 *>
 125 *> \param[in] LDB
 126 *> \verbatim
 127 *>          LDB is INTEGER
 128 *>          The leading dimension of the array B.  LDB >= max(1,N).
 129 *> \endverbatim
 130 *>
 131 *> \param[out] W
 132 *> \verbatim
 133 *>          W is DOUBLE PRECISION array, dimension (N)
 134 *>          If INFO = 0, the eigenvalues in ascending order.
 135 *> \endverbatim
 136 *>
 137 *> \param[out] WORK
 138 *> \verbatim
 139 *>          WORK is COMPLEX*16 array, dimension (MAX(1,LWORK))
 140 *>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
 141 *> \endverbatim
 142 *>
 143 *> \param[in] LWORK
 144 *> \verbatim
 145 *>          LWORK is INTEGER
 146 *>          The length of the array WORK.
 147 *>          If N <= 1,                LWORK >= 1.
 148 *>          If JOBZ  = 'N' and N > 1, LWORK >= N + 1.
 149 *>          If JOBZ  = 'V' and N > 1, LWORK >= 2*N + N**2.
 150 *>
 151 *>          If LWORK = -1, then a workspace query is assumed; the routine
 152 *>          only calculates the optimal sizes of the WORK, RWORK and
 153 *>          IWORK arrays, returns these values as the first entries of
 154 *>          the WORK, RWORK and IWORK arrays, and no error message
 155 *>          related to LWORK or LRWORK or LIWORK is issued by XERBLA.
 156 *> \endverbatim
 157 *>
 158 *> \param[out] RWORK
 159 *> \verbatim
 160 *>          RWORK is DOUBLE PRECISION array, dimension (MAX(1,LRWORK))
 161 *>          On exit, if INFO = 0, RWORK(1) returns the optimal LRWORK.
 162 *> \endverbatim
 163 *>
 164 *> \param[in] LRWORK
 165 *> \verbatim
 166 *>          LRWORK is INTEGER
 167 *>          The dimension of the array RWORK.
 168 *>          If N <= 1,                LRWORK >= 1.
 169 *>          If JOBZ  = 'N' and N > 1, LRWORK >= N.
 170 *>          If JOBZ  = 'V' and N > 1, LRWORK >= 1 + 5*N + 2*N**2.
 171 *>
 172 *>          If LRWORK = -1, then a workspace query is assumed; the
 173 *>          routine only calculates the optimal sizes of the WORK, RWORK
 174 *>          and IWORK arrays, returns these values as the first entries
 175 *>          of the WORK, RWORK and IWORK arrays, and no error message
 176 *>          related to LWORK or LRWORK or LIWORK is issued by XERBLA.
 177 *> \endverbatim
 178 *>
 179 *> \param[out] IWORK
 180 *> \verbatim
 181 *>          IWORK is INTEGER array, dimension (MAX(1,LIWORK))
 182 *>          On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK.
 183 *> \endverbatim
 184 *>
 185 *> \param[in] LIWORK
 186 *> \verbatim
 187 *>          LIWORK is INTEGER
 188 *>          The dimension of the array IWORK.
 189 *>          If N <= 1,                LIWORK >= 1.
 190 *>          If JOBZ  = 'N' and N > 1, LIWORK >= 1.
 191 *>          If JOBZ  = 'V' and N > 1, LIWORK >= 3 + 5*N.
 192 *>
 193 *>          If LIWORK = -1, then a workspace query is assumed; the
 194 *>          routine only calculates the optimal sizes of the WORK, RWORK
 195 *>          and IWORK arrays, returns these values as the first entries
 196 *>          of the WORK, RWORK and IWORK arrays, and no error message
 197 *>          related to LWORK or LRWORK or LIWORK is issued by XERBLA.
 198 *> \endverbatim
 199 *>
 200 *> \param[out] INFO
 201 *> \verbatim
 202 *>          INFO is INTEGER
 203 *>          = 0:  successful exit
 204 *>          < 0:  if INFO = -i, the i-th argument had an illegal value
 205 *>          > 0:  ZPOTRF or ZHEEVD returned an error code:
 206 *>             <= N:  if INFO = i and JOBZ = 'N', then the algorithm
 207 *>                    failed to converge; i off-diagonal elements of an
 208 *>                    intermediate tridiagonal form did not converge to
 209 *>                    zero;
 210 *>                    if INFO = i and JOBZ = 'V', then the algorithm
 211 *>                    failed to compute an eigenvalue while working on
 212 *>                    the submatrix lying in rows and columns INFO/(N+1)
 213 *>                    through mod(INFO,N+1);
 214 *>             > N:   if INFO = N + i, for 1 <= i <= N, then the leading
 215 *>                    minor of order i of B is not positive definite.
 216 *>                    The factorization of B could not be completed and
 217 *>                    no eigenvalues or eigenvectors were computed.
 218 *> \endverbatim
 219 *
 220 *  Authors:
 221 *  ========
 222 *
 223 *> \author Univ. of Tennessee
 224 *> \author Univ. of California Berkeley
 225 *> \author Univ. of Colorado Denver
 226 *> \author NAG Ltd.
 227 *
 228 *> \date November 2015
 229 *
 230 *> \ingroup complex16HEeigen
 231 *
 232 *> \par Further Details:
 233 *  =====================
 234 *>
 235 *> \verbatim
 236 *>
 237 *>  Modified so that no backsubstitution is performed if ZHEEVD fails to
 238 *>  converge (NEIG in old code could be greater than N causing out of
 239 *>  bounds reference to A - reported by Ralf Meyer).  Also corrected the
 240 *>  description of INFO and the test on ITYPE. Sven, 16 Feb 05.
 241 *> \endverbatim
 242 *
 243 *> \par Contributors:
 244 *  ==================
 245 *>
 246 *>     Mark Fahey, Department of Mathematics, Univ. of Kentucky, USA
 247 *>
 248 *  =====================================================================
 249       SUBROUTINE ZHEGVD( ITYPE, JOBZ, UPLO, N, A, LDA, B, LDB, W, WORK,
 250      $                   LWORK, RWORK, LRWORK, IWORK, LIWORK, INFO )
 251 *
 252 *  -- LAPACK driver routine (version 3.6.0) --
 253 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 254 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 255 *     November 2015
 256 *
 257 *     .. Scalar Arguments ..
 258       CHARACTER          JOBZ, UPLO
 259       INTEGER            INFO, ITYPE, LDA, LDB, LIWORK, LRWORK, LWORK, N
 260 *     ..
 261 *     .. Array Arguments ..
 262       INTEGER            IWORK( * )
 263       DOUBLE PRECISION   RWORK( * ), W( * )
 264       COMPLEX*16         A( LDA, * ), B( LDB, * ), WORK( * )
 265 *     ..
 266 *
 267 *  =====================================================================
 268 *
 269 *     .. Parameters ..
 270       COMPLEX*16         CONE
 271       PARAMETER          ( CONE = ( 1.0D+0, 0.0D+0 ) )
 272 *     ..
 273 *     .. Local Scalars ..
 274       LOGICAL            LQUERY, UPPER, WANTZ
 275       CHARACTER          TRANS
 276       INTEGER            LIOPT, LIWMIN, LOPT, LROPT, LRWMIN, LWMIN
 277 *     ..
 278 *     .. External Functions ..
 279       LOGICAL            LSAME
 280       EXTERNAL           LSAME
 281 *     ..
 282 *     .. External Subroutines ..
 283       EXTERNAL           XERBLA, ZHEEVD, ZHEGST, ZPOTRF, ZTRMM, ZTRSM
 284 *     ..
 285 *     .. Intrinsic Functions ..
 286       INTRINSIC          DBLE, MAX
 287 *     ..
 288 *     .. Executable Statements ..
 289 *
 290 *     Test the input parameters.
 291 *
 292       WANTZ = LSAME( JOBZ, 'V' )
 293       UPPER = LSAME( UPLO, 'U' )
 294       LQUERY = ( LWORK.EQ.-1 .OR. LRWORK.EQ.-1 .OR. LIWORK.EQ.-1 )
 295 *
 296       INFO = 0
 297       IF( N.LE.1 ) THEN
 298          LWMIN = 1
 299          LRWMIN = 1
 300          LIWMIN = 1
 301       ELSE IF( WANTZ ) THEN
 302          LWMIN = 2*N + N*N
 303          LRWMIN = 1 + 5*N + 2*N*N
 304          LIWMIN = 3 + 5*N
 305       ELSE
 306          LWMIN = N + 1
 307          LRWMIN = N
 308          LIWMIN = 1
 309       END IF
 310       LOPT = LWMIN
 311       LROPT = LRWMIN
 312       LIOPT = LIWMIN
 313       IF( ITYPE.LT.1 .OR. ITYPE.GT.3 ) THEN
 314          INFO = -1
 315       ELSE IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN
 316          INFO = -2
 317       ELSE IF( .NOT.( UPPER .OR. LSAME( UPLO, 'L' ) ) ) THEN
 318          INFO = -3
 319       ELSE IF( N.LT.0 ) THEN
 320          INFO = -4
 321       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
 322          INFO = -6
 323       ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
 324          INFO = -8
 325       END IF
 326 *
 327       IF( INFO.EQ.0 ) THEN
 328          WORK( 1 ) = LOPT
 329          RWORK( 1 ) = LROPT
 330          IWORK( 1 ) = LIOPT
 331 *
 332          IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN
 333             INFO = -11
 334          ELSE IF( LRWORK.LT.LRWMIN .AND. .NOT.LQUERY ) THEN
 335             INFO = -13
 336          ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT.LQUERY ) THEN
 337             INFO = -15
 338          END IF
 339       END IF
 340 *
 341       IF( INFO.NE.0 ) THEN
 342          CALL XERBLA( 'ZHEGVD', -INFO )
 343          RETURN
 344       ELSE IF( LQUERY ) THEN
 345          RETURN
 346       END IF
 347 *
 348 *     Quick return if possible
 349 *
 350       IF( N.EQ.0 )
 351      $   RETURN
 352 *
 353 *     Form a Cholesky factorization of B.
 354 *
 355       CALL ZPOTRF( UPLO, N, B, LDB, INFO )
 356       IF( INFO.NE.0 ) THEN
 357          INFO = N + INFO
 358          RETURN
 359       END IF
 360 *
 361 *     Transform problem to standard eigenvalue problem and solve.
 362 *
 363       CALL ZHEGST( ITYPE, UPLO, N, A, LDA, B, LDB, INFO )
 364       CALL ZHEEVD( JOBZ, UPLO, N, A, LDA, W, WORK, LWORK, RWORK, LRWORK,
 365      $             IWORK, LIWORK, INFO )
 366       LOPT = MAX( DBLE( LOPT ), DBLE( WORK( 1 ) ) )
 367       LROPT = MAX( DBLE( LROPT ), DBLE( RWORK( 1 ) ) )
 368       LIOPT = MAX( DBLE( LIOPT ), DBLE( IWORK( 1 ) ) )
 369 *
 370       IF( WANTZ .AND. INFO.EQ.0 ) THEN
 371 *
 372 *        Backtransform eigenvectors to the original problem.
 373 *
 374          IF( ITYPE.EQ.1 .OR. ITYPE.EQ.2 ) THEN
 375 *
 376 *           For A*x=(lambda)*B*x and A*B*x=(lambda)*x;
 377 *           backtransform eigenvectors: x = inv(L)**H *y or inv(U)*y
 378 *
 379             IF( UPPER ) THEN
 380                TRANS = 'N'
 381             ELSE
 382                TRANS = 'C'
 383             END IF
 384 *
 385             CALL ZTRSM( 'Left', UPLO, TRANS, 'Non-unit', N, N, CONE,
 386      $                  B, LDB, A, LDA )
 387 *
 388          ELSE IF( ITYPE.EQ.3 ) THEN
 389 *
 390 *           For B*A*x=(lambda)*x;
 391 *           backtransform eigenvectors: x = L*y or U**H *y
 392 *
 393             IF( UPPER ) THEN
 394                TRANS = 'C'
 395             ELSE
 396                TRANS = 'N'
 397             END IF
 398 *
 399             CALL ZTRMM( 'Left', UPLO, TRANS, 'Non-unit', N, N, CONE,
 400      $                  B, LDB, A, LDA )
 401          END IF
 402       END IF
 403 *
 404       WORK( 1 ) = LOPT
 405       RWORK( 1 ) = LROPT
 406       IWORK( 1 ) = LIOPT
 407 *
 408       RETURN
 409 *
 410 *     End of ZHEGVD
 411 *
 412       END