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_zuncsd (f08rn)

## Purpose

nag_lapack_zuncsd (f08rn) computes the CS decomposition of a complex $m$ by $m$ unitary matrix $X$, partitioned into a $2$ by $2$ array of submatrices.

## Syntax

[x11, x12, x21, x22, theta, u1, u2, v1t, v2t, info] = f08rn(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_zuncsd(m, p, q, x11, x12, x21, x22, 'jobu1', jobu1, 'jobu2', jobu2, 'jobv1t', jobv1t, 'jobv2t', jobv2t, 'trans', trans, 'signs', signs)

## Description

The $m$ by $m$ unitary 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 a unitary matrix containing the $p$ by $p$ unitary matrix ${U}_{1}$ and the $\left(m-p\right)$ by $\left(m-p\right)$ unitary matrix ${U}_{2}$;
 $V= V1 0 0 V2$
is a unitary matrix containing the $q$ by $q$ unitary matrix ${V}_{1}$ and the $\left(m-q\right)$ by $\left(m-q\right)$ unitary 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 unitary 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)$ – complex 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 unitary matrix $X$ whose CSD is desired.
5:     $\mathrm{x12}\left(\mathit{ldx12},:\right)$ – complex 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 unitary matrix $X$ whose CSD is desired.
6:     $\mathrm{x21}\left(\mathit{ldx21},:\right)$ – complex 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 unitary matrix $X$ whose CSD is desired.
7:     $\mathrm{x22}\left(\mathit{ldx22},:\right)$ – complex 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 unitary 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)$ – complex 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 unitary matrix used in a simultaneous bidiagonalization process.
2:     $\mathrm{x12}\left(\mathit{ldx12},:\right)$ – complex 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 unitary matrix used in a simultaneous bidiagonalization process.
3:     $\mathrm{x21}\left(\mathit{ldx21},:\right)$ – complex 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 unitary matrix used in a simultaneous bidiagonalization process.
4:     $\mathrm{x22}\left(\mathit{ldx22},:\right)$ – complex 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 unitary 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)$ – complex 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$ unitary matrix ${U}_{1}$.
7:     $\mathrm{u2}\left(\mathit{ldu2},:\right)$ – complex 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$ unitary matrix ${U}_{2}$.
8:     $\mathrm{v1t}\left(\mathit{ldv1t},:\right)$ – complex 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$ unitary matrix ${{V}_{1}}^{\mathrm{H}}$.
9:     $\mathrm{v2t}\left(\mathit{ldv2t},:\right)$ – complex 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$ unitary matrix ${{V}_{2}}^{\mathrm{H}}$.
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 real analogue of this function is nag_lapack_dorcsd (f08ra).

## Example

This example finds the full CS decomposition of a unitary $6$ by $6$ matrix $X$ partitioned in $3$ by $3$ blocks.
The decomposition is performed both on submatrices of the unitary matrix $X$ and on separated partition matrices. Code is also provided to perform a recombining check if required.
```function f08rn_example

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

% Find the full CS decomposition of X:
m = int64(6);
p = int64(2);
q = int64(4);

xr = [-1.3038E-02 -1.4039E-01  2.5177E-01 -5.0956E-02 -4.5947E-02 -5.2773E-02;
4.2764E-01  8.6298E-02 -3.2188E-01  1.1979E-01 -8.0311E-02 -3.8117E-02;
-3.2595E-01  3.8163E-01  1.3231E-01 -5.0671E-01  5.9714E-02 -1.3850E-01;
1.5906E-01 -2.8207E-01  2.1598E-01 -4.0163E-01 -4.6443E-02 -3.7354E-01;
-1.7210E-01 -5.0942E-01  3.6488E-02  1.9271E-01  5.7843E-01 -1.8815E-02;
-2.6336E-01 -1.0861E-01  1.0906E-01 -8.8159E-02  1.5763E-02  6.5007E-01];
xi = [-3.2595E-01 -2.6167E-01 -7.9789E-01 -2.1750E-01  1.4052E-04 -2.2492E-01;
-6.2582E-01 -3.8174E-02  1.6112E-01  1.6319E-01 -4.3605E-01 -2.1907E-01;
1.6428E-01 -1.8219E-01 -1.4565E-02  1.8615E-01 -5.8974E-01 -9.0941E-02;
-5.2151E-03  1.9732E-01  1.8813E-01  2.6787E-01  3.0864E-01 -5.5148E-01;
-1.3038E-02 -5.0319E-01  2.0316E-01  1.5574E-01 -1.2439E-01 -5.5686E-02;
-2.4772E-01  2.8474E-01 -1.2712E-01  5.6169E-01  4.7130E-02  4.9173E-03];
x = xr + i*xi;

% 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] = ...
f08rn(...
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');

```
```f08rn example results

Components of CS factorization of X:
Theta
0.1973
0.5386

U1
-0.2084 - 0.9384i   0.1565 + 0.2269i
-0.0490 + 0.2713i   0.1972 + 0.9408i

U2
0.1599 - 0.2574i  -0.0334 - 0.6270i   0.0573 + 0.6677i   0.2378 + 0.0912i
0.3473 + 0.6047i  -0.3155 - 0.0122i   0.0703 - 0.1039i   0.4497 + 0.4427i
0.3373 - 0.0221i   0.5536 - 0.1556i   0.6869 - 0.2398i  -0.1173 + 0.1096i
-0.5483 + 0.0839i   0.4188 + 0.0043i   0.0576 - 0.0505i   0.7043 - 0.1224i

V1
0.1202 + 0.0302i  -0.6762 + 0.6685i  -0.0150 - 0.1654i  -0.1297 + 0.1903i
0.2654 + 0.1007i  -0.1168 + 0.1140i  -0.3900 + 0.7376i  -0.0721 - 0.4376i
0.7707 - 0.4915i  -0.0624 - 0.1777i  -0.0175 - 0.0930i   0.3209 + 0.1306i
0.2581 + 0.0438i   0.1396 + 0.1200i   0.1007 - 0.5071i  -0.3019 - 0.7342i

V2
0.5354 + 0.0000i   0.8446 + 0.0000i
-0.8393 + 0.0941i   0.5321 - 0.0596i

```