LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
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 2015
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.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  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,
267  $ cuncsd2by1
268 * ..
269 * .. Intrinsic Functions ..
270  INTRINSIC cmplx, cos, max, min, REAL, SIN
271 * ..
272 * .. Executable Statements ..
273 *
274  ulp = slamch( 'Precision' )
275  ulpinv = realone / ulp
276 *
277 * The first half of the routine checks the 2-by-2 CSD
278 *
279  CALL claset( 'Full', m, m, zero, one, work, ldx )
280  CALL cherk( 'Upper', 'Conjugate transpose', m, m, -realone,
281  $ x, ldx, realone, work, ldx )
282  IF (m.GT.0) THEN
283  eps2 = max( ulp,
284  $ clange( '1', m, m, work, ldx, rwork ) / REAL( 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 clacpy( 'Full', m, m, x, ldx, xf, ldx )
293 *
294 * Compute the CSD
295 *
296  CALL cuncsd( '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 clacpy( 'Full', m, m, x, ldx, xf, ldx )
304 *
305  CALL cgemm( 'No transpose', 'Conjugate transpose', p, q, q, one,
306  $ xf, ldx, v1t, ldv1t, zero, work, ldx )
307 *
308  CALL cgemm( '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) - cmplx( cos(theta(i)),
317  $ 0.0e0 )
318  END DO
319 *
320  CALL cgemm( 'No transpose', 'Conjugate transpose', p, m-q, m-q,
321  $ one, xf(1,q+1), ldx, v2t, ldv2t, zero, work, ldx )
322 *
323  CALL cgemm( '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  $ cmplx( sin(theta(r-i+1)), 0.0e0 )
333  END DO
334 *
335  CALL cgemm( 'No transpose', 'Conjugate transpose', m-p, q, q, one,
336  $ xf(p+1,1), ldx, v1t, ldv1t, zero, work, ldx )
337 *
338  CALL cgemm( '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  $ cmplx( sin(theta(r-i+1)), 0.0e0 )
348  END DO
349 *
350  CALL cgemm( '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 cgemm( '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  $ cmplx( cos(theta(i)), 0.0e0 )
363  END DO
364 *
365 * Compute norm( U1'*X11*V1 - D11 ) / ( MAX(1,P,Q)*EPS2 ) .
366 *
367  resid = clange( '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 = clange( '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 = clange( '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 = clange( '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 claset( 'Full', p, p, zero, one, work, ldu1 )
388  CALL cherk( '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 = clanhe( '1', 'Upper', p, work, ldu1, rwork )
394  result( 5 ) = ( resid / REAL(MAX(1,P)) ) / ulp
395 *
396 * Compute I - U2'*U2
397 *
398  CALL claset( 'Full', m-p, m-p, zero, one, work, ldu2 )
399  CALL cherk( '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 = clanhe( '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 claset( 'Full', q, q, zero, one, work, ldv1t )
410  CALL cherk( '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 = clanhe( '1', 'Upper', q, work, ldv1t, rwork )
416  result( 7 ) = ( resid / REAL(MAX(1,Q)) ) / ulp
417 *
418 * Compute I - V2T*V2T'
419 *
420  CALL claset( 'Full', m-q, m-q, zero, one, work, ldv2t )
421  CALL cherk( '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 = clanhe( '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 claset( 'Full', q, q, zero, one, work, ldx )
446  CALL cherk( 'Upper', 'Conjugate transpose', q, m, -realone,
447  $ x, ldx, realone, work, ldx )
448  IF (m.GT.0) THEN
449  eps2 = max( ulp,
450  $ clange( '1', q, q, work, ldx, rwork ) / REAL( 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 clacpy( 'Full', m, q, x, ldx, xf, ldx )
459 *
460 * Compute the CSD
461 *
462  CALL cuncsd2by1( '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 cgemm( 'No transpose', 'Conjugate transpose', p, q, q, one,
469  $ x, ldx, v1t, ldv1t, zero, work, ldx )
470 *
471  CALL cgemm( '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) - cmplx( cos(theta(i)),
480  $ 0.0e0 )
481  END DO
482 *
483  CALL cgemm( 'No transpose', 'Conjugate transpose', m-p, q, q, one,
484  $ x(p+1,1), ldx, v1t, ldv1t, zero, work, ldx )
485 *
486  CALL cgemm( '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  $ cmplx( sin(theta(r-i+1)), 0.0e0 )
496  END DO
497 *
498 * Compute norm( U1'*X11*V1 - D11 ) / ( MAX(1,P,Q)*EPS2 ) .
499 *
500  resid = clange( '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 = clange( '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 claset( 'Full', p, p, zero, one, work, ldu1 )
511  CALL cherk( 'Upper', 'Conjugate transpose', p, p, -realone,
512  $ u1, ldu1, realone, work, ldu1 )
513 *
514 * Compute norm( I - U1'*U1 ) / ( MAX(1,P) * ULP ) .
515 *
516  resid = clanhe( '1', 'Upper', p, work, ldu1, rwork )
517  result( 12 ) = ( resid / REAL(MAX(1,P)) ) / ulp
518 *
519 * Compute I - U2'*U2
520 *
521  CALL claset( 'Full', m-p, m-p, zero, one, work, ldu2 )
522  CALL cherk( '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 = clanhe( '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 claset( 'Full', q, q, zero, one, work, ldv1t )
533  CALL cherk( '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 = clanhe( '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 CCSDTS
558 *
559  END
cherk
subroutine cherk(UPLO, TRANS, N, K, ALPHA, A, LDA, BETA, C, LDC)
CHERK
Definition: cherk.f:175
claset
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:108
ccsdts
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:231
cuncsd2by1
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:256
clacpy
subroutine clacpy(UPLO, M, N, A, LDA, B, LDB)
CLACPY copies all or part of one two-dimensional array to another.
Definition: clacpy.f:105
cgemm
subroutine cgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
CGEMM
Definition: cgemm.f:189
cuncsd
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:322