001:       SUBROUTINE SLARRR( N, D, E, INFO )
002: *
003: *  -- LAPACK auxiliary routine (version 3.2) --
004: *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
005: *     November 2006
006: *
007: *     .. Scalar Arguments ..
008:       INTEGER            N, INFO
009: *     ..
010: *     .. Array Arguments ..
011:       REAL               D( * ), E( * )
012: *     ..
013: *
014: *
015: *  Purpose
016: *  =======
017: *
018: *  Perform tests to decide whether the symmetric tridiagonal matrix T
019: *  warrants expensive computations which guarantee high relative accuracy
020: *  in the eigenvalues.
021: *
022: *  Arguments
023: *  =========
024: *
025: *  N       (input) INTEGER
026: *          The order of the matrix. N > 0.
027: *
028: *  D       (input) REAL             array, dimension (N)
029: *          The N diagonal elements of the tridiagonal matrix T.
030: *
031: *  E       (input/output) REAL             array, dimension (N)
032: *          On entry, the first (N-1) entries contain the subdiagonal
033: *          elements of the tridiagonal matrix T; E(N) is set to ZERO.
034: *
035: *  INFO    (output) INTEGER
036: *          INFO = 0(default) : the matrix warrants computations preserving
037: *                              relative accuracy.
038: *          INFO = 1          : the matrix warrants computations guaranteeing
039: *                              only absolute accuracy.
040: *
041: *  Further Details
042: *  ===============
043: *
044: *  Based on contributions by
045: *     Beresford Parlett, University of California, Berkeley, USA
046: *     Jim Demmel, University of California, Berkeley, USA
047: *     Inderjit Dhillon, University of Texas, Austin, USA
048: *     Osni Marques, LBNL/NERSC, USA
049: *     Christof Voemel, University of California, Berkeley, USA
050: *
051: *  =====================================================================
052: *
053: *     .. Parameters ..
054:       REAL               ZERO, RELCOND
055:       PARAMETER          ( ZERO = 0.0E0,
056:      $                     RELCOND = 0.999E0 )
057: *     ..
058: *     .. Local Scalars ..
059:       INTEGER            I
060:       LOGICAL            YESREL
061:       REAL               EPS, SAFMIN, SMLNUM, RMIN, TMP, TMP2,
062:      $          OFFDIG, OFFDIG2
063: 
064: *     ..
065: *     .. External Functions ..
066:       REAL               SLAMCH
067:       EXTERNAL           SLAMCH
068: *     ..
069: *     .. Intrinsic Functions ..
070:       INTRINSIC          ABS
071: *     ..
072: *     .. Executable Statements ..
073: *
074: *     As a default, do NOT go for relative-accuracy preserving computations.
075:       INFO = 1
076: 
077:       SAFMIN = SLAMCH( 'Safe minimum' )
078:       EPS = SLAMCH( 'Precision' )
079:       SMLNUM = SAFMIN / EPS
080:       RMIN = SQRT( SMLNUM )
081: 
082: *     Tests for relative accuracy
083: *
084: *     Test for scaled diagonal dominance
085: *     Scale the diagonal entries to one and check whether the sum of the
086: *     off-diagonals is less than one
087: *
088: *     The sdd relative error bounds have a 1/(1- 2*x) factor in them,
089: *     x = max(OFFDIG + OFFDIG2), so when x is close to 1/2, no relative
090: *     accuracy is promised.  In the notation of the code fragment below,
091: *     1/(1 - (OFFDIG + OFFDIG2)) is the condition number.
092: *     We don't think it is worth going into "sdd mode" unless the relative
093: *     condition number is reasonable, not 1/macheps.
094: *     The threshold should be compatible with other thresholds used in the
095: *     code. We set  OFFDIG + OFFDIG2 <= .999 =: RELCOND, it corresponds
096: *     to losing at most 3 decimal digits: 1 / (1 - (OFFDIG + OFFDIG2)) <= 1000
097: *     instead of the current OFFDIG + OFFDIG2 < 1
098: *
099:       YESREL = .TRUE.
100:       OFFDIG = ZERO
101:       TMP = SQRT(ABS(D(1)))
102:       IF (TMP.LT.RMIN) YESREL = .FALSE.
103:       IF(.NOT.YESREL) GOTO 11
104:       DO 10 I = 2, N
105:          TMP2 = SQRT(ABS(D(I)))
106:          IF (TMP2.LT.RMIN) YESREL = .FALSE.
107:          IF(.NOT.YESREL) GOTO 11
108:          OFFDIG2 = ABS(E(I-1))/(TMP*TMP2)
109:          IF(OFFDIG+OFFDIG2.GE.RELCOND) YESREL = .FALSE.
110:          IF(.NOT.YESREL) GOTO 11
111:          TMP = TMP2
112:          OFFDIG = OFFDIG2
113:  10   CONTINUE
114:  11   CONTINUE
115: 
116:       IF( YESREL ) THEN
117:          INFO = 0
118:          RETURN
119:       ELSE
120:       ENDIF
121: *
122: 
123: *
124: *     *** MORE TO BE IMPLEMENTED ***
125: *
126: 
127: *
128: *     Test if the lower bidiagonal matrix L from T = L D L^T
129: *     (zero shift facto) is well conditioned
130: *
131: 
132: *
133: *     Test if the upper bidiagonal matrix U from T = U D U^T
134: *     (zero shift facto) is well conditioned.
135: *     In this case, the matrix needs to be flipped and, at the end
136: *     of the eigenvector computation, the flip needs to be applied
137: *     to the computed eigenvectors (and the support)
138: *
139: 
140: *
141:       RETURN
142: *
143: *     END OF SLARRR
144: *
145:       END
146: