LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
slarrk.f
Go to the documentation of this file.
00001 *> \brief \b SLARRK
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *> \htmlonly
00009 *> Download SLARRK + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slarrk.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slarrk.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slarrk.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       SUBROUTINE SLARRK( N, IW, GL, GU,
00022 *                           D, E2, PIVMIN, RELTOL, W, WERR, INFO)
00023 * 
00024 *       .. Scalar Arguments ..
00025 *       INTEGER   INFO, IW, N
00026 *       REAL                PIVMIN, RELTOL, GL, GU, W, WERR
00027 *       ..
00028 *       .. Array Arguments ..
00029 *       REAL               D( * ), E2( * )
00030 *       ..
00031 *  
00032 *
00033 *> \par Purpose:
00034 *  =============
00035 *>
00036 *> \verbatim
00037 *>
00038 *> SLARRK computes one eigenvalue of a symmetric tridiagonal
00039 *> matrix T to suitable accuracy. This is an auxiliary code to be
00040 *> called from SSTEMR.
00041 *>
00042 *> To avoid overflow, the matrix must be scaled so that its
00043 *> largest element is no greater than overflow**(1/2) * underflow**(1/4) in absolute value, and for greatest
00044 *> accuracy, it should not be much smaller than that.
00045 *>
00046 *> See W. Kahan "Accurate Eigenvalues of a Symmetric Tridiagonal
00047 *> Matrix", Report CS41, Computer Science Dept., Stanford
00048 *> University, July 21, 1966.
00049 *> \endverbatim
00050 *
00051 *  Arguments:
00052 *  ==========
00053 *
00054 *> \param[in] N
00055 *> \verbatim
00056 *>          N is INTEGER
00057 *>          The order of the tridiagonal matrix T.  N >= 0.
00058 *> \endverbatim
00059 *>
00060 *> \param[in] IW
00061 *> \verbatim
00062 *>          IW is INTEGER
00063 *>          The index of the eigenvalues to be returned.
00064 *> \endverbatim
00065 *>
00066 *> \param[in] GL
00067 *> \verbatim
00068 *>          GL is REAL
00069 *> \endverbatim
00070 *>
00071 *> \param[in] GU
00072 *> \verbatim
00073 *>          GU is REAL
00074 *>          An upper and a lower bound on the eigenvalue.
00075 *> \endverbatim
00076 *>
00077 *> \param[in] D
00078 *> \verbatim
00079 *>          D is REAL array, dimension (N)
00080 *>          The n diagonal elements of the tridiagonal matrix T.
00081 *> \endverbatim
00082 *>
00083 *> \param[in] E2
00084 *> \verbatim
00085 *>          E2 is REAL array, dimension (N-1)
00086 *>          The (n-1) squared off-diagonal elements of the tridiagonal matrix T.
00087 *> \endverbatim
00088 *>
00089 *> \param[in] PIVMIN
00090 *> \verbatim
00091 *>          PIVMIN is REAL
00092 *>          The minimum pivot allowed in the Sturm sequence for T.
00093 *> \endverbatim
00094 *>
00095 *> \param[in] RELTOL
00096 *> \verbatim
00097 *>          RELTOL is REAL
00098 *>          The minimum relative width of an interval.  When an interval
00099 *>          is narrower than RELTOL times the larger (in
00100 *>          magnitude) endpoint, then it is considered to be
00101 *>          sufficiently small, i.e., converged.  Note: this should
00102 *>          always be at least radix*machine epsilon.
00103 *> \endverbatim
00104 *>
00105 *> \param[out] W
00106 *> \verbatim
00107 *>          W is REAL
00108 *> \endverbatim
00109 *>
00110 *> \param[out] WERR
00111 *> \verbatim
00112 *>          WERR is REAL
00113 *>          The error bound on the corresponding eigenvalue approximation
00114 *>          in W.
00115 *> \endverbatim
00116 *>
00117 *> \param[out] INFO
00118 *> \verbatim
00119 *>          INFO is INTEGER
00120 *>          = 0:       Eigenvalue converged
00121 *>          = -1:      Eigenvalue did NOT converge
00122 *> \endverbatim
00123 *
00124 *> \par Internal Parameters:
00125 *  =========================
00126 *>
00127 *> \verbatim
00128 *>  FUDGE   REAL            , default = 2
00129 *>          A "fudge factor" to widen the Gershgorin intervals.
00130 *> \endverbatim
00131 *
00132 *  Authors:
00133 *  ========
00134 *
00135 *> \author Univ. of Tennessee 
00136 *> \author Univ. of California Berkeley 
00137 *> \author Univ. of Colorado Denver 
00138 *> \author NAG Ltd. 
00139 *
00140 *> \date November 2011
00141 *
00142 *> \ingroup auxOTHERauxiliary
00143 *
00144 *  =====================================================================
00145       SUBROUTINE SLARRK( N, IW, GL, GU,
00146      $                    D, E2, PIVMIN, RELTOL, W, WERR, INFO)
00147 *
00148 *  -- LAPACK auxiliary routine (version 3.4.0) --
00149 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00150 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00151 *     November 2011
00152 *
00153 *     .. Scalar Arguments ..
00154       INTEGER   INFO, IW, N
00155       REAL                PIVMIN, RELTOL, GL, GU, W, WERR
00156 *     ..
00157 *     .. Array Arguments ..
00158       REAL               D( * ), E2( * )
00159 *     ..
00160 *
00161 *  =====================================================================
00162 *
00163 *     .. Parameters ..
00164       REAL               FUDGE, HALF, TWO, ZERO
00165       PARAMETER          ( HALF = 0.5E0, TWO = 2.0E0,
00166      $                     FUDGE = TWO, ZERO = 0.0E0 )
00167 *     ..
00168 *     .. Local Scalars ..
00169       INTEGER   I, IT, ITMAX, NEGCNT
00170       REAL               ATOLI, EPS, LEFT, MID, RIGHT, RTOLI, TMP1,
00171      $                   TMP2, TNORM
00172 *     ..
00173 *     .. External Functions ..
00174       REAL               SLAMCH
00175       EXTERNAL   SLAMCH
00176 *     ..
00177 *     .. Intrinsic Functions ..
00178       INTRINSIC          ABS, INT, LOG, MAX
00179 *     ..
00180 *     .. Executable Statements ..
00181 *
00182 *     Get machine constants
00183       EPS = SLAMCH( 'P' )
00184 
00185       TNORM = MAX( ABS( GL ), ABS( GU ) )
00186       RTOLI = RELTOL
00187       ATOLI = FUDGE*TWO*PIVMIN
00188 
00189       ITMAX = INT( ( LOG( TNORM+PIVMIN )-LOG( PIVMIN ) ) /
00190      $           LOG( TWO ) ) + 2
00191 
00192       INFO = -1
00193 
00194       LEFT = GL - FUDGE*TNORM*EPS*N - FUDGE*TWO*PIVMIN
00195       RIGHT = GU + FUDGE*TNORM*EPS*N + FUDGE*TWO*PIVMIN
00196       IT = 0
00197 
00198  10   CONTINUE
00199 *
00200 *     Check if interval converged or maximum number of iterations reached
00201 *
00202       TMP1 = ABS( RIGHT - LEFT )
00203       TMP2 = MAX( ABS(RIGHT), ABS(LEFT) )
00204       IF( TMP1.LT.MAX( ATOLI, PIVMIN, RTOLI*TMP2 ) ) THEN
00205          INFO = 0
00206          GOTO 30
00207       ENDIF
00208       IF(IT.GT.ITMAX)
00209      $   GOTO 30
00210 
00211 *
00212 *     Count number of negative pivots for mid-point
00213 *
00214       IT = IT + 1
00215       MID = HALF * (LEFT + RIGHT)
00216       NEGCNT = 0
00217       TMP1 = D( 1 ) - MID
00218       IF( ABS( TMP1 ).LT.PIVMIN )
00219      $   TMP1 = -PIVMIN
00220       IF( TMP1.LE.ZERO )
00221      $   NEGCNT = NEGCNT + 1
00222 *
00223       DO 20 I = 2, N
00224          TMP1 = D( I ) - E2( I-1 ) / TMP1 - MID
00225          IF( ABS( TMP1 ).LT.PIVMIN )
00226      $      TMP1 = -PIVMIN
00227          IF( TMP1.LE.ZERO )
00228      $      NEGCNT = NEGCNT + 1
00229  20   CONTINUE
00230 
00231       IF(NEGCNT.GE.IW) THEN
00232          RIGHT = MID
00233       ELSE
00234          LEFT = MID
00235       ENDIF
00236       GOTO 10
00237 
00238  30   CONTINUE
00239 *
00240 *     Converged or maximum number of iterations reached
00241 *
00242       W = HALF * (LEFT + RIGHT)
00243       WERR = HALF * ABS( RIGHT - LEFT )
00244 
00245       RETURN
00246 *
00247 *     End of SLARRK
00248 *
00249       END
 All Files Functions