NAME

SYNOPSIS

  SUBROUTINE ZGGSVD( JOBU, JOBV, JOBQ, M, N, P, K, L, A, LDA, B, LDB, 
 *      ALPHA, BETA, U, LDU, V, LDV, Q, LDQ, WORK, WORK2, IWORK3, INFO)
  CHARACTER * 1 JOBU, JOBV, JOBQ
  DOUBLE COMPLEX A(LDA,*), B(LDB,*), U(LDU,*), V(LDV,*), Q(LDQ,*), WORK(*)
  INTEGER M, N, P, K, L, LDA, LDB, LDU, LDV, LDQ, INFO
  INTEGER IWORK3(*)
  DOUBLE PRECISION ALPHA(*), BETA(*), WORK2(*)
 
  SUBROUTINE ZGGSVD_64( JOBU, JOBV, JOBQ, M, N, P, K, L, A, LDA, B, 
 *      LDB, ALPHA, BETA, U, LDU, V, LDV, Q, LDQ, WORK, WORK2, IWORK3, 
 *      INFO)
  CHARACTER * 1 JOBU, JOBV, JOBQ
  DOUBLE COMPLEX A(LDA,*), B(LDB,*), U(LDU,*), V(LDV,*), Q(LDQ,*), WORK(*)
  INTEGER*8 M, N, P, K, L, LDA, LDB, LDU, LDV, LDQ, INFO
  INTEGER*8 IWORK3(*)
  DOUBLE PRECISION ALPHA(*), BETA(*), WORK2(*)
 

F95 INTERFACE

  SUBROUTINE GGSVD( JOBU, JOBV, JOBQ, [M], [N], [P], K, L, A, [LDA], 
 *       B, [LDB], ALPHA, BETA, U, [LDU], V, [LDV], Q, [LDQ], [WORK], 
 *       [WORK2], IWORK3, [INFO])
  CHARACTER(LEN=1) :: JOBU, JOBV, JOBQ
  COMPLEX(8), DIMENSION(:) :: WORK
  COMPLEX(8), DIMENSION(:,:) :: A, B, U, V, Q
  INTEGER :: M, N, P, K, L, LDA, LDB, LDU, LDV, LDQ, INFO
  INTEGER, DIMENSION(:) :: IWORK3
  REAL(8), DIMENSION(:) :: ALPHA, BETA, WORK2
 
  SUBROUTINE GGSVD_64( JOBU, JOBV, JOBQ, [M], [N], [P], K, L, A, [LDA], 
 *       B, [LDB], ALPHA, BETA, U, [LDU], V, [LDV], Q, [LDQ], [WORK], 
 *       [WORK2], IWORK3, [INFO])
  CHARACTER(LEN=1) :: JOBU, JOBV, JOBQ
  COMPLEX(8), DIMENSION(:) :: WORK
  COMPLEX(8), DIMENSION(:,:) :: A, B, U, V, Q
  INTEGER(8) :: M, N, P, K, L, LDA, LDB, LDU, LDV, LDQ, INFO
  INTEGER(8), DIMENSION(:) :: IWORK3
  REAL(8), DIMENSION(:) :: ALPHA, BETA, WORK2
 

C INTERFACE

void zggsvd(char jobu, char jobv, char jobq, int m, int n, int p, int *k, int *l, doublecomplex *a, int lda, doublecomplex *b, int ldb, double *alpha, double *beta, doublecomplex *u, int ldu, doublecomplex *v, int ldv, doublecomplex *q, int ldq, int *iwork3, int *info);

void zggsvd_64(char jobu, char jobv, char jobq, long m, long n, long p, long *k, long *l, doublecomplex *a, long lda, doublecomplex *b, long ldb, double *alpha, double *beta, doublecomplex *u, long ldu, doublecomplex *v, long ldv, doublecomplex *q, long ldq, long *iwork3, long *info);

PURPOSE

      U'*A*Q = D1*( 0 R ),    V'*B*Q = D2*( 0 R )

where U, V and Q are unitary matrices, and Z' means the conjugate transpose of Z. Let K+L = the effective numerical rank of the matrix (A',B')', 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) ``diagonal'' 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'.

If ( A',B')' 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'*A x = lambda* B'*B x.

In some literature, the GSVD of A and B is presented in the form U'*A*X = ( 0 D1 ), V'*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) )

ARGUMENTS

* JOBU (input)

* JOBV (input)

* JOBQ (input)

* M (input)

The number of rows of the matrix A. M >= 0.

* N (input)

The number of columns of the matrices A and B. N >= 0.

* P (input)

The number of rows of the matrix B. P >= 0.

* K (output)

On exit, K and L specify the dimension of the subblocks described in Purpose. K + L = effective numerical rank of (A',B')'.

* L (output)

On exit, K and L specify the dimension of the subblocks described in Purpose. K + L = effective numerical rank of (A',B')'.

* A (input/output)

On entry, the M-by-N matrix A. On exit, A contains the triangular matrix R, or part of R. See Purpose for details.

* LDA (input)

The leading dimension of the array A. LDA >= max(1,M).

* B (input/output)

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.

* LDB (input)

The leading dimension of the array B. LDB >= max(1,P).

* ALPHA (output)

On exit, ALPHA and BETA contain the generalized singular value pairs of A and B; ALPHA(1:K) = 1, 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

* BETA (output)

See description of ALPHA.

* U (output)

If JOBU = 'U', U contains the M-by-M unitary matrix U. If JOBU = 'N', U is not referenced.

* LDU (input)

The leading dimension of the array U. LDU >= max(1,M) if JOBU = 'U'; LDU >= 1 otherwise.

* V (output)

If JOBV = 'V', V contains the P-by-P unitary matrix V. If JOBV = 'N', V is not referenced.

* LDV (input)

The leading dimension of the array V. LDV >= max(1,P) if JOBV = 'V'; LDV >= 1 otherwise.

* Q (output)

If JOBQ = 'Q', Q contains the N-by-N unitary matrix Q. If JOBQ = 'N', Q is not referenced.

* LDQ (input)

The leading dimension of the array Q. LDQ >= max(1,N) if JOBQ = 'Q'; LDQ >= 1 otherwise.

* WORK (workspace)

dimension(MAX(3*N,M,P)+N)

* WORK2 (workspace)

* IWORK3 (output)

On exit, IWORK3 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(IWORK3(I)) endfor such that ALPHA(1) >= ALPHA(2) >= ... >= ALPHA(N).

* INFO (output)