SRC/zhpgvx.f

   1 *> \brief \b ZHPGVX
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download ZHPGVX + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zhpgvx.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zhpgvx.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zhpgvx.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE ZHPGVX( ITYPE, JOBZ, RANGE, UPLO, N, AP, BP, VL, VU,
  22 *                          IL, IU, ABSTOL, M, W, Z, LDZ, WORK, RWORK,
  23 *                          IWORK, IFAIL, INFO )
  24 *
  25 *       .. Scalar Arguments ..
  26 *       CHARACTER          JOBZ, RANGE, UPLO
  27 *       INTEGER            IL, INFO, ITYPE, IU, LDZ, M, N
  28 *       DOUBLE PRECISION   ABSTOL, VL, VU
  29 *       ..
  30 *       .. Array Arguments ..
  31 *       INTEGER            IFAIL( * ), IWORK( * )
  32 *       DOUBLE PRECISION   RWORK( * ), W( * )
  33 *       COMPLEX*16         AP( * ), BP( * ), WORK( * ), Z( LDZ, * )
  34 *       ..
  35 *
  36 *
  37 *> \par Purpose:
  38 *  =============
  39 *>
  40 *> \verbatim
  41 *>
  42 *> ZHPGVX computes selected eigenvalues and, optionally, eigenvectors
  43 *> of a complex generalized Hermitian-definite eigenproblem, of the form
  44 *> A*x=(lambda)*B*x,  A*Bx=(lambda)*x,  or B*A*x=(lambda)*x.  Here A and
  45 *> B are assumed to be Hermitian, stored in packed format, and B is also
  46 *> positive definite.  Eigenvalues and eigenvectors can be selected by
  47 *> specifying either a range of values or a range of indices for the
  48 *> desired eigenvalues.
  49 *> \endverbatim
  50 *
  51 *  Arguments:
  52 *  ==========
  53 *
  54 *> \param[in] ITYPE
  55 *> \verbatim
  56 *>          ITYPE is INTEGER
  57 *>          Specifies the problem type to be solved:
  58 *>          = 1:  A*x = (lambda)*B*x
  59 *>          = 2:  A*B*x = (lambda)*x
  60 *>          = 3:  B*A*x = (lambda)*x
  61 *> \endverbatim
  62 *>
  63 *> \param[in] JOBZ
  64 *> \verbatim
  65 *>          JOBZ is CHARACTER*1
  66 *>          = 'N':  Compute eigenvalues only;
  67 *>          = 'V':  Compute eigenvalues and eigenvectors.
  68 *> \endverbatim
  69 *>
  70 *> \param[in] RANGE
  71 *> \verbatim
  72 *>          RANGE is CHARACTER*1
  73 *>          = 'A': all eigenvalues will be found;
  74 *>          = 'V': all eigenvalues in the half-open interval (VL,VU]
  75 *>                 will be found;
  76 *>          = 'I': the IL-th through IU-th eigenvalues will be found.
  77 *> \endverbatim
  78 *>
  79 *> \param[in] UPLO
  80 *> \verbatim
  81 *>          UPLO is CHARACTER*1
  82 *>          = 'U':  Upper triangles of A and B are stored;
  83 *>          = 'L':  Lower triangles of A and B are stored.
  84 *> \endverbatim
  85 *>
  86 *> \param[in] N
  87 *> \verbatim
  88 *>          N is INTEGER
  89 *>          The order of the matrices A and B.  N >= 0.
  90 *> \endverbatim
  91 *>
  92 *> \param[in,out] AP
  93 *> \verbatim
  94 *>          AP is COMPLEX*16 array, dimension (N*(N+1)/2)
  95 *>          On entry, the upper or lower triangle of the Hermitian matrix
  96 *>          A, packed columnwise in a linear array.  The j-th column of A
  97 *>          is stored in the array AP as follows:
  98 *>          if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
  99 *>          if UPLO = 'L', AP(i + (j-1)*(2*n-j)/2) = A(i,j) for j<=i<=n.
 100 *>
 101 *>          On exit, the contents of AP are destroyed.
 102 *> \endverbatim
 103 *>
 104 *> \param[in,out] BP
 105 *> \verbatim
 106 *>          BP is COMPLEX*16 array, dimension (N*(N+1)/2)
 107 *>          On entry, the upper or lower triangle of the Hermitian matrix
 108 *>          B, packed columnwise in a linear array.  The j-th column of B
 109 *>          is stored in the array BP as follows:
 110 *>          if UPLO = 'U', BP(i + (j-1)*j/2) = B(i,j) for 1<=i<=j;
 111 *>          if UPLO = 'L', BP(i + (j-1)*(2*n-j)/2) = B(i,j) for j<=i<=n.
 112 *>
 113 *>          On exit, the triangular factor U or L from the Cholesky
 114 *>          factorization B = U**H*U or B = L*L**H, in the same storage
 115 *>          format as B.
 116 *> \endverbatim
 117 *>
 118 *> \param[in] VL
 119 *> \verbatim
 120 *>          VL is DOUBLE PRECISION
 121 *>
 122 *>          If RANGE='V', the lower bound of the interval to
 123 *>          be searched for eigenvalues. VL < VU.
 124 *>          Not referenced if RANGE = 'A' or 'I'.
 125 *> \endverbatim
 126 *>
 127 *> \param[in] VU
 128 *> \verbatim
 129 *>          VU is DOUBLE PRECISION
 130 *>
 131 *>          If RANGE='V', the upper bound of the interval to
 132 *>          be searched for eigenvalues. VL < VU.
 133 *>          Not referenced if RANGE = 'A' or 'I'.
 134 *> \endverbatim
 135 *>
 136 *> \param[in] IL
 137 *> \verbatim
 138 *>          IL is INTEGER
 139 *>
 140 *>          If RANGE='I', the index of the
 141 *>          smallest eigenvalue to be returned.
 142 *>          1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
 143 *>          Not referenced if RANGE = 'A' or 'V'.
 144 *> \endverbatim
 145 *>
 146 *> \param[in] IU
 147 *> \verbatim
 148 *>          IU is INTEGER
 149 *>
 150 *>          If RANGE='I', the index of the
 151 *>          largest eigenvalue to be returned.
 152 *>          1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
 153 *>          Not referenced if RANGE = 'A' or 'V'.
 154 *> \endverbatim
 155 *>
 156 *> \param[in] ABSTOL
 157 *> \verbatim
 158 *>          ABSTOL is DOUBLE PRECISION
 159 *>          The absolute error tolerance for the eigenvalues.
 160 *>          An approximate eigenvalue is accepted as converged
 161 *>          when it is determined to lie in an interval [a,b]
 162 *>          of width less than or equal to
 163 *>
 164 *>                  ABSTOL + EPS *   max( |a|,|b| ) ,
 165 *>
 166 *>          where EPS is the machine precision.  If ABSTOL is less than
 167 *>          or equal to zero, then  EPS*|T|  will be used in its place,
 168 *>          where |T| is the 1-norm of the tridiagonal matrix obtained
 169 *>          by reducing AP to tridiagonal form.
 170 *>
 171 *>          Eigenvalues will be computed most accurately when ABSTOL is
 172 *>          set to twice the underflow threshold 2*DLAMCH('S'), not zero.
 173 *>          If this routine returns with INFO>0, indicating that some
 174 *>          eigenvectors did not converge, try setting ABSTOL to
 175 *>          2*DLAMCH('S').
 176 *> \endverbatim
 177 *>
 178 *> \param[out] M
 179 *> \verbatim
 180 *>          M is INTEGER
 181 *>          The total number of eigenvalues found.  0 <= M <= N.
 182 *>          If RANGE = 'A', M = N, and if RANGE = 'I', M = IU-IL+1.
 183 *> \endverbatim
 184 *>
 185 *> \param[out] W
 186 *> \verbatim
 187 *>          W is DOUBLE PRECISION array, dimension (N)
 188 *>          On normal exit, the first M elements contain the selected
 189 *>          eigenvalues in ascending order.
 190 *> \endverbatim
 191 *>
 192 *> \param[out] Z
 193 *> \verbatim
 194 *>          Z is COMPLEX*16 array, dimension (LDZ, N)
 195 *>          If JOBZ = 'N', then Z is not referenced.
 196 *>          If JOBZ = 'V', then if INFO = 0, the first M columns of Z
 197 *>          contain the orthonormal eigenvectors of the matrix A
 198 *>          corresponding to the selected eigenvalues, with the i-th
 199 *>          column of Z holding the eigenvector associated with W(i).
 200 *>          The eigenvectors are normalized as follows:
 201 *>          if ITYPE = 1 or 2, Z**H*B*Z = I;
 202 *>          if ITYPE = 3, Z**H*inv(B)*Z = I.
 203 *>
 204 *>          If an eigenvector fails to converge, then that column of Z
 205 *>          contains the latest approximation to the eigenvector, and the
 206 *>          index of the eigenvector is returned in IFAIL.
 207 *>          Note: the user must ensure that at least max(1,M) columns are
 208 *>          supplied in the array Z; if RANGE = 'V', the exact value of M
 209 *>          is not known in advance and an upper bound must be used.
 210 *> \endverbatim
 211 *>
 212 *> \param[in] LDZ
 213 *> \verbatim
 214 *>          LDZ is INTEGER
 215 *>          The leading dimension of the array Z.  LDZ >= 1, and if
 216 *>          JOBZ = 'V', LDZ >= max(1,N).
 217 *> \endverbatim
 218 *>
 219 *> \param[out] WORK
 220 *> \verbatim
 221 *>          WORK is COMPLEX*16 array, dimension (2*N)
 222 *> \endverbatim
 223 *>
 224 *> \param[out] RWORK
 225 *> \verbatim
 226 *>          RWORK is DOUBLE PRECISION array, dimension (7*N)
 227 *> \endverbatim
 228 *>
 229 *> \param[out] IWORK
 230 *> \verbatim
 231 *>          IWORK is INTEGER array, dimension (5*N)
 232 *> \endverbatim
 233 *>
 234 *> \param[out] IFAIL
 235 *> \verbatim
 236 *>          IFAIL is INTEGER array, dimension (N)
 237 *>          If JOBZ = 'V', then if INFO = 0, the first M elements of
 238 *>          IFAIL are zero.  If INFO > 0, then IFAIL contains the
 239 *>          indices of the eigenvectors that failed to converge.
 240 *>          If JOBZ = 'N', then IFAIL is not referenced.
 241 *> \endverbatim
 242 *>
 243 *> \param[out] INFO
 244 *> \verbatim
 245 *>          INFO is INTEGER
 246 *>          = 0:  successful exit
 247 *>          < 0:  if INFO = -i, the i-th argument had an illegal value
 248 *>          > 0:  ZPPTRF or ZHPEVX returned an error code:
 249 *>             <= N:  if INFO = i, ZHPEVX failed to converge;
 250 *>                    i eigenvectors failed to converge.  Their indices
 251 *>                    are stored in array IFAIL.
 252 *>             > N:   if INFO = N + i, for 1 <= i <= n, then the leading
 253 *>                    minor of order i of B is not positive definite.
 254 *>                    The factorization of B could not be completed and
 255 *>                    no eigenvalues or eigenvectors were computed.
 256 *> \endverbatim
 257 *
 258 *  Authors:
 259 *  ========
 260 *
 261 *> \author Univ. of Tennessee
 262 *> \author Univ. of California Berkeley
 263 *> \author Univ. of Colorado Denver
 264 *> \author NAG Ltd.
 265 *
 266 *> \date June 2016
 267 *
 268 *> \ingroup complex16OTHEReigen
 269 *
 270 *> \par Contributors:
 271 *  ==================
 272 *>
 273 *>     Mark Fahey, Department of Mathematics, Univ. of Kentucky, USA
 274 *
 275 *  =====================================================================
 276       SUBROUTINE ZHPGVX( ITYPE, JOBZ, RANGE, UPLO, N, AP, BP, VL, VU,
 277      $                   IL, IU, ABSTOL, M, W, Z, LDZ, WORK, RWORK,
 278      $                   IWORK, IFAIL, INFO )
 279 *
 280 *  -- LAPACK driver routine (version 3.6.1) --
 281 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 282 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 283 *     June 2016
 284 *
 285 *     .. Scalar Arguments ..
 286       CHARACTER          JOBZ, RANGE, UPLO
 287       INTEGER            IL, INFO, ITYPE, IU, LDZ, M, N
 288       DOUBLE PRECISION   ABSTOL, VL, VU
 289 *     ..
 290 *     .. Array Arguments ..
 291       INTEGER            IFAIL( * ), IWORK( * )
 292       DOUBLE PRECISION   RWORK( * ), W( * )
 293       COMPLEX*16         AP( * ), BP( * ), WORK( * ), Z( LDZ, * )
 294 *     ..
 295 *
 296 *  =====================================================================
 297 *
 298 *     .. Local Scalars ..
 299       LOGICAL            ALLEIG, INDEIG, UPPER, VALEIG, WANTZ
 300       CHARACTER          TRANS
 301       INTEGER            J
 302 *     ..
 303 *     .. External Functions ..
 304       LOGICAL            LSAME
 305       EXTERNAL           LSAME
 306 *     ..
 307 *     .. External Subroutines ..
 308       EXTERNAL           XERBLA, ZHPEVX, ZHPGST, ZPPTRF, ZTPMV, ZTPSV
 309 *     ..
 310 *     .. Intrinsic Functions ..
 311       INTRINSIC          MIN
 312 *     ..
 313 *     .. Executable Statements ..
 314 *
 315 *     Test the input parameters.
 316 *
 317       WANTZ = LSAME( JOBZ, 'V' )
 318       UPPER = LSAME( UPLO, 'U' )
 319       ALLEIG = LSAME( RANGE, 'A' )
 320       VALEIG = LSAME( RANGE, 'V' )
 321       INDEIG = LSAME( RANGE, 'I' )
 322 *
 323       INFO = 0
 324       IF( ITYPE.LT.1 .OR. ITYPE.GT.3 ) THEN
 325          INFO = -1
 326       ELSE IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN
 327          INFO = -2
 328       ELSE IF( .NOT.( ALLEIG .OR. VALEIG .OR. INDEIG ) ) THEN
 329          INFO = -3
 330       ELSE IF( .NOT.( UPPER .OR. LSAME( UPLO, 'L' ) ) ) THEN
 331          INFO = -4
 332       ELSE IF( N.LT.0 ) THEN
 333          INFO = -5
 334       ELSE
 335          IF( VALEIG ) THEN
 336             IF( N.GT.0 .AND. VU.LE.VL ) THEN
 337                INFO = -9
 338             END IF
 339          ELSE IF( INDEIG ) THEN
 340             IF( IL.LT.1 ) THEN
 341                INFO = -10
 342             ELSE IF( IU.LT.MIN( N, IL ) .OR. IU.GT.N ) THEN
 343                INFO = -11
 344             END IF
 345          END IF
 346       END IF
 347       IF( INFO.EQ.0 ) THEN
 348          IF( LDZ.LT.1 .OR. ( WANTZ .AND. LDZ.LT.N ) ) THEN
 349             INFO = -16
 350          END IF
 351       END IF
 352 *
 353       IF( INFO.NE.0 ) THEN
 354          CALL XERBLA( 'ZHPGVX', -INFO )
 355          RETURN
 356       END IF
 357 *
 358 *     Quick return if possible
 359 *
 360       IF( N.EQ.0 )
 361      $   RETURN
 362 *
 363 *     Form a Cholesky factorization of B.
 364 *
 365       CALL ZPPTRF( UPLO, N, BP, INFO )
 366       IF( INFO.NE.0 ) THEN
 367          INFO = N + INFO
 368          RETURN
 369       END IF
 370 *
 371 *     Transform problem to standard eigenvalue problem and solve.
 372 *
 373       CALL ZHPGST( ITYPE, UPLO, N, AP, BP, INFO )
 374       CALL ZHPEVX( JOBZ, RANGE, UPLO, N, AP, VL, VU, IL, IU, ABSTOL, M,
 375      $             W, Z, LDZ, WORK, RWORK, IWORK, IFAIL, INFO )
 376 *
 377       IF( WANTZ ) THEN
 378 *
 379 *        Backtransform eigenvectors to the original problem.
 380 *
 381          IF( INFO.GT.0 )
 382      $      M = INFO - 1
 383          IF( ITYPE.EQ.1 .OR. ITYPE.EQ.2 ) THEN
 384 *
 385 *           For A*x=(lambda)*B*x and A*B*x=(lambda)*x;
 386 *           backtransform eigenvectors: x = inv(L)**H *y or inv(U)*y
 387 *
 388             IF( UPPER ) THEN
 389                TRANS = 'N'
 390             ELSE
 391                TRANS = 'C'
 392             END IF
 393 *
 394             DO 10 J = 1, M
 395                CALL ZTPSV( UPLO, TRANS, 'Non-unit', N, BP, Z( 1, J ),
 396      $                     1 )
 397    10       CONTINUE
 398 *
 399          ELSE IF( ITYPE.EQ.3 ) THEN
 400 *
 401 *           For B*A*x=(lambda)*x;
 402 *           backtransform eigenvectors: x = L*y or U**H *y
 403 *
 404             IF( UPPER ) THEN
 405                TRANS = 'C'
 406             ELSE
 407                TRANS = 'N'
 408             END IF
 409 *
 410             DO 20 J = 1, M
 411                CALL ZTPMV( UPLO, TRANS, 'Non-unit', N, BP, Z( 1, J ),
 412      $                     1 )
 413    20       CONTINUE
 414          END IF
 415       END IF
 416 *
 417       RETURN
 418 *
 419 *     End of ZHPGVX
 420 *
 421       END