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