dgsvts.f

Go to the documentation of this file.

*> \brief \b DGSVTS
*
*  =========== DOCUMENTATION ===========
*
* Online html documentation available at 
*            http://www.netlib.org/lapack/explore-html/ 
*
*  Definition:
*  ===========
*
*       SUBROUTINE DGSVTS( M, P, N, A, AF, LDA, B, BF, LDB, U, LDU, V,
*                          LDV, Q, LDQ, ALPHA, BETA, R, LDR, IWORK, WORK,
*                          LWORK, RWORK, RESULT )
* 
*       .. Scalar Arguments ..
*       INTEGER            LDA, LDB, LDQ, LDR, LDU, LDV, LWORK, M, N, P
*       ..
*       .. Array Arguments ..
*       INTEGER            IWORK( * )
*       DOUBLE PRECISION   A( LDA, * ), AF( LDA, * ), ALPHA( * ),
*      $                   B( LDB, * ), BETA( * ), BF( LDB, * ),
*      $                   Q( LDQ, * ), R( LDR, * ), RESULT( 6 ),
*      $                   RWORK( * ), U( LDU, * ), V( LDV, * ),
*      $                   WORK( LWORK )
*       ..
*  
*
*> \par Purpose:
*  =============
*>
*> \verbatim
*>
*> DGSVTS tests DGGSVD, which computes the GSVD of an M-by-N matrix A
*> and a P-by-N matrix B:
*>              U'*A*Q = D1*R and V'*B*Q = D2*R.
*> \endverbatim
*
*  Arguments:
*  ==========
*
*> \param[in] M
*> \verbatim
*>          M is INTEGER
*>          The number of rows of the matrix A.  M >= 0.
*> \endverbatim
*>
*> \param[in] P
*> \verbatim
*>          P is INTEGER
*>          The number of rows of the matrix B.  P >= 0.
*> \endverbatim
*>
*> \param[in] N
*> \verbatim
*>          N is INTEGER
*>          The number of columns of the matrices A and B.  N >= 0.
*> \endverbatim
*>
*> \param[in] A
*> \verbatim
*>          A is DOUBLE PRECISION array, dimension (LDA,M)
*>          The M-by-N matrix A.
*> \endverbatim
*>
*> \param[out] AF
*> \verbatim
*>          AF is DOUBLE PRECISION array, dimension (LDA,N)
*>          Details of the GSVD of A and B, as returned by DGGSVD,
*>          see DGGSVD for further details.
*> \endverbatim
*>
*> \param[in] LDA
*> \verbatim
*>          LDA is INTEGER
*>          The leading dimension of the arrays A and AF.
*>          LDA >= max( 1,M ).
*> \endverbatim
*>
*> \param[in] B
*> \verbatim
*>          B is DOUBLE PRECISION array, dimension (LDB,P)
*>          On entry, the P-by-N matrix B.
*> \endverbatim
*>
*> \param[out] BF
*> \verbatim
*>          BF is DOUBLE PRECISION array, dimension (LDB,N)
*>          Details of the GSVD of A and B, as returned by DGGSVD,
*>          see DGGSVD for further details.
*> \endverbatim
*>
*> \param[in] LDB
*> \verbatim
*>          LDB is INTEGER
*>          The leading dimension of the arrays B and BF.
*>          LDB >= max(1,P).
*> \endverbatim
*>
*> \param[out] U
*> \verbatim
*>          U is DOUBLE PRECISION array, dimension(LDU,M)
*>          The M by M orthogonal matrix U.
*> \endverbatim
*>
*> \param[in] LDU
*> \verbatim
*>          LDU is INTEGER
*>          The leading dimension of the array U. LDU >= max(1,M).
*> \endverbatim
*>
*> \param[out] V
*> \verbatim
*>          V is DOUBLE PRECISION array, dimension(LDV,M)
*>          The P by P orthogonal matrix V.
*> \endverbatim
*>
*> \param[in] LDV
*> \verbatim
*>          LDV is INTEGER
*>          The leading dimension of the array V. LDV >= max(1,P).
*> \endverbatim
*>
*> \param[out] Q
*> \verbatim
*>          Q is DOUBLE PRECISION array, dimension(LDQ,N)
*>          The N by N orthogonal matrix Q.
*> \endverbatim
*>
*> \param[in] LDQ
*> \verbatim
*>          LDQ is INTEGER
*>          The leading dimension of the array Q. LDQ >= max(1,N).
*> \endverbatim
*>
*> \param[out] ALPHA
*> \verbatim
*>          ALPHA is DOUBLE PRECISION array, dimension (N)
*> \endverbatim
*>
*> \param[out] BETA
*> \verbatim
*>          BETA is DOUBLE PRECISION array, dimension (N)
*>
*>          The generalized singular value pairs of A and B, the
*>          ``diagonal'' matrices D1 and D2 are constructed from
*>          ALPHA and BETA, see subroutine DGGSVD for details.
*> \endverbatim
*>
*> \param[out] R
*> \verbatim
*>          R is DOUBLE PRECISION array, dimension(LDQ,N)
*>          The upper triangular matrix R.
*> \endverbatim
*>
*> \param[in] LDR
*> \verbatim
*>          LDR is INTEGER
*>          The leading dimension of the array R. LDR >= max(1,N).
*> \endverbatim
*>
*> \param[out] IWORK
*> \verbatim
*>          IWORK is INTEGER array, dimension (N)
*> \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,
*>          LWORK >= max(M,P,N)*max(M,P,N).
*> \endverbatim
*>
*> \param[out] RWORK
*> \verbatim
*>          RWORK is DOUBLE PRECISION array, dimension (max(M,P,N))
*> \endverbatim
*>
*> \param[out] RESULT
*> \verbatim
*>          RESULT is DOUBLE PRECISION array, dimension (6)
*>          The test ratios:
*>          RESULT(1) = norm( U'*A*Q - D1*R ) / ( MAX(M,N)*norm(A)*ULP)
*>          RESULT(2) = norm( V'*B*Q - D2*R ) / ( MAX(P,N)*norm(B)*ULP)
*>          RESULT(3) = norm( I - U'*U ) / ( M*ULP )
*>          RESULT(4) = norm( I - V'*V ) / ( P*ULP )
*>          RESULT(5) = norm( I - Q'*Q ) / ( N*ULP )
*>          RESULT(6) = 0        if ALPHA is in decreasing order;
*>                    = ULPINV   otherwise.
*> \endverbatim
*
*  Authors:
*  ========
*
*> \author Univ. of Tennessee 
*> \author Univ. of California Berkeley 
*> \author Univ. of Colorado Denver 
*> \author NAG Ltd. 
*
*> \date November 2011
*
*> \ingroup double_eig
*
*  =====================================================================
      SUBROUTINE dgsvts( M, P, N, A, AF, LDA, B, BF, LDB, U, LDU, V,
     $                   ldv, q, ldq, alpha, beta, r, ldr, iwork, work,
     $                   lwork, rwork, result )
*
*  -- LAPACK test routine (version 3.4.0) --
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*     November 2011
*
*     .. Scalar Arguments ..
      INTEGER            lda, ldb, ldq, ldr, ldu, ldv, lwork, m, n, p
*     ..
*     .. Array Arguments ..
      INTEGER            iwork( * )
      DOUBLE PRECISION   a( lda, * ), af( lda, * ), alpha( * ),
     $                   b( ldb, * ), beta( * ), bf( ldb, * ),
     $                   q( ldq, * ), r( ldr, * ), result( 6 ),
     $                   rwork( * ), u( ldu, * ), v( ldv, * ),
     $                   work( lwork )
*     ..
*
*  =====================================================================
*
*     .. Parameters ..
      DOUBLE PRECISION   zero, one
      parameter( zero = 0.0d+0, one = 1.0d+0 )
*     ..
*     .. Local Scalars ..
      INTEGER            i, info, j, k, l
      DOUBLE PRECISION   anorm, bnorm, resid, temp, ulp, ulpinv, unfl
*     ..
*     .. External Functions ..
      DOUBLE PRECISION   dlamch, dlange, dlansy
      EXTERNAL           dlamch, dlange, dlansy
*     ..
*     .. External Subroutines ..
      EXTERNAL           dcopy, dgemm, dggsvd, dlacpy, dlaset, dsyrk
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          dble, max, min
*     ..
*     .. Executable Statements ..
*
      ulp = dlamch( 'Precision' )
      ulpinv = one / ulp
      unfl = dlamch( 'Safe minimum' )
*
*     Copy the matrix A to the array AF.
*
      CALL dlacpy( 'Full', m, n, a, lda, af, lda )
      CALL dlacpy( 'Full', p, n, b, ldb, bf, ldb )
*
      anorm = max( dlange( '1', m, n, a, lda, rwork ), unfl )
      bnorm = max( dlange( '1', p, n, b, ldb, rwork ), unfl )
*
*     Factorize the matrices A and B in the arrays AF and BF.
*
      CALL dggsvd( 'U', 'V', 'Q', m, n, p, k, l, af, lda, bf, ldb,
     $             alpha, beta, u, ldu, v, ldv, q, ldq, work, iwork,
     $             info )
*
*     Copy R
*
      DO 20 i = 1, min( k+l, m )
         DO 10 j = i, k + l
            r( i, j ) = af( i, n-k-l+j )
   10    continue
   20 continue
*
      IF( m-k-l.LT.0 ) THEN
         DO 40 i = m + 1, k + l
            DO 30 j = i, k + l
               r( i, j ) = bf( i-k, n-k-l+j )
   30       continue
   40    continue
      END IF
*
*     Compute A:= U'*A*Q - D1*R
*
      CALL dgemm( 'No transpose', 'No transpose', m, n, n, one, a, lda,
     $            q, ldq, zero, work, lda )
*
      CALL dgemm( 'Transpose', 'No transpose', m, n, m, one, u, ldu,
     $            work, lda, zero, a, lda )
*
      DO 60 i = 1, k
         DO 50 j = i, k + l
            a( i, n-k-l+j ) = a( i, n-k-l+j ) - r( i, j )
   50    continue
   60 continue
*
      DO 80 i = k + 1, min( k+l, m )
         DO 70 j = i, k + l
            a( i, n-k-l+j ) = a( i, n-k-l+j ) - alpha( i )*r( i, j )
   70    continue
   80 continue
*
*     Compute norm( U'*A*Q - D1*R ) / ( MAX(1,M,N)*norm(A)*ULP ) .
*
      resid = dlange( '1', m, n, a, lda, rwork )
*
      IF( anorm.GT.zero ) THEN
         result( 1 ) = ( ( resid / dble( max( 1, m, n ) ) ) / anorm ) /
     $                 ulp
      ELSE
         result( 1 ) = zero
      END IF
*
*     Compute B := V'*B*Q - D2*R
*
      CALL dgemm( 'No transpose', 'No transpose', p, n, n, one, b, ldb,
     $            q, ldq, zero, work, ldb )
*
      CALL dgemm( 'Transpose', 'No transpose', p, n, p, one, v, ldv,
     $            work, ldb, zero, b, ldb )
*
      DO 100 i = 1, l
         DO 90 j = i, l
            b( i, n-l+j ) = b( i, n-l+j ) - beta( k+i )*r( k+i, k+j )
   90    continue
  100 continue
*
*     Compute norm( V'*B*Q - D2*R ) / ( MAX(P,N)*norm(B)*ULP ) .
*
      resid = dlange( '1', p, n, b, ldb, rwork )
      IF( bnorm.GT.zero ) THEN
         result( 2 ) = ( ( resid / dble( max( 1, p, n ) ) ) / bnorm ) /
     $                 ulp
      ELSE
         result( 2 ) = zero
      END IF
*
*     Compute I - U'*U
*
      CALL dlaset( 'Full', m, m, zero, one, work, ldq )
      CALL dsyrk( 'Upper', 'Transpose', m, m, -one, u, ldu, one, work,
     $            ldu )
*
*     Compute norm( I - U'*U ) / ( M * ULP ) .
*
      resid = dlansy( '1', 'Upper', m, work, ldu, rwork )
      result( 3 ) = ( resid / dble( max( 1, m ) ) ) / ulp
*
*     Compute I - V'*V
*
      CALL dlaset( 'Full', p, p, zero, one, work, ldv )
      CALL dsyrk( 'Upper', 'Transpose', p, p, -one, v, ldv, one, work,
     $            ldv )
*
*     Compute norm( I - V'*V ) / ( P * ULP ) .
*
      resid = dlansy( '1', 'Upper', p, work, ldv, rwork )
      result( 4 ) = ( resid / dble( max( 1, p ) ) ) / ulp
*
*     Compute I - Q'*Q
*
      CALL dlaset( 'Full', n, n, zero, one, work, ldq )
      CALL dsyrk( 'Upper', 'Transpose', n, n, -one, q, ldq, one, work,
     $            ldq )
*
*     Compute norm( I - Q'*Q ) / ( N * ULP ) .
*
      resid = dlansy( '1', 'Upper', n, work, ldq, rwork )
      result( 5 ) = ( resid / dble( max( 1, n ) ) ) / ulp
*
*     Check sorting
*
      CALL dcopy( n, alpha, 1, work, 1 )
      DO 110 i = k + 1, min( k+l, m )
         j = iwork( i )
         IF( i.NE.j ) THEN
            temp = work( i )
            work( i ) = work( j )
            work( j ) = temp
         END IF
  110 continue
*
      result( 6 ) = zero
      DO 120 i = k + 1, min( k+l, m ) - 1
         IF( work( i ).LT.work( i+1 ) )
     $      result( 6 ) = ulpinv
  120 continue
*
      return
*
*     End of DGSVTS
*
      END