LAPACK 3.11.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
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*> \ingroup realOTHERsing
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 sggsvd( 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 REAL 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 REAL ANORM, BNORM, SMAX, TEMP, TOLA, TOLB, ULP, UNFL
356* ..
357* .. External Functions ..
358 LOGICAL LSAME
359 REAL SLAMCH, SLANGE
360 EXTERNAL lsame, slamch, slange
361* ..
362* .. External Subroutines ..
363 EXTERNAL scopy, sggsvp, stgsja, 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( 'SGGSVD', -info )
402 RETURN
403 END IF
404*
405* Compute the Frobenius norm of matrices A and B
406*
407 anorm = slange( '1', m, n, a, lda, work )
408 bnorm = slange( '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 = slamch( 'Precision' )
414 unfl = slamch( '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 sggsvp( 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 stgsja( 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 scopy( 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 SGGSVD
460*
461 END
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
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:378
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:256
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY
Definition: scopy.f:82
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:334