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_zggsvp (f08vs)

## Purpose

nag_lapack_zggsvp (f08vs) uses unitary transformations to simultaneously reduce the $m$ by $n$ matrix $A$ and the $p$ by $n$ matrix $B$ to upper triangular form. This factorization is usually used as a preprocessing step for computing the generalized singular value decomposition (GSVD).

## Syntax

[a, b, k, l, u, v, q, info] = f08vs(jobu, jobv, jobq, a, b, tola, tolb, 'm', m, 'p', p, 'n', n)
[a, b, k, l, u, v, q, info] = nag_lapack_zggsvp(jobu, jobv, jobq, a, b, tola, tolb, 'm', m, 'p', p, 'n', n)

## Description

nag_lapack_zggsvp (f08vs) computes unitary matrices $U$, $V$ and $Q$ such that
where the $k$ by $k$ matrix ${A}_{12}$ and $l$ by $l$ matrix ${B}_{13}$ are nonsingular upper triangular; ${A}_{23}$ is $l$ by $l$ upper triangular if $m-k-l\ge 0$ and is $\left(m-k\right)$ by $l$ upper trapezoidal otherwise. $\left(k+l\right)$ is the effective numerical rank of the $\left(m+p\right)$ by $n$ matrix ${\left(\begin{array}{cc}{A}^{\mathrm{H}}& {B}^{\mathrm{H}}\end{array}\right)}^{\mathrm{H}}$.
This decomposition is usually used as the preprocessing step for computing the Generalized Singular Value Decomposition (GSVD), see function nag_lapack_zggsvd (f08vn).

## 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:     $\mathrm{jobu}$ – string (length ≥ 1)
If ${\mathbf{jobu}}=\text{'U'}$, the unitary matrix $U$ is computed.
If ${\mathbf{jobu}}=\text{'N'}$, $U$ is not computed.
Constraint: ${\mathbf{jobu}}=\text{'U'}$ or $\text{'N'}$.
2:     $\mathrm{jobv}$ – string (length ≥ 1)
If ${\mathbf{jobv}}=\text{'V'}$, the unitary matrix $V$ is computed.
If ${\mathbf{jobv}}=\text{'N'}$, $V$ is not computed.
Constraint: ${\mathbf{jobv}}=\text{'V'}$ or $\text{'N'}$.
3:     $\mathrm{jobq}$ – string (length ≥ 1)
If ${\mathbf{jobq}}=\text{'Q'}$, the unitary matrix $Q$ is computed.
If ${\mathbf{jobq}}=\text{'N'}$, $Q$ is not computed.
Constraint: ${\mathbf{jobq}}=\text{'Q'}$ or $\text{'N'}$.
4:     $\mathrm{a}\left(\mathit{lda},:\right)$ – complex array
The first dimension of the array a must be at least $\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{m}}\right)$.
The second dimension of the array a must be at least $\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{n}}\right)$.
The $m$ by $n$ matrix $A$.
5:     $\mathrm{b}\left(\mathit{ldb},:\right)$ – complex array
The first dimension of the array b must be at least $\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{p}}\right)$.
The second dimension of the array b must be at least $\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{n}}\right)$.
The $p$ by $n$ matrix $B$.
6:     $\mathrm{tola}$ – double scalar
7:     $\mathrm{tolb}$ – double scalar
tola and tolb are the thresholds to determine the effective numerical rank of matrix $B$ and a subblock of $A$. Generally, they are set to
 $tola=maxm,nAε, tolb=maxp,nBε,$
where $\epsilon$ is the machine precision.
The size of tola and tolb may affect the size of backward errors of the decomposition.

### Optional Input Parameters

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

### Output Parameters

1:     $\mathrm{a}\left(\mathit{lda},:\right)$ – complex array
The first dimension of the array a will be $\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{m}}\right)$.
The second dimension of the array a will be $\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{n}}\right)$.
Contains the triangular (or trapezoidal) matrix described in Description.
2:     $\mathrm{b}\left(\mathit{ldb},:\right)$ – complex array
The first dimension of the array b will be $\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{p}}\right)$.
The second dimension of the array b will be $\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{n}}\right)$.
Contains the triangular matrix described in Description.
3:     $\mathrm{k}$int64int32nag_int scalar
4:     $\mathrm{l}$int64int32nag_int scalar
k and l specify the dimension of the subblocks $k$ and $l$ as described in Description; $\left(k+l\right)$ is the effective numerical rank of ${\left(\begin{array}{cc}{{\mathbf{a}}}^{\mathrm{T}}& {{\mathbf{b}}}^{\mathrm{T}}\end{array}\right)}^{\mathrm{T}}$.
5:     $\mathrm{u}\left(\mathit{ldu},:\right)$ – complex array
The first dimension, $\mathit{ldu}$, of the array u will be
• if ${\mathbf{jobu}}=\text{'U'}$, $\mathit{ldu}=\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{m}}\right)$;
• otherwise $\mathit{ldu}=1$.
The second dimension of the array u will be $\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{m}}\right)$ if ${\mathbf{jobu}}=\text{'U'}$ and $1$ otherwise.
If ${\mathbf{jobu}}=\text{'U'}$, u contains the unitary matrix $U$.
If ${\mathbf{jobu}}=\text{'N'}$, u is not referenced.
6:     $\mathrm{v}\left(\mathit{ldv},:\right)$ – complex array
The first dimension, $\mathit{ldv}$, of the array v will be
• if ${\mathbf{jobv}}=\text{'V'}$, $\mathit{ldv}=\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{p}}\right)$;
• otherwise $\mathit{ldv}=1$.
The second dimension of the array v will be $\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{p}}\right)$ if ${\mathbf{jobv}}=\text{'V'}$ and $1$ otherwise.
If ${\mathbf{jobv}}=\text{'V'}$, v contains the unitary matrix $V$.
If ${\mathbf{jobv}}=\text{'N'}$, v is not referenced.
7:     $\mathrm{q}\left(\mathit{ldq},:\right)$ – complex array
The first dimension, $\mathit{ldq}$, of the array q will be
• if ${\mathbf{jobq}}=\text{'Q'}$, $\mathit{ldq}=\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{n}}\right)$;
• otherwise $\mathit{ldq}=1$.
The second dimension of the array q will be $\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{n}}\right)$ if ${\mathbf{jobq}}=\text{'Q'}$ and $1$ otherwise.
If ${\mathbf{jobq}}=\text{'Q'}$, q contains the unitary matrix $Q$.
If ${\mathbf{jobq}}=\text{'N'}$, q is not referenced.
8:     $\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}}=-i$
If ${\mathbf{info}}=-i$, parameter $i$ had an illegal value on entry. The parameters are numbered as follows:
1: jobu, 2: jobv, 3: jobq, 4: m, 5: p, 6: n, 7: a, 8: lda, 9: b, 10: ldb, 11: tola, 12: tolb, 13: k, 14: l, 15: u, 16: ldu, 17: v, 18: ldv, 19: q, 20: ldq, 21: iwork, 22: rwork, 23: tau, 24: work, 25: 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.

## Accuracy

The computed factorization is nearly the exact factorization for nearby matrices $\left(A+E\right)$ and $\left(B+F\right)$, where
 $E2 = OεA2 and F2= OεB2,$
and $\epsilon$ is the machine precision.

The real analogue of this function is nag_lapack_dggsvp (f08ve).

## Example

This example finds the generalized factorization
 $A = UΣ1 0 S QH , B= VΣ2 0 T QH ,$
of the matrix pair $\left(\begin{array}{cc}A& B\end{array}\right)$, where
 $A = 0.96-0.81i -0.03+0.96i -0.91+2.06i -0.05+0.41i -0.98+1.98i -1.20+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.60i 0.37+0.38i 0.19-0.54i -0.98-0.36i 0.22-0.20i 0.83+0.51i 0.20+0.01i -0.17-0.46i 1.47+1.59i 1.08-0.28i 0.20-0.12i -0.07+1.23i 0.26+0.26i$
and
 $B = 10-10 010-1 .$
```function f08vs_example

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

m = 6;
n = 4;
a = [ 0.96-0.81i -0.03+0.96i -0.91+2.06i -0.05+0.41i;
-0.98+1.98i -1.20+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.60i;
0.37+0.38i  0.19-0.54i -0.98-0.36i  0.22-0.20i;
0.83+0.51i  0.20+0.01i -0.17-0.46i  1.47+1.59i;
1.08-0.28i  0.20-0.12i -0.07+1.23i  0.26+0.26i];
p = 2;
b = complex([ 1 0 -1  0;
0 1  0 -1]);

% Reduce A and B to upper triangular form S = U^H A Q, T = V^H B Q
tola = max(m,n)*norm(a,1)*x02aj;
tolb = max(p,n)*norm(a,1)*x02aj;
[S, T, k, l, U, V, Q, info] = ...
f08vs( ...
'U', 'V', 'Q', a, b, tola, tolb);

% Compute singular values
[S, T, alpha, beta, U, V, Q, ncycle, info] = ...
f08ys( ...
'U', 'V', 'Q', k, l, S, T, tola, tolb, U, V, Q);

fprintf('Number of infinite generalized singular values = %3d\n',k);
fprintf('Number of   finite generalized singular values = %3d\n',l);
fprintf('Effective rank of the matrix pair  (A^T B^T)^T = %3d\n',k+l);
fprintf('Number of cycles of the Kogbetliantz method    = %3d\n\n',ncycle);
disp('Finite generalized singular values');
disp(alpha(k+1:k+l)./beta(k+1:k+l));

```
```f08vs example results

Number of infinite generalized singular values =   2
Number of   finite generalized singular values =   2
Effective rank of the matrix pair  (A^T B^T)^T =   4
Number of cycles of the Kogbetliantz method    =   2

Finite generalized singular values
2.0720
1.1058

```