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_dggsvd (f08va)

## Purpose

nag_lapack_dggsvd (f08va) computes the generalized singular value decomposition (GSVD) of an $m$ by $n$ real matrix $A$ and a $p$ by $n$ real matrix $B$.

## Syntax

[k, l, a, b, alpha, beta, u, v, q, iwork, info] = f08va(jobu, jobv, jobq, a, b, 'm', m, 'n', n, 'p', p)
[k, l, a, b, alpha, beta, u, v, q, iwork, info] = nag_lapack_dggsvd(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:
 At Mark 23: iwork was made an output parameter

## Description

The generalized singular value decomposition is given by
 $UT A Q = D1 0 R , VT B Q = D2 0 R ,$
where $U$, $V$ and $Q$ are orthogonal matrices. Let $\left(k+l\right)$ be the effective numerical rank of the matrix $\left(\begin{array}{c}A\\ B\end{array}\right)$, then $R$ is a $\left(k+l\right)$ by $\left(k+l\right)$ nonsingular upper triangular matrix, ${D}_{1}$ and ${D}_{2}$ are $m$ by $\left(k+l\right)$ and $p$ by $\left(k+l\right)$ ‘diagonal’ matrices structured as follows:
if $m-k-l\ge 0$,
 $D1= klkI0l0Cm-k-l00()$
 $D2= kll0Sp-l00()$
 $0R = n-k-lklk0R11R12l00R22()$
where
 $C = diagαk+1,…,αk+l ,$
 $S = diagβk+1,…,βk+l ,$
and
 $C2 + S2 = I .$
$R$ is stored as a submatrix of $A$ with elements ${R}_{ij}$ stored as ${A}_{i,n-k-l+j}$ on exit.
If $m-k-l<0$,
 $D1= km-kk+l-mkI00m-k0C0()$
 $D2= km-kk+l-mm-k0S0k+l-m00Ip-l000()$
 $0R = n-k-lkm-kk+l-mk0R11R12R13m-k00R22R23k+l-m000R33()$
where
 $C = diagαk+1,…,αm ,$
 $S = diagβk+1,…,βm ,$
and
 $C2 + S2 = I .$
$\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$ with ${R}_{ij}$ stored as ${A}_{i,n-k-l+j}$, and ${R}_{33}$ is stored as a submatrix of $B$ with ${\left({R}_{33}\right)}_{ij}$ stored as ${B}_{m-k+i,n+m-k-l+j}$.
The function computes $C$, $S$, $R$ and, optionally, the orthogonal transformation matrices $U$, $V$ and $Q$.
In particular, if $B$ is an $n$ by $n$ nonsingular matrix, then the GSVD of $A$ and $B$ implicitly gives the SVD of $A{B}^{-1}$:
 $A B-1 = U D1 D2-1 VT .$
If $\left(\begin{array}{c}A\\ B\end{array}\right)$ has orthonormal columns, then the GSVD of $A$ and $B$ is also equal to the CS decomposition of $A$ and $B$. Furthermore, the GSVD can be used to derive the solution of the eigenvalue problem:
 $AT Ax=λ BT Bx .$
In some literature, the GSVD of $A$ and $B$ is presented in the form
 $UT A X = 0D1 , VT B X = 0D2 ,$
where $U$ and $V$ are orthogonal and $X$ is nonsingular, and ${D}_{1}$ and ${D}_{2}$ are ‘diagonal’. The former GSVD form can be converted to the latter form by taking the nonsingular matrix $X$ as
 $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:     $\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 $\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)$ – 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$.

### 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 $\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 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 $\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 $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}=\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 $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}=\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 $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}=\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 $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

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

## Example

This example finds the generalized singular value decomposition
 $A = U Σ1 0R QT , B = V Σ2 0R QT ,$
where
 $A = 1 2 3 3 2 1 4 5 6 7 8 8 and B = -2 -3 3 4 6 5 ,$
together with estimates for the condition number of $R$ and the error bound for the computed generalized singular values.
The example program assumes that $m\ge n$, and would need slight modification if this is not the case.
```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
```