 LAPACK 3.11.0 LAPACK: Linear Algebra PACKage
Searching...
No Matches

## ◆ sggsvd()

 subroutine sggsvd ( character JOBU, character JOBV, character JOBQ, integer M, integer N, integer P, integer K, integer L, real, dimension( lda, * ) A, integer LDA, real, dimension( ldb, * ) B, integer LDB, real, dimension( * ) ALPHA, real, dimension( * ) BETA, real, dimension( ldu, * ) U, integer LDU, real, dimension( ldv, * ) V, integer LDV, real, dimension( ldq, * ) Q, integer LDQ, real, dimension( * ) WORK, integer, dimension( * ) IWORK, integer INFO )

SGGSVD computes the singular value decomposition (SVD) for OTHER matrices

Purpose:
``` This routine is deprecated and has been replaced by routine SGGSVD3.

SGGSVD computes the generalized singular value decomposition (GSVD)
of an M-by-N real matrix A and P-by-N real matrix B:

U**T*A*Q = D1*( 0 R ),    V**T*B*Q = D2*( 0 R )

where U, V and Q are orthogonal matrices.
Let K+L = the effective numerical rank of the matrix (A**T,B**T)**T,
then R is a K+L-by-K+L nonsingular upper triangular matrix, D1 and
D2 are M-by-(K+L) and P-by-(K+L) "diagonal" matrices and of the
following structures, respectively:

If M-K-L >= 0,

K  L
D1 =     K ( I  0 )
L ( 0  C )
M-K-L ( 0  0 )

K  L
D2 =   L ( 0  S )
P-L ( 0  0 )

N-K-L  K    L
( 0 R ) = K (  0   R11  R12 )
L (  0    0   R22 )

where

C = diag( ALPHA(K+1), ... , ALPHA(K+L) ),
S = diag( BETA(K+1),  ... , BETA(K+L) ),
C**2 + S**2 = I.

R is stored in A(1:K+L,N-K-L+1:N) on exit.

If M-K-L < 0,

K M-K K+L-M
D1 =   K ( I  0    0   )
M-K ( 0  C    0   )

K M-K K+L-M
D2 =   M-K ( 0  S    0  )
K+L-M ( 0  0    I  )
P-L ( 0  0    0  )

N-K-L  K   M-K  K+L-M
( 0 R ) =     K ( 0    R11  R12  R13  )
M-K ( 0     0   R22  R23  )
K+L-M ( 0     0    0   R33  )

where

C = diag( ALPHA(K+1), ... , ALPHA(M) ),
S = diag( BETA(K+1),  ... , BETA(M) ),
C**2 + S**2 = I.

(R11 R12 R13 ) is stored in A(1:M, N-K-L+1:N), and R33 is stored
( 0  R22 R23 )
in B(M-K+1:L,N+M-K-L+1:N) on exit.

The routine computes C, S, R, and optionally the orthogonal
transformation matrices U, V and Q.

In particular, if B is an N-by-N nonsingular matrix, then the GSVD of
A and B implicitly gives the SVD of A*inv(B):
A*inv(B) = U*(D1*inv(D2))*V**T.
If ( A**T,B**T)**T  has orthonormal columns, then the GSVD of A and B is
also equal to the CS decomposition of A and B. Furthermore, the GSVD
can be used to derive the solution of the eigenvalue problem:
A**T*A x = lambda* B**T*B x.
In some literature, the GSVD of A and B is presented in the form
U**T*A*X = ( 0 D1 ),   V**T*B*X = ( 0 D2 )
where U and V are orthogonal and X is nonsingular, D1 and D2 are
``diagonal''.  The former GSVD form can be converted to the latter
form by taking the nonsingular matrix X as

X = Q*( I   0    )
( 0 inv(R) ).```
Parameters
 [in] JOBU ``` JOBU is CHARACTER*1 = 'U': Orthogonal matrix U is computed; = 'N': U is not computed.``` [in] JOBV ``` JOBV is CHARACTER*1 = 'V': Orthogonal matrix V is computed; = 'N': V is not computed.``` [in] JOBQ ``` JOBQ is CHARACTER*1 = 'Q': Orthogonal matrix Q is computed; = 'N': Q is not computed.``` [in] M ``` M is INTEGER The number of rows of the matrix A. M >= 0.``` [in] N ``` N is INTEGER The number of columns of the matrices A and B. N >= 0.``` [in] P ``` P is INTEGER The number of rows of the matrix B. P >= 0.``` [out] K ` K is INTEGER` [out] L ``` L is INTEGER On exit, K and L specify the dimension of the subblocks described in Purpose. K + L = effective numerical rank of (A**T,B**T)**T.``` [in,out] A ``` A is REAL array, dimension (LDA,N) On entry, the M-by-N matrix A. On exit, A contains the triangular matrix R, or part of R. See Purpose for details.``` [in] LDA ``` LDA is INTEGER The leading dimension of the array A. LDA >= max(1,M).``` [in,out] B ``` B is REAL array, dimension (LDB,N) On entry, the P-by-N matrix B. On exit, B contains the triangular matrix R if M-K-L < 0. See Purpose for details.``` [in] LDB ``` LDB is INTEGER The leading dimension of the array B. LDB >= max(1,P).``` [out] ALPHA ` ALPHA is REAL array, dimension (N)` [out] BETA ``` BETA is REAL array, dimension (N) On exit, ALPHA and BETA contain the generalized singular value pairs of A and B; ALPHA(1:K) = 1, BETA(1:K) = 0, and if M-K-L >= 0, ALPHA(K+1:K+L) = C, BETA(K+1:K+L) = S, or if M-K-L < 0, ALPHA(K+1:M)=C, ALPHA(M+1:K+L)=0 BETA(K+1:M) =S, BETA(M+1:K+L) =1 and ALPHA(K+L+1:N) = 0 BETA(K+L+1:N) = 0``` [out] U ``` U is REAL array, dimension (LDU,M) If JOBU = 'U', U contains the M-by-M orthogonal matrix U. If JOBU = 'N', U is not referenced.``` [in] LDU ``` LDU is INTEGER The leading dimension of the array U. LDU >= max(1,M) if JOBU = 'U'; LDU >= 1 otherwise.``` [out] V ``` V is REAL array, dimension (LDV,P) If JOBV = 'V', V contains the P-by-P orthogonal matrix V. If JOBV = 'N', V is not referenced.``` [in] LDV ``` LDV is INTEGER The leading dimension of the array V. LDV >= max(1,P) if JOBV = 'V'; LDV >= 1 otherwise.``` [out] Q ``` Q is REAL array, dimension (LDQ,N) If JOBQ = 'Q', Q contains the N-by-N orthogonal matrix Q. If JOBQ = 'N', Q is not referenced.``` [in] LDQ ``` LDQ is INTEGER The leading dimension of the array Q. LDQ >= max(1,N) if JOBQ = 'Q'; LDQ >= 1 otherwise.``` [out] WORK ``` WORK is REAL array, dimension (max(3*N,M,P)+N)``` [out] IWORK ``` IWORK is INTEGER array, dimension (N) On exit, IWORK stores the sorting information. More precisely, the following loop will sort ALPHA for I = K+1, min(M,K+L) swap ALPHA(I) and ALPHA(IWORK(I)) endfor such that ALPHA(1) >= ALPHA(2) >= ... >= ALPHA(N).``` [out] INFO ``` INFO is INTEGER = 0: successful exit < 0: if INFO = -i, the i-th argument had an illegal value. > 0: if INFO = 1, the Jacobi-type procedure failed to converge. For further details, see subroutine STGSJA.```
Internal Parameters:
```  TOLA    REAL
TOLB    REAL
TOLA and TOLB are the thresholds to determine the effective
rank of (A**T,B**T)**T. Generally, they are set to
TOLA = MAX(M,N)*norm(A)*MACHEPS,
TOLB = MAX(P,N)*norm(B)*MACHEPS.
The size of TOLA and TOLB may affect the size of backward
errors of the decomposition.```
Contributors:
Ming Gu and Huan Ren, Computer Science Division, University of California at Berkeley, USA

Definition at line 331 of file sggsvd.f.

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*
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
real function slange(NORM, M, N, A, LDA, WORK)
SLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
Definition: slange.f:114
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
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:68
Here is the call graph for this function: