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.
2: *
3: * =========== DOCUMENTATION ===========
4: *
5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
7: *
8: *> \htmlonly
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">
15: *> [TXT]</a>
16: *> \endhtmlonly
17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE DLARRR( N, D, E, INFO )
22: *
23: * .. Scalar Arguments ..
24: * INTEGER N, INFO
25: * ..
26: * .. Array Arguments ..
27: * DOUBLE PRECISION D( * ), E( * )
28: * ..
29: *
30: *
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: *
76: *> \author Univ. of Tennessee
77: *> \author Univ. of California Berkeley
78: *> \author Univ. of Colorado Denver
79: *> \author NAG Ltd.
80: *
81: *> \date September 2012
82: *
83: *> \ingroup auxOTHERauxiliary
84: *
85: *> \par Contributors:
86: * ==================
87: *>
88: *> Beresford Parlett, University of California, Berkeley, USA \n
89: *> Jim Demmel, University of California, Berkeley, USA \n
90: *> Inderjit Dhillon, University of Texas, Austin, USA \n
91: *> Osni Marques, LBNL/NERSC, USA \n
92: *> Christof Voemel, University of California, Berkeley, USA
93: *
94: * =====================================================================
95: SUBROUTINE DLARRR( N, D, E, INFO )
96: *
97: * -- LAPACK auxiliary routine (version 3.4.2) --
98: * -- LAPACK is a software package provided by Univ. of Tennessee, --
99: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
100: * September 2012
101: *
102: * .. Scalar Arguments ..
103: INTEGER N, INFO
104: * ..
105: * .. Array Arguments ..
106: DOUBLE PRECISION D( * ), E( * )
107: * ..
108: *
109: *
110: * =====================================================================
111: *
112: * .. Parameters ..
113: DOUBLE PRECISION ZERO, RELCOND
114: PARAMETER ( ZERO = 0.0D0,
115: $ RELCOND = 0.999D0 )
116: * ..
117: * .. Local Scalars ..
118: INTEGER I
119: LOGICAL YESREL
120: DOUBLE PRECISION EPS, SAFMIN, SMLNUM, RMIN, TMP, TMP2,
121: $ OFFDIG, OFFDIG2
122:
123: * ..
124: * .. External Functions ..
125: DOUBLE PRECISION DLAMCH
126: EXTERNAL DLAMCH
127: * ..
128: * .. Intrinsic Functions ..
129: INTRINSIC ABS
130: * ..
131: * .. Executable Statements ..
132: *
133: * As a default, do NOT go for relative-accuracy preserving computations.
134: INFO = 1
135:
136: SAFMIN = DLAMCH( 'Safe minimum' )
137: EPS = DLAMCH( 'Precision' )
138: SMLNUM = SAFMIN / EPS
139: RMIN = SQRT( SMLNUM )
140:
141: * Tests for relative accuracy
142: *
143: * Test for scaled diagonal dominance
144: * Scale the diagonal entries to one and check whether the sum of the
145: * off-diagonals is less than one
146: *
147: * The sdd relative error bounds have a 1/(1- 2*x) factor in them,
148: * x = max(OFFDIG + OFFDIG2), so when x is close to 1/2, no relative
149: * accuracy is promised. In the notation of the code fragment below,
150: * 1/(1 - (OFFDIG + OFFDIG2)) is the condition number.
151: * We don't think it is worth going into "sdd mode" unless the relative
152: * condition number is reasonable, not 1/macheps.
153: * The threshold should be compatible with other thresholds used in the
154: * code. We set OFFDIG + OFFDIG2 <= .999 =: RELCOND, it corresponds
155: * to losing at most 3 decimal digits: 1 / (1 - (OFFDIG + OFFDIG2)) <= 1000
156: * instead of the current OFFDIG + OFFDIG2 < 1
157: *
158: YESREL = .TRUE.
159: OFFDIG = ZERO
160: TMP = SQRT(ABS(D(1)))
161: IF (TMP.LT.RMIN) YESREL = .FALSE.
162: IF(.NOT.YESREL) GOTO 11
163: DO 10 I = 2, N
164: TMP2 = SQRT(ABS(D(I)))
165: IF (TMP2.LT.RMIN) YESREL = .FALSE.
166: IF(.NOT.YESREL) GOTO 11
167: OFFDIG2 = ABS(E(I-1))/(TMP*TMP2)
168: IF(OFFDIG+OFFDIG2.GE.RELCOND) YESREL = .FALSE.
169: IF(.NOT.YESREL) GOTO 11
170: TMP = TMP2
171: OFFDIG = OFFDIG2
172: 10 CONTINUE
173: 11 CONTINUE
174:
175: IF( YESREL ) THEN
176: INFO = 0
177: RETURN
178: ELSE
179: ENDIF
180: *
181:
182: *
183: * *** MORE TO BE IMPLEMENTED ***
184: *
185:
186: *
187: * Test if the lower bidiagonal matrix L from T = L D L^T
188: * (zero shift facto) is well conditioned
189: *
190:
191: *
192: * Test if the upper bidiagonal matrix U from T = U D U^T
193: * (zero shift facto) is well conditioned.
194: * In this case, the matrix needs to be flipped and, at the end
195: * of the eigenvector computation, the flip needs to be applied
196: * to the computed eigenvectors (and the support)
197: *
198:
199: *
200: RETURN
201: *
202: * END OF DLARRR
203: *
204: END
CVSweb interface <joel.bertrand@systella.fr>