SRC/dsbgvx.f

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