SRC/dspevx.f

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