Annotation of rpl/lapack/lapack/dlarrk.f, revision 1.8

1.8     ! bertrand    1: *> \brief \b DLARRK
        !             2: *
        !             3: *  =========== DOCUMENTATION ===========
        !             4: *
        !             5: * Online html documentation available at 
        !             6: *            http://www.netlib.org/lapack/explore-html/ 
        !             7: *
        !             8: *> \htmlonly
        !             9: *> Download DLARRK + dependencies 
        !            10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlarrk.f"> 
        !            11: *> [TGZ]</a> 
        !            12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlarrk.f"> 
        !            13: *> [ZIP]</a> 
        !            14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlarrk.f"> 
        !            15: *> [TXT]</a>
        !            16: *> \endhtmlonly 
        !            17: *
        !            18: *  Definition:
        !            19: *  ===========
        !            20: *
        !            21: *       SUBROUTINE DLARRK( N, IW, GL, GU,
        !            22: *                           D, E2, PIVMIN, RELTOL, W, WERR, INFO)
        !            23: * 
        !            24: *       .. Scalar Arguments ..
        !            25: *       INTEGER   INFO, IW, N
        !            26: *       DOUBLE PRECISION    PIVMIN, RELTOL, GL, GU, W, WERR
        !            27: *       ..
        !            28: *       .. Array Arguments ..
        !            29: *       DOUBLE PRECISION   D( * ), E2( * )
        !            30: *       ..
        !            31: *  
        !            32: *
        !            33: *> \par Purpose:
        !            34: *  =============
        !            35: *>
        !            36: *> \verbatim
        !            37: *>
        !            38: *> DLARRK computes one eigenvalue of a symmetric tridiagonal
        !            39: *> matrix T to suitable accuracy. This is an auxiliary code to be
        !            40: *> called from DSTEMR.
        !            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 DOUBLE PRECISION
        !            69: *> \endverbatim
        !            70: *>
        !            71: *> \param[in] GU
        !            72: *> \verbatim
        !            73: *>          GU is DOUBLE PRECISION
        !            74: *>          An upper and a lower bound on the eigenvalue.
        !            75: *> \endverbatim
        !            76: *>
        !            77: *> \param[in] D
        !            78: *> \verbatim
        !            79: *>          D is DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION
        !            92: *>          The minimum pivot allowed in the Sturm sequence for T.
        !            93: *> \endverbatim
        !            94: *>
        !            95: *> \param[in] RELTOL
        !            96: *> \verbatim
        !            97: *>          RELTOL is DOUBLE PRECISION
        !            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 DOUBLE PRECISION
        !           108: *> \endverbatim
        !           109: *>
        !           110: *> \param[out] WERR
        !           111: *> \verbatim
        !           112: *>          WERR is DOUBLE PRECISION
        !           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   DOUBLE PRECISION, 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 November 2011
        !           141: *
        !           142: *> \ingroup auxOTHERauxiliary
        !           143: *
        !           144: *  =====================================================================
1.1       bertrand  145:       SUBROUTINE DLARRK( N, IW, GL, GU,
                    146:      $                    D, E2, PIVMIN, RELTOL, W, WERR, INFO)
                    147: *
1.8     ! bertrand  148: *  -- LAPACK auxiliary routine (version 3.4.0) --
1.1       bertrand  149: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
                    150: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.8     ! bertrand  151: *     November 2011
1.1       bertrand  152: *
                    153: *     .. Scalar Arguments ..
                    154:       INTEGER   INFO, IW, N
                    155:       DOUBLE PRECISION    PIVMIN, RELTOL, GL, GU, W, WERR
                    156: *     ..
                    157: *     .. Array Arguments ..
                    158:       DOUBLE PRECISION   D( * ), E2( * )
                    159: *     ..
                    160: *
                    161: *  =====================================================================
                    162: *
                    163: *     .. Parameters ..
                    164:       DOUBLE PRECISION   FUDGE, HALF, TWO, ZERO
                    165:       PARAMETER          ( HALF = 0.5D0, TWO = 2.0D0,
                    166:      $                     FUDGE = TWO, ZERO = 0.0D0 )
                    167: *     ..
                    168: *     .. Local Scalars ..
                    169:       INTEGER   I, IT, ITMAX, NEGCNT
                    170:       DOUBLE PRECISION   ATOLI, EPS, LEFT, MID, RIGHT, RTOLI, TMP1,
                    171:      $                   TMP2, TNORM
                    172: *     ..
                    173: *     .. External Functions ..
                    174:       DOUBLE PRECISION   DLAMCH
                    175:       EXTERNAL   DLAMCH
                    176: *     ..
                    177: *     .. Intrinsic Functions ..
                    178:       INTRINSIC          ABS, INT, LOG, MAX
                    179: *     ..
                    180: *     .. Executable Statements ..
                    181: *
                    182: *     Get machine constants
                    183:       EPS = DLAMCH( '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 DLARRK
                    248: *
                    249:       END

CVSweb interface <joel.bertrand@systella.fr>