1 SUBROUTINE SSTECH( N, A, B, EIG, TOL, WORK, INFO )
2 *
3 * -- LAPACK test routine (version 3.1) --
4 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
5 * November 2006
6 *
7 * .. Scalar Arguments ..
8 INTEGER INFO, N
9 REAL TOL
10 * ..
11 * .. Array Arguments ..
12 REAL A( * ), B( * ), EIG( * ), WORK( * )
13 * ..
14 *
15 * Purpose
16 * =======
17 *
18 * Let T be the tridiagonal matrix with diagonal entries A(1) ,...,
19 * A(N) and offdiagonal entries B(1) ,..., B(N-1)). SSTECH checks to
20 * see if EIG(1) ,..., EIG(N) are indeed accurate eigenvalues of T.
21 * It does this by expanding each EIG(I) into an interval
22 * [SVD(I) - EPS, SVD(I) + EPS], merging overlapping intervals if
23 * any, and using Sturm sequences to count and verify whether each
24 * resulting interval has the correct number of eigenvalues (using
25 * SSTECT). Here EPS = TOL*MACHEPS*MAXEIG, where MACHEPS is the
26 * machine precision and MAXEIG is the absolute value of the largest
27 * eigenvalue. If each interval contains the correct number of
28 * eigenvalues, INFO = 0 is returned, otherwise INFO is the index of
29 * the first eigenvalue in the first bad interval.
30 *
31 * Arguments
32 * =========
33 *
34 * N (input) INTEGER
35 * The dimension of the tridiagonal matrix T.
36 *
37 * A (input) REAL array, dimension (N)
38 * The diagonal entries of the tridiagonal matrix T.
39 *
40 * B (input) REAL array, dimension (N-1)
41 * The offdiagonal entries of the tridiagonal matrix T.
42 *
43 * EIG (input) REAL array, dimension (N)
44 * The purported eigenvalues to be checked.
45 *
46 * TOL (input) REAL
47 * Error tolerance for checking, a multiple of the
48 * machine precision.
49 *
50 * WORK (workspace) REAL array, dimension (N)
51 *
52 * INFO (output) INTEGER
53 * 0 if the eigenvalues are all correct (to within
54 * 1 +- TOL*MACHEPS*MAXEIG)
55 * >0 if the interval containing the INFO-th eigenvalue
56 * contains the incorrect number of eigenvalues.
57 *
58 * =====================================================================
59 *
60 * .. Parameters ..
61 REAL ZERO
62 PARAMETER ( ZERO = 0.0E0 )
63 * ..
64 * .. Local Scalars ..
65 INTEGER BPNT, COUNT, I, ISUB, J, NUML, NUMU, TPNT
66 REAL EMIN, EPS, LOWER, MX, TUPPR, UNFLEP, UPPER
67 * ..
68 * .. External Functions ..
69 REAL SLAMCH
70 EXTERNAL SLAMCH
71 * ..
72 * .. External Subroutines ..
73 EXTERNAL SSTECT
74 * ..
75 * .. Intrinsic Functions ..
76 INTRINSIC ABS, MAX
77 * ..
78 * .. Executable Statements ..
79 *
80 * Check input parameters
81 *
82 INFO = 0
83 IF( N.EQ.0 )
84 $ RETURN
85 IF( N.LT.0 ) THEN
86 INFO = -1
87 RETURN
88 END IF
89 IF( TOL.LT.ZERO ) THEN
90 INFO = -5
91 RETURN
92 END IF
93 *
94 * Get machine constants
95 *
96 EPS = SLAMCH( 'Epsilon' )*SLAMCH( 'Base' )
97 UNFLEP = SLAMCH( 'Safe minimum' ) / EPS
98 EPS = TOL*EPS
99 *
100 * Compute maximum absolute eigenvalue, error tolerance
101 *
102 MX = ABS( EIG( 1 ) )
103 DO 10 I = 2, N
104 MX = MAX( MX, ABS( EIG( I ) ) )
105 10 CONTINUE
106 EPS = MAX( EPS*MX, UNFLEP )
107 *
108 * Sort eigenvalues from EIG into WORK
109 *
110 DO 20 I = 1, N
111 WORK( I ) = EIG( I )
112 20 CONTINUE
113 DO 40 I = 1, N - 1
114 ISUB = 1
115 EMIN = WORK( 1 )
116 DO 30 J = 2, N + 1 - I
117 IF( WORK( J ).LT.EMIN ) THEN
118 ISUB = J
119 EMIN = WORK( J )
120 END IF
121 30 CONTINUE
122 IF( ISUB.NE.N+1-I ) THEN
123 WORK( ISUB ) = WORK( N+1-I )
124 WORK( N+1-I ) = EMIN
125 END IF
126 40 CONTINUE
127 *
128 * TPNT points to singular value at right endpoint of interval
129 * BPNT points to singular value at left endpoint of interval
130 *
131 TPNT = 1
132 BPNT = 1
133 *
134 * Begin loop over all intervals
135 *
136 50 CONTINUE
137 UPPER = WORK( TPNT ) + EPS
138 LOWER = WORK( BPNT ) - EPS
139 *
140 * Begin loop merging overlapping intervals
141 *
142 60 CONTINUE
143 IF( BPNT.EQ.N )
144 $ GO TO 70
145 TUPPR = WORK( BPNT+1 ) + EPS
146 IF( TUPPR.LT.LOWER )
147 $ GO TO 70
148 *
149 * Merge
150 *
151 BPNT = BPNT + 1
152 LOWER = WORK( BPNT ) - EPS
153 GO TO 60
154 70 CONTINUE
155 *
156 * Count singular values in interval [ LOWER, UPPER ]
157 *
158 CALL SSTECT( N, A, B, LOWER, NUML )
159 CALL SSTECT( N, A, B, UPPER, NUMU )
160 COUNT = NUMU - NUML
161 IF( COUNT.NE.BPNT-TPNT+1 ) THEN
162 *
163 * Wrong number of singular values in interval
164 *
165 INFO = TPNT
166 GO TO 80
167 END IF
168 TPNT = BPNT + 1
169 BPNT = TPNT
170 IF( TPNT.LE.N )
171 $ GO TO 50
172 80 CONTINUE
173 RETURN
174 *
175 * End of SSTECH
176 *
177 END
2 *
3 * -- LAPACK test routine (version 3.1) --
4 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
5 * November 2006
6 *
7 * .. Scalar Arguments ..
8 INTEGER INFO, N
9 REAL TOL
10 * ..
11 * .. Array Arguments ..
12 REAL A( * ), B( * ), EIG( * ), WORK( * )
13 * ..
14 *
15 * Purpose
16 * =======
17 *
18 * Let T be the tridiagonal matrix with diagonal entries A(1) ,...,
19 * A(N) and offdiagonal entries B(1) ,..., B(N-1)). SSTECH checks to
20 * see if EIG(1) ,..., EIG(N) are indeed accurate eigenvalues of T.
21 * It does this by expanding each EIG(I) into an interval
22 * [SVD(I) - EPS, SVD(I) + EPS], merging overlapping intervals if
23 * any, and using Sturm sequences to count and verify whether each
24 * resulting interval has the correct number of eigenvalues (using
25 * SSTECT). Here EPS = TOL*MACHEPS*MAXEIG, where MACHEPS is the
26 * machine precision and MAXEIG is the absolute value of the largest
27 * eigenvalue. If each interval contains the correct number of
28 * eigenvalues, INFO = 0 is returned, otherwise INFO is the index of
29 * the first eigenvalue in the first bad interval.
30 *
31 * Arguments
32 * =========
33 *
34 * N (input) INTEGER
35 * The dimension of the tridiagonal matrix T.
36 *
37 * A (input) REAL array, dimension (N)
38 * The diagonal entries of the tridiagonal matrix T.
39 *
40 * B (input) REAL array, dimension (N-1)
41 * The offdiagonal entries of the tridiagonal matrix T.
42 *
43 * EIG (input) REAL array, dimension (N)
44 * The purported eigenvalues to be checked.
45 *
46 * TOL (input) REAL
47 * Error tolerance for checking, a multiple of the
48 * machine precision.
49 *
50 * WORK (workspace) REAL array, dimension (N)
51 *
52 * INFO (output) INTEGER
53 * 0 if the eigenvalues are all correct (to within
54 * 1 +- TOL*MACHEPS*MAXEIG)
55 * >0 if the interval containing the INFO-th eigenvalue
56 * contains the incorrect number of eigenvalues.
57 *
58 * =====================================================================
59 *
60 * .. Parameters ..
61 REAL ZERO
62 PARAMETER ( ZERO = 0.0E0 )
63 * ..
64 * .. Local Scalars ..
65 INTEGER BPNT, COUNT, I, ISUB, J, NUML, NUMU, TPNT
66 REAL EMIN, EPS, LOWER, MX, TUPPR, UNFLEP, UPPER
67 * ..
68 * .. External Functions ..
69 REAL SLAMCH
70 EXTERNAL SLAMCH
71 * ..
72 * .. External Subroutines ..
73 EXTERNAL SSTECT
74 * ..
75 * .. Intrinsic Functions ..
76 INTRINSIC ABS, MAX
77 * ..
78 * .. Executable Statements ..
79 *
80 * Check input parameters
81 *
82 INFO = 0
83 IF( N.EQ.0 )
84 $ RETURN
85 IF( N.LT.0 ) THEN
86 INFO = -1
87 RETURN
88 END IF
89 IF( TOL.LT.ZERO ) THEN
90 INFO = -5
91 RETURN
92 END IF
93 *
94 * Get machine constants
95 *
96 EPS = SLAMCH( 'Epsilon' )*SLAMCH( 'Base' )
97 UNFLEP = SLAMCH( 'Safe minimum' ) / EPS
98 EPS = TOL*EPS
99 *
100 * Compute maximum absolute eigenvalue, error tolerance
101 *
102 MX = ABS( EIG( 1 ) )
103 DO 10 I = 2, N
104 MX = MAX( MX, ABS( EIG( I ) ) )
105 10 CONTINUE
106 EPS = MAX( EPS*MX, UNFLEP )
107 *
108 * Sort eigenvalues from EIG into WORK
109 *
110 DO 20 I = 1, N
111 WORK( I ) = EIG( I )
112 20 CONTINUE
113 DO 40 I = 1, N - 1
114 ISUB = 1
115 EMIN = WORK( 1 )
116 DO 30 J = 2, N + 1 - I
117 IF( WORK( J ).LT.EMIN ) THEN
118 ISUB = J
119 EMIN = WORK( J )
120 END IF
121 30 CONTINUE
122 IF( ISUB.NE.N+1-I ) THEN
123 WORK( ISUB ) = WORK( N+1-I )
124 WORK( N+1-I ) = EMIN
125 END IF
126 40 CONTINUE
127 *
128 * TPNT points to singular value at right endpoint of interval
129 * BPNT points to singular value at left endpoint of interval
130 *
131 TPNT = 1
132 BPNT = 1
133 *
134 * Begin loop over all intervals
135 *
136 50 CONTINUE
137 UPPER = WORK( TPNT ) + EPS
138 LOWER = WORK( BPNT ) - EPS
139 *
140 * Begin loop merging overlapping intervals
141 *
142 60 CONTINUE
143 IF( BPNT.EQ.N )
144 $ GO TO 70
145 TUPPR = WORK( BPNT+1 ) + EPS
146 IF( TUPPR.LT.LOWER )
147 $ GO TO 70
148 *
149 * Merge
150 *
151 BPNT = BPNT + 1
152 LOWER = WORK( BPNT ) - EPS
153 GO TO 60
154 70 CONTINUE
155 *
156 * Count singular values in interval [ LOWER, UPPER ]
157 *
158 CALL SSTECT( N, A, B, LOWER, NUML )
159 CALL SSTECT( N, A, B, UPPER, NUMU )
160 COUNT = NUMU - NUML
161 IF( COUNT.NE.BPNT-TPNT+1 ) THEN
162 *
163 * Wrong number of singular values in interval
164 *
165 INFO = TPNT
166 GO TO 80
167 END IF
168 TPNT = BPNT + 1
169 BPNT = TPNT
170 IF( TPNT.LE.N )
171 $ GO TO 50
172 80 CONTINUE
173 RETURN
174 *
175 * End of SSTECH
176 *
177 END