SRC/zhbgvd.f

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