SRC/slarrk.f

   1 *> \brief \b SLARRK computes one eigenvalue of a symmetric tridiagonal matrix T to suitable accuracy.
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download SLARRK + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slarrk.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slarrk.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slarrk.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE SLARRK( N, IW, GL, GU,
  22 *                           D, E2, PIVMIN, RELTOL, W, WERR, INFO)
  23 *
  24 *       .. Scalar Arguments ..
  25 *       INTEGER   INFO, IW, N
  26 *       REAL                PIVMIN, RELTOL, GL, GU, W, WERR
  27 *       ..
  28 *       .. Array Arguments ..
  29 *       REAL               D( * ), E2( * )
  30 *       ..
  31 *
  32 *
  33 *> \par Purpose:
  34 *  =============
  35 *>
  36 *> \verbatim
  37 *>
  38 *> SLARRK computes one eigenvalue of a symmetric tridiagonal
  39 *> matrix T to suitable accuracy. This is an auxiliary code to be
  40 *> called from SSTEMR.
  41 *>
  42 *> To avoid overflow, the matrix must be scaled so that its
  43 *> largest element is no greater than overflow**(1/2) * underflow**(1/4) in absolute value, and for greatest
  44 *> accuracy, it should not be much smaller than that.
  45 *>
  46 *> See W. Kahan "Accurate Eigenvalues of a Symmetric Tridiagonal
  47 *> Matrix", Report CS41, Computer Science Dept., Stanford
  48 *> University, July 21, 1966.
  49 *> \endverbatim
  50 *
  51 *  Arguments:
  52 *  ==========
  53 *
  54 *> \param[in] N
  55 *> \verbatim
  56 *>          N is INTEGER
  57 *>          The order of the tridiagonal matrix T.  N >= 0.
  58 *> \endverbatim
  59 *>
  60 *> \param[in] IW
  61 *> \verbatim
  62 *>          IW is INTEGER
  63 *>          The index of the eigenvalues to be returned.
  64 *> \endverbatim
  65 *>
  66 *> \param[in] GL
  67 *> \verbatim
  68 *>          GL is REAL
  69 *> \endverbatim
  70 *>
  71 *> \param[in] GU
  72 *> \verbatim
  73 *>          GU is REAL
  74 *>          An upper and a lower bound on the eigenvalue.
  75 *> \endverbatim
  76 *>
  77 *> \param[in] D
  78 *> \verbatim
  79 *>          D is REAL array, dimension (N)
  80 *>          The n diagonal elements of the tridiagonal matrix T.
  81 *> \endverbatim
  82 *>
  83 *> \param[in] E2
  84 *> \verbatim
  85 *>          E2 is REAL array, dimension (N-1)
  86 *>          The (n-1) squared off-diagonal elements of the tridiagonal matrix T.
  87 *> \endverbatim
  88 *>
  89 *> \param[in] PIVMIN
  90 *> \verbatim
  91 *>          PIVMIN is REAL
  92 *>          The minimum pivot allowed in the Sturm sequence for T.
  93 *> \endverbatim
  94 *>
  95 *> \param[in] RELTOL
  96 *> \verbatim
  97 *>          RELTOL is REAL
  98 *>          The minimum relative width of an interval.  When an interval
  99 *>          is narrower than RELTOL times the larger (in
 100 *>          magnitude) endpoint, then it is considered to be
 101 *>          sufficiently small, i.e., converged.  Note: this should
 102 *>          always be at least radix*machine epsilon.
 103 *> \endverbatim
 104 *>
 105 *> \param[out] W
 106 *> \verbatim
 107 *>          W is REAL
 108 *> \endverbatim
 109 *>
 110 *> \param[out] WERR
 111 *> \verbatim
 112 *>          WERR is REAL
 113 *>          The error bound on the corresponding eigenvalue approximation
 114 *>          in W.
 115 *> \endverbatim
 116 *>
 117 *> \param[out] INFO
 118 *> \verbatim
 119 *>          INFO is INTEGER
 120 *>          = 0:       Eigenvalue converged
 121 *>          = -1:      Eigenvalue did NOT converge
 122 *> \endverbatim
 123 *
 124 *> \par Internal Parameters:
 125 *  =========================
 126 *>
 127 *> \verbatim
 128 *>  FUDGE   REAL            , default = 2
 129 *>          A "fudge factor" to widen the Gershgorin intervals.
 130 *> \endverbatim
 131 *
 132 *  Authors:
 133 *  ========
 134 *
 135 *> \author Univ. of Tennessee
 136 *> \author Univ. of California Berkeley
 137 *> \author Univ. of Colorado Denver
 138 *> \author NAG Ltd.
 139 *
 140 *> \date September 2012
 141 *
 142 *> \ingroup OTHERauxiliary
 143 *
 144 *  =====================================================================
 145       SUBROUTINE SLARRK( N, IW, GL, GU,
 146      $                    D, E2, PIVMIN, RELTOL, W, WERR, INFO)
 147 *
 148 *  -- LAPACK auxiliary routine (version 3.4.2) --
 149 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 150 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 151 *     September 2012
 152 *
 153 *     .. Scalar Arguments ..
 154       INTEGER   INFO, IW, N
 155       REAL                PIVMIN, RELTOL, GL, GU, W, WERR
 156 *     ..
 157 *     .. Array Arguments ..
 158       REAL               D( * ), E2( * )
 159 *     ..
 160 *
 161 *  =====================================================================
 162 *
 163 *     .. Parameters ..
 164       REAL               FUDGE, HALF, TWO, ZERO
 165       PARAMETER          ( HALF = 0.5E0, TWO = 2.0E0,
 166      $                     FUDGE = TWO, ZERO = 0.0E0 )
 167 *     ..
 168 *     .. Local Scalars ..
 169       INTEGER   I, IT, ITMAX, NEGCNT
 170       REAL               ATOLI, EPS, LEFT, MID, RIGHT, RTOLI, TMP1,
 171      $                   TMP2, TNORM
 172 *     ..
 173 *     .. External Functions ..
 174       REAL               SLAMCH
 175       EXTERNAL   SLAMCH
 176 *     ..
 177 *     .. Intrinsic Functions ..
 178       INTRINSIC          ABS, INT, LOG, MAX
 179 *     ..
 180 *     .. Executable Statements ..
 181 *
 182 *     Get machine constants
 183       EPS = SLAMCH( 'P' )
 184
 185       TNORM = MAX( ABS( GL ), ABS( GU ) )
 186       RTOLI = RELTOL
 187       ATOLI = FUDGE*TWO*PIVMIN
 188
 189       ITMAX = INT( ( LOG( TNORM+PIVMIN )-LOG( PIVMIN ) ) /
 190      $           LOG( TWO ) ) + 2
 191
 192       INFO = -1
 193
 194       LEFT = GL - FUDGE*TNORM*EPS*N - FUDGE*TWO*PIVMIN
 195       RIGHT = GU + FUDGE*TNORM*EPS*N + FUDGE*TWO*PIVMIN
 196       IT = 0
 197
 198  10   CONTINUE
 199 *
 200 *     Check if interval converged or maximum number of iterations reached
 201 *
 202       TMP1 = ABS( RIGHT - LEFT )
 203       TMP2 = MAX( ABS(RIGHT), ABS(LEFT) )
 204       IF( TMP1.LT.MAX( ATOLI, PIVMIN, RTOLI*TMP2 ) ) THEN
 205          INFO = 0
 206          GOTO 30
 207       ENDIF
 208       IF(IT.GT.ITMAX)
 209      $   GOTO 30
 210
 211 *
 212 *     Count number of negative pivots for mid-point
 213 *
 214       IT = IT + 1
 215       MID = HALF * (LEFT + RIGHT)
 216       NEGCNT = 0
 217       TMP1 = D( 1 ) - MID
 218       IF( ABS( TMP1 ).LT.PIVMIN )
 219      $   TMP1 = -PIVMIN
 220       IF( TMP1.LE.ZERO )
 221      $   NEGCNT = NEGCNT + 1
 222 *
 223       DO 20 I = 2, N
 224          TMP1 = D( I ) - E2( I-1 ) / TMP1 - MID
 225          IF( ABS( TMP1 ).LT.PIVMIN )
 226      $      TMP1 = -PIVMIN
 227          IF( TMP1.LE.ZERO )
 228      $      NEGCNT = NEGCNT + 1
 229  20   CONTINUE
 230
 231       IF(NEGCNT.GE.IW) THEN
 232          RIGHT = MID
 233       ELSE
 234          LEFT = MID
 235       ENDIF
 236       GOTO 10
 237
 238  30   CONTINUE
 239 *
 240 *     Converged or maximum number of iterations reached
 241 *
 242       W = HALF * (LEFT + RIGHT)
 243       WERR = HALF * ABS( RIGHT - LEFT )
 244
 245       RETURN
 246 *
 247 *     End of SLARRK
 248 *
 249       END