Annotation of rpl/lapack/lapack/dlarrr.f, revision 1.19

1.11      bertrand    1: *> \brief \b DLARRR performs tests to decide whether the symmetric tridiagonal matrix T warrants expensive computations which guarantee high relative accuracy in the eigenvalues.
1.8       bertrand    2: *
                      3: *  =========== DOCUMENTATION ===========
                      4: *
1.15      bertrand    5: * Online html documentation available at
                      6: *            http://www.netlib.org/lapack/explore-html/
1.8       bertrand    7: *
                      8: *> \htmlonly
1.15      bertrand    9: *> Download DLARRR + dependencies
                     10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlarrr.f">
                     11: *> [TGZ]</a>
                     12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlarrr.f">
                     13: *> [ZIP]</a>
                     14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlarrr.f">
1.8       bertrand   15: *> [TXT]</a>
1.15      bertrand   16: *> \endhtmlonly
1.8       bertrand   17: *
                     18: *  Definition:
                     19: *  ===========
                     20: *
                     21: *       SUBROUTINE DLARRR( N, D, E, INFO )
1.15      bertrand   22: *
1.8       bertrand   23: *       .. Scalar Arguments ..
                     24: *       INTEGER            N, INFO
                     25: *       ..
                     26: *       .. Array Arguments ..
                     27: *       DOUBLE PRECISION   D( * ), E( * )
                     28: *       ..
1.15      bertrand   29: *
                     30: *
1.8       bertrand   31: *
                     32: *> \par Purpose:
                     33: *  =============
                     34: *>
                     35: *> \verbatim
                     36: *>
                     37: *> Perform tests to decide whether the symmetric tridiagonal matrix T
                     38: *> warrants expensive computations which guarantee high relative accuracy
                     39: *> in the eigenvalues.
                     40: *> \endverbatim
                     41: *
                     42: *  Arguments:
                     43: *  ==========
                     44: *
                     45: *> \param[in] N
                     46: *> \verbatim
                     47: *>          N is INTEGER
                     48: *>          The order of the matrix. N > 0.
                     49: *> \endverbatim
                     50: *>
                     51: *> \param[in] D
                     52: *> \verbatim
                     53: *>          D is DOUBLE PRECISION array, dimension (N)
                     54: *>          The N diagonal elements of the tridiagonal matrix T.
                     55: *> \endverbatim
                     56: *>
                     57: *> \param[in,out] E
                     58: *> \verbatim
                     59: *>          E is DOUBLE PRECISION array, dimension (N)
                     60: *>          On entry, the first (N-1) entries contain the subdiagonal
                     61: *>          elements of the tridiagonal matrix T; E(N) is set to ZERO.
                     62: *> \endverbatim
                     63: *>
                     64: *> \param[out] INFO
                     65: *> \verbatim
                     66: *>          INFO is INTEGER
                     67: *>          INFO = 0(default) : the matrix warrants computations preserving
                     68: *>                              relative accuracy.
                     69: *>          INFO = 1          : the matrix warrants computations guaranteeing
                     70: *>                              only absolute accuracy.
                     71: *> \endverbatim
                     72: *
                     73: *  Authors:
                     74: *  ========
                     75: *
1.15      bertrand   76: *> \author Univ. of Tennessee
                     77: *> \author Univ. of California Berkeley
                     78: *> \author Univ. of Colorado Denver
                     79: *> \author NAG Ltd.
1.8       bertrand   80: *
1.15      bertrand   81: *> \ingroup OTHERauxiliary
1.8       bertrand   82: *
                     83: *> \par Contributors:
                     84: *  ==================
                     85: *>
                     86: *> Beresford Parlett, University of California, Berkeley, USA \n
                     87: *> Jim Demmel, University of California, Berkeley, USA \n
                     88: *> Inderjit Dhillon, University of Texas, Austin, USA \n
                     89: *> Osni Marques, LBNL/NERSC, USA \n
                     90: *> Christof Voemel, University of California, Berkeley, USA
                     91: *
                     92: *  =====================================================================
1.1       bertrand   93:       SUBROUTINE DLARRR( N, D, E, INFO )
                     94: *
1.19    ! bertrand   95: *  -- LAPACK auxiliary routine --
1.1       bertrand   96: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
                     97: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
                     98: *
                     99: *     .. Scalar Arguments ..
                    100:       INTEGER            N, INFO
                    101: *     ..
                    102: *     .. Array Arguments ..
                    103:       DOUBLE PRECISION   D( * ), E( * )
                    104: *     ..
                    105: *
                    106: *
                    107: *  =====================================================================
                    108: *
                    109: *     .. Parameters ..
                    110:       DOUBLE PRECISION   ZERO, RELCOND
                    111:       PARAMETER          ( ZERO = 0.0D0,
                    112:      $                     RELCOND = 0.999D0 )
                    113: *     ..
                    114: *     .. Local Scalars ..
                    115:       INTEGER            I
                    116:       LOGICAL            YESREL
                    117:       DOUBLE PRECISION   EPS, SAFMIN, SMLNUM, RMIN, TMP, TMP2,
                    118:      $          OFFDIG, OFFDIG2
                    119: 
                    120: *     ..
                    121: *     .. External Functions ..
                    122:       DOUBLE PRECISION   DLAMCH
                    123:       EXTERNAL           DLAMCH
                    124: *     ..
                    125: *     .. Intrinsic Functions ..
                    126:       INTRINSIC          ABS
                    127: *     ..
                    128: *     .. Executable Statements ..
                    129: *
1.17      bertrand  130: *     Quick return if possible
                    131: *
                    132:       IF( N.LE.0 ) THEN
                    133:          INFO = 0
                    134:          RETURN
                    135:       END IF
                    136: *
1.1       bertrand  137: *     As a default, do NOT go for relative-accuracy preserving computations.
                    138:       INFO = 1
                    139: 
                    140:       SAFMIN = DLAMCH( 'Safe minimum' )
                    141:       EPS = DLAMCH( 'Precision' )
                    142:       SMLNUM = SAFMIN / EPS
                    143:       RMIN = SQRT( SMLNUM )
                    144: 
                    145: *     Tests for relative accuracy
                    146: *
                    147: *     Test for scaled diagonal dominance
                    148: *     Scale the diagonal entries to one and check whether the sum of the
                    149: *     off-diagonals is less than one
                    150: *
                    151: *     The sdd relative error bounds have a 1/(1- 2*x) factor in them,
                    152: *     x = max(OFFDIG + OFFDIG2), so when x is close to 1/2, no relative
                    153: *     accuracy is promised.  In the notation of the code fragment below,
                    154: *     1/(1 - (OFFDIG + OFFDIG2)) is the condition number.
                    155: *     We don't think it is worth going into "sdd mode" unless the relative
                    156: *     condition number is reasonable, not 1/macheps.
                    157: *     The threshold should be compatible with other thresholds used in the
                    158: *     code. We set  OFFDIG + OFFDIG2 <= .999 =: RELCOND, it corresponds
                    159: *     to losing at most 3 decimal digits: 1 / (1 - (OFFDIG + OFFDIG2)) <= 1000
                    160: *     instead of the current OFFDIG + OFFDIG2 < 1
                    161: *
                    162:       YESREL = .TRUE.
                    163:       OFFDIG = ZERO
                    164:       TMP = SQRT(ABS(D(1)))
                    165:       IF (TMP.LT.RMIN) YESREL = .FALSE.
                    166:       IF(.NOT.YESREL) GOTO 11
                    167:       DO 10 I = 2, N
                    168:          TMP2 = SQRT(ABS(D(I)))
                    169:          IF (TMP2.LT.RMIN) YESREL = .FALSE.
                    170:          IF(.NOT.YESREL) GOTO 11
                    171:          OFFDIG2 = ABS(E(I-1))/(TMP*TMP2)
                    172:          IF(OFFDIG+OFFDIG2.GE.RELCOND) YESREL = .FALSE.
                    173:          IF(.NOT.YESREL) GOTO 11
                    174:          TMP = TMP2
                    175:          OFFDIG = OFFDIG2
                    176:  10   CONTINUE
                    177:  11   CONTINUE
                    178: 
                    179:       IF( YESREL ) THEN
                    180:          INFO = 0
                    181:          RETURN
                    182:       ELSE
                    183:       ENDIF
                    184: *
                    185: 
                    186: *
                    187: *     *** MORE TO BE IMPLEMENTED ***
                    188: *
                    189: 
                    190: *
                    191: *     Test if the lower bidiagonal matrix L from T = L D L^T
                    192: *     (zero shift facto) is well conditioned
                    193: *
                    194: 
                    195: *
                    196: *     Test if the upper bidiagonal matrix U from T = U D U^T
                    197: *     (zero shift facto) is well conditioned.
                    198: *     In this case, the matrix needs to be flipped and, at the end
                    199: *     of the eigenvector computation, the flip needs to be applied
                    200: *     to the computed eigenvectors (and the support)
                    201: *
                    202: 
                    203: *
                    204:       RETURN
                    205: *
1.19    ! bertrand  206: *     End of DLARRR
1.1       bertrand  207: *
                    208:       END

CVSweb interface <joel.bertrand@systella.fr>