SRC/dsbgvd.f

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