1 SUBROUTINE ZCSDTS( 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.3.0) --
8 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
9 * November 2010
10 *
11 * Adapted to ZCSDTS
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 DOUBLE PRECISION RESULT( 9 ), RWORK( * ), THETA( * )
20 COMPLEX*16 U1( LDU1, * ), U2( LDU2, * ), V1T( LDV1T, * ),
21 $ V2T( LDV2T, * ), WORK( LWORK ), X( LDX, * ),
22 $ XF( LDX, * )
23 * ..
24 *
25 * Purpose
26 * =======
27 *
28 * ZCSDTS tests ZUNCSD, which, given an M-by-M partitioned unitary
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) COMPLEX*16 array, dimension (LDX,M)
60 * The M-by-M matrix X.
61 *
62 * XF (output) COMPLEX*16 array, dimension (LDX,M)
63 * Details of the CSD of X, as returned by ZUNCSD;
64 * see ZUNCSD 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) COMPLEX*16 array, dimension(LDU1,P)
71 * The P-by-P unitary matrix U1.
72 *
73 * LDU1 (input) INTEGER
74 * The leading dimension of the array U1. LDU >= max(1,P).
75 *
76 * U2 (output) COMPLEX*16 array, dimension(LDU2,M-P)
77 * The (M-P)-by-(M-P) unitary matrix U2.
78 *
79 * LDU2 (input) INTEGER
80 * The leading dimension of the array U2. LDU >= max(1,M-P).
81 *
82 * V1T (output) COMPLEX*16 array, dimension(LDV1T,Q)
83 * The Q-by-Q unitary matrix V1T.
84 *
85 * LDV1T (input) INTEGER
86 * The leading dimension of the array V1T. LDV1T >=
87 * max(1,Q).
88 *
89 * V2T (output) COMPLEX*16 array, dimension(LDV2T,M-Q)
90 * The (M-Q)-by-(M-Q) unitary 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) DOUBLE PRECISION 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 ZUNCSD for
99 * details.
100 *
101 * IWORK (workspace) INTEGER array, dimension (M)
102 *
103 * WORK (workspace) COMPLEX*16 array, dimension (LWORK)
104 *
105 * LWORK (input) INTEGER
106 * The dimension of the array WORK
107 *
108 * RWORK (workspace) DOUBLE PRECISION array
109 *
110 * RESULT (output) DOUBLE PRECISION 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 DOUBLE PRECISION PIOVER2, REALONE, REALZERO
129 PARAMETER ( PIOVER2 = 1.57079632679489662D0,
130 $ REALONE = 1.0D0, REALZERO = 0.0D0 )
131 COMPLEX*16 ZERO, ONE
132 PARAMETER ( ZERO = (0.0D0,0.0D0), ONE = (1.0D0,0.0D0) )
133 * ..
134 * .. Local Scalars ..
135 INTEGER I, INFO, R
136 DOUBLE PRECISION EPS2, RESID, ULP, ULPINV
137 * ..
138 * .. External Functions ..
139 DOUBLE PRECISION DLAMCH, ZLANGE, ZLANHE
140 EXTERNAL DLAMCH, ZLANGE, ZLANHE
141 * ..
142 * .. External Subroutines ..
143 EXTERNAL ZGEMM, ZLACPY, ZLASET, ZUNCSD, ZHERK
144 * ..
145 * .. Intrinsic Functions ..
146 INTRINSIC REAL, MAX, MIN
147 * ..
148 * .. Executable Statements ..
149 *
150 ULP = DLAMCH( 'Precision' )
151 ULPINV = REALONE / ULP
152 CALL ZLASET( 'Full', M, M, ZERO, ONE, WORK, LDX )
153 CALL ZHERK( 'Upper', 'Conjugate transpose', M, M, -ONE, X, LDX,
154 $ ONE, WORK, LDX )
155 EPS2 = MAX( ULP,
156 $ ZLANGE( '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 ZLACPY( 'Full', M, M, X, LDX, XF, LDX )
162 *
163 * Compute the CSD
164 *
165 CALL ZUNCSD( '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, RWORK, 17*(R+2), IWORK, INFO )
169 *
170 * Compute X := diag(U1,U2)'*X*diag(V1,V2) - [D11 D12; D21 D22]
171 *
172 CALL ZGEMM( 'No transpose', 'Conjugate transpose', P, Q, Q, ONE,
173 $ X, LDX, V1T, LDV1T, ZERO, WORK, LDX )
174 *
175 CALL ZGEMM( '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) - DCMPLX( COS(THETA(I)),
184 $ 0.0D0 )
185 END DO
186 *
187 CALL ZGEMM( 'No transpose', 'Conjugate transpose', P, M-Q, M-Q,
188 $ ONE, X(1,Q+1), LDX, V2T, LDV2T, ZERO, WORK, LDX )
189 *
190 CALL ZGEMM( 'Conjugate transpose', 'No transpose', P, M-Q, P,
191 $ ONE, U1, LDU1, WORK, LDX, ZERO, X(1,Q+1), LDX )
192 *
193 DO I = 1, MIN(P,M-Q)-R
194 X(P-I+1,M-I+1) = X(P-I+1,M-I+1) + ONE
195 END DO
196 DO I = 1, R
197 X(P-(MIN(P,M-Q)-R)+1-I,M-(MIN(P,M-Q)-R)+1-I) =
198 $ X(P-(MIN(P,M-Q)-R)+1-I,M-(MIN(P,M-Q)-R)+1-I) +
199 $ DCMPLX( SIN(THETA(R-I+1)), 0.0D0 )
200 END DO
201 *
202 CALL ZGEMM( 'No transpose', 'Conjugate transpose', M-P, Q, Q, ONE,
203 $ X(P+1,1), LDX, V1T, LDV1T, ZERO, WORK, LDX )
204 *
205 CALL ZGEMM( 'Conjugate transpose', 'No transpose', M-P, Q, M-P,
206 $ ONE, U2, LDU2, WORK, LDX, ZERO, X(P+1,1), LDX )
207 *
208 DO I = 1, MIN(M-P,Q)-R
209 X(M-I+1,Q-I+1) = X(M-I+1,Q-I+1) - ONE
210 END DO
211 DO I = 1, R
212 X(M-(MIN(M-P,Q)-R)+1-I,Q-(MIN(M-P,Q)-R)+1-I) =
213 $ X(M-(MIN(M-P,Q)-R)+1-I,Q-(MIN(M-P,Q)-R)+1-I) -
214 $ DCMPLX( SIN(THETA(R-I+1)), 0.0D0 )
215 END DO
216 *
217 CALL ZGEMM( 'No transpose', 'Conjugate transpose', M-P, M-Q, M-Q,
218 $ ONE, X(P+1,Q+1), LDX, V2T, LDV2T, ZERO, WORK, LDX )
219 *
220 CALL ZGEMM( 'Conjugate transpose', 'No transpose', M-P, M-Q, M-P,
221 $ ONE, U2, LDU2, WORK, LDX, ZERO, X(P+1,Q+1), LDX )
222 *
223 DO I = 1, MIN(M-P,M-Q)-R
224 X(P+I,Q+I) = X(P+I,Q+I) - ONE
225 END DO
226 DO I = 1, R
227 X(P+(MIN(M-P,M-Q)-R)+I,Q+(MIN(M-P,M-Q)-R)+I) =
228 $ X(P+(MIN(M-P,M-Q)-R)+I,Q+(MIN(M-P,M-Q)-R)+I) -
229 $ DCMPLX( COS(THETA(I)), 0.0D0 )
230 END DO
231 *
232 * Compute norm( U1'*X11*V1 - D11 ) / ( MAX(1,P,Q)*EPS2 ) .
233 *
234 RESID = ZLANGE( '1', P, Q, X, LDX, RWORK )
235 RESULT( 1 ) = ( RESID / REAL(MAX(1,P,Q)) ) / EPS2
236 *
237 * Compute norm( U1'*X12*V2 - D12 ) / ( MAX(1,P,M-Q)*EPS2 ) .
238 *
239 RESID = ZLANGE( '1', P, M-Q, X(1,Q+1), LDX, RWORK )
240 RESULT( 2 ) = ( RESID / REAL(MAX(1,P,M-Q)) ) / EPS2
241 *
242 * Compute norm( U2'*X21*V1 - D21 ) / ( MAX(1,M-P,Q)*EPS2 ) .
243 *
244 RESID = ZLANGE( '1', M-P, Q, X(P+1,1), LDX, RWORK )
245 RESULT( 3 ) = ( RESID / REAL(MAX(1,M-P,Q)) ) / EPS2
246 *
247 * Compute norm( U2'*X22*V2 - D22 ) / ( MAX(1,M-P,M-Q)*EPS2 ) .
248 *
249 RESID = ZLANGE( '1', M-P, M-Q, X(P+1,Q+1), LDX, RWORK )
250 RESULT( 4 ) = ( RESID / REAL(MAX(1,M-P,M-Q)) ) / EPS2
251 *
252 * Compute I - U1'*U1
253 *
254 CALL ZLASET( 'Full', P, P, ZERO, ONE, WORK, LDU1 )
255 CALL ZHERK( 'Upper', 'Conjugate transpose', P, P, -ONE, U1, LDU1,
256 $ ONE, WORK, LDU1 )
257 *
258 * Compute norm( I - U'*U ) / ( MAX(1,P) * ULP ) .
259 *
260 RESID = ZLANHE( '1', 'Upper', P, WORK, LDU1, RWORK )
261 RESULT( 5 ) = ( RESID / REAL(MAX(1,P)) ) / ULP
262 *
263 * Compute I - U2'*U2
264 *
265 CALL ZLASET( 'Full', M-P, M-P, ZERO, ONE, WORK, LDU2 )
266 CALL ZHERK( 'Upper', 'Conjugate transpose', M-P, M-P, -ONE, U2,
267 $ LDU2, ONE, WORK, LDU2 )
268 *
269 * Compute norm( I - U2'*U2 ) / ( MAX(1,M-P) * ULP ) .
270 *
271 RESID = ZLANHE( '1', 'Upper', M-P, WORK, LDU2, RWORK )
272 RESULT( 6 ) = ( RESID / REAL(MAX(1,M-P)) ) / ULP
273 *
274 * Compute I - V1T*V1T'
275 *
276 CALL ZLASET( 'Full', Q, Q, ZERO, ONE, WORK, LDV1T )
277 CALL ZHERK( 'Upper', 'No transpose', Q, Q, -ONE, V1T, LDV1T, ONE,
278 $ WORK, LDV1T )
279 *
280 * Compute norm( I - V1T*V1T' ) / ( MAX(1,Q) * ULP ) .
281 *
282 RESID = ZLANHE( '1', 'Upper', Q, WORK, LDV1T, RWORK )
283 RESULT( 7 ) = ( RESID / REAL(MAX(1,Q)) ) / ULP
284 *
285 * Compute I - V2T*V2T'
286 *
287 CALL ZLASET( 'Full', M-Q, M-Q, ZERO, ONE, WORK, LDV2T )
288 CALL ZHERK( 'Upper', 'No transpose', M-Q, M-Q, -ONE, V2T, LDV2T,
289 $ ONE, WORK, LDV2T )
290 *
291 * Compute norm( I - V2T*V2T' ) / ( MAX(1,M-Q) * ULP ) .
292 *
293 RESID = ZLANHE( '1', 'Upper', M-Q, WORK, LDV2T, RWORK )
294 RESULT( 8 ) = ( RESID / REAL(MAX(1,M-Q)) ) / ULP
295 *
296 * Check sorting
297 *
298 RESULT(9) = REALZERO
299 DO I = 1, R
300 IF( THETA(I).LT.REALZERO .OR. THETA(I).GT.PIOVER2 ) THEN
301 RESULT(9) = ULPINV
302 END IF
303 IF( I.GT.1) THEN
304 IF ( THETA(I).LT.THETA(I-1) ) THEN
305 RESULT(9) = ULPINV
306 END IF
307 END IF
308 END DO
309 *
310 RETURN
311 *
312 * End of ZCSDTS
313 *
314 END
315
2 $ LDV1T, V2T, LDV2T, THETA, IWORK, WORK, LWORK,
3 $ RWORK, RESULT )
4 IMPLICIT NONE
5 *
6 * Originally xGSVTS
7 * -- LAPACK test routine (version 3.3.0) --
8 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
9 * November 2010
10 *
11 * Adapted to ZCSDTS
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 DOUBLE PRECISION RESULT( 9 ), RWORK( * ), THETA( * )
20 COMPLEX*16 U1( LDU1, * ), U2( LDU2, * ), V1T( LDV1T, * ),
21 $ V2T( LDV2T, * ), WORK( LWORK ), X( LDX, * ),
22 $ XF( LDX, * )
23 * ..
24 *
25 * Purpose
26 * =======
27 *
28 * ZCSDTS tests ZUNCSD, which, given an M-by-M partitioned unitary
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) COMPLEX*16 array, dimension (LDX,M)
60 * The M-by-M matrix X.
61 *
62 * XF (output) COMPLEX*16 array, dimension (LDX,M)
63 * Details of the CSD of X, as returned by ZUNCSD;
64 * see ZUNCSD 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) COMPLEX*16 array, dimension(LDU1,P)
71 * The P-by-P unitary matrix U1.
72 *
73 * LDU1 (input) INTEGER
74 * The leading dimension of the array U1. LDU >= max(1,P).
75 *
76 * U2 (output) COMPLEX*16 array, dimension(LDU2,M-P)
77 * The (M-P)-by-(M-P) unitary matrix U2.
78 *
79 * LDU2 (input) INTEGER
80 * The leading dimension of the array U2. LDU >= max(1,M-P).
81 *
82 * V1T (output) COMPLEX*16 array, dimension(LDV1T,Q)
83 * The Q-by-Q unitary matrix V1T.
84 *
85 * LDV1T (input) INTEGER
86 * The leading dimension of the array V1T. LDV1T >=
87 * max(1,Q).
88 *
89 * V2T (output) COMPLEX*16 array, dimension(LDV2T,M-Q)
90 * The (M-Q)-by-(M-Q) unitary 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) DOUBLE PRECISION 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 ZUNCSD for
99 * details.
100 *
101 * IWORK (workspace) INTEGER array, dimension (M)
102 *
103 * WORK (workspace) COMPLEX*16 array, dimension (LWORK)
104 *
105 * LWORK (input) INTEGER
106 * The dimension of the array WORK
107 *
108 * RWORK (workspace) DOUBLE PRECISION array
109 *
110 * RESULT (output) DOUBLE PRECISION 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 DOUBLE PRECISION PIOVER2, REALONE, REALZERO
129 PARAMETER ( PIOVER2 = 1.57079632679489662D0,
130 $ REALONE = 1.0D0, REALZERO = 0.0D0 )
131 COMPLEX*16 ZERO, ONE
132 PARAMETER ( ZERO = (0.0D0,0.0D0), ONE = (1.0D0,0.0D0) )
133 * ..
134 * .. Local Scalars ..
135 INTEGER I, INFO, R
136 DOUBLE PRECISION EPS2, RESID, ULP, ULPINV
137 * ..
138 * .. External Functions ..
139 DOUBLE PRECISION DLAMCH, ZLANGE, ZLANHE
140 EXTERNAL DLAMCH, ZLANGE, ZLANHE
141 * ..
142 * .. External Subroutines ..
143 EXTERNAL ZGEMM, ZLACPY, ZLASET, ZUNCSD, ZHERK
144 * ..
145 * .. Intrinsic Functions ..
146 INTRINSIC REAL, MAX, MIN
147 * ..
148 * .. Executable Statements ..
149 *
150 ULP = DLAMCH( 'Precision' )
151 ULPINV = REALONE / ULP
152 CALL ZLASET( 'Full', M, M, ZERO, ONE, WORK, LDX )
153 CALL ZHERK( 'Upper', 'Conjugate transpose', M, M, -ONE, X, LDX,
154 $ ONE, WORK, LDX )
155 EPS2 = MAX( ULP,
156 $ ZLANGE( '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 ZLACPY( 'Full', M, M, X, LDX, XF, LDX )
162 *
163 * Compute the CSD
164 *
165 CALL ZUNCSD( '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, RWORK, 17*(R+2), IWORK, INFO )
169 *
170 * Compute X := diag(U1,U2)'*X*diag(V1,V2) - [D11 D12; D21 D22]
171 *
172 CALL ZGEMM( 'No transpose', 'Conjugate transpose', P, Q, Q, ONE,
173 $ X, LDX, V1T, LDV1T, ZERO, WORK, LDX )
174 *
175 CALL ZGEMM( '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) - DCMPLX( COS(THETA(I)),
184 $ 0.0D0 )
185 END DO
186 *
187 CALL ZGEMM( 'No transpose', 'Conjugate transpose', P, M-Q, M-Q,
188 $ ONE, X(1,Q+1), LDX, V2T, LDV2T, ZERO, WORK, LDX )
189 *
190 CALL ZGEMM( 'Conjugate transpose', 'No transpose', P, M-Q, P,
191 $ ONE, U1, LDU1, WORK, LDX, ZERO, X(1,Q+1), LDX )
192 *
193 DO I = 1, MIN(P,M-Q)-R
194 X(P-I+1,M-I+1) = X(P-I+1,M-I+1) + ONE
195 END DO
196 DO I = 1, R
197 X(P-(MIN(P,M-Q)-R)+1-I,M-(MIN(P,M-Q)-R)+1-I) =
198 $ X(P-(MIN(P,M-Q)-R)+1-I,M-(MIN(P,M-Q)-R)+1-I) +
199 $ DCMPLX( SIN(THETA(R-I+1)), 0.0D0 )
200 END DO
201 *
202 CALL ZGEMM( 'No transpose', 'Conjugate transpose', M-P, Q, Q, ONE,
203 $ X(P+1,1), LDX, V1T, LDV1T, ZERO, WORK, LDX )
204 *
205 CALL ZGEMM( 'Conjugate transpose', 'No transpose', M-P, Q, M-P,
206 $ ONE, U2, LDU2, WORK, LDX, ZERO, X(P+1,1), LDX )
207 *
208 DO I = 1, MIN(M-P,Q)-R
209 X(M-I+1,Q-I+1) = X(M-I+1,Q-I+1) - ONE
210 END DO
211 DO I = 1, R
212 X(M-(MIN(M-P,Q)-R)+1-I,Q-(MIN(M-P,Q)-R)+1-I) =
213 $ X(M-(MIN(M-P,Q)-R)+1-I,Q-(MIN(M-P,Q)-R)+1-I) -
214 $ DCMPLX( SIN(THETA(R-I+1)), 0.0D0 )
215 END DO
216 *
217 CALL ZGEMM( 'No transpose', 'Conjugate transpose', M-P, M-Q, M-Q,
218 $ ONE, X(P+1,Q+1), LDX, V2T, LDV2T, ZERO, WORK, LDX )
219 *
220 CALL ZGEMM( 'Conjugate transpose', 'No transpose', M-P, M-Q, M-P,
221 $ ONE, U2, LDU2, WORK, LDX, ZERO, X(P+1,Q+1), LDX )
222 *
223 DO I = 1, MIN(M-P,M-Q)-R
224 X(P+I,Q+I) = X(P+I,Q+I) - ONE
225 END DO
226 DO I = 1, R
227 X(P+(MIN(M-P,M-Q)-R)+I,Q+(MIN(M-P,M-Q)-R)+I) =
228 $ X(P+(MIN(M-P,M-Q)-R)+I,Q+(MIN(M-P,M-Q)-R)+I) -
229 $ DCMPLX( COS(THETA(I)), 0.0D0 )
230 END DO
231 *
232 * Compute norm( U1'*X11*V1 - D11 ) / ( MAX(1,P,Q)*EPS2 ) .
233 *
234 RESID = ZLANGE( '1', P, Q, X, LDX, RWORK )
235 RESULT( 1 ) = ( RESID / REAL(MAX(1,P,Q)) ) / EPS2
236 *
237 * Compute norm( U1'*X12*V2 - D12 ) / ( MAX(1,P,M-Q)*EPS2 ) .
238 *
239 RESID = ZLANGE( '1', P, M-Q, X(1,Q+1), LDX, RWORK )
240 RESULT( 2 ) = ( RESID / REAL(MAX(1,P,M-Q)) ) / EPS2
241 *
242 * Compute norm( U2'*X21*V1 - D21 ) / ( MAX(1,M-P,Q)*EPS2 ) .
243 *
244 RESID = ZLANGE( '1', M-P, Q, X(P+1,1), LDX, RWORK )
245 RESULT( 3 ) = ( RESID / REAL(MAX(1,M-P,Q)) ) / EPS2
246 *
247 * Compute norm( U2'*X22*V2 - D22 ) / ( MAX(1,M-P,M-Q)*EPS2 ) .
248 *
249 RESID = ZLANGE( '1', M-P, M-Q, X(P+1,Q+1), LDX, RWORK )
250 RESULT( 4 ) = ( RESID / REAL(MAX(1,M-P,M-Q)) ) / EPS2
251 *
252 * Compute I - U1'*U1
253 *
254 CALL ZLASET( 'Full', P, P, ZERO, ONE, WORK, LDU1 )
255 CALL ZHERK( 'Upper', 'Conjugate transpose', P, P, -ONE, U1, LDU1,
256 $ ONE, WORK, LDU1 )
257 *
258 * Compute norm( I - U'*U ) / ( MAX(1,P) * ULP ) .
259 *
260 RESID = ZLANHE( '1', 'Upper', P, WORK, LDU1, RWORK )
261 RESULT( 5 ) = ( RESID / REAL(MAX(1,P)) ) / ULP
262 *
263 * Compute I - U2'*U2
264 *
265 CALL ZLASET( 'Full', M-P, M-P, ZERO, ONE, WORK, LDU2 )
266 CALL ZHERK( 'Upper', 'Conjugate transpose', M-P, M-P, -ONE, U2,
267 $ LDU2, ONE, WORK, LDU2 )
268 *
269 * Compute norm( I - U2'*U2 ) / ( MAX(1,M-P) * ULP ) .
270 *
271 RESID = ZLANHE( '1', 'Upper', M-P, WORK, LDU2, RWORK )
272 RESULT( 6 ) = ( RESID / REAL(MAX(1,M-P)) ) / ULP
273 *
274 * Compute I - V1T*V1T'
275 *
276 CALL ZLASET( 'Full', Q, Q, ZERO, ONE, WORK, LDV1T )
277 CALL ZHERK( 'Upper', 'No transpose', Q, Q, -ONE, V1T, LDV1T, ONE,
278 $ WORK, LDV1T )
279 *
280 * Compute norm( I - V1T*V1T' ) / ( MAX(1,Q) * ULP ) .
281 *
282 RESID = ZLANHE( '1', 'Upper', Q, WORK, LDV1T, RWORK )
283 RESULT( 7 ) = ( RESID / REAL(MAX(1,Q)) ) / ULP
284 *
285 * Compute I - V2T*V2T'
286 *
287 CALL ZLASET( 'Full', M-Q, M-Q, ZERO, ONE, WORK, LDV2T )
288 CALL ZHERK( 'Upper', 'No transpose', M-Q, M-Q, -ONE, V2T, LDV2T,
289 $ ONE, WORK, LDV2T )
290 *
291 * Compute norm( I - V2T*V2T' ) / ( MAX(1,M-Q) * ULP ) .
292 *
293 RESID = ZLANHE( '1', 'Upper', M-Q, WORK, LDV2T, RWORK )
294 RESULT( 8 ) = ( RESID / REAL(MAX(1,M-Q)) ) / ULP
295 *
296 * Check sorting
297 *
298 RESULT(9) = REALZERO
299 DO I = 1, R
300 IF( THETA(I).LT.REALZERO .OR. THETA(I).GT.PIOVER2 ) THEN
301 RESULT(9) = ULPINV
302 END IF
303 IF( I.GT.1) THEN
304 IF ( THETA(I).LT.THETA(I-1) ) THEN
305 RESULT(9) = ULPINV
306 END IF
307 END IF
308 END DO
309 *
310 RETURN
311 *
312 * End of ZCSDTS
313 *
314 END
315