LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
subroutine ccsdts ( integer  M,
integer  P,
integer  Q,
complex, dimension( ldx, * )  X,
complex, dimension( ldx, * )  XF,
integer  LDX,
complex, dimension( ldu1, * )  U1,
integer  LDU1,
complex, dimension( ldu2, * )  U2,
integer  LDU2,
complex, dimension( ldv1t, * )  V1T,
integer  LDV1T,
complex, dimension( ldv2t, * )  V2T,
integer  LDV2T,
real, dimension( * )  THETA,
integer, dimension( * )  IWORK,
complex, dimension( lwork )  WORK,
integer  LWORK,
real, dimension( * )  RWORK,
real, dimension( 15 )  RESULT 
)

CCSDTS

Purpose: 
 CCSDTS tests CUNCSD, which, given an M-by-M partitioned unitary
 matrix X,
              Q  M-Q
       X = [ X11 X12 ] P   ,
           [ X21 X22 ] M-P

 computes the CSD

       [ U1    ]**T * [ X11 X12 ] * [ V1    ]
       [    U2 ]      [ X21 X22 ]   [    V2 ]

                             [  I  0  0 |  0  0  0 ]
                             [  0  C  0 |  0 -S  0 ]
                             [  0  0  0 |  0  0 -I ]
                           = [---------------------] = [ D11 D12 ] .
                             [  0  0  0 |  I  0  0 ]   [ D21 D22 ]
                             [  0  S  0 |  0  C  0 ]
                             [  0  0  I |  0  0  0 ]

 and also SORCSD2BY1, which, given
          Q
       [ X11 ] P   ,
       [ X21 ] M-P

 computes the 2-by-1 CSD

                                     [  I  0  0 ]
                                     [  0  C  0 ]
                                     [  0  0  0 ]
       [ U1    ]**T * [ X11 ] * V1 = [----------] = [ D11 ] ,
       [    U2 ]      [ X21 ]        [  0  0  0 ]   [ D21 ]
                                     [  0  S  0 ]
                                     [  0  0  I ]
Parameters
[in]M
          M is INTEGER
          The number of rows of the matrix X.  M >= 0.
[in]P
          P is INTEGER
          The number of rows of the matrix X11.  P >= 0.
[in]Q
          Q is INTEGER
          The number of columns of the matrix X11.  Q >= 0.
[in]X
          X is COMPLEX array, dimension (LDX,M)
          The M-by-M matrix X.
[out]XF
          XF is COMPLEX array, dimension (LDX,M)
          Details of the CSD of X, as returned by CUNCSD;
          see CUNCSD for further details.
[in]LDX
          LDX is INTEGER
          The leading dimension of the arrays X and XF.
          LDX >= max( 1,M ).
[out]U1
          U1 is COMPLEX array, dimension(LDU1,P)
          The P-by-P unitary matrix U1.
[in]LDU1
          LDU1 is INTEGER
          The leading dimension of the array U1. LDU >= max(1,P).
[out]U2
          U2 is COMPLEX array, dimension(LDU2,M-P)
          The (M-P)-by-(M-P) unitary matrix U2.
[in]LDU2
          LDU2 is INTEGER
          The leading dimension of the array U2. LDU >= max(1,M-P).
[out]V1T
          V1T is COMPLEX array, dimension(LDV1T,Q)
          The Q-by-Q unitary matrix V1T.
[in]LDV1T
          LDV1T is INTEGER
          The leading dimension of the array V1T. LDV1T >=
          max(1,Q).
[out]V2T
          V2T is COMPLEX array, dimension(LDV2T,M-Q)
          The (M-Q)-by-(M-Q) unitary matrix V2T.
[in]LDV2T
          LDV2T is INTEGER
          The leading dimension of the array V2T. LDV2T >=
          max(1,M-Q).
[out]THETA
          THETA is REAL array, dimension MIN(P,M-P,Q,M-Q)
          The CS values of X; the essentially diagonal matrices C and
          S are constructed from THETA; see subroutine CUNCSD for
          details.
[out]IWORK
          IWORK is INTEGER array, dimension (M)
[out]WORK
          WORK is COMPLEX array, dimension (LWORK)
[in]LWORK
          LWORK is INTEGER
          The dimension of the array WORK
[out]RWORK
          RWORK is REAL array
[out]RESULT
          RESULT is REAL array, dimension (15)
          The test ratios:
          First, the 2-by-2 CSD:
          RESULT(1) = norm( U1'*X11*V1 - D11 ) / ( MAX(1,P,Q)*EPS2 )
          RESULT(2) = norm( U1'*X12*V2 - D12 ) / ( MAX(1,P,M-Q)*EPS2 )
          RESULT(3) = norm( U2'*X21*V1 - D21 ) / ( MAX(1,M-P,Q)*EPS2 )
          RESULT(4) = norm( U2'*X22*V2 - D22 ) / ( MAX(1,M-P,M-Q)*EPS2 )
          RESULT(5) = norm( I - U1'*U1 ) / ( MAX(1,P)*ULP )
          RESULT(6) = norm( I - U2'*U2 ) / ( MAX(1,M-P)*ULP )
          RESULT(7) = norm( I - V1T'*V1T ) / ( MAX(1,Q)*ULP )
          RESULT(8) = norm( I - V2T'*V2T ) / ( MAX(1,M-Q)*ULP )
          RESULT(9) = 0        if THETA is in increasing order and
                               all angles are in [0,pi/2];
                    = ULPINV   otherwise.
          Then, the 2-by-1 CSD:
          RESULT(10) = norm( U1'*X11*V1 - D11 ) / ( MAX(1,P,Q)*EPS2 )
          RESULT(11) = norm( U2'*X21*V1 - D21 ) / ( MAX(1,M-P,Q)*EPS2 )
          RESULT(12) = norm( I - U1'*U1 ) / ( MAX(1,P)*ULP )
          RESULT(13) = norm( I - U2'*U2 ) / ( MAX(1,M-P)*ULP )
          RESULT(14) = norm( I - V1T'*V1T ) / ( MAX(1,Q)*ULP )
          RESULT(15) = 0        if THETA is in increasing order and
                                all angles are in [0,pi/2];
                     = ULPINV   otherwise.
          ( EPS2 = MAX( norm( I - X'*X ) / M, ULP ). )
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
November 2015

Definition at line 231 of file ccsdts.f.

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 *
clanhe
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:126
cherk
subroutine cherk(UPLO, TRANS, N, K, ALPHA, A, LDA, BETA, C, LDC)
CHERK
Definition: cherk.f:175
clange
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:117
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
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
slamch
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:69
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

Here is the call graph for this function:

Here is the caller graph for this function: