Commit Brian Sutton new CS Decomposition routines.

All the routines from the SRC folder have been updated to integrate the current Doxygen layout.
Everything seems to be fine, all tests passed without problem.

1 files changed, 176 insertions, 34 deletions

diff --git a/TESTING/EIG/scsdts.f b/TESTING/EIG/scsdts.f
index 74b32ead..a326f356 100644
--- a/TESTING/EIG/scsdts.f
+++ b/TESTING/EIG/scsdts.f
@@ -17,7 +17,7 @@
 *       ..
 *       .. Array Arguments ..
 *       INTEGER            IWORK( * )
-*       REAL               RESULT( 9 ), RWORK( * ), THETA( * )
+*       REAL               RESULT( 15 ), RWORK( * ), THETA( * )
 *       REAL               U1( LDU1, * ), U2( LDU2, * ), V1T( LDV1T, * ),
 *      $                   V2T( LDV2T, * ), WORK( LWORK ), X( LDX, * ),
 *      $                   XF( LDX, * )
@@ -47,6 +47,21 @@
 *>                             [  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 ]
 *> \endverbatim
 *
 *  Arguments:
@@ -171,8 +186,9 @@
 *>
 *> \param[out] RESULT
 *> \verbatim
-*>          RESULT is REAL array, dimension (9)
+*>          RESULT is REAL array, dimension (15)
 *>          The test ratios:
+*>          First, the 2-by-2 CSD:
 *>          RESULT(1) = norm( U1'*X11*V1 - D11 ) / ( MAX(1,P,Q)*EPS2 )
 *>          RESULT(2) = norm( U1'*X12*V2 - D12 ) / ( MAX(1,P,M-Q)*EPS2 )
 *>          RESULT(3) = norm( U2'*X21*V1 - D21 ) / ( MAX(1,M-P,Q)*EPS2 )
@@ -184,6 +200,15 @@
 *>          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 ). )
 *> \endverbatim
 *
@@ -214,7 +239,7 @@
 *     ..
 *     .. Array Arguments ..
       INTEGER            IWORK( * )
-      REAL               RESULT( 9 ), RWORK( * ), THETA( * )
+      REAL               RESULT( 15 ), RWORK( * ), THETA( * )
       REAL               U1( LDU1, * ), U2( LDU2, * ), V1T( LDV1T, * ),
      $                   V2T( LDV2T, * ), WORK( LWORK ), X( LDX, * ),
      $                   XF( LDX, * )
@@ -238,15 +263,19 @@
       EXTERNAL           SLAMCH, SLANGE, SLANSY
 *     ..
 *     .. External Subroutines ..
-      EXTERNAL           SGEMM, SLACPY, SLASET, SORCSD, SSYRK
+      EXTERNAL           SGEMM, SLACPY, SLASET, SORCSD, SORCSD2BY1,
+     $                   SSYRK
 *     ..
 *     .. Intrinsic Functions ..
-      INTRINSIC          REAL, MAX, MIN
+      INTRINSIC          COS, MAX, MIN, REAL, SIN
 *     ..
 *     .. Executable Statements ..
 *
       ULP = SLAMCH( 'Precision' )
       ULPINV = REALONE / ULP
+*
+*     The first half of the routine checks the 2-by-2 CSD
+*
       CALL SLASET( 'Full', M, M, ZERO, ONE, WORK, LDX )
       CALL SSYRK( 'Upper', 'Conjugate transpose', M, M, -ONE, X, LDX,
      $            ONE, WORK, LDX )
@@ -269,85 +298,87 @@
      $             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]
+*     Compute XF := diag(U1,U2)'*X*diag(V1,V2) - [D11 D12; D21 D22]
+*
+      CALL SLACPY( 'Full', M, M, X, LDX, XF, LDX )
 *
       CALL SGEMM( 'No transpose', 'Conjugate transpose', P, Q, Q, ONE,
-     $            X, LDX, V1T, LDV1T, ZERO, WORK, LDX )
+     $            XF, LDX, V1T, LDV1T, ZERO, WORK, LDX )
 *
       CALL SGEMM( 'Conjugate transpose', 'No transpose', P, Q, P, ONE,
-     $            U1, LDU1, WORK, LDX, ZERO, X, LDX )
+     $            U1, LDU1, WORK, LDX, ZERO, XF, LDX )
 *
       DO I = 1, MIN(P,Q)-R
-         X(I,I) = X(I,I) - ONE
+         XF(I,I) = XF(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))
+         XF(MIN(P,Q)-R+I,MIN(P,Q)-R+I) =
+     $           XF(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 )
+     $            ONE, XF(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 )
+     $            ONE, U1, LDU1, WORK, LDX, ZERO, XF(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
+         XF(P-I+1,M-I+1) = XF(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) +
+         XF(P-(MIN(P,M-Q)-R)+1-I,M-(MIN(P,M-Q)-R)+1-I) =
+     $      XF(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 )
+     $            XF(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 )
+     $            ONE, U2, LDU2, WORK, LDX, ZERO, XF(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
+         XF(M-I+1,Q-I+1) = XF(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) -
+         XF(M-(MIN(M-P,Q)-R)+1-I,Q-(MIN(M-P,Q)-R)+1-I) =
+     $             XF(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 )
+     $            ONE, XF(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 )
+     $            ONE, U2, LDU2, WORK, LDX, ZERO, XF(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
+         XF(P+I,Q+I) = XF(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) -
+         XF(P+(MIN(M-P,M-Q)-R)+I,Q+(MIN(M-P,M-Q)-R)+I) =
+     $      XF(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 )
+      RESID = SLANGE( '1', P, Q, XF, 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 )
+      RESID = SLANGE( '1', P, M-Q, XF(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 )
+      RESID = SLANGE( '1', M-P, Q, XF(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 )
+      RESID = SLANGE( '1', M-P, M-Q, XF(P+1,Q+1), LDX, RWORK )
       RESULT( 4 ) = ( RESID / REAL(MAX(1,M-P,M-Q)) ) / EPS2
 *
 *     Compute I - U1'*U1
@@ -396,14 +427,125 @@
 *
 *     Check sorting
 *
-      RESULT(9) = REALZERO
+      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 ) THEN
+            IF ( THETA(I).LT.THETA(I-1) ) THEN
+               RESULT( 9 ) = ULPINV
+            END IF
+         END IF
+      END DO
+*
+*     The second half of the routine checks the 2-by-1 CSD
+*
+      CALL SLASET( 'Full', Q, Q, ZERO, ONE, WORK, LDX )
+      CALL SSYRK( 'Upper', 'Conjugate transpose', Q, M, -ONE, X, LDX,
+     $            ONE, WORK, LDX )
+      IF (M.GT.0) THEN
+         EPS2 = MAX( ULP, 
+     $               SLANGE( '1', Q, Q, WORK, LDX, RWORK ) / REAL( M ) )
+      ELSE
+         EPS2 = ULP
+      END IF
+      R = MIN( P, M-P, Q, M-Q )
+*
+*     Copy the matrix [X11;X21] to the array XF.
+*
+      CALL SLACPY( 'Full', M, Q, X, LDX, XF, LDX )
+*
+*     Compute the CSD
+*
+      CALL SORCSD2BY1( 'Y', 'Y', 'Y', M, P, Q, XF(1,1), LDX, XF(P+1,1),
+     $                 LDX, THETA, U1, LDU1, U2, LDU2, V1T, LDV1T, WORK,
+     $                 LWORK, IWORK, INFO )
+*
+*     Compute [X11;X21] := diag(U1,U2)'*[X11;X21]*V1 - [D11;D21]
+*
+      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', 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
+*
+*     Compute norm( U1'*X11*V1 - D11 ) / ( MAX(1,P,Q)*EPS2 ) .
+*
+      RESID = SLANGE( '1', P, Q, X, LDX, RWORK )
+      RESULT( 10 ) = ( RESID / REAL(MAX(1,P,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( 11 ) = ( RESID / REAL(MAX(1,M-P,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 - U1'*U1 ) / ( MAX(1,P) * ULP ) .
+*
+      RESID = SLANSY( '1', 'Upper', P, WORK, LDU1, RWORK )
+      RESULT( 12 ) = ( 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( 13 ) = ( 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( 14 ) = ( RESID / REAL(MAX(1,Q)) ) / ULP
+*
+*     Check sorting
+*
+      RESULT( 15 ) = REALZERO
       DO I = 1, R
          IF( THETA(I).LT.REALZERO .OR. THETA(I).GT.PIOVER2 ) THEN
-            RESULT(9) = ULPINV
+            RESULT( 15 ) = ULPINV
          END IF
-         IF( I.GT.1) THEN
+         IF( I.GT.1 ) THEN
             IF ( THETA(I).LT.THETA(I-1) ) THEN
-               RESULT(9) = ULPINV
+               RESULT( 15 ) = ULPINV
             END IF
          END IF
       END DO