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 ]

M

          M is INTEGER


          The number of rows of the matrix X.  M >= 0.

P

          P is INTEGER


          The number of rows of the matrix X11.  P >= 0.

Q

          Q is INTEGER


          The number of columns of the matrix X11.  Q >= 0.

X

          X is COMPLEX array, dimension (LDX,M)


          The M-by-M matrix X.

XF

          XF is COMPLEX array, dimension (LDX,M)


          Details of the CSD of X, as returned by CUNCSD;


          see CUNCSD for further details.

LDX

          LDX is INTEGER


          The leading dimension of the arrays X and XF.


          LDX >= max( 1,M ).

U1

          U1 is COMPLEX array, dimension(LDU1,P)


          The P-by-P unitary matrix U1.

LDU1

          LDU1 is INTEGER


          The leading dimension of the array U1. LDU >= max(1,P).

U2

          U2 is COMPLEX array, dimension(LDU2,M-P)


          The (M-P)-by-(M-P) unitary matrix U2.

LDU2

          LDU2 is INTEGER


          The leading dimension of the array U2. LDU >= max(1,M-P).

V1T

          V1T is COMPLEX array, dimension(LDV1T,Q)


          The Q-by-Q unitary matrix V1T.

LDV1T

          LDV1T is INTEGER


          The leading dimension of the array V1T. LDV1T >=


          max(1,Q).

V2T

          V2T is COMPLEX array, dimension(LDV2T,M-Q)


          The (M-Q)-by-(M-Q) unitary matrix V2T.

LDV2T

          LDV2T is INTEGER


          The leading dimension of the array V2T. LDV2T >=


          max(1,M-Q).

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.

IWORK

          IWORK is INTEGER array, dimension (M)

WORK

          WORK is COMPLEX array, dimension (LWORK)

LWORK

          LWORK is INTEGER


          The dimension of the array WORK

RWORK

          RWORK is REAL array

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 ). )

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.