NAG Toolbox: nag_lapack_dggsvd (f08va)

Purpose

Syntax

Description

References

Parameters

Compulsory Input Parameters

1:     $\mathrm{jobu}$ – string (length ≥ 1)

If $\mathbf{jobu}=\text{'U'}$, the orthogonal 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 orthogonal 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 orthogonal 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)$ – 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$.

5:     $\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$.

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{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$.

3:     $\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$.

Output Parameters

1:     $\mathrm{k}$ – int64int32nag_int scalar

2:     $\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}{c}A\\ B\end{array}\right)$.

3:     $\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)$.

Contains the triangular matrix $R$, or part of $R$. See Description for details.

4:     $\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)$.

Contains the triangular matrix $R$ if $m-k-l<0$. See Description for details.

5:     $\mathrm{alpha}\left(\mathbf{n}\right)$ – double array

See the description of beta.

6:     $\mathrm{beta}\left(\mathbf{n}\right)$ – double array

alpha and beta contain the generalized singular value pairs of $A$ and $B$, ${\alpha }_i$ and ${\beta }_i$;

• $\mathbf{alpha}\left(1:\mathbf{k}\right)=1$,
• $\mathbf{beta}\left(1:\mathbf{k}\right)=0$,

and if $m-k-l\ge 0$,

• $\mathbf{alpha}\left(\mathbf{k}+1:\mathbf{k}+\mathbf{l}\right)=C$,
• $\mathbf{beta}\left(\mathbf{k}+1:\mathbf{k}+\mathbf{l}\right)=S$,

or if $m-k-l<0$,

• $\mathbf{alpha}\left(\mathbf{k}+1:\mathbf{m}\right)=C$,
• $\mathbf{alpha}\left(\mathbf{m}+1:\mathbf{k}+\mathbf{l}\right)=0$,
• $\mathbf{beta}\left(\mathbf{k}+1:\mathbf{m}\right)=S$,
• $\mathbf{beta}\left(\mathbf{m}+1:\mathbf{k}+\mathbf{l}\right)=1$, and
• $\mathbf{alpha}\left(\mathbf{k}+\mathbf{l}+1:\mathbf{n}\right)=0$,
• $\mathbf{beta}\left(\mathbf{k}+\mathbf{l}+1:\mathbf{n}\right)=0$.

The notation $\mathbf{alpha}\left(\mathbf{k}:\mathbf{n}\right)$ above refers to consecutive elements $\mathbf{alpha}\left(\mathit{i}\right)$, for $\mathit{i}=\mathbf{k},\dots ,\mathbf{n}$.

7:     $\mathrm{u}\left(\mathit{ldu},:\right)$ – double array

The first dimension, $\mathit{ldu}$, of the array u will be

• if $\mathbf{jobu}=\text{'U'}$, $\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'}$ and $1$ otherwise.

If $\mathbf{jobu}=\text{'U'}$, u contains the $m$ by $m$ orthogonal matrix $U$.

If $\mathbf{jobu}=\text{'N'}$, u is not referenced.

8:     $\mathrm{v}\left(\mathit{ldv},:\right)$ – double array

The first dimension, $\mathit{ldv}$, of the array v will be

• if $\mathbf{jobv}=\text{'V'}$, $\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'}$ and $1$ otherwise.

If $\mathbf{jobv}=\text{'V'}$, v contains the $p$ by $p$ orthogonal matrix $V$.

If $\mathbf{jobv}=\text{'N'}$, v is not referenced.

9:     $\mathrm{q}\left(\mathit{ldq},:\right)$ – double array

The first dimension, $\mathit{ldq}$, of the array q will be

• if $\mathbf{jobq}=\text{'Q'}$, $\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'}$ and $1$ otherwise.

If $\mathbf{jobq}=\text{'Q'}$, q contains the $n$ by $n$ orthogonal matrix $Q$.

If $\mathbf{jobq}=\text{'N'}$, q is not referenced.

10:   $\mathrm{iwork}\left(\mathbf{n}\right)$ – 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 $\mathbf{alpha}\left(1\right)\ge \mathbf{alpha}\left(2\right)\ge \cdots \ge \mathbf{alpha}\left(\mathbf{n}\right)$.

11:   $\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: 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: iwork, 23: 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$

If $\mathbf{info}=1$, the Jacobi-type procedure failed to converge.

Accuracy

Further Comments

Example

Open in the MATLAB editor: f08va_example

function f08va_example


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

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

% Compute the GSVD of matrix pair A, B
% (a, b) = (U*D1, V*D2)*(0 R)*(Q^T), m>=n
jobu = 'U';
jobv = 'V';
jobq = 'Q';
[k, l, DR, b, alpha, beta, U, V, Q, iwork, info] = ...
  f08va( ...
	 jobu, jobv, jobq, a, b);

irank = k + l;
fprintf('\nInfinite generalized singular values = %5d\n', k);
fprintf('  Finite generalized singular values = %5d\n', l);
fprintf('  Numerical rank of (a'' b'')''         = %5d\n', irank);

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

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');
R = DR(1:n,n-(irank-1):n);
disp(R);

% Estimate the reciprocal condition number of R
[rcond, info] = f07tg( ...
		       'Inf-norm','Upper','Non-unit', R);

fprintf('Condition number for R:\n%11.1e\n', 1/rcond);

% So long as irank = n, compute error bound for singular values
if (irank==n)
  serrbd = x02aj/rcond;
  fprintf('\nError bound for the generalized singular values\n%11.1e\n', ...
          serrbd);
else
  fprintf('\n(A'' B'')'' is not of full rank\n');
end

f08va example results


Infinite generalized singular values =     1
  Finite generalized singular values =     2
  Numerical rank of (a' b')'         =     3

Finite generalized singular values
    1.3151    0.0802

Orthogonal matrix U
   -0.1348    0.5252   -0.2092    0.8137
    0.6742   -0.5221   -0.3889    0.3487
    0.2697    0.5276   -0.6578   -0.4650
    0.6742    0.4161    0.6101    0.0000

Orthogonal matrix V
    0.3554   -0.9347
    0.9347    0.3554

Orthogonal matrix Q
   -0.8321   -0.0946   -0.5466
    0.5547   -0.1419   -0.8199
         0   -0.9853    0.1706

Non singular upper triangular matrix R
   -2.0569   -9.0121   -9.3705
         0  -10.8822   -7.2688
         0         0   -6.0405

Condition number for R:
    2.4e+01

Error bound for the generalized singular values
    2.6e-15