LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
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( 9 ), 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 *> \endverbatim
51 *
52 * Arguments:
53 * ==========
54 *
55 *> \param[in] M
56 *> \verbatim
57 *> M is INTEGER
58 *> The number of rows of the matrix X. M >= 0.
59 *> \endverbatim
60 *>
61 *> \param[in] P
62 *> \verbatim
63 *> P is INTEGER
64 *> The number of rows of the matrix X11. P >= 0.
65 *> \endverbatim
66 *>
67 *> \param[in] Q
68 *> \verbatim
69 *> Q is INTEGER
70 *> The number of columns of the matrix X11. Q >= 0.
71 *> \endverbatim
72 *>
73 *> \param[in] X
74 *> \verbatim
75 *> X is COMPLEX*16 array, dimension (LDX,M)
76 *> The M-by-M matrix X.
77 *> \endverbatim
78 *>
79 *> \param[out] XF
80 *> \verbatim
81 *> XF is COMPLEX*16 array, dimension (LDX,M)
82 *> Details of the CSD of X, as returned by ZUNCSD;
83 *> see ZUNCSD for further details.
84 *> \endverbatim
85 *>
86 *> \param[in] LDX
87 *> \verbatim
88 *> LDX is INTEGER
89 *> The leading dimension of the arrays X and XF.
90 *> LDX >= max( 1,M ).
91 *> \endverbatim
92 *>
93 *> \param[out] U1
94 *> \verbatim
95 *> U1 is COMPLEX*16 array, dimension(LDU1,P)
96 *> The P-by-P unitary matrix U1.
97 *> \endverbatim
98 *>
99 *> \param[in] LDU1
100 *> \verbatim
101 *> LDU1 is INTEGER
102 *> The leading dimension of the array U1. LDU >= max(1,P).
103 *> \endverbatim
104 *>
105 *> \param[out] U2
106 *> \verbatim
107 *> U2 is COMPLEX*16 array, dimension(LDU2,M-P)
108 *> The (M-P)-by-(M-P) unitary matrix U2.
109 *> \endverbatim
110 *>
111 *> \param[in] LDU2
112 *> \verbatim
113 *> LDU2 is INTEGER
114 *> The leading dimension of the array U2. LDU >= max(1,M-P).
115 *> \endverbatim
116 *>
117 *> \param[out] V1T
118 *> \verbatim
119 *> V1T is COMPLEX*16 array, dimension(LDV1T,Q)
120 *> The Q-by-Q unitary matrix V1T.
121 *> \endverbatim
122 *>
123 *> \param[in] LDV1T
124 *> \verbatim
125 *> LDV1T is INTEGER
126 *> The leading dimension of the array V1T. LDV1T >=
127 *> max(1,Q).
128 *> \endverbatim
129 *>
130 *> \param[out] V2T
131 *> \verbatim
132 *> V2T is COMPLEX*16 array, dimension(LDV2T,M-Q)
133 *> The (M-Q)-by-(M-Q) unitary matrix V2T.
134 *> \endverbatim
135 *>
136 *> \param[in] LDV2T
137 *> \verbatim
138 *> LDV2T is INTEGER
139 *> The leading dimension of the array V2T. LDV2T >=
140 *> max(1,M-Q).
141 *> \endverbatim
142 *>
143 *> \param[out] THETA
144 *> \verbatim
145 *> THETA is DOUBLE PRECISION array, dimension MIN(P,M-P,Q,M-Q)
146 *> The CS values of X; the essentially diagonal matrices C and
147 *> S are constructed from THETA; see subroutine ZUNCSD for
148 *> details.
149 *> \endverbatim
150 *>
151 *> \param[out] IWORK
152 *> \verbatim
153 *> IWORK is INTEGER array, dimension (M)
154 *> \endverbatim
155 *>
156 *> \param[out] WORK
157 *> \verbatim
158 *> WORK is COMPLEX*16 array, dimension (LWORK)
159 *> \endverbatim
160 *>
161 *> \param[in] LWORK
162 *> \verbatim
163 *> LWORK is INTEGER
164 *> The dimension of the array WORK
165 *> \endverbatim
166 *>
167 *> \param[out] RWORK
168 *> \verbatim
169 *> RWORK is DOUBLE PRECISION array
170 *> \endverbatim
171 *>
172 *> \param[out] RESULT
173 *> \verbatim
174 *> RESULT is DOUBLE PRECISION array, dimension (9)
175 *> The test ratios:
176 *> RESULT(1) = norm( U1'*X11*V1 - D11 ) / ( MAX(1,P,Q)*EPS2 )
177 *> RESULT(2) = norm( U1'*X12*V2 - D12 ) / ( MAX(1,P,M-Q)*EPS2 )
178 *> RESULT(3) = norm( U2'*X21*V1 - D21 ) / ( MAX(1,M-P,Q)*EPS2 )
179 *> RESULT(4) = norm( U2'*X22*V2 - D22 ) / ( MAX(1,M-P,M-Q)*EPS2 )
180 *> RESULT(5) = norm( I - U1'*U1 ) / ( MAX(1,P)*ULP )
181 *> RESULT(6) = norm( I - U2'*U2 ) / ( MAX(1,M-P)*ULP )
182 *> RESULT(7) = norm( I - V1T'*V1T ) / ( MAX(1,Q)*ULP )
183 *> RESULT(8) = norm( I - V2T'*V2T ) / ( MAX(1,M-Q)*ULP )
184 *> RESULT(9) = 0 if THETA is in increasing order and
185 *> all angles are in [0,pi/2];
186 *> = ULPINV otherwise.
187 *> ( EPS2 = MAX( norm( I - X'*X ) / M, ULP ). )
188 *> \endverbatim
189 *
190 * Authors:
191 * ========
192 *
193 *> \author Univ. of Tennessee
194 *> \author Univ. of California Berkeley
195 *> \author Univ. of Colorado Denver
196 *> \author NAG Ltd.
197 *
198 *> \date November 2011
199 *
200 *> \ingroup complex16_eig
201 *
202 * =====================================================================
203  SUBROUTINE zcsdts( M, P, Q, X, XF, LDX, U1, LDU1, U2, LDU2, V1T,
204  $ ldv1t, v2t, ldv2t, theta, iwork, work, lwork,
205  $ rwork, result )
206 *
207 * -- LAPACK test routine (version 3.4.0) --
208 * -- LAPACK is a software package provided by Univ. of Tennessee, --
209 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
210 * November 2011
211 *
212 * .. Scalar Arguments ..
213  INTEGER ldx, ldu1, ldu2, ldv1t, ldv2t, lwork, m, p, q
214 * ..
215 * .. Array Arguments ..
216  INTEGER iwork( * )
217  DOUBLE PRECISION result( 9 ), rwork( * ), theta( * )
218  COMPLEX*16 u1( ldu1, * ), u2( ldu2, * ), v1t( ldv1t, * ),
219  $ v2t( ldv2t, * ), work( lwork ), x( ldx, * ),
220  $ xf( ldx, * )
221 * ..
222 *
223 * =====================================================================
224 *
225 * .. Parameters ..
226  DOUBLE PRECISION piover2, realone, realzero
227  parameter( piover2 = 1.57079632679489662d0,
228  $ realone = 1.0d0, realzero = 0.0d0 )
229  COMPLEX*16 zero, one
230  parameter( zero = (0.0d0,0.0d0), one = (1.0d0,0.0d0) )
231 * ..
232 * .. Local Scalars ..
233  INTEGER i, info, r
234  DOUBLE PRECISION eps2, resid, ulp, ulpinv
235 * ..
236 * .. External Functions ..
237  DOUBLE PRECISION dlamch, zlange, zlanhe
238  EXTERNAL dlamch, zlange, zlanhe
239 * ..
240 * .. External Subroutines ..
241  EXTERNAL zgemm, zlacpy, zlaset, zuncsd, zherk
242 * ..
243 * .. Intrinsic Functions ..
244  INTRINSIC REAL, max, min
245 * ..
246 * .. Executable Statements ..
247 *
248  ulp = dlamch( 'Precision' )
249  ulpinv = realone / ulp
250  CALL zlaset( 'Full', m, m, zero, one, work, ldx )
251  CALL zherk( 'Upper', 'Conjugate transpose', m, m, -realone,
252  $ x, ldx, realone, work, ldx )
253  IF (m.GT.0) THEN
254  eps2 = max( ulp,
255  $ zlange( '1', m, m, work, ldx, rwork ) / dble( m ) )
256  ELSE
257  eps2 = ulp
258  END IF
259  r = min( p, m-p, q, m-q )
260 *
261 * Copy the matrix X to the array XF.
262 *
263  CALL zlacpy( 'Full', m, m, x, ldx, xf, ldx )
264 *
265 * Compute the CSD
266 *
267  CALL zuncsd( 'Y', 'Y', 'Y', 'Y', 'N', 'D', m, p, q, xf(1,1), ldx,
268  $ xf(1,q+1), ldx, xf(p+1,1), ldx, xf(p+1,q+1), ldx,
269  $ theta, u1, ldu1, u2, ldu2, v1t, ldv1t, v2t, ldv2t,
270  $ work, lwork, rwork, 17*(r+2), iwork, info )
271 *
272 * Compute X := diag(U1,U2)'*X*diag(V1,V2) - [D11 D12; D21 D22]
273 *
274  CALL zgemm( 'No transpose', 'Conjugate transpose', p, q, q, one,
275  $ x, ldx, v1t, ldv1t, zero, work, ldx )
276 *
277  CALL zgemm( 'Conjugate transpose', 'No transpose', p, q, p, one,
278  $ u1, ldu1, work, ldx, zero, x, ldx )
279 *
280  DO i = 1, min(p,q)-r
281  x(i,i) = x(i,i) - one
282  END DO
283  DO i = 1, r
284  x(min(p,q)-r+i,min(p,q)-r+i) =
285  $ x(min(p,q)-r+i,min(p,q)-r+i) - dcmplx( cos(theta(i)),
286  $ 0.0d0 )
287  END DO
288 *
289  CALL zgemm( 'No transpose', 'Conjugate transpose', p, m-q, m-q,
290  $ one, x(1,q+1), ldx, v2t, ldv2t, zero, work, ldx )
291 *
292  CALL zgemm( 'Conjugate transpose', 'No transpose', p, m-q, p,
293  $ one, u1, ldu1, work, ldx, zero, x(1,q+1), ldx )
294 *
295  DO i = 1, min(p,m-q)-r
296  x(p-i+1,m-i+1) = x(p-i+1,m-i+1) + one
297  END DO
298  DO i = 1, r
299  x(p-(min(p,m-q)-r)+1-i,m-(min(p,m-q)-r)+1-i) =
300  $ x(p-(min(p,m-q)-r)+1-i,m-(min(p,m-q)-r)+1-i) +
301  $ dcmplx( sin(theta(r-i+1)), 0.0d0 )
302  END DO
303 *
304  CALL zgemm( 'No transpose', 'Conjugate transpose', m-p, q, q, one,
305  $ x(p+1,1), ldx, v1t, ldv1t, zero, work, ldx )
306 *
307  CALL zgemm( 'Conjugate transpose', 'No transpose', m-p, q, m-p,
308  $ one, u2, ldu2, work, ldx, zero, x(p+1,1), ldx )
309 *
310  DO i = 1, min(m-p,q)-r
311  x(m-i+1,q-i+1) = x(m-i+1,q-i+1) - one
312  END DO
313  DO i = 1, r
314  x(m-(min(m-p,q)-r)+1-i,q-(min(m-p,q)-r)+1-i) =
315  $ x(m-(min(m-p,q)-r)+1-i,q-(min(m-p,q)-r)+1-i) -
316  $ dcmplx( sin(theta(r-i+1)), 0.0d0 )
317  END DO
318 *
319  CALL zgemm( 'No transpose', 'Conjugate transpose', m-p, m-q, m-q,
320  $ one, x(p+1,q+1), ldx, v2t, ldv2t, zero, work, ldx )
321 *
322  CALL zgemm( 'Conjugate transpose', 'No transpose', m-p, m-q, m-p,
323  $ one, u2, ldu2, work, ldx, zero, x(p+1,q+1), ldx )
324 *
325  DO i = 1, min(m-p,m-q)-r
326  x(p+i,q+i) = x(p+i,q+i) - one
327  END DO
328  DO i = 1, r
329  x(p+(min(m-p,m-q)-r)+i,q+(min(m-p,m-q)-r)+i) =
330  $ x(p+(min(m-p,m-q)-r)+i,q+(min(m-p,m-q)-r)+i) -
331  $ dcmplx( cos(theta(i)), 0.0d0 )
332  END DO
333 *
334 * Compute norm( U1'*X11*V1 - D11 ) / ( MAX(1,P,Q)*EPS2 ) .
335 *
336  resid = zlange( '1', p, q, x, ldx, rwork )
337  result( 1 ) = ( resid / REAL(MAX(1,P,Q)) ) / eps2
338 *
339 * Compute norm( U1'*X12*V2 - D12 ) / ( MAX(1,P,M-Q)*EPS2 ) .
340 *
341  resid = zlange( '1', p, m-q, x(1,q+1), ldx, rwork )
342  result( 2 ) = ( resid / REAL(MAX(1,P,M-Q)) ) / eps2
343 *
344 * Compute norm( U2'*X21*V1 - D21 ) / ( MAX(1,M-P,Q)*EPS2 ) .
345 *
346  resid = zlange( '1', m-p, q, x(p+1,1), ldx, rwork )
347  result( 3 ) = ( resid / REAL(MAX(1,M-P,Q)) ) / eps2
348 *
349 * Compute norm( U2'*X22*V2 - D22 ) / ( MAX(1,M-P,M-Q)*EPS2 ) .
350 *
351  resid = zlange( '1', m-p, m-q, x(p+1,q+1), ldx, rwork )
352  result( 4 ) = ( resid / REAL(MAX(1,M-P,M-Q)) ) / eps2
353 *
354 * Compute I - U1'*U1
355 *
356  CALL zlaset( 'Full', p, p, zero, one, work, ldu1 )
357  CALL zherk( 'Upper', 'Conjugate transpose', p, p, -realone,
358  $ u1, ldu1, realone, work, ldu1 )
359 *
360 * Compute norm( I - U'*U ) / ( MAX(1,P) * ULP ) .
361 *
362  resid = zlanhe( '1', 'Upper', p, work, ldu1, rwork )
363  result( 5 ) = ( resid / REAL(MAX(1,P)) ) / ulp
364 *
365 * Compute I - U2'*U2
366 *
367  CALL zlaset( 'Full', m-p, m-p, zero, one, work, ldu2 )
368  CALL zherk( 'Upper', 'Conjugate transpose', m-p, m-p, -realone,
369  $ u2, ldu2, realone, work, ldu2 )
370 *
371 * Compute norm( I - U2'*U2 ) / ( MAX(1,M-P) * ULP ) .
372 *
373  resid = zlanhe( '1', 'Upper', m-p, work, ldu2, rwork )
374  result( 6 ) = ( resid / REAL(MAX(1,M-P)) ) / ulp
375 *
376 * Compute I - V1T*V1T'
377 *
378  CALL zlaset( 'Full', q, q, zero, one, work, ldv1t )
379  CALL zherk( 'Upper', 'No transpose', q, q, -realone,
380  $ v1t, ldv1t, realone, work, ldv1t )
381 *
382 * Compute norm( I - V1T*V1T' ) / ( MAX(1,Q) * ULP ) .
383 *
384  resid = zlanhe( '1', 'Upper', q, work, ldv1t, rwork )
385  result( 7 ) = ( resid / REAL(MAX(1,Q)) ) / ulp
386 *
387 * Compute I - V2T*V2T'
388 *
389  CALL zlaset( 'Full', m-q, m-q, zero, one, work, ldv2t )
390  CALL zherk( 'Upper', 'No transpose', m-q, m-q, -realone,
391  $ v2t, ldv2t, realone, work, ldv2t )
392 *
393 * Compute norm( I - V2T*V2T' ) / ( MAX(1,M-Q) * ULP ) .
394 *
395  resid = zlanhe( '1', 'Upper', m-q, work, ldv2t, rwork )
396  result( 8 ) = ( resid / REAL(MAX(1,M-Q)) ) / ulp
397 *
398 * Check sorting
399 *
400  result(9) = realzero
401  DO i = 1, r
402  IF( theta(i).LT.realzero .OR. theta(i).GT.piover2 ) THEN
403  result(9) = ulpinv
404  END IF
405  IF( i.GT.1) THEN
406  IF ( theta(i).LT.theta(i-1) ) THEN
407  result(9) = ulpinv
408  END IF
409  END IF
410  END DO
411 *
412  return
413 *
414 * End of ZCSDTS
415 *
416  END
417