NAG Library Function Document

nag_zggsvd (f08vnc)

nag_zggsvd (Nag_OrderType order, Nag_ComputeUType jobu, Nag_ComputeVType jobv, Nag_ComputeQType jobq, Integer m, Integer n, Integer p, Integer *k, Integer *l, Complex a[], Integer pda, Complex b[], Integer pdb, double alpha[], double beta[], Complex u[], Integer pdu, Complex v[], Integer pdv, Complex q[], Integer pdq, Integer iwork[], NagError *fail)

3 Description

The generalized singular value decomposition is given by

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

where

U

V

and

Q

are unitary 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 function computes

C

S

R

and, optionally, the unitary 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 \times B^{- 1}

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

(\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^{H} A x = λ B^{H} B x .

In some literature, the GSVD of

A

and

B

is presented in the form

U^{H} A X = (\begin{matrix} 0 & D_{1} \end{matrix}), V^{H} 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 Arguments

1: order – Nag_OrderTypeInput

On entry: the order argument specifies the two-dimensional storage scheme being used, i.e., row-major ordering or column-major ordering. C language defined storage is specified by

order = Nag_RowMajor

. See Section 3.2.1.3 in the Essential Introduction for a more detailed explanation of the use of this argument.

Constraint:

order = Nag_RowMajor

or Nag_ColMajor.

2: jobu – Nag_ComputeUTypeInput

On entry: if

jobu = Nag_AllU

, the unitary matrix

U

is computed.

jobu = Nag_NotU

U

is not computed.

Constraint:

jobu = Nag_AllU

Nag_NotU

3: jobv – Nag_ComputeVTypeInput

On entry: if

jobv = Nag_ComputeV

, the unitary matrix

V

is computed.

jobv = Nag_NotV

V

is not computed.

Constraint:

jobv = Nag_ComputeV

Nag_NotV

4: jobq – Nag_ComputeQTypeInput

On entry: if

jobq = Nag_ComputeQ

, the unitary matrix

Q

is computed.

jobq = Nag_NotQ

Q

is not computed.

Constraint:

jobq = Nag_ComputeQ

Nag_NotQ

5: m – IntegerInput

On entry:

m

, the number of rows of the matrix

A

Constraint:

m \geq 0

6: n – IntegerInput

On entry:

n

, the number of columns of the matrices

A

and

B

Constraint:

n \geq 0

7: p – IntegerInput

On entry:

p

, the number of rows of the matrix

B

Constraint:

p \geq 0

8: k – Integer *Output

9: l – Integer *Output

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})

10: a[ $\dim$ ] – ComplexInput/Output

Note: the dimension, dim, of the array a must be at least

$\max (1, pda \times n)$ when $order = Nag_ColMajor$ ;
$\max (1, m \times pda)$ when $order = Nag_RowMajor$ .

The

(i, j)

th element of the matrix

A

is stored in

$a [(j - 1) \times pda + i - 1]$ when $order = Nag_ColMajor$ ;
$a [(i - 1) \times pda + j - 1]$ when $order = Nag_RowMajor$ .

On entry: the

m

n

matrix

A

On exit: contains the triangular matrix

R

, or part of

R

. See Section 3 for details.

11: pda – IntegerInput

On entry: the stride separating row or column elements (depending on the value of order) in the array a.

Constraints:

if $order = Nag_ColMajor$ , $pda \geq \max (1, m)$ ;
if $order = Nag_RowMajor$ , $pda \geq \max (1, n)$ .

12: b[ $\dim$ ] – ComplexInput/Output

Note: the dimension, dim, of the array b must be at least

$\max (1, pdb \times n)$ when $order = Nag_ColMajor$ ;
$\max (1, p \times pdb)$ when $order = Nag_RowMajor$ .

The

(i, j)

th element of the matrix

B

is stored in

$b [(j - 1) \times pdb + i - 1]$ when $order = Nag_ColMajor$ ;
$b [(i - 1) \times pdb + j - 1]$ when $order = Nag_RowMajor$ .

On entry: the

p

n

matrix

B

On exit: contains the triangular matrix

R

m - k - l < 0

. See Section 3 for details.

13: pdb – IntegerInput

On entry: the stride separating row or column elements (depending on the value of order) in the array b.

Constraints:

if $order = Nag_ColMajor$ , $pdb \geq \max (1, p)$ ;
if $order = Nag_RowMajor$ , $pdb \geq \max (1, n)$ .

14: alpha[n] – doubleOutput

On exit: see the description of beta.

15: beta[n] – doubleOutput

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 - 1]

, for

i = k, \dots, n

16: u[ $\dim$ ] – ComplexOutput

Note: the dimension, dim, of the array u must be at least

$\max (1, pdu \times m)$ when $jobu = Nag_AllU$ ;
$1$ otherwise.

The

(i, j)

th element of the matrix

U

is stored in

$u [(j - 1) \times pdu + i - 1]$ when $order = Nag_ColMajor$ ;
$u [(i - 1) \times pdu + j - 1]$ when $order = Nag_RowMajor$ .

On exit: if

jobu = Nag_AllU

, u contains the

m

m

unitary matrix

U

jobu = Nag_NotU

, u is not referenced.

17: pdu – IntegerInput

On entry: the stride separating row or column elements (depending on the value of order) in the array u.

Constraints:

if $jobu = Nag_AllU$ , $pdu \geq \max (1, m)$ ;
otherwise $pdu \geq 1$ .

18: v[ $\dim$ ] – ComplexOutput

Note: the dimension, dim, of the array v must be at least

$\max (1, pdv \times p)$ when $jobv = Nag_ComputeV$ ;
$1$ otherwise.

The

(i, j)

th element of the matrix

V

is stored in

$v [(j - 1) \times pdv + i - 1]$ when $order = Nag_ColMajor$ ;
$v [(i - 1) \times pdv + j - 1]$ when $order = Nag_RowMajor$ .

On exit: if

jobv = Nag_ComputeV

, v contains the

p

p

unitary matrix

V

jobv = Nag_NotV

, v is not referenced.

19: pdv – IntegerInput

On entry: the stride separating row or column elements (depending on the value of order) in the array v.

Constraints:

if $jobv = Nag_ComputeV$ , $pdv \geq \max (1, p)$ ;
otherwise $pdv \geq 1$ .

20: q[ $\dim$ ] – ComplexOutput

Note: the dimension, dim, of the array q must be at least

$\max (1, pdq \times n)$ when $jobq = Nag_ComputeQ$ ;
$1$ otherwise.

The

(i, j)

th element of the matrix

Q

is stored in

$q [(j - 1) \times pdq + i - 1]$ when $order = Nag_ColMajor$ ;
$q [(i - 1) \times pdq + j - 1]$ when $order = Nag_RowMajor$ .

On exit: if

jobq = Nag_ComputeQ

, q contains the

n

n

unitary matrix

Q

jobq = Nag_NotQ

, q is not referenced.

21: pdq – IntegerInput

On entry: the stride separating row or column elements (depending on the value of order) in the array q.

Constraints:

if $jobq = Nag_ComputeQ$ , $pdq \geq \max (1, n)$ ;
otherwise $pdq \geq 1$ .

22: iwork[n] – IntegerOutput

On exit: stores the sorting information. More precisely, the following loop will sort alpha such that

alpha [0] \geq alpha [1] \geq \dots \geq alpha [n - 1]

23: fail – NagError *Input/Output

The NAG error argument (see Section 3.6 in the Essential Introduction).

6 Error Indicators and Warnings

NE_ALLOC_FAIL: Dynamic memory allocation failed.
NE_BAD_PARAM: On entry, argument $〈value〉$ had an illegal value.
NE_CONVERGENCE: The Jacobi-type procedure failed to converge.
NE_ENUM_INT_2: On entry, $jobq = 〈value〉$ , $pdq = 〈value〉$ , $n = 〈value〉$ .
Constraint: if $jobq = Nag_ComputeQ$ , $pdq \geq \max (1, n)$ ;
otherwise $pdq \geq 1$ .; On entry, $jobu = 〈value〉$ , $pdu = 〈value〉$ and $m = 〈value〉$ .
Constraint: if $jobu = Nag_AllU$ , $pdu \geq \max (1, m)$ ;
otherwise $pdu \geq 1$ .; On entry, $jobv = 〈value〉$ , $pdv = 〈value〉$ and $p = 〈value〉$ .
Constraint: if $jobv = Nag_ComputeV$ , $pdv \geq \max (1, p)$ ;
otherwise $pdv \geq 1$ .; On entry, $pdq = 〈value〉$ , $jobq = 〈value〉$ and $n = 〈value〉$ .
Constraint: if $jobq = Nag_ComputeQ$ , $pdq \geq \max (1, n)$ ;
otherwise $pdq \geq 1$ .
NE_INT: On entry, $m = 〈value〉$ .
Constraint: $m \geq 0$ .; On entry, $n = 〈value〉$ .
Constraint: $n \geq 0$ .; On entry, $p = 〈value〉$ .
Constraint: $p \geq 0$ .; On entry, $pda = 〈value〉$ .
Constraint: $pda > 0$ .; On entry, $pdb = 〈value〉$ .
Constraint: $pdb > 0$ .; On entry, $pdq = 〈value〉$ .
Constraint: $pdq > 0$ .; On entry, $pdu = 〈value〉$ .
Constraint: $pdu > 0$ .; On entry, $pdv = 〈value〉$ .
Constraint: $pdv > 0$ .
NE_INT_2: On entry, $pda = 〈value〉$ and $m = 〈value〉$ .
Constraint: $pda \geq \max (1, m)$ .; On entry, $pda = 〈value〉$ and $n = 〈value〉$ .
Constraint: $pda \geq \max (1, n)$ .; On entry, $pdb = 〈value〉$ and $n = 〈value〉$ .
Constraint: $pdb \geq \max (1, n)$ .; On entry, $pdb = 〈value〉$ and $p = 〈value〉$ .
Constraint: $pdb \geq \max (1, p)$ .
NE_INTERNAL_ERROR: An internal error has occurred in this function. Check the function call and any array sizes. If the call is correct then please contact NAG for assistance.

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 diagonal elements of the matrix

R

are real.

The real analogue of this function is nag_dggsvd (f08vac).

9 Example

This example finds the generalized singular value decomposition

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

where

A = (\begin{array}{r} 0.96 - 0.81 i & - 0.03 + 0.96 i & - 0.91 + 2.06 i & - 0.05 + 0.41 i \\ - 0.98 + 1.98 i & - 1.20 + 0.19 i & - 0.66 + 0.42 i & - 0.81 + 0.56 i \\ 0.62 - 0.46 i & 1.01 + 0.02 i & 0.63 - 0.17 i & - 1.11 + 0.60 i \\ 0.37 + 0.38 i & 0.19 - 0.54 i & - 0.98 - 0.36 i & 0.22 - 0.20 i \\ 0.83 + 0.51 i & 0.20 + 0.01 i & - 0.17 - 0.46 i & 1.47 + 1.59 i \\ 1.08 - 0.28 i & 0.20 - 0.12 i & - 0.07 + 1.23 i & 0.26 + 0.26 i \end{array})

and

B = (\begin{matrix} 1 & 0 & - 1 & 0 \\ 0 & 1 & 0 & - 1 \end{matrix}),

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 Function Documentnag_zggsvd (f08vnc)

+− Contents

1 Purpose

2 Specification

3 Description

4 References

5 Arguments

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 Function Document

nag_zggsvd (f08vnc)