ZCSDTS tests ZUNCSD, 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*16 array, dimension (LDX,M)


          The M-by-M matrix X.

XF

          XF is COMPLEX*16 array, dimension (LDX,M)


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


          see ZUNCSD for further details.

LDX

          LDX is INTEGER


          The leading dimension of the arrays X and XF.


          LDX >= max( 1,M ).

U1

          U1 is COMPLEX*16 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*16 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*16 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*16 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 DOUBLE PRECISION 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 ZUNCSD for


          details.

IWORK

          IWORK is INTEGER array, dimension (M)

WORK

          WORK is COMPLEX*16 array, dimension (LWORK)

LWORK

          LWORK is INTEGER


          The dimension of the array WORK

RWORK

          RWORK is DOUBLE PRECISION array

RESULT

          RESULT is DOUBLE PRECISION 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.