 LAPACK  3.8.0 LAPACK: Linear Algebra PACKage

## ◆ 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.```
Date
December 2016
Contributors:
Ming Gu and Huan Ren, Computer Science Division, University of California at Berkeley, USA

Definition at line 336 of file sggsvd.f.

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 *
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY
Definition: scopy.f:84
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:69
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
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
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:116
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
Here is the call graph for this function: