LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
dcsdts.f
Go to the documentation of this file.
1 *> \brief \b DCSDTS
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 DCSDTS( 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 * DOUBLE PRECISION 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 *> DCSDTS tests DORCSD, which, given an M-by-M partitioned orthogonal
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 DOUBLE PRECISION array, dimension (LDX,M)
76 *> The M-by-M matrix X.
77 *> \endverbatim
78 *>
79 *> \param[out] XF
80 *> \verbatim
81 *> XF is DOUBLE PRECISION array, dimension (LDX,M)
82 *> Details of the CSD of X, as returned by DORCSD;
83 *> see DORCSD 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 DOUBLE PRECISION array, dimension(LDU1,P)
96 *> The P-by-P orthogonal 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 DOUBLE PRECISION array, dimension(LDU2,M-P)
108 *> The (M-P)-by-(M-P) orthogonal 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 DOUBLE PRECISION array, dimension(LDV1T,Q)
120 *> The Q-by-Q orthogonal 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 DOUBLE PRECISION array, dimension(LDV2T,M-Q)
133 *> The (M-Q)-by-(M-Q) orthogonal 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 DORCSD 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 DOUBLE PRECISION 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 double_eig
201 *
202 * =====================================================================
203  SUBROUTINE dcsdts( 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  DOUBLE PRECISION 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  DOUBLE PRECISION zero, one
230  parameter( zero = 0.0d0, one = 1.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, dlange, dlansy
238  EXTERNAL dlamch, dlange, dlansy
239 * ..
240 * .. External Subroutines ..
241  EXTERNAL dgemm, dlacpy, dlaset, dorcsd, dsyrk
242 * ..
243 * .. Intrinsic Functions ..
244  INTRINSIC dble, max, min
245 * ..
246 * .. Executable Statements ..
247 *
248  ulp = dlamch( 'Precision' )
249  ulpinv = realone / ulp
250  CALL dlaset( 'Full', m, m, zero, one, work, ldx )
251  CALL dsyrk( 'Upper', 'Conjugate transpose', m, m, -one, x, ldx,
252  $ one, work, ldx )
253  IF (m.GT.0) THEN
254  eps2 = max( ulp,
255  $ dlange( '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 dlacpy( 'Full', m, m, x, ldx, xf, ldx )
264 *
265 * Compute the CSD
266 *
267  CALL dorcsd( '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, iwork, info )
271 *
272 * Compute X := diag(U1,U2)'*X*diag(V1,V2) - [D11 D12; D21 D22]
273 *
274  CALL dgemm( 'No transpose', 'Conjugate transpose', p, q, q, one,
275  $ x, ldx, v1t, ldv1t, zero, work, ldx )
276 *
277  CALL dgemm( '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) - cos(theta(i))
286  END DO
287 *
288  CALL dgemm( 'No transpose', 'Conjugate transpose', p, m-q, m-q,
289  $ one, x(1,q+1), ldx, v2t, ldv2t, zero, work, ldx )
290 *
291  CALL dgemm( 'Conjugate transpose', 'No transpose', p, m-q, p,
292  $ one, u1, ldu1, work, ldx, zero, x(1,q+1), ldx )
293 *
294  DO i = 1, min(p,m-q)-r
295  x(p-i+1,m-i+1) = x(p-i+1,m-i+1) + one
296  END DO
297  DO i = 1, r
298  x(p-(min(p,m-q)-r)+1-i,m-(min(p,m-q)-r)+1-i) =
299  $ x(p-(min(p,m-q)-r)+1-i,m-(min(p,m-q)-r)+1-i) +
300  $ sin(theta(r-i+1))
301  END DO
302 *
303  CALL dgemm( 'No transpose', 'Conjugate transpose', m-p, q, q, one,
304  $ x(p+1,1), ldx, v1t, ldv1t, zero, work, ldx )
305 *
306  CALL dgemm( 'Conjugate transpose', 'No transpose', m-p, q, m-p,
307  $ one, u2, ldu2, work, ldx, zero, x(p+1,1), ldx )
308 *
309  DO i = 1, min(m-p,q)-r
310  x(m-i+1,q-i+1) = x(m-i+1,q-i+1) - one
311  END DO
312  DO i = 1, r
313  x(m-(min(m-p,q)-r)+1-i,q-(min(m-p,q)-r)+1-i) =
314  $ x(m-(min(m-p,q)-r)+1-i,q-(min(m-p,q)-r)+1-i) -
315  $ sin(theta(r-i+1))
316  END DO
317 *
318  CALL dgemm( 'No transpose', 'Conjugate transpose', m-p, m-q, m-q,
319  $ one, x(p+1,q+1), ldx, v2t, ldv2t, zero, work, ldx )
320 *
321  CALL dgemm( 'Conjugate transpose', 'No transpose', m-p, m-q, m-p,
322  $ one, u2, ldu2, work, ldx, zero, x(p+1,q+1), ldx )
323 *
324  DO i = 1, min(m-p,m-q)-r
325  x(p+i,q+i) = x(p+i,q+i) - one
326  END DO
327  DO i = 1, r
328  x(p+(min(m-p,m-q)-r)+i,q+(min(m-p,m-q)-r)+i) =
329  $ x(p+(min(m-p,m-q)-r)+i,q+(min(m-p,m-q)-r)+i) -
330  $ cos(theta(i))
331  END DO
332 *
333 * Compute norm( U1'*X11*V1 - D11 ) / ( MAX(1,P,Q)*EPS2 ) .
334 *
335  resid = dlange( '1', p, q, x, ldx, rwork )
336  result( 1 ) = ( resid / dble(max(1,p,q)) ) / eps2
337 *
338 * Compute norm( U1'*X12*V2 - D12 ) / ( MAX(1,P,M-Q)*EPS2 ) .
339 *
340  resid = dlange( '1', p, m-q, x(1,q+1), ldx, rwork )
341  result( 2 ) = ( resid / dble(max(1,p,m-q)) ) / eps2
342 *
343 * Compute norm( U2'*X21*V1 - D21 ) / ( MAX(1,M-P,Q)*EPS2 ) .
344 *
345  resid = dlange( '1', m-p, q, x(p+1,1), ldx, rwork )
346  result( 3 ) = ( resid / dble(max(1,m-p,q)) ) / eps2
347 *
348 * Compute norm( U2'*X22*V2 - D22 ) / ( MAX(1,M-P,M-Q)*EPS2 ) .
349 *
350  resid = dlange( '1', m-p, m-q, x(p+1,q+1), ldx, rwork )
351  result( 4 ) = ( resid / dble(max(1,m-p,m-q)) ) / eps2
352 *
353 * Compute I - U1'*U1
354 *
355  CALL dlaset( 'Full', p, p, zero, one, work, ldu1 )
356  CALL dsyrk( 'Upper', 'Conjugate transpose', p, p, -one, u1, ldu1,
357  $ one, work, ldu1 )
358 *
359 * Compute norm( I - U'*U ) / ( MAX(1,P) * ULP ) .
360 *
361  resid = dlansy( '1', 'Upper', p, work, ldu1, rwork )
362  result( 5 ) = ( resid / dble(max(1,p)) ) / ulp
363 *
364 * Compute I - U2'*U2
365 *
366  CALL dlaset( 'Full', m-p, m-p, zero, one, work, ldu2 )
367  CALL dsyrk( 'Upper', 'Conjugate transpose', m-p, m-p, -one, u2,
368  $ ldu2, one, work, ldu2 )
369 *
370 * Compute norm( I - U2'*U2 ) / ( MAX(1,M-P) * ULP ) .
371 *
372  resid = dlansy( '1', 'Upper', m-p, work, ldu2, rwork )
373  result( 6 ) = ( resid / dble(max(1,m-p)) ) / ulp
374 *
375 * Compute I - V1T*V1T'
376 *
377  CALL dlaset( 'Full', q, q, zero, one, work, ldv1t )
378  CALL dsyrk( 'Upper', 'No transpose', q, q, -one, v1t, ldv1t, one,
379  $ work, ldv1t )
380 *
381 * Compute norm( I - V1T*V1T' ) / ( MAX(1,Q) * ULP ) .
382 *
383  resid = dlansy( '1', 'Upper', q, work, ldv1t, rwork )
384  result( 7 ) = ( resid / dble(max(1,q)) ) / ulp
385 *
386 * Compute I - V2T*V2T'
387 *
388  CALL dlaset( 'Full', m-q, m-q, zero, one, work, ldv2t )
389  CALL dsyrk( 'Upper', 'No transpose', m-q, m-q, -one, v2t, ldv2t,
390  $ one, work, ldv2t )
391 *
392 * Compute norm( I - V2T*V2T' ) / ( MAX(1,M-Q) * ULP ) .
393 *
394  resid = dlansy( '1', 'Upper', m-q, work, ldv2t, rwork )
395  result( 8 ) = ( resid / dble(max(1,m-q)) ) / ulp
396 *
397 * Check sorting
398 *
399  result(9) = realzero
400  DO i = 1, r
401  IF( theta(i).LT.realzero .OR. theta(i).GT.piover2 ) THEN
402  result(9) = ulpinv
403  END IF
404  IF( i.GT.1) THEN
405  IF ( theta(i).LT.theta(i-1) ) THEN
406  result(9) = ulpinv
407  END IF
408  END IF
409  END DO
410 *
411  return
412 *
413 * End of DCSDTS
414 *
415  END
416