001:       SUBROUTINE SLARRR( N, D, E, INFO )
002: *
003: *  -- LAPACK auxiliary routine (version 3.2) --
004: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
005: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
006: *     November 2006
007: *
008: *     .. Scalar Arguments ..
009:       INTEGER            N, INFO
010: *     ..
011: *     .. Array Arguments ..
012:       REAL               D( * ), E( * )
013: *     ..
014: *
015: *
016: *  Purpose
017: *  =======
018: *
019: *  Perform tests to decide whether the symmetric tridiagonal matrix T
020: *  warrants expensive computations which guarantee high relative accuracy
021: *  in the eigenvalues.
022: *
023: *  Arguments
024: *  =========
025: *
026: *  N       (input) INTEGER
027: *          The order of the matrix. N > 0.
028: *
029: *  D       (input) REAL             array, dimension (N)
030: *          The N diagonal elements of the tridiagonal matrix T.
031: *
032: *  E       (input/output) REAL             array, dimension (N)
033: *          On entry, the first (N-1) entries contain the subdiagonal
034: *          elements of the tridiagonal matrix T; E(N) is set to ZERO.
035: *
036: *  INFO    (output) INTEGER
037: *          INFO = 0(default) : the matrix warrants computations preserving
038: *                              relative accuracy.
039: *          INFO = 1          : the matrix warrants computations guaranteeing
040: *                              only absolute accuracy.
041: *
042: *  Further Details
043: *  ===============
044: *
045: *  Based on contributions by
046: *     Beresford Parlett, University of California, Berkeley, USA
047: *     Jim Demmel, University of California, Berkeley, USA
048: *     Inderjit Dhillon, University of Texas, Austin, USA
049: *     Osni Marques, LBNL/NERSC, USA
050: *     Christof Voemel, University of California, Berkeley, USA
051: *
052: *  =====================================================================
053: *
054: *     .. Parameters ..
055:       REAL               ZERO, RELCOND
056:       PARAMETER          ( ZERO = 0.0E0,
057:      $                     RELCOND = 0.999E0 )
058: *     ..
059: *     .. Local Scalars ..
060:       INTEGER            I
061:       LOGICAL            YESREL
062:       REAL               EPS, SAFMIN, SMLNUM, RMIN, TMP, TMP2,
063:      $          OFFDIG, OFFDIG2
064: 
065: *     ..
066: *     .. External Functions ..
067:       REAL               SLAMCH
068:       EXTERNAL           SLAMCH
069: *     ..
070: *     .. Intrinsic Functions ..
071:       INTRINSIC          ABS
072: *     ..
073: *     .. Executable Statements ..
074: *
075: *     As a default, do NOT go for relative-accuracy preserving computations.
076:       INFO = 1
077: 
078:       SAFMIN = SLAMCH( 'Safe minimum' )
079:       EPS = SLAMCH( 'Precision' )
080:       SMLNUM = SAFMIN / EPS
081:       RMIN = SQRT( SMLNUM )
082: 
083: *     Tests for relative accuracy
084: *
085: *     Test for scaled diagonal dominance
086: *     Scale the diagonal entries to one and check whether the sum of the
087: *     off-diagonals is less than one
088: *
089: *     The sdd relative error bounds have a 1/(1- 2*x) factor in them,
090: *     x = max(OFFDIG + OFFDIG2), so when x is close to 1/2, no relative
091: *     accuracy is promised.  In the notation of the code fragment below,
092: *     1/(1 - (OFFDIG + OFFDIG2)) is the condition number.
093: *     We don't think it is worth going into "sdd mode" unless the relative
094: *     condition number is reasonable, not 1/macheps.
095: *     The threshold should be compatible with other thresholds used in the
096: *     code. We set  OFFDIG + OFFDIG2 <= .999 =: RELCOND, it corresponds
097: *     to losing at most 3 decimal digits: 1 / (1 - (OFFDIG + OFFDIG2)) <= 1000
098: *     instead of the current OFFDIG + OFFDIG2 < 1
099: *
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 SLARRR
145: *
146:       END
147: