LAPACK: SRC/dstevx.f Source File

Go to the documentation of this file.
00001       SUBROUTINE DSTEVX( JOBZ, RANGE, N, D, E, VL, VU, IL, IU, ABSTOL,
00002      $                   M, W, Z, LDZ, WORK, IWORK, IFAIL, INFO )
00003 *
00004 *  -- LAPACK driver routine (version 3.2) --
00005 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00006 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00007 *     November 2006
00008 *
00009 *     .. Scalar Arguments ..
00010       CHARACTER          JOBZ, RANGE
00011       INTEGER            IL, INFO, IU, LDZ, M, N
00012       DOUBLE PRECISION   ABSTOL, VL, VU
00013 *     ..
00014 *     .. Array Arguments ..
00015       INTEGER            IFAIL( * ), IWORK( * )
00016       DOUBLE PRECISION   D( * ), E( * ), W( * ), WORK( * ), Z( LDZ, * )
00017 *     ..
00018 *
00019 *  Purpose
00020 *  =======
00021 *
00022 *  DSTEVX computes selected eigenvalues and, optionally, eigenvectors
00023 *  of a real symmetric tridiagonal matrix A.  Eigenvalues and
00024 *  eigenvectors can be selected by specifying either a range of values
00025 *  or a range of indices for the desired eigenvalues.
00026 *
00027 *  Arguments
00028 *  =========
00029 *
00030 *  JOBZ    (input) CHARACTER*1
00031 *          = 'N':  Compute eigenvalues only;
00032 *          = 'V':  Compute eigenvalues and eigenvectors.
00033 *
00034 *  RANGE   (input) CHARACTER*1
00035 *          = 'A': all eigenvalues will be found.
00036 *          = 'V': all eigenvalues in the half-open interval (VL,VU]
00037 *                 will be found.
00038 *          = 'I': the IL-th through IU-th eigenvalues will be found.
00039 *
00040 *  N       (input) INTEGER
00041 *          The order of the matrix.  N >= 0.
00042 *
00043 *  D       (input/output) DOUBLE PRECISION array, dimension (N)
00044 *          On entry, the n diagonal elements of the tridiagonal matrix
00045 *          A.
00046 *          On exit, D may be multiplied by a constant factor chosen
00047 *          to avoid over/underflow in computing the eigenvalues.
00048 *
00049 *  E       (input/output) DOUBLE PRECISION array, dimension (max(1,N-1))
00050 *          On entry, the (n-1) subdiagonal elements of the tridiagonal
00051 *          matrix A in elements 1 to N-1 of E.
00052 *          On exit, E may be multiplied by a constant factor chosen
00053 *          to avoid over/underflow in computing the eigenvalues.
00054 *
00055 *  VL      (input) DOUBLE PRECISION
00056 *  VU      (input) DOUBLE PRECISION
00057 *          If RANGE='V', the lower and upper bounds of the interval to
00058 *          be searched for eigenvalues. VL < VU.
00059 *          Not referenced if RANGE = 'A' or 'I'.
00060 *
00061 *  IL      (input) INTEGER
00062 *  IU      (input) INTEGER
00063 *          If RANGE='I', the indices (in ascending order) of the
00064 *          smallest and largest eigenvalues to be returned.
00065 *          1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
00066 *          Not referenced if RANGE = 'A' or 'V'.
00067 *
00068 *  ABSTOL  (input) DOUBLE PRECISION
00069 *          The absolute error tolerance for the eigenvalues.
00070 *          An approximate eigenvalue is accepted as converged
00071 *          when it is determined to lie in an interval [a,b]
00072 *          of width less than or equal to
00073 *
00074 *                  ABSTOL + EPS *   max( |a|,|b| ) ,
00075 *
00076 *          where EPS is the machine precision.  If ABSTOL is less
00077 *          than or equal to zero, then  EPS*|T|  will be used in
00078 *          its place, where |T| is the 1-norm of the tridiagonal
00079 *          matrix.
00080 *
00081 *          Eigenvalues will be computed most accurately when ABSTOL is
00082 *          set to twice the underflow threshold 2*DLAMCH('S'), not zero.
00083 *          If this routine returns with INFO>0, indicating that some
00084 *          eigenvectors did not converge, try setting ABSTOL to
00085 *          2*DLAMCH('S').
00086 *
00087 *          See "Computing Small Singular Values of Bidiagonal Matrices
00088 *          with Guaranteed High Relative Accuracy," by Demmel and
00089 *          Kahan, LAPACK Working Note #3.
00090 *
00091 *  M       (output) INTEGER
00092 *          The total number of eigenvalues found.  0 <= M <= N.
00093 *          If RANGE = 'A', M = N, and if RANGE = 'I', M = IU-IL+1.
00094 *
00095 *  W       (output) DOUBLE PRECISION array, dimension (N)
00096 *          The first M elements contain the selected eigenvalues in
00097 *          ascending order.
00098 *
00099 *  Z       (output) DOUBLE PRECISION array, dimension (LDZ, max(1,M) )
00100 *          If JOBZ = 'V', then if INFO = 0, the first M columns of Z
00101 *          contain the orthonormal eigenvectors of the matrix A
00102 *          corresponding to the selected eigenvalues, with the i-th
00103 *          column of Z holding the eigenvector associated with W(i).
00104 *          If an eigenvector fails to converge (INFO > 0), then that
00105 *          column of Z contains the latest approximation to the
00106 *          eigenvector, and the index of the eigenvector is returned
00107 *          in IFAIL.  If JOBZ = 'N', then Z is not referenced.
00108 *          Note: the user must ensure that at least max(1,M) columns are
00109 *          supplied in the array Z; if RANGE = 'V', the exact value of M
00110 *          is not known in advance and an upper bound must be used.
00111 *
00112 *  LDZ     (input) INTEGER
00113 *          The leading dimension of the array Z.  LDZ >= 1, and if
00114 *          JOBZ = 'V', LDZ >= max(1,N).
00115 *
00116 *  WORK    (workspace) DOUBLE PRECISION array, dimension (5*N)
00117 *
00118 *  IWORK   (workspace) INTEGER array, dimension (5*N)
00119 *
00120 *  IFAIL   (output) INTEGER array, dimension (N)
00121 *          If JOBZ = 'V', then if INFO = 0, the first M elements of
00122 *          IFAIL are zero.  If INFO > 0, then IFAIL contains the
00123 *          indices of the eigenvectors that failed to converge.
00124 *          If JOBZ = 'N', then IFAIL is not referenced.
00125 *
00126 *  INFO    (output) INTEGER
00127 *          = 0:  successful exit
00128 *          < 0:  if INFO = -i, the i-th argument had an illegal value
00129 *          > 0:  if INFO = i, then i eigenvectors failed to converge.
00130 *                Their indices are stored in array IFAIL.
00131 *
00132 *  =====================================================================
00133 *
00134 *     .. Parameters ..
00135       DOUBLE PRECISION   ZERO, ONE
00136       PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0 )
00137 *     ..
00138 *     .. Local Scalars ..
00139       LOGICAL            ALLEIG, INDEIG, TEST, VALEIG, WANTZ
00140       CHARACTER          ORDER
00141       INTEGER            I, IMAX, INDIBL, INDISP, INDIWO, INDWRK,
00142      $                   ISCALE, ITMP1, J, JJ, NSPLIT
00143       DOUBLE PRECISION   BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA, SMLNUM,
00144      $                   TMP1, TNRM, VLL, VUU
00145 *     ..
00146 *     .. External Functions ..
00147       LOGICAL            LSAME
00148       DOUBLE PRECISION   DLAMCH, DLANST
00149       EXTERNAL           LSAME, DLAMCH, DLANST
00150 *     ..
00151 *     .. External Subroutines ..
00152       EXTERNAL           DCOPY, DSCAL, DSTEBZ, DSTEIN, DSTEQR, DSTERF,
00153      $                   DSWAP, XERBLA
00154 *     ..
00155 *     .. Intrinsic Functions ..
00156       INTRINSIC          MAX, MIN, SQRT
00157 *     ..
00158 *     .. Executable Statements ..
00159 *
00160 *     Test the input parameters.
00161 *
00162       WANTZ = LSAME( JOBZ, 'V' )
00163       ALLEIG = LSAME( RANGE, 'A' )
00164       VALEIG = LSAME( RANGE, 'V' )
00165       INDEIG = LSAME( RANGE, 'I' )
00166 *
00167       INFO = 0
00168       IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN
00169          INFO = -1
00170       ELSE IF( .NOT.( ALLEIG .OR. VALEIG .OR. INDEIG ) ) THEN
00171          INFO = -2
00172       ELSE IF( N.LT.0 ) THEN
00173          INFO = -3
00174       ELSE
00175          IF( VALEIG ) THEN
00176             IF( N.GT.0 .AND. VU.LE.VL )
00177      $         INFO = -7
00178          ELSE IF( INDEIG ) THEN
00179             IF( IL.LT.1 .OR. IL.GT.MAX( 1, N ) ) THEN
00180                INFO = -8
00181             ELSE IF( IU.LT.MIN( N, IL ) .OR. IU.GT.N ) THEN
00182                INFO = -9
00183             END IF
00184          END IF
00185       END IF
00186       IF( INFO.EQ.0 ) THEN
00187          IF( LDZ.LT.1 .OR. ( WANTZ .AND. LDZ.LT.N ) )
00188      $      INFO = -14
00189       END IF
00190 *
00191       IF( INFO.NE.0 ) THEN
00192          CALL XERBLA( 'DSTEVX', -INFO )
00193          RETURN
00194       END IF
00195 *
00196 *     Quick return if possible
00197 *
00198       M = 0
00199       IF( N.EQ.0 )
00200      $   RETURN
00201 *
00202       IF( N.EQ.1 ) THEN
00203          IF( ALLEIG .OR. INDEIG ) THEN
00204             M = 1
00205             W( 1 ) = D( 1 )
00206          ELSE
00207             IF( VL.LT.D( 1 ) .AND. VU.GE.D( 1 ) ) THEN
00208                M = 1
00209                W( 1 ) = D( 1 )
00210             END IF
00211          END IF
00212          IF( WANTZ )
00213      $      Z( 1, 1 ) = ONE
00214          RETURN
00215       END IF
00216 *
00217 *     Get machine constants.
00218 *
00219       SAFMIN = DLAMCH( 'Safe minimum' )
00220       EPS = DLAMCH( 'Precision' )
00221       SMLNUM = SAFMIN / EPS
00222       BIGNUM = ONE / SMLNUM
00223       RMIN = SQRT( SMLNUM )
00224       RMAX = MIN( SQRT( BIGNUM ), ONE / SQRT( SQRT( SAFMIN ) ) )
00225 *
00226 *     Scale matrix to allowable range, if necessary.
00227 *
00228       ISCALE = 0
00229       IF( VALEIG ) THEN
00230          VLL = VL
00231          VUU = VU
00232       ELSE
00233          VLL = ZERO
00234          VUU = ZERO
00235       END IF
00236       TNRM = DLANST( 'M', N, D, E )
00237       IF( TNRM.GT.ZERO .AND. TNRM.LT.RMIN ) THEN
00238          ISCALE = 1
00239          SIGMA = RMIN / TNRM
00240       ELSE IF( TNRM.GT.RMAX ) THEN
00241          ISCALE = 1
00242          SIGMA = RMAX / TNRM
00243       END IF
00244       IF( ISCALE.EQ.1 ) THEN
00245          CALL DSCAL( N, SIGMA, D, 1 )
00246          CALL DSCAL( N-1, SIGMA, E( 1 ), 1 )
00247          IF( VALEIG ) THEN
00248             VLL = VL*SIGMA
00249             VUU = VU*SIGMA
00250          END IF
00251       END IF
00252 *
00253 *     If all eigenvalues are desired and ABSTOL is less than zero, then
00254 *     call DSTERF or SSTEQR.  If this fails for some eigenvalue, then
00255 *     try DSTEBZ.
00256 *
00257       TEST = .FALSE.
00258       IF( INDEIG ) THEN
00259          IF( IL.EQ.1 .AND. IU.EQ.N ) THEN
00260             TEST = .TRUE.
00261          END IF
00262       END IF
00263       IF( ( ALLEIG .OR. TEST ) .AND. ( ABSTOL.LE.ZERO ) ) THEN
00264          CALL DCOPY( N, D, 1, W, 1 )
00265          CALL DCOPY( N-1, E( 1 ), 1, WORK( 1 ), 1 )
00266          INDWRK = N + 1
00267          IF( .NOT.WANTZ ) THEN
00268             CALL DSTERF( N, W, WORK, INFO )
00269          ELSE
00270             CALL DSTEQR( 'I', N, W, WORK, Z, LDZ, WORK( INDWRK ), INFO )
00271             IF( INFO.EQ.0 ) THEN
00272                DO 10 I = 1, N
00273                   IFAIL( I ) = 0
00274    10          CONTINUE
00275             END IF
00276          END IF
00277          IF( INFO.EQ.0 ) THEN
00278             M = N
00279             GO TO 20
00280          END IF
00281          INFO = 0
00282       END IF
00283 *
00284 *     Otherwise, call DSTEBZ and, if eigenvectors are desired, SSTEIN.
00285 *
00286       IF( WANTZ ) THEN
00287          ORDER = 'B'
00288       ELSE
00289          ORDER = 'E'
00290       END IF
00291       INDWRK = 1
00292       INDIBL = 1
00293       INDISP = INDIBL + N
00294       INDIWO = INDISP + N
00295       CALL DSTEBZ( RANGE, ORDER, N, VLL, VUU, IL, IU, ABSTOL, D, E, M,
00296      $             NSPLIT, W, IWORK( INDIBL ), IWORK( INDISP ),
00297      $             WORK( INDWRK ), IWORK( INDIWO ), INFO )
00298 *
00299       IF( WANTZ ) THEN
00300          CALL DSTEIN( N, D, E, M, W, IWORK( INDIBL ), IWORK( INDISP ),
00301      $                Z, LDZ, WORK( INDWRK ), IWORK( INDIWO ), IFAIL,
00302      $                INFO )
00303       END IF
00304 *
00305 *     If matrix was scaled, then rescale eigenvalues appropriately.
00306 *
00307    20 CONTINUE
00308       IF( ISCALE.EQ.1 ) THEN
00309          IF( INFO.EQ.0 ) THEN
00310             IMAX = M
00311          ELSE
00312             IMAX = INFO - 1
00313          END IF
00314          CALL DSCAL( IMAX, ONE / SIGMA, W, 1 )
00315       END IF
00316 *
00317 *     If eigenvalues are not in order, then sort them, along with
00318 *     eigenvectors.
00319 *
00320       IF( WANTZ ) THEN
00321          DO 40 J = 1, M - 1
00322             I = 0
00323             TMP1 = W( J )
00324             DO 30 JJ = J + 1, M
00325                IF( W( JJ ).LT.TMP1 ) THEN
00326                   I = JJ
00327                   TMP1 = W( JJ )
00328                END IF
00329    30       CONTINUE
00330 *
00331             IF( I.NE.0 ) THEN
00332                ITMP1 = IWORK( INDIBL+I-1 )
00333                W( I ) = W( J )
00334                IWORK( INDIBL+I-1 ) = IWORK( INDIBL+J-1 )
00335                W( J ) = TMP1
00336                IWORK( INDIBL+J-1 ) = ITMP1
00337                CALL DSWAP( N, Z( 1, I ), 1, Z( 1, J ), 1 )
00338                IF( INFO.NE.0 ) THEN
00339                   ITMP1 = IFAIL( I )
00340                   IFAIL( I ) = IFAIL( J )
00341                   IFAIL( J ) = ITMP1
00342                END IF
00343             END IF
00344    40    CONTINUE
00345       END IF
00346 *
00347       RETURN
00348 *
00349 *     End of DSTEVX
00350 *
00351       END