?ggsvd

Computes the generalized singular value decomposition of a pair of general rectangular matrices.

Syntax

FORTRAN 77:

call sggsvd(jobu, jobv, jobq, m, n, p, k, l, a, lda, b, ldb, alpha, beta, u, ldu, v, ldv, q, ldq, work, iwork, info)

call dggsvd(jobu, jobv, jobq, m, n, p, k, l, a, lda, b, ldb, alpha, beta, u, ldu, v, ldv, q, ldq, work, iwork, info)

call cggsvd(jobu, jobv, jobq, m, n, p, k, l, a, lda, b, ldb, alpha, beta, u, ldu, v, ldv, q, ldq, work, rwork, iwork, info)

call zggsvd(jobu, jobv, jobq, m, n, p, k, l, a, lda, b, ldb, alpha, beta, u, ldu, v, ldv, q, ldq, work, rwork, iwork, info)

FORTRAN 95:

call ggsvd(a, b, alpha, beta [, k] [,l] [,u] [,v] [,q] [,iwork] [,info])

C:

lapack_int LAPACKE_sggsvd( int matrix_layout, char jobu, char jobv, char jobq, lapack_int m, lapack_int n, lapack_int p, lapack_int* k, lapack_int* l, float* a, lapack_int lda, float* b, lapack_int ldb, float* alpha, float* beta, float* u, lapack_int ldu, float* v, lapack_int ldv, float* q, lapack_int ldq, lapack_int* iwork );

lapack_int LAPACKE_dggsvd( int matrix_layout, char jobu, char jobv, char jobq, lapack_int m, lapack_int n, lapack_int p, lapack_int* k, lapack_int* l, double* a, lapack_int lda, double* b, lapack_int ldb, double* alpha, double* beta, double* u, lapack_int ldu, double* v, lapack_int ldv, double* q, lapack_int ldq, lapack_int* iwork );

lapack_int LAPACKE_cggsvd( int matrix_layout, char jobu, char jobv, char jobq, lapack_int m, lapack_int n, lapack_int p, lapack_int* k, lapack_int* l, lapack_complex_float* a, lapack_int lda, lapack_complex_float* b, lapack_int ldb, float* alpha, float* beta, lapack_complex_float* u, lapack_int ldu, lapack_complex_float* v, lapack_int ldv, lapack_complex_float* q, lapack_int ldq, lapack_int* iwork );

lapack_int LAPACKE_zggsvd( int matrix_layout, char jobu, char jobv, char jobq, lapack_int m, lapack_int n, lapack_int p, lapack_int* k, lapack_int* l, lapack_complex_double* a, lapack_int lda, lapack_complex_double* b, lapack_int ldb, double* alpha, double* beta, lapack_complex_double* u, lapack_int ldu, lapack_complex_double* v, lapack_int ldv, lapack_complex_double* q, lapack_int ldq, lapack_int* iwork );

Include Files

  • Fortran: mkl.fi
  • Fortran 95: lapack.f90
  • C: mkl.h

Description

The routine computes the generalized singular value decomposition (GSVD) of an m-by-n real/complex matrix A and p-by-n real/complex matrix B:

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

where U, V and Q are orthogonal/unitary matrices and U', V' mean transpose/conjugate transpose of U and V respectively.

Let k+l = the effective numerical rank of the matrix (A', B')', 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,


Equation


Equation


Equation

where

C = diag(alpha(K+1),..., alpha(K+l))

S = diag(beta(K+1),...,beta(K+l))

C2 + S2 = I

R is stored in a(1:k+l, n-k-l+1:n ) on exit.

If m-k-l < 0,


Equation


Equation


Equation

where

C = diag(alpha(K+1),..., alpha(m)),

S = diag(beta(K+1),...,beta(m)),

C2 + S2 = I

On exit, Equation is stored in a(1:m, n-k-l+1:n ) and R33 is stored in b(m-k+1:l, n+m-k-l+1:n ).

The routine computes C, S, R, and optionally the orthogonal/unitary 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*B-1:

A*B-1 = U*(D1*D2-1)*V'.

If (A', B')' 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'**A*x = λ*B'*B*x.

Input Parameters

The data types are given for the Fortran interface. A <datatype> placeholder, if present, is used for the C interface data types in the C interface section above. See C Interface Conventions for the C interface principal conventions and type definitions.

jobu

CHARACTER*1. Must be 'U' or 'N'.

If jobu = 'U', orthogonal/unitary matrix U is computed.

If jobu = 'N', U is not computed.

jobv

CHARACTER*1. Must be 'V' or 'N'.

If jobv = 'V', orthogonal/unitary matrix V is computed.

If jobv = 'N', V is not computed.

jobq

CHARACTER*1. Must be 'Q' or 'N'.

If jobq = 'Q', orthogonal/unitary matrix Q is computed.

If jobq = 'N', Q is not computed.

m

INTEGER. The number of rows of the matrix A (m 0).

n

INTEGER. The number of columns of the matrices A and B (n 0).

p

INTEGER. The number of rows of the matrix B (p 0).

a, b, work

REAL for sggsvd

DOUBLE PRECISION for dggsvd

COMPLEX for cggsvd

DOUBLE COMPLEX for zggsvd.

Arrays:

a(lda,*) contains the m-by-n matrix A.

The second dimension of a must be at least max(1, n).

b(ldb,*) contains the p-by-n matrix B.

The second dimension of b must be at least max(1, n).

work(*) is a workspace array.

The dimension of work must be at least max(3n, m, p)+n.

lda

INTEGER. The leading dimension of a; at least max(1, m).

ldb

INTEGER. The leading dimension of b; at least max(1, p).

ldu

INTEGER. The leading dimension of the array u .

ldu max(1, m) if jobu = 'U'; ldu 1 otherwise.

ldv

INTEGER. The leading dimension of the array v .

ldv max(1, p) if jobv = 'V'; ldv 1 otherwise.

ldq

INTEGER. The leading dimension of the array q .

ldq max(1, n) if jobq = 'Q'; ldq 1 otherwise.

iwork

INTEGER.

Workspace array, DIMENSION at least max(1, n).

rwork

REAL for cggsvd DOUBLE PRECISION for zggsvd.

Workspace array, DIMENSION at least max(1, 2n). Used in complex flavors only.

Output Parameters

k, l

INTEGER. On exit, k and l specify the dimension of the subblocks. The sum k+l is equal to the effective numerical rank of (A', B')'.

a

On exit, a contains the triangular matrix R or part of R.

b

On exit, b contains part of the triangular matrix R if m-k-l < 0.

alpha, beta

REAL for single-precision flavors

DOUBLE PRECISION for double-precision flavors.

Arrays, DIMENSION at least max(1, n) each.

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.

u, v, q

REAL for sggsvd

DOUBLE PRECISION for dggsvd

COMPLEX for cggsvd

DOUBLE COMPLEX for zggsvd.

Arrays:

u(ldu,*); the second dimension of u must be at least max(1, m).

If jobu = 'U', u contains the m-by-m orthogonal/unitary matrix U.

If jobu = 'N', u is not referenced.

v(ldv,*); the second dimension of v must be at least max(1, p).

If jobv = 'V', v contains the p-by-p orthogonal/unitary matrix V.

If jobv = 'N', v is not referenced.

q(ldq,*); the second dimension of q must be at least max(1, n).

If jobq = 'Q', q contains the n-by-n orthogonal/unitary matrix Q.

If jobq = 'N', q is not referenced.

iwork

On exit, iwork stores the sorting information.

info

INTEGER.

If info = 0, the execution is successful.

If info = -i, the i-th parameter had an illegal value.

If info = 1, the Jacobi-type procedure failed to converge. For further details, see subroutine tgsja.

Fortran 95 Interface Notes

Routines in Fortran 95 interface have fewer arguments in the calling sequence than their FORTRAN 77 counterparts. For general conventions applied to skip redundant or restorable arguments, see Fortran 95 Interface Conventions.

Specific details for the routine ggsvd interface are the following:

a

Holds the matrix A of size (m, n).

b

Holds the matrix B of size (p, n).

alpha

Holds the vector of length n.

beta

Holds the vector of length n.

u

Holds the matrix U of size (m, m).

v

Holds the matrix V of size (p, p).

q

Holds the matrix Q of size (n, n).

iwork

Holds the vector of length n.

jobu

Restored based on the presence of the argument u as follows:

jobu = 'U', if u is present, jobu = 'N', if u is omitted.

jobv

Restored based on the presence of the argument v as follows:

jobz = 'V', if v is present,

jobz = 'N', if v is omitted.

jobq

Restored based on the presence of the argument q as follows:

jobz = 'Q', if q is present,

jobz = 'N', if q is omitted.

For more complete information about compiler optimizations, see our Optimization Notice.