LAPACK  3.10.1 LAPACK: Linear Algebra PACKage

## ◆ dggsvd()

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

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

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

DGGSVD 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION array, dimension (N)` [out] BETA ``` BETA is DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DTGSJA.```
Internal Parameters:
```  TOLA    DOUBLE PRECISION
TOLB    DOUBLE PRECISION
TOLA and TOLB are the thresholds to determine the effective
rank of (A',B')**T. Generally, they are set to
TOLA = MAX(M,N)*norm(A)*MAZHEPS,
TOLB = MAX(P,N)*norm(B)*MAZHEPS.
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 dggsvd.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  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 *
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:69
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
Definition: dcopy.f:82
double precision function dlange(NORM, M, N, A, LDA, WORK)
DLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
Definition: dlange.f:114
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
Here is the call graph for this function: