dcsdts.f

Go to the documentation of this file.

*> \brief \b DCSDTS
*
*  =========== DOCUMENTATION ===========
*
* Online html documentation available at
*            http://www.netlib.org/lapack/explore-html/
*
*  Definition:
*  ===========
*
*       SUBROUTINE DCSDTS( M, P, Q, X, XF, LDX, U1, LDU1, U2, LDU2, V1T,
*                          LDV1T, V2T, LDV2T, THETA, IWORK, WORK, LWORK,
*                          RWORK, RESULT )
*
*       .. Scalar Arguments ..
*       INTEGER            LDX, LDU1, LDU2, LDV1T, LDV2T, LWORK, M, P, Q
*       ..
*       .. Array Arguments ..
*       INTEGER            IWORK( * )
*       DOUBLE PRECISION   RESULT( 15 ), RWORK( * ), THETA( * )
*       DOUBLE PRECISION   U1( LDU1, * ), U2( LDU2, * ), V1T( LDV1T, * ),
*      $                   V2T( LDV2T, * ), WORK( LWORK ), X( LDX, * ),
*      $                   XF( LDX, * )
*       ..
*
*
*> \par Purpose:
*  =============
*>
*> \verbatim
*>
*> DCSDTS tests DORCSD, which, given an M-by-M partitioned orthogonal
*> matrix X,
*>              Q  M-Q
*>       X = [ X11 X12 ] P   ,
*>           [ X21 X22 ] M-P
*>
*> computes the CSD
*>
*>       [ U1    ]**T * [ X11 X12 ] * [ V1    ]
*>       [    U2 ]      [ X21 X22 ]   [    V2 ]
*>
*>                             [  I  0  0 |  0  0  0 ]
*>                             [  0  C  0 |  0 -S  0 ]
*>                             [  0  0  0 |  0  0 -I ]
*>                           = [---------------------] = [ D11 D12 ] ,
*>                             [  0  0  0 |  I  0  0 ]   [ D21 D22 ]
*>                             [  0  S  0 |  0  C  0 ]
*>                             [  0  0  I |  0  0  0 ]
*>
*> and also DORCSD2BY1, which, given
*>          Q
*>       [ X11 ] P   ,
*>       [ X21 ] M-P
*>
*> computes the 2-by-1 CSD
*>
*>                                     [  I  0  0 ]
*>                                     [  0  C  0 ]
*>                                     [  0  0  0 ]
*>       [ U1    ]**T * [ X11 ] * V1 = [----------] = [ D11 ] ,
*>       [    U2 ]      [ X21 ]        [  0  0  0 ]   [ D21 ]
*>                                     [  0  S  0 ]
*>                                     [  0  0  I ]
*> \endverbatim
*
*  Arguments:
*  ==========
*
*> \param[in] M
*> \verbatim
*>          M is INTEGER
*>          The number of rows of the matrix X.  M >= 0.
*> \endverbatim
*>
*> \param[in] P
*> \verbatim
*>          P is INTEGER
*>          The number of rows of the matrix X11.  P >= 0.
*> \endverbatim
*>
*> \param[in] Q
*> \verbatim
*>          Q is INTEGER
*>          The number of columns of the matrix X11.  Q >= 0.
*> \endverbatim
*>
*> \param[in] X
*> \verbatim
*>          X is DOUBLE PRECISION array, dimension (LDX,M)
*>          The M-by-M matrix X.
*> \endverbatim
*>
*> \param[out] XF
*> \verbatim
*>          XF is DOUBLE PRECISION array, dimension (LDX,M)
*>          Details of the CSD of X, as returned by DORCSD;
*>          see DORCSD for further details.
*> \endverbatim
*>
*> \param[in] LDX
*> \verbatim
*>          LDX is INTEGER
*>          The leading dimension of the arrays X and XF.
*>          LDX >= max( 1,M ).
*> \endverbatim
*>
*> \param[out] U1
*> \verbatim
*>          U1 is DOUBLE PRECISION array, dimension(LDU1,P)
*>          The P-by-P orthogonal matrix U1.
*> \endverbatim
*>
*> \param[in] LDU1
*> \verbatim
*>          LDU1 is INTEGER
*>          The leading dimension of the array U1. LDU >= max(1,P).
*> \endverbatim
*>
*> \param[out] U2
*> \verbatim
*>          U2 is DOUBLE PRECISION array, dimension(LDU2,M-P)
*>          The (M-P)-by-(M-P) orthogonal matrix U2.
*> \endverbatim
*>
*> \param[in] LDU2
*> \verbatim
*>          LDU2 is INTEGER
*>          The leading dimension of the array U2. LDU >= max(1,M-P).
*> \endverbatim
*>
*> \param[out] V1T
*> \verbatim
*>          V1T is DOUBLE PRECISION array, dimension(LDV1T,Q)
*>          The Q-by-Q orthogonal matrix V1T.
*> \endverbatim
*>
*> \param[in] LDV1T
*> \verbatim
*>          LDV1T is INTEGER
*>          The leading dimension of the array V1T. LDV1T >=
*>          max(1,Q).
*> \endverbatim
*>
*> \param[out] V2T
*> \verbatim
*>          V2T is DOUBLE PRECISION array, dimension(LDV2T,M-Q)
*>          The (M-Q)-by-(M-Q) orthogonal matrix V2T.
*> \endverbatim
*>
*> \param[in] LDV2T
*> \verbatim
*>          LDV2T is INTEGER
*>          The leading dimension of the array V2T. LDV2T >=
*>          max(1,M-Q).
*> \endverbatim
*>
*> \param[out] THETA
*> \verbatim
*>          THETA is DOUBLE PRECISION array, dimension MIN(P,M-P,Q,M-Q)
*>          The CS values of X; the essentially diagonal matrices C and
*>          S are constructed from THETA; see subroutine DORCSD for
*>          details.
*> \endverbatim
*>
*> \param[out] IWORK
*> \verbatim
*>          IWORK is INTEGER array, dimension (M)
*> \endverbatim
*>
*> \param[out] WORK
*> \verbatim
*>          WORK is DOUBLE PRECISION array, dimension (LWORK)
*> \endverbatim
*>
*> \param[in] LWORK
*> \verbatim
*>          LWORK is INTEGER
*>          The dimension of the array WORK
*> \endverbatim
*>
*> \param[out] RWORK
*> \verbatim
*>          RWORK is DOUBLE PRECISION array
*> \endverbatim
*>
*> \param[out] RESULT
*> \verbatim
*>          RESULT is DOUBLE PRECISION array, dimension (15)
*>          The test ratios:
*>          First, the 2-by-2 CSD:
*>          RESULT(1) = norm( U1'*X11*V1 - D11 ) / ( MAX(1,P,Q)*EPS2 )
*>          RESULT(2) = norm( U1'*X12*V2 - D12 ) / ( MAX(1,P,M-Q)*EPS2 )
*>          RESULT(3) = norm( U2'*X21*V1 - D21 ) / ( MAX(1,M-P,Q)*EPS2 )
*>          RESULT(4) = norm( U2'*X22*V2 - D22 ) / ( MAX(1,M-P,M-Q)*EPS2 )
*>          RESULT(5) = norm( I - U1'*U1 ) / ( MAX(1,P)*ULP )
*>          RESULT(6) = norm( I - U2'*U2 ) / ( MAX(1,M-P)*ULP )
*>          RESULT(7) = norm( I - V1T'*V1T ) / ( MAX(1,Q)*ULP )
*>          RESULT(8) = norm( I - V2T'*V2T ) / ( MAX(1,M-Q)*ULP )
*>          RESULT(9) = 0        if THETA is in increasing order and
*>                               all angles are in [0,pi/2];
*>                    = ULPINV   otherwise.
*>          Then, the 2-by-1 CSD:
*>          RESULT(10) = norm( U1'*X11*V1 - D11 ) / ( MAX(1,P,Q)*EPS2 )
*>          RESULT(11) = norm( U2'*X21*V1 - D21 ) / ( MAX(1,M-P,Q)*EPS2 )
*>          RESULT(12) = norm( I - U1'*U1 ) / ( MAX(1,P)*ULP )
*>          RESULT(13) = norm( I - U2'*U2 ) / ( MAX(1,M-P)*ULP )
*>          RESULT(14) = norm( I - V1T'*V1T ) / ( MAX(1,Q)*ULP )
*>          RESULT(15) = 0        if THETA is in increasing order and
*>                                all angles are in [0,pi/2];
*>                     = ULPINV   otherwise.
*>          ( EPS2 = MAX( norm( I - X'*X ) / M, ULP ). )
*> \endverbatim
*
*  Authors:
*  ========
*
*> \author Univ. of Tennessee
*> \author Univ. of California Berkeley
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
*> \ingroup double_eig
*
*  =====================================================================
      SUBROUTINE dcsdts( M, P, Q, X, XF, LDX, U1, LDU1, U2, LDU2, V1T,
     $                   LDV1T, V2T, LDV2T, THETA, IWORK, WORK, LWORK,
     $                   RWORK, RESULT )
*
*  -- LAPACK test routine --
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*
*     .. Scalar Arguments ..
      INTEGER            LDX, LDU1, LDU2, LDV1T, LDV2T, LWORK, M, P, Q
*     ..
*     .. Array Arguments ..
      INTEGER            IWORK( * )
      DOUBLE PRECISION   RESULT( 15 ), RWORK( * ), THETA( * )
      DOUBLE PRECISION   U1( LDU1, * ), U2( LDU2, * ), V1T( LDV1T, * ),
     $                   v2t( ldv2t, * ), work( lwork ), x( ldx, * ),
     $                   xf( ldx, * )
*     ..
*
*  =====================================================================
*
*     .. Parameters ..
      DOUBLE PRECISION   REALONE, REALZERO
      PARAMETER          ( REALONE = 1.0d0, realzero = 0.0d0 )
      DOUBLE PRECISION   ZERO, ONE
      parameter( zero = 0.0d0, one = 1.0d0 )
      DOUBLE PRECISION   PIOVER2
      parameter( piover2 = 1.57079632679489661923132169163975144210d0 )
*     ..
*     .. 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, dorcsd2by1,
     $                   dsyrk
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          cos, dble, max, min, sin
*     ..
*     .. Executable Statements ..
*
      ulp = dlamch( 'Precision' )
      ulpinv = realone / ulp
*
*     The first half of the routine checks the 2-by-2 CSD
*
      CALL dlaset( 'Full', m, m, zero, one, work, ldx )
      CALL dsyrk( 'upper', 'conjugate transpose', M, M, -ONE, X, LDX,
     $            ONE, WORK, LDX )
.GT.      IF (M0) THEN
         EPS2 = MAX( ULP,
     $               DLANGE( '1', M, M, WORK, LDX, RWORK ) / DBLE( M ) )
      ELSE
         EPS2 = ULP
      END IF
      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 XF := diag(U1,U2)'*X*diag(V1,V2) - [D11 D12; D21 D22]
*
      CALL dlacpy( 'Full', m, m, x, ldx, xf, ldx )
*
      CALL dgemm( 'No transpose', 'Conjugate transpose', p, q, q, one,
     $            xf, ldx, v1t, ldv1t, zero, work, ldx )
*
      CALL dgemm( 'Conjugate transpose', 'No transpose', p, q, p, one,
     $            u1, ldu1, work, ldx, zero, xf, ldx )
*
      DO i = 1, min(p,q)-r
         xf(i,i) = xf(i,i) - one
      END DO
      DO i = 1, r
         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 dgemm( 'No transpose', 'Conjugate transpose', p, m-q, m-q,
     $            one, xf(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, xf(1,q+1), ldx )
*
      DO i = 1, min(p,m-q)-r
         xf(p-i+1,m-i+1) = xf(p-i+1,m-i+1) + one
      END DO
      DO i = 1, r
         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 dgemm( 'No transpose', 'Conjugate transpose', m-p, q, q, one,
     $            xf(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, xf(p+1,1), ldx )
*
      DO i = 1, min(m-p,q)-r
         xf(m-i+1,q-i+1) = xf(m-i+1,q-i+1) - one
      END DO
      DO i = 1, r
         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 dgemm( 'No transpose', 'conjugate transpose', M-P, M-Q, M-Q,
     $            ONE, XF(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, XF(P+1,Q+1), LDX )
*
      DO I = 1, MIN(M-P,M-Q)-R
         XF(P+I,Q+I) = XF(P+I,Q+I) - ONE
      END DO
      DO I = 1, R
         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 = DLANGE( '1', P, Q, XF, 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, XF(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, XF(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, XF(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 ) 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 dlaset( 'Full', q, q, zero, one, work, ldx )
      CALL dsyrk( 'Upper', 'Conjugate transpose', q, m, -one, x, ldx,
     $            one, work, ldx )
      IF( m.GT.0 ) THEN
         eps2 = max( ulp,
     $               dlange( '1', q, q, work, ldx, rwork ) / dble( m ) )
      ELSE
         eps2 = ulp
      END IF
      r = min( p, m-p, q, m-q )
*
*     Copy the matrix [ X11; X21 ] to the array XF.
*
      CALL dlacpy( 'Full', m, q, x, ldx, xf, ldx )
*
*     Compute the CSD
*
      CALL dorcsd2by1( '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 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', 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
*
*     Compute norm( U1'*X11*V1 - D11 ) / ( MAX(1,P,Q)*EPS2 ) .
*
      resid = dlange( '1', p, q, x, ldx, rwork )
      result( 10 ) = ( resid / dble(max(1,p,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( 11 ) = ( resid / dble(max(1,m-p,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 - U1'*U1 ) / ( MAX(1,P) * ULP ) .
*
      resid = dlansy( '1', 'Upper', p, work, ldu1, rwork )
      result( 12 ) = ( 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( 13 ) = ( 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( 14 ) = ( RESID / DBLE(MAX(1,Q)) ) / ULP
*
*     Check sorting
*
      RESULT( 15 ) = REALZERO
      DO I = 1, R
.LT..OR..GT.         IF( THETA(I)REALZERO  THETA(I)PIOVER2 ) THEN
            RESULT( 15 ) = ULPINV
         END IF
.GT.         IF( I1 ) THEN
.LT.            IF ( THETA(I)THETA(I-1) ) THEN
               RESULT( 15 ) = ULPINV
            END IF
         END IF
      END DO
*
      RETURN
*
*     End of DCSDTS
*
      END