zggsvd.f File Reference

Go to the source code of this file.

Functions/Subroutines
subroutine	zggsvd (jobu, jobv, jobq, m, n, p, k, l, a, lda, b, ldb, alpha, beta, u, ldu, v, ldv, q, ldq, work, rwork, iwork, info)
	ZGGSVD computes the singular value decomposition (SVD) for OTHER matrices

Function/Subroutine Documentation

zggsvd()

subroutine zggsvd	(	character	jobu,
		character	jobv,
		character	jobq,
		integer	m,
		integer	n,
		integer	p,
		integer	k,
		integer	l,
		complex*16, dimension( lda, * )	a,
		integer	lda,
		complex*16, dimension( ldb, * )	b,
		integer	ldb,
		double precision, dimension( * )	alpha,
		double precision, dimension( * )	beta,
		complex*16, dimension( ldu, * )	u,
		integer	ldu,
		complex*16, dimension( ldv, * )	v,
		integer	ldv,
		complex*16, dimension( ldq, * )	q,
		integer	ldq,
		complex*16, dimension( * )	work,
		double precision, dimension( * )	rwork,
		integer, dimension( * )	iwork,
		integer	info )

ZGGSVD computes the singular value decomposition (SVD) for OTHER matrices

Download ZGGSVD + dependencies [TGZ] [ZIP] [TXT]

Purpose:

!>
!> This routine is deprecated and has been replaced by routine ZGGSVD3.
!>
!> ZGGSVD computes the generalized singular value decomposition (GSVD)
!> of an M-by-N complex matrix A and P-by-N complex matrix B:
!>
!>       U**H*A*Q = D1*( 0 R ),    V**H*B*Q = D2*( 0 R )
!>
!> where U, V and Q are unitary matrices.
!> Let K+L = the effective numerical rank of the
!> matrix (A**H,B**H)**H, then R is a (K+L)-by-(K+L) nonsingular upper
!> triangular matrix, D1 and D2 are M-by-(K+L) and P-by-(K+L) 
!> matrices and of the following structures, respectively:
!>
!> If M-K-L >= 0,
!>
!>                     K  L
!>        D1 =     K ( I  0 )
!>                 L ( 0  C )
!>             M-K-L ( 0  0 )
!>
!>                   K  L
!>        D2 =   L ( 0  S )
!>             P-L ( 0  0 )
!>
!>                 N-K-L  K    L
!>   ( 0 R ) = K (  0   R11  R12 )
!>             L (  0    0   R22 )
!> where
!>
!>   C = diag( ALPHA(K+1), ... , ALPHA(K+L) ),
!>   S = diag( BETA(K+1),  ... , BETA(K+L) ),
!>   C**2 + S**2 = I.
!>
!>   R is stored in A(1:K+L,N-K-L+1:N) on exit.
!>
!> If M-K-L < 0,
!>
!>                   K M-K K+L-M
!>        D1 =   K ( I  0    0   )
!>             M-K ( 0  C    0   )
!>
!>                     K M-K K+L-M
!>        D2 =   M-K ( 0  S    0  )
!>             K+L-M ( 0  0    I  )
!>               P-L ( 0  0    0  )
!>
!>                    N-K-L  K   M-K  K+L-M
!>   ( 0 R ) =     K ( 0    R11  R12  R13  )
!>               M-K ( 0     0   R22  R23  )
!>             K+L-M ( 0     0    0   R33  )
!>
!> where
!>
!>   C = diag( ALPHA(K+1), ... , ALPHA(M) ),
!>   S = diag( BETA(K+1),  ... , BETA(M) ),
!>   C**2 + S**2 = I.
!>
!>   (R11 R12 R13 ) is stored in A(1:M, N-K-L+1:N), and R33 is stored
!>   ( 0  R22 R23 )
!>   in B(M-K+1:L,N+M-K-L+1:N) on exit.
!>
!> The routine computes C, S, R, and optionally the unitary
!> transformation matrices U, V and Q.
!>
!> In particular, if B is an N-by-N nonsingular matrix, then the GSVD of
!> A and B implicitly gives the SVD of A*inv(B):
!>                      A*inv(B) = U*(D1*inv(D2))*V**H.
!> If ( A**H,B**H)**H has orthnormal columns, then the GSVD of A and B is also
!> equal to the CS decomposition of A and B. Furthermore, the GSVD can
!> be used to derive the solution of the eigenvalue problem:
!>                      A**H*A x = lambda* B**H*B x.
!> In some literature, the GSVD of A and B is presented in the form
!>                  U**H*A*X = ( 0 D1 ),   V**H*B*X = ( 0 D2 )
!> where U and V are orthogonal and X is nonsingular, and D1 and D2 are
!> ``diagonal''.  The former GSVD form can be converted to the latter
!> form by taking the nonsingular matrix X as
!>
!>                       X = Q*(  I   0    )
!>                             (  0 inv(R) )
!> 

Parameters

[in]	JOBU	!> JOBU is CHARACTER*1 !> = 'U': Unitary matrix U is computed; !> = 'N': U is not computed. !>
[in]	JOBV	!> JOBV is CHARACTER*1 !> = 'V': Unitary matrix V is computed; !> = 'N': V is not computed. !>
[in]	JOBQ	!> JOBQ is CHARACTER*1 !> = 'Q': Unitary matrix Q is computed; !> = 'N': Q is not computed. !>
[in]	M	!> M is INTEGER !> The number of rows of the matrix A. M >= 0. !>
[in]	N	!> N is INTEGER !> The number of columns of the matrices A and B. N >= 0. !>
[in]	P	!> P is INTEGER !> The number of rows of the matrix B. P >= 0. !>
[out]	K	!> K is INTEGER !>
[out]	L	!> L is INTEGER !> !> On exit, K and L specify the dimension of the subblocks !> described in Purpose. !> K + L = effective numerical rank of (A**H,B**H)**H. !>
[in,out]	A	!> A is COMPLEX*16 array, dimension (LDA,N) !> On entry, the M-by-N matrix A. !> On exit, A contains the triangular matrix R, or part of R. !> See Purpose for details. !>
[in]	LDA	!> LDA is INTEGER !> The leading dimension of the array A. LDA >= max(1,M). !>
[in,out]	B	!> B is COMPLEX*16 array, dimension (LDB,N) !> On entry, the P-by-N matrix B. !> On exit, B contains part of the triangular matrix R if !> M-K-L < 0. See Purpose for details. !>
[in]	LDB	!> LDB is INTEGER !> The leading dimension of the array B. LDB >= max(1,P). !>
[out]	ALPHA	!> ALPHA is DOUBLE PRECISION array, dimension (N) !>
[out]	BETA	!> BETA is DOUBLE PRECISION array, dimension (N) !> !> On exit, ALPHA and BETA contain the generalized singular !> value pairs of A and B; !> ALPHA(1:K) = 1, !> BETA(1:K) = 0, !> and if M-K-L >= 0, !> ALPHA(K+1:K+L) = C, !> BETA(K+1:K+L) = S, !> or if M-K-L < 0, !> ALPHA(K+1:M)=C, ALPHA(M+1:K+L)=0 !> BETA(K+1:M) =S, BETA(M+1:K+L) =1 !> and !> ALPHA(K+L+1:N) = 0 !> BETA(K+L+1:N) = 0 !>
[out]	U	!> U is COMPLEX*16 array, dimension (LDU,M) !> If JOBU = 'U', U contains the M-by-M unitary matrix U. !> If JOBU = 'N', U is not referenced. !>
[in]	LDU	!> LDU is INTEGER !> The leading dimension of the array U. LDU >= max(1,M) if !> JOBU = 'U'; LDU >= 1 otherwise. !>
[out]	V	!> V is COMPLEX*16 array, dimension (LDV,P) !> If JOBV = 'V', V contains the P-by-P unitary matrix V. !> If JOBV = 'N', V is not referenced. !>
[in]	LDV	!> LDV is INTEGER !> The leading dimension of the array V. LDV >= max(1,P) if !> JOBV = 'V'; LDV >= 1 otherwise. !>
[out]	Q	!> Q is COMPLEX*16 array, dimension (LDQ,N) !> If JOBQ = 'Q', Q contains the N-by-N unitary matrix Q. !> If JOBQ = 'N', Q is not referenced. !>
[in]	LDQ	!> LDQ is INTEGER !> The leading dimension of the array Q. LDQ >= max(1,N) if !> JOBQ = 'Q'; LDQ >= 1 otherwise. !>
[out]	WORK	!> WORK is COMPLEX*16 array, dimension (max(3*N,M,P)+N) !>
[out]	RWORK	!> RWORK is DOUBLE PRECISION array, dimension (2*N) !>
[out]	IWORK	!> IWORK is INTEGER array, dimension (N) !> On exit, IWORK stores the sorting information. More !> precisely, the following loop will sort ALPHA !> for I = K+1, min(M,K+L) !> swap ALPHA(I) and ALPHA(IWORK(I)) !> endfor !> such that ALPHA(1) >= ALPHA(2) >= ... >= ALPHA(N). !>
[out]	INFO	!> INFO is INTEGER !> = 0: successful exit. !> < 0: if INFO = -i, the i-th argument had an illegal value. !> > 0: if INFO = 1, the Jacobi-type procedure failed to !> converge. For further details, see subroutine ZTGSJA. !>

Internal Parameters:

!>  TOLA    DOUBLE PRECISION
!>  TOLB    DOUBLE PRECISION
!>          TOLA and TOLB are the thresholds to determine the effective
!>          rank of (A**H,B**H)**H. Generally, they are set to
!>                   TOLA = MAX(M,N)*norm(A)*MAZHEPS,
!>                   TOLB = MAX(P,N)*norm(B)*MAZHEPS.
!>          The size of TOLA and TOLB may affect the size of backward
!>          errors of the decomposition.
!> 

Author

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

Contributors:

Ming Gu and Huan Ren, Computer Science Division, University of California at Berkeley, USA

Definition at line 334 of file zggsvd.f.

*
*  -- LAPACK driver routine --
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*
*     .. Scalar Arguments ..
      CHARACTER          JOBQ, JOBU, JOBV
      INTEGER            INFO, K, L, LDA, LDB, LDQ, LDU, LDV, M, N, P
*     ..
*     .. Array Arguments ..
      INTEGER            IWORK( * )
      DOUBLE PRECISION   ALPHA( * ), BETA( * ), RWORK( * )
      COMPLEX*16         A( LDA, * ), B( LDB, * ), Q( LDQ, * ),
     $                   U( LDU, * ), V( LDV, * ), WORK( * )
*     ..
*
*  =====================================================================
*
*     .. Local Scalars ..
      LOGICAL            WANTQ, WANTU, WANTV
      INTEGER            I, IBND, ISUB, J, NCYCLE
      DOUBLE PRECISION   ANORM, BNORM, SMAX, TEMP, TOLA, TOLB, ULP, UNFL
*     ..
*     .. External Functions ..
      LOGICAL            LSAME
      DOUBLE PRECISION   DLAMCH, ZLANGE
      EXTERNAL           lsame, dlamch, zlange
*     ..
*     .. External Subroutines ..
      EXTERNAL           dcopy, xerbla, zggsvp, ztgsja
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          max, min
*     ..
*     .. Executable Statements ..
*
*     Decode and test the input parameters
*
      wantu = lsame( jobu, 'U' )
      wantv = lsame( jobv, 'V' )
      wantq = lsame( jobq, 'Q' )
*
      info = 0
      IF( .NOT.( wantu .OR. lsame( jobu, 'N' ) ) ) THEN
         info = -1
      ELSE IF( .NOT.( wantv .OR. lsame( jobv, 'N' ) ) ) THEN
         info = -2
      ELSE IF( .NOT.( wantq .OR. lsame( jobq, 'N' ) ) ) THEN
         info = -3
      ELSE IF( m.LT.0 ) THEN
         info = -4
      ELSE IF( n.LT.0 ) THEN
         info = -5
      ELSE IF( p.LT.0 ) THEN
         info = -6
      ELSE IF( lda.LT.max( 1, m ) ) THEN
         info = -10
      ELSE IF( ldb.LT.max( 1, p ) ) THEN
         info = -12
      ELSE IF( ldu.LT.1 .OR. ( wantu .AND. ldu.LT.m ) ) THEN
         info = -16
      ELSE IF( ldv.LT.1 .OR. ( wantv .AND. ldv.LT.p ) ) THEN
         info = -18
      ELSE IF( ldq.LT.1 .OR. ( wantq .AND. ldq.LT.n ) ) THEN
         info = -20
      END IF
      IF( info.NE.0 ) THEN
         CALL xerbla( 'ZGGSVD', -info )
         RETURN
      END IF
*
*     Compute the Frobenius norm of matrices A and B
*
      anorm = zlange( '1', m, n, a, lda, rwork )
      bnorm = zlange( '1', p, n, b, ldb, rwork )
*
*     Get machine precision and set up threshold for determining
*     the effective numerical rank of the matrices A and B.
*
      ulp = dlamch( 'Precision' )
      unfl = dlamch( 'Safe Minimum' )
      tola = max( m, n )*max( anorm, unfl )*ulp
      tolb = max( p, n )*max( bnorm, unfl )*ulp
*
      CALL zggsvp( jobu, jobv, jobq, m, p, n, a, lda, b, ldb, tola,
     $             tolb, k, l, u, ldu, v, ldv, q, ldq, iwork, rwork,
     $             work, work( n+1 ), info )
*
*     Compute the GSVD of two upper "triangular" matrices
*
      CALL ztgsja( jobu, jobv, jobq, m, p, n, k, l, a, lda, b, ldb,
     $             tola, tolb, alpha, beta, u, ldu, v, ldv, q, ldq,
     $             work, ncycle, info )
*
*     Sort the singular values and store the pivot indices in IWORK
*     Copy ALPHA to RWORK, then sort ALPHA in RWORK
*
      CALL dcopy( n, alpha, 1, rwork, 1 )
      ibnd = min( l, m-k )
      DO 20 i = 1, ibnd
*
*        Scan for largest ALPHA(K+I)
*
         isub = i
         smax = rwork( k+i )
         DO 10 j = i + 1, ibnd
            temp = rwork( k+j )
            IF( temp.GT.smax ) THEN
               isub = j
               smax = temp
            END IF
   10    CONTINUE
         IF( isub.NE.i ) THEN
            rwork( k+isub ) = rwork( k+i )
            rwork( k+i ) = smax
            iwork( k+i ) = k + isub
         ELSE
            iwork( k+i ) = k + i
         END IF
   20 CONTINUE
*
      RETURN
*
*     End of ZGGSVD
*