Return to dlarrr.f CVS log

Up to [local] / rpl / lapack / lapack

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

1.1     ! bertrand    1:       SUBROUTINE DLARRR( N, D, E, INFO )
        !             2: *
        !             3: *  -- LAPACK auxiliary routine (version 3.2) --
        !             4: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
        !             5: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
        !             6: *     November 2006
        !             7: *
        !             8: *     .. Scalar Arguments ..
        !             9:       INTEGER            N, INFO
        !            10: *     ..
        !            11: *     .. Array Arguments ..
        !            12:       DOUBLE PRECISION   D( * ), E( * )
        !            13: *     ..
        !            14: *
        !            15: *
        !            16: *  Purpose
        !            17: *  =======
        !            18: *
        !            19: *  Perform tests to decide whether the symmetric tridiagonal matrix T
        !            20: *  warrants expensive computations which guarantee high relative accuracy
        !            21: *  in the eigenvalues.
        !            22: *
        !            23: *  Arguments
        !            24: *  =========
        !            25: *
        !            26: *  N       (input) INTEGER
        !            27: *          The order of the matrix. N > 0.
        !            28: *
        !            29: *  D       (input) DOUBLE PRECISION array, dimension (N)
        !            30: *          The N diagonal elements of the tridiagonal matrix T.
        !            31: *
        !            32: *  E       (input/output) DOUBLE PRECISION array, dimension (N)
        !            33: *          On entry, the first (N-1) entries contain the subdiagonal
        !            34: *          elements of the tridiagonal matrix T; E(N) is set to ZERO.
        !            35: *
        !            36: *  INFO    (output) INTEGER
        !            37: *          INFO = 0(default) : the matrix warrants computations preserving
        !            38: *                              relative accuracy.
        !            39: *          INFO = 1          : the matrix warrants computations guaranteeing
        !            40: *                              only absolute accuracy.
        !            41: *
        !            42: *  Further Details
        !            43: *  ===============
        !            44: *
        !            45: *  Based on contributions by
        !            46: *     Beresford Parlett, University of California, Berkeley, USA
        !            47: *     Jim Demmel, University of California, Berkeley, USA
        !            48: *     Inderjit Dhillon, University of Texas, Austin, USA
        !            49: *     Osni Marques, LBNL/NERSC, USA
        !            50: *     Christof Voemel, University of California, Berkeley, USA
        !            51: *
        !            52: *  =====================================================================
        !            53: *
        !            54: *     .. Parameters ..
        !            55:       DOUBLE PRECISION   ZERO, RELCOND
        !            56:       PARAMETER          ( ZERO = 0.0D0,
        !            57:      $                     RELCOND = 0.999D0 )
        !            58: *     ..
        !            59: *     .. Local Scalars ..
        !            60:       INTEGER            I
        !            61:       LOGICAL            YESREL
        !            62:       DOUBLE PRECISION   EPS, SAFMIN, SMLNUM, RMIN, TMP, TMP2,
        !            63:      $          OFFDIG, OFFDIG2
        !            64: 
        !            65: *     ..
        !            66: *     .. External Functions ..
        !            67:       DOUBLE PRECISION   DLAMCH
        !            68:       EXTERNAL           DLAMCH
        !            69: *     ..
        !            70: *     .. Intrinsic Functions ..
        !            71:       INTRINSIC          ABS
        !            72: *     ..
        !            73: *     .. Executable Statements ..
        !            74: *
        !            75: *     As a default, do NOT go for relative-accuracy preserving computations.
        !            76:       INFO = 1
        !            77: 
        !            78:       SAFMIN = DLAMCH( 'Safe minimum' )
        !            79:       EPS = DLAMCH( 'Precision' )
        !            80:       SMLNUM = SAFMIN / EPS
        !            81:       RMIN = SQRT( SMLNUM )
        !            82: 
        !            83: *     Tests for relative accuracy
        !            84: *
        !            85: *     Test for scaled diagonal dominance
        !            86: *     Scale the diagonal entries to one and check whether the sum of the
        !            87: *     off-diagonals is less than one
        !            88: *
        !            89: *     The sdd relative error bounds have a 1/(1- 2*x) factor in them,
        !            90: *     x = max(OFFDIG + OFFDIG2), so when x is close to 1/2, no relative
        !            91: *     accuracy is promised.  In the notation of the code fragment below,
        !            92: *     1/(1 - (OFFDIG + OFFDIG2)) is the condition number.
        !            93: *     We don't think it is worth going into "sdd mode" unless the relative
        !            94: *     condition number is reasonable, not 1/macheps.
        !            95: *     The threshold should be compatible with other thresholds used in the
        !            96: *     code. We set  OFFDIG + OFFDIG2 <= .999 =: RELCOND, it corresponds
        !            97: *     to losing at most 3 decimal digits: 1 / (1 - (OFFDIG + OFFDIG2)) <= 1000
        !            98: *     instead of the current OFFDIG + OFFDIG2 < 1
        !            99: *
        !           100:       YESREL = .TRUE.
        !           101:       OFFDIG = ZERO
        !           102:       TMP = SQRT(ABS(D(1)))
        !           103:       IF (TMP.LT.RMIN) YESREL = .FALSE.
        !           104:       IF(.NOT.YESREL) GOTO 11
        !           105:       DO 10 I = 2, N
        !           106:          TMP2 = SQRT(ABS(D(I)))
        !           107:          IF (TMP2.LT.RMIN) YESREL = .FALSE.
        !           108:          IF(.NOT.YESREL) GOTO 11
        !           109:          OFFDIG2 = ABS(E(I-1))/(TMP*TMP2)
        !           110:          IF(OFFDIG+OFFDIG2.GE.RELCOND) YESREL = .FALSE.
        !           111:          IF(.NOT.YESREL) GOTO 11
        !           112:          TMP = TMP2
        !           113:          OFFDIG = OFFDIG2
        !           114:  10   CONTINUE
        !           115:  11   CONTINUE
        !           116: 
        !           117:       IF( YESREL ) THEN
        !           118:          INFO = 0
        !           119:          RETURN
        !           120:       ELSE
        !           121:       ENDIF
        !           122: *
        !           123: 
        !           124: *
        !           125: *     *** MORE TO BE IMPLEMENTED ***
        !           126: *
        !           127: 
        !           128: *
        !           129: *     Test if the lower bidiagonal matrix L from T = L D L^T
        !           130: *     (zero shift facto) is well conditioned
        !           131: *
        !           132: 
        !           133: *
        !           134: *     Test if the upper bidiagonal matrix U from T = U D U^T
        !           135: *     (zero shift facto) is well conditioned.
        !           136: *     In this case, the matrix needs to be flipped and, at the end
        !           137: *     of the eigenvector computation, the flip needs to be applied
        !           138: *     to the computed eigenvectors (and the support)
        !           139: *
        !           140: 
        !           141: *
        !           142:       RETURN
        !           143: *
        !           144: *     END OF DLARRR
        !           145: *
        !           146:       END