1DGGSVP ‐ orthogonal matrices U, V and Q such that N‐K‐L K L
2U'*A*Q = K ( 0 A12 A13 ) if M‐K‐L >= 0 SUBROUTINE DGGSVP( JOBU,
3JOBV, JOBQ, M, P, N, A, LDA, B, LDB, TOLA, TOLB, K, L, U, LDU, V,
4LDV, Q, LDQ, IWORK, TAU, WORK, INFO )
5 CHARACTER JOBQ, JOBU, JOBV
6 INTEGER INFO, K, L, LDA, LDB, LDQ, LDU, LDV, M, N, P
7 DOUBLE PRECISION TOLA, TOLB
8 INTEGER IWORK( * )
9 DOUBLE PRECISION A( LDA, * ), B( LDB, * ), Q( LDQ, * ), TAU(
10* ), U( LDU, * ), V( LDV, * ), WORK( * ) DGGSVP computes orthogo‐
11nal matrices U, V and Q such that
12 L ( 0 0 A23 )
13 M‐K‐L ( 0 0 0 )
14 N‐K‐L K L
16 = K ( 0 A12 A13 ) if M‐K‐L < 0;
17 M‐K ( 0 0 A23 )
18 N‐K‐L K L
20 V'*B*Q = L ( 0 0 B13 )
21 P‐L ( 0 0 0 )
22where the K‐by‐K matrix A12 and L‐by‐L matrix B13 are nonsingular
24upper triangular; A23 is L‐by‐L upper triangular if M‐K‐L >= 0,
25otherwise A23 is (M‐K)‐by‐L upper trapezoidal. K+L = the effec‐
26tive numerical rank of the (M+P)‐by‐N matrix (A',B')'. Z' de‐
27notes the transpose of Z.
28This decomposition is the preprocessing step for computing the
30Generalized Singular Value Decomposition (GSVD), see subroutine
31DGGSVD.
32JOBU (input) CHARACTER*1 = 'U': Orthogonal matrix U is com‐
34puted;
35= 'N': U is not computed. JOBV (input) CHARACTER*1
36= 'V': Orthogonal matrix V is computed;
37= 'N': V is not computed. JOBQ (input) CHARACTER*1
38= 'Q': Orthogonal matrix Q is computed;
39= 'N': Q is not computed. M (input) INTEGER The number of
40rows of the matrix A. M >= 0. P (input) INTEGER The num‐
41ber of rows of the matrix B. P >= 0. N (input) INTEGER
42The number of columns of the matrices A and B. N >= 0. A
43(input/output) DOUBLE PRECISION array, dimension (LDA,N) On en‐
44try, the M‐by‐N matrix A. On exit, A contains the triangular (or
45trapezoidal) matrix described in the Purpose section. LDA
46(input) INTEGER The leading dimension of the array A. LDA >=
47max(1,M). B (input/output) DOUBLE PRECISION array, dimen‐
48sion (LDB,N) On entry, the P‐by‐N matrix B. On exit, B contains
49the triangular matrix described in the Purpose section. LDB
50(input) INTEGER The leading dimension of the array B. LDB >=
51max(1,P). TOLA (input) DOUBLE PRECISION TOLB (input) DOU‐
52BLE PRECISION TOLA and TOLB are the thresholds to determine the
53effective numerical rank of matrix B and a subblock of A. Gener‐
54ally, they are set to TOLA = MAX(M,N)*norm(A)*MAZHEPS, TOLB =
55MAX(P,N)*norm(B)*MAZHEPS. The size of TOLA and TOLB may affect
56the size of backward errors of the decomposition. K (out‐
57put) INTEGER L (output) INTEGER On exit, K and L specify
58the dimension of the subblocks described in Purpose. K + L = ef‐
59fective numerical rank of (A',B')'. U (output) DOUBLE PRE‐
60CISION array, dimension (LDU,M) If JOBU = 'U', U contains the or‐
61thogonal matrix U. If JOBU = 'N', U is not referenced. LDU
62(input) INTEGER The leading dimension of the array U. LDU >=
63max(1,M) if JOBU = 'U'; LDU >= 1 otherwise. V (output)
64DOUBLE PRECISION array, dimension (LDV,M) If JOBV = 'V', V con‐
65tains the orthogonal matrix V. If JOBV = 'N', V is not refer‐
66enced. LDV (input) INTEGER The leading dimension of the ar‐
67ray V. LDV >= max(1,P) if JOBV = 'V'; LDV >= 1 otherwise. Q
68(output) DOUBLE PRECISION array, dimension (LDQ,N) If JOBQ = 'Q',
69Q contains the orthogonal matrix Q. If JOBQ = 'N', Q is not ref‐
70erenced. LDQ (input) INTEGER The leading dimension of the
71array Q. LDQ >= max(1,N) if JOBQ = 'Q'; LDQ >= 1 otherwise.
72IWORK (workspace) INTEGER array, dimension (N) TAU
73(workspace) DOUBLE PRECISION array, dimension (N) WORK
74(workspace) DOUBLE PRECISION array, dimension (max(3*N,M,P)) INFO
75(output) INTEGER = 0: successful exit
76< 0: if INFO = ‐i, the i‐th argument had an illegal value. The
77subroutine uses LAPACK subroutine DGEQPF for the QR factorization
78with column pivoting to detect the effective numerical rank of
79the a matrix. It may be replaced by a better rank determination
80strategy.
81