NAME

sggsvd - compute the generalized singular value decomposition (GSVD) of an M-by-N real matrix A and P-by-N real matrix B

SYNOPSIS

  SUBROUTINE SGGSVD( JOBU, JOBV, JOBQ, M, N, P, K, L, A, LDA, B, LDB, 
 *      ALPHA, BETA, U, LDU, V, LDV, Q, LDQ, WORK, IWORK3, INFO)
  CHARACTER * 1 JOBU, JOBV, JOBQ
  INTEGER M, N, P, K, L, LDA, LDB, LDU, LDV, LDQ, INFO
  INTEGER IWORK3(*)
  REAL A(LDA,*), B(LDB,*), ALPHA(*), BETA(*), U(LDU,*), V(LDV,*), Q(LDQ,*), WORK(*)
  SUBROUTINE SGGSVD_64( JOBU, JOBV, JOBQ, M, N, P, K, L, A, LDA, B, 
 *      LDB, ALPHA, BETA, U, LDU, V, LDV, Q, LDQ, WORK, IWORK3, INFO)
  CHARACTER * 1 JOBU, JOBV, JOBQ
  INTEGER*8 M, N, P, K, L, LDA, LDB, LDU, LDV, LDQ, INFO
  INTEGER*8 IWORK3(*)
  REAL A(LDA,*), B(LDB,*), ALPHA(*), BETA(*), U(LDU,*), V(LDV,*), Q(LDQ,*), WORK(*)

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], 
 *       IWORK3, [INFO])
  CHARACTER(LEN=1) :: JOBU, JOBV, JOBQ
  INTEGER :: M, N, P, K, L, LDA, LDB, LDU, LDV, LDQ, INFO
  INTEGER, DIMENSION(:) :: IWORK3
  REAL, DIMENSION(:) :: ALPHA, BETA, WORK
  REAL, DIMENSION(:,:) :: A, B, U, V, Q
  SUBROUTINE GGSVD_64( JOBU, JOBV, JOBQ, [M], [N], [P], K, L, A, [LDA], 
 *       B, [LDB], ALPHA, BETA, U, [LDU], V, [LDV], Q, [LDQ], [WORK], 
 *       IWORK3, [INFO])
  CHARACTER(LEN=1) :: JOBU, JOBV, JOBQ
  INTEGER(8) :: M, N, P, K, L, LDA, LDB, LDU, LDV, LDQ, INFO
  INTEGER(8), DIMENSION(:) :: IWORK3
  REAL, DIMENSION(:) :: ALPHA, BETA, WORK
  REAL, DIMENSION(:,:) :: A, B, U, V, Q

C INTERFACE

#include <sunperf.h>

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

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

PURPOSE

sggsvd computes the generalized singular value decomposition (GSVD) of an M-by-N real matrix A and P-by-N real matrix B:

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

where U, V and Q are orthogonal matrices, and Z' is the 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 orthogonal 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 orthonormal 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, 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