Previous Next Contents Generated Index Doc Set Home


Generalized Singular Value Decomposition of General Matrices

The subroutines described in this section compute the matrices C, S, and R from the generalized singular value decomposition of the M×N general matrix A and the P×N general matrix B. Optionally, these subroutines may also compute the orthogonal transformation matrices U, V, and Q.

If the numerical rank of the matrix [AT BT]T is equal to (K+L), and (M - K - L) 0 then the generalized SVD is defined to be:

and

where U, V, and Q are orthogonal or unitary matrices, R is nonsingular upper triangular, [C]2 + [S]2 = [I], and

On exit, R is stored in A(1:K+L, N-K-L+1:N).

If the numerical rank of the matrix [AT BT]T is equal to (K+L), and (M - K - L) < 0 then the generalized SVD is defined to be:

and

where U, V, and Q are orthogonal or unitary matrices, R is nonsingular upper triangular, [C]2 + [S]2 = [I], and

On exit, R(1:M, K+L) is stored in A(1:M, N-K-L+1:N) and R(M+1:K+L, M+1:K+L) is stored in B(M-K+1:L, N+M-K-L+1:N).

Calling Sequence

CALL DGGSVD 

(JOBU, JOBV, JOBQ, M, N, P, K, L, DA, LDA, DB, LDB, 
DALPHA, DBETA, DU, LDU, DV, LDV, DQ, LDQ, DWORK, 
IWORK3, INFO)

CALL SGGSVD 

(JOBU, JOBV, JOBQ, M, N, P, K, L, SA, LDA, SB, LDB, 
SALPHA, SBETA, SU, LDU, SV, LDV, SQ, LDQ, SWORK, 
IWORK3, INFO)

CALL ZGGSVD 

(JOBU, JOBV, JOBQ, M, N, P, K, L, ZA, LDA, ZB, LDB, 
DALPHA, DBETA, ZU, LDU, ZV, LDV, ZQ, LDQ, ZWORK, 
DWORK2, IWORK3, INFO)

CALL CGGSVD 

(JOBU, JOBV, JOBQ, M, N, P, K, L, CA, LDA, CB, LDB, 
SALPHA, SBETA, CU, LDU, CV, LDV, CQ, LDQ, CWORK, 
SWORK2, IWORK3, INFO)
















void dggsvd 







(char jobu, char jobv, char jobq, int m, int n, int p, 
int *k, int *l, double *da, int lda, double *db, int 
ldb, double *dalpha, double *dbeta, double *du, int 
ldu, double *dv, int ldv, double *dq, int ldq, int 
*info)










void sggsvd 







(char jobu, char jobv, char jobq, int m, int n, int p, 
int *k, int *l, float *sa, int lda, float *sb, int ldb, 
float *salpha, float *sbeta, float *su, int ldu, float 
*sv, int ldv, float *q, int ldq, int *info)










void zggsvd 







(char jobu, char jobv, char jobq, int m, int n, int p, 
int *k, int *l, doublecomplex *za, int lda, 
doublecomplex *zb, int ldb, double *dalpha, double 
*dbeta, doublecomplex *zu, int ldu, doublecomplex *zv, 
int ldv, doublecomplex *zq, int ldq, int *info)










void cggsvd 







(char jobu, char jobv, char jobq, int m, int n, int p, 
int *k, int *l, complex *ca, int lda, complex *cb, int 
ldb, float *salpha, float *sbeta, complex *cu, int ldu, 
complex *cv, int ldv, complex *cq, int ldq, int *info)










 














 Arguments







JOBU





Indicates whether or not to compute the matrix U.  Legal values of JOBU are shown below.  Any values not shown below are illegal.














'U' or 'u'





Matrix U is computed.














'N' or 'n'





Matrix U is not computed.








JOBV





Indicates whether or not to compute the matrix V.  Legal values of JOBV are shown below.  Any values not shown below are illegal.














'V' or 'v'





Matrix V is computed.














'N' or 'n'





Matrix V is not computed.








JOBQ





Indicates whether or not to compute the matrix Q.  Legal values of JOBQ are shown below.  Any values not shown below are illegal.














'Q' or 'q'





Matrix Q is computed.














'N' or 'n'





Matrix Q is not computed.








M





Number of rows of the matrix A.  M  0.








N





Number of columns of the matrices A and B.  N  0.








P





Number of rows of the matrix B.  P  0.








K, L





On exit, K and L specify the dimension of the sub-blocks described above.  K+L = effective numerical rank of [AT BT]T.








xA





On entry, the M×N matrix A.
On exit, the triangular matrix R, or part of R.  If M-K-L  0 then R is stored in A(1:K+L, N-K-L+1:N).  If M-K-L < 0 then R(1:M, K+L) is stored in A(1:M, N-K-L+1:N) and R(M+1:K+L, M+1:K+L) is stored in B(M-K+1:L, N+M-K-L+1:N).








LDA





Leading dimension of the array A as specified in a dimension or type statement.  LDA  max(1, M).








xB





On entry, the P×N matrix B.
On exit, the triangular matrix R if necessary.  If M-K-L < 0 then R(1:M, K+L) is stored in A(1:M, N-K-L+1:N) and R(M+1:K+L, M+1:K+L) is stored in B(M-K+1:L, N+M-K-L+1:N).








LDB





Leading dimension of the array B as specified in a dimension or type statement.  LDB  max(1, P).








xALPHA, xBETA














On exit, ALPHA and BETA contain the generalized singular value pairs of A and B:














if M-K-L  0




















ALPHA(1:K) = ONE	ALPHA(K+1:K+L) = C



BETA(1:K) = ZERO	BETA(K+1:K+L) = S














if M-K-L < 0




















ALPHA(1:K) = ONE	ALPHA(K+1:M) = C	ALPHA(M+1:K+L) = ZERO



BETA(1:K) = ZERO	BETA(K+1:M) = S	BETA(M+1:K+L) = ONE














and




















ALPHA(K+L+1:N) = ZERO	 



BETA(K+L+1:N) = ZERO








xU





On exit, if JOBU = 'U' or 'u', then U contains the M×M matrix U; otherwise U is not used.








LDU





Leading dimension of the array U as specified in a dimension or type statement.  LDU  max(1, M).








xV





On exit, if JOBV = 'V' or 'v', then V contains the P×P matrix V; otherwise V is not used.








LDV





Leading dimension of the array V as specified in a dimension or type statement.  LDV  max(1, P).








xQ





On exit, if JOBQ = 'Q' or 'q', then Q contains the N×N matrix Q; otherwise Q is not used.








LDQ





Leading dimension of the array Q as specified in a dimension or type statement.  LDQ  max(1, N).








xWORK





Scratch array with a dimension of max(M, 3 × N, P) + N.








xWORK2





Scratch array with a dimension of 2 × N for complex subroutines.








IWORK3





Scratch array with a dimension of N.








INFO





On exit:




















INFO = 0





Subroutine completed normally.














INFO < 0





The ith argument, where i = |INFO|, had an illegal value.














INFO = 1





Convergence failed.








 














 Sample Program
















      PROGRAM TEST










      IMPLICIT NONE










C










      INTEGER           LDA, LDB, LDIWRK, LDWORK, N










      PARAMETER        (N = 3)










      PARAMETER        (LDA = N)










      PARAMETER        (LDB = N)










      PARAMETER        (LDIWRK = N)










      PARAMETER        (LDWORK = 5 * N)










C










      DOUBLE PRECISION  A(LDA,N), ALPHA(N), BETA(N), B(LDB,N)










      DOUBLE PRECISION  TEMP, WORK(LDWORK)










      INTEGER           ICOL, INFO, IROW, IWORK(LDIWRK), K, L










C










      EXTERNAL          DGGSVD










      INTRINSIC         ABS










C










C     Initialize the array A to store the matrix A shown below.










C     Initialize the array B to store the matrix B shown below.










C










C         2   4  12       2   4  12










C     A = 0  12  20   B = 0   6  10










C         0   0  24       0   0   8










C










      DATA A / 2.0D0, 0.0D0, 0.0D0, 4.0D0, 1.2D1, 0.0D0,










     $         1.2D1, 2.0D1, 2.4D1 /










      DATA B / 2.0D0, 0.0D0, 0.0D0, 4.0D0, 6.0D0, 0.0D0,










     $         1.2D1, 1.0D1, 8.0D0 /










C










C     Print the initial values of A and B.










C










      PRINT 1000










      PRINT 1010, ((A(IROW,ICOL), ICOL = 1, N), IROW = 1, N)










      PRINT 1020










      PRINT 1010, ((B(IROW,ICOL), ICOL = 1, N), IROW = 1, N)










C










C     Compute the generalized singular values of (A,B)










C










      CALL DGGSVD ('NO COMPUTE U', 'NO COMPUTE V', 'NO COMPUTE Q',










     $             N, N, N, K, L, A, LDA, B, LDB, ALPHA, BETA,










     $             TEMP, 1, TEMP, 1, TEMP, 1, WORK, IWORK, INFO)










      IF (INFO .LT. 0) THEN










        PRINT 1030, ABS(INFO)










        STOP 1










      ELSE IF (INFO .GT. 0) THEN










        PRINT 1040










        STOP 2










      END IF










      PRINT 1050










      PRINT 1060, ALPHA










      PRINT 1070










      PRINT 1080, BETA










      PRINT 1090










      DO 100, ICOL = 1, N










        PRINT 1100, (ALPHA(ICOL) / BETA(ICOL))










  100 CONTINUE










C










 1000 FORMAT (1X, 'A:')










 1010 FORMAT (3(3X, F4.1))










 1020 FORMAT (/1X, 'B:')










 1030 FORMAT (1X, 'Illegal value for argument #', I1,










     $        ' in DGGSVD.')










 1040 FORMAT (1X, 'Singular value decomposition failed',










     $        1X, 'to converge.')










 1050 FORMAT (/1X, 'ALPHA:')










 1060 FORMAT (1X, 10(F8.4))










 1070 FORMAT (/1X, 'BETA:')










 1080 FORMAT (1X, 10(F8.4))










 1090 FORMAT (/1X, 'ALPHA(i) / BETA(i):')










 1100 FORMAT (1X, F8.4)










C










      END










 










 














 Sample Output







 










 A:










    2.0    4.0   12.0










    0.0   12.0   20.0










    0.0    0.0   24.0



















 B:










    2.0    4.0   12.0










    0.0    6.0   10.0










    0.0    0.0    8.0



















 ALPHA:










   0.8944  0.9487  0.7071



















 BETA:










   0.4472  0.3162  0.7071



















 ALPHA(i) / BETA(i):










   2.0000










   3.0000










   1.0000



















 























Previous


Next

Contents

Generated Index

Doc Set

Home