1 SUBROUTINE SCSDTS( M, P, Q, X, XF, LDX, U1, LDU1, U2, LDU2, V1T,
2 $ LDV1T, V2T, LDV2T, THETA, IWORK, WORK, LWORK,
3 $ RWORK, RESULT )
4 IMPLICIT NONE
5 *
6 * Originally xGSVTS
7 * -- LAPACK test routine (version 3.1) --
8 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
9 * November 2006
10 *
11 * Adapted to SCSDTS
12 * July 2010
13 *
14 * .. Scalar Arguments ..
15 INTEGER LDX, LDU1, LDU2, LDV1T, LDV2T, LWORK, M, P, Q
16 * ..
17 * .. Array Arguments ..
18 INTEGER IWORK( * )
19 REAL RESULT( 9 ), RWORK( * ), THETA( * )
20 REAL U1( LDU1, * ), U2( LDU2, * ), V1T( LDV1T, * ),
21 $ V2T( LDV2T, * ), WORK( LWORK ), X( LDX, * ),
22 $ XF( LDX, * )
23 * ..
24 *
25 * Purpose
26 * =======
27 *
28 * SCSDTS tests SORCSD, which, given an M-by-M partitioned orthogonal
29 * matrix X,
30 * Q M-Q
31 * X = [ X11 X12 ] P ,
32 * [ X21 X22 ] M-P
33 *
34 * computes the CSD
35 *
36 * [ U1 ]**T * [ X11 X12 ] * [ V1 ]
37 * [ U2 ] [ X21 X22 ] [ V2 ]
38 *
39 * [ I 0 0 | 0 0 0 ]
40 * [ 0 C 0 | 0 -S 0 ]
41 * [ 0 0 0 | 0 0 -I ]
42 * = [---------------------] = [ D11 D12 ] .
43 * [ 0 0 0 | I 0 0 ] [ D21 D22 ]
44 * [ 0 S 0 | 0 C 0 ]
45 * [ 0 0 I | 0 0 0 ]
46 *
47 * Arguments
48 * =========
49 *
50 * M (input) INTEGER
51 * The number of rows of the matrix X. M >= 0.
52 *
53 * P (input) INTEGER
54 * The number of rows of the matrix X11. P >= 0.
55 *
56 * Q (input) INTEGER
57 * The number of columns of the matrix X11. Q >= 0.
58 *
59 * X (input) REAL array, dimension (LDX,M)
60 * The M-by-M matrix X.
61 *
62 * XF (output) REAL array, dimension (LDX,M)
63 * Details of the CSD of X, as returned by SORCSD;
64 * see SORCSD for further details.
65 *
66 * LDX (input) INTEGER
67 * The leading dimension of the arrays X and XF.
68 * LDX >= max( 1,M ).
69 *
70 * U1 (output) REAL array, dimension(LDU1,P)
71 * The P-by-P orthogonal matrix U1.
72 *
73 * LDU1 (input) INTEGER
74 * The leading dimension of the array U1. LDU >= max(1,P).
75 *
76 * U2 (output) REAL array, dimension(LDU2,M-P)
77 * The (M-P)-by-(M-P) orthogonal matrix U2.
78 *
79 * LDU2 (input) INTEGER
80 * The leading dimension of the array U2. LDU >= max(1,M-P).
81 *
82 * V1T (output) REAL array, dimension(LDV1T,Q)
83 * The Q-by-Q orthogonal matrix V1T.
84 *
85 * LDV1T (input) INTEGER
86 * The leading dimension of the array V1T. LDV1T >=
87 * max(1,Q).
88 *
89 * V2T (output) REAL array, dimension(LDV2T,M-Q)
90 * The (M-Q)-by-(M-Q) orthogonal matrix V2T.
91 *
92 * LDV2T (input) INTEGER
93 * The leading dimension of the array V2T. LDV2T >=
94 * max(1,M-Q).
95 *
96 * THETA (output) REAL array, dimension MIN(P,M-P,Q,M-Q)
97 * The CS values of X; the essentially diagonal matrices C and
98 * S are constructed from THETA; see subroutine SORCSD for
99 * details.
100 *
101 * IWORK (workspace) INTEGER array, dimension (M)
102 *
103 * WORK (workspace) REAL array, dimension (LWORK)
104 *
105 * LWORK (input) INTEGER
106 * The dimension of the array WORK
107 *
108 * RWORK (workspace) REAL array
109 *
110 * RESULT (output) REAL array, dimension (9)
111 * The test ratios:
112 * RESULT(1) = norm( U1'*X11*V1 - D11 ) / ( MAX(1,P,Q)*EPS2 )
113 * RESULT(2) = norm( U1'*X12*V2 - D12 ) / ( MAX(1,P,M-Q)*EPS2 )
114 * RESULT(3) = norm( U2'*X21*V1 - D21 ) / ( MAX(1,M-P,Q)*EPS2 )
115 * RESULT(4) = norm( U2'*X22*V2 - D22 ) / ( MAX(1,M-P,M-Q)*EPS2 )
116 * RESULT(5) = norm( I - U1'*U1 ) / ( MAX(1,P)*ULP )
117 * RESULT(6) = norm( I - U2'*U2 ) / ( MAX(1,M-P)*ULP )
118 * RESULT(7) = norm( I - V1T'*V1T ) / ( MAX(1,Q)*ULP )
119 * RESULT(8) = norm( I - V2T'*V2T ) / ( MAX(1,M-Q)*ULP )
120 * RESULT(9) = 0 if THETA is in increasing order and
121 * all angles are in [0,pi/2];
122 * = ULPINV otherwise.
123 * ( EPS2 = MAX( norm( I - X'*X ) / M, ULP ). )
124 *
125 * =====================================================================
126 *
127 * .. Parameters ..
128 REAL PIOVER2, REALONE, REALZERO
129 PARAMETER ( PIOVER2 = 1.57079632679489662E0,
130 $ REALONE = 1.0E0, REALZERO = 0.0E0 )
131 REAL ZERO, ONE
132 PARAMETER ( ZERO = 0.0E0, ONE = 1.0E0 )
133 * ..
134 * .. Local Scalars ..
135 INTEGER I, INFO, R
136 REAL EPS2, RESID, ULP, ULPINV
137 * ..
138 * .. External Functions ..
139 REAL SLAMCH, SLANGE, SLANSY
140 EXTERNAL SLAMCH, SLANGE, SLANSY
141 * ..
142 * .. External Subroutines ..
143 EXTERNAL SGEMM, SLACPY, SLASET, SORCSD, SSYRK
144 * ..
145 * .. Intrinsic Functions ..
146 INTRINSIC REAL, MAX, MIN
147 * ..
148 * .. Executable Statements ..
149 *
150 ULP = SLAMCH( 'Precision' )
151 ULPINV = REALONE / ULP
152 CALL SLASET( 'Full', M, M, ZERO, ONE, WORK, LDX )
153 CALL SSYRK( 'Upper', 'Conjugate transpose', M, M, -ONE, X, LDX,
154 $ ONE, WORK, LDX )
155 EPS2 = MAX( ULP,
156 $ SLANGE( '1', M, M, WORK, LDX, RWORK ) / REAL( M ) )
157 R = MIN( P, M-P, Q, M-Q )
158 *
159 * Copy the matrix X to the array XF.
160 *
161 CALL SLACPY( 'Full', M, M, X, LDX, XF, LDX )
162 *
163 * Compute the CSD
164 *
165 CALL SORCSD( 'Y', 'Y', 'Y', 'Y', 'N', 'D', M, P, Q, XF(1,1), LDX,
166 $ XF(1,Q+1), LDX, XF(P+1,1), LDX, XF(P+1,Q+1), LDX,
167 $ THETA, U1, LDU1, U2, LDU2, V1T, LDV1T, V2T, LDV2T,
168 $ WORK, LWORK, IWORK, INFO )
169 *
170 * Compute X := diag(U1,U2)'*X*diag(V1,V2) - [D11 D12; D21 D22]
171 *
172 CALL SGEMM( 'No transpose', 'Conjugate transpose', P, Q, Q, ONE,
173 $ X, LDX, V1T, LDV1T, ZERO, WORK, LDX )
174 *
175 CALL SGEMM( 'Conjugate transpose', 'No transpose', P, Q, P, ONE,
176 $ U1, LDU1, WORK, LDX, ZERO, X, LDX )
177 *
178 DO I = 1, MIN(P,Q)-R
179 X(I,I) = X(I,I) - ONE
180 END DO
181 DO I = 1, R
182 X(MIN(P,Q)-R+I,MIN(P,Q)-R+I) =
183 $ X(MIN(P,Q)-R+I,MIN(P,Q)-R+I) - COS(THETA(I))
184 END DO
185 *
186 CALL SGEMM( 'No transpose', 'Conjugate transpose', P, M-Q, M-Q,
187 $ ONE, X(1,Q+1), LDX, V2T, LDV2T, ZERO, WORK, LDX )
188 *
189 CALL SGEMM( 'Conjugate transpose', 'No transpose', P, M-Q, P,
190 $ ONE, U1, LDU1, WORK, LDX, ZERO, X(1,Q+1), LDX )
191 *
192 DO I = 1, MIN(P,M-Q)-R
193 X(P-I+1,M-I+1) = X(P-I+1,M-I+1) + ONE
194 END DO
195 DO I = 1, R
196 X(P-(MIN(P,M-Q)-R)+1-I,M-(MIN(P,M-Q)-R)+1-I) =
197 $ X(P-(MIN(P,M-Q)-R)+1-I,M-(MIN(P,M-Q)-R)+1-I) +
198 $ SIN(THETA(R-I+1))
199 END DO
200 *
201 CALL SGEMM( 'No transpose', 'Conjugate transpose', M-P, Q, Q, ONE,
202 $ X(P+1,1), LDX, V1T, LDV1T, ZERO, WORK, LDX )
203 *
204 CALL SGEMM( 'Conjugate transpose', 'No transpose', M-P, Q, M-P,
205 $ ONE, U2, LDU2, WORK, LDX, ZERO, X(P+1,1), LDX )
206 *
207 DO I = 1, MIN(M-P,Q)-R
208 X(M-I+1,Q-I+1) = X(M-I+1,Q-I+1) - ONE
209 END DO
210 DO I = 1, R
211 X(M-(MIN(M-P,Q)-R)+1-I,Q-(MIN(M-P,Q)-R)+1-I) =
212 $ X(M-(MIN(M-P,Q)-R)+1-I,Q-(MIN(M-P,Q)-R)+1-I) -
213 $ SIN(THETA(R-I+1))
214 END DO
215 *
216 CALL SGEMM( 'No transpose', 'Conjugate transpose', M-P, M-Q, M-Q,
217 $ ONE, X(P+1,Q+1), LDX, V2T, LDV2T, ZERO, WORK, LDX )
218 *
219 CALL SGEMM( 'Conjugate transpose', 'No transpose', M-P, M-Q, M-P,
220 $ ONE, U2, LDU2, WORK, LDX, ZERO, X(P+1,Q+1), LDX )
221 *
222 DO I = 1, MIN(M-P,M-Q)-R
223 X(P+I,Q+I) = X(P+I,Q+I) - ONE
224 END DO
225 DO I = 1, R
226 X(P+(MIN(M-P,M-Q)-R)+I,Q+(MIN(M-P,M-Q)-R)+I) =
227 $ X(P+(MIN(M-P,M-Q)-R)+I,Q+(MIN(M-P,M-Q)-R)+I) -
228 $ COS(THETA(I))
229 END DO
230 *
231 * Compute norm( U1'*X11*V1 - D11 ) / ( MAX(1,P,Q)*EPS2 ) .
232 *
233 RESID = SLANGE( '1', P, Q, X, LDX, RWORK )
234 RESULT( 1 ) = ( RESID / REAL(MAX(1,P,Q)) ) / EPS2
235 *
236 * Compute norm( U1'*X12*V2 - D12 ) / ( MAX(1,P,M-Q)*EPS2 ) .
237 *
238 RESID = SLANGE( '1', P, M-Q, X(1,Q+1), LDX, RWORK )
239 RESULT( 2 ) = ( RESID / REAL(MAX(1,P,M-Q)) ) / EPS2
240 *
241 * Compute norm( U2'*X21*V1 - D21 ) / ( MAX(1,M-P,Q)*EPS2 ) .
242 *
243 RESID = SLANGE( '1', M-P, Q, X(P+1,1), LDX, RWORK )
244 RESULT( 3 ) = ( RESID / REAL(MAX(1,M-P,Q)) ) / EPS2
245 *
246 * Compute norm( U2'*X22*V2 - D22 ) / ( MAX(1,M-P,M-Q)*EPS2 ) .
247 *
248 RESID = SLANGE( '1', M-P, M-Q, X(P+1,Q+1), LDX, RWORK )
249 RESULT( 4 ) = ( RESID / REAL(MAX(1,M-P,M-Q)) ) / EPS2
250 *
251 * Compute I - U1'*U1
252 *
253 CALL SLASET( 'Full', P, P, ZERO, ONE, WORK, LDU1 )
254 CALL SSYRK( 'Upper', 'Conjugate transpose', P, P, -ONE, U1, LDU1,
255 $ ONE, WORK, LDU1 )
256 *
257 * Compute norm( I - U'*U ) / ( MAX(1,P) * ULP ) .
258 *
259 RESID = SLANSY( '1', 'Upper', P, WORK, LDU1, RWORK )
260 RESULT( 5 ) = ( RESID / REAL(MAX(1,P)) ) / ULP
261 *
262 * Compute I - U2'*U2
263 *
264 CALL SLASET( 'Full', M-P, M-P, ZERO, ONE, WORK, LDU2 )
265 CALL SSYRK( 'Upper', 'Conjugate transpose', M-P, M-P, -ONE, U2,
266 $ LDU2, ONE, WORK, LDU2 )
267 *
268 * Compute norm( I - U2'*U2 ) / ( MAX(1,M-P) * ULP ) .
269 *
270 RESID = SLANSY( '1', 'Upper', M-P, WORK, LDU2, RWORK )
271 RESULT( 6 ) = ( RESID / REAL(MAX(1,M-P)) ) / ULP
272 *
273 * Compute I - V1T*V1T'
274 *
275 CALL SLASET( 'Full', Q, Q, ZERO, ONE, WORK, LDV1T )
276 CALL SSYRK( 'Upper', 'No transpose', Q, Q, -ONE, V1T, LDV1T, ONE,
277 $ WORK, LDV1T )
278 *
279 * Compute norm( I - V1T*V1T' ) / ( MAX(1,Q) * ULP ) .
280 *
281 RESID = SLANSY( '1', 'Upper', Q, WORK, LDV1T, RWORK )
282 RESULT( 7 ) = ( RESID / REAL(MAX(1,Q)) ) / ULP
283 *
284 * Compute I - V2T*V2T'
285 *
286 CALL SLASET( 'Full', M-Q, M-Q, ZERO, ONE, WORK, LDV2T )
287 CALL SSYRK( 'Upper', 'No transpose', M-Q, M-Q, -ONE, V2T, LDV2T,
288 $ ONE, WORK, LDV2T )
289 *
290 * Compute norm( I - V2T*V2T' ) / ( MAX(1,M-Q) * ULP ) .
291 *
292 RESID = SLANSY( '1', 'Upper', M-Q, WORK, LDV2T, RWORK )
293 RESULT( 8 ) = ( RESID / REAL(MAX(1,M-Q)) ) / ULP
294 *
295 * Check sorting
296 *
297 RESULT(9) = REALZERO
298 DO I = 1, R
299 IF( THETA(I).LT.REALZERO .OR. THETA(I).GT.PIOVER2 ) THEN
300 RESULT(9) = ULPINV
301 END IF
302 IF( I.GT.1) THEN
303 IF ( THETA(I).LT.THETA(I-1) ) THEN
304 RESULT(9) = ULPINV
305 END IF
306 END IF
307 END DO
308 *
309 RETURN
310 *
311 * End of SCSDTS
312 *
313 END
314
2 $ LDV1T, V2T, LDV2T, THETA, IWORK, WORK, LWORK,
3 $ RWORK, RESULT )
4 IMPLICIT NONE
5 *
6 * Originally xGSVTS
7 * -- LAPACK test routine (version 3.1) --
8 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
9 * November 2006
10 *
11 * Adapted to SCSDTS
12 * July 2010
13 *
14 * .. Scalar Arguments ..
15 INTEGER LDX, LDU1, LDU2, LDV1T, LDV2T, LWORK, M, P, Q
16 * ..
17 * .. Array Arguments ..
18 INTEGER IWORK( * )
19 REAL RESULT( 9 ), RWORK( * ), THETA( * )
20 REAL U1( LDU1, * ), U2( LDU2, * ), V1T( LDV1T, * ),
21 $ V2T( LDV2T, * ), WORK( LWORK ), X( LDX, * ),
22 $ XF( LDX, * )
23 * ..
24 *
25 * Purpose
26 * =======
27 *
28 * SCSDTS tests SORCSD, which, given an M-by-M partitioned orthogonal
29 * matrix X,
30 * Q M-Q
31 * X = [ X11 X12 ] P ,
32 * [ X21 X22 ] M-P
33 *
34 * computes the CSD
35 *
36 * [ U1 ]**T * [ X11 X12 ] * [ V1 ]
37 * [ U2 ] [ X21 X22 ] [ V2 ]
38 *
39 * [ I 0 0 | 0 0 0 ]
40 * [ 0 C 0 | 0 -S 0 ]
41 * [ 0 0 0 | 0 0 -I ]
42 * = [---------------------] = [ D11 D12 ] .
43 * [ 0 0 0 | I 0 0 ] [ D21 D22 ]
44 * [ 0 S 0 | 0 C 0 ]
45 * [ 0 0 I | 0 0 0 ]
46 *
47 * Arguments
48 * =========
49 *
50 * M (input) INTEGER
51 * The number of rows of the matrix X. M >= 0.
52 *
53 * P (input) INTEGER
54 * The number of rows of the matrix X11. P >= 0.
55 *
56 * Q (input) INTEGER
57 * The number of columns of the matrix X11. Q >= 0.
58 *
59 * X (input) REAL array, dimension (LDX,M)
60 * The M-by-M matrix X.
61 *
62 * XF (output) REAL array, dimension (LDX,M)
63 * Details of the CSD of X, as returned by SORCSD;
64 * see SORCSD for further details.
65 *
66 * LDX (input) INTEGER
67 * The leading dimension of the arrays X and XF.
68 * LDX >= max( 1,M ).
69 *
70 * U1 (output) REAL array, dimension(LDU1,P)
71 * The P-by-P orthogonal matrix U1.
72 *
73 * LDU1 (input) INTEGER
74 * The leading dimension of the array U1. LDU >= max(1,P).
75 *
76 * U2 (output) REAL array, dimension(LDU2,M-P)
77 * The (M-P)-by-(M-P) orthogonal matrix U2.
78 *
79 * LDU2 (input) INTEGER
80 * The leading dimension of the array U2. LDU >= max(1,M-P).
81 *
82 * V1T (output) REAL array, dimension(LDV1T,Q)
83 * The Q-by-Q orthogonal matrix V1T.
84 *
85 * LDV1T (input) INTEGER
86 * The leading dimension of the array V1T. LDV1T >=
87 * max(1,Q).
88 *
89 * V2T (output) REAL array, dimension(LDV2T,M-Q)
90 * The (M-Q)-by-(M-Q) orthogonal matrix V2T.
91 *
92 * LDV2T (input) INTEGER
93 * The leading dimension of the array V2T. LDV2T >=
94 * max(1,M-Q).
95 *
96 * THETA (output) REAL array, dimension MIN(P,M-P,Q,M-Q)
97 * The CS values of X; the essentially diagonal matrices C and
98 * S are constructed from THETA; see subroutine SORCSD for
99 * details.
100 *
101 * IWORK (workspace) INTEGER array, dimension (M)
102 *
103 * WORK (workspace) REAL array, dimension (LWORK)
104 *
105 * LWORK (input) INTEGER
106 * The dimension of the array WORK
107 *
108 * RWORK (workspace) REAL array
109 *
110 * RESULT (output) REAL array, dimension (9)
111 * The test ratios:
112 * RESULT(1) = norm( U1'*X11*V1 - D11 ) / ( MAX(1,P,Q)*EPS2 )
113 * RESULT(2) = norm( U1'*X12*V2 - D12 ) / ( MAX(1,P,M-Q)*EPS2 )
114 * RESULT(3) = norm( U2'*X21*V1 - D21 ) / ( MAX(1,M-P,Q)*EPS2 )
115 * RESULT(4) = norm( U2'*X22*V2 - D22 ) / ( MAX(1,M-P,M-Q)*EPS2 )
116 * RESULT(5) = norm( I - U1'*U1 ) / ( MAX(1,P)*ULP )
117 * RESULT(6) = norm( I - U2'*U2 ) / ( MAX(1,M-P)*ULP )
118 * RESULT(7) = norm( I - V1T'*V1T ) / ( MAX(1,Q)*ULP )
119 * RESULT(8) = norm( I - V2T'*V2T ) / ( MAX(1,M-Q)*ULP )
120 * RESULT(9) = 0 if THETA is in increasing order and
121 * all angles are in [0,pi/2];
122 * = ULPINV otherwise.
123 * ( EPS2 = MAX( norm( I - X'*X ) / M, ULP ). )
124 *
125 * =====================================================================
126 *
127 * .. Parameters ..
128 REAL PIOVER2, REALONE, REALZERO
129 PARAMETER ( PIOVER2 = 1.57079632679489662E0,
130 $ REALONE = 1.0E0, REALZERO = 0.0E0 )
131 REAL ZERO, ONE
132 PARAMETER ( ZERO = 0.0E0, ONE = 1.0E0 )
133 * ..
134 * .. Local Scalars ..
135 INTEGER I, INFO, R
136 REAL EPS2, RESID, ULP, ULPINV
137 * ..
138 * .. External Functions ..
139 REAL SLAMCH, SLANGE, SLANSY
140 EXTERNAL SLAMCH, SLANGE, SLANSY
141 * ..
142 * .. External Subroutines ..
143 EXTERNAL SGEMM, SLACPY, SLASET, SORCSD, SSYRK
144 * ..
145 * .. Intrinsic Functions ..
146 INTRINSIC REAL, MAX, MIN
147 * ..
148 * .. Executable Statements ..
149 *
150 ULP = SLAMCH( 'Precision' )
151 ULPINV = REALONE / ULP
152 CALL SLASET( 'Full', M, M, ZERO, ONE, WORK, LDX )
153 CALL SSYRK( 'Upper', 'Conjugate transpose', M, M, -ONE, X, LDX,
154 $ ONE, WORK, LDX )
155 EPS2 = MAX( ULP,
156 $ SLANGE( '1', M, M, WORK, LDX, RWORK ) / REAL( M ) )
157 R = MIN( P, M-P, Q, M-Q )
158 *
159 * Copy the matrix X to the array XF.
160 *
161 CALL SLACPY( 'Full', M, M, X, LDX, XF, LDX )
162 *
163 * Compute the CSD
164 *
165 CALL SORCSD( 'Y', 'Y', 'Y', 'Y', 'N', 'D', M, P, Q, XF(1,1), LDX,
166 $ XF(1,Q+1), LDX, XF(P+1,1), LDX, XF(P+1,Q+1), LDX,
167 $ THETA, U1, LDU1, U2, LDU2, V1T, LDV1T, V2T, LDV2T,
168 $ WORK, LWORK, IWORK, INFO )
169 *
170 * Compute X := diag(U1,U2)'*X*diag(V1,V2) - [D11 D12; D21 D22]
171 *
172 CALL SGEMM( 'No transpose', 'Conjugate transpose', P, Q, Q, ONE,
173 $ X, LDX, V1T, LDV1T, ZERO, WORK, LDX )
174 *
175 CALL SGEMM( 'Conjugate transpose', 'No transpose', P, Q, P, ONE,
176 $ U1, LDU1, WORK, LDX, ZERO, X, LDX )
177 *
178 DO I = 1, MIN(P,Q)-R
179 X(I,I) = X(I,I) - ONE
180 END DO
181 DO I = 1, R
182 X(MIN(P,Q)-R+I,MIN(P,Q)-R+I) =
183 $ X(MIN(P,Q)-R+I,MIN(P,Q)-R+I) - COS(THETA(I))
184 END DO
185 *
186 CALL SGEMM( 'No transpose', 'Conjugate transpose', P, M-Q, M-Q,
187 $ ONE, X(1,Q+1), LDX, V2T, LDV2T, ZERO, WORK, LDX )
188 *
189 CALL SGEMM( 'Conjugate transpose', 'No transpose', P, M-Q, P,
190 $ ONE, U1, LDU1, WORK, LDX, ZERO, X(1,Q+1), LDX )
191 *
192 DO I = 1, MIN(P,M-Q)-R
193 X(P-I+1,M-I+1) = X(P-I+1,M-I+1) + ONE
194 END DO
195 DO I = 1, R
196 X(P-(MIN(P,M-Q)-R)+1-I,M-(MIN(P,M-Q)-R)+1-I) =
197 $ X(P-(MIN(P,M-Q)-R)+1-I,M-(MIN(P,M-Q)-R)+1-I) +
198 $ SIN(THETA(R-I+1))
199 END DO
200 *
201 CALL SGEMM( 'No transpose', 'Conjugate transpose', M-P, Q, Q, ONE,
202 $ X(P+1,1), LDX, V1T, LDV1T, ZERO, WORK, LDX )
203 *
204 CALL SGEMM( 'Conjugate transpose', 'No transpose', M-P, Q, M-P,
205 $ ONE, U2, LDU2, WORK, LDX, ZERO, X(P+1,1), LDX )
206 *
207 DO I = 1, MIN(M-P,Q)-R
208 X(M-I+1,Q-I+1) = X(M-I+1,Q-I+1) - ONE
209 END DO
210 DO I = 1, R
211 X(M-(MIN(M-P,Q)-R)+1-I,Q-(MIN(M-P,Q)-R)+1-I) =
212 $ X(M-(MIN(M-P,Q)-R)+1-I,Q-(MIN(M-P,Q)-R)+1-I) -
213 $ SIN(THETA(R-I+1))
214 END DO
215 *
216 CALL SGEMM( 'No transpose', 'Conjugate transpose', M-P, M-Q, M-Q,
217 $ ONE, X(P+1,Q+1), LDX, V2T, LDV2T, ZERO, WORK, LDX )
218 *
219 CALL SGEMM( 'Conjugate transpose', 'No transpose', M-P, M-Q, M-P,
220 $ ONE, U2, LDU2, WORK, LDX, ZERO, X(P+1,Q+1), LDX )
221 *
222 DO I = 1, MIN(M-P,M-Q)-R
223 X(P+I,Q+I) = X(P+I,Q+I) - ONE
224 END DO
225 DO I = 1, R
226 X(P+(MIN(M-P,M-Q)-R)+I,Q+(MIN(M-P,M-Q)-R)+I) =
227 $ X(P+(MIN(M-P,M-Q)-R)+I,Q+(MIN(M-P,M-Q)-R)+I) -
228 $ COS(THETA(I))
229 END DO
230 *
231 * Compute norm( U1'*X11*V1 - D11 ) / ( MAX(1,P,Q)*EPS2 ) .
232 *
233 RESID = SLANGE( '1', P, Q, X, LDX, RWORK )
234 RESULT( 1 ) = ( RESID / REAL(MAX(1,P,Q)) ) / EPS2
235 *
236 * Compute norm( U1'*X12*V2 - D12 ) / ( MAX(1,P,M-Q)*EPS2 ) .
237 *
238 RESID = SLANGE( '1', P, M-Q, X(1,Q+1), LDX, RWORK )
239 RESULT( 2 ) = ( RESID / REAL(MAX(1,P,M-Q)) ) / EPS2
240 *
241 * Compute norm( U2'*X21*V1 - D21 ) / ( MAX(1,M-P,Q)*EPS2 ) .
242 *
243 RESID = SLANGE( '1', M-P, Q, X(P+1,1), LDX, RWORK )
244 RESULT( 3 ) = ( RESID / REAL(MAX(1,M-P,Q)) ) / EPS2
245 *
246 * Compute norm( U2'*X22*V2 - D22 ) / ( MAX(1,M-P,M-Q)*EPS2 ) .
247 *
248 RESID = SLANGE( '1', M-P, M-Q, X(P+1,Q+1), LDX, RWORK )
249 RESULT( 4 ) = ( RESID / REAL(MAX(1,M-P,M-Q)) ) / EPS2
250 *
251 * Compute I - U1'*U1
252 *
253 CALL SLASET( 'Full', P, P, ZERO, ONE, WORK, LDU1 )
254 CALL SSYRK( 'Upper', 'Conjugate transpose', P, P, -ONE, U1, LDU1,
255 $ ONE, WORK, LDU1 )
256 *
257 * Compute norm( I - U'*U ) / ( MAX(1,P) * ULP ) .
258 *
259 RESID = SLANSY( '1', 'Upper', P, WORK, LDU1, RWORK )
260 RESULT( 5 ) = ( RESID / REAL(MAX(1,P)) ) / ULP
261 *
262 * Compute I - U2'*U2
263 *
264 CALL SLASET( 'Full', M-P, M-P, ZERO, ONE, WORK, LDU2 )
265 CALL SSYRK( 'Upper', 'Conjugate transpose', M-P, M-P, -ONE, U2,
266 $ LDU2, ONE, WORK, LDU2 )
267 *
268 * Compute norm( I - U2'*U2 ) / ( MAX(1,M-P) * ULP ) .
269 *
270 RESID = SLANSY( '1', 'Upper', M-P, WORK, LDU2, RWORK )
271 RESULT( 6 ) = ( RESID / REAL(MAX(1,M-P)) ) / ULP
272 *
273 * Compute I - V1T*V1T'
274 *
275 CALL SLASET( 'Full', Q, Q, ZERO, ONE, WORK, LDV1T )
276 CALL SSYRK( 'Upper', 'No transpose', Q, Q, -ONE, V1T, LDV1T, ONE,
277 $ WORK, LDV1T )
278 *
279 * Compute norm( I - V1T*V1T' ) / ( MAX(1,Q) * ULP ) .
280 *
281 RESID = SLANSY( '1', 'Upper', Q, WORK, LDV1T, RWORK )
282 RESULT( 7 ) = ( RESID / REAL(MAX(1,Q)) ) / ULP
283 *
284 * Compute I - V2T*V2T'
285 *
286 CALL SLASET( 'Full', M-Q, M-Q, ZERO, ONE, WORK, LDV2T )
287 CALL SSYRK( 'Upper', 'No transpose', M-Q, M-Q, -ONE, V2T, LDV2T,
288 $ ONE, WORK, LDV2T )
289 *
290 * Compute norm( I - V2T*V2T' ) / ( MAX(1,M-Q) * ULP ) .
291 *
292 RESID = SLANSY( '1', 'Upper', M-Q, WORK, LDV2T, RWORK )
293 RESULT( 8 ) = ( RESID / REAL(MAX(1,M-Q)) ) / ULP
294 *
295 * Check sorting
296 *
297 RESULT(9) = REALZERO
298 DO I = 1, R
299 IF( THETA(I).LT.REALZERO .OR. THETA(I).GT.PIOVER2 ) THEN
300 RESULT(9) = ULPINV
301 END IF
302 IF( I.GT.1) THEN
303 IF ( THETA(I).LT.THETA(I-1) ) THEN
304 RESULT(9) = ULPINV
305 END IF
306 END IF
307 END DO
308 *
309 RETURN
310 *
311 * End of SCSDTS
312 *
313 END
314