SRC/zhbgvx.f

   1 *> \brief \b ZHBGVX
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download ZHBGVX + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zhbgvx.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zhbgvx.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zhbgvx.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE ZHBGVX( JOBZ, RANGE, UPLO, N, KA, KB, AB, LDAB, BB,
  22 *                          LDBB, Q, LDQ, VL, VU, IL, IU, ABSTOL, M, W, Z,
  23 *                          LDZ, WORK, RWORK, IWORK, IFAIL, INFO )
  24 *
  25 *       .. Scalar Arguments ..
  26 *       CHARACTER          JOBZ, RANGE, UPLO
  27 *       INTEGER            IL, INFO, IU, KA, KB, LDAB, LDBB, LDQ, LDZ, M,
  28 *      $                   N
  29 *       DOUBLE PRECISION   ABSTOL, VL, VU
  30 *       ..
  31 *       .. Array Arguments ..
  32 *       INTEGER            IFAIL( * ), IWORK( * )
  33 *       DOUBLE PRECISION   RWORK( * ), W( * )
  34 *       COMPLEX*16         AB( LDAB, * ), BB( LDBB, * ), Q( LDQ, * ),
  35 *      $                   WORK( * ), Z( LDZ, * )
  36 *       ..
  37 *
  38 *
  39 *> \par Purpose:
  40 *  =============
  41 *>
  42 *> \verbatim
  43 *>
  44 *> ZHBGVX computes all the eigenvalues, and optionally, the eigenvectors
  45 *> of a complex generalized Hermitian-definite banded eigenproblem, of
  46 *> the form A*x=(lambda)*B*x. Here A and B are assumed to be Hermitian
  47 *> and banded, and B is also positive definite.  Eigenvalues and
  48 *> eigenvectors can be selected by specifying either all eigenvalues,
  49 *> a range of values or a range of indices for the desired eigenvalues.
  50 *> \endverbatim
  51 *
  52 *  Arguments:
  53 *  ==========
  54 *
  55 *> \param[in] JOBZ
  56 *> \verbatim
  57 *>          JOBZ is CHARACTER*1
  58 *>          = 'N':  Compute eigenvalues only;
  59 *>          = 'V':  Compute eigenvalues and eigenvectors.
  60 *> \endverbatim
  61 *>
  62 *> \param[in] RANGE
  63 *> \verbatim
  64 *>          RANGE is CHARACTER*1
  65 *>          = 'A': all eigenvalues will be found;
  66 *>          = 'V': all eigenvalues in the half-open interval (VL,VU]
  67 *>                 will be found;
  68 *>          = 'I': the IL-th through IU-th eigenvalues will be found.
  69 *> \endverbatim
  70 *>
  71 *> \param[in] UPLO
  72 *> \verbatim
  73 *>          UPLO is CHARACTER*1
  74 *>          = 'U':  Upper triangles of A and B are stored;
  75 *>          = 'L':  Lower triangles of A and B are stored.
  76 *> \endverbatim
  77 *>
  78 *> \param[in] N
  79 *> \verbatim
  80 *>          N is INTEGER
  81 *>          The order of the matrices A and B.  N >= 0.
  82 *> \endverbatim
  83 *>
  84 *> \param[in] KA
  85 *> \verbatim
  86 *>          KA is INTEGER
  87 *>          The number of superdiagonals of the matrix A if UPLO = 'U',
  88 *>          or the number of subdiagonals if UPLO = 'L'. KA >= 0.
  89 *> \endverbatim
  90 *>
  91 *> \param[in] KB
  92 *> \verbatim
  93 *>          KB is INTEGER
  94 *>          The number of superdiagonals of the matrix B if UPLO = 'U',
  95 *>          or the number of subdiagonals if UPLO = 'L'. KB >= 0.
  96 *> \endverbatim
  97 *>
  98 *> \param[in,out] AB
  99 *> \verbatim
 100 *>          AB is COMPLEX*16 array, dimension (LDAB, N)
 101 *>          On entry, the upper or lower triangle of the Hermitian band
 102 *>          matrix A, stored in the first ka+1 rows of the array.  The
 103 *>          j-th column of A is stored in the j-th column of the array AB
 104 *>          as follows:
 105 *>          if UPLO = 'U', AB(ka+1+i-j,j) = A(i,j) for max(1,j-ka)<=i<=j;
 106 *>          if UPLO = 'L', AB(1+i-j,j)    = A(i,j) for j<=i<=min(n,j+ka).
 107 *>
 108 *>          On exit, the contents of AB are destroyed.
 109 *> \endverbatim
 110 *>
 111 *> \param[in] LDAB
 112 *> \verbatim
 113 *>          LDAB is INTEGER
 114 *>          The leading dimension of the array AB.  LDAB >= KA+1.
 115 *> \endverbatim
 116 *>
 117 *> \param[in,out] BB
 118 *> \verbatim
 119 *>          BB is COMPLEX*16 array, dimension (LDBB, N)
 120 *>          On entry, the upper or lower triangle of the Hermitian band
 121 *>          matrix B, stored in the first kb+1 rows of the array.  The
 122 *>          j-th column of B is stored in the j-th column of the array BB
 123 *>          as follows:
 124 *>          if UPLO = 'U', BB(kb+1+i-j,j) = B(i,j) for max(1,j-kb)<=i<=j;
 125 *>          if UPLO = 'L', BB(1+i-j,j)    = B(i,j) for j<=i<=min(n,j+kb).
 126 *>
 127 *>          On exit, the factor S from the split Cholesky factorization
 128 *>          B = S**H*S, as returned by ZPBSTF.
 129 *> \endverbatim
 130 *>
 131 *> \param[in] LDBB
 132 *> \verbatim
 133 *>          LDBB is INTEGER
 134 *>          The leading dimension of the array BB.  LDBB >= KB+1.
 135 *> \endverbatim
 136 *>
 137 *> \param[out] Q
 138 *> \verbatim
 139 *>          Q is COMPLEX*16 array, dimension (LDQ, N)
 140 *>          If JOBZ = 'V', the n-by-n matrix used in the reduction of
 141 *>          A*x = (lambda)*B*x to standard form, i.e. C*x = (lambda)*x,
 142 *>          and consequently C to tridiagonal form.
 143 *>          If JOBZ = 'N', the array Q is not referenced.
 144 *> \endverbatim
 145 *>
 146 *> \param[in] LDQ
 147 *> \verbatim
 148 *>          LDQ is INTEGER
 149 *>          The leading dimension of the array Q.  If JOBZ = 'N',
 150 *>          LDQ >= 1. If JOBZ = 'V', LDQ >= max(1,N).
 151 *> \endverbatim
 152 *>
 153 *> \param[in] VL
 154 *> \verbatim
 155 *>          VL is DOUBLE PRECISION
 156 *>
 157 *>          If RANGE='V', the lower bound of the interval to
 158 *>          be searched for eigenvalues. VL < VU.
 159 *>          Not referenced if RANGE = 'A' or 'I'.
 160 *> \endverbatim
 161 *>
 162 *> \param[in] VU
 163 *> \verbatim
 164 *>          VU is DOUBLE PRECISION
 165 *>
 166 *>          If RANGE='V', the upper bound of the interval to
 167 *>          be searched for eigenvalues. VL < VU.
 168 *>          Not referenced if RANGE = 'A' or 'I'.
 169 *> \endverbatim
 170 *>
 171 *> \param[in] IL
 172 *> \verbatim
 173 *>          IL is INTEGER
 174 *>
 175 *>          If RANGE='I', the index of the
 176 *>          smallest eigenvalue to be returned.
 177 *>          1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
 178 *>          Not referenced if RANGE = 'A' or 'V'.
 179 *> \endverbatim
 180 *>
 181 *> \param[in] IU
 182 *> \verbatim
 183 *>          IU is INTEGER
 184 *>
 185 *>          If RANGE='I', the index of the
 186 *>          largest eigenvalue to be returned.
 187 *>          1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
 188 *>          Not referenced if RANGE = 'A' or 'V'.
 189 *> \endverbatim
 190 *>
 191 *> \param[in] ABSTOL
 192 *> \verbatim
 193 *>          ABSTOL is DOUBLE PRECISION
 194 *>          The absolute error tolerance for the eigenvalues.
 195 *>          An approximate eigenvalue is accepted as converged
 196 *>          when it is determined to lie in an interval [a,b]
 197 *>          of width less than or equal to
 198 *>
 199 *>                  ABSTOL + EPS *   max( |a|,|b| ) ,
 200 *>
 201 *>          where EPS is the machine precision.  If ABSTOL is less than
 202 *>          or equal to zero, then  EPS*|T|  will be used in its place,
 203 *>          where |T| is the 1-norm of the tridiagonal matrix obtained
 204 *>          by reducing AP to tridiagonal form.
 205 *>
 206 *>          Eigenvalues will be computed most accurately when ABSTOL is
 207 *>          set to twice the underflow threshold 2*DLAMCH('S'), not zero.
 208 *>          If this routine returns with INFO>0, indicating that some
 209 *>          eigenvectors did not converge, try setting ABSTOL to
 210 *>          2*DLAMCH('S').
 211 *> \endverbatim
 212 *>
 213 *> \param[out] M
 214 *> \verbatim
 215 *>          M is INTEGER
 216 *>          The total number of eigenvalues found.  0 <= M <= N.
 217 *>          If RANGE = 'A', M = N, and if RANGE = 'I', M = IU-IL+1.
 218 *> \endverbatim
 219 *>
 220 *> \param[out] W
 221 *> \verbatim
 222 *>          W is DOUBLE PRECISION array, dimension (N)
 223 *>          If INFO = 0, the eigenvalues in ascending order.
 224 *> \endverbatim
 225 *>
 226 *> \param[out] Z
 227 *> \verbatim
 228 *>          Z is COMPLEX*16 array, dimension (LDZ, N)
 229 *>          If JOBZ = 'V', then if INFO = 0, Z contains the matrix Z of
 230 *>          eigenvectors, with the i-th column of Z holding the
 231 *>          eigenvector associated with W(i). The eigenvectors are
 232 *>          normalized so that Z**H*B*Z = I.
 233 *>          If JOBZ = 'N', then Z is not referenced.
 234 *> \endverbatim
 235 *>
 236 *> \param[in] LDZ
 237 *> \verbatim
 238 *>          LDZ is INTEGER
 239 *>          The leading dimension of the array Z.  LDZ >= 1, and if
 240 *>          JOBZ = 'V', LDZ >= N.
 241 *> \endverbatim
 242 *>
 243 *> \param[out] WORK
 244 *> \verbatim
 245 *>          WORK is COMPLEX*16 array, dimension (N)
 246 *> \endverbatim
 247 *>
 248 *> \param[out] RWORK
 249 *> \verbatim
 250 *>          RWORK is DOUBLE PRECISION array, dimension (7*N)
 251 *> \endverbatim
 252 *>
 253 *> \param[out] IWORK
 254 *> \verbatim
 255 *>          IWORK is INTEGER array, dimension (5*N)
 256 *> \endverbatim
 257 *>
 258 *> \param[out] IFAIL
 259 *> \verbatim
 260 *>          IFAIL is INTEGER array, dimension (N)
 261 *>          If JOBZ = 'V', then if INFO = 0, the first M elements of
 262 *>          IFAIL are zero.  If INFO > 0, then IFAIL contains the
 263 *>          indices of the eigenvectors that failed to converge.
 264 *>          If JOBZ = 'N', then IFAIL is not referenced.
 265 *> \endverbatim
 266 *>
 267 *> \param[out] INFO
 268 *> \verbatim
 269 *>          INFO is INTEGER
 270 *>          = 0:  successful exit
 271 *>          < 0:  if INFO = -i, the i-th argument had an illegal value
 272 *>          > 0:  if INFO = i, and i is:
 273 *>             <= N:  then i eigenvectors failed to converge.  Their
 274 *>                    indices are stored in array IFAIL.
 275 *>             > N:   if INFO = N + i, for 1 <= i <= N, then ZPBSTF
 276 *>                    returned INFO = i: B is not positive definite.
 277 *>                    The factorization of B could not be completed and
 278 *>                    no eigenvalues or eigenvectors were computed.
 279 *> \endverbatim
 280 *
 281 *  Authors:
 282 *  ========
 283 *
 284 *> \author Univ. of Tennessee
 285 *> \author Univ. of California Berkeley
 286 *> \author Univ. of Colorado Denver
 287 *> \author NAG Ltd.
 288 *
 289 *> \date June 2016
 290 *
 291 *> \ingroup complex16OTHEReigen
 292 *
 293 *> \par Contributors:
 294 *  ==================
 295 *>
 296 *>     Mark Fahey, Department of Mathematics, Univ. of Kentucky, USA
 297 *
 298 *  =====================================================================
 299       SUBROUTINE ZHBGVX( JOBZ, RANGE, UPLO, N, KA, KB, AB, LDAB, BB,
 300      $                   LDBB, Q, LDQ, VL, VU, IL, IU, ABSTOL, M, W, Z,
 301      $                   LDZ, WORK, RWORK, IWORK, IFAIL, INFO )
 302 *
 303 *  -- LAPACK driver routine (version 3.6.1) --
 304 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 305 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 306 *     June 2016
 307 *
 308 *     .. Scalar Arguments ..
 309       CHARACTER          JOBZ, RANGE, UPLO
 310       INTEGER            IL, INFO, IU, KA, KB, LDAB, LDBB, LDQ, LDZ, M,
 311      $                   N
 312       DOUBLE PRECISION   ABSTOL, VL, VU
 313 *     ..
 314 *     .. Array Arguments ..
 315       INTEGER            IFAIL( * ), IWORK( * )
 316       DOUBLE PRECISION   RWORK( * ), W( * )
 317       COMPLEX*16         AB( LDAB, * ), BB( LDBB, * ), Q( LDQ, * ),
 318      $                   WORK( * ), Z( LDZ, * )
 319 *     ..
 320 *
 321 *  =====================================================================
 322 *
 323 *     .. Parameters ..
 324       DOUBLE PRECISION   ZERO
 325       PARAMETER          ( ZERO = 0.0D+0 )
 326       COMPLEX*16         CZERO, CONE
 327       PARAMETER          ( CZERO = ( 0.0D+0, 0.0D+0 ),
 328      $                   CONE = ( 1.0D+0, 0.0D+0 ) )
 329 *     ..
 330 *     .. Local Scalars ..
 331       LOGICAL            ALLEIG, INDEIG, TEST, UPPER, VALEIG, WANTZ
 332       CHARACTER          ORDER, VECT
 333       INTEGER            I, IINFO, INDD, INDE, INDEE, INDIBL, INDISP,
 334      $                   INDIWK, INDRWK, INDWRK, ITMP1, J, JJ, NSPLIT
 335       DOUBLE PRECISION   TMP1
 336 *     ..
 337 *     .. External Functions ..
 338       LOGICAL            LSAME
 339       EXTERNAL           LSAME
 340 *     ..
 341 *     .. External Subroutines ..
 342       EXTERNAL           DCOPY, DSTEBZ, DSTERF, XERBLA, ZCOPY, ZGEMV,
 343      $                   ZHBGST, ZHBTRD, ZLACPY, ZPBSTF, ZSTEIN, ZSTEQR,
 344      $                   ZSWAP
 345 *     ..
 346 *     .. Intrinsic Functions ..
 347       INTRINSIC          MIN
 348 *     ..
 349 *     .. Executable Statements ..
 350 *
 351 *     Test the input parameters.
 352 *
 353       WANTZ = LSAME( JOBZ, 'V' )
 354       UPPER = LSAME( UPLO, 'U' )
 355       ALLEIG = LSAME( RANGE, 'A' )
 356       VALEIG = LSAME( RANGE, 'V' )
 357       INDEIG = LSAME( RANGE, 'I' )
 358 *
 359       INFO = 0
 360       IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN
 361          INFO = -1
 362       ELSE IF( .NOT.( ALLEIG .OR. VALEIG .OR. INDEIG ) ) THEN
 363          INFO = -2
 364       ELSE IF( .NOT.( UPPER .OR. LSAME( UPLO, 'L' ) ) ) THEN
 365          INFO = -3
 366       ELSE IF( N.LT.0 ) THEN
 367          INFO = -4
 368       ELSE IF( KA.LT.0 ) THEN
 369          INFO = -5
 370       ELSE IF( KB.LT.0 .OR. KB.GT.KA ) THEN
 371          INFO = -6
 372       ELSE IF( LDAB.LT.KA+1 ) THEN
 373          INFO = -8
 374       ELSE IF( LDBB.LT.KB+1 ) THEN
 375          INFO = -10
 376       ELSE IF( LDQ.LT.1 .OR. ( WANTZ .AND. LDQ.LT.N ) ) THEN
 377          INFO = -12
 378       ELSE
 379          IF( VALEIG ) THEN
 380             IF( N.GT.0 .AND. VU.LE.VL )
 381      $         INFO = -14
 382          ELSE IF( INDEIG ) THEN
 383             IF( IL.LT.1 .OR. IL.GT.MAX( 1, N ) ) THEN
 384                INFO = -15
 385             ELSE IF ( IU.LT.MIN( N, IL ) .OR. IU.GT.N ) THEN
 386                INFO = -16
 387             END IF
 388          END IF
 389       END IF
 390       IF( INFO.EQ.0) THEN
 391          IF( LDZ.LT.1 .OR. ( WANTZ .AND. LDZ.LT.N ) ) THEN
 392             INFO = -21
 393          END IF
 394       END IF
 395 *
 396       IF( INFO.NE.0 ) THEN
 397          CALL XERBLA( 'ZHBGVX', -INFO )
 398          RETURN
 399       END IF
 400 *
 401 *     Quick return if possible
 402 *
 403       M = 0
 404       IF( N.EQ.0 )
 405      $   RETURN
 406 *
 407 *     Form a split Cholesky factorization of B.
 408 *
 409       CALL ZPBSTF( UPLO, N, KB, BB, LDBB, INFO )
 410       IF( INFO.NE.0 ) THEN
 411          INFO = N + INFO
 412          RETURN
 413       END IF
 414 *
 415 *     Transform problem to standard eigenvalue problem.
 416 *
 417       CALL ZHBGST( JOBZ, UPLO, N, KA, KB, AB, LDAB, BB, LDBB, Q, LDQ,
 418      $             WORK, RWORK, IINFO )
 419 *
 420 *     Solve the standard eigenvalue problem.
 421 *     Reduce Hermitian band matrix to tridiagonal form.
 422 *
 423       INDD = 1
 424       INDE = INDD + N
 425       INDRWK = INDE + N
 426       INDWRK = 1
 427       IF( WANTZ ) THEN
 428          VECT = 'U'
 429       ELSE
 430          VECT = 'N'
 431       END IF
 432       CALL ZHBTRD( VECT, UPLO, N, KA, AB, LDAB, RWORK( INDD ),
 433      $             RWORK( INDE ), Q, LDQ, WORK( INDWRK ), IINFO )
 434 *
 435 *     If all eigenvalues are desired and ABSTOL is less than or equal
 436 *     to zero, then call DSTERF or ZSTEQR.  If this fails for some
 437 *     eigenvalue, then try DSTEBZ.
 438 *
 439       TEST = .FALSE.
 440       IF( INDEIG ) THEN
 441          IF( IL.EQ.1 .AND. IU.EQ.N ) THEN
 442             TEST = .TRUE.
 443          END IF
 444       END IF
 445       IF( ( ALLEIG .OR. TEST ) .AND. ( ABSTOL.LE.ZERO ) ) THEN
 446          CALL DCOPY( N, RWORK( INDD ), 1, W, 1 )
 447          INDEE = INDRWK + 2*N
 448          CALL DCOPY( N-1, RWORK( INDE ), 1, RWORK( INDEE ), 1 )
 449          IF( .NOT.WANTZ ) THEN
 450             CALL DSTERF( N, W, RWORK( INDEE ), INFO )
 451          ELSE
 452             CALL ZLACPY( 'A', N, N, Q, LDQ, Z, LDZ )
 453             CALL ZSTEQR( JOBZ, N, W, RWORK( INDEE ), Z, LDZ,
 454      $                   RWORK( INDRWK ), INFO )
 455             IF( INFO.EQ.0 ) THEN
 456                DO 10 I = 1, N
 457                   IFAIL( I ) = 0
 458    10          CONTINUE
 459             END IF
 460          END IF
 461          IF( INFO.EQ.0 ) THEN
 462             M = N
 463             GO TO 30
 464          END IF
 465          INFO = 0
 466       END IF
 467 *
 468 *     Otherwise, call DSTEBZ and, if eigenvectors are desired,
 469 *     call ZSTEIN.
 470 *
 471       IF( WANTZ ) THEN
 472          ORDER = 'B'
 473       ELSE
 474          ORDER = 'E'
 475       END IF
 476       INDIBL = 1
 477       INDISP = INDIBL + N
 478       INDIWK = INDISP + N
 479       CALL DSTEBZ( RANGE, ORDER, N, VL, VU, IL, IU, ABSTOL,
 480      $             RWORK( INDD ), RWORK( INDE ), M, NSPLIT, W,
 481      $             IWORK( INDIBL ), IWORK( INDISP ), RWORK( INDRWK ),
 482      $             IWORK( INDIWK ), INFO )
 483 *
 484       IF( WANTZ ) THEN
 485          CALL ZSTEIN( N, RWORK( INDD ), RWORK( INDE ), M, W,
 486      $                IWORK( INDIBL ), IWORK( INDISP ), Z, LDZ,
 487      $                RWORK( INDRWK ), IWORK( INDIWK ), IFAIL, INFO )
 488 *
 489 *        Apply unitary matrix used in reduction to tridiagonal
 490 *        form to eigenvectors returned by ZSTEIN.
 491 *
 492          DO 20 J = 1, M
 493             CALL ZCOPY( N, Z( 1, J ), 1, WORK( 1 ), 1 )
 494             CALL ZGEMV( 'N', N, N, CONE, Q, LDQ, WORK, 1, CZERO,
 495      $                  Z( 1, J ), 1 )
 496    20    CONTINUE
 497       END IF
 498 *
 499    30 CONTINUE
 500 *
 501 *     If eigenvalues are not in order, then sort them, along with
 502 *     eigenvectors.
 503 *
 504       IF( WANTZ ) THEN
 505          DO 50 J = 1, M - 1
 506             I = 0
 507             TMP1 = W( J )
 508             DO 40 JJ = J + 1, M
 509                IF( W( JJ ).LT.TMP1 ) THEN
 510                   I = JJ
 511                   TMP1 = W( JJ )
 512                END IF
 513    40       CONTINUE
 514 *
 515             IF( I.NE.0 ) THEN
 516                ITMP1 = IWORK( INDIBL+I-1 )
 517                W( I ) = W( J )
 518                IWORK( INDIBL+I-1 ) = IWORK( INDIBL+J-1 )
 519                W( J ) = TMP1
 520                IWORK( INDIBL+J-1 ) = ITMP1
 521                CALL ZSWAP( N, Z( 1, I ), 1, Z( 1, J ), 1 )
 522                IF( INFO.NE.0 ) THEN
 523                   ITMP1 = IFAIL( I )
 524                   IFAIL( I ) = IFAIL( J )
 525                   IFAIL( J ) = ITMP1
 526                END IF
 527             END IF
 528    50    CONTINUE
 529       END IF
 530 *
 531       RETURN
 532 *
 533 *     End of ZHBGVX
 534 *
 535       END