1 SUBROUTINE DLARRK( N, IW, GL, GU,
2 $ D, E2, PIVMIN, RELTOL, W, WERR, INFO)
3 IMPLICIT NONE
4 *
5 * -- LAPACK auxiliary routine (version 3.2) --
6 * -- LAPACK is a software package provided by Univ. of Tennessee, --
7 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
8 * November 2006
9 *
10 * .. Scalar Arguments ..
11 INTEGER INFO, IW, N
12 DOUBLE PRECISION PIVMIN, RELTOL, GL, GU, W, WERR
13 * ..
14 * .. Array Arguments ..
15 DOUBLE PRECISION D( * ), E2( * )
16 * ..
17 *
18 * Purpose
19 * =======
20 *
21 * DLARRK computes one eigenvalue of a symmetric tridiagonal
22 * matrix T to suitable accuracy. This is an auxiliary code to be
23 * called from DSTEMR.
24 *
25 * To avoid overflow, the matrix must be scaled so that its
26 * largest element is no greater than overflow**(1/2) *
27 * underflow**(1/4) in absolute value, and for greatest
28 * accuracy, it should not be much smaller than that.
29 *
30 * See W. Kahan "Accurate Eigenvalues of a Symmetric Tridiagonal
31 * Matrix", Report CS41, Computer Science Dept., Stanford
32 * University, July 21, 1966.
33 *
34 * Arguments
35 * =========
36 *
37 * N (input) INTEGER
38 * The order of the tridiagonal matrix T. N >= 0.
39 *
40 * IW (input) INTEGER
41 * The index of the eigenvalues to be returned.
42 *
43 * GL (input) DOUBLE PRECISION
44 * GU (input) DOUBLE PRECISION
45 * An upper and a lower bound on the eigenvalue.
46 *
47 * D (input) DOUBLE PRECISION array, dimension (N)
48 * The n diagonal elements of the tridiagonal matrix T.
49 *
50 * E2 (input) DOUBLE PRECISION array, dimension (N-1)
51 * The (n-1) squared off-diagonal elements of the tridiagonal matrix T.
52 *
53 * PIVMIN (input) DOUBLE PRECISION
54 * The minimum pivot allowed in the Sturm sequence for T.
55 *
56 * RELTOL (input) DOUBLE PRECISION
57 * The minimum relative width of an interval. When an interval
58 * is narrower than RELTOL times the larger (in
59 * magnitude) endpoint, then it is considered to be
60 * sufficiently small, i.e., converged. Note: this should
61 * always be at least radix*machine epsilon.
62 *
63 * W (output) DOUBLE PRECISION
64 *
65 * WERR (output) DOUBLE PRECISION
66 * The error bound on the corresponding eigenvalue approximation
67 * in W.
68 *
69 * INFO (output) INTEGER
70 * = 0: Eigenvalue converged
71 * = -1: Eigenvalue did NOT converge
72 *
73 * Internal Parameters
74 * ===================
75 *
76 * FUDGE DOUBLE PRECISION, default = 2
77 * A "fudge factor" to widen the Gershgorin intervals.
78 *
79 * =====================================================================
80 *
81 * .. Parameters ..
82 DOUBLE PRECISION FUDGE, HALF, TWO, ZERO
83 PARAMETER ( HALF = 0.5D0, TWO = 2.0D0,
84 $ FUDGE = TWO, ZERO = 0.0D0 )
85 * ..
86 * .. Local Scalars ..
87 INTEGER I, IT, ITMAX, NEGCNT
88 DOUBLE PRECISION ATOLI, EPS, LEFT, MID, RIGHT, RTOLI, TMP1,
89 $ TMP2, TNORM
90 * ..
91 * .. External Functions ..
92 DOUBLE PRECISION DLAMCH
93 EXTERNAL DLAMCH
94 * ..
95 * .. Intrinsic Functions ..
96 INTRINSIC ABS, INT, LOG, MAX
97 * ..
98 * .. Executable Statements ..
99 *
100 * Get machine constants
101 EPS = DLAMCH( 'P' )
102
103 TNORM = MAX( ABS( GL ), ABS( GU ) )
104 RTOLI = RELTOL
105 ATOLI = FUDGE*TWO*PIVMIN
106
107 ITMAX = INT( ( LOG( TNORM+PIVMIN )-LOG( PIVMIN ) ) /
108 $ LOG( TWO ) ) + 2
109
110 INFO = -1
111
112 LEFT = GL - FUDGE*TNORM*EPS*N - FUDGE*TWO*PIVMIN
113 RIGHT = GU + FUDGE*TNORM*EPS*N + FUDGE*TWO*PIVMIN
114 IT = 0
115
116 10 CONTINUE
117 *
118 * Check if interval converged or maximum number of iterations reached
119 *
120 TMP1 = ABS( RIGHT - LEFT )
121 TMP2 = MAX( ABS(RIGHT), ABS(LEFT) )
122 IF( TMP1.LT.MAX( ATOLI, PIVMIN, RTOLI*TMP2 ) ) THEN
123 INFO = 0
124 GOTO 30
125 ENDIF
126 IF(IT.GT.ITMAX)
127 $ GOTO 30
128
129 *
130 * Count number of negative pivots for mid-point
131 *
132 IT = IT + 1
133 MID = HALF * (LEFT + RIGHT)
134 NEGCNT = 0
135 TMP1 = D( 1 ) - MID
136 IF( ABS( TMP1 ).LT.PIVMIN )
137 $ TMP1 = -PIVMIN
138 IF( TMP1.LE.ZERO )
139 $ NEGCNT = NEGCNT + 1
140 *
141 DO 20 I = 2, N
142 TMP1 = D( I ) - E2( I-1 ) / TMP1 - MID
143 IF( ABS( TMP1 ).LT.PIVMIN )
144 $ TMP1 = -PIVMIN
145 IF( TMP1.LE.ZERO )
146 $ NEGCNT = NEGCNT + 1
147 20 CONTINUE
148
149 IF(NEGCNT.GE.IW) THEN
150 RIGHT = MID
151 ELSE
152 LEFT = MID
153 ENDIF
154 GOTO 10
155
156 30 CONTINUE
157 *
158 * Converged or maximum number of iterations reached
159 *
160 W = HALF * (LEFT + RIGHT)
161 WERR = HALF * ABS( RIGHT - LEFT )
162
163 RETURN
164 *
165 * End of DLARRK
166 *
167 END
2 $ D, E2, PIVMIN, RELTOL, W, WERR, INFO)
3 IMPLICIT NONE
4 *
5 * -- LAPACK auxiliary routine (version 3.2) --
6 * -- LAPACK is a software package provided by Univ. of Tennessee, --
7 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
8 * November 2006
9 *
10 * .. Scalar Arguments ..
11 INTEGER INFO, IW, N
12 DOUBLE PRECISION PIVMIN, RELTOL, GL, GU, W, WERR
13 * ..
14 * .. Array Arguments ..
15 DOUBLE PRECISION D( * ), E2( * )
16 * ..
17 *
18 * Purpose
19 * =======
20 *
21 * DLARRK computes one eigenvalue of a symmetric tridiagonal
22 * matrix T to suitable accuracy. This is an auxiliary code to be
23 * called from DSTEMR.
24 *
25 * To avoid overflow, the matrix must be scaled so that its
26 * largest element is no greater than overflow**(1/2) *
27 * underflow**(1/4) in absolute value, and for greatest
28 * accuracy, it should not be much smaller than that.
29 *
30 * See W. Kahan "Accurate Eigenvalues of a Symmetric Tridiagonal
31 * Matrix", Report CS41, Computer Science Dept., Stanford
32 * University, July 21, 1966.
33 *
34 * Arguments
35 * =========
36 *
37 * N (input) INTEGER
38 * The order of the tridiagonal matrix T. N >= 0.
39 *
40 * IW (input) INTEGER
41 * The index of the eigenvalues to be returned.
42 *
43 * GL (input) DOUBLE PRECISION
44 * GU (input) DOUBLE PRECISION
45 * An upper and a lower bound on the eigenvalue.
46 *
47 * D (input) DOUBLE PRECISION array, dimension (N)
48 * The n diagonal elements of the tridiagonal matrix T.
49 *
50 * E2 (input) DOUBLE PRECISION array, dimension (N-1)
51 * The (n-1) squared off-diagonal elements of the tridiagonal matrix T.
52 *
53 * PIVMIN (input) DOUBLE PRECISION
54 * The minimum pivot allowed in the Sturm sequence for T.
55 *
56 * RELTOL (input) DOUBLE PRECISION
57 * The minimum relative width of an interval. When an interval
58 * is narrower than RELTOL times the larger (in
59 * magnitude) endpoint, then it is considered to be
60 * sufficiently small, i.e., converged. Note: this should
61 * always be at least radix*machine epsilon.
62 *
63 * W (output) DOUBLE PRECISION
64 *
65 * WERR (output) DOUBLE PRECISION
66 * The error bound on the corresponding eigenvalue approximation
67 * in W.
68 *
69 * INFO (output) INTEGER
70 * = 0: Eigenvalue converged
71 * = -1: Eigenvalue did NOT converge
72 *
73 * Internal Parameters
74 * ===================
75 *
76 * FUDGE DOUBLE PRECISION, default = 2
77 * A "fudge factor" to widen the Gershgorin intervals.
78 *
79 * =====================================================================
80 *
81 * .. Parameters ..
82 DOUBLE PRECISION FUDGE, HALF, TWO, ZERO
83 PARAMETER ( HALF = 0.5D0, TWO = 2.0D0,
84 $ FUDGE = TWO, ZERO = 0.0D0 )
85 * ..
86 * .. Local Scalars ..
87 INTEGER I, IT, ITMAX, NEGCNT
88 DOUBLE PRECISION ATOLI, EPS, LEFT, MID, RIGHT, RTOLI, TMP1,
89 $ TMP2, TNORM
90 * ..
91 * .. External Functions ..
92 DOUBLE PRECISION DLAMCH
93 EXTERNAL DLAMCH
94 * ..
95 * .. Intrinsic Functions ..
96 INTRINSIC ABS, INT, LOG, MAX
97 * ..
98 * .. Executable Statements ..
99 *
100 * Get machine constants
101 EPS = DLAMCH( 'P' )
102
103 TNORM = MAX( ABS( GL ), ABS( GU ) )
104 RTOLI = RELTOL
105 ATOLI = FUDGE*TWO*PIVMIN
106
107 ITMAX = INT( ( LOG( TNORM+PIVMIN )-LOG( PIVMIN ) ) /
108 $ LOG( TWO ) ) + 2
109
110 INFO = -1
111
112 LEFT = GL - FUDGE*TNORM*EPS*N - FUDGE*TWO*PIVMIN
113 RIGHT = GU + FUDGE*TNORM*EPS*N + FUDGE*TWO*PIVMIN
114 IT = 0
115
116 10 CONTINUE
117 *
118 * Check if interval converged or maximum number of iterations reached
119 *
120 TMP1 = ABS( RIGHT - LEFT )
121 TMP2 = MAX( ABS(RIGHT), ABS(LEFT) )
122 IF( TMP1.LT.MAX( ATOLI, PIVMIN, RTOLI*TMP2 ) ) THEN
123 INFO = 0
124 GOTO 30
125 ENDIF
126 IF(IT.GT.ITMAX)
127 $ GOTO 30
128
129 *
130 * Count number of negative pivots for mid-point
131 *
132 IT = IT + 1
133 MID = HALF * (LEFT + RIGHT)
134 NEGCNT = 0
135 TMP1 = D( 1 ) - MID
136 IF( ABS( TMP1 ).LT.PIVMIN )
137 $ TMP1 = -PIVMIN
138 IF( TMP1.LE.ZERO )
139 $ NEGCNT = NEGCNT + 1
140 *
141 DO 20 I = 2, N
142 TMP1 = D( I ) - E2( I-1 ) / TMP1 - MID
143 IF( ABS( TMP1 ).LT.PIVMIN )
144 $ TMP1 = -PIVMIN
145 IF( TMP1.LE.ZERO )
146 $ NEGCNT = NEGCNT + 1
147 20 CONTINUE
148
149 IF(NEGCNT.GE.IW) THEN
150 RIGHT = MID
151 ELSE
152 LEFT = MID
153 ENDIF
154 GOTO 10
155
156 30 CONTINUE
157 *
158 * Converged or maximum number of iterations reached
159 *
160 W = HALF * (LEFT + RIGHT)
161 WERR = HALF * ABS( RIGHT - LEFT )
162
163 RETURN
164 *
165 * End of DLARRK
166 *
167 END