NAG Toolbox: nag_lapack_dtgsja (f08ye)

Purpose

Syntax

Description

References

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 $\max \left(1,\mathbf{m}\right)$.

The second dimension of the array a must be at least $\max \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 $\max \left(1,\mathbf{p}\right)$.

The second dimension of the array b must be at least $\max \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

$$\begin{array}{c}\mathbf{tola}=\max \left(\mathbf{m},\mathbf{n}\right)‖A‖\epsilon \text{,}\\ \mathbf{tolb}=\max \left(\mathbf{p},\mathbf{n}\right)‖B‖\epsilon \text{,}\end{array}$$

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 \max \left(1,\mathbf{m}\right)$;
• otherwise $\mathit{ldu}\ge 1$.

The second dimension of the array u must be at least $\max \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 \max \left(1,\mathbf{p}\right)$;
• otherwise $\mathit{ldv}\ge 1$.

The second dimension of the array v must be at least $\max \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 \max \left(1,\mathbf{n}\right)$;
• otherwise $\mathit{ldq}\ge 1$.

The second dimension of the array q must be at least $\max \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 $\max \left(1,\mathbf{m}\right)$.

The second dimension of the array a will be $\max \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 $\max \left(1,\mathbf{p}\right)$.

The second dimension of the array b will be $\max \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<n$, $\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}=\max \left(1,\mathbf{m}\right)$;
• otherwise $\mathit{ldu}=1$.

The second dimension of the array u will be $\max \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}=\max \left(1,\mathbf{p}\right)$;
• otherwise $\mathit{ldv}=1$.

The second dimension of the array v will be $\max \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}=\max \left(1,\mathbf{n}\right)$;
• otherwise $\mathit{ldq}=1$.

The second dimension of the array q will be $\max \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

Further Comments

Example

Open in the MATLAB editor: f08ye_example

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