Previous: stgevc Up: ../lapack-s.html Next: stpcon

NAME STGSJA - compute the generalized singular value decomposi- tion (GSVD) of two real upper ``triangular (or tra- pezoidal)'' matrices A and B SYNOPSIS SUBROUTINE STGSJA( 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 ) CHARACTER JOBQ, JOBU, JOBV INTEGER INFO, K, L, LDA, LDB, LDQ, LDU, LDV, M, N, NCYCLE, P REAL TOLA, TOLB REAL ALPHA( * ), BETA( * ), A( LDA, * ), B( LDB, * ), Q( LDQ, * ), U( LDU, * ), V( LDV, * ), WORK( * ) PURPOSE STGSJA computes the generalized singular value decomposition (GSVD) of two real upper ``triangular (or trapezoidal)'' matrices A and B. On entry, it is assumed that matrices A and B have the fol- lowing forms, which may be obtained by the preprocessing subroutine SGGSVP for two general M-by-N matrix A and P-by-N matrix B: If M-K-L >= 0 A = ( 0 A12 A13 ) K , B = ( 0 0 B13 ) L ( 0 0 A23 ) L ( 0 0 0 ) P-L ( 0 0 0 ) M-K-L N-K-L K L N-K-L K L if M-K-L < 0 A = ( 0 A12 A13 ) K , B = ( 0 0 B13 ) L ( 0 0 A23 ) M-K ( 0 0 0 ) P-L N-K-L K L N-K-L K L where K-by-K matrix A12 and L-by-L matrix B13 are nonsingu- lar upper triangular. A23 is L-by-L upper triangular if M- K-L > 0, otherwise A23 is L-by-(M-K) upper trapezoidal. On exit, U'*A*Q = D1*( 0 R ), V'*B*Q = D2*( 0 R ), where U, V and Q are orthogonal matrices, Z' denotes the transpose of Z, R is a nonsingular upper triangular matrix, and D1 and D2 are ``diagonal'' matrices, which are of the following structures: If M-K-L >= 0, U'*A*Q = D1*( 0 R ) = K ( I 0 ) * ( 0 R11 R12 ) K L ( 0 C ) ( 0 0 R22 ) L M-K-L ( 0 0 ) N-K-L K L K L V'*B*Q = D2*( 0 R ) = L ( 0 S ) * ( 0 R11 R12 ) K P-L ( 0 0 ) ( 0 0 R22 ) L K L N-K-L K L where C = diag( ALPHA(K+1), ... , ALPHA(K+L) ), S = diag( BETA(K+1), ... , BETA(K+L) ), C**2 + S**2 = I. The nonsingular triangular matrix R = ( R11 R12 ) is stored ( 0 R22 ) in A(1:K+L,N-K-L+1:N) on exit. If M-K-L < 0, U'*A*Q = D1*( 0 R ) = K ( I 0 0 ) * ( 0 R11 R12 R13 ) K M-K ( 0 C 0 ) ( 0 0 R22 R23 ) M-K K M-K K+L-M ( 0 0 0 R33 ) K+L-M N-K-L K M-K K+L-M V'*B*Q = D2*( 0 R ) = M-K ( 0 S 0 ) * ( 0 R11 R12 R13 ) K K+L-M ( 0 0 I ) ( 0 0 R22 R23 ) M-K P-L ( 0 0 0 ) ( 0 0 0 R33 ) K+L-M K M-K K+L-M N-K-L K M-K K+L-M where C = diag( ALPHA(K+1), ... , ALPHA(M) ), S = diag( BETA(K+1), ... , BETA(M) ), C**2 + S**2 = I. R = ( R11 R12 R13 ) is a nonsingular upper triangular matrix, the ( 0 R22 R23 ) ( 0 0 R33 ) first M rows of R are stored in A(1:M, N-K-L+1:N) and R33 is stored in B(M-K+1:L,N+M-K-L+1:N) on exit. The computations of the orthogonal transformation matrices U, V and Q are optional and may also be applied to the input orthogonal matrices U, V and Q. ARGUMENTS JOBU (input) CHARACTER*1 = 'U': U is overwritten on the input orthogonal matrix U; = 'I': U is initialized to the identity matrix; = 'N': U is not computed. JOBV (input) CHARACTER*1 = 'V': V is overwritten on the input orthogonal matrix V; = 'I': V is initialized to the identity matrix; = 'N': V is not computed. JOBQ (input) CHARACTER*1 = 'Q': Q is overwritten on the input orthogonal matrix Q; = 'I': Q is initialized to the identity matrix; = 'N': Q is not computed. M (input) INTEGER The number of rows of the matrix A. M >= 0. P (input) INTEGER The number of rows of the matrix B. P >= 0. N (input) INTEGER The number of columns of the matrices A and B. N >= 0. K (input) INTEGER L (input) INTEGER K and L specify the sub- blocks in the input matrices A and B: A23 = A(K+1:MIN(K+L,M),N-L+1:N) and B13 = B(1:L,,N- L+1:N) of A and B, whose GSVD is going to be com- puted by STGSJA. See Further details. A (input/output) REAL array, dimension (LDA,N) On entry, the M-by-N matrix A. On exit, A(N- K+1:N,1:MIN(K+L,M) ) contains the triangular matrix R or part of R. See Purpose for details. LDA (input) INTEGER The leading dimension of the array A. LDA >= max(1,M). B (input/output) REAL array, dimension (LDB,N) On entry, the P-by-N matrix B. On exit, if neces- sary, B(M-K+1:L,N+M-K-L+1:N) contains a part of R. See Purpose for details. LDB (input) INTEGER The leading dimension of the array B. LDB >= max(1,P). TOLA (input) REAL TOLB (input) REAL TOLA and TOLB are the conver- gence criteria for the Jacobi- Kogbetliantz itera- tion procedure. Generally, they are the same as used in the preprocessing step, say TOLA = MAX(M,N)*norm(A)*MACHEPS, TOLB = MAX(P,N)*norm(B)*MACHEPS. ALPHA (output) REAL array, dimension (N) BETA (output) REAL array, dimension (N) 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) = diag(C), BETA(1:K) = ZERO, BETA(K+1:K+L) = diag(S), and 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. Furthermore, if K+L < N, ALPHA(K+L+1:N) = ZERO BETA(K+L+1:N) = ZERO. U (input/output) REAL array, dimension (LDU,M) On entry, if JOBU = 'U', U contains the orthogonal matrix U, On exit, if JOBU = 'U', U is overwritten on the input orthogonal matrix U. If JOBU = 'I', U is first set to the identity matrix. If JOBU = 'N', U is not referenced. LDU (input) INTEGER The leading dimension of the array U. LDU >= max(1,M). V (input/output) REAL array, dimension (LDV,P) On entry, if JOBV = 'V', V contains the orthogonal matrix V. On exit, if JOBV = 'V', V is overwritten on the input orthogonal matrix V. If JOBV = 'I', U is first set to the identity matrix. If JOBV = 'N', V is not referenced. LDV (input) INTEGER The leading dimension of the array V. LDV >= max(1,P). Q (input/output) REAL array, dimension (LDQ,N) On entry, if JOBQ = 'Q', Q contains the orthogonal matrix Q, On exit, if JOBQ = 'Q', Q is overwritten on the input orthogonal matrix Q. If JOBQ = 'I', Q is first set to the identity matrix. If JOBQ = 'N', Q is not referenced. LDQ (input) INTEGER The leading dimension of the array Q. LDQ >= MAX(1,N). WORK (workspace) REAL array, dimension (2*N) NCYCLE (output) INTEGER The number of cycles required for convergence. INFO (output) INTEGER = 0: successful exit < 0: if INFO = -i, the i-th argument had an illegal value. = 1: the procedure does not converge after MAXIT cycles. PARAMETERS MAXIT INTEGER MAXIT specifies the total loops that the iterative procedure may take. If after MAXIT cycles, the rou- tine fails to converge, we return INFO = 1. Further Details =============== STGSJA essentially uses a variant of Kogbetliantz algorithm to reduce min(L,M-K)-by-L triangular (or trapezoidal) matrix A23 and L-by-L matrix B13 to the form: U1'*A13*Q1 = C1*R1; V1'*B13*Q1 = S1*R1, where U1, V1 and Q1 are orthogonal matrix, and Z' is the transpose of Z. C1 and S1 are diagonal matrices satisfying C1**2 + S1**2 = I, and R1 is an L-by-L nonsingular upper triangular matrix.