LAPACK  3.5.0
LAPACK: Linear Algebra PACKage
 All Classes Files Functions Variables Typedefs Macros
ccsdts.f
Go to the documentation of this file.
1 *> \brief \b CCSDTS
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 CCSDTS( 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 * REAL RESULT( 15 ), RWORK( * ), THETA( * )
21 * COMPLEX 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 *> CCSDTS tests CUNCSD, 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 array, dimension (LDX,M)
91 *> The M-by-M matrix X.
92 *> \endverbatim
93 *>
94 *> \param[out] XF
95 *> \verbatim
96 *> XF is COMPLEX array, dimension (LDX,M)
97 *> Details of the CSD of X, as returned by CUNCSD;
98 *> see CUNCSD 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 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 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 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 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 REAL 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 CUNCSD 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 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 REAL array
185 *> \endverbatim
186 *>
187 *> \param[out] RESULT
188 *> \verbatim
189 *> RESULT is REAL 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 2011
224 *
225 *> \ingroup complex_eig
226 *
227 * =====================================================================
228  SUBROUTINE ccsdts( 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.4.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 2011
236 *
237 * .. Scalar Arguments ..
238  INTEGER ldx, ldu1, ldu2, ldv1t, ldv2t, lwork, m, p, q
239 * ..
240 * .. Array Arguments ..
241  INTEGER iwork( * )
242  REAL result( 15 ), rwork( * ), theta( * )
243  COMPLEX u1( ldu1, * ), u2( ldu2, * ), v1t( ldv1t, * ),
244  $ v2t( ldv2t, * ), work( lwork ), x( ldx, * ),
245  $ xf( ldx, * )
246 * ..
247 *
248 * =====================================================================
249 *
250 * .. Parameters ..
251  REAL piover2, realone, realzero
252  parameter( piover2 = 1.57079632679489662e0,
253  $ realone = 1.0e0, realzero = 0.0e0 )
254  COMPLEX zero, one
255  parameter( zero = (0.0e0,0.0e0), one = (1.0e0,0.0e0) )
256 * ..
257 * .. Local Scalars ..
258  INTEGER i, info, r
259  REAL eps2, resid, ulp, ulpinv
260 * ..
261 * .. External Functions ..
262  REAL slamch, clange, clanhe
263  EXTERNAL slamch, clange, clanhe
264 * ..
265 * .. External Subroutines ..
266  EXTERNAL cgemm, cherk, clacpy, claset, cuncsd, cuncsd2by1
267 * ..
268 * .. Intrinsic Functions ..
269  INTRINSIC cmplx, cos, max, min, REAL, sin
270 * ..
271 * .. Executable Statements ..
272 *
273  ulp = slamch( 'Precision' )
274  ulpinv = realone / ulp
275 *
276 * The first half of the routine checks the 2-by-2 CSD
277 *
278  CALL claset( 'Full', m, m, zero, one, work, ldx )
279  CALL cherk( 'Upper', 'Conjugate transpose', m, m, -realone,
280  $ x, ldx, realone, work, ldx )
281  IF (m.GT.0) THEN
282  eps2 = max( ulp,
283  $ clange( '1', m, m, work, ldx, rwork ) / REAL( M ) )
284  ELSE
285  eps2 = ulp
286  END IF
287  r = min( p, m-p, q, m-q )
288 *
289 * Copy the matrix X to the array XF.
290 *
291  CALL clacpy( 'Full', m, m, x, ldx, xf, ldx )
292 *
293 * Compute the CSD
294 *
295  CALL cuncsd( 'Y', 'Y', 'Y', 'Y', 'N', 'D', m, p, q, xf(1,1), ldx,
296  $ xf(1,q+1), ldx, xf(p+1,1), ldx, xf(p+1,q+1), ldx,
297  $ theta, u1, ldu1, u2, ldu2, v1t, ldv1t, v2t, ldv2t,
298  $ work, lwork, rwork, 17*(r+2), iwork, info )
299 *
300 * Compute XF := diag(U1,U2)'*X*diag(V1,V2) - [D11 D12; D21 D22]
301 *
302  CALL clacpy( 'Full', m, m, x, ldx, xf, ldx )
303 *
304  CALL cgemm( 'No transpose', 'Conjugate transpose', p, q, q, one,
305  $ xf, ldx, v1t, ldv1t, zero, work, ldx )
306 *
307  CALL cgemm( 'Conjugate transpose', 'No transpose', p, q, p, one,
308  $ u1, ldu1, work, ldx, zero, xf, ldx )
309 *
310  DO i = 1, min(p,q)-r
311  xf(i,i) = xf(i,i) - one
312  END DO
313  DO i = 1, r
314  xf(min(p,q)-r+i,min(p,q)-r+i) =
315  $ xf(min(p,q)-r+i,min(p,q)-r+i) - cmplx( cos(theta(i)),
316  $ 0.0e0 )
317  END DO
318 *
319  CALL cgemm( 'No transpose', 'Conjugate transpose', p, m-q, m-q,
320  $ one, xf(1,q+1), ldx, v2t, ldv2t, zero, work, ldx )
321 *
322  CALL cgemm( 'Conjugate transpose', 'No transpose', p, m-q, p,
323  $ one, u1, ldu1, work, ldx, zero, xf(1,q+1), ldx )
324 *
325  DO i = 1, min(p,m-q)-r
326  xf(p-i+1,m-i+1) = xf(p-i+1,m-i+1) + one
327  END DO
328  DO i = 1, r
329  xf(p-(min(p,m-q)-r)+1-i,m-(min(p,m-q)-r)+1-i) =
330  $ xf(p-(min(p,m-q)-r)+1-i,m-(min(p,m-q)-r)+1-i) +
331  $ cmplx( sin(theta(r-i+1)), 0.0e0 )
332  END DO
333 *
334  CALL cgemm( 'No transpose', 'Conjugate transpose', m-p, q, q, one,
335  $ xf(p+1,1), ldx, v1t, ldv1t, zero, work, ldx )
336 *
337  CALL cgemm( 'Conjugate transpose', 'No transpose', m-p, q, m-p,
338  $ one, u2, ldu2, work, ldx, zero, xf(p+1,1), ldx )
339 *
340  DO i = 1, min(m-p,q)-r
341  xf(m-i+1,q-i+1) = xf(m-i+1,q-i+1) - one
342  END DO
343  DO i = 1, r
344  xf(m-(min(m-p,q)-r)+1-i,q-(min(m-p,q)-r)+1-i) =
345  $ xf(m-(min(m-p,q)-r)+1-i,q-(min(m-p,q)-r)+1-i) -
346  $ cmplx( sin(theta(r-i+1)), 0.0e0 )
347  END DO
348 *
349  CALL cgemm( 'No transpose', 'Conjugate transpose', m-p, m-q, m-q,
350  $ one, xf(p+1,q+1), ldx, v2t, ldv2t, zero, work, ldx )
351 *
352  CALL cgemm( 'Conjugate transpose', 'No transpose', m-p, m-q, m-p,
353  $ one, u2, ldu2, work, ldx, zero, xf(p+1,q+1), ldx )
354 *
355  DO i = 1, min(m-p,m-q)-r
356  xf(p+i,q+i) = xf(p+i,q+i) - one
357  END DO
358  DO i = 1, r
359  xf(p+(min(m-p,m-q)-r)+i,q+(min(m-p,m-q)-r)+i) =
360  $ xf(p+(min(m-p,m-q)-r)+i,q+(min(m-p,m-q)-r)+i) -
361  $ cmplx( cos(theta(i)), 0.0e0 )
362  END DO
363 *
364 * Compute norm( U1'*X11*V1 - D11 ) / ( MAX(1,P,Q)*EPS2 ) .
365 *
366  resid = clange( '1', p, q, xf, ldx, rwork )
367  result( 1 ) = ( resid / REAL(MAX(1,P,Q)) ) / eps2
368 *
369 * Compute norm( U1'*X12*V2 - D12 ) / ( MAX(1,P,M-Q)*EPS2 ) .
370 *
371  resid = clange( '1', p, m-q, xf(1,q+1), ldx, rwork )
372  result( 2 ) = ( resid / REAL(MAX(1,P,M-Q)) ) / eps2
373 *
374 * Compute norm( U2'*X21*V1 - D21 ) / ( MAX(1,M-P,Q)*EPS2 ) .
375 *
376  resid = clange( '1', m-p, q, xf(p+1,1), ldx, rwork )
377  result( 3 ) = ( resid / REAL(MAX(1,M-P,Q)) ) / eps2
378 *
379 * Compute norm( U2'*X22*V2 - D22 ) / ( MAX(1,M-P,M-Q)*EPS2 ) .
380 *
381  resid = clange( '1', m-p, m-q, xf(p+1,q+1), ldx, rwork )
382  result( 4 ) = ( resid / REAL(MAX(1,M-P,M-Q)) ) / eps2
383 *
384 * Compute I - U1'*U1
385 *
386  CALL claset( 'Full', p, p, zero, one, work, ldu1 )
387  CALL cherk( 'Upper', 'Conjugate transpose', p, p, -realone,
388  $ u1, ldu1, realone, work, ldu1 )
389 *
390 * Compute norm( I - U'*U ) / ( MAX(1,P) * ULP ) .
391 *
392  resid = clanhe( '1', 'Upper', p, work, ldu1, rwork )
393  result( 5 ) = ( resid / REAL(MAX(1,P)) ) / ulp
394 *
395 * Compute I - U2'*U2
396 *
397  CALL claset( 'Full', m-p, m-p, zero, one, work, ldu2 )
398  CALL cherk( 'Upper', 'Conjugate transpose', m-p, m-p, -realone,
399  $ u2, ldu2, realone, work, ldu2 )
400 *
401 * Compute norm( I - U2'*U2 ) / ( MAX(1,M-P) * ULP ) .
402 *
403  resid = clanhe( '1', 'Upper', m-p, work, ldu2, rwork )
404  result( 6 ) = ( resid / REAL(MAX(1,M-P)) ) / ulp
405 *
406 * Compute I - V1T*V1T'
407 *
408  CALL claset( 'Full', q, q, zero, one, work, ldv1t )
409  CALL cherk( 'Upper', 'No transpose', q, q, -realone,
410  $ v1t, ldv1t, realone, work, ldv1t )
411 *
412 * Compute norm( I - V1T*V1T' ) / ( MAX(1,Q) * ULP ) .
413 *
414  resid = clanhe( '1', 'Upper', q, work, ldv1t, rwork )
415  result( 7 ) = ( resid / REAL(MAX(1,Q)) ) / ulp
416 *
417 * Compute I - V2T*V2T'
418 *
419  CALL claset( 'Full', m-q, m-q, zero, one, work, ldv2t )
420  CALL cherk( 'Upper', 'No transpose', m-q, m-q, -realone,
421  $ v2t, ldv2t, realone, work, ldv2t )
422 *
423 * Compute norm( I - V2T*V2T' ) / ( MAX(1,M-Q) * ULP ) .
424 *
425  resid = clanhe( '1', 'Upper', m-q, work, ldv2t, rwork )
426  result( 8 ) = ( resid / REAL(MAX(1,M-Q)) ) / ulp
427 *
428 * Check sorting
429 *
430  result( 9 ) = realzero
431  DO i = 1, r
432  IF( theta(i).LT.realzero .OR. theta(i).GT.piover2 ) THEN
433  result( 9 ) = ulpinv
434  END IF
435  IF( i.GT.1) THEN
436  IF ( theta(i).LT.theta(i-1) ) THEN
437  result( 9 ) = ulpinv
438  END IF
439  END IF
440  END DO
441 *
442 * The second half of the routine checks the 2-by-1 CSD
443 *
444  CALL claset( 'Full', q, q, zero, one, work, ldx )
445  CALL cherk( 'Upper', 'Conjugate transpose', q, m, -realone,
446  $ x, ldx, realone, work, ldx )
447  IF (m.GT.0) THEN
448  eps2 = max( ulp,
449  $ clange( '1', q, q, work, ldx, rwork ) / REAL( M ) )
450  ELSE
451  eps2 = ulp
452  END IF
453  r = min( p, m-p, q, m-q )
454 *
455 * Copy the matrix X to the array XF.
456 *
457  CALL clacpy( 'Full', m, q, x, ldx, xf, ldx )
458 *
459 * Compute the CSD
460 *
461  CALL cuncsd2by1( 'Y', 'Y', 'Y', m, p, q, xf(1,1), ldx, xf(p+1,1),
462  $ ldx, theta, u1, ldu1, u2, ldu2, v1t, ldv1t, work,
463  $ lwork, rwork, 17*(r+2), iwork, info )
464 *
465 * Compute [X11;X21] := diag(U1,U2)'*[X11;X21]*V1 - [D11;D21]
466 *
467  CALL cgemm( 'No transpose', 'Conjugate transpose', p, q, q, one,
468  $ x, ldx, v1t, ldv1t, zero, work, ldx )
469 *
470  CALL cgemm( 'Conjugate transpose', 'No transpose', p, q, p, one,
471  $ u1, ldu1, work, ldx, zero, x, ldx )
472 *
473  DO i = 1, min(p,q)-r
474  x(i,i) = x(i,i) - one
475  END DO
476  DO i = 1, r
477  x(min(p,q)-r+i,min(p,q)-r+i) =
478  $ x(min(p,q)-r+i,min(p,q)-r+i) - cmplx( cos(theta(i)),
479  $ 0.0e0 )
480  END DO
481 *
482  CALL cgemm( 'No transpose', 'Conjugate transpose', m-p, q, q, one,
483  $ x(p+1,1), ldx, v1t, ldv1t, zero, work, ldx )
484 *
485  CALL cgemm( 'Conjugate transpose', 'No transpose', m-p, q, m-p,
486  $ one, u2, ldu2, work, ldx, zero, x(p+1,1), ldx )
487 *
488  DO i = 1, min(m-p,q)-r
489  x(m-i+1,q-i+1) = x(m-i+1,q-i+1) - one
490  END DO
491  DO i = 1, r
492  x(m-(min(m-p,q)-r)+1-i,q-(min(m-p,q)-r)+1-i) =
493  $ x(m-(min(m-p,q)-r)+1-i,q-(min(m-p,q)-r)+1-i) -
494  $ cmplx( sin(theta(r-i+1)), 0.0e0 )
495  END DO
496 *
497 * Compute norm( U1'*X11*V1 - D11 ) / ( MAX(1,P,Q)*EPS2 ) .
498 *
499  resid = clange( '1', p, q, x, ldx, rwork )
500  result( 10 ) = ( resid / REAL(MAX(1,P,Q)) ) / eps2
501 *
502 * Compute norm( U2'*X21*V1 - D21 ) / ( MAX(1,M-P,Q)*EPS2 ) .
503 *
504  resid = clange( '1', m-p, q, x(p+1,1), ldx, rwork )
505  result( 11 ) = ( resid / REAL(MAX(1,M-P,Q)) ) / eps2
506 *
507 * Compute I - U1'*U1
508 *
509  CALL claset( 'Full', p, p, zero, one, work, ldu1 )
510  CALL cherk( 'Upper', 'Conjugate transpose', p, p, -realone,
511  $ u1, ldu1, realone, work, ldu1 )
512 *
513 * Compute norm( I - U1'*U1 ) / ( MAX(1,P) * ULP ) .
514 *
515  resid = clanhe( '1', 'Upper', p, work, ldu1, rwork )
516  result( 12 ) = ( resid / REAL(MAX(1,P)) ) / ulp
517 *
518 * Compute I - U2'*U2
519 *
520  CALL claset( 'Full', m-p, m-p, zero, one, work, ldu2 )
521  CALL cherk( 'Upper', 'Conjugate transpose', m-p, m-p, -realone,
522  $ u2, ldu2, realone, work, ldu2 )
523 *
524 * Compute norm( I - U2'*U2 ) / ( MAX(1,M-P) * ULP ) .
525 *
526  resid = clanhe( '1', 'Upper', m-p, work, ldu2, rwork )
527  result( 13 ) = ( resid / REAL(MAX(1,M-P)) ) / ulp
528 *
529 * Compute I - V1T*V1T'
530 *
531  CALL claset( 'Full', q, q, zero, one, work, ldv1t )
532  CALL cherk( 'Upper', 'No transpose', q, q, -realone,
533  $ v1t, ldv1t, realone, work, ldv1t )
534 *
535 * Compute norm( I - V1T*V1T' ) / ( MAX(1,Q) * ULP ) .
536 *
537  resid = clanhe( '1', 'Upper', q, work, ldv1t, rwork )
538  result( 14 ) = ( resid / REAL(MAX(1,Q)) ) / ulp
539 *
540 * Check sorting
541 *
542  result( 15 ) = realzero
543  DO i = 1, r
544  IF( theta(i).LT.realzero .OR. theta(i).GT.piover2 ) THEN
545  result( 15 ) = ulpinv
546  END IF
547  IF( i.GT.1) THEN
548  IF ( theta(i).LT.theta(i-1) ) THEN
549  result( 15 ) = ulpinv
550  END IF
551  END IF
552  END DO
553 *
554  RETURN
555 *
556 * End of CCSDTS
557 *
558  END
subroutine claset(UPLO, M, N, ALPHA, BETA, A, LDA)
CLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
Definition: claset.f:107
recursive subroutine cuncsd(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)
CUNCSD
Definition: cuncsd.f:316
REAL function clanhe(NORM, UPLO, N, A, LDA, WORK)
CLANHE returns the value of the 1-norm, or the Frobenius norm, or the infinity norm, or the element of largest absolute value of a complex Hermitian matrix.
Definition: clanhe.f:125
REAL function slamch(CMACH)
SLAMCH
Definition: slamch.f:68
subroutine ccsdts(M, P, Q, X, XF, LDX, U1, LDU1, U2, LDU2, V1T, LDV1T, V2T, LDV2T, THETA, IWORK, WORK, LWORK, RWORK, RESULT)
CCSDTS
Definition: ccsdts.f:228
subroutine cuncsd2by1(JOBU1, JOBU2, JOBV1T, M, P, Q, X11, LDX11, X21, LDX21, THETA, U1, LDU1, U2, LDU2, V1T, LDV1T, WORK, LWORK, RWORK, LRWORK, IWORK, INFO)
CUNCSD2BY1
Definition: cuncsd2by1.f:260
REAL function clange(NORM, M, N, A, LDA, WORK)
CLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
Definition: clange.f:116
subroutine clacpy(UPLO, M, N, A, LDA, B, LDB)
CLACPY copies all or part of one two-dimensional array to another.
Definition: clacpy.f:104
ELF f x
Definition: testslamch:1
subroutine cgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
CGEMM
Definition: cgemm.f:188
subroutine cherk(UPLO, TRANS, N, K, ALPHA, A, LDA, BETA, C, LDC)
CHERK
Definition: cherk.f:174