!>
!> DCSDTS tests DORCSD, which, given an M-by-M partitioned orthogonal
!> 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 DORCSD2BY1, 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 DOUBLE PRECISION array, dimension (LDX,M)
!>          The M-by-M matrix X.
!> 

XF

!>          XF is DOUBLE PRECISION array, dimension (LDX,M)
!>          Details of the CSD of X, as returned by DORCSD;
!>          see DORCSD for further details.
!> 

LDX

!>          LDX is INTEGER
!>          The leading dimension of the arrays X and XF.
!>          LDX >= max( 1,M ).
!> 

U1

!>          U1 is DOUBLE PRECISION array, dimension(LDU1,P)
!>          The P-by-P orthogonal matrix U1.
!> 

LDU1

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

U2

!>          U2 is DOUBLE PRECISION array, dimension(LDU2,M-P)
!>          The (M-P)-by-(M-P) orthogonal matrix U2.
!> 

LDU2

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

V1T

!>          V1T is DOUBLE PRECISION array, dimension(LDV1T,Q)
!>          The Q-by-Q orthogonal matrix V1T.
!> 

LDV1T

!>          LDV1T is INTEGER
!>          The leading dimension of the array V1T. LDV1T >=
!>          max(1,Q).
!> 

V2T

!>          V2T is DOUBLE PRECISION array, dimension(LDV2T,M-Q)
!>          The (M-Q)-by-(M-Q) orthogonal 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 DORCSD for
!>          details.
!> 

IWORK

!>          IWORK is INTEGER array, dimension (M)
!> 

WORK

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