LAPACK  3.10.1
LAPACK: Linear Algebra PACKage
zcsdts.f
Go to the documentation of this file.
1 *> \brief \b ZCSDTS
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 * Definition:
9 * ===========
10 *
11 * SUBROUTINE ZCSDTS( M, P, Q, X, XF, LDX, U1, LDU1, U2, LDU2, V1T,
12 * LDV1T, V2T, LDV2T, THETA, IWORK, WORK, LWORK,
13 * RWORK, RESULT )
14 *
15 * .. Scalar Arguments ..
16 * INTEGER LDX, LDU1, LDU2, LDV1T, LDV2T, LWORK, M, P, Q
17 * ..
18 * .. Array Arguments ..
19 * INTEGER IWORK( * )
20 * DOUBLE PRECISION RESULT( 15 ), RWORK( * ), THETA( * )
21 * COMPLEX*16 U1( LDU1, * ), U2( LDU2, * ), V1T( LDV1T, * ),
22 * $ V2T( LDV2T, * ), WORK( LWORK ), X( LDX, * ),
23 * $ XF( LDX, * )
24 * ..
25 *
26 *
27 *> \par Purpose:
28 * =============
29 *>
30 *> \verbatim
31 *>
32 *> ZCSDTS tests ZUNCSD, which, given an M-by-M partitioned unitary
33 *> matrix X,
34 *> Q M-Q
35 *> X = [ X11 X12 ] P ,
36 *> [ X21 X22 ] M-P
37 *>
38 *> computes the CSD
39 *>
40 *> [ U1 ]**T * [ X11 X12 ] * [ V1 ]
41 *> [ U2 ] [ X21 X22 ] [ V2 ]
42 *>
43 *> [ I 0 0 | 0 0 0 ]
44 *> [ 0 C 0 | 0 -S 0 ]
45 *> [ 0 0 0 | 0 0 -I ]
46 *> = [---------------------] = [ D11 D12 ] .
47 *> [ 0 0 0 | I 0 0 ] [ D21 D22 ]
48 *> [ 0 S 0 | 0 C 0 ]
49 *> [ 0 0 I | 0 0 0 ]
50 *>
51 *> and also SORCSD2BY1, which, given
52 *> Q
53 *> [ X11 ] P ,
54 *> [ X21 ] M-P
55 *>
56 *> computes the 2-by-1 CSD
57 *>
58 *> [ I 0 0 ]
59 *> [ 0 C 0 ]
60 *> [ 0 0 0 ]
61 *> [ U1 ]**T * [ X11 ] * V1 = [----------] = [ D11 ] ,
62 *> [ U2 ] [ X21 ] [ 0 0 0 ] [ D21 ]
63 *> [ 0 S 0 ]
64 *> [ 0 0 I ]
65 *> \endverbatim
66 *
67 * Arguments:
68 * ==========
69 *
70 *> \param[in] M
71 *> \verbatim
72 *> M is INTEGER
73 *> The number of rows of the matrix X. M >= 0.
74 *> \endverbatim
75 *>
76 *> \param[in] P
77 *> \verbatim
78 *> P is INTEGER
79 *> The number of rows of the matrix X11. P >= 0.
80 *> \endverbatim
81 *>
82 *> \param[in] Q
83 *> \verbatim
84 *> Q is INTEGER
85 *> The number of columns of the matrix X11. Q >= 0.
86 *> \endverbatim
87 *>
88 *> \param[in] X
89 *> \verbatim
90 *> X is COMPLEX*16 array, dimension (LDX,M)
91 *> The M-by-M matrix X.
92 *> \endverbatim
93 *>
94 *> \param[out] XF
95 *> \verbatim
96 *> XF is COMPLEX*16 array, dimension (LDX,M)
97 *> Details of the CSD of X, as returned by ZUNCSD;
98 *> see ZUNCSD for further details.
99 *> \endverbatim
100 *>
101 *> \param[in] LDX
102 *> \verbatim
103 *> LDX is INTEGER
104 *> The leading dimension of the arrays X and XF.
105 *> LDX >= max( 1,M ).
106 *> \endverbatim
107 *>
108 *> \param[out] U1
109 *> \verbatim
110 *> U1 is COMPLEX*16 array, dimension(LDU1,P)
111 *> The P-by-P unitary matrix U1.
112 *> \endverbatim
113 *>
114 *> \param[in] LDU1
115 *> \verbatim
116 *> LDU1 is INTEGER
117 *> The leading dimension of the array U1. LDU >= max(1,P).
118 *> \endverbatim
119 *>
120 *> \param[out] U2
121 *> \verbatim
122 *> U2 is COMPLEX*16 array, dimension(LDU2,M-P)
123 *> The (M-P)-by-(M-P) unitary matrix U2.
124 *> \endverbatim
125 *>
126 *> \param[in] LDU2
127 *> \verbatim
128 *> LDU2 is INTEGER
129 *> The leading dimension of the array U2. LDU >= max(1,M-P).
130 *> \endverbatim
131 *>
132 *> \param[out] V1T
133 *> \verbatim
134 *> V1T is COMPLEX*16 array, dimension(LDV1T,Q)
135 *> The Q-by-Q unitary matrix V1T.
136 *> \endverbatim
137 *>
138 *> \param[in] LDV1T
139 *> \verbatim
140 *> LDV1T is INTEGER
141 *> The leading dimension of the array V1T. LDV1T >=
142 *> max(1,Q).
143 *> \endverbatim
144 *>
145 *> \param[out] V2T
146 *> \verbatim
147 *> V2T is COMPLEX*16 array, dimension(LDV2T,M-Q)
148 *> The (M-Q)-by-(M-Q) unitary matrix V2T.
149 *> \endverbatim
150 *>
151 *> \param[in] LDV2T
152 *> \verbatim
153 *> LDV2T is INTEGER
154 *> The leading dimension of the array V2T. LDV2T >=
155 *> max(1,M-Q).
156 *> \endverbatim
157 *>
158 *> \param[out] THETA
159 *> \verbatim
160 *> THETA is DOUBLE PRECISION array, dimension MIN(P,M-P,Q,M-Q)
161 *> The CS values of X; the essentially diagonal matrices C and
162 *> S are constructed from THETA; see subroutine ZUNCSD for
163 *> details.
164 *> \endverbatim
165 *>
166 *> \param[out] IWORK
167 *> \verbatim
168 *> IWORK is INTEGER array, dimension (M)
169 *> \endverbatim
170 *>
171 *> \param[out] WORK
172 *> \verbatim
173 *> WORK is COMPLEX*16 array, dimension (LWORK)
174 *> \endverbatim
175 *>
176 *> \param[in] LWORK
177 *> \verbatim
178 *> LWORK is INTEGER
179 *> The dimension of the array WORK
180 *> \endverbatim
181 *>
182 *> \param[out] RWORK
183 *> \verbatim
184 *> RWORK is DOUBLE PRECISION array
185 *> \endverbatim
186 *>
187 *> \param[out] RESULT
188 *> \verbatim
189 *> RESULT is DOUBLE PRECISION array, dimension (15)
190 *> The test ratios:
191 *> First, the 2-by-2 CSD:
192 *> RESULT(1) = norm( U1'*X11*V1 - D11 ) / ( MAX(1,P,Q)*EPS2 )
193 *> RESULT(2) = norm( U1'*X12*V2 - D12 ) / ( MAX(1,P,M-Q)*EPS2 )
194 *> RESULT(3) = norm( U2'*X21*V1 - D21 ) / ( MAX(1,M-P,Q)*EPS2 )
195 *> RESULT(4) = norm( U2'*X22*V2 - D22 ) / ( MAX(1,M-P,M-Q)*EPS2 )
196 *> RESULT(5) = norm( I - U1'*U1 ) / ( MAX(1,P)*ULP )
197 *> RESULT(6) = norm( I - U2'*U2 ) / ( MAX(1,M-P)*ULP )
198 *> RESULT(7) = norm( I - V1T'*V1T ) / ( MAX(1,Q)*ULP )
199 *> RESULT(8) = norm( I - V2T'*V2T ) / ( MAX(1,M-Q)*ULP )
200 *> RESULT(9) = 0 if THETA is in increasing order and
201 *> all angles are in [0,pi/2];
202 *> = ULPINV otherwise.
203 *> Then, the 2-by-1 CSD:
204 *> RESULT(10) = norm( U1'*X11*V1 - D11 ) / ( MAX(1,P,Q)*EPS2 )
205 *> RESULT(11) = norm( U2'*X21*V1 - D21 ) / ( MAX(1,M-P,Q)*EPS2 )
206 *> RESULT(12) = norm( I - U1'*U1 ) / ( MAX(1,P)*ULP )
207 *> RESULT(13) = norm( I - U2'*U2 ) / ( MAX(1,M-P)*ULP )
208 *> RESULT(14) = norm( I - V1T'*V1T ) / ( MAX(1,Q)*ULP )
209 *> RESULT(15) = 0 if THETA is in increasing order and
210 *> all angles are in [0,pi/2];
211 *> = ULPINV otherwise.
212 *> ( EPS2 = MAX( norm( I - X'*X ) / M, ULP ). )
213 *> \endverbatim
214 *
215 * Authors:
216 * ========
217 *
218 *> \author Univ. of Tennessee
219 *> \author Univ. of California Berkeley
220 *> \author Univ. of Colorado Denver
221 *> \author NAG Ltd.
222 *
223 *> \ingroup complex16_eig
224 *
225 * =====================================================================
226  SUBROUTINE zcsdts( M, P, Q, X, XF, LDX, U1, LDU1, U2, LDU2, V1T,
227  $ LDV1T, V2T, LDV2T, THETA, IWORK, WORK, LWORK,
228  $ RWORK, RESULT )
229 *
230 * -- LAPACK test routine --
231 * -- LAPACK is a software package provided by Univ. of Tennessee, --
232 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
233 *
234 * .. Scalar Arguments ..
235  INTEGER LDX, LDU1, LDU2, LDV1T, LDV2T, LWORK, M, P, Q
236 * ..
237 * .. Array Arguments ..
238  INTEGER IWORK( * )
239  DOUBLE PRECISION RESULT( 15 ), RWORK( * ), THETA( * )
240  COMPLEX*16 U1( LDU1, * ), U2( LDU2, * ), V1T( LDV1T, * ),
241  $ v2t( ldv2t, * ), work( lwork ), x( ldx, * ),
242  $ xf( ldx, * )
243 * ..
244 *
245 * =====================================================================
246 *
247 * .. Parameters ..
248  DOUBLE PRECISION REALONE, REALZERO
249  PARAMETER ( REALONE = 1.0d0, realzero = 0.0d0 )
250  COMPLEX*16 ZERO, ONE
251  parameter( zero = (0.0d0,0.0d0), one = (1.0d0,0.0d0) )
252  DOUBLE PRECISION PIOVER2
253  parameter( piover2 = 1.57079632679489661923132169163975144210d0 )
254 * ..
255 * .. Local Scalars ..
256  INTEGER I, INFO, R
257  DOUBLE PRECISION EPS2, RESID, ULP, ULPINV
258 * ..
259 * .. External Functions ..
260  DOUBLE PRECISION DLAMCH, ZLANGE, ZLANHE
261  EXTERNAL DLAMCH, ZLANGE, ZLANHE
262 * ..
263 * .. External Subroutines ..
264  EXTERNAL zgemm, zherk, zlacpy, zlaset, zuncsd,
265  $ zuncsd2by1
266 * ..
267 * .. Intrinsic Functions ..
268  INTRINSIC cos, dble, dcmplx, max, min, real, sin
269 * ..
270 * .. Executable Statements ..
271 *
272  ulp = dlamch( 'Precision' )
273  ulpinv = realone / ulp
274 *
275 * The first half of the routine checks the 2-by-2 CSD
276 *
277  CALL zlaset( 'Full', m, m, zero, one, work, ldx )
278  CALL zherk( 'Upper', 'Conjugate transpose', m, m, -realone,
279  $ x, ldx, realone, work, ldx )
280  IF (m.GT.0) THEN
281  eps2 = max( ulp,
282  $ zlange( '1', m, m, work, ldx, rwork ) / dble( m ) )
283  ELSE
284  eps2 = ulp
285  END IF
286  r = min( p, m-p, q, m-q )
287 *
288 * Copy the matrix X to the array XF.
289 *
290  CALL zlacpy( 'Full', m, m, x, ldx, xf, ldx )
291 *
292 * Compute the CSD
293 *
294  CALL zuncsd( 'Y', 'Y', 'Y', 'Y', 'N', 'D', m, p, q, xf(1,1), ldx,
295  $ xf(1,q+1), ldx, xf(p+1,1), ldx, xf(p+1,q+1), ldx,
296  $ theta, u1, ldu1, u2, ldu2, v1t, ldv1t, v2t, ldv2t,
297  $ work, lwork, rwork, 17*(r+2), iwork, info )
298 *
299 * Compute XF := diag(U1,U2)'*X*diag(V1,V2) - [D11 D12; D21 D22]
300 *
301  CALL zlacpy( 'Full', m, m, x, ldx, xf, ldx )
302 *
303  CALL zgemm( 'No transpose', 'Conjugate transpose', p, q, q, one,
304  $ xf, ldx, v1t, ldv1t, zero, work, ldx )
305 *
306  CALL zgemm( 'Conjugate transpose', 'No transpose', p, q, p, one,
307  $ u1, ldu1, work, ldx, zero, xf, ldx )
308 *
309  DO i = 1, min(p,q)-r
310  xf(i,i) = xf(i,i) - one
311  END DO
312  DO i = 1, r
313  xf(min(p,q)-r+i,min(p,q)-r+i) =
314  $ xf(min(p,q)-r+i,min(p,q)-r+i) - dcmplx( cos(theta(i)),
315  $ 0.0d0 )
316  END DO
317 *
318  CALL zgemm( 'No transpose', 'Conjugate transpose', p, m-q, m-q,
319  $ one, xf(1,q+1), ldx, v2t, ldv2t, zero, work, ldx )
320 *
321  CALL zgemm( 'Conjugate transpose', 'No transpose', p, m-q, p,
322  $ one, u1, ldu1, work, ldx, zero, xf(1,q+1), ldx )
323 *
324  DO i = 1, min(p,m-q)-r
325  xf(p-i+1,m-i+1) = xf(p-i+1,m-i+1) + one
326  END DO
327  DO i = 1, r
328  xf(p-(min(p,m-q)-r)+1-i,m-(min(p,m-q)-r)+1-i) =
329  $ xf(p-(min(p,m-q)-r)+1-i,m-(min(p,m-q)-r)+1-i) +
330  $ dcmplx( sin(theta(r-i+1)), 0.0d0 )
331  END DO
332 *
333  CALL zgemm( 'No transpose', 'Conjugate transpose', m-p, q, q, one,
334  $ xf(p+1,1), ldx, v1t, ldv1t, zero, work, ldx )
335 *
336  CALL zgemm( 'Conjugate transpose', 'No transpose', m-p, q, m-p,
337  $ one, u2, ldu2, work, ldx, zero, xf(p+1,1), ldx )
338 *
339  DO i = 1, min(m-p,q)-r
340  xf(m-i+1,q-i+1) = xf(m-i+1,q-i+1) - one
341  END DO
342  DO i = 1, r
343  xf(m-(min(m-p,q)-r)+1-i,q-(min(m-p,q)-r)+1-i) =
344  $ xf(m-(min(m-p,q)-r)+1-i,q-(min(m-p,q)-r)+1-i) -
345  $ dcmplx( sin(theta(r-i+1)), 0.0d0 )
346  END DO
347 *
348  CALL zgemm( 'No transpose', 'Conjugate transpose', m-p, m-q, m-q,
349  $ one, xf(p+1,q+1), ldx, v2t, ldv2t, zero, work, ldx )
350 *
351  CALL zgemm( 'Conjugate transpose', 'No transpose', m-p, m-q, m-p,
352  $ one, u2, ldu2, work, ldx, zero, xf(p+1,q+1), ldx )
353 *
354  DO i = 1, min(m-p,m-q)-r
355  xf(p+i,q+i) = xf(p+i,q+i) - one
356  END DO
357  DO i = 1, r
358  xf(p+(min(m-p,m-q)-r)+i,q+(min(m-p,m-q)-r)+i) =
359  $ xf(p+(min(m-p,m-q)-r)+i,q+(min(m-p,m-q)-r)+i) -
360  $ dcmplx( cos(theta(i)), 0.0d0 )
361  END DO
362 *
363 * Compute norm( U1'*X11*V1 - D11 ) / ( MAX(1,P,Q)*EPS2 ) .
364 *
365  resid = zlange( '1', p, q, xf, ldx, rwork )
366  result( 1 ) = ( resid / real(max(1,p,q)) ) / eps2
367 *
368 * Compute norm( U1'*X12*V2 - D12 ) / ( MAX(1,P,M-Q)*EPS2 ) .
369 *
370  resid = zlange( '1', p, m-q, xf(1,q+1), ldx, rwork )
371  result( 2 ) = ( resid / real(max(1,p,m-q)) ) / eps2
372 *
373 * Compute norm( U2'*X21*V1 - D21 ) / ( MAX(1,M-P,Q)*EPS2 ) .
374 *
375  resid = zlange( '1', m-p, q, xf(p+1,1), ldx, rwork )
376  result( 3 ) = ( resid / real(max(1,m-p,q)) ) / eps2
377 *
378 * Compute norm( U2'*X22*V2 - D22 ) / ( MAX(1,M-P,M-Q)*EPS2 ) .
379 *
380  resid = zlange( '1', m-p, m-q, xf(p+1,q+1), ldx, rwork )
381  result( 4 ) = ( resid / real(max(1,m-p,m-q)) ) / eps2
382 *
383 * Compute I - U1'*U1
384 *
385  CALL zlaset( 'Full', p, p, zero, one, work, ldu1 )
386  CALL zherk( 'Upper', 'Conjugate transpose', p, p, -realone,
387  $ u1, ldu1, realone, work, ldu1 )
388 *
389 * Compute norm( I - U'*U ) / ( MAX(1,P) * ULP ) .
390 *
391  resid = zlanhe( '1', 'Upper', p, work, ldu1, rwork )
392  result( 5 ) = ( resid / real(max(1,p)) ) / ulp
393 *
394 * Compute I - U2'*U2
395 *
396  CALL zlaset( 'Full', m-p, m-p, zero, one, work, ldu2 )
397  CALL zherk( 'Upper', 'Conjugate transpose', m-p, m-p, -realone,
398  $ u2, ldu2, realone, work, ldu2 )
399 *
400 * Compute norm( I - U2'*U2 ) / ( MAX(1,M-P) * ULP ) .
401 *
402  resid = zlanhe( '1', 'Upper', m-p, work, ldu2, rwork )
403  result( 6 ) = ( resid / real(max(1,m-p)) ) / ulp
404 *
405 * Compute I - V1T*V1T'
406 *
407  CALL zlaset( 'Full', q, q, zero, one, work, ldv1t )
408  CALL zherk( 'Upper', 'No transpose', q, q, -realone,
409  $ v1t, ldv1t, realone, work, ldv1t )
410 *
411 * Compute norm( I - V1T*V1T' ) / ( MAX(1,Q) * ULP ) .
412 *
413  resid = zlanhe( '1', 'Upper', q, work, ldv1t, rwork )
414  result( 7 ) = ( resid / real(max(1,q)) ) / ulp
415 *
416 * Compute I - V2T*V2T'
417 *
418  CALL zlaset( 'Full', m-q, m-q, zero, one, work, ldv2t )
419  CALL zherk( 'Upper', 'No transpose', m-q, m-q, -realone,
420  $ v2t, ldv2t, realone, work, ldv2t )
421 *
422 * Compute norm( I - V2T*V2T' ) / ( MAX(1,M-Q) * ULP ) .
423 *
424  resid = zlanhe( '1', 'Upper', m-q, work, ldv2t, rwork )
425  result( 8 ) = ( resid / real(max(1,m-q)) ) / ulp
426 *
427 * Check sorting
428 *
429  result( 9 ) = realzero
430  DO i = 1, r
431  IF( theta(i).LT.realzero .OR. theta(i).GT.piover2 ) THEN
432  result( 9 ) = ulpinv
433  END IF
434  IF( i.GT.1) THEN
435  IF ( theta(i).LT.theta(i-1) ) THEN
436  result( 9 ) = ulpinv
437  END IF
438  END IF
439  END DO
440 *
441 * The second half of the routine checks the 2-by-1 CSD
442 *
443  CALL zlaset( 'Full', q, q, zero, one, work, ldx )
444  CALL zherk( 'Upper', 'Conjugate transpose', q, m, -realone,
445  $ x, ldx, realone, work, ldx )
446  IF (m.GT.0) THEN
447  eps2 = max( ulp,
448  $ zlange( '1', q, q, work, ldx, rwork ) / dble( m ) )
449  ELSE
450  eps2 = ulp
451  END IF
452  r = min( p, m-p, q, m-q )
453 *
454 * Copy the matrix X to the array XF.
455 *
456  CALL zlacpy( 'Full', m, m, x, ldx, xf, ldx )
457 *
458 * Compute the CSD
459 *
460  CALL zuncsd2by1( 'Y', 'Y', 'Y', m, p, q, xf(1,1), ldx, xf(p+1,1),
461  $ ldx, theta, u1, ldu1, u2, ldu2, v1t, ldv1t, work,
462  $ lwork, rwork, 17*(r+2), iwork, info )
463 *
464 * Compute [X11;X21] := diag(U1,U2)'*[X11;X21]*V1 - [D11;D21]
465 *
466  CALL zgemm( 'No transpose', 'Conjugate transpose', p, q, q, one,
467  $ x, ldx, v1t, ldv1t, zero, work, ldx )
468 *
469  CALL zgemm( 'Conjugate transpose', 'No transpose', p, q, p, one,
470  $ u1, ldu1, work, ldx, zero, x, ldx )
471 *
472  DO i = 1, min(p,q)-r
473  x(i,i) = x(i,i) - one
474  END DO
475  DO i = 1, r
476  x(min(p,q)-r+i,min(p,q)-r+i) =
477  $ x(min(p,q)-r+i,min(p,q)-r+i) - dcmplx( cos(theta(i)),
478  $ 0.0d0 )
479  END DO
480 *
481  CALL zgemm( 'No transpose', 'Conjugate transpose', m-p, q, q, one,
482  $ x(p+1,1), ldx, v1t, ldv1t, zero, work, ldx )
483 *
484  CALL zgemm( 'Conjugate transpose', 'No transpose', m-p, q, m-p,
485  $ one, u2, ldu2, work, ldx, zero, x(p+1,1), ldx )
486 *
487  DO i = 1, min(m-p,q)-r
488  x(m-i+1,q-i+1) = x(m-i+1,q-i+1) - one
489  END DO
490  DO i = 1, r
491  x(m-(min(m-p,q)-r)+1-i,q-(min(m-p,q)-r)+1-i) =
492  $ x(m-(min(m-p,q)-r)+1-i,q-(min(m-p,q)-r)+1-i) -
493  $ dcmplx( sin(theta(r-i+1)), 0.0d0 )
494  END DO
495 *
496 * Compute norm( U1'*X11*V1 - D11 ) / ( MAX(1,P,Q)*EPS2 ) .
497 *
498  resid = zlange( '1', p, q, x, ldx, rwork )
499  result( 10 ) = ( resid / real(max(1,p,q)) ) / eps2
500 *
501 * Compute norm( U2'*X21*V1 - D21 ) / ( MAX(1,M-P,Q)*EPS2 ) .
502 *
503  resid = zlange( '1', m-p, q, x(p+1,1), ldx, rwork )
504  result( 11 ) = ( resid / real(max(1,m-p,q)) ) / eps2
505 *
506 * Compute I - U1'*U1
507 *
508  CALL zlaset( 'Full', p, p, zero, one, work, ldu1 )
509  CALL zherk( 'Upper', 'Conjugate transpose', p, p, -realone,
510  $ u1, ldu1, realone, work, ldu1 )
511 *
512 * Compute norm( I - U'*U ) / ( MAX(1,P) * ULP ) .
513 *
514  resid = zlanhe( '1', 'Upper', p, work, ldu1, rwork )
515  result( 12 ) = ( resid / real(max(1,p)) ) / ulp
516 *
517 * Compute I - U2'*U2
518 *
519  CALL zlaset( 'Full', m-p, m-p, zero, one, work, ldu2 )
520  CALL zherk( 'Upper', 'Conjugate transpose', m-p, m-p, -realone,
521  $ u2, ldu2, realone, work, ldu2 )
522 *
523 * Compute norm( I - U2'*U2 ) / ( MAX(1,M-P) * ULP ) .
524 *
525  resid = zlanhe( '1', 'Upper', m-p, work, ldu2, rwork )
526  result( 13 ) = ( resid / real(max(1,m-p)) ) / ulp
527 *
528 * Compute I - V1T*V1T'
529 *
530  CALL zlaset( 'Full', q, q, zero, one, work, ldv1t )
531  CALL zherk( 'Upper', 'No transpose', q, q, -realone,
532  $ v1t, ldv1t, realone, work, ldv1t )
533 *
534 * Compute norm( I - V1T*V1T' ) / ( MAX(1,Q) * ULP ) .
535 *
536  resid = zlanhe( '1', 'Upper', q, work, ldv1t, rwork )
537  result( 14 ) = ( resid / real(max(1,q)) ) / ulp
538 *
539 * Check sorting
540 *
541  result( 15 ) = realzero
542  DO i = 1, r
543  IF( theta(i).LT.realzero .OR. theta(i).GT.piover2 ) THEN
544  result( 15 ) = ulpinv
545  END IF
546  IF( i.GT.1) THEN
547  IF ( theta(i).LT.theta(i-1) ) THEN
548  result( 15 ) = ulpinv
549  END IF
550  END IF
551  END DO
552 *
553  RETURN
554 *
555 * End of ZCSDTS
556 *
557  END
558 
subroutine zgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
ZGEMM
Definition: zgemm.f:187
subroutine zherk(UPLO, TRANS, N, K, ALPHA, A, LDA, BETA, C, LDC)
ZHERK
Definition: zherk.f:173
subroutine zcsdts(M, P, Q, X, XF, LDX, U1, LDU1, U2, LDU2, V1T, LDV1T, V2T, LDV2T, THETA, IWORK, WORK, LWORK, RWORK, RESULT)
ZCSDTS
Definition: zcsdts.f:229
subroutine zlacpy(UPLO, M, N, A, LDA, B, LDB)
ZLACPY copies all or part of one two-dimensional array to another.
Definition: zlacpy.f:103
subroutine zlaset(UPLO, M, N, ALPHA, BETA, A, LDA)
ZLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition: zlaset.f:106
recursive subroutine zuncsd(JOBU1, JOBU2, JOBV1T, JOBV2T, TRANS, SIGNS, M, P, Q, X11, LDX11, X12, LDX12, X21, LDX21, X22, LDX22, THETA, U1, LDU1, U2, LDU2, V1T, LDV1T, V2T, LDV2T, WORK, LWORK, RWORK, LRWORK, IWORK, INFO)
ZUNCSD
Definition: zuncsd.f:320
subroutine zuncsd2by1(JOBU1, JOBU2, JOBV1T, M, P, Q, X11, LDX11, X21, LDX21, THETA, U1, LDU1, U2, LDU2, V1T, LDV1T, WORK, LWORK, RWORK, LRWORK, IWORK, INFO)
ZUNCSD2BY1
Definition: zuncsd2by1.f:256