SRC/zhegv.f

   1 *> \brief \b ZHEGV
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download ZHEGV + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zhegv.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zhegv.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zhegv.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE ZHEGV( ITYPE, JOBZ, UPLO, N, A, LDA, B, LDB, W, WORK,
  22 *                         LWORK, RWORK, INFO )
  23 *
  24 *       .. Scalar Arguments ..
  25 *       CHARACTER          JOBZ, UPLO
  26 *       INTEGER            INFO, ITYPE, LDA, LDB, LWORK, N
  27 *       ..
  28 *       .. Array Arguments ..
  29 *       DOUBLE PRECISION   RWORK( * ), W( * )
  30 *       COMPLEX*16         A( LDA, * ), B( LDB, * ), WORK( * )
  31 *       ..
  32 *
  33 *
  34 *> \par Purpose:
  35 *  =============
  36 *>
  37 *> \verbatim
  38 *>
  39 *> ZHEGV computes all the eigenvalues, and optionally, the eigenvectors
  40 *> of a complex generalized Hermitian-definite eigenproblem, of the form
  41 *> A*x=(lambda)*B*x,  A*Bx=(lambda)*x,  or B*A*x=(lambda)*x.
  42 *> Here A and B are assumed to be Hermitian and B is also
  43 *> positive definite.
  44 *> \endverbatim
  45 *
  46 *  Arguments:
  47 *  ==========
  48 *
  49 *> \param[in] ITYPE
  50 *> \verbatim
  51 *>          ITYPE is INTEGER
  52 *>          Specifies the problem type to be solved:
  53 *>          = 1:  A*x = (lambda)*B*x
  54 *>          = 2:  A*B*x = (lambda)*x
  55 *>          = 3:  B*A*x = (lambda)*x
  56 *> \endverbatim
  57 *>
  58 *> \param[in] JOBZ
  59 *> \verbatim
  60 *>          JOBZ is CHARACTER*1
  61 *>          = 'N':  Compute eigenvalues only;
  62 *>          = 'V':  Compute eigenvalues and eigenvectors.
  63 *> \endverbatim
  64 *>
  65 *> \param[in] UPLO
  66 *> \verbatim
  67 *>          UPLO is CHARACTER*1
  68 *>          = 'U':  Upper triangles of A and B are stored;
  69 *>          = 'L':  Lower triangles of A and B are stored.
  70 *> \endverbatim
  71 *>
  72 *> \param[in] N
  73 *> \verbatim
  74 *>          N is INTEGER
  75 *>          The order of the matrices A and B.  N >= 0.
  76 *> \endverbatim
  77 *>
  78 *> \param[in,out] A
  79 *> \verbatim
  80 *>          A is COMPLEX*16 array, dimension (LDA, N)
  81 *>          On entry, the Hermitian matrix A.  If UPLO = 'U', the
  82 *>          leading N-by-N upper triangular part of A contains the
  83 *>          upper triangular part of the matrix A.  If UPLO = 'L',
  84 *>          the leading N-by-N lower triangular part of A contains
  85 *>          the lower triangular part of the matrix A.
  86 *>
  87 *>          On exit, if JOBZ = 'V', then if INFO = 0, A contains the
  88 *>          matrix Z of eigenvectors.  The eigenvectors are normalized
  89 *>          as follows:
  90 *>          if ITYPE = 1 or 2, Z**H*B*Z = I;
  91 *>          if ITYPE = 3, Z**H*inv(B)*Z = I.
  92 *>          If JOBZ = 'N', then on exit the upper triangle (if UPLO='U')
  93 *>          or the lower triangle (if UPLO='L') of A, including the
  94 *>          diagonal, is destroyed.
  95 *> \endverbatim
  96 *>
  97 *> \param[in] LDA
  98 *> \verbatim
  99 *>          LDA is INTEGER
 100 *>          The leading dimension of the array A.  LDA >= max(1,N).
 101 *> \endverbatim
 102 *>
 103 *> \param[in,out] B
 104 *> \verbatim
 105 *>          B is COMPLEX*16 array, dimension (LDB, N)
 106 *>          On entry, the Hermitian positive definite matrix B.
 107 *>          If UPLO = 'U', the leading N-by-N upper triangular part of B
 108 *>          contains the upper triangular part of the matrix B.
 109 *>          If UPLO = 'L', the leading N-by-N lower triangular part of B
 110 *>          contains the lower triangular part of the matrix B.
 111 *>
 112 *>          On exit, if INFO <= N, the part of B containing the matrix is
 113 *>          overwritten by the triangular factor U or L from the Cholesky
 114 *>          factorization B = U**H*U or B = L*L**H.
 115 *> \endverbatim
 116 *>
 117 *> \param[in] LDB
 118 *> \verbatim
 119 *>          LDB is INTEGER
 120 *>          The leading dimension of the array B.  LDB >= max(1,N).
 121 *> \endverbatim
 122 *>
 123 *> \param[out] W
 124 *> \verbatim
 125 *>          W is DOUBLE PRECISION array, dimension (N)
 126 *>          If INFO = 0, the eigenvalues in ascending order.
 127 *> \endverbatim
 128 *>
 129 *> \param[out] WORK
 130 *> \verbatim
 131 *>          WORK is COMPLEX*16 array, dimension (MAX(1,LWORK))
 132 *>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
 133 *> \endverbatim
 134 *>
 135 *> \param[in] LWORK
 136 *> \verbatim
 137 *>          LWORK is INTEGER
 138 *>          The length of the array WORK.  LWORK >= max(1,2*N-1).
 139 *>          For optimal efficiency, LWORK >= (NB+1)*N,
 140 *>          where NB is the blocksize for ZHETRD returned by ILAENV.
 141 *>
 142 *>          If LWORK = -1, then a workspace query is assumed; the routine
 143 *>          only calculates the optimal size of the WORK array, returns
 144 *>          this value as the first entry of the WORK array, and no error
 145 *>          message related to LWORK is issued by XERBLA.
 146 *> \endverbatim
 147 *>
 148 *> \param[out] RWORK
 149 *> \verbatim
 150 *>          RWORK is DOUBLE PRECISION array, dimension (max(1, 3*N-2))
 151 *> \endverbatim
 152 *>
 153 *> \param[out] INFO
 154 *> \verbatim
 155 *>          INFO is INTEGER
 156 *>          = 0:  successful exit
 157 *>          < 0:  if INFO = -i, the i-th argument had an illegal value
 158 *>          > 0:  ZPOTRF or ZHEEV returned an error code:
 159 *>             <= N:  if INFO = i, ZHEEV failed to converge;
 160 *>                    i off-diagonal elements of an intermediate
 161 *>                    tridiagonal form did not converge to zero;
 162 *>             > N:   if INFO = N + i, for 1 <= i <= N, then the leading
 163 *>                    minor of order i of B is not positive definite.
 164 *>                    The factorization of B could not be completed and
 165 *>                    no eigenvalues or eigenvectors were computed.
 166 *> \endverbatim
 167 *
 168 *  Authors:
 169 *  ========
 170 *
 171 *> \author Univ. of Tennessee
 172 *> \author Univ. of California Berkeley
 173 *> \author Univ. of Colorado Denver
 174 *> \author NAG Ltd.
 175 *
 176 *> \date November 2015
 177 *
 178 *> \ingroup complex16HEeigen
 179 *
 180 *  =====================================================================
 181       SUBROUTINE ZHEGV( ITYPE, JOBZ, UPLO, N, A, LDA, B, LDB, W, WORK,
 182      $                  LWORK, RWORK, INFO )
 183 *
 184 *  -- LAPACK driver routine (version 3.6.0) --
 185 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 186 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 187 *     November 2015
 188 *
 189 *     .. Scalar Arguments ..
 190       CHARACTER          JOBZ, UPLO
 191       INTEGER            INFO, ITYPE, LDA, LDB, LWORK, N
 192 *     ..
 193 *     .. Array Arguments ..
 194       DOUBLE PRECISION   RWORK( * ), W( * )
 195       COMPLEX*16         A( LDA, * ), B( LDB, * ), WORK( * )
 196 *     ..
 197 *
 198 *  =====================================================================
 199 *
 200 *     .. Parameters ..
 201       COMPLEX*16         ONE
 202       PARAMETER          ( ONE = ( 1.0D+0, 0.0D+0 ) )
 203 *     ..
 204 *     .. Local Scalars ..
 205       LOGICAL            LQUERY, UPPER, WANTZ
 206       CHARACTER          TRANS
 207       INTEGER            LWKOPT, NB, NEIG
 208 *     ..
 209 *     .. External Functions ..
 210       LOGICAL            LSAME
 211       INTEGER            ILAENV
 212       EXTERNAL           LSAME, ILAENV
 213 *     ..
 214 *     .. External Subroutines ..
 215       EXTERNAL           XERBLA, ZHEEV, ZHEGST, ZPOTRF, ZTRMM, ZTRSM
 216 *     ..
 217 *     .. Intrinsic Functions ..
 218       INTRINSIC          MAX
 219 *     ..
 220 *     .. Executable Statements ..
 221 *
 222 *     Test the input parameters.
 223 *
 224       WANTZ = LSAME( JOBZ, 'V' )
 225       UPPER = LSAME( UPLO, 'U' )
 226       LQUERY = ( LWORK.EQ.-1 )
 227 *
 228       INFO = 0
 229       IF( ITYPE.LT.1 .OR. ITYPE.GT.3 ) THEN
 230          INFO = -1
 231       ELSE IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN
 232          INFO = -2
 233       ELSE IF( .NOT.( UPPER .OR. LSAME( UPLO, 'L' ) ) ) THEN
 234          INFO = -3
 235       ELSE IF( N.LT.0 ) THEN
 236          INFO = -4
 237       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
 238          INFO = -6
 239       ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
 240          INFO = -8
 241       END IF
 242 *
 243       IF( INFO.EQ.0 ) THEN
 244          NB = ILAENV( 1, 'ZHETRD', UPLO, N, -1, -1, -1 )
 245          LWKOPT = MAX( 1, ( NB + 1 )*N )
 246          WORK( 1 ) = LWKOPT
 247 *
 248          IF( LWORK.LT.MAX( 1, 2*N - 1 ) .AND. .NOT.LQUERY ) THEN
 249             INFO = -11
 250          END IF
 251       END IF
 252 *
 253       IF( INFO.NE.0 ) THEN
 254          CALL XERBLA( 'ZHEGV ', -INFO )
 255          RETURN
 256       ELSE IF( LQUERY ) THEN
 257          RETURN
 258       END IF
 259 *
 260 *     Quick return if possible
 261 *
 262       IF( N.EQ.0 )
 263      $   RETURN
 264 *
 265 *     Form a Cholesky factorization of B.
 266 *
 267       CALL ZPOTRF( UPLO, N, B, LDB, INFO )
 268       IF( INFO.NE.0 ) THEN
 269          INFO = N + INFO
 270          RETURN
 271       END IF
 272 *
 273 *     Transform problem to standard eigenvalue problem and solve.
 274 *
 275       CALL ZHEGST( ITYPE, UPLO, N, A, LDA, B, LDB, INFO )
 276       CALL ZHEEV( JOBZ, UPLO, N, A, LDA, W, WORK, LWORK, RWORK, INFO )
 277 *
 278       IF( WANTZ ) THEN
 279 *
 280 *        Backtransform eigenvectors to the original problem.
 281 *
 282          NEIG = N
 283          IF( INFO.GT.0 )
 284      $      NEIG = INFO - 1
 285          IF( ITYPE.EQ.1 .OR. ITYPE.EQ.2 ) THEN
 286 *
 287 *           For A*x=(lambda)*B*x and A*B*x=(lambda)*x;
 288 *           backtransform eigenvectors: x = inv(L)**H *y or inv(U)*y
 289 *
 290             IF( UPPER ) THEN
 291                TRANS = 'N'
 292             ELSE
 293                TRANS = 'C'
 294             END IF
 295 *
 296             CALL ZTRSM( 'Left', UPLO, TRANS, 'Non-unit', N, NEIG, ONE,
 297      $                  B, LDB, A, LDA )
 298 *
 299          ELSE IF( ITYPE.EQ.3 ) THEN
 300 *
 301 *           For B*A*x=(lambda)*x;
 302 *           backtransform eigenvectors: x = L*y or U**H *y
 303 *
 304             IF( UPPER ) THEN
 305                TRANS = 'C'
 306             ELSE
 307                TRANS = 'N'
 308             END IF
 309 *
 310             CALL ZTRMM( 'Left', UPLO, TRANS, 'Non-unit', N, NEIG, ONE,
 311      $                  B, LDB, A, LDA )
 312          END IF
 313       END IF
 314 *
 315       WORK( 1 ) = LWKOPT
 316 *
 317       RETURN
 318 *
 319 *     End of ZHEGV
 320 *
 321       END