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