1: SUBROUTINE DGGSVD( JOBU, JOBV, JOBQ, M, N, P, K, L, A, LDA, B,
2: $ LDB, ALPHA, BETA, U, LDU, V, LDV, Q, LDQ, WORK,
3: $ IWORK, INFO )
4: *
5: * -- LAPACK driver routine (version 3.2) --
6: * -- LAPACK is a software package provided by Univ. of Tennessee, --
7: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
8: * November 2006
9: *
10: * .. Scalar Arguments ..
11: CHARACTER JOBQ, JOBU, JOBV
12: INTEGER INFO, K, L, LDA, LDB, LDQ, LDU, LDV, M, N, P
13: * ..
14: * .. Array Arguments ..
15: INTEGER IWORK( * )
16: DOUBLE PRECISION A( LDA, * ), ALPHA( * ), B( LDB, * ),
17: $ BETA( * ), Q( LDQ, * ), U( LDU, * ),
18: $ V( LDV, * ), WORK( * )
19: * ..
20: *
21: * Purpose
22: * =======
23: *
24: * DGGSVD computes the generalized singular value decomposition (GSVD)
25: * of an M-by-N real matrix A and P-by-N real matrix B:
26: *
27: * U'*A*Q = D1*( 0 R ), V'*B*Q = D2*( 0 R )
28: *
29: * where U, V and Q are orthogonal matrices, and Z' is the transpose
30: * of Z. Let K+L = the effective numerical rank of the matrix (A',B')',
31: * then R is a K+L-by-K+L nonsingular upper triangular matrix, D1 and
32: * D2 are M-by-(K+L) and P-by-(K+L) "diagonal" matrices and of the
33: * following structures, respectively:
34: *
35: * If M-K-L >= 0,
36: *
37: * K L
38: * D1 = K ( I 0 )
39: * L ( 0 C )
40: * M-K-L ( 0 0 )
41: *
42: * K L
43: * D2 = L ( 0 S )
44: * P-L ( 0 0 )
45: *
46: * N-K-L K L
47: * ( 0 R ) = K ( 0 R11 R12 )
48: * L ( 0 0 R22 )
49: *
50: * where
51: *
52: * C = diag( ALPHA(K+1), ... , ALPHA(K+L) ),
53: * S = diag( BETA(K+1), ... , BETA(K+L) ),
54: * C**2 + S**2 = I.
55: *
56: * R is stored in A(1:K+L,N-K-L+1:N) on exit.
57: *
58: * If M-K-L < 0,
59: *
60: * K M-K K+L-M
61: * D1 = K ( I 0 0 )
62: * M-K ( 0 C 0 )
63: *
64: * K M-K K+L-M
65: * D2 = M-K ( 0 S 0 )
66: * K+L-M ( 0 0 I )
67: * P-L ( 0 0 0 )
68: *
69: * N-K-L K M-K K+L-M
70: * ( 0 R ) = K ( 0 R11 R12 R13 )
71: * M-K ( 0 0 R22 R23 )
72: * K+L-M ( 0 0 0 R33 )
73: *
74: * where
75: *
76: * C = diag( ALPHA(K+1), ... , ALPHA(M) ),
77: * S = diag( BETA(K+1), ... , BETA(M) ),
78: * C**2 + S**2 = I.
79: *
80: * (R11 R12 R13 ) is stored in A(1:M, N-K-L+1:N), and R33 is stored
81: * ( 0 R22 R23 )
82: * in B(M-K+1:L,N+M-K-L+1:N) on exit.
83: *
84: * The routine computes C, S, R, and optionally the orthogonal
85: * transformation matrices U, V and Q.
86: *
87: * In particular, if B is an N-by-N nonsingular matrix, then the GSVD of
88: * A and B implicitly gives the SVD of A*inv(B):
89: * A*inv(B) = U*(D1*inv(D2))*V'.
90: * If ( A',B')' has orthonormal columns, then the GSVD of A and B is
91: * also equal to the CS decomposition of A and B. Furthermore, the GSVD
92: * can be used to derive the solution of the eigenvalue problem:
93: * A'*A x = lambda* B'*B x.
94: * In some literature, the GSVD of A and B is presented in the form
95: * U'*A*X = ( 0 D1 ), V'*B*X = ( 0 D2 )
96: * where U and V are orthogonal and X is nonsingular, D1 and D2 are
97: * ``diagonal''. The former GSVD form can be converted to the latter
98: * form by taking the nonsingular matrix X as
99: *
100: * X = Q*( I 0 )
101: * ( 0 inv(R) ).
102: *
103: * Arguments
104: * =========
105: *
106: * JOBU (input) CHARACTER*1
107: * = 'U': Orthogonal matrix U is computed;
108: * = 'N': U is not computed.
109: *
110: * JOBV (input) CHARACTER*1
111: * = 'V': Orthogonal matrix V is computed;
112: * = 'N': V is not computed.
113: *
114: * JOBQ (input) CHARACTER*1
115: * = 'Q': Orthogonal matrix Q is computed;
116: * = 'N': Q is not computed.
117: *
118: * M (input) INTEGER
119: * The number of rows of the matrix A. M >= 0.
120: *
121: * N (input) INTEGER
122: * The number of columns of the matrices A and B. N >= 0.
123: *
124: * P (input) INTEGER
125: * The number of rows of the matrix B. P >= 0.
126: *
127: * K (output) INTEGER
128: * L (output) INTEGER
129: * On exit, K and L specify the dimension of the subblocks
130: * described in the Purpose section.
131: * K + L = effective numerical rank of (A',B')'.
132: *
133: * A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
134: * On entry, the M-by-N matrix A.
135: * On exit, A contains the triangular matrix R, or part of R.
136: * See Purpose for details.
137: *
138: * LDA (input) INTEGER
139: * The leading dimension of the array A. LDA >= max(1,M).
140: *
141: * B (input/output) DOUBLE PRECISION array, dimension (LDB,N)
142: * On entry, the P-by-N matrix B.
143: * On exit, B contains the triangular matrix R if M-K-L < 0.
144: * See Purpose for details.
145: *
146: * LDB (input) INTEGER
147: * The leading dimension of the array B. LDB >= max(1,P).
148: *
149: * ALPHA (output) DOUBLE PRECISION array, dimension (N)
150: * BETA (output) DOUBLE PRECISION array, dimension (N)
151: * On exit, ALPHA and BETA contain the generalized singular
152: * value pairs of A and B;
153: * ALPHA(1:K) = 1,
154: * BETA(1:K) = 0,
155: * and if M-K-L >= 0,
156: * ALPHA(K+1:K+L) = C,
157: * BETA(K+1:K+L) = S,
158: * or if M-K-L < 0,
159: * ALPHA(K+1:M)=C, ALPHA(M+1:K+L)=0
160: * BETA(K+1:M) =S, BETA(M+1:K+L) =1
161: * and
162: * ALPHA(K+L+1:N) = 0
163: * BETA(K+L+1:N) = 0
164: *
165: * U (output) DOUBLE PRECISION array, dimension (LDU,M)
166: * If JOBU = 'U', U contains the M-by-M orthogonal matrix U.
167: * If JOBU = 'N', U is not referenced.
168: *
169: * LDU (input) INTEGER
170: * The leading dimension of the array U. LDU >= max(1,M) if
171: * JOBU = 'U'; LDU >= 1 otherwise.
172: *
173: * V (output) DOUBLE PRECISION array, dimension (LDV,P)
174: * If JOBV = 'V', V contains the P-by-P orthogonal matrix V.
175: * If JOBV = 'N', V is not referenced.
176: *
177: * LDV (input) INTEGER
178: * The leading dimension of the array V. LDV >= max(1,P) if
179: * JOBV = 'V'; LDV >= 1 otherwise.
180: *
181: * Q (output) DOUBLE PRECISION array, dimension (LDQ,N)
182: * If JOBQ = 'Q', Q contains the N-by-N orthogonal matrix Q.
183: * If JOBQ = 'N', Q is not referenced.
184: *
185: * LDQ (input) INTEGER
186: * The leading dimension of the array Q. LDQ >= max(1,N) if
187: * JOBQ = 'Q'; LDQ >= 1 otherwise.
188: *
189: * WORK (workspace) DOUBLE PRECISION array,
190: * dimension (max(3*N,M,P)+N)
191: *
192: * IWORK (workspace/output) INTEGER array, dimension (N)
193: * On exit, IWORK stores the sorting information. More
194: * precisely, the following loop will sort ALPHA
195: * for I = K+1, min(M,K+L)
196: * swap ALPHA(I) and ALPHA(IWORK(I))
197: * endfor
198: * such that ALPHA(1) >= ALPHA(2) >= ... >= ALPHA(N).
199: *
200: * INFO (output) INTEGER
201: * = 0: successful exit
202: * < 0: if INFO = -i, the i-th argument had an illegal value.
203: * > 0: if INFO = 1, the Jacobi-type procedure failed to
204: * converge. For further details, see subroutine DTGSJA.
205: *
206: * Internal Parameters
207: * ===================
208: *
209: * TOLA DOUBLE PRECISION
210: * TOLB DOUBLE PRECISION
211: * TOLA and TOLB are the thresholds to determine the effective
212: * rank of (A',B')'. Generally, they are set to
213: * TOLA = MAX(M,N)*norm(A)*MAZHEPS,
214: * TOLB = MAX(P,N)*norm(B)*MAZHEPS.
215: * The size of TOLA and TOLB may affect the size of backward
216: * errors of the decomposition.
217: *
218: * Further Details
219: * ===============
220: *
221: * 2-96 Based on modifications by
222: * Ming Gu and Huan Ren, Computer Science Division, University of
223: * California at Berkeley, USA
224: *
225: * =====================================================================
226: *
227: * .. Local Scalars ..
228: LOGICAL WANTQ, WANTU, WANTV
229: INTEGER I, IBND, ISUB, J, NCYCLE
230: DOUBLE PRECISION ANORM, BNORM, SMAX, TEMP, TOLA, TOLB, ULP, UNFL
231: * ..
232: * .. External Functions ..
233: LOGICAL LSAME
234: DOUBLE PRECISION DLAMCH, DLANGE
235: EXTERNAL LSAME, DLAMCH, DLANGE
236: * ..
237: * .. External Subroutines ..
238: EXTERNAL DCOPY, DGGSVP, DTGSJA, XERBLA
239: * ..
240: * .. Intrinsic Functions ..
241: INTRINSIC MAX, MIN
242: * ..
243: * .. Executable Statements ..
244: *
245: * Test the input parameters
246: *
247: WANTU = LSAME( JOBU, 'U' )
248: WANTV = LSAME( JOBV, 'V' )
249: WANTQ = LSAME( JOBQ, 'Q' )
250: *
251: INFO = 0
252: IF( .NOT.( WANTU .OR. LSAME( JOBU, 'N' ) ) ) THEN
253: INFO = -1
254: ELSE IF( .NOT.( WANTV .OR. LSAME( JOBV, 'N' ) ) ) THEN
255: INFO = -2
256: ELSE IF( .NOT.( WANTQ .OR. LSAME( JOBQ, 'N' ) ) ) THEN
257: INFO = -3
258: ELSE IF( M.LT.0 ) THEN
259: INFO = -4
260: ELSE IF( N.LT.0 ) THEN
261: INFO = -5
262: ELSE IF( P.LT.0 ) THEN
263: INFO = -6
264: ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
265: INFO = -10
266: ELSE IF( LDB.LT.MAX( 1, P ) ) THEN
267: INFO = -12
268: ELSE IF( LDU.LT.1 .OR. ( WANTU .AND. LDU.LT.M ) ) THEN
269: INFO = -16
270: ELSE IF( LDV.LT.1 .OR. ( WANTV .AND. LDV.LT.P ) ) THEN
271: INFO = -18
272: ELSE IF( LDQ.LT.1 .OR. ( WANTQ .AND. LDQ.LT.N ) ) THEN
273: INFO = -20
274: END IF
275: IF( INFO.NE.0 ) THEN
276: CALL XERBLA( 'DGGSVD', -INFO )
277: RETURN
278: END IF
279: *
280: * Compute the Frobenius norm of matrices A and B
281: *
282: ANORM = DLANGE( '1', M, N, A, LDA, WORK )
283: BNORM = DLANGE( '1', P, N, B, LDB, WORK )
284: *
285: * Get machine precision and set up threshold for determining
286: * the effective numerical rank of the matrices A and B.
287: *
288: ULP = DLAMCH( 'Precision' )
289: UNFL = DLAMCH( 'Safe Minimum' )
290: TOLA = MAX( M, N )*MAX( ANORM, UNFL )*ULP
291: TOLB = MAX( P, N )*MAX( BNORM, UNFL )*ULP
292: *
293: * Preprocessing
294: *
295: CALL DGGSVP( JOBU, JOBV, JOBQ, M, P, N, A, LDA, B, LDB, TOLA,
296: $ TOLB, K, L, U, LDU, V, LDV, Q, LDQ, IWORK, WORK,
297: $ WORK( N+1 ), INFO )
298: *
299: * Compute the GSVD of two upper "triangular" matrices
300: *
301: CALL DTGSJA( JOBU, JOBV, JOBQ, M, P, N, K, L, A, LDA, B, LDB,
302: $ TOLA, TOLB, ALPHA, BETA, U, LDU, V, LDV, Q, LDQ,
303: $ WORK, NCYCLE, INFO )
304: *
305: * Sort the singular values and store the pivot indices in IWORK
306: * Copy ALPHA to WORK, then sort ALPHA in WORK
307: *
308: CALL DCOPY( N, ALPHA, 1, WORK, 1 )
309: IBND = MIN( L, M-K )
310: DO 20 I = 1, IBND
311: *
312: * Scan for largest ALPHA(K+I)
313: *
314: ISUB = I
315: SMAX = WORK( K+I )
316: DO 10 J = I + 1, IBND
317: TEMP = WORK( K+J )
318: IF( TEMP.GT.SMAX ) THEN
319: ISUB = J
320: SMAX = TEMP
321: END IF
322: 10 CONTINUE
323: IF( ISUB.NE.I ) THEN
324: WORK( K+ISUB ) = WORK( K+I )
325: WORK( K+I ) = SMAX
326: IWORK( K+I ) = K + ISUB
327: ELSE
328: IWORK( K+I ) = K + I
329: END IF
330: 20 CONTINUE
331: *
332: RETURN
333: *
334: * End of DGGSVD
335: *
336: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dggsvd.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack