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
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