LAPACK  3.6.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 *> \date November 2015
224 *
225 *> \ingroup complex16_eig
226 *
227 * =====================================================================
228  SUBROUTINE zcsdts( M, P, Q, X, XF, LDX, U1, LDU1, U2, LDU2, V1T,
229  $ ldv1t, v2t, ldv2t, theta, iwork, work, lwork,
230  $ rwork, result )
231 *
232 * -- LAPACK test routine (version 3.6.0) --
233 * -- LAPACK is a software package provided by Univ. of Tennessee, --
234 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
235 * November 2015
236 *
237 * .. Scalar Arguments ..
238  INTEGER LDX, LDU1, LDU2, LDV1T, LDV2T, LWORK, M, P, Q
239 * ..
240 * .. Array Arguments ..
241  INTEGER IWORK( * )
242  DOUBLE PRECISION RESULT( 15 ), RWORK( * ), THETA( * )
243  COMPLEX*16 U1( ldu1, * ), U2( ldu2, * ), V1T( ldv1t, * ),
244  $ v2t( ldv2t, * ), work( lwork ), x( ldx, * ),
245  $ xf( ldx, * )
246 * ..
247 *
248 * =====================================================================
249 *
250 * .. Parameters ..
251  DOUBLE PRECISION PIOVER2, REALONE, REALZERO
252  parameter ( piover2 = 1.57079632679489662d0,
253  $ realone = 1.0d0, realzero = 0.0d0 )
254  COMPLEX*16 ZERO, ONE
255  parameter ( zero = (0.0d0,0.0d0), one = (1.0d0,0.0d0) )
256 * ..
257 * .. Local Scalars ..
258  INTEGER I, INFO, R
259  DOUBLE PRECISION EPS2, RESID, ULP, ULPINV
260 * ..
261 * .. External Functions ..
262  DOUBLE PRECISION DLAMCH, ZLANGE, ZLANHE
263  EXTERNAL dlamch, zlange, zlanhe
264 * ..
265 * .. External Subroutines ..
266  EXTERNAL zgemm, zherk, zlacpy, zlaset, zuncsd,
267  $ zuncsd2by1
268 * ..
269 * .. Intrinsic Functions ..
270  INTRINSIC cos, dble, dcmplx, max, min, REAL, SIN
271 * ..
272 * .. Executable Statements ..
273 *
274  ulp = dlamch( 'Precision' )
275  ulpinv = realone / ulp
276 *
277 * The first half of the routine checks the 2-by-2 CSD
278 *
279  CALL zlaset( 'Full', m, m, zero, one, work, ldx )
280  CALL zherk( 'Upper', 'Conjugate transpose', m, m, -realone,
281  $ x, ldx, realone, work, ldx )
282  IF (m.GT.0) THEN
283  eps2 = max( ulp,
284  $ zlange( '1', m, m, work, ldx, rwork ) / dble( m ) )
285  ELSE
286  eps2 = ulp
287  END IF
288  r = min( p, m-p, q, m-q )
289 *
290 * Copy the matrix X to the array XF.
291 *
292  CALL zlacpy( 'Full', m, m, x, ldx, xf, ldx )
293 *
294 * Compute the CSD
295 *
296  CALL zuncsd( 'Y', 'Y', 'Y', 'Y', 'N', 'D', m, p, q, xf(1,1), ldx,
297  $ xf(1,q+1), ldx, xf(p+1,1), ldx, xf(p+1,q+1), ldx,
298  $ theta, u1, ldu1, u2, ldu2, v1t, ldv1t, v2t, ldv2t,
299  $ work, lwork, rwork, 17*(r+2), iwork, info )
300 *
301 * Compute XF := diag(U1,U2)'*X*diag(V1,V2) - [D11 D12; D21 D22]
302 *
303  CALL zlacpy( 'Full', m, m, x, ldx, xf, ldx )
304 *
305  CALL zgemm( 'No transpose', 'Conjugate transpose', p, q, q, one,
306  $ xf, ldx, v1t, ldv1t, zero, work, ldx )
307 *
308  CALL zgemm( 'Conjugate transpose', 'No transpose', p, q, p, one,
309  $ u1, ldu1, work, ldx, zero, xf, ldx )
310 *
311  DO i = 1, min(p,q)-r
312  xf(i,i) = xf(i,i) - one
313  END DO
314  DO i = 1, r
315  xf(min(p,q)-r+i,min(p,q)-r+i) =
316  $ xf(min(p,q)-r+i,min(p,q)-r+i) - dcmplx( cos(theta(i)),
317  $ 0.0d0 )
318  END DO
319 *
320  CALL zgemm( 'No transpose', 'Conjugate transpose', p, m-q, m-q,
321  $ one, xf(1,q+1), ldx, v2t, ldv2t, zero, work, ldx )
322 *
323  CALL zgemm( 'Conjugate transpose', 'No transpose', p, m-q, p,
324  $ one, u1, ldu1, work, ldx, zero, xf(1,q+1), ldx )
325 *
326  DO i = 1, min(p,m-q)-r
327  xf(p-i+1,m-i+1) = xf(p-i+1,m-i+1) + one
328  END DO
329  DO i = 1, r
330  xf(p-(min(p,m-q)-r)+1-i,m-(min(p,m-q)-r)+1-i) =
331  $ xf(p-(min(p,m-q)-r)+1-i,m-(min(p,m-q)-r)+1-i) +
332  $ dcmplx( sin(theta(r-i+1)), 0.0d0 )
333  END DO
334 *
335  CALL zgemm( 'No transpose', 'Conjugate transpose', m-p, q, q, one,
336  $ xf(p+1,1), ldx, v1t, ldv1t, zero, work, ldx )
337 *
338  CALL zgemm( 'Conjugate transpose', 'No transpose', m-p, q, m-p,
339  $ one, u2, ldu2, work, ldx, zero, xf(p+1,1), ldx )
340 *
341  DO i = 1, min(m-p,q)-r
342  xf(m-i+1,q-i+1) = xf(m-i+1,q-i+1) - one
343  END DO
344  DO i = 1, r
345  xf(m-(min(m-p,q)-r)+1-i,q-(min(m-p,q)-r)+1-i) =
346  $ xf(m-(min(m-p,q)-r)+1-i,q-(min(m-p,q)-r)+1-i) -
347  $ dcmplx( sin(theta(r-i+1)), 0.0d0 )
348  END DO
349 *
350  CALL zgemm( 'No transpose', 'Conjugate transpose', m-p, m-q, m-q,
351  $ one, xf(p+1,q+1), ldx, v2t, ldv2t, zero, work, ldx )
352 *
353  CALL zgemm( 'Conjugate transpose', 'No transpose', m-p, m-q, m-p,
354  $ one, u2, ldu2, work, ldx, zero, xf(p+1,q+1), ldx )
355 *
356  DO i = 1, min(m-p,m-q)-r
357  xf(p+i,q+i) = xf(p+i,q+i) - one
358  END DO
359  DO i = 1, r
360  xf(p+(min(m-p,m-q)-r)+i,q+(min(m-p,m-q)-r)+i) =
361  $ xf(p+(min(m-p,m-q)-r)+i,q+(min(m-p,m-q)-r)+i) -
362  $ dcmplx( cos(theta(i)), 0.0d0 )
363  END DO
364 *
365 * Compute norm( U1'*X11*V1 - D11 ) / ( MAX(1,P,Q)*EPS2 ) .
366 *
367  resid = zlange( '1', p, q, xf, ldx, rwork )
368  result( 1 ) = ( resid / REAL(MAX(1,P,Q)) ) / eps2
369 *
370 * Compute norm( U1'*X12*V2 - D12 ) / ( MAX(1,P,M-Q)*EPS2 ) .
371 *
372  resid = zlange( '1', p, m-q, xf(1,q+1), ldx, rwork )
373  result( 2 ) = ( resid / REAL(MAX(1,P,M-Q)) ) / eps2
374 *
375 * Compute norm( U2'*X21*V1 - D21 ) / ( MAX(1,M-P,Q)*EPS2 ) .
376 *
377  resid = zlange( '1', m-p, q, xf(p+1,1), ldx, rwork )
378  result( 3 ) = ( resid / REAL(MAX(1,M-P,Q)) ) / eps2
379 *
380 * Compute norm( U2'*X22*V2 - D22 ) / ( MAX(1,M-P,M-Q)*EPS2 ) .
381 *
382  resid = zlange( '1', m-p, m-q, xf(p+1,q+1), ldx, rwork )
383  result( 4 ) = ( resid / REAL(MAX(1,M-P,M-Q)) ) / eps2
384 *
385 * Compute I - U1'*U1
386 *
387  CALL zlaset( 'Full', p, p, zero, one, work, ldu1 )
388  CALL zherk( 'Upper', 'Conjugate transpose', p, p, -realone,
389  $ u1, ldu1, realone, work, ldu1 )
390 *
391 * Compute norm( I - U'*U ) / ( MAX(1,P) * ULP ) .
392 *
393  resid = zlanhe( '1', 'Upper', p, work, ldu1, rwork )
394  result( 5 ) = ( resid / REAL(MAX(1,P)) ) / ulp
395 *
396 * Compute I - U2'*U2
397 *
398  CALL zlaset( 'Full', m-p, m-p, zero, one, work, ldu2 )
399  CALL zherk( 'Upper', 'Conjugate transpose', m-p, m-p, -realone,
400  $ u2, ldu2, realone, work, ldu2 )
401 *
402 * Compute norm( I - U2'*U2 ) / ( MAX(1,M-P) * ULP ) .
403 *
404  resid = zlanhe( '1', 'Upper', m-p, work, ldu2, rwork )
405  result( 6 ) = ( resid / REAL(MAX(1,M-P)) ) / ulp
406 *
407 * Compute I - V1T*V1T'
408 *
409  CALL zlaset( 'Full', q, q, zero, one, work, ldv1t )
410  CALL zherk( 'Upper', 'No transpose', q, q, -realone,
411  $ v1t, ldv1t, realone, work, ldv1t )
412 *
413 * Compute norm( I - V1T*V1T' ) / ( MAX(1,Q) * ULP ) .
414 *
415  resid = zlanhe( '1', 'Upper', q, work, ldv1t, rwork )
416  result( 7 ) = ( resid / REAL(MAX(1,Q)) ) / ulp
417 *
418 * Compute I - V2T*V2T'
419 *
420  CALL zlaset( 'Full', m-q, m-q, zero, one, work, ldv2t )
421  CALL zherk( 'Upper', 'No transpose', m-q, m-q, -realone,
422  $ v2t, ldv2t, realone, work, ldv2t )
423 *
424 * Compute norm( I - V2T*V2T' ) / ( MAX(1,M-Q) * ULP ) .
425 *
426  resid = zlanhe( '1', 'Upper', m-q, work, ldv2t, rwork )
427  result( 8 ) = ( resid / REAL(MAX(1,M-Q)) ) / ulp
428 *
429 * Check sorting
430 *
431  result( 9 ) = realzero
432  DO i = 1, r
433  IF( theta(i).LT.realzero .OR. theta(i).GT.piover2 ) THEN
434  result( 9 ) = ulpinv
435  END IF
436  IF( i.GT.1) THEN
437  IF ( theta(i).LT.theta(i-1) ) THEN
438  result( 9 ) = ulpinv
439  END IF
440  END IF
441  END DO
442 *
443 * The second half of the routine checks the 2-by-1 CSD
444 *
445  CALL zlaset( 'Full', q, q, zero, one, work, ldx )
446  CALL zherk( 'Upper', 'Conjugate transpose', q, m, -realone,
447  $ x, ldx, realone, work, ldx )
448  IF (m.GT.0) THEN
449  eps2 = max( ulp,
450  $ zlange( '1', q, q, work, ldx, rwork ) / dble( m ) )
451  ELSE
452  eps2 = ulp
453  END IF
454  r = min( p, m-p, q, m-q )
455 *
456 * Copy the matrix X to the array XF.
457 *
458  CALL zlacpy( 'Full', m, m, x, ldx, xf, ldx )
459 *
460 * Compute the CSD
461 *
462  CALL zuncsd2by1( 'Y', 'Y', 'Y', m, p, q, xf(1,1), ldx, xf(p+1,1),
463  $ ldx, theta, u1, ldu1, u2, ldu2, v1t, ldv1t, work,
464  $ lwork, rwork, 17*(r+2), iwork, info )
465 *
466 * Compute [X11;X21] := diag(U1,U2)'*[X11;X21]*V1 - [D11;D21]
467 *
468  CALL zgemm( 'No transpose', 'Conjugate transpose', p, q, q, one,
469  $ x, ldx, v1t, ldv1t, zero, work, ldx )
470 *
471  CALL zgemm( 'Conjugate transpose', 'No transpose', p, q, p, one,
472  $ u1, ldu1, work, ldx, zero, x, ldx )
473 *
474  DO i = 1, min(p,q)-r
475  x(i,i) = x(i,i) - one
476  END DO
477  DO i = 1, r
478  x(min(p,q)-r+i,min(p,q)-r+i) =
479  $ x(min(p,q)-r+i,min(p,q)-r+i) - dcmplx( cos(theta(i)),
480  $ 0.0d0 )
481  END DO
482 *
483  CALL zgemm( 'No transpose', 'Conjugate transpose', m-p, q, q, one,
484  $ x(p+1,1), ldx, v1t, ldv1t, zero, work, ldx )
485 *
486  CALL zgemm( 'Conjugate transpose', 'No transpose', m-p, q, m-p,
487  $ one, u2, ldu2, work, ldx, zero, x(p+1,1), ldx )
488 *
489  DO i = 1, min(m-p,q)-r
490  x(m-i+1,q-i+1) = x(m-i+1,q-i+1) - one
491  END DO
492  DO i = 1, r
493  x(m-(min(m-p,q)-r)+1-i,q-(min(m-p,q)-r)+1-i) =
494  $ x(m-(min(m-p,q)-r)+1-i,q-(min(m-p,q)-r)+1-i) -
495  $ dcmplx( sin(theta(r-i+1)), 0.0d0 )
496  END DO
497 *
498 * Compute norm( U1'*X11*V1 - D11 ) / ( MAX(1,P,Q)*EPS2 ) .
499 *
500  resid = zlange( '1', p, q, x, ldx, rwork )
501  result( 10 ) = ( resid / REAL(MAX(1,P,Q)) ) / eps2
502 *
503 * Compute norm( U2'*X21*V1 - D21 ) / ( MAX(1,M-P,Q)*EPS2 ) .
504 *
505  resid = zlange( '1', m-p, q, x(p+1,1), ldx, rwork )
506  result( 11 ) = ( resid / REAL(MAX(1,M-P,Q)) ) / eps2
507 *
508 * Compute I - U1'*U1
509 *
510  CALL zlaset( 'Full', p, p, zero, one, work, ldu1 )
511  CALL zherk( 'Upper', 'Conjugate transpose', p, p, -realone,
512  $ u1, ldu1, realone, work, ldu1 )
513 *
514 * Compute norm( I - U'*U ) / ( MAX(1,P) * ULP ) .
515 *
516  resid = zlanhe( '1', 'Upper', p, work, ldu1, rwork )
517  result( 12 ) = ( resid / REAL(MAX(1,P)) ) / ulp
518 *
519 * Compute I - U2'*U2
520 *
521  CALL zlaset( 'Full', m-p, m-p, zero, one, work, ldu2 )
522  CALL zherk( 'Upper', 'Conjugate transpose', m-p, m-p, -realone,
523  $ u2, ldu2, realone, work, ldu2 )
524 *
525 * Compute norm( I - U2'*U2 ) / ( MAX(1,M-P) * ULP ) .
526 *
527  resid = zlanhe( '1', 'Upper', m-p, work, ldu2, rwork )
528  result( 13 ) = ( resid / REAL(MAX(1,M-P)) ) / ulp
529 *
530 * Compute I - V1T*V1T'
531 *
532  CALL zlaset( 'Full', q, q, zero, one, work, ldv1t )
533  CALL zherk( 'Upper', 'No transpose', q, q, -realone,
534  $ v1t, ldv1t, realone, work, ldv1t )
535 *
536 * Compute norm( I - V1T*V1T' ) / ( MAX(1,Q) * ULP ) .
537 *
538  resid = zlanhe( '1', 'Upper', q, work, ldv1t, rwork )
539  result( 14 ) = ( resid / REAL(MAX(1,Q)) ) / ulp
540 *
541 * Check sorting
542 *
543  result( 15 ) = realzero
544  DO i = 1, r
545  IF( theta(i).LT.realzero .OR. theta(i).GT.piover2 ) THEN
546  result( 15 ) = ulpinv
547  END IF
548  IF( i.GT.1) THEN
549  IF ( theta(i).LT.theta(i-1) ) THEN
550  result( 15 ) = ulpinv
551  END IF
552  END IF
553  END DO
554 *
555  RETURN
556 *
557 * End of ZCSDTS
558 *
559  END
560 
zcsdts
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:231
zlacpy
subroutine zlacpy(UPLO, M, N, A, LDA, B, LDB)
ZLACPY copies all or part of one two-dimensional array to another.
Definition: zlacpy.f:105
zgemm
subroutine zgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
ZGEMM
Definition: zgemm.f:189
zlaset
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:108
zuncsd
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:322
zuncsd2by1
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:255
zherk
subroutine zherk(UPLO, TRANS, N, K, ALPHA, A, LDA, BETA, C, LDC)
ZHERK
Definition: zherk.f:175