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_zggsvd (f08vn)

## Purpose

nag_lapack_zggsvd (f08vn) computes the generalized singular value decomposition (GSVD) of an m$m$ by n$n$ complex matrix A$A$ and a p$p$ by n$n$ complex matrix B$B$.

## Syntax

[k, l, a, b, alpha, beta, u, v, q, iwork, info] = f08vn(jobu, jobv, jobq, a, b, 'm', m, 'n', n, 'p', p)
[k, l, a, b, alpha, beta, u, v, q, iwork, info] = nag_lapack_zggsvd(jobu, jobv, jobq, a, b, 'm', m, 'n', n, 'p', p)
Note: the interface to this routine has changed since earlier releases of the toolbox:
Mark 23: iwork is now an output parameter
.

## Description

The generalized singular value decomposition is given by
UH A Q = D1
 ( 0 R )
,   VH B Q = D2
 ( 0 R )
,
$UH A Q = D1 0 R , VH B Q = D2 0 R ,$
where U$U$, V$V$ and Q$Q$ are unitary matrices. Let (k + l)$\left(k+l\right)$ be the effective numerical rank of the matrix
 ( A ) B
$\left(\begin{array}{c}A\\ B\end{array}\right)$, then R$R$ is a (k + l)$\left(k+l\right)$ by (k + l)$\left(k+l\right)$ nonsingular upper triangular matrix, D1${D}_{1}$ and D2${D}_{2}$ are m$m$ by (k + l)$\left(k+l\right)$ and p$p$ by (k + l)$\left(k+l\right)$ ‘diagonal’ matrices structured as follows:
if mkl0$m-k-l\ge 0$,
D1 =
 k l k ( I 0 ) l 0 C m − k − l 0 0
$D1= klk(I0) l 0 C m-k-l 0 0$
D2 =
 k l l ( 0 S ) p − l 0 0
$D2= kll(0S) p-l 0 0$
 ( 0 R )
=
 n − k − l k l k ( 0 R11 R12 ) l 0 0 R22
$0R = n-k-lklk(0R11R12) l 0 0 R22$
where
 C = diag(αk + 1, … ,αk + l) , $C = diag(αk+1,…,αk+l) ,$
 S = diag(βk + 1, … ,βk + l) , $S = diag(βk+1,…,βk+l) ,$
and
 C2 + S2 = I . $C2 + S2 = I .$
R$R$ is stored as a submatrix of A$A$ with elements Rij${R}_{ij}$ stored as Ai,nkl + j${A}_{i,n-k-l+j}$ on exit.
If mkl < 0 $m-k-l<0$,
D1 =
 k m − k k + l − m k ( I 0 0 ) m − k 0 C 0
$D1= km-kk+l-mk(I00) m-k 0 C 0$
D2 =
 k m − k k + l − m m − k ( 0 S 0 ) k + l − m 0 0 I p − l 0 0 0
$D2= km-kk+l-mm-k(0S0) k+l-m 0 0 I p-l 0 0 0$
 ( 0 R )
=
 n − k − l k m − k k + l − m k ( 0 R11 R12 R13 ) m − k 0 0 R22 R23 k + l − m 0 0 0 R33
$0R = n-k-lkm-kk+l-mk(0R11R12R13) m-k 0 0 R22 R23 k+l-m 0 0 0 R33$
where
 C = diag(αk + 1, … ,αm) , $C = diag(αk+1,…,αm) ,$
 S = diag(βk + 1, … ,βm) , $S = diag(βk+1,…,βm) ,$
and
 C2 + S2 = I . $C2 + S2 = I .$
 ( R11 R12 R13 ) 0 R22 R23
$\left(\begin{array}{ccc}{R}_{11}& {R}_{12}& {R}_{13}\\ 0& {R}_{22}& {R}_{23}\end{array}\right)$ is stored as a submatrix of A$A$ with Rij${R}_{ij}$ stored as Ai,nkl + j${A}_{i,n-k-l+j}$, and R33 ${R}_{33}$ is stored as a submatrix of B$B$ with (R33)ij${\left({R}_{33}\right)}_{ij}$ stored as Bmk + i,n + mkl + j${B}_{m-k+i,n+m-k-l+j}$.
The function computes C$C$, S$S$, R$R$ and, optionally, the unitary transformation matrices U$U$, V$V$ and Q$Q$.
In particular, if B$B$ is an n$n$ by n$n$ nonsingular matrix, then the GSVD of A$A$ and B$B$ implicitly gives the SVD of A × B1$A×{B}^{-1}$:
 A B − 1 = U (D1D2 − 1) VH . $A B-1 = U ( D1 D2-1 ) VH .$
If
 ( A ) B
$\left(\begin{array}{c}A\\ B\end{array}\right)$ has orthonormal columns, then the GSVD of A$A$ and B$B$ is also equal to the CS decomposition of A$A$ and B$B$. Furthermore, the GSVD can be used to derive the solution of the eigenvalue problem:
 AH Ax = λ BH Bx . $AH Ax=λ BH Bx .$
In some literature, the GSVD of A$A$ and B$B$ is presented in the form
UH A X =
 ( 0 D1 )
,   VH B X =
 ( 0 D2 )
,
$UH A X = 0D1 , VH B X = 0D2 ,$
where U$U$ and V$V$ are orthogonal and X$X$ is nonsingular, and D1${D}_{1}$ and D2${D}_{2}$ are ‘diagonal’. The former GSVD form can be converted to the latter form by taking the nonsingular matrix X$X$ as
X = Q
 ( I 0 ) 0 R − 1
.
$X = Q I 0 0 R-1 .$

## 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

## Parameters

### Compulsory Input Parameters

1:     jobu – string (length ≥ 1)
If jobu = 'U'${\mathbf{jobu}}=\text{'U'}$, the unitary matrix U$U$ is computed.
If jobu = 'N'${\mathbf{jobu}}=\text{'N'}$, U$U$ is not computed.
Constraint: jobu = 'U'${\mathbf{jobu}}=\text{'U'}$ or 'N'$\text{'N'}$.
2:     jobv – string (length ≥ 1)
If jobv = 'V'${\mathbf{jobv}}=\text{'V'}$, the unitary matrix V$V$ is computed.
If jobv = 'N'${\mathbf{jobv}}=\text{'N'}$, V$V$ is not computed.
Constraint: jobv = 'V'${\mathbf{jobv}}=\text{'V'}$ or 'N'$\text{'N'}$.
3:     jobq – string (length ≥ 1)
If jobq = 'Q'${\mathbf{jobq}}=\text{'Q'}$, the unitary matrix Q$Q$ is computed.
If jobq = 'N'${\mathbf{jobq}}=\text{'N'}$, Q$Q$ is not computed.
Constraint: jobq = 'Q'${\mathbf{jobq}}=\text{'Q'}$ or 'N'$\text{'N'}$.
4:     a(lda, : $:$) – complex array
The first dimension of the array a must be at least max (1,m)$\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{m}}\right)$
The second dimension of the array must be at least max (1,n)$\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{n}}\right)$
The m$m$ by n$n$ matrix A$A$.
5:     b(ldb, : $:$) – complex array
The first dimension of the array b must be at least max (1,p)$\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{p}}\right)$
The second dimension of the array must be at least max (1,n)$\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{n}}\right)$
The p$p$ by n$n$ matrix B$B$.

### Optional Input Parameters

1:     m – int64int32nag_int scalar
Default: The first dimension of the array a.
m$m$, the number of rows of the matrix A$A$.
Constraint: m0${\mathbf{m}}\ge 0$.
2:     n – int64int32nag_int scalar
Default: The second dimension of the array a.
n$n$, the number of columns of the matrices A$A$ and B$B$.
Constraint: n0${\mathbf{n}}\ge 0$.
3:     p – int64int32nag_int scalar
Default: The first dimension of the array b.
p$p$, the number of rows of the matrix B$B$.
Constraint: p0${\mathbf{p}}\ge 0$.

### Input Parameters Omitted from the MATLAB Interface

lda ldb ldu ldv ldq work rwork

### Output Parameters

1:     k – int64int32nag_int scalar
2:     l – int64int32nag_int scalar
k and l specify the dimension of the subblocks k$k$ and l$l$ as described in Section [Description]; (k + l)$\left(k+l\right)$ is the effective numerical rank of
 ( A ) B
$\left(\begin{array}{c}A\\ B\end{array}\right)$.
3:     a(lda, : $:$) – complex array
The first dimension of the array a will be max (1,m)$\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{m}}\right)$
The second dimension of the array will be max (1,n)$\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{n}}\right)$
ldamax (1,m)$\mathit{lda}\ge \mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{m}}\right)$.
Contains the triangular matrix R$R$, or part of R$R$. See Section [Description] for details.
4:     b(ldb, : $:$) – complex array
The first dimension of the array b will be max (1,p)$\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{p}}\right)$
The second dimension of the array will be max (1,n)$\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{n}}\right)$
ldbmax (1,p)$\mathit{ldb}\ge \mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{p}}\right)$.
Contains the triangular matrix R$R$ if mkl < 0$m-k-l<0$. See Section [Description] for details.
5:     alpha(n) – double array
See the description of beta.
6:     beta(n) – double array
alpha and beta contain the generalized singular value pairs of A$A$ and B$B$, αi ${\alpha }_{i}$ and βi ${\beta }_{i}$;
• alpha(1 : k) = 1 ${\mathbf{alpha}}\left(1:{\mathbf{k}}\right)=1$,
• beta(1 : k) = 0 ${\mathbf{beta}}\left(1:{\mathbf{k}}\right)=0$,
and if mkl0 $m-k-l\ge 0$,
• alpha(k + 1 : k + l) = C ${\mathbf{alpha}}\left({\mathbf{k}}+1:{\mathbf{k}}+{\mathbf{l}}\right)=C$,
• beta(k + 1 : k + l) = S ${\mathbf{beta}}\left({\mathbf{k}}+1:{\mathbf{k}}+{\mathbf{l}}\right)=S$,
or if mkl < 0 $m-k-l<0$,
• alpha(k + 1 : m) = C ${\mathbf{alpha}}\left({\mathbf{k}}+1:{\mathbf{m}}\right)=C$,
• alpha(m + 1 : k + l) = 0 ${\mathbf{alpha}}\left({\mathbf{m}}+1:{\mathbf{k}}+{\mathbf{l}}\right)=0$,
• beta(k + 1 : m) = S ${\mathbf{beta}}\left({\mathbf{k}}+1:{\mathbf{m}}\right)=S$,
• beta(m + 1 : k + l) = 1 ${\mathbf{beta}}\left({\mathbf{m}}+1:{\mathbf{k}}+{\mathbf{l}}\right)=1$, and
• alpha(k + l + 1 : n) = 0 ${\mathbf{alpha}}\left({\mathbf{k}}+{\mathbf{l}}+1:{\mathbf{n}}\right)=0$,
• beta(k + l + 1 : n) = 0 ${\mathbf{beta}}\left({\mathbf{k}}+{\mathbf{l}}+1:{\mathbf{n}}\right)=0$.
The notation alpha(k : n)${\mathbf{alpha}}\left({\mathbf{k}}:{\mathbf{n}}\right)$ above refers to consecutive elements alpha(i)${\mathbf{alpha}}\left(\mathit{i}\right)$, for i = k,,n$\mathit{i}={\mathbf{k}},\dots ,{\mathbf{n}}$.
7:     u(ldu, : $:$) – complex array
The first dimension, ldu, of the array u will be
• if jobu = 'U'${\mathbf{jobu}}=\text{'U'}$, ldu max (1,m) $\mathit{ldu}\ge \mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{m}}\right)$;
• otherwise ldu1$\mathit{ldu}\ge 1$.
The second dimension of the array will be max (1,m)$\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{m}}\right)$ if jobu = 'U'${\mathbf{jobu}}=\text{'U'}$, and at least 1$1$ otherwise
If jobu = 'U'${\mathbf{jobu}}=\text{'U'}$, u contains the m$m$ by m$m$ unitary matrix U$U$.
If jobu = 'N'${\mathbf{jobu}}=\text{'N'}$, u is not referenced.
8:     v(ldv, : $:$) – complex array
The first dimension, ldv, of the array v will be
• if jobv = 'V'${\mathbf{jobv}}=\text{'V'}$, ldv max (1,p) $\mathit{ldv}\ge \mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{p}}\right)$;
• otherwise ldv1$\mathit{ldv}\ge 1$.
The second dimension of the array will be max (1,p)$\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{p}}\right)$ if jobv = 'V'${\mathbf{jobv}}=\text{'V'}$, and at least 1$1$ otherwise
If jobv = 'V'${\mathbf{jobv}}=\text{'V'}$, v contains the p$p$ by p$p$ unitary matrix V$V$.
If jobv = 'N'${\mathbf{jobv}}=\text{'N'}$, v is not referenced.
9:     q(ldq, : $:$) – complex array
The first dimension, ldq, of the array q will be
• if jobq = 'Q'${\mathbf{jobq}}=\text{'Q'}$, ldq max (1,n) $\mathit{ldq}\ge \mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{n}}\right)$;
• otherwise ldq1$\mathit{ldq}\ge 1$.
The second dimension of the array will be max (1,n)$\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{n}}\right)$ if jobq = 'Q'${\mathbf{jobq}}=\text{'Q'}$, and at least 1$1$ otherwise
If jobq = 'Q'${\mathbf{jobq}}=\text{'Q'}$, q contains the n$n$ by n$n$ unitary matrix Q$Q$.
If jobq = 'N'${\mathbf{jobq}}=\text{'N'}$, q is not referenced.
10:   iwork(n) – int64int32nag_int array
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)) end `
such that alpha(1)alpha(2)alpha(n)${\mathbf{alpha}}\left(1\right)\ge {\mathbf{alpha}}\left(2\right)\ge \cdots \ge {\mathbf{alpha}}\left({\mathbf{n}}\right)$.
11:   info – int64int32nag_int scalar
info = 0${\mathbf{info}}=0$ unless the function detects an error (see Section [Error Indicators and Warnings]).

## Error Indicators and Warnings

info = i${\mathbf{info}}=-i$
If info = i${\mathbf{info}}=-i$, parameter i$i$ had an illegal value on entry. The parameters are numbered as follows:
1: jobu, 2: jobv, 3: jobq, 4: m, 5: n, 6: p, 7: k, 8: l, 9: a, 10: lda, 11: b, 12: ldb, 13: alpha, 14: beta, 15: u, 16: ldu, 17: v, 18: ldv, 19: q, 20: ldq, 21: work, 22: rwork, 23: iwork, 24: info.
It is possible that info refers to a parameter that is omitted from the MATLAB interface. This usually indicates that an error in one of the other input parameters has caused an incorrect value to be inferred.
INFO = 1${\mathbf{INFO}}=1$
If info = 1${\mathbf{info}}=1$, the Jacobi-type procedure failed to converge.

## Accuracy

The computed generalized singular value decomposition is nearly the exact generalized singular value decomposition for nearby matrices (A + E) $\left(A+E\right)$ and (B + F) $\left(B+F\right)$, where
 ‖E‖2 = O(ε) ‖A‖2 ​ and ​ ‖F‖2 = O(ε) ‖B‖2 , $‖E‖2 = O(ε) ‖A‖2 ​ and ​ ‖F‖2 = O(ε) ‖B‖2 ,$
and ε $\epsilon$ is the machine precision. See Section 4.12 of Anderson et al. (1999) for further details.

The diagonal elements of the matrix R$R$ are real.
The real analogue of this function is nag_lapack_dggsvd (f08va).

## Example

```function nag_lapack_zggsvd_example
jobu = 'U';
jobv = 'V';
jobq = 'Q';
a = [ 0.96 - 0.81i,  -0.03 + 0.96i,  -0.91 + 2.06i, ...
-0.05 + 0.41i;
-0.98 + 1.98i,  -1.2 + 0.19i,  -0.66 + 0.42i, ...
-0.81 + 0.56i;
0.62 - 0.46i,  1.01 + 0.02i,  0.63 - 0.17i, ...
-1.11 + 0.6i;
0.37 + 0.38i,  0.19 - 0.54i,  -0.98 - 0.36i, ...
0.22 - 0.2i;
0.83 + 0.51i,  0.2 + 0.01i,  -0.17 - 0.46i, ...
1.47 + 1.59i;
1.08 - 0.28i,  0.2 - 0.12i,  -0.07 + 1.23i, ...
0.26 + 0.26i];
b = [complex(1),  0 + 0i,  -1 + 0i,  0 + 0i;
0 + 0i,  1 + 0i,  0 + 0i,  -1 + 0i];
% Compute the generalized singular value decomposition of (a, b)
% (a = U*D1*(0 R)*(Q^H), b = V*D2*(0 R)*(Q^H), m>=n)
[k, l, a, b, alpha, beta, u, v, q, iwork, info] = nag_lapack_zggsvd(jobu, jobv, jobq, a, b);

if (info == 0)
kk = double(k);
irank = kk + double(l);
fprintf('\nNumber of infinite generalized singular values (k)\n%5d\n', k);
fprintf('Number of finite generalized singular values (l)\n%5d\n', l);
fprintf('Numerical rank of (a^H b^H)^H (k+l)\n%5d\n', irank);

fprintf('\nFinite generalized singular values\n');
disp(transpose((alpha(kk+1:irank)./beta(kk+1:irank))));

fprintf('Orthogonal matrix u\n');
disp(u);

fprintf('Orthogonal matrix v\n');
disp(v);

fprintf('Orthogonal matrix q\n');
disp(q);

fprintf('Non singular upper triangular matrix R\n');
disp(triu(a(1:4,4-irank+1:4)));

% Estimate the reciprocal condition number of R
[rcond, info] = nag_lapack_ztrcon('Infinity-norm','Upper','Non-unit', a(1:4,4-irank+1:4));

fprintf('Estimate of reciprocal condition number for R\n%11.1e\n', rcond);

% So long as irank = n, get the machine precision, eps, from nag_machine_precision and compute
% the approximate error bound for the computed generalized singular values
if (irank==4)
serrbd = nag_machine_precision/rcond;
fprintf('\nError estimate for the generalized singular values\n%11.1e\n', serrbd);
else
fprintf('\n(A^H B^H)^H is not of full rank\n');
end
end
```
```

Number of infinite generalized singular values (k)
2
Number of finite generalized singular values (l)
2
Numerical rank of (a^H b^H)^H (k+l)
4

Finite generalized singular values
2.0720    1.1058

Orthogonal matrix u
Columns 1 through 5

-0.0130 - 0.3259i  -0.1404 - 0.2617i   0.2518 - 0.7979i  -0.0510 - 0.2175i  -0.0459 + 0.0001i
0.4276 - 0.6258i   0.0863 - 0.0382i  -0.3219 + 0.1611i   0.1198 + 0.1632i  -0.0803 - 0.4360i
-0.3259 + 0.1643i   0.3816 - 0.1822i   0.1323 - 0.0146i  -0.5067 + 0.1862i   0.0597 - 0.5897i
0.1591 - 0.0052i  -0.2821 + 0.1973i   0.2160 + 0.1881i  -0.4016 + 0.2679i  -0.0464 + 0.3086i
-0.1721 - 0.0130i  -0.5094 - 0.5032i   0.0365 + 0.2032i   0.1927 + 0.1557i   0.5784 - 0.1244i
-0.2634 - 0.2477i  -0.1086 + 0.2847i   0.1091 - 0.1271i  -0.0882 + 0.5617i   0.0158 + 0.0471i

Column 6

-0.0528 - 0.2249i
-0.0381 - 0.2191i
-0.1385 - 0.0909i
-0.3735 - 0.5515i
-0.0188 - 0.0557i
0.6501 + 0.0049i

Orthogonal matrix v
0.9893 + 0.0000i  -0.1146 + 0.0903i
-0.1146 - 0.0903i  -0.9893 + 0.0000i

Orthogonal matrix q
0.7071 + 0.0000i   0.0000 + 0.0000i   0.6995 - 0.0000i   0.0810 - 0.0638i
0.0000 + 0.0000i   0.7071 + 0.0000i  -0.0810 - 0.0638i   0.6995 + 0.0000i
0.7071 + 0.0000i   0.0000 + 0.0000i  -0.6995 + 0.0000i  -0.0810 + 0.0638i
0.0000 + 0.0000i   0.7071 + 0.0000i   0.0810 + 0.0638i  -0.6995 - 0.0000i

Non singular upper triangular matrix R
-2.7118 + 0.0000i  -1.4390 - 1.0315i  -0.0769 + 1.3613i  -0.2814 - 0.0324i
0.0000 + 0.0000i  -1.8583 + 0.0000i  -1.0760 + 0.0310i   1.3292 + 0.3677i
0.0000 + 0.0000i   0.0000 + 0.0000i   3.2537 + 0.0000i  -0.0000 - 0.0000i
0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i  -2.1084 + 0.0000i

Estimate of reciprocal condition number for R
1.3e-01

Error estimate for the generalized singular values
8.3e-16

```
```function f08vn_example
jobu = 'U';
jobv = 'V';
jobq = 'Q';
a = [ 0.96 - 0.81i,  -0.03 + 0.96i,  -0.91 + 2.06i, ...
-0.05 + 0.41i;
-0.98 + 1.98i,  -1.2 + 0.19i,  -0.66 + 0.42i, ...
-0.81 + 0.56i;
0.62 - 0.46i,  1.01 + 0.02i,  0.63 - 0.17i, ...
-1.11 + 0.6i;
0.37 + 0.38i,  0.19 - 0.54i,  -0.98 - 0.36i, ...
0.22 - 0.2i;
0.83 + 0.51i,  0.2 + 0.01i,  -0.17 - 0.46i, ...
1.47 + 1.59i;
1.08 - 0.28i,  0.2 - 0.12i,  -0.07 + 1.23i, ...
0.26 + 0.26i];
b = [complex(1),  0 + 0i,  -1 + 0i,  0 + 0i;
0 + 0i,  1 + 0i,  0 + 0i,  -1 + 0i];
% Compute the generalized singular value decomposition of (a, b)
% (a = U*D1*(0 R)*(Q^H), b = V*D2*(0 R)*(Q^H), m>=n)
[k, l, a, b, alpha, beta, u, v, q, iwork, info] = f08vn(jobu, jobv, jobq, a, b);

if (info == 0)
kk = double(k);
irank = kk + double(l);
fprintf('\nNumber of infinite generalized singular values (k)\n%5d\n', k);
fprintf('Number of finite generalized singular values (l)\n%5d\n', l);
fprintf('Numerical rank of (a^H b^H)^H (k+l)\n%5d\n', irank);

fprintf('\nFinite generalized singular values\n');
disp(transpose((alpha(kk+1:irank)./beta(kk+1:irank))));

fprintf('Orthogonal matrix u\n');
disp(u);

fprintf('Orthogonal matrix v\n');
disp(v);

fprintf('Orthogonal matrix q\n');
disp(q);

fprintf('Non singular upper triangular matrix R\n');
disp(triu(a(1:4,4-irank+1:4)));

% Estimate the reciprocal condition number of R
[rcond, info] = f07tu('Infinity-norm','Upper','Non-unit', a(1:4,4-irank+1:4));

fprintf('Estimate of reciprocal condition number for R\n%11.1e\n', rcond);

% So long as irank = n, get the machine precision, eps, from x02aj and compute
% the approximate error bound for the computed generalized singular values
if (irank==4)
serrbd = x02aj/rcond;
fprintf('\nError estimate for the generalized singular values\n%11.1e\n', serrbd);
else
fprintf('\n(A^H B^H)^H is not of full rank\n');
end
end
```
```

Number of infinite generalized singular values (k)
2
Number of finite generalized singular values (l)
2
Numerical rank of (a^H b^H)^H (k+l)
4

Finite generalized singular values
2.0720    1.1058

Orthogonal matrix u
Columns 1 through 5

-0.0130 - 0.3259i  -0.1404 - 0.2617i   0.2518 - 0.7979i  -0.0510 - 0.2175i  -0.0459 + 0.0001i
0.4276 - 0.6258i   0.0863 - 0.0382i  -0.3219 + 0.1611i   0.1198 + 0.1632i  -0.0803 - 0.4360i
-0.3259 + 0.1643i   0.3816 - 0.1822i   0.1323 - 0.0146i  -0.5067 + 0.1862i   0.0597 - 0.5897i
0.1591 - 0.0052i  -0.2821 + 0.1973i   0.2160 + 0.1881i  -0.4016 + 0.2679i  -0.0464 + 0.3086i
-0.1721 - 0.0130i  -0.5094 - 0.5032i   0.0365 + 0.2032i   0.1927 + 0.1557i   0.5784 - 0.1244i
-0.2634 - 0.2477i  -0.1086 + 0.2847i   0.1091 - 0.1271i  -0.0882 + 0.5617i   0.0158 + 0.0471i

Column 6

-0.0528 - 0.2249i
-0.0381 - 0.2191i
-0.1385 - 0.0909i
-0.3735 - 0.5515i
-0.0188 - 0.0557i
0.6501 + 0.0049i

Orthogonal matrix v
0.9893 + 0.0000i  -0.1146 + 0.0903i
-0.1146 - 0.0903i  -0.9893 + 0.0000i

Orthogonal matrix q
0.7071 + 0.0000i   0.0000 + 0.0000i   0.6995 - 0.0000i   0.0810 - 0.0638i
0.0000 + 0.0000i   0.7071 + 0.0000i  -0.0810 - 0.0638i   0.6995 + 0.0000i
0.7071 + 0.0000i   0.0000 + 0.0000i  -0.6995 + 0.0000i  -0.0810 + 0.0638i
0.0000 + 0.0000i   0.7071 + 0.0000i   0.0810 + 0.0638i  -0.6995 - 0.0000i

Non singular upper triangular matrix R
-2.7118 + 0.0000i  -1.4390 - 1.0315i  -0.0769 + 1.3613i  -0.2814 - 0.0324i
0.0000 + 0.0000i  -1.8583 + 0.0000i  -1.0760 + 0.0310i   1.3292 + 0.3677i
0.0000 + 0.0000i   0.0000 + 0.0000i   3.2537 + 0.0000i  -0.0000 - 0.0000i
0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i  -2.1084 + 0.0000i

Estimate of reciprocal condition number for R
1.3e-01

Error estimate for the generalized singular values
8.3e-16

```