Added CS decomposition test files to TESTING

9 files changed, 2444 insertions, 0 deletions

diff --git a/TESTING/EIG/cckcsd.f b/TESTING/EIG/cckcsd.f
new file mode 100644
index 00000000..47315354
--- /dev/null
+++ b/TESTING/EIG/cckcsd.f
@@ -0,0 +1,296 @@
+      SUBROUTINE CCKCSD( NM, MVAL, PVAL, QVAL, NMATS, ISEED, THRESH,
+     $                   MMAX, X, XF, U1, U2, V1T, V2T, THETA, IWORK,
+     $                   WORK, RWORK, NIN, NOUT, INFO )
+      IMPLICIT NONE
+*
+*     Originally CCKGSV
+*  -- LAPACK test routine (version 3.1) --
+*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
+*     November 2006
+*
+*     Adapted to CCKCSD by Brian Sutton
+*     July 2010
+*
+*     .. Scalar Arguments ..
+      INTEGER            INFO, NIN, NM, NMATS, MMAX, NOUT
+      REAL               THRESH
+*     ..
+*     .. Array Arguments ..
+      INTEGER            ISEED( 4 ), IWORK( * ), MVAL( * ), PVAL( * ),
+     $                   QVAL( * )
+      REAL               RWORK( * ), THETA( * )
+      COMPLEX            U1( * ), U2( * ), V1T( * ), V2T( * ),
+     $                   WORK( * ), X( * ), XF( * )
+*     ..
+*
+*  Purpose
+*  =======
+*
+*  CCKCSD tests CUNCSD:
+*         the CSD for an M-by-M unitary matrix X partitioned as
+*         [ X11 X12; X21 X22 ]. X11 is P-by-Q.
+*
+*  Arguments
+*  =========
+*
+*  NM      (input) INTEGER
+*          The number of values of M contained in the vector MVAL.
+*
+*  MVAL    (input) INTEGER array, dimension (NM)
+*          The values of the matrix row dimension M.
+*
+*  PVAL    (input) INTEGER array, dimension (NM)
+*          The values of the matrix row dimension P.
+*
+*  QVAL    (input) INTEGER array, dimension (NM)
+*          The values of the matrix column dimension Q.
+*
+*  NMATS   (input) INTEGER
+*          The number of matrix types to be tested for each combination
+*          of matrix dimensions.  If NMATS >= NTYPES (the maximum
+*          number of matrix types), then all the different types are
+*          generated for testing.  If NMATS < NTYPES, another input line
+*          is read to get the numbers of the matrix types to be used.
+*
+*  ISEED   (input/output) INTEGER array, dimension (4)
+*          On entry, the seed of the random number generator.  The array
+*          elements should be between 0 and 4095, otherwise they will be
+*          reduced mod 4096, and ISEED(4) must be odd.
+*          On exit, the next seed in the random number sequence after
+*          all the test matrices have been generated.
+*
+*  THRESH  (input) REAL
+*          The threshold value for the test ratios.  A result is
+*          included in the output file if RESULT >= THRESH.  To have
+*          every test ratio printed, use THRESH = 0.
+*
+*  MMAX    (input) INTEGER
+*          The maximum value permitted for M, used in dimensioning the
+*          work arrays.
+*
+*  X       (workspace) COMPLEX array, dimension (MMAX*MMAX)
+*
+*  XF      (workspace) COMPLEX array, dimension (MMAX*MMAX)
+*
+*  U1      (workspace) COMPLEX array, dimension (MMAX*MMAX)
+*
+*  U2      (workspace) COMPLEX array, dimension (MMAX*MMAX)
+*
+*  V1T     (workspace) COMPLEX array, dimension (MMAX*MMAX)
+*
+*  V2T     (workspace) COMPLEX array, dimension (MMAX*MMAX)
+*
+*  THETA   (workspace) REAL array, dimension (MMAX)
+*
+*  IWORK   (workspace) INTEGER array, dimension (MMAX)
+*
+*  WORK    (workspace) COMPLEX array
+*
+*  RWORK   (workspace) REAL array
+*
+*  NIN     (input) INTEGER
+*          The unit number for input.
+*
+*  NOUT    (input) INTEGER
+*          The unit number for output.
+*
+*  INFO    (output) INTEGER
+*          = 0 :  successful exit
+*          > 0 :  If CLAROR returns an error code, the absolute value
+*                 of it is returned.
+*
+*  =====================================================================
+*
+*     .. Parameters ..
+      INTEGER            NTESTS
+      PARAMETER          ( NTESTS = 9 )
+      INTEGER            NTYPES
+      PARAMETER          ( NTYPES = 4 )
+      REAL               GAPDIGIT, ORTH, PIOVER2, TEN
+      PARAMETER          ( GAPDIGIT = 10.0E0, ORTH = 1.0E-4,
+     $                     PIOVER2 = 1.57079632679489662E0,
+     $                     TEN = 10.0D0 )
+      COMPLEX            ZERO, ONE
+      PARAMETER          ( ZERO = (0.0E0,0.0E0), ONE = (1.0E0,0.0E0) )
+*     ..
+*     .. Local Scalars ..
+      LOGICAL            FIRSTT
+      CHARACTER*3        PATH
+      INTEGER            I, IINFO, IM, IMAT, J, LDU1, LDU2, LDV1T,
+     $                   LDV2T, LDX, LWORK, M, NFAIL, NRUN, NT, P, Q, R
+*     ..
+*     .. Local Arrays ..
+      LOGICAL            DOTYPE( NTYPES )
+      DOUBLE PRECISION   RESULT( NTESTS )
+*     ..
+*     .. External Subroutines ..
+      EXTERNAL           ALAHDG, ALAREQ, ALASUM, CCSDTS, CLACSG, CLAROR,
+     $                   CLASET
+*     ..
+*     .. Intrinsic Functions ..
+      INTRINSIC          ABS, COS, MIN, SIN
+*     ..
+*     .. External Functions ..
+      REAL               CLANGE, SLARND
+      EXTERNAL           CLANGE, SLARND
+*     ..
+*     .. Executable Statements ..
+*
+*     Initialize constants and the random number seed.
+*
+      PATH( 1: 3 ) = 'CSD'
+      INFO = 0
+      NRUN = 0
+      NFAIL = 0
+      FIRSTT = .TRUE.
+      CALL ALAREQ( PATH, NMATS, DOTYPE, NTYPES, NIN, NOUT )
+      LDX = MMAX
+      LDU1 = MMAX
+      LDU2 = MMAX
+      LDV1T = MMAX
+      LDV2T = MMAX
+      LWORK = MMAX*MMAX
+*
+*     Do for each value of M in MVAL.
+*
+      DO 30 IM = 1, NM
+         M = MVAL( IM )
+         P = PVAL( IM )
+         Q = QVAL( IM )
+*
+         DO 20 IMAT = 1, NTYPES
+*
+*           Do the tests only if DOTYPE( IMAT ) is true.
+*
+            IF( .NOT.DOTYPE( IMAT ) )
+     $         GO TO 20
+*
+*           Generate X
+*
+            IF( IMAT.EQ.1 ) THEN
+               CALL CLAROR( 'L', 'I', M, M, X, LDX, ISEED, WORK, IINFO )
+               IF( M .NE. 0 .AND. IINFO .NE. 0 ) THEN
+                  WRITE( NOUT, FMT = 9999 ) M, IINFO
+                  INFO = ABS( IINFO )
+                  GO TO 20
+               END IF
+            ELSE IF( IMAT.EQ.2 ) THEN
+               R = MIN( P, M-P, Q, M-Q )
+               DO I = 1, R
+                  THETA(I) = PIOVER2 * SLARND( 1, ISEED )
+               END DO
+               CALL CLACSG( M, P, Q, THETA, ISEED, X, LDX, WORK )
+               DO I = 1, M
+                  DO J = 1, M
+                     X(I+(J-1)*LDX) = X(I+(J-1)*LDX) +
+     $                                ORTH*SLARND(2,ISEED)
+                  END DO
+               END DO
+            ELSE IF( IMAT.EQ.3 ) THEN
+               R = MIN( P, M-P, Q, M-Q )
+               DO I = 1, R+1
+                  THETA(I) = TEN**(-SLARND(1,ISEED)*GAPDIGIT)
+               END DO
+               DO I = 2, R+1
+                  THETA(I) = THETA(I-1) + THETA(I)
+               END DO
+               DO I = 1, R
+                  THETA(I) = PIOVER2 * THETA(I) / THETA(R+1)
+               END DO
+               CALL CLACSG( M, P, Q, THETA, ISEED, X, LDX, WORK )
+            ELSE
+            END IF
+*
+            NT = 9
+*
+            CALL CCSDTS( M, P, Q, X, XF, LDX, U1, LDU1, U2, LDU2, V1T,
+     $                   LDV1T, V2T, LDV2T, THETA, IWORK, WORK, LWORK,
+     $                   RWORK, RESULT )
+*
+*           Print information about the tests that did not
+*           pass the threshold.
+*
+            DO 10 I = 1, NT
+               IF( RESULT( I ).GE.THRESH ) THEN
+                  IF( NFAIL.EQ.0 .AND. FIRSTT ) THEN
+                     FIRSTT = .FALSE.
+                     CALL ALAHDG( NOUT, PATH )
+                  END IF
+                  WRITE( NOUT, FMT = 9998 )M, P, Q, IMAT, I,
+     $               RESULT( I )
+                  NFAIL = NFAIL + 1
+               END IF
+   10       CONTINUE
+            NRUN = NRUN + NT
+   20    CONTINUE
+   30 CONTINUE
+*
+*     Print a summary of the results.
+*
+      CALL ALASUM( PATH, NOUT, NFAIL, NRUN, 0 )
+*
+ 9999 FORMAT( ' CLAROR in CCKCSD: M = ', I5, ', INFO = ', I15 )
+ 9998 FORMAT( ' M=', I4, ' P=', I4, ', Q=', I4, ', type ', I2,
+     $      ', test ', I2, ', ratio=', G13.6 )
+      RETURN
+*
+*     End of CCKCSD
+*
+      END
+*
+*
+*
+      SUBROUTINE CLACSG( M, P, Q, THETA, ISEED, X, LDX, WORK )
+      IMPLICIT NONE
+*
+      INTEGER            LDX, M, P, Q
+      INTEGER            ISEED( 4 )
+      REAL               THETA( * )
+      COMPLEX            WORK( * ), X( LDX, * )
+*
+      COMPLEX            ONE, ZERO
+      PARAMETER          ( ONE = (1.0E0,0.0E0), ZERO = (0.0E0,0.0E0) )
+*
+      INTEGER            I, INFO, R
+*
+      R = MIN( P, M-P, Q, M-Q )
+*
+      CALL CLASET( 'Full', M, M, ZERO, ZERO, X, LDX )
+*
+      DO I = 1, MIN(P,Q)-R
+         X(I,I) = ONE
+      END DO
+      DO I = 1, R
+         X(MIN(P,Q)-R+I,MIN(P,Q)-R+I) = COMPLEX( COS(THETA(I)), 0.0E0 )
+      END DO
+      DO I = 1, MIN(P,M-Q)-R
+         X(P-I+1,M-I+1) = -ONE
+      END DO
+      DO I = 1, R
+         X(P-(MIN(P,M-Q)-R)+1-I,M-(MIN(P,M-Q)-R)+1-I) =
+     $      COMPLEX( -SIN(THETA(R-I+1)), 0.0E0 )
+      END DO
+      DO I = 1, MIN(M-P,Q)-R
+         X(M-I+1,Q-I+1) = ONE
+      END DO
+      DO I = 1, R
+         X(M-(MIN(M-P,Q)-R)+1-I,Q-(MIN(M-P,Q)-R)+1-I) =
+     $      COMPLEX( SIN(THETA(R-I+1)), 0.0E0 )
+      END DO
+      DO I = 1, MIN(M-P,M-Q)-R
+         X(P+I,Q+I) = ONE
+      END DO
+      DO I = 1, R
+         X(P+(MIN(M-P,M-Q)-R)+I,Q+(MIN(M-P,M-Q)-R)+I) =
+     $      COMPLEX( COS(THETA(I)), 0.0E0 )
+      END DO
+      CALL CLAROR( 'Left', 'No init', P, M, X, LDX, ISEED, WORK, INFO )
+      CALL CLAROR( 'Left', 'No init', M-P, M, X(P+1,1), LDX,
+     $             ISEED, WORK, INFO )
+      CALL CLAROR( 'Right', 'No init', M, Q, X, LDX, ISEED,
+     $             WORK, INFO )
+      CALL CLAROR( 'Right', 'No init', M, M-Q,
+     $             X(1,Q+1), LDX, ISEED, WORK, INFO )
+*
+      END
+
diff --git a/TESTING/EIG/ccsdts.f b/TESTING/EIG/ccsdts.f
new file mode 100644
index 00000000..e9be18b4
--- /dev/null
+++ b/TESTING/EIG/ccsdts.f
@@ -0,0 +1,313 @@
+      SUBROUTINE CCSDTS( M, P, Q, X, XF, LDX, U1, LDU1, U2, LDU2, V1T,
+     $                   LDV1T, V2T, LDV2T, THETA, IWORK, WORK, LWORK,
+     $                   RWORK, RESULT )
+      IMPLICIT NONE
+*
+*     Originally xGSVTS
+*  -- LAPACK test routine (version 3.1) --
+*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
+*     November 2006
+*
+*     Adapted to CCSDTS by Brian Sutton
+*     July 2010
+*
+*     .. Scalar Arguments ..
+      INTEGER            LDX, LDU1, LDU2, LDV1T, LDV2T, LWORK, M, P, Q
+*     ..
+*     .. Array Arguments ..
+      INTEGER            IWORK( * )
+      REAL               RESULT( 9 ), RWORK( * ), THETA( * )
+      COMPLEX            U1( LDU1, * ), U2( LDU2, * ), V1T( LDV1T, * ),
+     $                   V2T( LDV2T, * ), WORK( LWORK ), X( LDX, * ),
+     $                   XF( LDX, * )
+*     ..
+*
+*  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 ]
+*
+*  Arguments
+*  =========
+*
+*  M       (input) INTEGER
+*          The number of rows of the matrix X.  M >= 0.
+*
+*  P       (input) INTEGER
+*          The number of rows of the matrix X11.  P >= 0.
+*
+*  Q       (input) INTEGER
+*          The number of columns of the matrix X11.  Q >= 0.
+*
+*  X       (input) COMPLEX array, dimension (LDX,M)
+*          The M-by-M matrix X.
+*
+*  XF      (output) COMPLEX array, dimension (LDX,M)
+*          Details of the CSD of X, as returned by CUNCSD;
+*          see CUNCSD for further details.
+*
+*  LDX     (input) INTEGER
+*          The leading dimension of the arrays X and XF.
+*          LDX >= max( 1,M ).
+*
+*  U1      (output) COMPLEX array, dimension(LDU1,P)
+*          The P-by-P unitary matrix U1.
+*
+*  LDU1    (input) INTEGER
+*          The leading dimension of the array U1. LDU >= max(1,P).
+*
+*  U2      (output) COMPLEX array, dimension(LDU2,M-P)
+*          The (M-P)-by-(M-P) unitary matrix U2.
+*
+*  LDU2    (input) INTEGER
+*          The leading dimension of the array U2. LDU >= max(1,M-P).
+*
+*  V1T     (output) COMPLEX array, dimension(LDV1T,Q)
+*          The Q-by-Q unitary matrix V1T.
+*
+*  LDV1T   (input) INTEGER
+*          The leading dimension of the array V1T. LDV1T >=
+*          max(1,Q).
+*
+*  V2T     (output) COMPLEX array, dimension(LDV2T,M-Q)
+*          The (M-Q)-by-(M-Q) unitary matrix V2T.
+*
+*  LDV2T   (input) INTEGER
+*          The leading dimension of the array V2T. LDV2T >=
+*          max(1,M-Q).
+*
+*  THETA   (output) 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 CUNCSD for
+*          details.
+*
+*  IWORK   (workspace) INTEGER array, dimension (M)
+*
+*  WORK    (workspace) COMPLEX array, dimension (LWORK)
+*
+*  LWORK   (input) INTEGER
+*          The dimension of the array WORK
+*
+*  RWORK   (workspace) REAL array
+*
+*  RESULT  (output) REAL array, dimension (9)
+*          The test ratios:
+*          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.
+*          ( EPS2 = MAX( norm( I - X'*X ) / M, ULP ). )
+*
+*  =====================================================================
+*
+*     .. Parameters ..
+      REAL               PIOVER2, REALONE, REALZERO
+      PARAMETER          ( PIOVER2 = 1.57079632679489662E0,
+     $                     REALONE = 1.0E0, REALZERO = 0.0E0 )
+      COMPLEX            ZERO, ONE
+      PARAMETER          ( ZERO = (0.0E0,0.0E0), ONE = (1.0E0,0.0E0) )
+*     ..
+*     .. Local Scalars ..
+      INTEGER            I, INFO, R
+      REAL               EPS2, RESID, ULP, ULPINV
+*     ..
+*     .. External Functions ..
+      REAL               SLAMCH, CLANGE, CLANHE
+      EXTERNAL           SLAMCH, CLANGE, CLANHE
+*     ..
+*     .. External Subroutines ..
+      EXTERNAL           CGEMM, CLACPY, CLASET, CUNCSD, CHERK
+*     ..
+*     .. Intrinsic Functions ..
+      INTRINSIC          REAL, MAX, MIN
+*     ..
+*     .. Executable Statements ..
+*
+      ULP = SLAMCH( 'Precision' )
+      ULPINV = REALONE / ULP
+      CALL CLASET( 'Full', M, M, ZERO, ONE, WORK, LDX )
+      CALL CHERK( 'Upper', 'Conjugate transpose', M, M, -ONE, X, LDX,
+     $            ONE, WORK, LDX )
+      EPS2 = MAX( ULP, 
+     $            CLANGE( '1', M, M, WORK, LDX, RWORK ) / REAL( M ) )
+      R = MIN( P, M-P, Q, M-Q )
+*
+*     Copy the matrix X to the array XF.
+*
+      CALL CLACPY( 'Full', M, M, X, LDX, XF, LDX )
+*
+*     Compute the CSD
+*
+      CALL CUNCSD( 'Y', 'Y', 'Y', 'Y', 'N', 'D', M, P, Q, XF(1,1), LDX,
+     $             XF(1,Q+1), LDX, XF(P+1,1), LDX, XF(P+1,Q+1), LDX,
+     $             THETA, U1, LDU1, U2, LDU2, V1T, LDV1T, V2T, LDV2T,
+     $             WORK, LWORK, RWORK, 17*(R+2), IWORK, INFO )
+*
+*     Compute X := diag(U1,U2)'*X*diag(V1,V2) - [D11 D12; D21 D22]
+*
+      CALL CGEMM( 'No transpose', 'Conjugate transpose', P, Q, Q, ONE,
+     $            X, LDX, V1T, LDV1T, ZERO, WORK, LDX )
+*
+      CALL CGEMM( 'Conjugate transpose', 'No transpose', P, Q, P, ONE,
+     $            U1, LDU1, WORK, LDX, ZERO, X, LDX )
+*
+      DO I = 1, MIN(P,Q)-R
+         X(I,I) = X(I,I) - ONE
+      END DO
+      DO I = 1, R
+         X(MIN(P,Q)-R+I,MIN(P,Q)-R+I) =
+     $           X(MIN(P,Q)-R+I,MIN(P,Q)-R+I) - COMPLEX( COS(THETA(I)),
+     $              0.0E0 )
+      END DO
+*
+      CALL CGEMM( 'No transpose', 'Conjugate transpose', P, M-Q, M-Q,
+     $            ONE, X(1,Q+1), LDX, V2T, LDV2T, ZERO, WORK, LDX )
+*
+      CALL CGEMM( 'Conjugate transpose', 'No transpose', P, M-Q, P,
+     $            ONE, U1, LDU1, WORK, LDX, ZERO, X(1,Q+1), LDX )
+*
+      DO I = 1, MIN(P,M-Q)-R
+         X(P-I+1,M-I+1) = X(P-I+1,M-I+1) + ONE
+      END DO
+      DO I = 1, R
+         X(P-(MIN(P,M-Q)-R)+1-I,M-(MIN(P,M-Q)-R)+1-I) =
+     $      X(P-(MIN(P,M-Q)-R)+1-I,M-(MIN(P,M-Q)-R)+1-I) +
+     $      COMPLEX( SIN(THETA(R-I+1)), 0.0E0 )
+      END DO
+*
+      CALL CGEMM( 'No transpose', 'Conjugate transpose', M-P, Q, Q, ONE,
+     $            X(P+1,1), LDX, V1T, LDV1T, ZERO, WORK, LDX )
+*
+      CALL CGEMM( 'Conjugate transpose', 'No transpose', M-P, Q, M-P,
+     $            ONE, U2, LDU2, WORK, LDX, ZERO, X(P+1,1), LDX )
+*
+      DO I = 1, MIN(M-P,Q)-R
+         X(M-I+1,Q-I+1) = X(M-I+1,Q-I+1) - ONE
+      END DO
+      DO I = 1, R
+         X(M-(MIN(M-P,Q)-R)+1-I,Q-(MIN(M-P,Q)-R)+1-I) =
+     $             X(M-(MIN(M-P,Q)-R)+1-I,Q-(MIN(M-P,Q)-R)+1-I) -
+     $             COMPLEX( SIN(THETA(R-I+1)), 0.0E0 )
+      END DO
+*
+      CALL CGEMM( 'No transpose', 'Conjugate transpose', M-P, M-Q, M-Q,
+     $            ONE, X(P+1,Q+1), LDX, V2T, LDV2T, ZERO, WORK, LDX )
+*
+      CALL CGEMM( 'Conjugate transpose', 'No transpose', M-P, M-Q, M-P,
+     $            ONE, U2, LDU2, WORK, LDX, ZERO, X(P+1,Q+1), LDX )
+*
+      DO I = 1, MIN(M-P,M-Q)-R
+         X(P+I,Q+I) = X(P+I,Q+I) - ONE
+      END DO
+      DO I = 1, R
+         X(P+(MIN(M-P,M-Q)-R)+I,Q+(MIN(M-P,M-Q)-R)+I) =
+     $      X(P+(MIN(M-P,M-Q)-R)+I,Q+(MIN(M-P,M-Q)-R)+I) -
+     $      COMPLEX( COS(THETA(I)), 0.0E0 )
+      END DO
+*
+*     Compute norm( U1'*X11*V1 - D11 ) / ( MAX(1,P,Q)*EPS2 ) .
+*
+      RESID = CLANGE( '1', P, Q, X, LDX, RWORK )
+      RESULT( 1 ) = ( RESID / REAL(MAX(1,P,Q)) ) / EPS2
+*
+*     Compute norm( U1'*X12*V2 - D12 ) / ( MAX(1,P,M-Q)*EPS2 ) .
+*
+      RESID = CLANGE( '1', P, M-Q, X(1,Q+1), LDX, RWORK )
+      RESULT( 2 ) = ( RESID / REAL(MAX(1,P,M-Q)) ) / EPS2
+*
+*     Compute norm( U2'*X21*V1 - D21 ) / ( MAX(1,M-P,Q)*EPS2 ) .
+*
+      RESID = CLANGE( '1', M-P, Q, X(P+1,1), LDX, RWORK )
+      RESULT( 3 ) = ( RESID / REAL(MAX(1,M-P,Q)) ) / EPS2
+*
+*     Compute norm( U2'*X22*V2 - D22 ) / ( MAX(1,M-P,M-Q)*EPS2 ) .
+*
+      RESID = CLANGE( '1', M-P, M-Q, X(P+1,Q+1), LDX, RWORK )
+      RESULT( 4 ) = ( RESID / REAL(MAX(1,M-P,M-Q)) ) / EPS2
+*
+*     Compute I - U1'*U1
+*
+      CALL CLASET( 'Full', P, P, ZERO, ONE, WORK, LDU1 )
+      CALL CHERK( 'Upper', 'Conjugate transpose', P, P, -ONE, U1, LDU1,
+     $            ONE, WORK, LDU1 )
+*
+*     Compute norm( I - U'*U ) / ( MAX(1,P) * ULP ) .
+*
+      RESID = CLANHE( '1', 'Upper', P, WORK, LDU1, RWORK )
+      RESULT( 5 ) = ( RESID / REAL(MAX(1,P)) ) / ULP
+*
+*     Compute I - U2'*U2
+*
+      CALL CLASET( 'Full', M-P, M-P, ZERO, ONE, WORK, LDU2 )
+      CALL CHERK( 'Upper', 'Conjugate transpose', M-P, M-P, -ONE, U2,
+     $            LDU2, ONE, WORK, LDU2 )
+*
+*     Compute norm( I - U2'*U2 ) / ( MAX(1,M-P) * ULP ) .
+*
+      RESID = CLANHE( '1', 'Upper', M-P, WORK, LDU2, RWORK )
+      RESULT( 6 ) = ( RESID / REAL(MAX(1,M-P)) ) / ULP
+*
+*     Compute I - V1T*V1T'
+*
+      CALL CLASET( 'Full', Q, Q, ZERO, ONE, WORK, LDV1T )
+      CALL CHERK( 'Upper', 'No transpose', Q, Q, -ONE, V1T, LDV1T, ONE,
+     $            WORK, LDV1T )
+*
+*     Compute norm( I - V1T*V1T' ) / ( MAX(1,Q) * ULP ) .
+*
+      RESID = CLANHE( '1', 'Upper', Q, WORK, LDV1T, RWORK )
+      RESULT( 7 ) = ( RESID / REAL(MAX(1,Q)) ) / ULP
+*
+*     Compute I - V2T*V2T'
+*
+      CALL CLASET( 'Full', M-Q, M-Q, ZERO, ONE, WORK, LDV2T )
+      CALL CHERK( 'Upper', 'No transpose', M-Q, M-Q, -ONE, V2T, LDV2T,
+     $            ONE, WORK, LDV2T )
+*
+*     Compute norm( I - V2T*V2T' ) / ( MAX(1,M-Q) * ULP ) .
+*
+      RESID = CLANHE( '1', 'Upper', M-Q, WORK, LDV2T, RWORK )
+      RESULT( 8 ) = ( RESID / REAL(MAX(1,M-Q)) ) / ULP
+*
+*     Check sorting
+*
+      RESULT(9) = REALZERO
+      DO I = 1, R
+         IF( THETA(I).LT.REALZERO .OR. THETA(I).GT.PIOVER2 ) THEN
+            RESULT(9) = ULPINV
+         END IF
+         IF( I.GT.1 .AND. THETA(I).LT.THETA(I-1) ) THEN
+            RESULT(9) = ULPINV
+         END IF
+      END DO
+*
+      RETURN
+*      
+*     End of CCSDTS
+*
+      END
+
diff --git a/TESTING/EIG/dckcsd.f b/TESTING/EIG/dckcsd.f
new file mode 100644
index 00000000..ad3cacc6
--- /dev/null
+++ b/TESTING/EIG/dckcsd.f
@@ -0,0 +1,296 @@
+      SUBROUTINE DCKCSD( NM, MVAL, PVAL, QVAL, NMATS, ISEED, THRESH,
+     $                   MMAX, X, XF, U1, U2, V1T, V2T, THETA, IWORK,
+     $                   WORK, RWORK, NIN, NOUT, INFO )
+      IMPLICIT NONE
+*
+*     Originally DCKGSV
+*  -- LAPACK test routine (version 3.1) --
+*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
+*     November 2006
+*
+*     Adapted to DCKCSD by Brian Sutton
+*     July 2010
+*
+*     .. Scalar Arguments ..
+      INTEGER            INFO, NIN, NM, NMATS, MMAX, NOUT
+      DOUBLE PRECISION   THRESH
+*     ..
+*     .. Array Arguments ..
+      INTEGER            ISEED( 4 ), IWORK( * ), MVAL( * ), PVAL( * ),
+     $                   QVAL( * )
+      DOUBLE PRECISION   RWORK( * ), THETA( * )
+      DOUBLE PRECISION   U1( * ), U2( * ), V1T( * ), V2T( * ),
+     $                   WORK( * ), X( * ), XF( * )
+*     ..
+*
+*  Purpose
+*  =======
+*
+*  DCKCSD tests DORCSD:
+*         the CSD for an M-by-M orthogonal matrix X partitioned as
+*         [ X11 X12; X21 X22 ]. X11 is P-by-Q.
+*
+*  Arguments
+*  =========
+*
+*  NM      (input) INTEGER
+*          The number of values of M contained in the vector MVAL.
+*
+*  MVAL    (input) INTEGER array, dimension (NM)
+*          The values of the matrix row dimension M.
+*
+*  PVAL    (input) INTEGER array, dimension (NM)
+*          The values of the matrix row dimension P.
+*
+*  QVAL    (input) INTEGER array, dimension (NM)
+*          The values of the matrix column dimension Q.
+*
+*  NMATS   (input) INTEGER
+*          The number of matrix types to be tested for each combination
+*          of matrix dimensions.  If NMATS >= NTYPES (the maximum
+*          number of matrix types), then all the different types are
+*          generated for testing.  If NMATS < NTYPES, another input line
+*          is read to get the numbers of the matrix types to be used.
+*
+*  ISEED   (input/output) INTEGER array, dimension (4)
+*          On entry, the seed of the random number generator.  The array
+*          elements should be between 0 and 4095, otherwise they will be
+*          reduced mod 4096, and ISEED(4) must be odd.
+*          On exit, the next seed in the random number sequence after
+*          all the test matrices have been generated.
+*
+*  THRESH  (input) DOUBLE PRECISION
+*          The threshold value for the test ratios.  A result is
+*          included in the output file if RESULT >= THRESH.  To have
+*          every test ratio printed, use THRESH = 0.
+*
+*  MMAX    (input) INTEGER
+*          The maximum value permitted for M, used in dimensioning the
+*          work arrays.
+*
+*  X       (workspace) DOUBLE PRECISION array, dimension (MMAX*MMAX)
+*
+*  XF      (workspace) DOUBLE PRECISION array, dimension (MMAX*MMAX)
+*
+*  U1      (workspace) DOUBLE PRECISION array, dimension (MMAX*MMAX)
+*
+*  U2      (workspace) DOUBLE PRECISION array, dimension (MMAX*MMAX)
+*
+*  V1T     (workspace) DOUBLE PRECISION array, dimension (MMAX*MMAX)
+*
+*  V2T     (workspace) DOUBLE PRECISION array, dimension (MMAX*MMAX)
+*
+*  THETA   (workspace) DOUBLE PRECISION array, dimension (MMAX)
+*
+*  IWORK   (workspace) INTEGER array, dimension (MMAX)
+*
+*  WORK    (workspace) DOUBLE PRECISION array
+*
+*  RWORK   (workspace) DOUBLE PRECISION array
+*
+*  NIN     (input) INTEGER
+*          The unit number for input.
+*
+*  NOUT    (input) INTEGER
+*          The unit number for output.
+*
+*  INFO    (output) INTEGER
+*          = 0 :  successful exit
+*          > 0 :  If DLAROR returns an error code, the absolute value
+*                 of it is returned.
+*
+*  =====================================================================
+*
+*     .. Parameters ..
+      INTEGER            NTESTS
+      PARAMETER          ( NTESTS = 9 )
+      INTEGER            NTYPES
+      PARAMETER          ( NTYPES = 4 )
+      DOUBLE PRECISION   GAPDIGIT, ORTH, PIOVER2, TEN
+      PARAMETER          ( GAPDIGIT = 18.0D0, ORTH = 1.0D-12,
+     $                     PIOVER2 = 1.57079632679489662D0,
+     $                     TEN = 10.0D0 )
+      DOUBLE PRECISION   ZERO, ONE
+      PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0 )
+*     ..
+*     .. Local Scalars ..
+      LOGICAL            FIRSTT
+      CHARACTER*3        PATH
+      INTEGER            I, IINFO, IM, IMAT, J, LDU1, LDU2, LDV1T,
+     $                   LDV2T, LDX, LWORK, M, NFAIL, NRUN, NT, P, Q, R
+*     ..
+*     .. Local Arrays ..
+      LOGICAL            DOTYPE( NTYPES )
+      DOUBLE PRECISION   RESULT( NTESTS )
+*     ..
+*     .. External Subroutines ..
+      EXTERNAL           ALAHDG, ALAREQ, ALASUM, DCSDTS, DLACSG, DLAROR,
+     $                   DLASET
+*     ..
+*     .. Intrinsic Functions ..
+      INTRINSIC          ABS, COS, MIN, SIN
+*     ..
+*     .. External Functions ..
+      DOUBLE PRECISION   DLANGE, DLARND
+      EXTERNAL           DLANGE, DLARND
+*     ..
+*     .. Executable Statements ..
+*
+*     Initialize constants and the random number seed.
+*
+      PATH( 1: 3 ) = 'CSD'
+      INFO = 0
+      NRUN = 0
+      NFAIL = 0
+      FIRSTT = .TRUE.
+      CALL ALAREQ( PATH, NMATS, DOTYPE, NTYPES, NIN, NOUT )
+      LDX = MMAX
+      LDU1 = MMAX
+      LDU2 = MMAX
+      LDV1T = MMAX
+      LDV2T = MMAX
+      LWORK = MMAX*MMAX
+*
+*     Do for each value of M in MVAL.
+*
+      DO 30 IM = 1, NM
+         M = MVAL( IM )
+         P = PVAL( IM )
+         Q = QVAL( IM )
+*
+         DO 20 IMAT = 1, NTYPES
+*
+*           Do the tests only if DOTYPE( IMAT ) is true.
+*
+            IF( .NOT.DOTYPE( IMAT ) )
+     $         GO TO 20
+*
+*           Generate X
+*
+            IF( IMAT.EQ.1 ) THEN
+               CALL DLAROR( 'L', 'I', M, M, X, LDX, ISEED, WORK, IINFO )
+               IF( M .NE. 0 .AND. IINFO .NE. 0 ) THEN
+                  WRITE( NOUT, FMT = 9999 ) M, IINFO
+                  INFO = ABS( IINFO )
+                  GO TO 20
+               END IF
+            ELSE IF( IMAT.EQ.2 ) THEN
+               R = MIN( P, M-P, Q, M-Q )
+               DO I = 1, R
+                  THETA(I) = PIOVER2 * DLARND( 1, ISEED )
+               END DO
+               CALL DLACSG( M, P, Q, THETA, ISEED, X, LDX, WORK )
+               DO I = 1, M
+                  DO J = 1, M
+                     X(I+(J-1)*LDX) = X(I+(J-1)*LDX) +
+     $                                ORTH*DLARND(2,ISEED)
+                  END DO
+               END DO
+            ELSE IF( IMAT.EQ.3 ) THEN
+               R = MIN( P, M-P, Q, M-Q )
+               DO I = 1, R+1
+                  THETA(I) = TEN**(-DLARND(1,ISEED)*GAPDIGIT)
+               END DO
+               DO I = 2, R+1
+                  THETA(I) = THETA(I-1) + THETA(I)
+               END DO
+               DO I = 1, R
+                  THETA(I) = PIOVER2 * THETA(I) / THETA(R+1)
+               END DO
+               CALL DLACSG( M, P, Q, THETA, ISEED, X, LDX, WORK )
+            ELSE
+            END IF
+*
+            NT = 9
+*
+            CALL DCSDTS( M, P, Q, X, XF, LDX, U1, LDU1, U2, LDU2, V1T,
+     $                   LDV1T, V2T, LDV2T, THETA, IWORK, WORK, LWORK,
+     $                   RWORK, RESULT )
+*
+*           Print information about the tests that did not
+*           pass the threshold.
+*
+            DO 10 I = 1, NT
+               IF( RESULT( I ).GE.THRESH ) THEN
+                  IF( NFAIL.EQ.0 .AND. FIRSTT ) THEN
+                     FIRSTT = .FALSE.
+                     CALL ALAHDG( NOUT, PATH )
+                  END IF
+                  WRITE( NOUT, FMT = 9998 )M, P, Q, IMAT, I,
+     $               RESULT( I )
+                  NFAIL = NFAIL + 1
+               END IF
+   10       CONTINUE
+            NRUN = NRUN + NT
+   20    CONTINUE
+   30 CONTINUE
+*
+*     Print a summary of the results.
+*
+      CALL ALASUM( PATH, NOUT, NFAIL, NRUN, 0 )
+*
+ 9999 FORMAT( ' DLAROR in DCKCSD: M = ', I5, ', INFO = ', I15 )
+ 9998 FORMAT( ' M=', I4, ' P=', I4, ', Q=', I4, ', type ', I2,
+     $      ', test ', I2, ', ratio=', G13.6 )
+      RETURN
+*
+*     End of DCKCSD
+*
+      END
+*
+*
+*
+      SUBROUTINE DLACSG( M, P, Q, THETA, ISEED, X, LDX, WORK )
+      IMPLICIT NONE
+*
+      INTEGER            LDX, M, P, Q
+      INTEGER            ISEED( 4 )
+      DOUBLE PRECISION   THETA( * )
+      DOUBLE PRECISION   WORK( * ), X( LDX, * )
+*
+      DOUBLE PRECISION   ONE, ZERO
+      PARAMETER          ( ONE = 1.0D0, ZERO = 0.0D0 )
+*
+      INTEGER            I, INFO, R
+*
+      R = MIN( P, M-P, Q, M-Q )
+*
+      CALL DLASET( 'Full', M, M, ZERO, ZERO, X, LDX )
+*
+      DO I = 1, MIN(P,Q)-R
+         X(I,I) = ONE
+      END DO
+      DO I = 1, R
+         X(MIN(P,Q)-R+I,MIN(P,Q)-R+I) = COS(THETA(I))
+      END DO
+      DO I = 1, MIN(P,M-Q)-R
+         X(P-I+1,M-I+1) = -ONE
+      END DO
+      DO I = 1, R
+         X(P-(MIN(P,M-Q)-R)+1-I,M-(MIN(P,M-Q)-R)+1-I) =
+     $      -SIN(THETA(R-I+1))
+      END DO
+      DO I = 1, MIN(M-P,Q)-R
+         X(M-I+1,Q-I+1) = ONE
+      END DO
+      DO I = 1, R
+         X(M-(MIN(M-P,Q)-R)+1-I,Q-(MIN(M-P,Q)-R)+1-I) =
+     $      SIN(THETA(R-I+1))
+      END DO
+      DO I = 1, MIN(M-P,M-Q)-R
+         X(P+I,Q+I) = ONE
+      END DO
+      DO I = 1, R
+         X(P+(MIN(M-P,M-Q)-R)+I,Q+(MIN(M-P,M-Q)-R)+I) =
+     $      COS(THETA(I))
+      END DO
+      CALL DLAROR( 'Left', 'No init', P, M, X, LDX, ISEED, WORK, INFO )
+      CALL DLAROR( 'Left', 'No init', M-P, M, X(P+1,1), LDX,
+     $             ISEED, WORK, INFO )
+      CALL DLAROR( 'Right', 'No init', M, Q, X, LDX, ISEED,
+     $             WORK, INFO )
+      CALL DLAROR( 'Right', 'No init', M, M-Q,
+     $             X(1,Q+1), LDX, ISEED, WORK, INFO )
+*
+      END
+
diff --git a/TESTING/EIG/dcsdts.f b/TESTING/EIG/dcsdts.f
new file mode 100644
index 00000000..22ef4bce
--- /dev/null
+++ b/TESTING/EIG/dcsdts.f
@@ -0,0 +1,312 @@
+      SUBROUTINE DCSDTS( M, P, Q, X, XF, LDX, U1, LDU1, U2, LDU2, V1T,
+     $                   LDV1T, V2T, LDV2T, THETA, IWORK, WORK, LWORK,
+     $                   RWORK, RESULT )
+      IMPLICIT NONE
+*
+*     Originally xGSVTS
+*  -- LAPACK test routine (version 3.1) --
+*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
+*     November 2006
+*
+*     Adapted to DCSDTS by Brian Sutton
+*     July 2010
+*
+*     .. Scalar Arguments ..
+      INTEGER            LDX, LDU1, LDU2, LDV1T, LDV2T, LWORK, M, P, Q
+*     ..
+*     .. Array Arguments ..
+      INTEGER            IWORK( * )
+      DOUBLE PRECISION   RESULT( 9 ), RWORK( * ), THETA( * )
+      DOUBLE PRECISION   U1( LDU1, * ), U2( LDU2, * ), V1T( LDV1T, * ),
+     $                   V2T( LDV2T, * ), WORK( LWORK ), X( LDX, * ),
+     $                   XF( LDX, * )
+*     ..
+*
+*  Purpose
+*  =======
+*
+*  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 ]
+*
+*  Arguments
+*  =========
+*
+*  M       (input) INTEGER
+*          The number of rows of the matrix X.  M >= 0.
+*
+*  P       (input) INTEGER
+*          The number of rows of the matrix X11.  P >= 0.
+*
+*  Q       (input) INTEGER
+*          The number of columns of the matrix X11.  Q >= 0.
+*
+*  X       (input) DOUBLE PRECISION array, dimension (LDX,M)
+*          The M-by-M matrix X.
+*
+*  XF      (output) DOUBLE PRECISION array, dimension (LDX,M)
+*          Details of the CSD of X, as returned by DORCSD;
+*          see DORCSD for further details.
+*
+*  LDX     (input) INTEGER
+*          The leading dimension of the arrays X and XF.
+*          LDX >= max( 1,M ).
+*
+*  U1      (output) DOUBLE PRECISION array, dimension(LDU1,P)
+*          The P-by-P orthogonal matrix U1.
+*
+*  LDU1    (input) INTEGER
+*          The leading dimension of the array U1. LDU >= max(1,P).
+*
+*  U2      (output) DOUBLE PRECISION array, dimension(LDU2,M-P)
+*          The (M-P)-by-(M-P) orthogonal matrix U2.
+*
+*  LDU2    (input) INTEGER
+*          The leading dimension of the array U2. LDU >= max(1,M-P).
+*
+*  V1T     (output) DOUBLE PRECISION array, dimension(LDV1T,Q)
+*          The Q-by-Q orthogonal matrix V1T.
+*
+*  LDV1T   (input) INTEGER
+*          The leading dimension of the array V1T. LDV1T >=
+*          max(1,Q).
+*
+*  V2T     (output) DOUBLE PRECISION array, dimension(LDV2T,M-Q)
+*          The (M-Q)-by-(M-Q) orthogonal matrix V2T.
+*
+*  LDV2T   (input) INTEGER
+*          The leading dimension of the array V2T. LDV2T >=
+*          max(1,M-Q).
+*
+*  THETA   (output) 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   (workspace) INTEGER array, dimension (M)
+*
+*  WORK    (workspace) DOUBLE PRECISION array, dimension (LWORK)
+*
+*  LWORK   (input) INTEGER
+*          The dimension of the array WORK
+*
+*  RWORK   (workspace) DOUBLE PRECISION array
+*
+*  RESULT  (output) DOUBLE PRECISION array, dimension (9)
+*          The test ratios:
+*          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.
+*          ( EPS2 = MAX( norm( I - X'*X ) / M, ULP ). )
+*
+*  =====================================================================
+*
+*     .. Parameters ..
+      DOUBLE PRECISION   PIOVER2, REALONE, REALZERO
+      PARAMETER          ( PIOVER2 = 1.57079632679489662D0,
+     $                     REALONE = 1.0D0, REALZERO = 0.0D0 )
+      DOUBLE PRECISION   ZERO, ONE
+      PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0 )
+*     ..
+*     .. Local Scalars ..
+      INTEGER            I, INFO, R
+      DOUBLE PRECISION   EPS2, RESID, ULP, ULPINV
+*     ..
+*     .. External Functions ..
+      DOUBLE PRECISION   DLAMCH, DLANGE, DLANSY
+      EXTERNAL           DLAMCH, DLANGE, DLANSY
+*     ..
+*     .. External Subroutines ..
+      EXTERNAL           DGEMM, DLACPY, DLASET, DORCSD, DSYRK
+*     ..
+*     .. Intrinsic Functions ..
+      INTRINSIC          DBLE, MAX, MIN
+*     ..
+*     .. Executable Statements ..
+*
+      ULP = DLAMCH( 'Precision' )
+      ULPINV = REALONE / ULP
+      CALL DLASET( 'Full', M, M, ZERO, ONE, WORK, LDX )
+      CALL DSYRK( 'Upper', 'Conjugate transpose', M, M, -ONE, X, LDX,
+     $            ONE, WORK, LDX )
+      EPS2 = MAX( ULP, 
+     $            DLANGE( '1', M, M, WORK, LDX, RWORK ) / DBLE( M ) )
+      R = MIN( P, M-P, Q, M-Q )
+*
+*     Copy the matrix X to the array XF.
+*
+      CALL DLACPY( 'Full', M, M, X, LDX, XF, LDX )
+*
+*     Compute the CSD
+*
+      CALL DORCSD( 'Y', 'Y', 'Y', 'Y', 'N', 'D', M, P, Q, XF(1,1), LDX,
+     $             XF(1,Q+1), LDX, XF(P+1,1), LDX, XF(P+1,Q+1), LDX,
+     $             THETA, U1, LDU1, U2, LDU2, V1T, LDV1T, V2T, LDV2T,
+     $             WORK, LWORK, IWORK, INFO )
+*
+*     Compute X := diag(U1,U2)'*X*diag(V1,V2) - [D11 D12; D21 D22]
+*
+      CALL DGEMM( 'No transpose', 'Conjugate transpose', P, Q, Q, ONE,
+     $            X, LDX, V1T, LDV1T, ZERO, WORK, LDX )
+*
+      CALL DGEMM( 'Conjugate transpose', 'No transpose', P, Q, P, ONE,
+     $            U1, LDU1, WORK, LDX, ZERO, X, LDX )
+*
+      DO I = 1, MIN(P,Q)-R
+         X(I,I) = X(I,I) - ONE
+      END DO
+      DO I = 1, R
+         X(MIN(P,Q)-R+I,MIN(P,Q)-R+I) =
+     $           X(MIN(P,Q)-R+I,MIN(P,Q)-R+I) - COS(THETA(I))
+      END DO
+*
+      CALL DGEMM( 'No transpose', 'Conjugate transpose', P, M-Q, M-Q,
+     $            ONE, X(1,Q+1), LDX, V2T, LDV2T, ZERO, WORK, LDX )
+*
+      CALL DGEMM( 'Conjugate transpose', 'No transpose', P, M-Q, P,
+     $            ONE, U1, LDU1, WORK, LDX, ZERO, X(1,Q+1), LDX )
+*
+      DO I = 1, MIN(P,M-Q)-R
+         X(P-I+1,M-I+1) = X(P-I+1,M-I+1) + ONE
+      END DO
+      DO I = 1, R
+         X(P-(MIN(P,M-Q)-R)+1-I,M-(MIN(P,M-Q)-R)+1-I) =
+     $      X(P-(MIN(P,M-Q)-R)+1-I,M-(MIN(P,M-Q)-R)+1-I) +
+     $      SIN(THETA(R-I+1))
+      END DO
+*
+      CALL DGEMM( 'No transpose', 'Conjugate transpose', M-P, Q, Q, ONE,
+     $            X(P+1,1), LDX, V1T, LDV1T, ZERO, WORK, LDX )
+*
+      CALL DGEMM( 'Conjugate transpose', 'No transpose', M-P, Q, M-P,
+     $            ONE, U2, LDU2, WORK, LDX, ZERO, X(P+1,1), LDX )
+*
+      DO I = 1, MIN(M-P,Q)-R
+         X(M-I+1,Q-I+1) = X(M-I+1,Q-I+1) - ONE
+      END DO
+      DO I = 1, R
+         X(M-(MIN(M-P,Q)-R)+1-I,Q-(MIN(M-P,Q)-R)+1-I) =
+     $             X(M-(MIN(M-P,Q)-R)+1-I,Q-(MIN(M-P,Q)-R)+1-I) -
+     $             SIN(THETA(R-I+1))
+      END DO
+*
+      CALL DGEMM( 'No transpose', 'Conjugate transpose', M-P, M-Q, M-Q,
+     $            ONE, X(P+1,Q+1), LDX, V2T, LDV2T, ZERO, WORK, LDX )
+*
+      CALL DGEMM( 'Conjugate transpose', 'No transpose', M-P, M-Q, M-P,
+     $            ONE, U2, LDU2, WORK, LDX, ZERO, X(P+1,Q+1), LDX )
+*
+      DO I = 1, MIN(M-P,M-Q)-R
+         X(P+I,Q+I) = X(P+I,Q+I) - ONE
+      END DO
+      DO I = 1, R
+         X(P+(MIN(M-P,M-Q)-R)+I,Q+(MIN(M-P,M-Q)-R)+I) =
+     $      X(P+(MIN(M-P,M-Q)-R)+I,Q+(MIN(M-P,M-Q)-R)+I) -
+     $      COS(THETA(I))
+      END DO
+*
+*     Compute norm( U1'*X11*V1 - D11 ) / ( MAX(1,P,Q)*EPS2 ) .
+*
+      RESID = DLANGE( '1', P, Q, X, LDX, RWORK )
+      RESULT( 1 ) = ( RESID / DBLE(MAX(1,P,Q)) ) / EPS2
+*
+*     Compute norm( U1'*X12*V2 - D12 ) / ( MAX(1,P,M-Q)*EPS2 ) .
+*
+      RESID = DLANGE( '1', P, M-Q, X(1,Q+1), LDX, RWORK )
+      RESULT( 2 ) = ( RESID / DBLE(MAX(1,P,M-Q)) ) / EPS2
+*
+*     Compute norm( U2'*X21*V1 - D21 ) / ( MAX(1,M-P,Q)*EPS2 ) .
+*
+      RESID = DLANGE( '1', M-P, Q, X(P+1,1), LDX, RWORK )
+      RESULT( 3 ) = ( RESID / DBLE(MAX(1,M-P,Q)) ) / EPS2
+*
+*     Compute norm( U2'*X22*V2 - D22 ) / ( MAX(1,M-P,M-Q)*EPS2 ) .
+*
+      RESID = DLANGE( '1', M-P, M-Q, X(P+1,Q+1), LDX, RWORK )
+      RESULT( 4 ) = ( RESID / DBLE(MAX(1,M-P,M-Q)) ) / EPS2
+*
+*     Compute I - U1'*U1
+*
+      CALL DLASET( 'Full', P, P, ZERO, ONE, WORK, LDU1 )
+      CALL DSYRK( 'Upper', 'Conjugate transpose', P, P, -ONE, U1, LDU1,
+     $            ONE, WORK, LDU1 )
+*
+*     Compute norm( I - U'*U ) / ( MAX(1,P) * ULP ) .
+*
+      RESID = DLANSY( '1', 'Upper', P, WORK, LDU1, RWORK )
+      RESULT( 5 ) = ( RESID / DBLE(MAX(1,P)) ) / ULP
+*
+*     Compute I - U2'*U2
+*
+      CALL DLASET( 'Full', M-P, M-P, ZERO, ONE, WORK, LDU2 )
+      CALL DSYRK( 'Upper', 'Conjugate transpose', M-P, M-P, -ONE, U2,
+     $            LDU2, ONE, WORK, LDU2 )
+*
+*     Compute norm( I - U2'*U2 ) / ( MAX(1,M-P) * ULP ) .
+*
+      RESID = DLANSY( '1', 'Upper', M-P, WORK, LDU2, RWORK )
+      RESULT( 6 ) = ( RESID / DBLE(MAX(1,M-P)) ) / ULP
+*
+*     Compute I - V1T*V1T'
+*
+      CALL DLASET( 'Full', Q, Q, ZERO, ONE, WORK, LDV1T )
+      CALL DSYRK( 'Upper', 'No transpose', Q, Q, -ONE, V1T, LDV1T, ONE,
+     $            WORK, LDV1T )
+*
+*     Compute norm( I - V1T*V1T' ) / ( MAX(1,Q) * ULP ) .
+*
+      RESID = DLANSY( '1', 'Upper', Q, WORK, LDV1T, RWORK )
+      RESULT( 7 ) = ( RESID / DBLE(MAX(1,Q)) ) / ULP
+*
+*     Compute I - V2T*V2T'
+*
+      CALL DLASET( 'Full', M-Q, M-Q, ZERO, ONE, WORK, LDV2T )
+      CALL DSYRK( 'Upper', 'No transpose', M-Q, M-Q, -ONE, V2T, LDV2T,
+     $            ONE, WORK, LDV2T )
+*
+*     Compute norm( I - V2T*V2T' ) / ( MAX(1,M-Q) * ULP ) .
+*
+      RESID = DLANSY( '1', 'Upper', M-Q, WORK, LDV2T, RWORK )
+      RESULT( 8 ) = ( RESID / DBLE(MAX(1,M-Q)) ) / ULP
+*
+*     Check sorting
+*
+      RESULT(9) = REALZERO
+      DO I = 1, R
+         IF( THETA(I).LT.REALZERO .OR. THETA(I).GT.PIOVER2 ) THEN
+            RESULT(9) = ULPINV
+         END IF
+         IF( I.GT.1 .AND. THETA(I).LT.THETA(I-1) ) THEN
+            RESULT(9) = ULPINV
+         END IF
+      END DO
+*
+      RETURN
+*      
+*     End of DCSDTS
+*
+      END
+
diff --git a/TESTING/EIG/sckcsd.f b/TESTING/EIG/sckcsd.f
new file mode 100644
index 00000000..e9d4646a
--- /dev/null
+++ b/TESTING/EIG/sckcsd.f
@@ -0,0 +1,296 @@
+      SUBROUTINE SCKCSD( NM, MVAL, PVAL, QVAL, NMATS, ISEED, THRESH,
+     $                   MMAX, X, XF, U1, U2, V1T, V2T, THETA, IWORK,
+     $                   WORK, RWORK, NIN, NOUT, INFO )
+      IMPLICIT NONE
+*
+*     Originally SCKGSV
+*  -- LAPACK test routine (version 3.1) --
+*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
+*     November 2006
+*
+*     Adapted to SCKCSD by Brian Sutton
+*     July 2010
+*
+*     .. Scalar Arguments ..
+      INTEGER            INFO, NIN, NM, NMATS, MMAX, NOUT
+      REAL               THRESH
+*     ..
+*     .. Array Arguments ..
+      INTEGER            ISEED( 4 ), IWORK( * ), MVAL( * ), PVAL( * ),
+     $                   QVAL( * )
+      REAL               RWORK( * ), THETA( * )
+      REAL               U1( * ), U2( * ), V1T( * ), V2T( * ),
+     $                   WORK( * ), X( * ), XF( * )
+*     ..
+*
+*  Purpose
+*  =======
+*
+*  SCKCSD tests SORCSD:
+*         the CSD for an M-by-M orthogonal matrix X partitioned as
+*         [ X11 X12; X21 X22 ]. X11 is P-by-Q.
+*
+*  Arguments
+*  =========
+*
+*  NM      (input) INTEGER
+*          The number of values of M contained in the vector MVAL.
+*
+*  MVAL    (input) INTEGER array, dimension (NM)
+*          The values of the matrix row dimension M.
+*
+*  PVAL    (input) INTEGER array, dimension (NM)
+*          The values of the matrix row dimension P.
+*
+*  QVAL    (input) INTEGER array, dimension (NM)
+*          The values of the matrix column dimension Q.
+*
+*  NMATS   (input) INTEGER
+*          The number of matrix types to be tested for each combination
+*          of matrix dimensions.  If NMATS >= NTYPES (the maximum
+*          number of matrix types), then all the different types are
+*          generated for testing.  If NMATS < NTYPES, another input line
+*          is read to get the numbers of the matrix types to be used.
+*
+*  ISEED   (input/output) INTEGER array, dimension (4)
+*          On entry, the seed of the random number generator.  The array
+*          elements should be between 0 and 4095, otherwise they will be
+*          reduced mod 4096, and ISEED(4) must be odd.
+*          On exit, the next seed in the random number sequence after
+*          all the test matrices have been generated.
+*
+*  THRESH  (input) REAL
+*          The threshold value for the test ratios.  A result is
+*          included in the output file if RESULT >= THRESH.  To have
+*          every test ratio printed, use THRESH = 0.
+*
+*  MMAX    (input) INTEGER
+*          The maximum value permitted for M, used in dimensioning the
+*          work arrays.
+*
+*  X       (workspace) REAL array, dimension (MMAX*MMAX)
+*
+*  XF      (workspace) REAL array, dimension (MMAX*MMAX)
+*
+*  U1      (workspace) REAL array, dimension (MMAX*MMAX)
+*
+*  U2      (workspace) REAL array, dimension (MMAX*MMAX)
+*
+*  V1T     (workspace) REAL array, dimension (MMAX*MMAX)
+*
+*  V2T     (workspace) REAL array, dimension (MMAX*MMAX)
+*
+*  THETA   (workspace) REAL array, dimension (MMAX)
+*
+*  IWORK   (workspace) INTEGER array, dimension (MMAX)
+*
+*  WORK    (workspace) REAL array
+*
+*  RWORK   (workspace) REAL array
+*
+*  NIN     (input) INTEGER
+*          The unit number for input.
+*
+*  NOUT    (input) INTEGER
+*          The unit number for output.
+*
+*  INFO    (output) INTEGER
+*          = 0 :  successful exit
+*          > 0 :  If SLAROR returns an error code, the absolute value
+*                 of it is returned.
+*
+*  =====================================================================
+*
+*     .. Parameters ..
+      INTEGER            NTESTS
+      PARAMETER          ( NTESTS = 9 )
+      INTEGER            NTYPES
+      PARAMETER          ( NTYPES = 4 )
+      REAL               GAPDIGIT, ORTH, PIOVER2, TEN
+      PARAMETER          ( GAPDIGIT = 10.0E0, ORTH = 1.0E-4,
+     $                     PIOVER2 = 1.57079632679489662E0,
+     $                     TEN = 10.0D0 )
+      REAL               ZERO, ONE
+      PARAMETER          ( ZERO = 0.0E0, ONE = 1.0E0 )
+*     ..
+*     .. Local Scalars ..
+      LOGICAL            FIRSTT
+      CHARACTER*3        PATH
+      INTEGER            I, IINFO, IM, IMAT, J, LDU1, LDU2, LDV1T,
+     $                   LDV2T, LDX, LWORK, M, NFAIL, NRUN, NT, P, Q, R
+*     ..
+*     .. Local Arrays ..
+      LOGICAL            DOTYPE( NTYPES )
+      DOUBLE PRECISION   RESULT( NTESTS )
+*     ..
+*     .. External Subroutines ..
+      EXTERNAL           ALAHDG, ALAREQ, ALASUM, SCSDTS, SLACSG, SLAROR,
+     $                   SLASET
+*     ..
+*     .. Intrinsic Functions ..
+      INTRINSIC          ABS, COS, MIN, SIN
+*     ..
+*     .. External Functions ..
+      REAL               SLANGE, SLARND
+      EXTERNAL           SLANGE, SLARND
+*     ..
+*     .. Executable Statements ..
+*
+*     Initialize constants and the random number seed.
+*
+      PATH( 1: 3 ) = 'CSD'
+      INFO = 0
+      NRUN = 0
+      NFAIL = 0
+      FIRSTT = .TRUE.
+      CALL ALAREQ( PATH, NMATS, DOTYPE, NTYPES, NIN, NOUT )
+      LDX = MMAX
+      LDU1 = MMAX
+      LDU2 = MMAX
+      LDV1T = MMAX
+      LDV2T = MMAX
+      LWORK = MMAX*MMAX
+*
+*     Do for each value of M in MVAL.
+*
+      DO 30 IM = 1, NM
+         M = MVAL( IM )
+         P = PVAL( IM )
+         Q = QVAL( IM )
+*
+         DO 20 IMAT = 1, NTYPES
+*
+*           Do the tests only if DOTYPE( IMAT ) is true.
+*
+            IF( .NOT.DOTYPE( IMAT ) )
+     $         GO TO 20
+*
+*           Generate X
+*
+            IF( IMAT.EQ.1 ) THEN
+               CALL SLAROR( 'L', 'I', M, M, X, LDX, ISEED, WORK, IINFO )
+               IF( M .NE. 0 .AND. IINFO .NE. 0 ) THEN
+                  WRITE( NOUT, FMT = 9999 ) M, IINFO
+                  INFO = ABS( IINFO )
+                  GO TO 20
+               END IF
+            ELSE IF( IMAT.EQ.2 ) THEN
+               R = MIN( P, M-P, Q, M-Q )
+               DO I = 1, R
+                  THETA(I) = PIOVER2 * SLARND( 1, ISEED )
+               END DO
+               CALL SLACSG( M, P, Q, THETA, ISEED, X, LDX, WORK )
+               DO I = 1, M
+                  DO J = 1, M
+                     X(I+(J-1)*LDX) = X(I+(J-1)*LDX) +
+     $                                ORTH*SLARND(2,ISEED)
+                  END DO
+               END DO
+            ELSE IF( IMAT.EQ.3 ) THEN
+               R = MIN( P, M-P, Q, M-Q )
+               DO I = 1, R+1
+                  THETA(I) = TEN**(-SLARND(1,ISEED)*GAPDIGIT)
+               END DO
+               DO I = 2, R+1
+                  THETA(I) = THETA(I-1) + THETA(I)
+               END DO
+               DO I = 1, R
+                  THETA(I) = PIOVER2 * THETA(I) / THETA(R+1)
+               END DO
+               CALL SLACSG( M, P, Q, THETA, ISEED, X, LDX, WORK )
+            ELSE
+            END IF
+*
+            NT = 9
+*
+            CALL SCSDTS( M, P, Q, X, XF, LDX, U1, LDU1, U2, LDU2, V1T,
+     $                   LDV1T, V2T, LDV2T, THETA, IWORK, WORK, LWORK,
+     $                   RWORK, RESULT )
+*
+*           Print information about the tests that did not
+*           pass the threshold.
+*
+            DO 10 I = 1, NT
+               IF( RESULT( I ).GE.THRESH ) THEN
+                  IF( NFAIL.EQ.0 .AND. FIRSTT ) THEN
+                     FIRSTT = .FALSE.
+                     CALL ALAHDG( NOUT, PATH )
+                  END IF
+                  WRITE( NOUT, FMT = 9998 )M, P, Q, IMAT, I,
+     $               RESULT( I )
+                  NFAIL = NFAIL + 1
+               END IF
+   10       CONTINUE
+            NRUN = NRUN + NT
+   20    CONTINUE
+   30 CONTINUE
+*
+*     Print a summary of the results.
+*
+      CALL ALASUM( PATH, NOUT, NFAIL, NRUN, 0 )
+*
+ 9999 FORMAT( ' SLAROR in SCKCSD: M = ', I5, ', INFO = ', I15 )
+ 9998 FORMAT( ' M=', I4, ' P=', I4, ', Q=', I4, ', type ', I2,
+     $      ', test ', I2, ', ratio=', G13.6 )
+      RETURN
+*
+*     End of SCKCSD
+*
+      END
+*
+*
+*
+      SUBROUTINE SLACSG( M, P, Q, THETA, ISEED, X, LDX, WORK )
+      IMPLICIT NONE
+*
+      INTEGER            LDX, M, P, Q
+      INTEGER            ISEED( 4 )
+      REAL               THETA( * )
+      REAL               WORK( * ), X( LDX, * )
+*
+      REAL               ONE, ZERO
+      PARAMETER          ( ONE = 1.0E0, ZERO = 0.0E0 )
+*
+      INTEGER            I, INFO, R
+*
+      R = MIN( P, M-P, Q, M-Q )
+*
+      CALL SLASET( 'Full', M, M, ZERO, ZERO, X, LDX )
+*
+      DO I = 1, MIN(P,Q)-R
+         X(I,I) = ONE
+      END DO
+      DO I = 1, R
+         X(MIN(P,Q)-R+I,MIN(P,Q)-R+I) = COS(THETA(I))
+      END DO
+      DO I = 1, MIN(P,M-Q)-R
+         X(P-I+1,M-I+1) = -ONE
+      END DO
+      DO I = 1, R
+         X(P-(MIN(P,M-Q)-R)+1-I,M-(MIN(P,M-Q)-R)+1-I) =
+     $      -SIN(THETA(R-I+1))
+      END DO
+      DO I = 1, MIN(M-P,Q)-R
+         X(M-I+1,Q-I+1) = ONE
+      END DO
+      DO I = 1, R
+         X(M-(MIN(M-P,Q)-R)+1-I,Q-(MIN(M-P,Q)-R)+1-I) =
+     $      SIN(THETA(R-I+1))
+      END DO
+      DO I = 1, MIN(M-P,M-Q)-R
+         X(P+I,Q+I) = ONE
+      END DO
+      DO I = 1, R
+         X(P+(MIN(M-P,M-Q)-R)+I,Q+(MIN(M-P,M-Q)-R)+I) =
+     $      COS(THETA(I))
+      END DO
+      CALL SLAROR( 'Left', 'No init', P, M, X, LDX, ISEED, WORK, INFO )
+      CALL SLAROR( 'Left', 'No init', M-P, M, X(P+1,1), LDX,
+     $             ISEED, WORK, INFO )
+      CALL SLAROR( 'Right', 'No init', M, Q, X, LDX, ISEED,
+     $             WORK, INFO )
+      CALL SLAROR( 'Right', 'No init', M, M-Q,
+     $             X(1,Q+1), LDX, ISEED, WORK, INFO )
+*
+      END
+
diff --git a/TESTING/EIG/scsdts.f b/TESTING/EIG/scsdts.f
new file mode 100644
index 00000000..6f1b8264
--- /dev/null
+++ b/TESTING/EIG/scsdts.f
@@ -0,0 +1,312 @@
+      SUBROUTINE SCSDTS( M, P, Q, X, XF, LDX, U1, LDU1, U2, LDU2, V1T,
+     $                   LDV1T, V2T, LDV2T, THETA, IWORK, WORK, LWORK,
+     $                   RWORK, RESULT )
+      IMPLICIT NONE
+*
+*     Originally xGSVTS
+*  -- LAPACK test routine (version 3.1) --
+*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
+*     November 2006
+*
+*     Adapted to SCSDTS by Brian Sutton
+*     July 2010
+*
+*     .. Scalar Arguments ..
+      INTEGER            LDX, LDU1, LDU2, LDV1T, LDV2T, LWORK, M, P, Q
+*     ..
+*     .. Array Arguments ..
+      INTEGER            IWORK( * )
+      REAL               RESULT( 9 ), RWORK( * ), THETA( * )
+      REAL               U1( LDU1, * ), U2( LDU2, * ), V1T( LDV1T, * ),
+     $                   V2T( LDV2T, * ), WORK( LWORK ), X( LDX, * ),
+     $                   XF( LDX, * )
+*     ..
+*
+*  Purpose
+*  =======
+*
+*  SCSDTS tests SORCSD, 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 ]
+*
+*  Arguments
+*  =========
+*
+*  M       (input) INTEGER
+*          The number of rows of the matrix X.  M >= 0.
+*
+*  P       (input) INTEGER
+*          The number of rows of the matrix X11.  P >= 0.
+*
+*  Q       (input) INTEGER
+*          The number of columns of the matrix X11.  Q >= 0.
+*
+*  X       (input) REAL array, dimension (LDX,M)
+*          The M-by-M matrix X.
+*
+*  XF      (output) REAL array, dimension (LDX,M)
+*          Details of the CSD of X, as returned by SORCSD;
+*          see SORCSD for further details.
+*
+*  LDX     (input) INTEGER
+*          The leading dimension of the arrays X and XF.
+*          LDX >= max( 1,M ).
+*
+*  U1      (output) REAL array, dimension(LDU1,P)
+*          The P-by-P orthogonal matrix U1.
+*
+*  LDU1    (input) INTEGER
+*          The leading dimension of the array U1. LDU >= max(1,P).
+*
+*  U2      (output) REAL array, dimension(LDU2,M-P)
+*          The (M-P)-by-(M-P) orthogonal matrix U2.
+*
+*  LDU2    (input) INTEGER
+*          The leading dimension of the array U2. LDU >= max(1,M-P).
+*
+*  V1T     (output) REAL array, dimension(LDV1T,Q)
+*          The Q-by-Q orthogonal matrix V1T.
+*
+*  LDV1T   (input) INTEGER
+*          The leading dimension of the array V1T. LDV1T >=
+*          max(1,Q).
+*
+*  V2T     (output) REAL array, dimension(LDV2T,M-Q)
+*          The (M-Q)-by-(M-Q) orthogonal matrix V2T.
+*
+*  LDV2T   (input) INTEGER
+*          The leading dimension of the array V2T. LDV2T >=
+*          max(1,M-Q).
+*
+*  THETA   (output) 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 SORCSD for
+*          details.
+*
+*  IWORK   (workspace) INTEGER array, dimension (M)
+*
+*  WORK    (workspace) REAL array, dimension (LWORK)
+*
+*  LWORK   (input) INTEGER
+*          The dimension of the array WORK
+*
+*  RWORK   (workspace) REAL array
+*
+*  RESULT  (output) REAL array, dimension (9)
+*          The test ratios:
+*          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.
+*          ( EPS2 = MAX( norm( I - X'*X ) / M, ULP ). )
+*
+*  =====================================================================
+*
+*     .. Parameters ..
+      REAL               PIOVER2, REALONE, REALZERO
+      PARAMETER          ( PIOVER2 = 1.57079632679489662E0,
+     $                     REALONE = 1.0E0, REALZERO = 0.0E0 )
+      REAL               ZERO, ONE
+      PARAMETER          ( ZERO = 0.0E0, ONE = 1.0E0 )
+*     ..
+*     .. Local Scalars ..
+      INTEGER            I, INFO, R
+      REAL               EPS2, RESID, ULP, ULPINV
+*     ..
+*     .. External Functions ..
+      REAL               SLAMCH, SLANGE, SLANSY
+      EXTERNAL           SLAMCH, SLANGE, SLANSY
+*     ..
+*     .. External Subroutines ..
+      EXTERNAL           SGEMM, SLACPY, SLASET, SORCSD, SSYRK
+*     ..
+*     .. Intrinsic Functions ..
+      INTRINSIC          REAL, MAX, MIN
+*     ..
+*     .. Executable Statements ..
+*
+      ULP = SLAMCH( 'Precision' )
+      ULPINV = REALONE / ULP
+      CALL SLASET( 'Full', M, M, ZERO, ONE, WORK, LDX )
+      CALL SSYRK( 'Upper', 'Conjugate transpose', M, M, -ONE, X, LDX,
+     $            ONE, WORK, LDX )
+      EPS2 = MAX( ULP, 
+     $            SLANGE( '1', M, M, WORK, LDX, RWORK ) / REAL( M ) )
+      R = MIN( P, M-P, Q, M-Q )
+*
+*     Copy the matrix X to the array XF.
+*
+      CALL SLACPY( 'Full', M, M, X, LDX, XF, LDX )
+*
+*     Compute the CSD
+*
+      CALL SORCSD( 'Y', 'Y', 'Y', 'Y', 'N', 'D', M, P, Q, XF(1,1), LDX,
+     $             XF(1,Q+1), LDX, XF(P+1,1), LDX, XF(P+1,Q+1), LDX,
+     $             THETA, U1, LDU1, U2, LDU2, V1T, LDV1T, V2T, LDV2T,
+     $             WORK, LWORK, IWORK, INFO )
+*
+*     Compute X := diag(U1,U2)'*X*diag(V1,V2) - [D11 D12; D21 D22]
+*
+      CALL SGEMM( 'No transpose', 'Conjugate transpose', P, Q, Q, ONE,
+     $            X, LDX, V1T, LDV1T, ZERO, WORK, LDX )
+*
+      CALL SGEMM( 'Conjugate transpose', 'No transpose', P, Q, P, ONE,
+     $            U1, LDU1, WORK, LDX, ZERO, X, LDX )
+*
+      DO I = 1, MIN(P,Q)-R
+         X(I,I) = X(I,I) - ONE
+      END DO
+      DO I = 1, R
+         X(MIN(P,Q)-R+I,MIN(P,Q)-R+I) =
+     $           X(MIN(P,Q)-R+I,MIN(P,Q)-R+I) - COS(THETA(I))
+      END DO
+*
+      CALL SGEMM( 'No transpose', 'Conjugate transpose', P, M-Q, M-Q,
+     $            ONE, X(1,Q+1), LDX, V2T, LDV2T, ZERO, WORK, LDX )
+*
+      CALL SGEMM( 'Conjugate transpose', 'No transpose', P, M-Q, P,
+     $            ONE, U1, LDU1, WORK, LDX, ZERO, X(1,Q+1), LDX )
+*
+      DO I = 1, MIN(P,M-Q)-R
+         X(P-I+1,M-I+1) = X(P-I+1,M-I+1) + ONE
+      END DO
+      DO I = 1, R
+         X(P-(MIN(P,M-Q)-R)+1-I,M-(MIN(P,M-Q)-R)+1-I) =
+     $      X(P-(MIN(P,M-Q)-R)+1-I,M-(MIN(P,M-Q)-R)+1-I) +
+     $      SIN(THETA(R-I+1))
+      END DO
+*
+      CALL SGEMM( 'No transpose', 'Conjugate transpose', M-P, Q, Q, ONE,
+     $            X(P+1,1), LDX, V1T, LDV1T, ZERO, WORK, LDX )
+*
+      CALL SGEMM( 'Conjugate transpose', 'No transpose', M-P, Q, M-P,
+     $            ONE, U2, LDU2, WORK, LDX, ZERO, X(P+1,1), LDX )
+*
+      DO I = 1, MIN(M-P,Q)-R
+         X(M-I+1,Q-I+1) = X(M-I+1,Q-I+1) - ONE
+      END DO
+      DO I = 1, R
+         X(M-(MIN(M-P,Q)-R)+1-I,Q-(MIN(M-P,Q)-R)+1-I) =
+     $             X(M-(MIN(M-P,Q)-R)+1-I,Q-(MIN(M-P,Q)-R)+1-I) -
+     $             SIN(THETA(R-I+1))
+      END DO
+*
+      CALL SGEMM( 'No transpose', 'Conjugate transpose', M-P, M-Q, M-Q,
+     $            ONE, X(P+1,Q+1), LDX, V2T, LDV2T, ZERO, WORK, LDX )
+*
+      CALL SGEMM( 'Conjugate transpose', 'No transpose', M-P, M-Q, M-P,
+     $            ONE, U2, LDU2, WORK, LDX, ZERO, X(P+1,Q+1), LDX )
+*
+      DO I = 1, MIN(M-P,M-Q)-R
+         X(P+I,Q+I) = X(P+I,Q+I) - ONE
+      END DO
+      DO I = 1, R
+         X(P+(MIN(M-P,M-Q)-R)+I,Q+(MIN(M-P,M-Q)-R)+I) =
+     $      X(P+(MIN(M-P,M-Q)-R)+I,Q+(MIN(M-P,M-Q)-R)+I) -
+     $      COS(THETA(I))
+      END DO
+*
+*     Compute norm( U1'*X11*V1 - D11 ) / ( MAX(1,P,Q)*EPS2 ) .
+*
+      RESID = SLANGE( '1', P, Q, X, LDX, RWORK )
+      RESULT( 1 ) = ( RESID / REAL(MAX(1,P,Q)) ) / EPS2
+*
+*     Compute norm( U1'*X12*V2 - D12 ) / ( MAX(1,P,M-Q)*EPS2 ) .
+*
+      RESID = SLANGE( '1', P, M-Q, X(1,Q+1), LDX, RWORK )
+      RESULT( 2 ) = ( RESID / REAL(MAX(1,P,M-Q)) ) / EPS2
+*
+*     Compute norm( U2'*X21*V1 - D21 ) / ( MAX(1,M-P,Q)*EPS2 ) .
+*
+      RESID = SLANGE( '1', M-P, Q, X(P+1,1), LDX, RWORK )
+      RESULT( 3 ) = ( RESID / REAL(MAX(1,M-P,Q)) ) / EPS2
+*
+*     Compute norm( U2'*X22*V2 - D22 ) / ( MAX(1,M-P,M-Q)*EPS2 ) .
+*
+      RESID = SLANGE( '1', M-P, M-Q, X(P+1,Q+1), LDX, RWORK )
+      RESULT( 4 ) = ( RESID / REAL(MAX(1,M-P,M-Q)) ) / EPS2
+*
+*     Compute I - U1'*U1
+*
+      CALL SLASET( 'Full', P, P, ZERO, ONE, WORK, LDU1 )
+      CALL SSYRK( 'Upper', 'Conjugate transpose', P, P, -ONE, U1, LDU1,
+     $            ONE, WORK, LDU1 )
+*
+*     Compute norm( I - U'*U ) / ( MAX(1,P) * ULP ) .
+*
+      RESID = SLANSY( '1', 'Upper', P, WORK, LDU1, RWORK )
+      RESULT( 5 ) = ( RESID / REAL(MAX(1,P)) ) / ULP
+*
+*     Compute I - U2'*U2
+*
+      CALL SLASET( 'Full', M-P, M-P, ZERO, ONE, WORK, LDU2 )
+      CALL SSYRK( 'Upper', 'Conjugate transpose', M-P, M-P, -ONE, U2,
+     $            LDU2, ONE, WORK, LDU2 )
+*
+*     Compute norm( I - U2'*U2 ) / ( MAX(1,M-P) * ULP ) .
+*
+      RESID = SLANSY( '1', 'Upper', M-P, WORK, LDU2, RWORK )
+      RESULT( 6 ) = ( RESID / REAL(MAX(1,M-P)) ) / ULP
+*
+*     Compute I - V1T*V1T'
+*
+      CALL SLASET( 'Full', Q, Q, ZERO, ONE, WORK, LDV1T )
+      CALL SSYRK( 'Upper', 'No transpose', Q, Q, -ONE, V1T, LDV1T, ONE,
+     $            WORK, LDV1T )
+*
+*     Compute norm( I - V1T*V1T' ) / ( MAX(1,Q) * ULP ) .
+*
+      RESID = SLANSY( '1', 'Upper', Q, WORK, LDV1T, RWORK )
+      RESULT( 7 ) = ( RESID / REAL(MAX(1,Q)) ) / ULP
+*
+*     Compute I - V2T*V2T'
+*
+      CALL SLASET( 'Full', M-Q, M-Q, ZERO, ONE, WORK, LDV2T )
+      CALL SSYRK( 'Upper', 'No transpose', M-Q, M-Q, -ONE, V2T, LDV2T,
+     $            ONE, WORK, LDV2T )
+*
+*     Compute norm( I - V2T*V2T' ) / ( MAX(1,M-Q) * ULP ) .
+*
+      RESID = SLANSY( '1', 'Upper', M-Q, WORK, LDV2T, RWORK )
+      RESULT( 8 ) = ( RESID / REAL(MAX(1,M-Q)) ) / ULP
+*
+*     Check sorting
+*
+      RESULT(9) = REALZERO
+      DO I = 1, R
+         IF( THETA(I).LT.REALZERO .OR. THETA(I).GT.PIOVER2 ) THEN
+            RESULT(9) = ULPINV
+         END IF
+         IF( I.GT.1 .AND. THETA(I).LT.THETA(I-1) ) THEN
+            RESULT(9) = ULPINV
+         END IF
+      END DO
+*
+      RETURN
+*      
+*     End of SCSDTS
+*
+      END
+
diff --git a/TESTING/EIG/zckcsd.f b/TESTING/EIG/zckcsd.f
new file mode 100644
index 00000000..4cc7b921
--- /dev/null
+++ b/TESTING/EIG/zckcsd.f
@@ -0,0 +1,296 @@
+      SUBROUTINE ZCKCSD( NM, MVAL, PVAL, QVAL, NMATS, ISEED, THRESH,
+     $                   MMAX, X, XF, U1, U2, V1T, V2T, THETA, IWORK,
+     $                   WORK, RWORK, NIN, NOUT, INFO )
+      IMPLICIT NONE
+*
+*     Originally ZCKGSV
+*  -- LAPACK test routine (version 3.1) --
+*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
+*     November 2006
+*
+*     Adapted to ZCKCSD by Brian Sutton
+*     July 2010
+*
+*     .. Scalar Arguments ..
+      INTEGER            INFO, NIN, NM, NMATS, MMAX, NOUT
+      DOUBLE PRECISION   THRESH
+*     ..
+*     .. Array Arguments ..
+      INTEGER            ISEED( 4 ), IWORK( * ), MVAL( * ), PVAL( * ),
+     $                   QVAL( * )
+      DOUBLE PRECISION   RWORK( * ), THETA( * )
+      COMPLEX*16         U1( * ), U2( * ), V1T( * ), V2T( * ),
+     $                   WORK( * ), X( * ), XF( * )
+*     ..
+*
+*  Purpose
+*  =======
+*
+*  ZCKCSD tests ZUNCSD:
+*         the CSD for an M-by-M unitary matrix X partitioned as
+*         [ X11 X12; X21 X22 ]. X11 is P-by-Q.
+*
+*  Arguments
+*  =========
+*
+*  NM      (input) INTEGER
+*          The number of values of M contained in the vector MVAL.
+*
+*  MVAL    (input) INTEGER array, dimension (NM)
+*          The values of the matrix row dimension M.
+*
+*  PVAL    (input) INTEGER array, dimension (NM)
+*          The values of the matrix row dimension P.
+*
+*  QVAL    (input) INTEGER array, dimension (NM)
+*          The values of the matrix column dimension Q.
+*
+*  NMATS   (input) INTEGER
+*          The number of matrix types to be tested for each combination
+*          of matrix dimensions.  If NMATS >= NTYPES (the maximum
+*          number of matrix types), then all the different types are
+*          generated for testing.  If NMATS < NTYPES, another input line
+*          is read to get the numbers of the matrix types to be used.
+*
+*  ISEED   (input/output) INTEGER array, dimension (4)
+*          On entry, the seed of the random number generator.  The array
+*          elements should be between 0 and 4095, otherwise they will be
+*          reduced mod 4096, and ISEED(4) must be odd.
+*          On exit, the next seed in the random number sequence after
+*          all the test matrices have been generated.
+*
+*  THRESH  (input) DOUBLE PRECISION
+*          The threshold value for the test ratios.  A result is
+*          included in the output file if RESULT >= THRESH.  To have
+*          every test ratio printed, use THRESH = 0.
+*
+*  MMAX    (input) INTEGER
+*          The maximum value permitted for M, used in dimensioning the
+*          work arrays.
+*
+*  X       (workspace) COMPLEX*16 array, dimension (MMAX*MMAX)
+*
+*  XF      (workspace) COMPLEX*16 array, dimension (MMAX*MMAX)
+*
+*  U1      (workspace) COMPLEX*16 array, dimension (MMAX*MMAX)
+*
+*  U2      (workspace) COMPLEX*16 array, dimension (MMAX*MMAX)
+*
+*  V1T     (workspace) COMPLEX*16 array, dimension (MMAX*MMAX)
+*
+*  V2T     (workspace) COMPLEX*16 array, dimension (MMAX*MMAX)
+*
+*  THETA   (workspace) DOUBLE PRECISION array, dimension (MMAX)
+*
+*  IWORK   (workspace) INTEGER array, dimension (MMAX)
+*
+*  WORK    (workspace) COMPLEX*16 array
+*
+*  RWORK   (workspace) DOUBLE PRECISION array
+*
+*  NIN     (input) INTEGER
+*          The unit number for input.
+*
+*  NOUT    (input) INTEGER
+*          The unit number for output.
+*
+*  INFO    (output) INTEGER
+*          = 0 :  successful exit
+*          > 0 :  If ZLAROR returns an error code, the absolute value
+*                 of it is returned.
+*
+*  =====================================================================
+*
+*     .. Parameters ..
+      INTEGER            NTESTS
+      PARAMETER          ( NTESTS = 9 )
+      INTEGER            NTYPES
+      PARAMETER          ( NTYPES = 4 )
+      DOUBLE PRECISION   GAPDIGIT, ORTH, PIOVER2, TEN
+      PARAMETER          ( GAPDIGIT = 18.0D0, ORTH = 1.0D-12,
+     $                     PIOVER2 = 1.57079632679489662D0,
+     $                     TEN = 10.0D0 )
+      COMPLEX*16         ZERO, ONE
+      PARAMETER          ( ZERO = (0.0D0,0.0D0), ONE = (1.0D0,0.0D0) )
+*     ..
+*     .. Local Scalars ..
+      LOGICAL            FIRSTT
+      CHARACTER*3        PATH
+      INTEGER            I, IINFO, IM, IMAT, J, LDU1, LDU2, LDV1T,
+     $                   LDV2T, LDX, LWORK, M, NFAIL, NRUN, NT, P, Q, R
+*     ..
+*     .. Local Arrays ..
+      LOGICAL            DOTYPE( NTYPES )
+      DOUBLE PRECISION   RESULT( NTESTS )
+*     ..
+*     .. External Subroutines ..
+      EXTERNAL           ALAHDG, ALAREQ, ALASUM, ZCSDTS, ZLACSG, ZLAROR,
+     $                   ZLASET
+*     ..
+*     .. Intrinsic Functions ..
+      INTRINSIC          ABS, COS, MIN, SIN
+*     ..
+*     .. External Functions ..
+      DOUBLE PRECISION   ZLANGE, DLARND
+      EXTERNAL           ZLANGE, DLARND
+*     ..
+*     .. Executable Statements ..
+*
+*     Initialize constants and the random number seed.
+*
+      PATH( 1: 3 ) = 'CSD'
+      INFO = 0
+      NRUN = 0
+      NFAIL = 0
+      FIRSTT = .TRUE.
+      CALL ALAREQ( PATH, NMATS, DOTYPE, NTYPES, NIN, NOUT )
+      LDX = MMAX
+      LDU1 = MMAX
+      LDU2 = MMAX
+      LDV1T = MMAX
+      LDV2T = MMAX
+      LWORK = MMAX*MMAX
+*
+*     Do for each value of M in MVAL.
+*
+      DO 30 IM = 1, NM
+         M = MVAL( IM )
+         P = PVAL( IM )
+         Q = QVAL( IM )
+*
+         DO 20 IMAT = 1, NTYPES
+*
+*           Do the tests only if DOTYPE( IMAT ) is true.
+*
+            IF( .NOT.DOTYPE( IMAT ) )
+     $         GO TO 20
+*
+*           Generate X
+*
+            IF( IMAT.EQ.1 ) THEN
+               CALL ZLAROR( 'L', 'I', M, M, X, LDX, ISEED, WORK, IINFO )
+               IF( M .NE. 0 .AND. IINFO .NE. 0 ) THEN
+                  WRITE( NOUT, FMT = 9999 ) M, IINFO
+                  INFO = ABS( IINFO )
+                  GO TO 20
+               END IF
+            ELSE IF( IMAT.EQ.2 ) THEN
+               R = MIN( P, M-P, Q, M-Q )
+               DO I = 1, R
+                  THETA(I) = PIOVER2 * DLARND( 1, ISEED )
+               END DO
+               CALL ZLACSG( M, P, Q, THETA, ISEED, X, LDX, WORK )
+               DO I = 1, M
+                  DO J = 1, M
+                     X(I+(J-1)*LDX) = X(I+(J-1)*LDX) +
+     $                                ORTH*DLARND(2,ISEED)
+                  END DO
+               END DO
+            ELSE IF( IMAT.EQ.3 ) THEN
+               R = MIN( P, M-P, Q, M-Q )
+               DO I = 1, R+1
+                  THETA(I) = TEN**(-DLARND(1,ISEED)*GAPDIGIT)
+               END DO
+               DO I = 2, R+1
+                  THETA(I) = THETA(I-1) + THETA(I)
+               END DO
+               DO I = 1, R
+                  THETA(I) = PIOVER2 * THETA(I) / THETA(R+1)
+               END DO
+               CALL ZLACSG( M, P, Q, THETA, ISEED, X, LDX, WORK )
+            ELSE
+            END IF
+*
+            NT = 9
+*
+            CALL ZCSDTS( M, P, Q, X, XF, LDX, U1, LDU1, U2, LDU2, V1T,
+     $                   LDV1T, V2T, LDV2T, THETA, IWORK, WORK, LWORK,
+     $                   RWORK, RESULT )
+*
+*           Print information about the tests that did not
+*           pass the threshold.
+*
+            DO 10 I = 1, NT
+               IF( RESULT( I ).GE.THRESH ) THEN
+                  IF( NFAIL.EQ.0 .AND. FIRSTT ) THEN
+                     FIRSTT = .FALSE.
+                     CALL ALAHDG( NOUT, PATH )
+                  END IF
+                  WRITE( NOUT, FMT = 9998 )M, P, Q, IMAT, I,
+     $               RESULT( I )
+                  NFAIL = NFAIL + 1
+               END IF
+   10       CONTINUE
+            NRUN = NRUN + NT
+   20    CONTINUE
+   30 CONTINUE
+*
+*     Print a summary of the results.
+*
+      CALL ALASUM( PATH, NOUT, NFAIL, NRUN, 0 )
+*
+ 9999 FORMAT( ' ZLAROR in ZCKCSD: M = ', I5, ', INFO = ', I15 )
+ 9998 FORMAT( ' M=', I4, ' P=', I4, ', Q=', I4, ', type ', I2,
+     $      ', test ', I2, ', ratio=', G13.6 )
+      RETURN
+*
+*     End of ZCKCSD
+*
+      END
+*
+*
+*
+      SUBROUTINE ZLACSG( M, P, Q, THETA, ISEED, X, LDX, WORK )
+      IMPLICIT NONE
+*
+      INTEGER            LDX, M, P, Q
+      INTEGER            ISEED( 4 )
+      DOUBLE PRECISION   THETA( * )
+      COMPLEX*16         WORK( * ), X( LDX, * )
+*
+      COMPLEX*16         ONE, ZERO
+      PARAMETER          ( ONE = (1.0D0,0.0D0), ZERO = (0.0D0,0.0D0) )
+*
+      INTEGER            I, INFO, R
+*
+      R = MIN( P, M-P, Q, M-Q )
+*
+      CALL ZLASET( 'Full', M, M, ZERO, ZERO, X, LDX )
+*
+      DO I = 1, MIN(P,Q)-R
+         X(I,I) = ONE
+      END DO
+      DO I = 1, R
+         X(MIN(P,Q)-R+I,MIN(P,Q)-R+I) = COMPLEX( COS(THETA(I)), 0.0D0 )
+      END DO
+      DO I = 1, MIN(P,M-Q)-R
+         X(P-I+1,M-I+1) = -ONE
+      END DO
+      DO I = 1, R
+         X(P-(MIN(P,M-Q)-R)+1-I,M-(MIN(P,M-Q)-R)+1-I) =
+     $      COMPLEX( -SIN(THETA(R-I+1)), 0.0D0 )
+      END DO
+      DO I = 1, MIN(M-P,Q)-R
+         X(M-I+1,Q-I+1) = ONE
+      END DO
+      DO I = 1, R
+         X(M-(MIN(M-P,Q)-R)+1-I,Q-(MIN(M-P,Q)-R)+1-I) =
+     $      COMPLEX( SIN(THETA(R-I+1)), 0.0D0 )
+      END DO
+      DO I = 1, MIN(M-P,M-Q)-R
+         X(P+I,Q+I) = ONE
+      END DO
+      DO I = 1, R
+         X(P+(MIN(M-P,M-Q)-R)+I,Q+(MIN(M-P,M-Q)-R)+I) =
+     $      COMPLEX( COS(THETA(I)), 0.0D0 )
+      END DO
+      CALL ZLAROR( 'Left', 'No init', P, M, X, LDX, ISEED, WORK, INFO )
+      CALL ZLAROR( 'Left', 'No init', M-P, M, X(P+1,1), LDX,
+     $             ISEED, WORK, INFO )
+      CALL ZLAROR( 'Right', 'No init', M, Q, X, LDX, ISEED,
+     $             WORK, INFO )
+      CALL ZLAROR( 'Right', 'No init', M, M-Q,
+     $             X(1,Q+1), LDX, ISEED, WORK, INFO )
+*
+      END
+
diff --git a/TESTING/EIG/zcsdts.f b/TESTING/EIG/zcsdts.f
new file mode 100644
index 00000000..ae9a48ce
--- /dev/null
+++ b/TESTING/EIG/zcsdts.f
@@ -0,0 +1,313 @@
+      SUBROUTINE ZCSDTS( M, P, Q, X, XF, LDX, U1, LDU1, U2, LDU2, V1T,
+     $                   LDV1T, V2T, LDV2T, THETA, IWORK, WORK, LWORK,
+     $                   RWORK, RESULT )
+      IMPLICIT NONE
+*
+*     Originally xGSVTS
+*  -- LAPACK test routine (version 3.1) --
+*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
+*     November 2006
+*
+*     Adapted to ZCSDTS by Brian Sutton
+*     July 2010
+*
+*     .. Scalar Arguments ..
+      INTEGER            LDX, LDU1, LDU2, LDV1T, LDV2T, LWORK, M, P, Q
+*     ..
+*     .. Array Arguments ..
+      INTEGER            IWORK( * )
+      DOUBLE PRECISION   RESULT( 9 ), RWORK( * ), THETA( * )
+      COMPLEX*16         U1( LDU1, * ), U2( LDU2, * ), V1T( LDV1T, * ),
+     $                   V2T( LDV2T, * ), WORK( LWORK ), X( LDX, * ),
+     $                   XF( LDX, * )
+*     ..
+*
+*  Purpose
+*  =======
+*
+*  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 ]
+*
+*  Arguments
+*  =========
+*
+*  M       (input) INTEGER
+*          The number of rows of the matrix X.  M >= 0.
+*
+*  P       (input) INTEGER
+*          The number of rows of the matrix X11.  P >= 0.
+*
+*  Q       (input) INTEGER
+*          The number of columns of the matrix X11.  Q >= 0.
+*
+*  X       (input) COMPLEX*16 array, dimension (LDX,M)
+*          The M-by-M matrix X.
+*
+*  XF      (output) COMPLEX*16 array, dimension (LDX,M)
+*          Details of the CSD of X, as returned by ZUNCSD;
+*          see ZUNCSD for further details.
+*
+*  LDX     (input) INTEGER
+*          The leading dimension of the arrays X and XF.
+*          LDX >= max( 1,M ).
+*
+*  U1      (output) COMPLEX*16 array, dimension(LDU1,P)
+*          The P-by-P unitary matrix U1.
+*
+*  LDU1    (input) INTEGER
+*          The leading dimension of the array U1. LDU >= max(1,P).
+*
+*  U2      (output) COMPLEX*16 array, dimension(LDU2,M-P)
+*          The (M-P)-by-(M-P) unitary matrix U2.
+*
+*  LDU2    (input) INTEGER
+*          The leading dimension of the array U2. LDU >= max(1,M-P).
+*
+*  V1T     (output) COMPLEX*16 array, dimension(LDV1T,Q)
+*          The Q-by-Q unitary matrix V1T.
+*
+*  LDV1T   (input) INTEGER
+*          The leading dimension of the array V1T. LDV1T >=
+*          max(1,Q).
+*
+*  V2T     (output) COMPLEX*16 array, dimension(LDV2T,M-Q)
+*          The (M-Q)-by-(M-Q) unitary matrix V2T.
+*
+*  LDV2T   (input) INTEGER
+*          The leading dimension of the array V2T. LDV2T >=
+*          max(1,M-Q).
+*
+*  THETA   (output) 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   (workspace) INTEGER array, dimension (M)
+*
+*  WORK    (workspace) COMPLEX*16 array, dimension (LWORK)
+*
+*  LWORK   (input) INTEGER
+*          The dimension of the array WORK
+*
+*  RWORK   (workspace) DOUBLE PRECISION array
+*
+*  RESULT  (output) DOUBLE PRECISION array, dimension (9)
+*          The test ratios:
+*          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.
+*          ( EPS2 = MAX( norm( I - X'*X ) / M, ULP ). )
+*
+*  =====================================================================
+*
+*     .. Parameters ..
+      DOUBLE PRECISION   PIOVER2, REALONE, REALZERO
+      PARAMETER          ( PIOVER2 = 1.57079632679489662D0,
+     $                     REALONE = 1.0D0, REALZERO = 0.0D0 )
+      COMPLEX*16         ZERO, ONE
+      PARAMETER          ( ZERO = (0.0D0,0.0D0), ONE = (1.0D0,0.0D0) )
+*     ..
+*     .. Local Scalars ..
+      INTEGER            I, INFO, R
+      DOUBLE PRECISION   EPS2, RESID, ULP, ULPINV
+*     ..
+*     .. External Functions ..
+      DOUBLE PRECISION   DLAMCH, ZLANGE, ZLANHE
+      EXTERNAL           DLAMCH, ZLANGE, ZLANHE
+*     ..
+*     .. External Subroutines ..
+      EXTERNAL           ZGEMM, ZLACPY, ZLASET, ZUNCSD, ZHERK
+*     ..
+*     .. Intrinsic Functions ..
+      INTRINSIC          REAL, MAX, MIN
+*     ..
+*     .. Executable Statements ..
+*
+      ULP = DLAMCH( 'Precision' )
+      ULPINV = REALONE / ULP
+      CALL ZLASET( 'Full', M, M, ZERO, ONE, WORK, LDX )
+      CALL ZHERK( 'Upper', 'Conjugate transpose', M, M, -ONE, X, LDX,
+     $            ONE, WORK, LDX )
+      EPS2 = MAX( ULP, 
+     $            ZLANGE( '1', M, M, WORK, LDX, RWORK ) / REAL( M ) )
+      R = MIN( P, M-P, Q, M-Q )
+*
+*     Copy the matrix X to the array XF.
+*
+      CALL ZLACPY( 'Full', M, M, X, LDX, XF, LDX )
+*
+*     Compute the CSD
+*
+      CALL ZUNCSD( 'Y', 'Y', 'Y', 'Y', 'N', 'D', M, P, Q, XF(1,1), LDX,
+     $             XF(1,Q+1), LDX, XF(P+1,1), LDX, XF(P+1,Q+1), LDX,
+     $             THETA, U1, LDU1, U2, LDU2, V1T, LDV1T, V2T, LDV2T,
+     $             WORK, LWORK, RWORK, 17*(R+2), IWORK, INFO )
+*
+*     Compute X := diag(U1,U2)'*X*diag(V1,V2) - [D11 D12; D21 D22]
+*
+      CALL ZGEMM( 'No transpose', 'Conjugate transpose', P, Q, Q, ONE,
+     $            X, LDX, V1T, LDV1T, ZERO, WORK, LDX )
+*
+      CALL ZGEMM( 'Conjugate transpose', 'No transpose', P, Q, P, ONE,
+     $            U1, LDU1, WORK, LDX, ZERO, X, LDX )
+*
+      DO I = 1, MIN(P,Q)-R
+         X(I,I) = X(I,I) - ONE
+      END DO
+      DO I = 1, R
+         X(MIN(P,Q)-R+I,MIN(P,Q)-R+I) =
+     $           X(MIN(P,Q)-R+I,MIN(P,Q)-R+I) - COMPLEX( COS(THETA(I)),
+     $              0.0D0 )
+      END DO
+*
+      CALL ZGEMM( 'No transpose', 'Conjugate transpose', P, M-Q, M-Q,
+     $            ONE, X(1,Q+1), LDX, V2T, LDV2T, ZERO, WORK, LDX )
+*
+      CALL ZGEMM( 'Conjugate transpose', 'No transpose', P, M-Q, P,
+     $            ONE, U1, LDU1, WORK, LDX, ZERO, X(1,Q+1), LDX )
+*
+      DO I = 1, MIN(P,M-Q)-R
+         X(P-I+1,M-I+1) = X(P-I+1,M-I+1) + ONE
+      END DO
+      DO I = 1, R
+         X(P-(MIN(P,M-Q)-R)+1-I,M-(MIN(P,M-Q)-R)+1-I) =
+     $      X(P-(MIN(P,M-Q)-R)+1-I,M-(MIN(P,M-Q)-R)+1-I) +
+     $      COMPLEX( SIN(THETA(R-I+1)), 0.0D0 )
+      END DO
+*
+      CALL ZGEMM( 'No transpose', 'Conjugate transpose', M-P, Q, Q, ONE,
+     $            X(P+1,1), LDX, V1T, LDV1T, ZERO, WORK, LDX )
+*
+      CALL ZGEMM( 'Conjugate transpose', 'No transpose', M-P, Q, M-P,
+     $            ONE, U2, LDU2, WORK, LDX, ZERO, X(P+1,1), LDX )
+*
+      DO I = 1, MIN(M-P,Q)-R
+         X(M-I+1,Q-I+1) = X(M-I+1,Q-I+1) - ONE
+      END DO
+      DO I = 1, R
+         X(M-(MIN(M-P,Q)-R)+1-I,Q-(MIN(M-P,Q)-R)+1-I) =
+     $             X(M-(MIN(M-P,Q)-R)+1-I,Q-(MIN(M-P,Q)-R)+1-I) -
+     $             COMPLEX( SIN(THETA(R-I+1)), 0.0D0 )
+      END DO
+*
+      CALL ZGEMM( 'No transpose', 'Conjugate transpose', M-P, M-Q, M-Q,
+     $            ONE, X(P+1,Q+1), LDX, V2T, LDV2T, ZERO, WORK, LDX )
+*
+      CALL ZGEMM( 'Conjugate transpose', 'No transpose', M-P, M-Q, M-P,
+     $            ONE, U2, LDU2, WORK, LDX, ZERO, X(P+1,Q+1), LDX )
+*
+      DO I = 1, MIN(M-P,M-Q)-R
+         X(P+I,Q+I) = X(P+I,Q+I) - ONE
+      END DO
+      DO I = 1, R
+         X(P+(MIN(M-P,M-Q)-R)+I,Q+(MIN(M-P,M-Q)-R)+I) =
+     $      X(P+(MIN(M-P,M-Q)-R)+I,Q+(MIN(M-P,M-Q)-R)+I) -
+     $      COMPLEX( COS(THETA(I)), 0.0D0 )
+      END DO
+*
+*     Compute norm( U1'*X11*V1 - D11 ) / ( MAX(1,P,Q)*EPS2 ) .
+*
+      RESID = ZLANGE( '1', P, Q, X, LDX, RWORK )
+      RESULT( 1 ) = ( RESID / REAL(MAX(1,P,Q)) ) / EPS2
+*
+*     Compute norm( U1'*X12*V2 - D12 ) / ( MAX(1,P,M-Q)*EPS2 ) .
+*
+      RESID = ZLANGE( '1', P, M-Q, X(1,Q+1), LDX, RWORK )
+      RESULT( 2 ) = ( RESID / REAL(MAX(1,P,M-Q)) ) / EPS2
+*
+*     Compute norm( U2'*X21*V1 - D21 ) / ( MAX(1,M-P,Q)*EPS2 ) .
+*
+      RESID = ZLANGE( '1', M-P, Q, X(P+1,1), LDX, RWORK )
+      RESULT( 3 ) = ( RESID / REAL(MAX(1,M-P,Q)) ) / EPS2
+*
+*     Compute norm( U2'*X22*V2 - D22 ) / ( MAX(1,M-P,M-Q)*EPS2 ) .
+*
+      RESID = ZLANGE( '1', M-P, M-Q, X(P+1,Q+1), LDX, RWORK )
+      RESULT( 4 ) = ( RESID / REAL(MAX(1,M-P,M-Q)) ) / EPS2
+*
+*     Compute I - U1'*U1
+*
+      CALL ZLASET( 'Full', P, P, ZERO, ONE, WORK, LDU1 )
+      CALL ZHERK( 'Upper', 'Conjugate transpose', P, P, -ONE, U1, LDU1,
+     $            ONE, WORK, LDU1 )
+*
+*     Compute norm( I - U'*U ) / ( MAX(1,P) * ULP ) .
+*
+      RESID = ZLANHE( '1', 'Upper', P, WORK, LDU1, RWORK )
+      RESULT( 5 ) = ( RESID / REAL(MAX(1,P)) ) / ULP
+*
+*     Compute I - U2'*U2
+*
+      CALL ZLASET( 'Full', M-P, M-P, ZERO, ONE, WORK, LDU2 )
+      CALL ZHERK( 'Upper', 'Conjugate transpose', M-P, M-P, -ONE, U2,
+     $            LDU2, ONE, WORK, LDU2 )
+*
+*     Compute norm( I - U2'*U2 ) / ( MAX(1,M-P) * ULP ) .
+*
+      RESID = ZLANHE( '1', 'Upper', M-P, WORK, LDU2, RWORK )
+      RESULT( 6 ) = ( RESID / REAL(MAX(1,M-P)) ) / ULP
+*
+*     Compute I - V1T*V1T'
+*
+      CALL ZLASET( 'Full', Q, Q, ZERO, ONE, WORK, LDV1T )
+      CALL ZHERK( 'Upper', 'No transpose', Q, Q, -ONE, V1T, LDV1T, ONE,
+     $            WORK, LDV1T )
+*
+*     Compute norm( I - V1T*V1T' ) / ( MAX(1,Q) * ULP ) .
+*
+      RESID = ZLANHE( '1', 'Upper', Q, WORK, LDV1T, RWORK )
+      RESULT( 7 ) = ( RESID / REAL(MAX(1,Q)) ) / ULP
+*
+*     Compute I - V2T*V2T'
+*
+      CALL ZLASET( 'Full', M-Q, M-Q, ZERO, ONE, WORK, LDV2T )
+      CALL ZHERK( 'Upper', 'No transpose', M-Q, M-Q, -ONE, V2T, LDV2T,
+     $            ONE, WORK, LDV2T )
+*
+*     Compute norm( I - V2T*V2T' ) / ( MAX(1,M-Q) * ULP ) .
+*
+      RESID = ZLANHE( '1', 'Upper', M-Q, WORK, LDV2T, RWORK )
+      RESULT( 8 ) = ( RESID / REAL(MAX(1,M-Q)) ) / ULP
+*
+*     Check sorting
+*
+      RESULT(9) = REALZERO
+      DO I = 1, R
+         IF( THETA(I).LT.REALZERO .OR. THETA(I).GT.PIOVER2 ) THEN
+            RESULT(9) = ULPINV
+         END IF
+         IF( I.GT.1 .AND. THETA(I).LT.THETA(I-1) ) THEN
+            RESULT(9) = ULPINV
+         END IF
+      END DO
+*
+      RETURN
+*      
+*     End of ZCSDTS
+*
+      END
+
diff --git a/TESTING/csd.in b/TESTING/csd.in
new file mode 100644
index 00000000..c2edcec5
--- /dev/null
+++ b/TESTING/csd.in
@@ -0,0 +1,10 @@
+CSD:  Data file for testing CS decomposition routines
+10                                   Number of values of M, P, Q
+0  10 10 10 10 21 24 30 22 32 55     Values of M (row and column dimension of unitary matrix)
+0  4   4 0  10 9  10 20 12 12 40     Values of P (row dimension of top-left block)
+0  0  10 4  4  15 12 8  20 8  20     Values of Q (column dimension of top-left block)
+10.0                                 Threshold value of test ratio
+T                                    Put T to test the error exits
+1                                    Code to interpret the seed
+CSD   3                              List types on next line if 0 < NTYPES < 4
+1 2 3