LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
sstevx.f
Go to the documentation of this file.
00001 *> \brief <b> SSTEVX computes the eigenvalues and, optionally, the left and/or right eigenvectors for OTHER matrices</b>
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *> \htmlonly
00009 *> Download SSTEVX + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/sstevx.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/sstevx.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/sstevx.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       SUBROUTINE SSTEVX( JOBZ, RANGE, N, D, E, VL, VU, IL, IU, ABSTOL,
00022 *                          M, W, Z, LDZ, WORK, IWORK, IFAIL, INFO )
00023 * 
00024 *       .. Scalar Arguments ..
00025 *       CHARACTER          JOBZ, RANGE
00026 *       INTEGER            IL, INFO, IU, LDZ, M, N
00027 *       REAL               ABSTOL, VL, VU
00028 *       ..
00029 *       .. Array Arguments ..
00030 *       INTEGER            IFAIL( * ), IWORK( * )
00031 *       REAL               D( * ), E( * ), W( * ), WORK( * ), Z( LDZ, * )
00032 *       ..
00033 *  
00034 *
00035 *> \par Purpose:
00036 *  =============
00037 *>
00038 *> \verbatim
00039 *>
00040 *> SSTEVX computes selected eigenvalues and, optionally, eigenvectors
00041 *> of a real symmetric tridiagonal matrix A.  Eigenvalues and
00042 *> eigenvectors can be selected by specifying either a range of values
00043 *> or a range of indices for the desired eigenvalues.
00044 *> \endverbatim
00045 *
00046 *  Arguments:
00047 *  ==========
00048 *
00049 *> \param[in] JOBZ
00050 *> \verbatim
00051 *>          JOBZ is CHARACTER*1
00052 *>          = 'N':  Compute eigenvalues only;
00053 *>          = 'V':  Compute eigenvalues and eigenvectors.
00054 *> \endverbatim
00055 *>
00056 *> \param[in] RANGE
00057 *> \verbatim
00058 *>          RANGE is CHARACTER*1
00059 *>          = 'A': all eigenvalues will be found.
00060 *>          = 'V': all eigenvalues in the half-open interval (VL,VU]
00061 *>                 will be found.
00062 *>          = 'I': the IL-th through IU-th eigenvalues will be found.
00063 *> \endverbatim
00064 *>
00065 *> \param[in] N
00066 *> \verbatim
00067 *>          N is INTEGER
00068 *>          The order of the matrix.  N >= 0.
00069 *> \endverbatim
00070 *>
00071 *> \param[in,out] D
00072 *> \verbatim
00073 *>          D is REAL array, dimension (N)
00074 *>          On entry, the n diagonal elements of the tridiagonal matrix
00075 *>          A.
00076 *>          On exit, D may be multiplied by a constant factor chosen
00077 *>          to avoid over/underflow in computing the eigenvalues.
00078 *> \endverbatim
00079 *>
00080 *> \param[in,out] E
00081 *> \verbatim
00082 *>          E is REAL array, dimension (max(1,N-1))
00083 *>          On entry, the (n-1) subdiagonal elements of the tridiagonal
00084 *>          matrix A in elements 1 to N-1 of E.
00085 *>          On exit, E may be multiplied by a constant factor chosen
00086 *>          to avoid over/underflow in computing the eigenvalues.
00087 *> \endverbatim
00088 *>
00089 *> \param[in] VL
00090 *> \verbatim
00091 *>          VL is REAL
00092 *> \endverbatim
00093 *>
00094 *> \param[in] VU
00095 *> \verbatim
00096 *>          VU is REAL
00097 *>          If RANGE='V', the lower and upper bounds of the interval to
00098 *>          be searched for eigenvalues. VL < VU.
00099 *>          Not referenced if RANGE = 'A' or 'I'.
00100 *> \endverbatim
00101 *>
00102 *> \param[in] IL
00103 *> \verbatim
00104 *>          IL is INTEGER
00105 *> \endverbatim
00106 *>
00107 *> \param[in] IU
00108 *> \verbatim
00109 *>          IU is INTEGER
00110 *>          If RANGE='I', the indices (in ascending order) of the
00111 *>          smallest and largest eigenvalues to be returned.
00112 *>          1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
00113 *>          Not referenced if RANGE = 'A' or 'V'.
00114 *> \endverbatim
00115 *>
00116 *> \param[in] ABSTOL
00117 *> \verbatim
00118 *>          ABSTOL is REAL
00119 *>          The absolute error tolerance for the eigenvalues.
00120 *>          An approximate eigenvalue is accepted as converged
00121 *>          when it is determined to lie in an interval [a,b]
00122 *>          of width less than or equal to
00123 *>
00124 *>                  ABSTOL + EPS *   max( |a|,|b| ) ,
00125 *>
00126 *>          where EPS is the machine precision.  If ABSTOL is less
00127 *>          than or equal to zero, then  EPS*|T|  will be used in
00128 *>          its place, where |T| is the 1-norm of the tridiagonal
00129 *>          matrix.
00130 *>
00131 *>          Eigenvalues will be computed most accurately when ABSTOL is
00132 *>          set to twice the underflow threshold 2*SLAMCH('S'), not zero.
00133 *>          If this routine returns with INFO>0, indicating that some
00134 *>          eigenvectors did not converge, try setting ABSTOL to
00135 *>          2*SLAMCH('S').
00136 *>
00137 *>          See "Computing Small Singular Values of Bidiagonal Matrices
00138 *>          with Guaranteed High Relative Accuracy," by Demmel and
00139 *>          Kahan, LAPACK Working Note #3.
00140 *> \endverbatim
00141 *>
00142 *> \param[out] M
00143 *> \verbatim
00144 *>          M is INTEGER
00145 *>          The total number of eigenvalues found.  0 <= M <= N.
00146 *>          If RANGE = 'A', M = N, and if RANGE = 'I', M = IU-IL+1.
00147 *> \endverbatim
00148 *>
00149 *> \param[out] W
00150 *> \verbatim
00151 *>          W is REAL array, dimension (N)
00152 *>          The first M elements contain the selected eigenvalues in
00153 *>          ascending order.
00154 *> \endverbatim
00155 *>
00156 *> \param[out] Z
00157 *> \verbatim
00158 *>          Z is REAL array, dimension (LDZ, max(1,M) )
00159 *>          If JOBZ = 'V', then if INFO = 0, the first M columns of Z
00160 *>          contain the orthonormal eigenvectors of the matrix A
00161 *>          corresponding to the selected eigenvalues, with the i-th
00162 *>          column of Z holding the eigenvector associated with W(i).
00163 *>          If an eigenvector fails to converge (INFO > 0), then that
00164 *>          column of Z contains the latest approximation to the
00165 *>          eigenvector, and the index of the eigenvector is returned
00166 *>          in IFAIL.  If JOBZ = 'N', then Z is not referenced.
00167 *>          Note: the user must ensure that at least max(1,M) columns are
00168 *>          supplied in the array Z; if RANGE = 'V', the exact value of M
00169 *>          is not known in advance and an upper bound must be used.
00170 *> \endverbatim
00171 *>
00172 *> \param[in] LDZ
00173 *> \verbatim
00174 *>          LDZ is INTEGER
00175 *>          The leading dimension of the array Z.  LDZ >= 1, and if
00176 *>          JOBZ = 'V', LDZ >= max(1,N).
00177 *> \endverbatim
00178 *>
00179 *> \param[out] WORK
00180 *> \verbatim
00181 *>          WORK is REAL array, dimension (5*N)
00182 *> \endverbatim
00183 *>
00184 *> \param[out] IWORK
00185 *> \verbatim
00186 *>          IWORK is INTEGER array, dimension (5*N)
00187 *> \endverbatim
00188 *>
00189 *> \param[out] IFAIL
00190 *> \verbatim
00191 *>          IFAIL is INTEGER array, dimension (N)
00192 *>          If JOBZ = 'V', then if INFO = 0, the first M elements of
00193 *>          IFAIL are zero.  If INFO > 0, then IFAIL contains the
00194 *>          indices of the eigenvectors that failed to converge.
00195 *>          If JOBZ = 'N', then IFAIL is not referenced.
00196 *> \endverbatim
00197 *>
00198 *> \param[out] INFO
00199 *> \verbatim
00200 *>          INFO is INTEGER
00201 *>          = 0:  successful exit
00202 *>          < 0:  if INFO = -i, the i-th argument had an illegal value
00203 *>          > 0:  if INFO = i, then i eigenvectors failed to converge.
00204 *>                Their indices are stored in array IFAIL.
00205 *> \endverbatim
00206 *
00207 *  Authors:
00208 *  ========
00209 *
00210 *> \author Univ. of Tennessee 
00211 *> \author Univ. of California Berkeley 
00212 *> \author Univ. of Colorado Denver 
00213 *> \author NAG Ltd. 
00214 *
00215 *> \date November 2011
00216 *
00217 *> \ingroup realOTHEReigen
00218 *
00219 *  =====================================================================
00220       SUBROUTINE SSTEVX( JOBZ, RANGE, N, D, E, VL, VU, IL, IU, ABSTOL,
00221      $                   M, W, Z, LDZ, WORK, IWORK, IFAIL, INFO )
00222 *
00223 *  -- LAPACK driver routine (version 3.4.0) --
00224 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00225 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00226 *     November 2011
00227 *
00228 *     .. Scalar Arguments ..
00229       CHARACTER          JOBZ, RANGE
00230       INTEGER            IL, INFO, IU, LDZ, M, N
00231       REAL               ABSTOL, VL, VU
00232 *     ..
00233 *     .. Array Arguments ..
00234       INTEGER            IFAIL( * ), IWORK( * )
00235       REAL               D( * ), E( * ), W( * ), WORK( * ), Z( LDZ, * )
00236 *     ..
00237 *
00238 *  =====================================================================
00239 *
00240 *     .. Parameters ..
00241       REAL               ZERO, ONE
00242       PARAMETER          ( ZERO = 0.0E0, ONE = 1.0E0 )
00243 *     ..
00244 *     .. Local Scalars ..
00245       LOGICAL            ALLEIG, INDEIG, TEST, VALEIG, WANTZ
00246       CHARACTER          ORDER
00247       INTEGER            I, IMAX, INDIBL, INDISP, INDIWO, INDWRK,
00248      $                   ISCALE, ITMP1, J, JJ, NSPLIT
00249       REAL               BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA, SMLNUM,
00250      $                   TMP1, TNRM, VLL, VUU
00251 *     ..
00252 *     .. External Functions ..
00253       LOGICAL            LSAME
00254       REAL               SLAMCH, SLANST
00255       EXTERNAL           LSAME, SLAMCH, SLANST
00256 *     ..
00257 *     .. External Subroutines ..
00258       EXTERNAL           SCOPY, SSCAL, SSTEBZ, SSTEIN, SSTEQR, SSTERF,
00259      $                   SSWAP, XERBLA
00260 *     ..
00261 *     .. Intrinsic Functions ..
00262       INTRINSIC          MAX, MIN, SQRT
00263 *     ..
00264 *     .. Executable Statements ..
00265 *
00266 *     Test the input parameters.
00267 *
00268       WANTZ = LSAME( JOBZ, 'V' )
00269       ALLEIG = LSAME( RANGE, 'A' )
00270       VALEIG = LSAME( RANGE, 'V' )
00271       INDEIG = LSAME( RANGE, 'I' )
00272 *
00273       INFO = 0
00274       IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN
00275          INFO = -1
00276       ELSE IF( .NOT.( ALLEIG .OR. VALEIG .OR. INDEIG ) ) THEN
00277          INFO = -2
00278       ELSE IF( N.LT.0 ) THEN
00279          INFO = -3
00280       ELSE
00281          IF( VALEIG ) THEN
00282             IF( N.GT.0 .AND. VU.LE.VL )
00283      $         INFO = -7
00284          ELSE IF( INDEIG ) THEN
00285             IF( IL.LT.1 .OR. IL.GT.MAX( 1, N ) ) THEN
00286                INFO = -8
00287             ELSE IF( IU.LT.MIN( N, IL ) .OR. IU.GT.N ) THEN
00288                INFO = -9
00289             END IF
00290          END IF
00291       END IF
00292       IF( INFO.EQ.0 ) THEN
00293          IF( LDZ.LT.1 .OR. ( WANTZ .AND. LDZ.LT.N ) )
00294      $      INFO = -14
00295       END IF
00296 *
00297       IF( INFO.NE.0 ) THEN
00298          CALL XERBLA( 'SSTEVX', -INFO )
00299          RETURN
00300       END IF
00301 *
00302 *     Quick return if possible
00303 *
00304       M = 0
00305       IF( N.EQ.0 )
00306      $   RETURN
00307 *
00308       IF( N.EQ.1 ) THEN
00309          IF( ALLEIG .OR. INDEIG ) THEN
00310             M = 1
00311             W( 1 ) = D( 1 )
00312          ELSE
00313             IF( VL.LT.D( 1 ) .AND. VU.GE.D( 1 ) ) THEN
00314                M = 1
00315                W( 1 ) = D( 1 )
00316             END IF
00317          END IF
00318          IF( WANTZ )
00319      $      Z( 1, 1 ) = ONE
00320          RETURN
00321       END IF
00322 *
00323 *     Get machine constants.
00324 *
00325       SAFMIN = SLAMCH( 'Safe minimum' )
00326       EPS = SLAMCH( 'Precision' )
00327       SMLNUM = SAFMIN / EPS
00328       BIGNUM = ONE / SMLNUM
00329       RMIN = SQRT( SMLNUM )
00330       RMAX = MIN( SQRT( BIGNUM ), ONE / SQRT( SQRT( SAFMIN ) ) )
00331 *
00332 *     Scale matrix to allowable range, if necessary.
00333 *
00334       ISCALE = 0
00335       IF ( VALEIG ) THEN
00336          VLL = VL
00337          VUU = VU
00338       ELSE
00339          VLL = ZERO
00340          VUU = ZERO
00341       ENDIF
00342       TNRM = SLANST( 'M', N, D, E )
00343       IF( TNRM.GT.ZERO .AND. TNRM.LT.RMIN ) THEN
00344          ISCALE = 1
00345          SIGMA = RMIN / TNRM
00346       ELSE IF( TNRM.GT.RMAX ) THEN
00347          ISCALE = 1
00348          SIGMA = RMAX / TNRM
00349       END IF
00350       IF( ISCALE.EQ.1 ) THEN
00351          CALL SSCAL( N, SIGMA, D, 1 )
00352          CALL SSCAL( N-1, SIGMA, E( 1 ), 1 )
00353          IF( VALEIG ) THEN
00354             VLL = VL*SIGMA
00355             VUU = VU*SIGMA
00356          END IF
00357       END IF
00358 *
00359 *     If all eigenvalues are desired and ABSTOL is less than zero, then
00360 *     call SSTERF or SSTEQR.  If this fails for some eigenvalue, then
00361 *     try SSTEBZ.
00362 *
00363       TEST = .FALSE.
00364       IF( INDEIG ) THEN
00365          IF( IL.EQ.1 .AND. IU.EQ.N ) THEN
00366             TEST = .TRUE.
00367          END IF
00368       END IF
00369       IF( ( ALLEIG .OR. TEST ) .AND. ( ABSTOL.LE.ZERO ) ) THEN
00370          CALL SCOPY( N, D, 1, W, 1 )
00371          CALL SCOPY( N-1, E( 1 ), 1, WORK( 1 ), 1 )
00372          INDWRK = N + 1
00373          IF( .NOT.WANTZ ) THEN
00374             CALL SSTERF( N, W, WORK, INFO )
00375          ELSE
00376             CALL SSTEQR( 'I', N, W, WORK, Z, LDZ, WORK( INDWRK ), INFO )
00377             IF( INFO.EQ.0 ) THEN
00378                DO 10 I = 1, N
00379                   IFAIL( I ) = 0
00380    10          CONTINUE
00381             END IF
00382          END IF
00383          IF( INFO.EQ.0 ) THEN
00384             M = N
00385             GO TO 20
00386          END IF
00387          INFO = 0
00388       END IF
00389 *
00390 *     Otherwise, call SSTEBZ and, if eigenvectors are desired, SSTEIN.
00391 *
00392       IF( WANTZ ) THEN
00393          ORDER = 'B'
00394       ELSE
00395          ORDER = 'E'
00396       END IF
00397       INDWRK = 1
00398       INDIBL = 1
00399       INDISP = INDIBL + N
00400       INDIWO = INDISP + N
00401       CALL SSTEBZ( RANGE, ORDER, N, VLL, VUU, IL, IU, ABSTOL, D, E, M,
00402      $             NSPLIT, W, IWORK( INDIBL ), IWORK( INDISP ),
00403      $             WORK( INDWRK ), IWORK( INDIWO ), INFO )
00404 *
00405       IF( WANTZ ) THEN
00406          CALL SSTEIN( N, D, E, M, W, IWORK( INDIBL ), IWORK( INDISP ),
00407      $                Z, LDZ, WORK( INDWRK ), IWORK( INDIWO ), IFAIL,
00408      $                INFO )
00409       END IF
00410 *
00411 *     If matrix was scaled, then rescale eigenvalues appropriately.
00412 *
00413    20 CONTINUE
00414       IF( ISCALE.EQ.1 ) THEN
00415          IF( INFO.EQ.0 ) THEN
00416             IMAX = M
00417          ELSE
00418             IMAX = INFO - 1
00419          END IF
00420          CALL SSCAL( IMAX, ONE / SIGMA, W, 1 )
00421       END IF
00422 *
00423 *     If eigenvalues are not in order, then sort them, along with
00424 *     eigenvectors.
00425 *
00426       IF( WANTZ ) THEN
00427          DO 40 J = 1, M - 1
00428             I = 0
00429             TMP1 = W( J )
00430             DO 30 JJ = J + 1, M
00431                IF( W( JJ ).LT.TMP1 ) THEN
00432                   I = JJ
00433                   TMP1 = W( JJ )
00434                END IF
00435    30       CONTINUE
00436 *
00437             IF( I.NE.0 ) THEN
00438                ITMP1 = IWORK( INDIBL+I-1 )
00439                W( I ) = W( J )
00440                IWORK( INDIBL+I-1 ) = IWORK( INDIBL+J-1 )
00441                W( J ) = TMP1
00442                IWORK( INDIBL+J-1 ) = ITMP1
00443                CALL SSWAP( N, Z( 1, I ), 1, Z( 1, J ), 1 )
00444                IF( INFO.NE.0 ) THEN
00445                   ITMP1 = IFAIL( I )
00446                   IFAIL( I ) = IFAIL( J )
00447                   IFAIL( J ) = ITMP1
00448                END IF
00449             END IF
00450    40    CONTINUE
00451       END IF
00452 *
00453       RETURN
00454 *
00455 *     End of SSTEVX
00456 *
00457       END
 All Files Functions