SRC/dspgvx.f

   1 *> \brief \b DSPGVX
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download DSPGVX + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dspgvx.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dspgvx.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dspgvx.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE DSPGVX( ITYPE, JOBZ, RANGE, UPLO, N, AP, BP, VL, VU,
  22 *                          IL, IU, ABSTOL, M, W, Z, LDZ, WORK, IWORK,
  23 *                          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   AP( * ), BP( * ), W( * ), WORK( * ),
  33 *      $                   Z( LDZ, * )
  34 *       ..
  35 *
  36 *
  37 *> \par Purpose:
  38 *  =============
  39 *>
  40 *> \verbatim
  41 *>
  42 *> DSPGVX computes selected eigenvalues, and optionally, eigenvectors
  43 *> of a real generalized symmetric-definite eigenproblem, of the form
  44 *> A*x=(lambda)*B*x,  A*Bx=(lambda)*x,  or B*A*x=(lambda)*x.  Here A
  45 *> and B are assumed to be symmetric, stored in packed storage, and B
  46 *> is also positive definite.  Eigenvalues and eigenvectors can be
  47 *> selected by specifying either a range of values or a range of indices
  48 *> for the 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 triangle of A and B are stored;
  83 *>          = 'L':  Lower triangle of A and B are stored.
  84 *> \endverbatim
  85 *>
  86 *> \param[in] N
  87 *> \verbatim
  88 *>          N is INTEGER
  89 *>          The order of the matrix pencil (A,B).  N >= 0.
  90 *> \endverbatim
  91 *>
  92 *> \param[in,out] AP
  93 *> \verbatim
  94 *>          AP is DOUBLE PRECISION array, dimension (N*(N+1)/2)
  95 *>          On entry, the upper or lower triangle of the symmetric 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 DOUBLE PRECISION array, dimension (N*(N+1)/2)
 107 *>          On entry, the upper or lower triangle of the symmetric 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**T*U or B = L*L**T, 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 A 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 DOUBLE PRECISION array, dimension (LDZ, max(1,M))
 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**T*B*Z = I;
 202 *>          if ITYPE = 3, Z**T*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 DOUBLE PRECISION array, dimension (8*N)
 222 *> \endverbatim
 223 *>
 224 *> \param[out] IWORK
 225 *> \verbatim
 226 *>          IWORK is INTEGER array, dimension (5*N)
 227 *> \endverbatim
 228 *>
 229 *> \param[out] IFAIL
 230 *> \verbatim
 231 *>          IFAIL is INTEGER array, dimension (N)
 232 *>          If JOBZ = 'V', then if INFO = 0, the first M elements of
 233 *>          IFAIL are zero.  If INFO > 0, then IFAIL contains the
 234 *>          indices of the eigenvectors that failed to converge.
 235 *>          If JOBZ = 'N', then IFAIL is not referenced.
 236 *> \endverbatim
 237 *>
 238 *> \param[out] INFO
 239 *> \verbatim
 240 *>          INFO is INTEGER
 241 *>          = 0:  successful exit
 242 *>          < 0:  if INFO = -i, the i-th argument had an illegal value
 243 *>          > 0:  DPPTRF or DSPEVX returned an error code:
 244 *>             <= N:  if INFO = i, DSPEVX failed to converge;
 245 *>                    i eigenvectors failed to converge.  Their indices
 246 *>                    are stored in array IFAIL.
 247 *>             > N:   if INFO = N + i, for 1 <= i <= N, then the leading
 248 *>                    minor of order i of B is not positive definite.
 249 *>                    The factorization of B could not be completed and
 250 *>                    no eigenvalues or eigenvectors were computed.
 251 *> \endverbatim
 252 *
 253 *  Authors:
 254 *  ========
 255 *
 256 *> \author Univ. of Tennessee
 257 *> \author Univ. of California Berkeley
 258 *> \author Univ. of Colorado Denver
 259 *> \author NAG Ltd.
 260 *
 261 *> \date June 2016
 262 *
 263 *> \ingroup doubleOTHEReigen
 264 *
 265 *> \par Contributors:
 266 *  ==================
 267 *>
 268 *>     Mark Fahey, Department of Mathematics, Univ. of Kentucky, USA
 269 *
 270 *  =====================================================================
 271       SUBROUTINE DSPGVX( ITYPE, JOBZ, RANGE, UPLO, N, AP, BP, VL, VU,
 272      $                   IL, IU, ABSTOL, M, W, Z, LDZ, WORK, IWORK,
 273      $                   IFAIL, INFO )
 274 *
 275 *  -- LAPACK driver routine (version 3.6.1) --
 276 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 277 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 278 *     June 2016
 279 *
 280 *     .. Scalar Arguments ..
 281       CHARACTER          JOBZ, RANGE, UPLO
 282       INTEGER            IL, INFO, ITYPE, IU, LDZ, M, N
 283       DOUBLE PRECISION   ABSTOL, VL, VU
 284 *     ..
 285 *     .. Array Arguments ..
 286       INTEGER            IFAIL( * ), IWORK( * )
 287       DOUBLE PRECISION   AP( * ), BP( * ), W( * ), WORK( * ),
 288      $                   Z( LDZ, * )
 289 *     ..
 290 *
 291 * =====================================================================
 292 *
 293 *     .. Local Scalars ..
 294       LOGICAL            ALLEIG, INDEIG, UPPER, VALEIG, WANTZ
 295       CHARACTER          TRANS
 296       INTEGER            J
 297 *     ..
 298 *     .. External Functions ..
 299       LOGICAL            LSAME
 300       EXTERNAL           LSAME
 301 *     ..
 302 *     .. External Subroutines ..
 303       EXTERNAL           DPPTRF, DSPEVX, DSPGST, DTPMV, DTPSV, XERBLA
 304 *     ..
 305 *     .. Intrinsic Functions ..
 306       INTRINSIC          MIN
 307 *     ..
 308 *     .. Executable Statements ..
 309 *
 310 *     Test the input parameters.
 311 *
 312       UPPER = LSAME( UPLO, 'U' )
 313       WANTZ = LSAME( JOBZ, 'V' )
 314       ALLEIG = LSAME( RANGE, 'A' )
 315       VALEIG = LSAME( RANGE, 'V' )
 316       INDEIG = LSAME( RANGE, 'I' )
 317 *
 318       INFO = 0
 319       IF( ITYPE.LT.1 .OR. ITYPE.GT.3 ) THEN
 320          INFO = -1
 321       ELSE IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN
 322          INFO = -2
 323       ELSE IF( .NOT.( ALLEIG .OR. VALEIG .OR. INDEIG ) ) THEN
 324          INFO = -3
 325       ELSE IF( .NOT.( UPPER .OR. LSAME( UPLO, 'L' ) ) ) THEN
 326          INFO = -4
 327       ELSE IF( N.LT.0 ) THEN
 328          INFO = -5
 329       ELSE
 330          IF( VALEIG ) THEN
 331             IF( N.GT.0 .AND. VU.LE.VL ) THEN
 332                INFO = -9
 333             END IF
 334          ELSE IF( INDEIG ) THEN
 335             IF( IL.LT.1 ) THEN
 336                INFO = -10
 337             ELSE IF( IU.LT.MIN( N, IL ) .OR. IU.GT.N ) THEN
 338                INFO = -11
 339             END IF
 340          END IF
 341       END IF
 342       IF( INFO.EQ.0 ) THEN
 343          IF( LDZ.LT.1 .OR. ( WANTZ .AND. LDZ.LT.N ) ) THEN
 344             INFO = -16
 345          END IF
 346       END IF
 347 *
 348       IF( INFO.NE.0 ) THEN
 349          CALL XERBLA( 'DSPGVX', -INFO )
 350          RETURN
 351       END IF
 352 *
 353 *     Quick return if possible
 354 *
 355       M = 0
 356       IF( N.EQ.0 )
 357      $   RETURN
 358 *
 359 *     Form a Cholesky factorization of B.
 360 *
 361       CALL DPPTRF( UPLO, N, BP, INFO )
 362       IF( INFO.NE.0 ) THEN
 363          INFO = N + INFO
 364          RETURN
 365       END IF
 366 *
 367 *     Transform problem to standard eigenvalue problem and solve.
 368 *
 369       CALL DSPGST( ITYPE, UPLO, N, AP, BP, INFO )
 370       CALL DSPEVX( JOBZ, RANGE, UPLO, N, AP, VL, VU, IL, IU, ABSTOL, M,
 371      $             W, Z, LDZ, WORK, IWORK, IFAIL, INFO )
 372 *
 373       IF( WANTZ ) THEN
 374 *
 375 *        Backtransform eigenvectors to the original problem.
 376 *
 377          IF( INFO.GT.0 )
 378      $      M = INFO - 1
 379          IF( ITYPE.EQ.1 .OR. ITYPE.EQ.2 ) THEN
 380 *
 381 *           For A*x=(lambda)*B*x and A*B*x=(lambda)*x;
 382 *           backtransform eigenvectors: x = inv(L)**T*y or inv(U)*y
 383 *
 384             IF( UPPER ) THEN
 385                TRANS = 'N'
 386             ELSE
 387                TRANS = 'T'
 388             END IF
 389 *
 390             DO 10 J = 1, M
 391                CALL DTPSV( UPLO, TRANS, 'Non-unit', N, BP, Z( 1, J ),
 392      $                     1 )
 393    10       CONTINUE
 394 *
 395          ELSE IF( ITYPE.EQ.3 ) THEN
 396 *
 397 *           For B*A*x=(lambda)*x;
 398 *           backtransform eigenvectors: x = L*y or U**T*y
 399 *
 400             IF( UPPER ) THEN
 401                TRANS = 'T'
 402             ELSE
 403                TRANS = 'N'
 404             END IF
 405 *
 406             DO 20 J = 1, M
 407                CALL DTPMV( UPLO, TRANS, 'Non-unit', N, BP, Z( 1, J ),
 408      $                     1 )
 409    20       CONTINUE
 410          END IF
 411       END IF
 412 *
 413       RETURN
 414 *
 415 *     End of DSPGVX
 416 *
 417       END