LAPACK  3.8.0
LAPACK: Linear Algebra PACKage
sggsvd.f
Go to the documentation of this file.
1 *> \brief <b> SGGSVD 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 SGGSVD + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/sggsvd.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/sggsvd.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/sggsvd.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE SGGSVD( 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 * REAL 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 SGGSVD3.
43 *>
44 *> SGGSVD 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 REAL 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 REAL 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 REAL array, dimension (N)
211 *> \endverbatim
212 *>
213 *> \param[out] BETA
214 *> \verbatim
215 *> BETA is REAL 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 REAL 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 REAL 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 REAL 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 REAL 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 STGSJA.
298 *> \endverbatim
299 *
300 *> \par Internal Parameters:
301 * =========================
302 *>
303 *> \verbatim
304 *> TOLA REAL
305 *> TOLB REAL
306 *> TOLA and TOLB are the thresholds to determine the effective
307 *> rank of (A**T,B**T)**T. Generally, they are set to
308 *> TOLA = MAX(M,N)*norm(A)*MACHEPS,
309 *> TOLB = MAX(P,N)*norm(B)*MACHEPS.
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 *> \date December 2016
323 *
324 *> \ingroup realOTHERsing
325 *
326 *> \par Contributors:
327 * ==================
328 *>
329 *> Ming Gu and Huan Ren, Computer Science Division, University of
330 *> California at Berkeley, USA
331 *>
332 * =====================================================================
333  SUBROUTINE sggsvd( JOBU, JOBV, JOBQ, M, N, P, K, L, A, LDA, B,
334  $ LDB, ALPHA, BETA, U, LDU, V, LDV, Q, LDQ, WORK,
335  $ IWORK, INFO )
336 *
337 * -- LAPACK driver routine (version 3.7.0) --
338 * -- LAPACK is a software package provided by Univ. of Tennessee, --
339 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
340 * December 2016
341 *
342 * .. Scalar Arguments ..
343  CHARACTER JOBQ, JOBU, JOBV
344  INTEGER INFO, K, L, LDA, LDB, LDQ, LDU, LDV, M, N, P
345 * ..
346 * .. Array Arguments ..
347  INTEGER IWORK( * )
348  REAL A( lda, * ), ALPHA( * ), B( ldb, * ),
349  $ beta( * ), q( ldq, * ), u( ldu, * ),
350  $ v( ldv, * ), work( * )
351 * ..
352 *
353 * =====================================================================
354 *
355 * .. Local Scalars ..
356  LOGICAL WANTQ, WANTU, WANTV
357  INTEGER I, IBND, ISUB, J, NCYCLE
358  REAL ANORM, BNORM, SMAX, TEMP, TOLA, TOLB, ULP, UNFL
359 * ..
360 * .. External Functions ..
361  LOGICAL LSAME
362  REAL SLAMCH, SLANGE
363  EXTERNAL lsame, slamch, slange
364 * ..
365 * .. External Subroutines ..
366  EXTERNAL scopy, sggsvp, stgsja, xerbla
367 * ..
368 * .. Intrinsic Functions ..
369  INTRINSIC max, min
370 * ..
371 * .. Executable Statements ..
372 *
373 * Test the input parameters
374 *
375  wantu = lsame( jobu, 'U' )
376  wantv = lsame( jobv, 'V' )
377  wantq = lsame( jobq, 'Q' )
378 *
379  info = 0
380  IF( .NOT.( wantu .OR. lsame( jobu, 'N' ) ) ) THEN
381  info = -1
382  ELSE IF( .NOT.( wantv .OR. lsame( jobv, 'N' ) ) ) THEN
383  info = -2
384  ELSE IF( .NOT.( wantq .OR. lsame( jobq, 'N' ) ) ) THEN
385  info = -3
386  ELSE IF( m.LT.0 ) THEN
387  info = -4
388  ELSE IF( n.LT.0 ) THEN
389  info = -5
390  ELSE IF( p.LT.0 ) THEN
391  info = -6
392  ELSE IF( lda.LT.max( 1, m ) ) THEN
393  info = -10
394  ELSE IF( ldb.LT.max( 1, p ) ) THEN
395  info = -12
396  ELSE IF( ldu.LT.1 .OR. ( wantu .AND. ldu.LT.m ) ) THEN
397  info = -16
398  ELSE IF( ldv.LT.1 .OR. ( wantv .AND. ldv.LT.p ) ) THEN
399  info = -18
400  ELSE IF( ldq.LT.1 .OR. ( wantq .AND. ldq.LT.n ) ) THEN
401  info = -20
402  END IF
403  IF( info.NE.0 ) THEN
404  CALL xerbla( 'SGGSVD', -info )
405  RETURN
406  END IF
407 *
408 * Compute the Frobenius norm of matrices A and B
409 *
410  anorm = slange( '1', m, n, a, lda, work )
411  bnorm = slange( '1', p, n, b, ldb, work )
412 *
413 * Get machine precision and set up threshold for determining
414 * the effective numerical rank of the matrices A and B.
415 *
416  ulp = slamch( 'Precision' )
417  unfl = slamch( 'Safe Minimum' )
418  tola = max( m, n )*max( anorm, unfl )*ulp
419  tolb = max( p, n )*max( bnorm, unfl )*ulp
420 *
421 * Preprocessing
422 *
423  CALL sggsvp( jobu, jobv, jobq, m, p, n, a, lda, b, ldb, tola,
424  $ tolb, k, l, u, ldu, v, ldv, q, ldq, iwork, work,
425  $ work( n+1 ), info )
426 *
427 * Compute the GSVD of two upper "triangular" matrices
428 *
429  CALL stgsja( jobu, jobv, jobq, m, p, n, k, l, a, lda, b, ldb,
430  $ tola, tolb, alpha, beta, u, ldu, v, ldv, q, ldq,
431  $ work, ncycle, info )
432 *
433 * Sort the singular values and store the pivot indices in IWORK
434 * Copy ALPHA to WORK, then sort ALPHA in WORK
435 *
436  CALL scopy( n, alpha, 1, work, 1 )
437  ibnd = min( l, m-k )
438  DO 20 i = 1, ibnd
439 *
440 * Scan for largest ALPHA(K+I)
441 *
442  isub = i
443  smax = work( k+i )
444  DO 10 j = i + 1, ibnd
445  temp = work( k+j )
446  IF( temp.GT.smax ) THEN
447  isub = j
448  smax = temp
449  END IF
450  10 CONTINUE
451  IF( isub.NE.i ) THEN
452  work( k+isub ) = work( k+i )
453  work( k+i ) = smax
454  iwork( k+i ) = k + isub
455  ELSE
456  iwork( k+i ) = k + i
457  END IF
458  20 CONTINUE
459 *
460  RETURN
461 *
462 * End of SGGSVD
463 *
464  END
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY
Definition: scopy.f:84
subroutine sggsvd(JOBU, JOBV, JOBQ, M, N, P, K, L, A, LDA, B, LDB, ALPHA, BETA, U, LDU, V, LDV, Q, LDQ, WORK, IWORK, INFO)
SGGSVD computes the singular value decomposition (SVD) for OTHER matrices
Definition: sggsvd.f:336
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)
STGSJA
Definition: stgsja.f:380
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine sggsvp(JOBU, JOBV, JOBQ, M, P, N, A, LDA, B, LDB, TOLA, TOLB, K, L, U, LDU, V, LDV, Q, LDQ, IWORK, TAU, WORK, INFO)
SGGSVP
Definition: sggsvp.f:258