LAPACK  3.10.0
LAPACK: Linear Algebra PACKage
dggsvd.f
Go to the documentation of this file.
1 *> \brief <b> DGGSVD computes the singular value decomposition (SVD) for OTHER matrices</b>
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download DGGSVD + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dggsvd.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dggsvd.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dggsvd.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE DGGSVD( JOBU, JOBV, JOBQ, M, N, P, K, L, A, LDA, B,
22 * LDB, ALPHA, BETA, U, LDU, V, LDV, Q, LDQ, WORK,
23 * IWORK, INFO )
24 *
25 * .. Scalar Arguments ..
26 * CHARACTER JOBQ, JOBU, JOBV
27 * INTEGER INFO, K, L, LDA, LDB, LDQ, LDU, LDV, M, N, P
28 * ..
29 * .. Array Arguments ..
30 * INTEGER IWORK( * )
31 * DOUBLE PRECISION A( LDA, * ), ALPHA( * ), B( LDB, * ),
32 * $ BETA( * ), Q( LDQ, * ), U( LDU, * ),
33 * $ V( LDV, * ), WORK( * )
34 * ..
35 *
36 *
37 *> \par Purpose:
38 * =============
39 *>
40 *> \verbatim
41 *>
42 *> This routine is deprecated and has been replaced by routine DGGSVD3.
43 *>
44 *> DGGSVD computes the generalized singular value decomposition (GSVD)
45 *> of an M-by-N real matrix A and P-by-N real matrix B:
46 *>
47 *> U**T*A*Q = D1*( 0 R ), V**T*B*Q = D2*( 0 R )
48 *>
49 *> where U, V and Q are orthogonal matrices.
50 *> Let K+L = the effective numerical rank of the matrix (A**T,B**T)**T,
51 *> then R is a K+L-by-K+L nonsingular upper triangular matrix, D1 and
52 *> D2 are M-by-(K+L) and P-by-(K+L) "diagonal" matrices and of the
53 *> following structures, respectively:
54 *>
55 *> If M-K-L >= 0,
56 *>
57 *> K L
58 *> D1 = K ( I 0 )
59 *> L ( 0 C )
60 *> M-K-L ( 0 0 )
61 *>
62 *> K L
63 *> D2 = L ( 0 S )
64 *> P-L ( 0 0 )
65 *>
66 *> N-K-L K L
67 *> ( 0 R ) = K ( 0 R11 R12 )
68 *> L ( 0 0 R22 )
69 *>
70 *> where
71 *>
72 *> C = diag( ALPHA(K+1), ... , ALPHA(K+L) ),
73 *> S = diag( BETA(K+1), ... , BETA(K+L) ),
74 *> C**2 + S**2 = I.
75 *>
76 *> R is stored in A(1:K+L,N-K-L+1:N) on exit.
77 *>
78 *> If M-K-L < 0,
79 *>
80 *> K M-K K+L-M
81 *> D1 = K ( I 0 0 )
82 *> M-K ( 0 C 0 )
83 *>
84 *> K M-K K+L-M
85 *> D2 = M-K ( 0 S 0 )
86 *> K+L-M ( 0 0 I )
87 *> P-L ( 0 0 0 )
88 *>
89 *> N-K-L K M-K K+L-M
90 *> ( 0 R ) = K ( 0 R11 R12 R13 )
91 *> M-K ( 0 0 R22 R23 )
92 *> K+L-M ( 0 0 0 R33 )
93 *>
94 *> where
95 *>
96 *> C = diag( ALPHA(K+1), ... , ALPHA(M) ),
97 *> S = diag( BETA(K+1), ... , BETA(M) ),
98 *> C**2 + S**2 = I.
99 *>
100 *> (R11 R12 R13 ) is stored in A(1:M, N-K-L+1:N), and R33 is stored
101 *> ( 0 R22 R23 )
102 *> in B(M-K+1:L,N+M-K-L+1:N) on exit.
103 *>
104 *> The routine computes C, S, R, and optionally the orthogonal
105 *> transformation matrices U, V and Q.
106 *>
107 *> In particular, if B is an N-by-N nonsingular matrix, then the GSVD of
108 *> A and B implicitly gives the SVD of A*inv(B):
109 *> A*inv(B) = U*(D1*inv(D2))*V**T.
110 *> If ( A**T,B**T)**T has orthonormal columns, then the GSVD of A and B is
111 *> also equal to the CS decomposition of A and B. Furthermore, the GSVD
112 *> can be used to derive the solution of the eigenvalue problem:
113 *> A**T*A x = lambda* B**T*B x.
114 *> In some literature, the GSVD of A and B is presented in the form
115 *> U**T*A*X = ( 0 D1 ), V**T*B*X = ( 0 D2 )
116 *> where U and V are orthogonal and X is nonsingular, D1 and D2 are
117 *> ``diagonal''. The former GSVD form can be converted to the latter
118 *> form by taking the nonsingular matrix X as
119 *>
120 *> X = Q*( I 0 )
121 *> ( 0 inv(R) ).
122 *> \endverbatim
123 *
124 * Arguments:
125 * ==========
126 *
127 *> \param[in] JOBU
128 *> \verbatim
129 *> JOBU is CHARACTER*1
130 *> = 'U': Orthogonal matrix U is computed;
131 *> = 'N': U is not computed.
132 *> \endverbatim
133 *>
134 *> \param[in] JOBV
135 *> \verbatim
136 *> JOBV is CHARACTER*1
137 *> = 'V': Orthogonal matrix V is computed;
138 *> = 'N': V is not computed.
139 *> \endverbatim
140 *>
141 *> \param[in] JOBQ
142 *> \verbatim
143 *> JOBQ is CHARACTER*1
144 *> = 'Q': Orthogonal matrix Q is computed;
145 *> = 'N': Q is not computed.
146 *> \endverbatim
147 *>
148 *> \param[in] M
149 *> \verbatim
150 *> M is INTEGER
151 *> The number of rows of the matrix A. M >= 0.
152 *> \endverbatim
153 *>
154 *> \param[in] N
155 *> \verbatim
156 *> N is INTEGER
157 *> The number of columns of the matrices A and B. N >= 0.
158 *> \endverbatim
159 *>
160 *> \param[in] P
161 *> \verbatim
162 *> P is INTEGER
163 *> The number of rows of the matrix B. P >= 0.
164 *> \endverbatim
165 *>
166 *> \param[out] K
167 *> \verbatim
168 *> K is INTEGER
169 *> \endverbatim
170 *>
171 *> \param[out] L
172 *> \verbatim
173 *> L is INTEGER
174 *>
175 *> On exit, K and L specify the dimension of the subblocks
176 *> described in Purpose.
177 *> K + L = effective numerical rank of (A**T,B**T)**T.
178 *> \endverbatim
179 *>
180 *> \param[in,out] A
181 *> \verbatim
182 *> A is DOUBLE PRECISION array, dimension (LDA,N)
183 *> On entry, the M-by-N matrix A.
184 *> On exit, A contains the triangular matrix R, or part of R.
185 *> See Purpose for details.
186 *> \endverbatim
187 *>
188 *> \param[in] LDA
189 *> \verbatim
190 *> LDA is INTEGER
191 *> The leading dimension of the array A. LDA >= max(1,M).
192 *> \endverbatim
193 *>
194 *> \param[in,out] B
195 *> \verbatim
196 *> B is DOUBLE PRECISION array, dimension (LDB,N)
197 *> On entry, the P-by-N matrix B.
198 *> On exit, B contains the triangular matrix R if M-K-L < 0.
199 *> See Purpose for details.
200 *> \endverbatim
201 *>
202 *> \param[in] LDB
203 *> \verbatim
204 *> LDB is INTEGER
205 *> The leading dimension of the array B. LDB >= max(1,P).
206 *> \endverbatim
207 *>
208 *> \param[out] ALPHA
209 *> \verbatim
210 *> ALPHA is DOUBLE PRECISION array, dimension (N)
211 *> \endverbatim
212 *>
213 *> \param[out] BETA
214 *> \verbatim
215 *> BETA is DOUBLE PRECISION array, dimension (N)
216 *>
217 *> On exit, ALPHA and BETA contain the generalized singular
218 *> value pairs of A and B;
219 *> ALPHA(1:K) = 1,
220 *> BETA(1:K) = 0,
221 *> and if M-K-L >= 0,
222 *> ALPHA(K+1:K+L) = C,
223 *> BETA(K+1:K+L) = S,
224 *> or if M-K-L < 0,
225 *> ALPHA(K+1:M)=C, ALPHA(M+1:K+L)=0
226 *> BETA(K+1:M) =S, BETA(M+1:K+L) =1
227 *> and
228 *> ALPHA(K+L+1:N) = 0
229 *> BETA(K+L+1:N) = 0
230 *> \endverbatim
231 *>
232 *> \param[out] U
233 *> \verbatim
234 *> U is DOUBLE PRECISION array, dimension (LDU,M)
235 *> If JOBU = 'U', U contains the M-by-M orthogonal matrix U.
236 *> If JOBU = 'N', U is not referenced.
237 *> \endverbatim
238 *>
239 *> \param[in] LDU
240 *> \verbatim
241 *> LDU is INTEGER
242 *> The leading dimension of the array U. LDU >= max(1,M) if
243 *> JOBU = 'U'; LDU >= 1 otherwise.
244 *> \endverbatim
245 *>
246 *> \param[out] V
247 *> \verbatim
248 *> V is DOUBLE PRECISION array, dimension (LDV,P)
249 *> If JOBV = 'V', V contains the P-by-P orthogonal matrix V.
250 *> If JOBV = 'N', V is not referenced.
251 *> \endverbatim
252 *>
253 *> \param[in] LDV
254 *> \verbatim
255 *> LDV is INTEGER
256 *> The leading dimension of the array V. LDV >= max(1,P) if
257 *> JOBV = 'V'; LDV >= 1 otherwise.
258 *> \endverbatim
259 *>
260 *> \param[out] Q
261 *> \verbatim
262 *> Q is DOUBLE PRECISION array, dimension (LDQ,N)
263 *> If JOBQ = 'Q', Q contains the N-by-N orthogonal matrix Q.
264 *> If JOBQ = 'N', Q is not referenced.
265 *> \endverbatim
266 *>
267 *> \param[in] LDQ
268 *> \verbatim
269 *> LDQ is INTEGER
270 *> The leading dimension of the array Q. LDQ >= max(1,N) if
271 *> JOBQ = 'Q'; LDQ >= 1 otherwise.
272 *> \endverbatim
273 *>
274 *> \param[out] WORK
275 *> \verbatim
276 *> WORK is DOUBLE PRECISION array,
277 *> dimension (max(3*N,M,P)+N)
278 *> \endverbatim
279 *>
280 *> \param[out] IWORK
281 *> \verbatim
282 *> IWORK is INTEGER array, dimension (N)
283 *> On exit, IWORK stores the sorting information. More
284 *> precisely, the following loop will sort ALPHA
285 *> for I = K+1, min(M,K+L)
286 *> swap ALPHA(I) and ALPHA(IWORK(I))
287 *> endfor
288 *> such that ALPHA(1) >= ALPHA(2) >= ... >= ALPHA(N).
289 *> \endverbatim
290 *>
291 *> \param[out] INFO
292 *> \verbatim
293 *> INFO is INTEGER
294 *> = 0: successful exit
295 *> < 0: if INFO = -i, the i-th argument had an illegal value.
296 *> > 0: if INFO = 1, the Jacobi-type procedure failed to
297 *> converge. For further details, see subroutine DTGSJA.
298 *> \endverbatim
299 *
300 *> \par Internal Parameters:
301 * =========================
302 *>
303 *> \verbatim
304 *> TOLA DOUBLE PRECISION
305 *> TOLB DOUBLE PRECISION
306 *> TOLA and TOLB are the thresholds to determine the effective
307 *> rank of (A',B')**T. Generally, they are set to
308 *> TOLA = MAX(M,N)*norm(A)*MAZHEPS,
309 *> TOLB = MAX(P,N)*norm(B)*MAZHEPS.
310 *> The size of TOLA and TOLB may affect the size of backward
311 *> errors of the decomposition.
312 *> \endverbatim
313 *
314 * Authors:
315 * ========
316 *
317 *> \author Univ. of Tennessee
318 *> \author Univ. of California Berkeley
319 *> \author Univ. of Colorado Denver
320 *> \author NAG Ltd.
321 *
322 *> \ingroup doubleOTHERsing
323 *
324 *> \par Contributors:
325 * ==================
326 *>
327 *> Ming Gu and Huan Ren, Computer Science Division, University of
328 *> California at Berkeley, USA
329 *>
330 * =====================================================================
331  SUBROUTINE dggsvd( JOBU, JOBV, JOBQ, M, N, P, K, L, A, LDA, B,
332  $ LDB, ALPHA, BETA, U, LDU, V, LDV, Q, LDQ, WORK,
333  $ IWORK, INFO )
334 *
335 * -- LAPACK driver routine --
336 * -- LAPACK is a software package provided by Univ. of Tennessee, --
337 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
338 *
339 * .. Scalar Arguments ..
340  CHARACTER JOBQ, JOBU, JOBV
341  INTEGER INFO, K, L, LDA, LDB, LDQ, LDU, LDV, M, N, P
342 * ..
343 * .. Array Arguments ..
344  INTEGER IWORK( * )
345  DOUBLE PRECISION A( LDA, * ), ALPHA( * ), B( LDB, * ),
346  $ beta( * ), q( ldq, * ), u( ldu, * ),
347  $ v( ldv, * ), work( * )
348 * ..
349 *
350 * =====================================================================
351 *
352 * .. Local Scalars ..
353  LOGICAL WANTQ, WANTU, WANTV
354  INTEGER I, IBND, ISUB, J, NCYCLE
355  DOUBLE PRECISION ANORM, BNORM, SMAX, TEMP, TOLA, TOLB, ULP, UNFL
356 * ..
357 * .. External Functions ..
358  LOGICAL LSAME
359  DOUBLE PRECISION DLAMCH, DLANGE
360  EXTERNAL lsame, dlamch, dlange
361 * ..
362 * .. External Subroutines ..
363  EXTERNAL dcopy, dggsvp, dtgsja, xerbla
364 * ..
365 * .. Intrinsic Functions ..
366  INTRINSIC max, min
367 * ..
368 * .. Executable Statements ..
369 *
370 * Test the input parameters
371 *
372  wantu = lsame( jobu, 'U' )
373  wantv = lsame( jobv, 'V' )
374  wantq = lsame( jobq, 'Q' )
375 *
376  info = 0
377  IF( .NOT.( wantu .OR. lsame( jobu, 'N' ) ) ) THEN
378  info = -1
379  ELSE IF( .NOT.( wantv .OR. lsame( jobv, 'N' ) ) ) THEN
380  info = -2
381  ELSE IF( .NOT.( wantq .OR. lsame( jobq, 'N' ) ) ) THEN
382  info = -3
383  ELSE IF( m.LT.0 ) THEN
384  info = -4
385  ELSE IF( n.LT.0 ) THEN
386  info = -5
387  ELSE IF( p.LT.0 ) THEN
388  info = -6
389  ELSE IF( lda.LT.max( 1, m ) ) THEN
390  info = -10
391  ELSE IF( ldb.LT.max( 1, p ) ) THEN
392  info = -12
393  ELSE IF( ldu.LT.1 .OR. ( wantu .AND. ldu.LT.m ) ) THEN
394  info = -16
395  ELSE IF( ldv.LT.1 .OR. ( wantv .AND. ldv.LT.p ) ) THEN
396  info = -18
397  ELSE IF( ldq.LT.1 .OR. ( wantq .AND. ldq.LT.n ) ) THEN
398  info = -20
399  END IF
400  IF( info.NE.0 ) THEN
401  CALL xerbla( 'DGGSVD', -info )
402  RETURN
403  END IF
404 *
405 * Compute the Frobenius norm of matrices A and B
406 *
407  anorm = dlange( '1', m, n, a, lda, work )
408  bnorm = dlange( '1', p, n, b, ldb, work )
409 *
410 * Get machine precision and set up threshold for determining
411 * the effective numerical rank of the matrices A and B.
412 *
413  ulp = dlamch( 'Precision' )
414  unfl = dlamch( 'Safe Minimum' )
415  tola = max( m, n )*max( anorm, unfl )*ulp
416  tolb = max( p, n )*max( bnorm, unfl )*ulp
417 *
418 * Preprocessing
419 *
420  CALL dggsvp( jobu, jobv, jobq, m, p, n, a, lda, b, ldb, tola,
421  $ tolb, k, l, u, ldu, v, ldv, q, ldq, iwork, work,
422  $ work( n+1 ), info )
423 *
424 * Compute the GSVD of two upper "triangular" matrices
425 *
426  CALL dtgsja( jobu, jobv, jobq, m, p, n, k, l, a, lda, b, ldb,
427  $ tola, tolb, alpha, beta, u, ldu, v, ldv, q, ldq,
428  $ work, ncycle, info )
429 *
430 * Sort the singular values and store the pivot indices in IWORK
431 * Copy ALPHA to WORK, then sort ALPHA in WORK
432 *
433  CALL dcopy( n, alpha, 1, work, 1 )
434  ibnd = min( l, m-k )
435  DO 20 i = 1, ibnd
436 *
437 * Scan for largest ALPHA(K+I)
438 *
439  isub = i
440  smax = work( k+i )
441  DO 10 j = i + 1, ibnd
442  temp = work( k+j )
443  IF( temp.GT.smax ) THEN
444  isub = j
445  smax = temp
446  END IF
447  10 CONTINUE
448  IF( isub.NE.i ) THEN
449  work( k+isub ) = work( k+i )
450  work( k+i ) = smax
451  iwork( k+i ) = k + isub
452  ELSE
453  iwork( k+i ) = k + i
454  END IF
455  20 CONTINUE
456 *
457  RETURN
458 *
459 * End of DGGSVD
460 *
461  END
subroutine dggsvd(JOBU, JOBV, JOBQ, M, N, P, K, L, A, LDA, B, LDB, ALPHA, BETA, U, LDU, V, LDV, Q, LDQ, WORK, IWORK, INFO)
DGGSVD computes the singular value decomposition (SVD) for OTHER matrices
Definition: dggsvd.f:334
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
Definition: dcopy.f:82
subroutine dtgsja(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)
DTGSJA
Definition: dtgsja.f:378
subroutine dggsvp(JOBU, JOBV, JOBQ, M, P, N, A, LDA, B, LDB, TOLA, TOLB, K, L, U, LDU, V, LDV, Q, LDQ, IWORK, TAU, WORK, INFO)
DGGSVP
Definition: dggsvp.f:256