1 SUBROUTINE SSVDCH( N, S, E, SVD, TOL, 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 E( * ), S( * ), SVD( * )
13 * ..
14 *
15 * Purpose
16 * =======
17 *
18 * SSVDCH checks to see if SVD(1) ,..., SVD(N) are accurate singular
19 * values of the bidiagonal matrix B with diagonal entries
20 * S(1) ,..., S(N) and superdiagonal entries E(1) ,..., E(N-1)).
21 * It does this by expanding each SVD(I) into an interval
22 * [SVD(I) * (1-EPS) , SVD(I) * (1+EPS)], merging overlapping intervals
23 * if any, and using Sturm sequences to count and verify whether each
24 * resulting interval has the correct number of singular values (using
25 * SSVDCT). Here EPS=TOL*MAX(N/10,1)*MACHEP, where MACHEP is the
26 * machine precision. The routine assumes the singular values are sorted
27 * with SVD(1) the largest and SVD(N) smallest. If each interval
28 * contains the correct number of singular values, INFO = 0 is returned,
29 * otherwise INFO is the index of the first singular value in the first
30 * bad interval.
31 *
32 * Arguments
33 * ==========
34 *
35 * N (input) INTEGER
36 * The dimension of the bidiagonal matrix B.
37 *
38 * S (input) REAL array, dimension (N)
39 * The diagonal entries of the bidiagonal matrix B.
40 *
41 * E (input) REAL array, dimension (N-1)
42 * The superdiagonal entries of the bidiagonal matrix B.
43 *
44 * SVD (input) REAL array, dimension (N)
45 * The computed singular values to be checked.
46 *
47 * TOL (input) REAL
48 * Error tolerance for checking, a multiplier of the
49 * machine precision.
50 *
51 * INFO (output) INTEGER
52 * =0 if the singular values are all correct (to within
53 * 1 +- TOL*MACHEPS)
54 * >0 if the interval containing the INFO-th singular value
55 * contains the incorrect number of singular values.
56 *
57 * =====================================================================
58 *
59 * .. Parameters ..
60 REAL ONE
61 PARAMETER ( ONE = 1.0E0 )
62 REAL ZERO
63 PARAMETER ( ZERO = 0.0E0 )
64 * ..
65 * .. Local Scalars ..
66 INTEGER BPNT, COUNT, NUML, NUMU, TPNT
67 REAL EPS, LOWER, OVFL, TUPPR, UNFL, UNFLEP, UPPER
68 * ..
69 * .. External Functions ..
70 REAL SLAMCH
71 EXTERNAL SLAMCH
72 * ..
73 * .. External Subroutines ..
74 EXTERNAL SSVDCT
75 * ..
76 * .. Intrinsic Functions ..
77 INTRINSIC MAX, SQRT
78 * ..
79 * .. Executable Statements ..
80 *
81 * Get machine constants
82 *
83 INFO = 0
84 IF( N.LE.0 )
85 $ RETURN
86 UNFL = SLAMCH( 'Safe minimum' )
87 OVFL = SLAMCH( 'Overflow' )
88 EPS = SLAMCH( 'Epsilon' )*SLAMCH( 'Base' )
89 *
90 * UNFLEP is chosen so that when an eigenvalue is multiplied by the
91 * scale factor sqrt(OVFL)*sqrt(sqrt(UNFL))/MX in SSVDCT, it exceeds
92 * sqrt(UNFL), which is the lower limit for SSVDCT.
93 *
94 UNFLEP = ( SQRT( SQRT( UNFL ) ) / SQRT( OVFL ) )*SVD( 1 ) +
95 $ UNFL / EPS
96 *
97 * The value of EPS works best when TOL .GE. 10.
98 *
99 EPS = TOL*MAX( N / 10, 1 )*EPS
100 *
101 * TPNT points to singular value at right endpoint of interval
102 * BPNT points to singular value at left endpoint of interval
103 *
104 TPNT = 1
105 BPNT = 1
106 *
107 * Begin loop over all intervals
108 *
109 10 CONTINUE
110 UPPER = ( ONE+EPS )*SVD( TPNT ) + UNFLEP
111 LOWER = ( ONE-EPS )*SVD( BPNT ) - UNFLEP
112 IF( LOWER.LE.UNFLEP )
113 $ LOWER = -UPPER
114 *
115 * Begin loop merging overlapping intervals
116 *
117 20 CONTINUE
118 IF( BPNT.EQ.N )
119 $ GO TO 30
120 TUPPR = ( ONE+EPS )*SVD( BPNT+1 ) + UNFLEP
121 IF( TUPPR.LT.LOWER )
122 $ GO TO 30
123 *
124 * Merge
125 *
126 BPNT = BPNT + 1
127 LOWER = ( ONE-EPS )*SVD( BPNT ) - UNFLEP
128 IF( LOWER.LE.UNFLEP )
129 $ LOWER = -UPPER
130 GO TO 20
131 30 CONTINUE
132 *
133 * Count singular values in interval [ LOWER, UPPER ]
134 *
135 CALL SSVDCT( N, S, E, LOWER, NUML )
136 CALL SSVDCT( N, S, E, UPPER, NUMU )
137 COUNT = NUMU - NUML
138 IF( LOWER.LT.ZERO )
139 $ COUNT = COUNT / 2
140 IF( COUNT.NE.BPNT-TPNT+1 ) THEN
141 *
142 * Wrong number of singular values in interval
143 *
144 INFO = TPNT
145 GO TO 40
146 END IF
147 TPNT = BPNT + 1
148 BPNT = TPNT
149 IF( TPNT.LE.N )
150 $ GO TO 10
151 40 CONTINUE
152 RETURN
153 *
154 * End of SSVDCH
155 *
156 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 E( * ), S( * ), SVD( * )
13 * ..
14 *
15 * Purpose
16 * =======
17 *
18 * SSVDCH checks to see if SVD(1) ,..., SVD(N) are accurate singular
19 * values of the bidiagonal matrix B with diagonal entries
20 * S(1) ,..., S(N) and superdiagonal entries E(1) ,..., E(N-1)).
21 * It does this by expanding each SVD(I) into an interval
22 * [SVD(I) * (1-EPS) , SVD(I) * (1+EPS)], merging overlapping intervals
23 * if any, and using Sturm sequences to count and verify whether each
24 * resulting interval has the correct number of singular values (using
25 * SSVDCT). Here EPS=TOL*MAX(N/10,1)*MACHEP, where MACHEP is the
26 * machine precision. The routine assumes the singular values are sorted
27 * with SVD(1) the largest and SVD(N) smallest. If each interval
28 * contains the correct number of singular values, INFO = 0 is returned,
29 * otherwise INFO is the index of the first singular value in the first
30 * bad interval.
31 *
32 * Arguments
33 * ==========
34 *
35 * N (input) INTEGER
36 * The dimension of the bidiagonal matrix B.
37 *
38 * S (input) REAL array, dimension (N)
39 * The diagonal entries of the bidiagonal matrix B.
40 *
41 * E (input) REAL array, dimension (N-1)
42 * The superdiagonal entries of the bidiagonal matrix B.
43 *
44 * SVD (input) REAL array, dimension (N)
45 * The computed singular values to be checked.
46 *
47 * TOL (input) REAL
48 * Error tolerance for checking, a multiplier of the
49 * machine precision.
50 *
51 * INFO (output) INTEGER
52 * =0 if the singular values are all correct (to within
53 * 1 +- TOL*MACHEPS)
54 * >0 if the interval containing the INFO-th singular value
55 * contains the incorrect number of singular values.
56 *
57 * =====================================================================
58 *
59 * .. Parameters ..
60 REAL ONE
61 PARAMETER ( ONE = 1.0E0 )
62 REAL ZERO
63 PARAMETER ( ZERO = 0.0E0 )
64 * ..
65 * .. Local Scalars ..
66 INTEGER BPNT, COUNT, NUML, NUMU, TPNT
67 REAL EPS, LOWER, OVFL, TUPPR, UNFL, UNFLEP, UPPER
68 * ..
69 * .. External Functions ..
70 REAL SLAMCH
71 EXTERNAL SLAMCH
72 * ..
73 * .. External Subroutines ..
74 EXTERNAL SSVDCT
75 * ..
76 * .. Intrinsic Functions ..
77 INTRINSIC MAX, SQRT
78 * ..
79 * .. Executable Statements ..
80 *
81 * Get machine constants
82 *
83 INFO = 0
84 IF( N.LE.0 )
85 $ RETURN
86 UNFL = SLAMCH( 'Safe minimum' )
87 OVFL = SLAMCH( 'Overflow' )
88 EPS = SLAMCH( 'Epsilon' )*SLAMCH( 'Base' )
89 *
90 * UNFLEP is chosen so that when an eigenvalue is multiplied by the
91 * scale factor sqrt(OVFL)*sqrt(sqrt(UNFL))/MX in SSVDCT, it exceeds
92 * sqrt(UNFL), which is the lower limit for SSVDCT.
93 *
94 UNFLEP = ( SQRT( SQRT( UNFL ) ) / SQRT( OVFL ) )*SVD( 1 ) +
95 $ UNFL / EPS
96 *
97 * The value of EPS works best when TOL .GE. 10.
98 *
99 EPS = TOL*MAX( N / 10, 1 )*EPS
100 *
101 * TPNT points to singular value at right endpoint of interval
102 * BPNT points to singular value at left endpoint of interval
103 *
104 TPNT = 1
105 BPNT = 1
106 *
107 * Begin loop over all intervals
108 *
109 10 CONTINUE
110 UPPER = ( ONE+EPS )*SVD( TPNT ) + UNFLEP
111 LOWER = ( ONE-EPS )*SVD( BPNT ) - UNFLEP
112 IF( LOWER.LE.UNFLEP )
113 $ LOWER = -UPPER
114 *
115 * Begin loop merging overlapping intervals
116 *
117 20 CONTINUE
118 IF( BPNT.EQ.N )
119 $ GO TO 30
120 TUPPR = ( ONE+EPS )*SVD( BPNT+1 ) + UNFLEP
121 IF( TUPPR.LT.LOWER )
122 $ GO TO 30
123 *
124 * Merge
125 *
126 BPNT = BPNT + 1
127 LOWER = ( ONE-EPS )*SVD( BPNT ) - UNFLEP
128 IF( LOWER.LE.UNFLEP )
129 $ LOWER = -UPPER
130 GO TO 20
131 30 CONTINUE
132 *
133 * Count singular values in interval [ LOWER, UPPER ]
134 *
135 CALL SSVDCT( N, S, E, LOWER, NUML )
136 CALL SSVDCT( N, S, E, UPPER, NUMU )
137 COUNT = NUMU - NUML
138 IF( LOWER.LT.ZERO )
139 $ COUNT = COUNT / 2
140 IF( COUNT.NE.BPNT-TPNT+1 ) THEN
141 *
142 * Wrong number of singular values in interval
143 *
144 INFO = TPNT
145 GO TO 40
146 END IF
147 TPNT = BPNT + 1
148 BPNT = TPNT
149 IF( TPNT.LE.N )
150 $ GO TO 10
151 40 CONTINUE
152 RETURN
153 *
154 * End of SSVDCH
155 *
156 END