NAG Library Function Document

nag_dorcsd (f08rac)

nag_dorcsd (Nag_OrderType order, Nag_ComputeUType jobu1, Nag_ComputeUType jobu2, Nag_ComputeVTType jobv1t, Nag_ComputeVTType jobv2t, Nag_SignsType signs, Integer m, Integer p, Integer q, double x11[], Integer pdx11, double x12[], Integer pdx12, double x21[], Integer pdx21, double x22[], Integer pdx22, double theta[], double u1[], Integer pdu1, double u2[], Integer pdu2, double v1t[], Integer pdv1t, double v2t[], Integer pdv2t, NagError *fail)

3 Description

The

m

m

orthogonal matrix

X

is partitioned as

X = (\begin{matrix} X_{11} & X_{12} \\ X_{21} & X_{22} \end{matrix})

where

X_{11}

is a

p

q

submatrix and the dimensions of the other submatrices

X_{12}

X_{21}

and

X_{22}

are such that

X

remains

m

m

The CS decomposition of

X

X = U Σ_{p} V^{T}

where

U

V

and

Σ_{p}

are

m

m

matrices, such that

U = (\begin{matrix} U_{1} & 0 \\ 0 & U_{2} \end{matrix})

is an orthogonal matrix containing the

p

p

orthogonal matrix

U_{1}

and the

(m - p)

(m - p)

orthogonal matrix

U_{2}

;

V = (\begin{matrix} V_{1} & 0 \\ 0 & V_{2} \end{matrix})

is an orthogonal matrix containing the

q

q

orthogonal matrix

V_{1}

and the

(m - q)

(m - q)

orthogonal matrix

V_{2}

; and

Σ_{p} = (\begin{matrix} I_{11} & 0 & 0 & 0 \\ C & 0 & 0 & - S \\ 0 & 0 & 0 & {- I}_{12} \\ 0 & 0 & I_{22} & 0 \\ 0 & S & C & 0 \\ 0 & I_{21} & 0 & 0 \end{matrix})

contains the

r

r

non-negative diagonal submatrices

C

and

S

satisfying

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

, where

r = \min (p, m - p, q, m - q)

and the top left partition is

p

q

The identity matrix

I_{11}

is of order

\min (p, q) - r

and vanishes if

\min (p, q) = r

The identity matrix

I_{12}

is of order

\min (p, m - q) - r

and vanishes if

\min (p, m - q) = r

The identity matrix

I_{21}

is of order

\min (m - p, q) - r

and vanishes if

\min (m - p, q) = r

The identity matrix

I_{22}

is of order

\min (m - p, m - q) - r

and vanishes if

\min (m - p, m - q) = r

In each of the four cases

r = p, q, m - p, m - q

at least two of the identity matrices vanish.

The indicated zeros represent augmentations by additional rows or columns (but not both) to the square diagonal matrices formed by

I_{i j}

and

C

S

Σ_{p}

does not need to be stored in full; it is sufficient to return only the values

θ_{i}

for

i = 1, 2, \dots, r

where

C_{i i} = \cos (θ_{i})

and

S_{i i} = \sin (θ_{i})

The algorithm used to perform the complete

C S

decomposition is described fully in Sutton (2009) including discussions of the stability and accuracy of the algorithm.

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 (2012) Matrix Computations (4th Edition) Johns Hopkins University Press, Baltimore

Sutton B D (2009) Computing the complete

C S

decomposition Numerical Algorithms (Volume 50) 1017–1398 Springer US 33–65 http://dx.doi.org/10.1007/s11075-008-9215-6

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

Nag_ColMajor

2: jobu1 – Nag_ComputeUTypeInput

On entry:

if $jobu1 = Nag_AllU$ , $U_{1}$ is computed;
if $jobu1 = Nag_NotU$ , $U_{1}$ is not computed.

Constraint:

jobu1 = Nag_AllU

Nag_NotU

3: jobu2 – Nag_ComputeUTypeInput

On entry:

if $jobu2 = Nag_AllU$ , $U_{2}$ is computed;
if $jobu2 = Nag_NotU$ , $U_{2}$ is not computed.

Constraint:

jobu2 = Nag_AllU

Nag_NotU

4: jobv1t – Nag_ComputeVTTypeInput

On entry:

if $jobv1t = Nag_AllVT$ , $V_{1}^{T}$ is computed;
if $jobv1t = Nag_NotVT$ , $V_{1}^{T}$ is not computed.

Constraint:

jobv1t = Nag_AllVT

Nag_NotVT

5: jobv2t – Nag_ComputeVTTypeInput

On entry:

if $jobv2t = Nag_AllVT$ , $V_{2}^{T}$ is computed;
if $jobv2t = Nag_NotVT$ , $V_{2}^{T}$ is not computed.

Constraint:

jobv2t = Nag_AllVT

Nag_NotVT

6: signs – Nag_SignsTypeInput

On entry:

if $signs = Nag_LowerMinus$ , the lower-left block is made nonpositive (the other convention);
if $signs = Nag_UpperMinus$ , the upper-right block is made nonpositive (the default convention).

Constraint:

signs = Nag_LowerMinus

Nag_UpperMinus

7: m – IntegerInput

On entry:

m

, the number of rows and columns in the orthogonal matrix

X

Constraint:

m \geq 0

8: p – IntegerInput

On entry:

p

, the number of rows in

X_{11}

and

X_{12}

Constraint:

0 \leq p \leq m

9: q – IntegerInput

On entry:

q

, the number of columns in

X_{11}

and

X_{21}

Constraint:

0 \leq q \leq m

10: x11[ $\dim$ ] – doubleInput/Output

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

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

The

(i, j)

th element of the matrix is stored in

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

On entry: the upper left partition of the orthogonal matrix

X

whose CSD is desired.

On exit: contains details of the orthogonal matrix used in a simultaneous bidiagonalization process.

11: pdx11 – IntegerInput

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

Constraints:

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

12: x12[ $\dim$ ] – doubleInput/Output

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

$\max (1, pdx12 \times p)$ when $order = Nag_RowMajor$ ;
$\max (1, pdx12 \times (m - q))$ when $order = Nag_ColMajor$ .

The

(i, j)

th element of the matrix is stored in

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

On entry: the upper right partition of the orthogonal matrix

X

whose CSD is desired.

On exit: contains details of the orthogonal matrix used in a simultaneous bidiagonalization process.

13: pdx12 – IntegerInput

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

Constraints:

if $order = Nag_RowMajor$ , $pdx12 \geq \max (1, m - q)$ ;
if $order = Nag_ColMajor$ , $pdx12 \geq \max (1, p)$ .

14: x21[ $\dim$ ] – doubleInput/Output

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

$\max (1, pdx21 \times (m - p))$ when $order = Nag_RowMajor$ ;
$\max (1, pdx21 \times q)$ when $order = Nag_ColMajor$ .

The

(i, j)

th element of the matrix is stored in

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

On entry: the lower left partition of the orthogonal matrix

X

whose CSD is desired.

On exit: contains details of the orthogonal matrix used in a simultaneous bidiagonalization process.

15: pdx21 – IntegerInput

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

Constraints:

if $order = Nag_RowMajor$ , $pdx21 \geq \max (1, q)$ ;
if $order = Nag_ColMajor$ , $pdx21 \geq \max (1, m - p)$ .

16: x22[ $\dim$ ] – doubleInput/Output

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

$\max (1, pdx22 \times (m - p))$ when $order = Nag_RowMajor$ ;
$\max (1, pdx22 \times (m - q))$ when $order = Nag_ColMajor$ .

The

(i, j)

th element of the matrix is stored in

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

On entry: the lower right partition of the orthogonal matrix

X

CSD is desired.

On exit: contains details of the orthogonal matrix used in a simultaneous bidiagonalization process.

17: pdx22 – IntegerInput

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

Constraints:

if $order = Nag_RowMajor$ , $pdx22 \geq \max (1, m - q)$ ;
if $order = Nag_ColMajor$ , $pdx22 \geq \max (1, m - p)$ .

18: theta[ $\dim$ ] – doubleOutput

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

\min (p, m - p, q, m - q)

On exit: the values

θ_{i}

for

i = 1, 2, \dots, r

where

r = \min (p, m - p, q, m - q)

. The diagonal submatrices

C

and

S

Σ_{p}

are constructed from these values as

$C = diag (\cos (theta [0]), \dots, \cos (theta [r - 1]))$ and
$S = diag (\sin (theta [0]), \dots, \sin (theta [r - 1]))$ .

19: u1[ $\dim$ ] – doubleOutput

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

$\max (1, pdu1 \times p)$ when $jobu1 = Nag_AllU$ ;
otherwise u1 may be NULL.

The

(i, j)

th element of the matrix is stored in

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

On exit: if

jobu1 = Nag_AllU

, u1 contains the

p

p

orthogonal matrix

U_{1}

20: pdu1 – IntegerInput

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

Constraint: if

jobu1 = Nag_AllU

pdu1 \geq \max (1, p)

21: u2[ $\dim$ ] – doubleOutput

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

$\max (1, pdu2 \times (m - p))$ when $jobu2 = Nag_AllU$ ;
otherwise u2 may be NULL.

The

(i, j)

th element of the matrix is stored in

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

On exit: if

jobu2 = Nag_AllU

, u2 contains the

m - p

m - p

orthogonal matrix

U_{2}

22: pdu2 – IntegerInput

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

Constraint: if

jobu2 = Nag_AllU

pdu2 \geq \max (1, m - p)

23: v1t[ $\dim$ ] – doubleOutput

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

$\max (1, pdv1t \times q)$ when $jobv1t = Nag_AllVT$ ;
otherwise v1t may be NULL.

The

(i, j)

th element of the matrix is stored in

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

On exit: if

jobv1t = Nag_AllVT

, v1t contains the

q

q

orthogonal matrix

{V_{1}}^{T}

24: pdv1t – IntegerInput

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

Constraint: if

jobv1t = Nag_AllVT

pdv1t \geq \max (1, q)

25: v2t[ $\dim$ ] – doubleOutput

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

$\max (1, pdv2t \times (m - q))$ when $jobv2t = Nag_AllVT$ ;
otherwise v2t may be NULL.

The

(i, j)

th element of the matrix is stored in

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

On exit: if

jobv2t = Nag_AllVT

, v2t contains the

m - q

m - q

orthogonal matrix

{V_{2}}^{T}

26: pdv2t – IntegerInput

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

Constraint: if

jobv2t = Nag_AllVT

pdv2t \geq \max (1, m - q)

27: 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 during an internal reduction to bidiagonal-block form. The process requires convergence to $\min (p, m - p, q, m - q)$ values, the value of $fail . errnum$ gives the number of converged values.
NE_ENUM_INT_2: On entry, $jobu1 = ⟨value⟩$ , $pdu1 = ⟨value⟩$ and $p = ⟨value⟩$ .
Constraint: if $jobu1 = Nag_AllU$ , $pdu1 \geq \max (1, p)$ .
On entry, $jobv1t = ⟨value⟩$ , $pdv1t = ⟨value⟩$ and $q = ⟨value⟩$ .
Constraint: if $jobv1t = Nag_AllVT$ , $pdv1t \geq \max (1, q)$ .
NE_ENUM_INT_3: On entry, $jobu2 = ⟨value⟩$ , $pdu2 = ⟨value⟩$ , $m = ⟨value⟩$ and $p = ⟨value⟩$ .
Constraint: if $jobu2 = Nag_AllU$ , $pdu2 \geq \max (1, m - p)$ .
On entry, $jobv2t = ⟨value⟩$ , $pdv2t = ⟨value⟩$ , $m = ⟨value⟩$ and $q = ⟨value⟩$ .
Constraint: if $jobv2t = Nag_AllVT$ , $pdv2t \geq \max (1, m - q)$ .
On entry, $pdx11 = ⟨value⟩$ , $p = ⟨value⟩$ , $q = ⟨value⟩$ and $order = ⟨value⟩$ .
Constraint: if $order = Nag_RowMajor$ , $pdx11 \geq \max (1, q)$ ;
if $order = Nag_ColMajor$ , $pdx11 \geq \max (1, p)$ .
NE_ENUM_INT_4: On entry, $pdx12 = ⟨value⟩$ , $m = ⟨value⟩$ , $p = ⟨value⟩$ , $q = ⟨value⟩$ and $order = ⟨value⟩$ .
Constraint: if $order = Nag_RowMajor$ , $pdx12 \geq \max (1, m - q)$ ;
if $order = Nag_ColMajor$ , $pdx12 \geq \max (1, p)$ .
On entry, $pdx21 = ⟨value⟩$ , $m = ⟨value⟩$ , $p = ⟨value⟩$ , $q = ⟨value⟩$ and $order = ⟨value⟩$ .
Constraint: if $order = Nag_RowMajor$ , $pdx21 \geq \max (1, q)$ ;
if $order = Nag_ColMajor$ , $pdx21 \geq \max (1, m - p)$ .
On entry, $pdx22 = ⟨value⟩$ , $m = ⟨value⟩$ , $p = ⟨value⟩$ , $q = ⟨value⟩$ and $order = ⟨value⟩$ .
Constraint: if $order = Nag_RowMajor$ , $pdx22 \geq \max (1, m - q)$ ;
if $order = Nag_ColMajor$ , $pdx22 \geq \max (1, m - p)$ .
NE_INT: On entry, $m = ⟨value⟩$ .
Constraint: $m \geq 0$ .
NE_INT_2: On entry, $m = ⟨value⟩$ and $p = ⟨value⟩$ .
Constraint: $0 \leq p \leq m$ .
On entry, $m = ⟨value⟩$ and $q = ⟨value⟩$ .
Constraint: $0 \leq q \leq m$ .
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

C S

decomposition is nearly the exact

C S

decomposition for the nearby matrix

(X + E)

, where

{‖E‖}_{2} = O (ε),

and

ε

is the machine precision.

8 Parallelism and Performance

nag_dorcsd (f08rac) is threaded by NAG for parallel execution in multithreaded implementations of the NAG Library.

nag_dorcsd (f08rac) makes calls to BLAS and/or LAPACK routines, which may be threaded within the vendor library used by this implementation. Consult the documentation for the vendor library for further information.

Please consult the Users' Note for your implementation for any additional implementation-specific information.

9 Further Comments

The total number of floating-point operations required to perform the full

C S

decomposition is approximately

2 m^{3}

The complex analogue of this function is nag_zuncsd (f08rnc).

10 Example

This example finds the full CS decomposition of

X = (\begin{array}{r} -0.13484 & 0.52524 & -0.20924 & 0.81373 \\ 0.67420 & -0.52213 & -0.38886 & 0.34874 \\ 0.26968 & 0.52757 & -0.65782 & -0.46499 \\ 0.67420 & 0.41615 & 0.61014 & 0.00000 \end{array})

partitioned in

2

2

blocks.

The decomposition is performed both on submatrices of the orthogonal matrix

X

and on separated partition matrices. Code is also provided to perform a recombining check if required.

NAG Library Function Documentnag_dorcsd (f08rac)

+− Contents

1 Purpose

2 Specification

3 Description

4 References

5 Arguments

6 Error Indicators and Warnings

7 Accuracy

8 Parallelism and Performance

9 Further Comments

10 Example

10.1 Program Text

10.2 Program Data

10.3 Program Results

NAG Library Function Document

nag_dorcsd (f08rac)