1 SUBROUTINE DSVDCH( 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 DOUBLE PRECISION TOL
10 * ..
11 * .. Array Arguments ..
12 DOUBLE PRECISION E( * ), S( * ), SVD( * )
13 * ..
14 *
15 * Purpose
16 * =======
17 *
18 * DSVDCH 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 * DSVDCT). Here EPS=TOL*MAX(N/10,1)*MAZHEP, 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) DOUBLE PRECISION array, dimension (N)
39 * The diagonal entries of the bidiagonal matrix B.
40 *
41 * E (input) DOUBLE PRECISION array, dimension (N-1)
42 * The superdiagonal entries of the bidiagonal matrix B.
43 *
44 * SVD (input) DOUBLE PRECISION array, dimension (N)
45 * The computed singular values to be checked.
46 *
47 * TOL (input) DOUBLE PRECISION
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*MAZHEPS)
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 DOUBLE PRECISION ONE
61 PARAMETER ( ONE = 1.0D0 )
62 DOUBLE PRECISION ZERO
63 PARAMETER ( ZERO = 0.0D0 )
64 * ..
65 * .. Local Scalars ..
66 INTEGER BPNT, COUNT, NUML, NUMU, TPNT
67 DOUBLE PRECISION EPS, LOWER, OVFL, TUPPR, UNFL, UNFLEP, UPPER
68 * ..
69 * .. External Functions ..
70 DOUBLE PRECISION DLAMCH
71 EXTERNAL DLAMCH
72 * ..
73 * .. External Subroutines ..
74 EXTERNAL DSVDCT
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 = DLAMCH( 'Safe minimum' )
87 OVFL = DLAMCH( 'Overflow' )
88 EPS = DLAMCH( 'Epsilon' )*DLAMCH( '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 DSVDCT, it exceeds
92 * sqrt(UNFL), which is the lower limit for DSVDCT.
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 DSVDCT( N, S, E, LOWER, NUML )
136 CALL DSVDCT( 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 DSVDCH
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 DOUBLE PRECISION TOL
10 * ..
11 * .. Array Arguments ..
12 DOUBLE PRECISION E( * ), S( * ), SVD( * )
13 * ..
14 *
15 * Purpose
16 * =======
17 *
18 * DSVDCH 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 * DSVDCT). Here EPS=TOL*MAX(N/10,1)*MAZHEP, 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) DOUBLE PRECISION array, dimension (N)
39 * The diagonal entries of the bidiagonal matrix B.
40 *
41 * E (input) DOUBLE PRECISION array, dimension (N-1)
42 * The superdiagonal entries of the bidiagonal matrix B.
43 *
44 * SVD (input) DOUBLE PRECISION array, dimension (N)
45 * The computed singular values to be checked.
46 *
47 * TOL (input) DOUBLE PRECISION
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*MAZHEPS)
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 DOUBLE PRECISION ONE
61 PARAMETER ( ONE = 1.0D0 )
62 DOUBLE PRECISION ZERO
63 PARAMETER ( ZERO = 0.0D0 )
64 * ..
65 * .. Local Scalars ..
66 INTEGER BPNT, COUNT, NUML, NUMU, TPNT
67 DOUBLE PRECISION EPS, LOWER, OVFL, TUPPR, UNFL, UNFLEP, UPPER
68 * ..
69 * .. External Functions ..
70 DOUBLE PRECISION DLAMCH
71 EXTERNAL DLAMCH
72 * ..
73 * .. External Subroutines ..
74 EXTERNAL DSVDCT
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 = DLAMCH( 'Safe minimum' )
87 OVFL = DLAMCH( 'Overflow' )
88 EPS = DLAMCH( 'Epsilon' )*DLAMCH( '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 DSVDCT, it exceeds
92 * sqrt(UNFL), which is the lower limit for DSVDCT.
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 DSVDCT( N, S, E, LOWER, NUML )
136 CALL DSVDCT( 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 DSVDCH
155 *
156 END