Integer type:  int32  int64  nag_int  show int32  show int32  show int64  show int64  show nag_int  show nag_int

Chapter Contents
Chapter Introduction
NAG Toolbox

# NAG Toolbox: nag_lapack_dorcsd (f08ra)

## Purpose

nag_lapack_dorcsd (f08ra) computes the CS decomposition of a real $m$ by $m$ orthogonal matrix $X$, partitioned into a $2$ by $2$ array of submatrices.

## Syntax

[x11, x12, x21, x22, theta, u1, u2, v1t, v2t, info] = f08ra(m, p, q, x11, x12, x21, x22, 'jobu1', jobu1, 'jobu2', jobu2, 'jobv1t', jobv1t, 'jobv2t', jobv2t, 'trans', trans, 'signs', signs)
[x11, x12, x21, x22, theta, u1, u2, v1t, v2t, info] = nag_lapack_dorcsd(m, p, q, x11, x12, x21, x22, 'jobu1', jobu1, 'jobu2', jobu2, 'jobv1t', jobv1t, 'jobv2t', jobv2t, 'trans', trans, 'signs', signs)

## Description

The $m$ by $m$ orthogonal matrix $X$ is partitioned as
 $X= X11 X12 X21 X22$
where ${X}_{11}$ is a $p$ by $q$ submatrix and the dimensions of the other submatrices ${X}_{12}$, ${X}_{21}$ and ${X}_{22}$ are such that $X$ remains $m$ by $m$.
The CS decomposition of $X$ is $X=U{\Sigma }_{p}{V}^{\mathrm{T}}$ where $U$, $V$ and ${\Sigma }_{p}$ are $m$ by $m$ matrices, such that
 $U= U1 0 0 U2$
is an orthogonal matrix containing the $p$ by $p$ orthogonal matrix ${U}_{1}$ and the $\left(m-p\right)$ by $\left(m-p\right)$ orthogonal matrix ${U}_{2}$;
 $V= V1 0 0 V2$
is an orthogonal matrix containing the $q$ by $q$ orthogonal matrix ${V}_{1}$ and the $\left(m-q\right)$ by $\left(m-q\right)$ orthogonal matrix ${V}_{2}$; and
 $Σp= I11 0 0 0 C 0 0 -S 0 0 0 -I12 0 0 I22 0 0 S C 0 0 I21 0 0$
contains the $r$ by $r$ non-negative diagonal submatrices $C$ and $S$ satisfying ${C}^{2}+{S}^{2}=I$, where $r=\mathrm{min}\phantom{\rule{0.125em}{0ex}}\left(p,m-p,q,m-q\right)$ and the top left partition is $p$ by $q$.
The identity matrix ${I}_{11}$ is of order $\mathrm{min}\phantom{\rule{0.125em}{0ex}}\left(p,q\right)-r$ and vanishes if $\mathrm{min}\phantom{\rule{0.125em}{0ex}}\left(p,q\right)=r$.
The identity matrix ${I}_{12}$ is of order $\mathrm{min}\phantom{\rule{0.125em}{0ex}}\left(p,m-q\right)-r$ and vanishes if $\mathrm{min}\phantom{\rule{0.125em}{0ex}}\left(p,m-q\right)=r$.
The identity matrix ${I}_{21}$ is of order $\mathrm{min}\phantom{\rule{0.125em}{0ex}}\left(m-p,q\right)-r$ and vanishes if $\mathrm{min}\phantom{\rule{0.125em}{0ex}}\left(m-p,q\right)=r$.
The identity matrix ${I}_{22}$ is of order $\mathrm{min}\phantom{\rule{0.125em}{0ex}}\left(m-p,m-q\right)-r$ and vanishes if $\mathrm{min}\phantom{\rule{0.125em}{0ex}}\left(m-p,m-q\right)=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}_{ij}$ and $C$ or $S$.
${\Sigma }_{p}$ does not need to be stored in full; it is sufficient to return only the values ${\theta }_{i}$ for $i=1,2,\dots ,r$ where ${C}_{ii}=\mathrm{cos}\left({\theta }_{i}\right)$ and ${S}_{ii}=\mathrm{sin}\left({\theta }_{i}\right)$.
The algorithm used to perform the complete $CS$ decomposition is described fully in Sutton (2009) including discussions of the stability and accuracy of the algorithm.

## 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 $CS$ decomposition Numerical Algorithms (Volume 50) 1017–1398 Springer US 33–65 http://dx.doi.org/10.1007/s11075-008-9215-6

## Parameters

### Compulsory Input Parameters

1:     $\mathrm{m}$int64int32nag_int scalar
$m$, the number of rows and columns in the orthogonal matrix $X$.
Constraint: ${\mathbf{m}}\ge 0$.
2:     $\mathrm{p}$int64int32nag_int scalar
$p$, the number of rows in ${X}_{11}$ and ${X}_{12}$.
Constraint: $0\le {\mathbf{p}}\le {\mathbf{m}}$.
3:     $\mathrm{q}$int64int32nag_int scalar
$q$, the number of columns in ${X}_{11}$ and ${X}_{21}$.
Constraint: $0\le {\mathbf{q}}\le {\mathbf{m}}$.
4:     $\mathrm{x11}\left(\mathit{ldx11},:\right)$ – double array
The first dimension, $\mathit{ldx11}$, of the array x11 must satisfy
• if ${\mathbf{trans}}=\text{'T'}$, $\mathit{ldx11}\ge \mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{q}}\right)$;
• otherwise $\mathit{ldx11}\ge \mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{p}}\right)$.
The second dimension of the array x11 must be at least $\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{p}}\right)$ if ${\mathbf{trans}}=\text{'T'}$, and at least ${\mathbf{q}}$ otherwise.
The upper left partition of the orthogonal matrix $X$ whose CSD is desired.
5:     $\mathrm{x12}\left(\mathit{ldx12},:\right)$ – double array
The first dimension, $\mathit{ldx12}$, of the array x12 must satisfy
• if ${\mathbf{trans}}=\text{'T'}$, $\mathit{ldx12}\ge \mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{m}}-{\mathbf{q}}\right)$;
• otherwise $\mathit{ldx12}\ge \mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{p}}\right)$.
The second dimension of the array x12 must be at least $\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{p}}\right)$ if ${\mathbf{trans}}=\text{'T'}$, and at least ${\mathbf{m}}-{\mathbf{q}}$ otherwise.
The upper right partition of the orthogonal matrix $X$ whose CSD is desired.
6:     $\mathrm{x21}\left(\mathit{ldx21},:\right)$ – double array
The first dimension, $\mathit{ldx21}$, of the array x21 must satisfy
• if ${\mathbf{trans}}=\text{'T'}$, $\mathit{ldx21}\ge \mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{q}}\right)$;
• otherwise $\mathit{ldx21}\ge \mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{m}}-{\mathbf{p}}\right)$.
The second dimension of the array x21 must be at least $\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{m}}-{\mathbf{p}}\right)$ if ${\mathbf{trans}}=\text{'T'}$, and at least ${\mathbf{q}}$ otherwise.
The lower left partition of the orthogonal matrix $X$ whose CSD is desired.
7:     $\mathrm{x22}\left(\mathit{ldx22},:\right)$ – double array
The first dimension, $\mathit{ldx22}$, of the array x22 must satisfy
• if ${\mathbf{trans}}=\text{'T'}$, $\mathit{ldx22}\ge \mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{m}}-{\mathbf{q}}\right)$;
• otherwise $\mathit{ldx22}\ge \mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{m}}-{\mathbf{p}}\right)$.
The second dimension of the array x22 must be at least $\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{m}}-{\mathbf{p}}\right)$ if ${\mathbf{trans}}=\text{'T'}$, and at least ${\mathbf{m}}-{\mathbf{q}}$ otherwise.
The lower right partition of the orthogonal matrix $X$ CSD is desired.

### Optional Input Parameters

1:     $\mathrm{jobu1}$ – string (length ≥ 1)
Default: $\text{'Y'}$
• if ${\mathbf{jobu1}}=\text{'Y'}$, ${U}_{1}$ is computed;
• otherwise, ${U}_{1}$ is not computed.
2:     $\mathrm{jobu2}$ – string (length ≥ 1)
Default: $\text{'Y'}$
• if ${\mathbf{jobu2}}=\text{'Y'}$, ${U}_{2}$ is computed;
• otherwise, ${U}_{2}$ is not computed.
3:     $\mathrm{jobv1t}$ – string (length ≥ 1)
Default: $\text{'Y'}$
• if ${\mathbf{jobv1t}}=\text{'Y'}$, ${V}_{1}^{\mathrm{T}}$ is computed;
• otherwise, ${V}_{1}^{\mathrm{T}}$ is not computed.
4:     $\mathrm{jobv2t}$ – string (length ≥ 1)
Default: $\text{'Y'}$
• if ${\mathbf{jobv2t}}=\text{'Y'}$, ${V}_{2}^{\mathrm{T}}$ is computed;
• otherwise, ${V}_{2}^{\mathrm{T}}$ is not computed.
5:     $\mathrm{trans}$ – string (length ≥ 1)
Default: $\text{'N'}$
• if ${\mathbf{trans}}=\text{'T'}$, $X$, ${U}_{1}$, ${U}_{2}$, ${V}_{1}^{\mathrm{T}}$ and ${V}_{2}^{\mathrm{T}}$ are stored in row-major order;
• otherwise, $X$, ${U}_{1}$, ${U}_{2}$, ${V}_{1}^{\mathrm{T}}$ and ${V}_{2}^{\mathrm{T}}$ are stored in column-major order.
6:     $\mathrm{signs}$ – string (length ≥ 1)
Default: $\text{'D'}$
• if ${\mathbf{signs}}=\text{'O'}$, the lower-left block is made nonpositive (the other convention);
• otherwise, the upper-right block is made nonpositive (the default convention).

### Output Parameters

1:     $\mathrm{x11}\left(\mathit{ldx11},:\right)$ – double array
The first dimension, $\mathit{ldx11}$, of the array x11 will be
• if ${\mathbf{trans}}=\text{'T'}$, $\mathit{ldx11}=\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{q}}\right)$;
• otherwise $\mathit{ldx11}=\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{p}}\right)$.
The second dimension of the array x11 will be $\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{p}}\right)$ if ${\mathbf{trans}}=\text{'T'}$ and ${\mathbf{q}}$ otherwise.
Contains details of the orthogonal matrix used in a simultaneous bidiagonalization process.
2:     $\mathrm{x12}\left(\mathit{ldx12},:\right)$ – double array
The first dimension, $\mathit{ldx12}$, of the array x12 will be
• if ${\mathbf{trans}}=\text{'T'}$, $\mathit{ldx12}=\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{m}}-{\mathbf{q}}\right)$;
• otherwise $\mathit{ldx12}=\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{p}}\right)$.
The second dimension of the array x12 will be $\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{p}}\right)$ if ${\mathbf{trans}}=\text{'T'}$ and ${\mathbf{m}}-{\mathbf{q}}$ otherwise.
Contains details of the orthogonal matrix used in a simultaneous bidiagonalization process.
3:     $\mathrm{x21}\left(\mathit{ldx21},:\right)$ – double array
The first dimension, $\mathit{ldx21}$, of the array x21 will be
• if ${\mathbf{trans}}=\text{'T'}$, $\mathit{ldx21}=\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{q}}\right)$;
• otherwise $\mathit{ldx21}=\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{m}}-{\mathbf{p}}\right)$.
The second dimension of the array x21 will be $\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{m}}-{\mathbf{p}}\right)$ if ${\mathbf{trans}}=\text{'T'}$ and ${\mathbf{q}}$ otherwise.
Contains details of the orthogonal matrix used in a simultaneous bidiagonalization process.
4:     $\mathrm{x22}\left(\mathit{ldx22},:\right)$ – double array
The first dimension, $\mathit{ldx22}$, of the array x22 will be
• if ${\mathbf{trans}}=\text{'T'}$, $\mathit{ldx22}=\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{m}}-{\mathbf{q}}\right)$;
• otherwise $\mathit{ldx22}=\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{m}}-{\mathbf{p}}\right)$.
The second dimension of the array x22 will be $\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{m}}-{\mathbf{p}}\right)$ if ${\mathbf{trans}}=\text{'T'}$ and ${\mathbf{m}}-{\mathbf{q}}$ otherwise.
Contains details of the orthogonal matrix used in a simultaneous bidiagonalization process.
5:     $\mathrm{theta}\left(\mathrm{min}\phantom{\rule{0.125em}{0ex}}\left({\mathbf{p}},{\mathbf{m}}-{\mathbf{p}},{\mathbf{q}},{\mathbf{m}}-{\mathbf{q}}\right)\right)$ – double array
The values ${\theta }_{i}$ for $i=1,2,\dots ,r$ where $r=\mathrm{min}\phantom{\rule{0.125em}{0ex}}\left(p,m-p,q,m-q\right)$. The diagonal submatrices $C$ and $S$ of ${\Sigma }_{p}$ are constructed from these values as
• $C=\mathrm{diag}\left(\mathrm{cos}\left({\mathbf{theta}}\left(1\right)\right),\dots ,\mathrm{cos}\left({\mathbf{theta}}\left(r\right)\right)\right)$ and
• $S=\mathrm{diag}\left(\mathrm{sin}\left({\mathbf{theta}}\left(1\right)\right),\dots ,\mathrm{sin}\left({\mathbf{theta}}\left(r\right)\right)\right)$.
6:     $\mathrm{u1}\left(\mathit{ldu1},:\right)$ – double array
The first dimension of the array u1 will be if ${\mathbf{jobu1}}=\text{'Y'}$, $\mathit{ldu1}=\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{p}}\right)$.
The second dimension of the array u1 will be $\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{p}}\right)$ if ${\mathbf{jobu1}}=\text{'Y'}$ and $1$ otherwise.
If ${\mathbf{jobu1}}=\text{'Y'}$, u1 contains the $p$ by $p$ orthogonal matrix ${U}_{1}$.
7:     $\mathrm{u2}\left(\mathit{ldu2},:\right)$ – double array
The first dimension of the array u2 will be if ${\mathbf{jobu2}}=\text{'Y'}$, $\mathit{ldu2}=\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{m}}-{\mathbf{p}}\right)$.
The second dimension of the array u2 will be $\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{m}}-{\mathbf{p}}\right)$ if ${\mathbf{jobu2}}=\text{'Y'}$ and $1$ otherwise.
If ${\mathbf{jobu2}}=\text{'Y'}$, u2 contains the $m-p$ by $m-p$ orthogonal matrix ${U}_{2}$.
8:     $\mathrm{v1t}\left(\mathit{ldv1t},:\right)$ – double array
The first dimension of the array v1t will be if ${\mathbf{jobv1t}}=\text{'Y'}$, $\mathit{ldv1t}=\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{q}}\right)$.
The second dimension of the array v1t will be $\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{q}}\right)$ if ${\mathbf{jobv1t}}=\text{'Y'}$ and $1$ otherwise.
If ${\mathbf{jobv1t}}=\text{'Y'}$, v1t contains the $q$ by $q$ orthogonal matrix ${{V}_{1}}^{\mathrm{T}}$.
9:     $\mathrm{v2t}\left(\mathit{ldv2t},:\right)$ – double array
The first dimension of the array v2t will be if ${\mathbf{jobv2t}}=\text{'Y'}$, $\mathit{ldv2t}=\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{m}}-{\mathbf{q}}\right)$.
The second dimension of the array v2t will be $\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{m}}-{\mathbf{q}}\right)$ if ${\mathbf{jobv2t}}=\text{'Y'}$ and $1$ otherwise.
If ${\mathbf{jobv2t}}=\text{'Y'}$, v2t contains the $m-q$ by $m-q$ orthogonal matrix ${{V}_{2}}^{\mathrm{T}}$.
10:   $\mathrm{info}$int64int32nag_int scalar
${\mathbf{info}}=0$ unless the function detects an error (see Error Indicators and Warnings).

## Error Indicators and Warnings

${\mathbf{info}}<0$
If ${\mathbf{info}}=-i$, argument $i$ had an illegal value. An explanatory message is output, and execution of the program is terminated.
${\mathbf{info}}>0$
The Jacobi-type procedure failed to converge during an internal reduction to bidiagonal-block form. The process requires convergence to $\mathrm{min}\phantom{\rule{0.125em}{0ex}}\left({\mathbf{p}},{\mathbf{m}}-{\mathbf{p}},{\mathbf{q}},{\mathbf{m}}-{\mathbf{q}}\right)$ values, the value of info gives the number of converged values.

## Accuracy

The computed $CS$ decomposition is nearly the exact $CS$ decomposition for the nearby matrix $\left(X+E\right)$, where
 $E2 = Oε ,$
and $\epsilon$ is the machine precision.

The total number of floating-point operations required to perform the full $CS$ decomposition is approximately $2{m}^{3}$.
The complex analogue of this function is nag_lapack_zuncsd (f08rn).

## Example

This example finds the full CS decomposition of
 $X = -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$
partitioned in $2$ by $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.
```function f08ra_example

fprintf('f08ra example results\n\n');

% Find the full CS decomposition of X:
m = int64(5);
p = int64(3);
q = int64(2);
x = [-0.7576  0.3697  0.3838  0.2126 -0.3112;
-0.4077 -0.1552 -0.1129  0.2676  0.8517;
-0.0488  0.7240 -0.6730 -0.1301  0.0602;
-0.2287  0.0088  0.2235 -0.9235  0.2120;
0.4530  0.5612  0.5806  0.1162  0.3595];

% Partition matrix
x11 = x(1:p,  1:q);   x12 = x(1:p,  q+1:m);
x21 = x(p+1:m,1:q);   x22 = x(p+1:m,q+1:m);

[x11, x12, x21, x22, theta, u1, u2, v1t, v2t, info] = ...
f08ra( ...
m, p, q, x11, x12, x21, x22);

disp('Components of CS factorization of X:');
disp('    Theta');
disp(theta);
disp('    U1');
disp(u1);
disp('    U2');
disp(u2);
disp('    V1');
disp(v1t');
disp('    V2');
disp(v2t');

```
```f08ra example results

Components of CS factorization of X:
Theta
0.1811
0.8255

U1
-0.8249   -0.3370   -0.4538
-0.2042   -0.5710    0.7952
-0.5271    0.7486    0.4022

U2
-0.9802   -0.1982
-0.1982    0.9802

V1
0.7461    0.6658
-0.6658    0.7461

V2
-0.3397    0.7738    0.5346
0.8967    0.4379   -0.0640
-0.2837    0.4576   -0.8427

```