NAG Library Routine Document

F08VAF (DGGSVD)

INTEGER	M, N, P, K, L, LDA, LDB, LDU, LDV, LDQ, IWORK(N), INFO
REAL (KIND=nag_wp)	A(LDA,), B(LDB,), ALPHA(N), BETA(N), U(LDU,), V(LDV,), Q(LDQ,), WORK(max(3N,M,P)+N)
CHARACTER(1)	JOBU, JOBV, JOBQ

The routine may be called by its LAPACK name dggsvd.

3 Description

The generalized singular value decomposition is given by

U^{T} A Q = D_{1} (\begin{matrix} 0 & R \end{matrix}), V^{T} B Q = D_{2} (\begin{matrix} 0 & R \end{matrix}),

where

U

V

and

Q

are orthogonal matrices. Let

(k + l)

be the effective numerical rank of the matrix

(\begin{matrix} A \\ B \end{matrix})

, then

R

is a

(k + l)

(k + l)

nonsingular upper triangular matrix,

D_{1}

and

D_{2}

are

m

(k + l)

and

p

(k + l)

‘diagonal’ matrices structured as follows:

m - k - l \geq 0

D_{1} = \begin{array}{rc} k & l \\ k & ( & I & 0 & ) \\ l & 0 & C \\ m - k - l & 0 & 0 \end{array}

D_{2} = \begin{array}{rc} k & l \\ l & ( & 0 & S & ) \\ p - l & 0 & 0 \end{array}

(\begin{matrix} 0 & R \end{matrix}) = \begin{array}{rc} n - k - l & k & l \\ k & ( & 0 & R_{11} & R_{12} & ) \\ l & 0 & 0 & R_{22} \end{array}

where

C = diag (α_{k + 1}, \dots, α_{k + l}),

S = diag (β_{k + 1}, \dots, β_{k + l}),

and

C^{2} + S^{2} = I .

R

is stored as a submatrix of

A

with elements

R_{i j}

stored as

A_{i, n - k - l + j}

on exit.

m - k - l < 0

D_{1} = \begin{array}{rc} k & m - k & k + l - m \\ k & ( & I & 0 & 0 & ) \\ m - k & 0 & C & 0 \end{array}

D_{2} = \begin{array}{rc} k & m - k & k + l - m \\ m - k & ( & 0 & S & 0 & ) \\ k + l - m & 0 & 0 & I \\ p - l & 0 & 0 & 0 \end{array}

(\begin{matrix} 0 & R \end{matrix}) = \begin{array}{rc} n - k - l & k & m - k & k + l - m \\ k & ( & 0 & R_{11} & R_{12} & R_{13} & ) \\ m - k & 0 & 0 & R_{22} & R_{23} \\ k + l - m & 0 & 0 & 0 & R_{33} \end{array}

where

C = diag (α_{k + 1}, \dots, α_{m}),

S = diag (β_{k + 1}, \dots, β_{m}),

and

C^{2} + S^{2} = I .

(\begin{matrix} R_{11} & R_{12} & R_{13} \\ 0 & R_{22} & R_{23} \end{matrix})

is stored as a submatrix of

A

with

R_{i j}

stored as

A_{i, n - k - l + j}

, and

R_{33}

is stored as a submatrix of

B

with

{(R_{33})}_{i j}

stored as

B_{m - k + i, n + m - k - l + j}

The routine computes

C

S

R

and, optionally, the orthogonal transformation matrices

U

V

and

Q

In particular, if

B

is an

n

n

nonsingular matrix, then the GSVD of

A

and

B

implicitly gives the SVD of

A B^{- 1}

A B^{- 1} = U (D_{1} D_{2}^{- 1}) V^{T} .

(\begin{matrix} A \\ B \end{matrix})

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 = λ B^{T} B x .

In some literature, the GSVD of

A

and

B

is presented in the form

U^{T} A X = (\begin{matrix} 0 & D_{1} \end{matrix}), V^{T} B X = (\begin{matrix} 0 & D_{2} \end{matrix}),

where

U

and

V

are orthogonal and

X

is nonsingular, and

D_{1}

and

D_{2}

are ‘diagonal’. The former GSVD form can be converted to the latter form by taking the nonsingular matrix

X

X = Q (\begin{matrix} I & 0 \\ 0 & R^{- 1} \end{matrix}) .

4 References

Anderson E, Bai Z, Bischof C, Blackford S, Demmel J, Dongarra J J, Du Croz J J, Greenbaum A, Hammarling S, McKenney A and Sorensen D (1999) LAPACK Users' Guide (3rd Edition) SIAM, Philadelphia http://www.netlib.org/lapack/lug

Golub G H and Van Loan C F (1996) Matrix Computations (3rd Edition) Johns Hopkins University Press, Baltimore

5 Parameters

1: JOBU – CHARACTER(1)Input

On entry: if

JOBU ='U'

, the orthogonal matrix

U

is computed.

JOBU ='N'

U

is not computed.

Constraint:

JOBU ='U'

'N'

2: JOBV – CHARACTER(1)Input

On entry: if

JOBV ='V'

, the orthogonal matrix

V

is computed.

JOBV ='N'

V

is not computed.

Constraint:

JOBV ='V'

'N'

3: JOBQ – CHARACTER(1)Input

On entry: if

JOBQ ='Q'

, the orthogonal matrix

Q

is computed.

JOBQ ='N'

Q

is not computed.

Constraint:

JOBQ ='Q'

'N'

4: M – INTEGERInput

On entry:

m

, the number of rows of the matrix

A

Constraint:

M \geq 0

5: N – INTEGERInput

On entry:

n

, the number of columns of the matrices

A

and

B

Constraint:

N \geq 0

6: P – INTEGERInput

On entry:

p

, the number of rows of the matrix

B

Constraint:

P \geq 0

7: K – INTEGEROutput

8: L – INTEGEROutput

On exit: K and L specify the dimension of the subblocks

k

and

l

as described in Section 3;

(k + l)

is the effective numerical rank of

(\begin{matrix} A \\ B \end{matrix})

9: A(LDA, $*$ ) – REAL (KIND=nag_wp) arrayInput/Output

Note: the second dimension of the array A must be at least

\max (1, N)

On entry: the

m

n

matrix

A

On exit: contains the triangular matrix

R

, or part of

R

. See Section 3 for details.

10: LDA – INTEGERInput

On entry: the first dimension of the array A as declared in the (sub)program from which F08VAF (DGGSVD) is called.

Constraint:

LDA \geq \max (1, M)

11: B(LDB, $*$ ) – REAL (KIND=nag_wp) arrayInput/Output

Note: the second dimension of the array B must be at least

\max (1, N)

On entry: the

p

n

matrix

B

On exit: contains the triangular matrix

R

m - k - l < 0

. See Section 3 for details.

12: LDB – INTEGERInput

On entry: the first dimension of the array B as declared in the (sub)program from which F08VAF (DGGSVD) is called.

Constraint:

LDB \geq \max (1, P)

13: ALPHA(N) – REAL (KIND=nag_wp) arrayOutput

On exit: see the description of BETA.

14: BETA(N) – REAL (KIND=nag_wp) arrayOutput

On exit: ALPHA and BETA contain the generalized singular value pairs of

A

and

B

α_{i}

and

β_{i}

;

$ALPHA (1 : K) = 1$ ,
$BETA (1 : K) = 0$ ,

and if

m - k - l \geq 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$ .

The notation

ALPHA (K : N)

above refers to consecutive elements

ALPHA (i)

, for

i = K, \dots, N

15: U(LDU, $*$ ) – REAL (KIND=nag_wp) arrayOutput

Note: the second dimension of the array U must be at least

\max (1, M)

JOBU ='U'

, and at least

1

otherwise.

On exit: if

JOBU ='U'

, U contains the

m

m

orthogonal matrix

U

JOBU ='N'

, U is not referenced.

16: LDU – INTEGERInput

On entry: the first dimension of the array U as declared in the (sub)program from which F08VAF (DGGSVD) is called.

Constraints:

if $JOBU ='U'$ , $LDU \geq \max (1, M)$ ;
otherwise $LDU \geq 1$ .

17: V(LDV, $*$ ) – REAL (KIND=nag_wp) arrayOutput

Note: the second dimension of the array V must be at least

\max (1, P)

JOBV ='V'

, and at least

1

otherwise.

On exit: if

JOBV ='V'

, V contains the

p

p

orthogonal matrix

V

JOBV ='N'

, V is not referenced.

18: LDV – INTEGERInput

On entry: the first dimension of the array V as declared in the (sub)program from which F08VAF (DGGSVD) is called.

Constraints:

if $JOBV ='V'$ , $LDV \geq \max (1, P)$ ;
otherwise $LDV \geq 1$ .

19: Q(LDQ, $*$ ) – REAL (KIND=nag_wp) arrayOutput

Note: the second dimension of the array Q must be at least

\max (1, N)

JOBQ ='Q'

, and at least

1

otherwise.

On exit: if

JOBQ ='Q'

, Q contains the

n

n

orthogonal matrix

Q

JOBQ ='N'

, Q is not referenced.

20: LDQ – INTEGERInput

On entry: the first dimension of the array Q as declared in the (sub)program from which F08VAF (DGGSVD) is called.

Constraints:

if $JOBQ ='Q'$ , $LDQ \geq \max (1, N)$ ;
otherwise $LDQ \geq 1$ .

21: WORK( $\max (3 \times N, M, P) + N$ ) – REAL (KIND=nag_wp) arrayWorkspace

22: IWORK(N) – INTEGER arrayOutput

On exit: 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) \geq ALPHA (2) \geq \dots \geq ALPHA (N)

23: INFO – INTEGEROutput

On exit:

INFO = 0

unless the routine detects an error (see Section 6).

6 Error Indicators and Warnings

Errors or warnings detected by the routine:

$INFO < 0$: If $INFO = - i$ , argument $i$ had an illegal value. An explanatory message is output, and execution of the program is terminated.

$INFO = 1$: If $INFO = 1$ , the Jacobi-type procedure failed to converge.

7 Accuracy

The computed generalized singular value decomposition is nearly the exact generalized singular value decomposition for nearby matrices

(A + E)

and

(B + F)

, where

{‖E‖}_{2} = O (ε) {‖A‖}_{2} ​ and ​ {‖F‖}_{2} = O (ε) {‖B‖}_{2},

and

ε

is the machine precision. See Section 4.12 of Anderson et al. (1999) for further details.

8 Further Comments

The complex analogue of this routine is F08VNF (ZGGSVD).

9 Example

This example finds the generalized singular value decomposition

A = U Σ_{1} (\begin{matrix} 0 & R \end{matrix}) Q^{T}, B = V Σ_{2} (\begin{matrix} 0 & R \end{matrix}) Q^{T},

where

A = (\begin{matrix} 1 & 2 & 3 \\ 3 & 2 & 1 \\ 4 & 5 & 6 \\ 7 & 8 & 8 \end{matrix}) and B = (\begin{array}{r} - 2 & - 3 & 3 \\ 4 & 6 & 5 \end{array}),

together with estimates for the condition number of

R

and the error bound for the computed generalized singular values.

The example program assumes that

m \geq n

, and would need slight modification if this is not the case.

NAG Library Routine DocumentF08VAF (DGGSVD)

+− Contents

1 Purpose

2 Specification

3 Description

4 References

5 Parameters

6 Error Indicators and Warnings

7 Accuracy

8 Further Comments

9 Example

9.1 Program Text

9.2 Program Data

9.3 Program Results

NAG Library Routine Document

F08VAF (DGGSVD)