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
CVSweb interface <joel.bertrand@systella.fr>