LAPACK  3.10.1
LAPACK: Linear Algebra PACKage
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( 15 ), 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 *>
51 *> and also DORCSD2BY1, 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 DOUBLE PRECISION array, dimension (LDX,M)
91 *> The M-by-M matrix X.
92 *> \endverbatim
93 *>
94 *> \param[out] XF
95 *> \verbatim
96 *> XF is DOUBLE PRECISION array, dimension (LDX,M)
97 *> Details of the CSD of X, as returned by DORCSD;
98 *> see DORCSD 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 DOUBLE PRECISION array, dimension(LDU1,P)
111 *> The P-by-P orthogonal 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 DOUBLE PRECISION array, dimension(LDU2,M-P)
123 *> The (M-P)-by-(M-P) orthogonal 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 DOUBLE PRECISION array, dimension(LDV1T,Q)
135 *> The Q-by-Q orthogonal 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 DOUBLE PRECISION array, dimension(LDV2T,M-Q)
148 *> The (M-Q)-by-(M-Q) orthogonal 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 DORCSD 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 DOUBLE PRECISION 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 double_eig
224 *
225 * =====================================================================
226  SUBROUTINE dcsdts( 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  DOUBLE PRECISION 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  DOUBLE PRECISION ZERO, ONE
251  parameter( zero = 0.0d0, one = 1.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, DLANGE, DLANSY
261  EXTERNAL DLAMCH, DLANGE, DLANSY
262 * ..
263 * .. External Subroutines ..
264  EXTERNAL dgemm, dlacpy, dlaset, dorcsd, dorcsd2by1,
265  $ dsyrk
266 * ..
267 * .. Intrinsic Functions ..
268  INTRINSIC cos, dble, max, min, 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 dlaset( 'Full', m, m, zero, one, work, ldx )
278  CALL dsyrk( 'Upper', 'Conjugate transpose', m, m, -one, x, ldx,
279  $ one, work, ldx )
280  IF (m.GT.0) THEN
281  eps2 = max( ulp,
282  $ dlange( '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 dlacpy( 'Full', m, m, x, ldx, xf, ldx )
291 *
292 * Compute the CSD
293 *
294  CALL dorcsd( '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, iwork, info )
298 *
299 * Compute XF := diag(U1,U2)'*X*diag(V1,V2) - [D11 D12; D21 D22]
300 *
301  CALL dlacpy( 'Full', m, m, x, ldx, xf, ldx )
302 *
303  CALL dgemm( 'No transpose', 'Conjugate transpose', p, q, q, one,
304  $ xf, ldx, v1t, ldv1t, zero, work, ldx )
305 *
306  CALL dgemm( '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) - cos(theta(i))
315  END DO
316 *
317  CALL dgemm( 'No transpose', 'Conjugate transpose', p, m-q, m-q,
318  $ one, xf(1,q+1), ldx, v2t, ldv2t, zero, work, ldx )
319 *
320  CALL dgemm( 'Conjugate transpose', 'No transpose', p, m-q, p,
321  $ one, u1, ldu1, work, ldx, zero, xf(1,q+1), ldx )
322 *
323  DO i = 1, min(p,m-q)-r
324  xf(p-i+1,m-i+1) = xf(p-i+1,m-i+1) + one
325  END DO
326  DO i = 1, r
327  xf(p-(min(p,m-q)-r)+1-i,m-(min(p,m-q)-r)+1-i) =
328  $ xf(p-(min(p,m-q)-r)+1-i,m-(min(p,m-q)-r)+1-i) +
329  $ sin(theta(r-i+1))
330  END DO
331 *
332  CALL dgemm( 'No transpose', 'Conjugate transpose', m-p, q, q, one,
333  $ xf(p+1,1), ldx, v1t, ldv1t, zero, work, ldx )
334 *
335  CALL dgemm( 'Conjugate transpose', 'No transpose', m-p, q, m-p,
336  $ one, u2, ldu2, work, ldx, zero, xf(p+1,1), ldx )
337 *
338  DO i = 1, min(m-p,q)-r
339  xf(m-i+1,q-i+1) = xf(m-i+1,q-i+1) - one
340  END DO
341  DO i = 1, r
342  xf(m-(min(m-p,q)-r)+1-i,q-(min(m-p,q)-r)+1-i) =
343  $ xf(m-(min(m-p,q)-r)+1-i,q-(min(m-p,q)-r)+1-i) -
344  $ sin(theta(r-i+1))
345  END DO
346 *
347  CALL dgemm( 'No transpose', 'Conjugate transpose', m-p, m-q, m-q,
348  $ one, xf(p+1,q+1), ldx, v2t, ldv2t, zero, work, ldx )
349 *
350  CALL dgemm( 'Conjugate transpose', 'No transpose', m-p, m-q, m-p,
351  $ one, u2, ldu2, work, ldx, zero, xf(p+1,q+1), ldx )
352 *
353  DO i = 1, min(m-p,m-q)-r
354  xf(p+i,q+i) = xf(p+i,q+i) - one
355  END DO
356  DO i = 1, r
357  xf(p+(min(m-p,m-q)-r)+i,q+(min(m-p,m-q)-r)+i) =
358  $ xf(p+(min(m-p,m-q)-r)+i,q+(min(m-p,m-q)-r)+i) -
359  $ cos(theta(i))
360  END DO
361 *
362 * Compute norm( U1'*X11*V1 - D11 ) / ( MAX(1,P,Q)*EPS2 ) .
363 *
364  resid = dlange( '1', p, q, xf, ldx, rwork )
365  result( 1 ) = ( resid / dble(max(1,p,q)) ) / eps2
366 *
367 * Compute norm( U1'*X12*V2 - D12 ) / ( MAX(1,P,M-Q)*EPS2 ) .
368 *
369  resid = dlange( '1', p, m-q, xf(1,q+1), ldx, rwork )
370  result( 2 ) = ( resid / dble(max(1,p,m-q)) ) / eps2
371 *
372 * Compute norm( U2'*X21*V1 - D21 ) / ( MAX(1,M-P,Q)*EPS2 ) .
373 *
374  resid = dlange( '1', m-p, q, xf(p+1,1), ldx, rwork )
375  result( 3 ) = ( resid / dble(max(1,m-p,q)) ) / eps2
376 *
377 * Compute norm( U2'*X22*V2 - D22 ) / ( MAX(1,M-P,M-Q)*EPS2 ) .
378 *
379  resid = dlange( '1', m-p, m-q, xf(p+1,q+1), ldx, rwork )
380  result( 4 ) = ( resid / dble(max(1,m-p,m-q)) ) / eps2
381 *
382 * Compute I - U1'*U1
383 *
384  CALL dlaset( 'Full', p, p, zero, one, work, ldu1 )
385  CALL dsyrk( 'Upper', 'Conjugate transpose', p, p, -one, u1, ldu1,
386  $ one, work, ldu1 )
387 *
388 * Compute norm( I - U'*U ) / ( MAX(1,P) * ULP ) .
389 *
390  resid = dlansy( '1', 'Upper', p, work, ldu1, rwork )
391  result( 5 ) = ( resid / dble(max(1,p)) ) / ulp
392 *
393 * Compute I - U2'*U2
394 *
395  CALL dlaset( 'Full', m-p, m-p, zero, one, work, ldu2 )
396  CALL dsyrk( 'Upper', 'Conjugate transpose', m-p, m-p, -one, u2,
397  $ ldu2, one, work, ldu2 )
398 *
399 * Compute norm( I - U2'*U2 ) / ( MAX(1,M-P) * ULP ) .
400 *
401  resid = dlansy( '1', 'Upper', m-p, work, ldu2, rwork )
402  result( 6 ) = ( resid / dble(max(1,m-p)) ) / ulp
403 *
404 * Compute I - V1T*V1T'
405 *
406  CALL dlaset( 'Full', q, q, zero, one, work, ldv1t )
407  CALL dsyrk( 'Upper', 'No transpose', q, q, -one, v1t, ldv1t, one,
408  $ work, ldv1t )
409 *
410 * Compute norm( I - V1T*V1T' ) / ( MAX(1,Q) * ULP ) .
411 *
412  resid = dlansy( '1', 'Upper', q, work, ldv1t, rwork )
413  result( 7 ) = ( resid / dble(max(1,q)) ) / ulp
414 *
415 * Compute I - V2T*V2T'
416 *
417  CALL dlaset( 'Full', m-q, m-q, zero, one, work, ldv2t )
418  CALL dsyrk( 'Upper', 'No transpose', m-q, m-q, -one, v2t, ldv2t,
419  $ one, work, ldv2t )
420 *
421 * Compute norm( I - V2T*V2T' ) / ( MAX(1,M-Q) * ULP ) .
422 *
423  resid = dlansy( '1', 'Upper', m-q, work, ldv2t, rwork )
424  result( 8 ) = ( resid / dble(max(1,m-q)) ) / ulp
425 *
426 * Check sorting
427 *
428  result( 9 ) = realzero
429  DO i = 1, r
430  IF( theta(i).LT.realzero .OR. theta(i).GT.piover2 ) THEN
431  result( 9 ) = ulpinv
432  END IF
433  IF( i.GT.1 ) THEN
434  IF ( theta(i).LT.theta(i-1) ) THEN
435  result( 9 ) = ulpinv
436  END IF
437  END IF
438  END DO
439 *
440 * The second half of the routine checks the 2-by-1 CSD
441 *
442  CALL dlaset( 'Full', q, q, zero, one, work, ldx )
443  CALL dsyrk( 'Upper', 'Conjugate transpose', q, m, -one, x, ldx,
444  $ one, work, ldx )
445  IF( m.GT.0 ) THEN
446  eps2 = max( ulp,
447  $ dlange( '1', q, q, work, ldx, rwork ) / dble( m ) )
448  ELSE
449  eps2 = ulp
450  END IF
451  r = min( p, m-p, q, m-q )
452 *
453 * Copy the matrix [ X11; X21 ] to the array XF.
454 *
455  CALL dlacpy( 'Full', m, q, x, ldx, xf, ldx )
456 *
457 * Compute the CSD
458 *
459  CALL dorcsd2by1( 'Y', 'Y', 'Y', m, p, q, xf(1,1), ldx, xf(p+1,1),
460  $ ldx, theta, u1, ldu1, u2, ldu2, v1t, ldv1t, work,
461  $ lwork, iwork, info )
462 *
463 * Compute [X11;X21] := diag(U1,U2)'*[X11;X21]*V1 - [D11;D21]
464 *
465  CALL dgemm( 'No transpose', 'Conjugate transpose', p, q, q, one,
466  $ x, ldx, v1t, ldv1t, zero, work, ldx )
467 *
468  CALL dgemm( 'Conjugate transpose', 'No transpose', p, q, p, one,
469  $ u1, ldu1, work, ldx, zero, x, ldx )
470 *
471  DO i = 1, min(p,q)-r
472  x(i,i) = x(i,i) - one
473  END DO
474  DO i = 1, r
475  x(min(p,q)-r+i,min(p,q)-r+i) =
476  $ x(min(p,q)-r+i,min(p,q)-r+i) - cos(theta(i))
477  END DO
478 *
479  CALL dgemm( 'No transpose', 'Conjugate transpose', m-p, q, q, one,
480  $ x(p+1,1), ldx, v1t, ldv1t, zero, work, ldx )
481 *
482  CALL dgemm( 'Conjugate transpose', 'No transpose', m-p, q, m-p,
483  $ one, u2, ldu2, work, ldx, zero, x(p+1,1), ldx )
484 *
485  DO i = 1, min(m-p,q)-r
486  x(m-i+1,q-i+1) = x(m-i+1,q-i+1) - one
487  END DO
488  DO i = 1, r
489  x(m-(min(m-p,q)-r)+1-i,q-(min(m-p,q)-r)+1-i) =
490  $ x(m-(min(m-p,q)-r)+1-i,q-(min(m-p,q)-r)+1-i) -
491  $ sin(theta(r-i+1))
492  END DO
493 *
494 * Compute norm( U1'*X11*V1 - D11 ) / ( MAX(1,P,Q)*EPS2 ) .
495 *
496  resid = dlange( '1', p, q, x, ldx, rwork )
497  result( 10 ) = ( resid / dble(max(1,p,q)) ) / eps2
498 *
499 * Compute norm( U2'*X21*V1 - D21 ) / ( MAX(1,M-P,Q)*EPS2 ) .
500 *
501  resid = dlange( '1', m-p, q, x(p+1,1), ldx, rwork )
502  result( 11 ) = ( resid / dble(max(1,m-p,q)) ) / eps2
503 *
504 * Compute I - U1'*U1
505 *
506  CALL dlaset( 'Full', p, p, zero, one, work, ldu1 )
507  CALL dsyrk( 'Upper', 'Conjugate transpose', p, p, -one, u1, ldu1,
508  $ one, work, ldu1 )
509 *
510 * Compute norm( I - U1'*U1 ) / ( MAX(1,P) * ULP ) .
511 *
512  resid = dlansy( '1', 'Upper', p, work, ldu1, rwork )
513  result( 12 ) = ( resid / dble(max(1,p)) ) / ulp
514 *
515 * Compute I - U2'*U2
516 *
517  CALL dlaset( 'Full', m-p, m-p, zero, one, work, ldu2 )
518  CALL dsyrk( 'Upper', 'Conjugate transpose', m-p, m-p, -one, u2,
519  $ ldu2, one, work, ldu2 )
520 *
521 * Compute norm( I - U2'*U2 ) / ( MAX(1,M-P) * ULP ) .
522 *
523  resid = dlansy( '1', 'Upper', m-p, work, ldu2, rwork )
524  result( 13 ) = ( resid / dble(max(1,m-p)) ) / ulp
525 *
526 * Compute I - V1T*V1T'
527 *
528  CALL dlaset( 'Full', q, q, zero, one, work, ldv1t )
529  CALL dsyrk( 'Upper', 'No transpose', q, q, -one, v1t, ldv1t, one,
530  $ work, ldv1t )
531 *
532 * Compute norm( I - V1T*V1T' ) / ( MAX(1,Q) * ULP ) .
533 *
534  resid = dlansy( '1', 'Upper', q, work, ldv1t, rwork )
535  result( 14 ) = ( resid / dble(max(1,q)) ) / ulp
536 *
537 * Check sorting
538 *
539  result( 15 ) = realzero
540  DO i = 1, r
541  IF( theta(i).LT.realzero .OR. theta(i).GT.piover2 ) THEN
542  result( 15 ) = ulpinv
543  END IF
544  IF( i.GT.1 ) THEN
545  IF ( theta(i).LT.theta(i-1) ) THEN
546  result( 15 ) = ulpinv
547  END IF
548  END IF
549  END DO
550 *
551  RETURN
552 *
553 * End of DCSDTS
554 *
555  END
556 
subroutine dlacpy(UPLO, M, N, A, LDA, B, LDB)
DLACPY copies all or part of one two-dimensional array to another.
Definition: dlacpy.f:103
subroutine dlaset(UPLO, M, N, ALPHA, BETA, A, LDA)
DLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition: dlaset.f:110
subroutine dsyrk(UPLO, TRANS, N, K, ALPHA, A, LDA, BETA, C, LDC)
DSYRK
Definition: dsyrk.f:169
subroutine dgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
DGEMM
Definition: dgemm.f:187
subroutine dcsdts(M, P, Q, X, XF, LDX, U1, LDU1, U2, LDU2, V1T, LDV1T, V2T, LDV2T, THETA, IWORK, WORK, LWORK, RWORK, RESULT)
DCSDTS
Definition: dcsdts.f:229
recursive subroutine dorcsd(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, IWORK, INFO)
DORCSD
Definition: dorcsd.f:300
subroutine dorcsd2by1(JOBU1, JOBU2, JOBV1T, M, P, Q, X11, LDX11, X21, LDX21, THETA, U1, LDU1, U2, LDU2, V1T, LDV1T, WORK, LWORK, IWORK, INFO)
DORCSD2BY1
Definition: dorcsd2by1.f:233