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_dtgsja (f08ye)

## Purpose

nag_lapack_dtgsja (f08ye) computes the generalized singular value decomposition (GSVD) of two real upper trapezoidal matrices $A$ and $B$, where $A$ is an $m$ by $n$ matrix and $B$ is a $p$ by $n$ matrix.
$A$ and $B$ are assumed to be in the form returned by nag_lapack_dggsvp (f08ve).

## Syntax

[a, b, alpha, beta, u, v, q, ncycle, info] = f08ye(jobu, jobv, jobq, k, l, a, b, tola, tolb, u, v, q, 'm', m, 'p', p, 'n', n)
[a, b, alpha, beta, u, v, q, ncycle, info] = nag_lapack_dtgsja(jobu, jobv, jobq, k, l, a, b, tola, tolb, u, v, q, 'm', m, 'p', p, 'n', n)

## Description

nag_lapack_dtgsja (f08ye) computes the GSVD of the matrices $A$ and $B$ which are assumed to have the form as returned by nag_lapack_dggsvp (f08ve)
 $A= n-k-lklk0A12A13l00A23m-k-l000() , if ​ m-k-l ≥ 0; n-k-lklk0A12A13m-k00A23() , if ​ m-k-l < 0 ; B= n-k-lkll00B13p-l000() ,$
where the $k$ by $k$ matrix ${A}_{12}$ and the $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.
nag_lapack_dtgsja (f08ye) computes orthogonal matrices $Q$, $U$ and $V$, diagonal matrices ${D}_{1}$ and ${D}_{2}$, and an upper triangular matrix $R$ such that
 $UTAQ = D1 0 R , VTBQ = D2 0 R .$
Optionally $Q$, $U$ and $V$ may or may not be computed, or they may be premultiplied by matrices ${Q}_{1}$, ${U}_{1}$ and ${V}_{1}$ respectively.
If $\left(m-k-l\right)\ge 0$ then ${D}_{1}$, ${D}_{2}$ and $R$ have the form
 $D1= klkI0l0Cm-k-l00() ,$
 $D2= kll0Sp-l00() ,$
 $R = klkR11R12l0R22() ,$
where $C=\mathrm{diag}\left({\alpha }_{k+1},,,\dots ,,,{\alpha }_{k+l}\right)\text{, }S=\mathrm{diag}\left({\beta }_{k+1},,,\dots ,,,{\beta }_{k+l}\right)$.
If $\left(m-k-l\right)<0$ then ${D}_{1}$, ${D}_{2}$ and $R$ have the form
 $D1= km-kk+l-mkI00m-k0C0() ,$
 $D2= km-kk+l-mm-k0S0k+l-m00Ip-l000() ,$
 $R = km-kk+l-mkR11R12R13m-k0R22R23k+l-m00R33() ,$
where $C=\mathrm{diag}\left({\alpha }_{k+1},,,\dots ,,,{\alpha }_{m}\right)\text{, }S=\mathrm{diag}\left({\beta }_{k+1},,,\dots ,,,{\beta }_{m}\right)$.
In both cases the diagonal matrix $C$ has non-negative diagonal elements, the diagonal matrix $S$ has positive diagonal elements, so that $S$ is nonsingular, and ${C}^{2}+{S}^{2}=1$. See Section 2.3.5.3 of Anderson et al. (1999) for further information.

## 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'}$, u must contain an orthogonal matrix ${U}_{1}$ on entry, and the product ${U}_{1}U$ is returned.
If ${\mathbf{jobu}}=\text{'I'}$, u is initialized to the unit matrix, and the orthogonal matrix $U$ is returned.
If ${\mathbf{jobu}}=\text{'N'}$, $U$ is not computed.
Constraint: ${\mathbf{jobu}}=\text{'U'}$, $\text{'I'}$ or $\text{'N'}$.
2:     $\mathrm{jobv}$ – string (length ≥ 1)
If ${\mathbf{jobv}}=\text{'V'}$, v must contain an orthogonal matrix ${V}_{1}$ on entry, and the product ${V}_{1}V$ is returned.
If ${\mathbf{jobv}}=\text{'I'}$, v is initialized to the unit matrix, and the orthogonal matrix $V$ is returned.
If ${\mathbf{jobv}}=\text{'N'}$, $V$ is not computed.
Constraint: ${\mathbf{jobv}}=\text{'V'}$, $\text{'I'}$ or $\text{'N'}$.
3:     $\mathrm{jobq}$ – string (length ≥ 1)
If ${\mathbf{jobq}}=\text{'Q'}$, q must contain an orthogonal matrix ${Q}_{1}$ on entry, and the product ${Q}_{1}Q$ is returned.
If ${\mathbf{jobq}}=\text{'I'}$, q is initialized to the unit matrix, and the orthogonal matrix $Q$ is returned.
If ${\mathbf{jobq}}=\text{'N'}$, $Q$ is not computed.
Constraint: ${\mathbf{jobq}}=\text{'Q'}$, $\text{'I'}$ or $\text{'N'}$.
4:     $\mathrm{k}$int64int32nag_int scalar
5:     $\mathrm{l}$int64int32nag_int scalar
k and l specify the sizes, $k$ and $l$, of the subblocks of $A$ and $B$, whose GSVD is to be computed by nag_lapack_dtgsja (f08ye).
6:     $\mathrm{a}\left(\mathit{lda},:\right)$ – double 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$.
7:     $\mathrm{b}\left(\mathit{ldb},:\right)$ – double 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$.
8:     $\mathrm{tola}$ – double scalar
9:     $\mathrm{tolb}$ – double scalar
tola and tolb are the convergence criteria for the Jacobi–Kogbetliantz iteration procedure. Generally, they should be the same as used in the preprocessing step performed by nag_lapack_zggsvp (f08vs), say
 $tola=maxm,nAε, tolb=maxp,nBε,$
where $\epsilon$ is the machine precision.
10:   $\mathrm{u}\left(\mathit{ldu},:\right)$ – double array
The first dimension, $\mathit{ldu}$, of the array u must satisfy
• if ${\mathbf{jobu}}=\text{'U'}$ or $\text{'I'}$, $\mathit{ldu}\ge \mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{m}}\right)$;
• otherwise $\mathit{ldu}\ge 1$.
The second dimension of the array u must be at least $\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{m}}\right)$ if ${\mathbf{jobu}}=\text{'U'}$ or $\text{'I'}$, and at least $1$ otherwise.
If ${\mathbf{jobu}}=\text{'U'}$, u must contain an $m$ by $m$ matrix ${U}_{1}$ (usually the orthogonal matrix returned by nag_lapack_dggsvp (f08ve)).
11:   $\mathrm{v}\left(\mathit{ldv},:\right)$ – double array
The first dimension, $\mathit{ldv}$, of the array v must satisfy
• if ${\mathbf{jobv}}=\text{'V'}$ or $\text{'I'}$, $\mathit{ldv}\ge \mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{p}}\right)$;
• otherwise $\mathit{ldv}\ge 1$.
The second dimension of the array v must be at least $\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{p}}\right)$ if ${\mathbf{jobv}}=\text{'V'}$ or $\text{'I'}$, and at least $1$ otherwise.
If ${\mathbf{jobv}}=\text{'V'}$, v must contain an $p$ by $p$ matrix ${V}_{1}$ (usually the orthogonal matrix returned by nag_lapack_dggsvp (f08ve)).
12:   $\mathrm{q}\left(\mathit{ldq},:\right)$ – double array
The first dimension, $\mathit{ldq}$, of the array q must satisfy
• if ${\mathbf{jobq}}=\text{'Q'}$ or $\text{'I'}$, $\mathit{ldq}\ge \mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{n}}\right)$;
• otherwise $\mathit{ldq}\ge 1$.
The second dimension of the array q must be at least $\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{n}}\right)$ if ${\mathbf{jobq}}=\text{'Q'}$ or $\text{'I'}$, and at least $1$ otherwise.
If ${\mathbf{jobq}}=\text{'Q'}$, q must contain an $n$ by $n$ matrix ${Q}_{1}$ (usually the orthogonal matrix returned by nag_lapack_dggsvp (f08ve)).

### 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)$ – double 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)$.
If $m-k-l\ge 0$, ${\mathbf{a}}\left(1:k+l,n-k-l+1:n\right)$ contains the $\left(k+l\right)$ by $\left(k+l\right)$ upper triangular matrix $R$.
If $m-k-l<0$, ${\mathbf{a}}\left(1:m,n-k-l+1:n\right)$ contains the first $m$ rows of the $\left(k+l\right)$ by $\left(k+l\right)$ upper triangular matrix $R$, and the submatrix ${R}_{33}$ is returned in ${\mathbf{b}}\left(m-k+1:l,n+m-k-l+1:n\right)$.
2:     $\mathrm{b}\left(\mathit{ldb},:\right)$ – double 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)$.
If $m-k-l<0$, ${\mathbf{b}}\left(m-k+1:l,n+m-k-l+1:n\right)$ contains the submatrix ${R}_{33}$ of $R$.
3:     $\mathrm{alpha}\left({\mathbf{n}}\right)$ – double array
See the description of beta.
4:     $\mathrm{beta}\left({\mathbf{n}}\right)$ – double array
alpha and beta contain the generalized singular value pairs of $A$ and $B$;
• ${\mathbf{alpha}}\left(\mathit{i}\right)=1$, ${\mathbf{beta}}\left(\mathit{i}\right)=0$, for $\mathit{i}=1,2,\dots ,k$, and
• if $m-k-l\ge 0$, ${\mathbf{alpha}}\left(\mathit{i}\right)={\alpha }_{\mathit{i}}$, ${\mathbf{beta}}\left(\mathit{i}\right)={\beta }_{\mathit{i}}$, for $\mathit{i}=k+1,\dots ,k+l$, or
• if $m-k-l<0$, ${\mathbf{alpha}}\left(\mathit{i}\right)={\alpha }_{\mathit{i}}$, ${\mathbf{beta}}\left(\mathit{i}\right)={\beta }_{\mathit{i}}$, for $\mathit{i}=k+1,\dots ,m$ and ${\mathbf{alpha}}\left(\mathit{i}\right)=0$, ${\mathbf{beta}}\left(\mathit{i}\right)=1$, for $\mathit{i}=m+1,\dots ,k+l$.
Furthermore, if $k+l, ${\mathbf{alpha}}\left(\mathit{i}\right)={\mathbf{beta}}\left(\mathit{i}\right)=0$, for $\mathit{i}=k+l+1,\dots ,n$.
5:     $\mathrm{u}\left(\mathit{ldu},:\right)$ – double array
The first dimension, $\mathit{ldu}$, of the array u will be
• if ${\mathbf{jobu}}=\text{'U'}$ or $\text{'I'}$, $\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'}$ or $\text{'I'}$ and $1$ otherwise.
If ${\mathbf{jobu}}=\text{'U'}$, u contains the product ${U}_{1}U$.
If ${\mathbf{jobu}}=\text{'I'}$, u contains the orthogonal matrix $U$.
If ${\mathbf{jobu}}=\text{'N'}$, u is not referenced.
6:     $\mathrm{v}\left(\mathit{ldv},:\right)$ – double array
The first dimension, $\mathit{ldv}$, of the array v will be
• if ${\mathbf{jobv}}=\text{'V'}$ or $\text{'I'}$, $\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'}$ or $\text{'I'}$ and $1$ otherwise.
If ${\mathbf{jobv}}=\text{'I'}$, v contains the orthogonal matrix $V$.
If ${\mathbf{jobv}}=\text{'V'}$, v contains the product ${V}_{1}V$.
If ${\mathbf{jobv}}=\text{'N'}$, v is not referenced.
7:     $\mathrm{q}\left(\mathit{ldq},:\right)$ – double array
The first dimension, $\mathit{ldq}$, of the array q will be
• if ${\mathbf{jobq}}=\text{'Q'}$ or $\text{'I'}$, $\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'}$ or $\text{'I'}$ and $1$ otherwise.
If ${\mathbf{jobq}}=\text{'I'}$, q contains the orthogonal matrix $Q$.
If ${\mathbf{jobq}}=\text{'Q'}$, q contains the product ${Q}_{1}Q$.
If ${\mathbf{jobq}}=\text{'N'}$, q is not referenced.
8:     $\mathrm{ncycle}$int64int32nag_int scalar
The number of cycles required for convergence.
9:     $\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: k, 8: l, 9: a, 10: lda, 11: b, 12: ldb, 13: tola, 14: tolb, 15: alpha, 16: beta, 17: u, 18: ldu, 19: v, 20: ldv, 21: q, 22: ldq, 23: work, 24: ncycle, 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.
${\mathbf{info}}=1$
The procedure does not converge after $40$ cycles.

## Accuracy

The computed generalized singular value decomposition is nearly the exact generalized singular value decomposition 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. See Section 4.12 of Anderson et al. (1999) for further details.

The complex analogue of this function is nag_lapack_ztgsja (f08ys).

## Example

This example finds the generalized singular value decomposition
 $A = UΣ1 0 R QT , B= VΣ2 0 R QT ,$
of the matrix pair $\left(A,B\right)$, where
 $A = 1 2 3 3 2 1 4 5 6 7 8 8 and B= -2 -3 3 4 6 5 .$
```function f08ye_example

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

% Generalized singular values of (A,B) where:
m = 4;
n = 3;
a = [ 1  2  3;
3  2  1;
4  5  6;
7  8  8];
p = 2;
b = [-2 -3  3;
4  6  5];

% Reduce A and B to upper triangular form S = U^T A Q, T = V^T 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] = ...
f08ve( ...
'U', 'V', 'Q', a, b, tola, tolb);

% Compute singular values
[S, T, alpha, beta, U, V, Q, ncycle, info] = ...
f08ye( ...
'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));

```
```f08ye example results

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

Finite generalized singular values
1.3151
0.0802

```