Integer type:  int32  int64  nag_int  show int32  show int32  show int64  show int64  show nag_int  show nag_int

# NAG Toolbox: nag_eigen_real_triang_svd (f02wu)

## Purpose

nag_eigen_real_triang_svd (f02wu) returns all, or part, of the singular value decomposition of a real upper triangular matrix.

## Syntax

[a, b, q, sv, work, ifail] = f02wu(a, b, wantq, wantp, 'n', n, 'ncolb', ncolb)
[a, b, q, sv, work, ifail] = nag_eigen_real_triang_svd(a, b, wantq, wantp, 'n', n, 'ncolb', ncolb)

## Description

The n$n$ by n$n$ upper triangular matrix R$R$ is factorized as
 R = QSPT, $R=QSPT,$
where Q$Q$ and P$P$ are n$n$ by n$n$ orthogonal matrices and S$S$ is an n$n$ by n$n$ diagonal matrix with non-negative diagonal elements, σ1,σ2,,σn${\sigma }_{1},{\sigma }_{2},\dots ,{\sigma }_{n}$, ordered such that
 σ1 ≥ σ2 ≥ … ≥ σn ≥ 0. $σ1≥σ2≥…≥σn≥0.$
The columns of Q$Q$ are the left-hand singular vectors of R$R$, the diagonal elements of S$S$ are the singular values of R$R$ and the columns of P$P$ are the right-hand singular vectors of R$R$.
Either or both of Q$Q$ and PT${P}^{\mathrm{T}}$ may be requested and the matrix C$C$ given by
 C = QTB, $C=QTB,$
where B$B$ is an n$n$ by ncolb$\mathit{ncolb}$ given matrix, may also be requested.
The function obtains the singular value decomposition by first reducing R$R$ to bidiagonal form by means of Givens plane rotations and then using the QR$QR$ algorithm to obtain the singular value decomposition of the bidiagonal form.
Good background descriptions to the singular value decomposition are given in Chan (1982), Dongarra et al. (1979), Golub and Van Loan (1996), Hammarling (1985) and Wilkinson (1978).
Note that if K$K$ is any orthogonal diagonal matrix so that
 KKT = I $KKT=I$
(that is the diagonal elements of K$K$ are + 1$+1$ or 1$-1$) then
 A = (QK)S(PK)T $A= (QK) S (PK) T$
is also a singular value decomposition of A$A$.

## References

Chan T F (1982) An improved algorithm for computing the singular value decomposition ACM Trans. Math. Software 8 72–83
Dongarra J J, Moler C B, Bunch J R and Stewart G W (1979) LINPACK Users' Guide SIAM, Philadelphia
Golub G H and Van Loan C F (1996) Matrix Computations (3rd Edition) Johns Hopkins University Press, Baltimore
Hammarling S (1985) The singular value decomposition in multivariate statistics SIGNUM Newsl. 20(3) 2–25
Wilkinson J H (1978) Singular Value Decomposition – Basic Aspects Numerical Software – Needs and Availability (ed D A H Jacobs) Academic Press

## Parameters

### Compulsory Input Parameters

1:     a(lda, : $:$) – double array
The first dimension of the array a must be at least max (1,n)$\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{n}}\right)$
The second dimension of the array must be at least max (1,n)$\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{n}}\right)$
The leading n$n$ by n$n$ upper triangular part of the array a must contain the upper triangular matrix R$R$.
2:     b(ldb, : $:$) – double array
The first dimension, ldb, of the array b must satisfy
• if ncolb > 0${\mathbf{ncolb}}>0$, ldbmax (1,n)$\mathit{ldb}\ge \mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{n}}\right)$;
• otherwise ldb1$\mathit{ldb}\ge 1$.
The second dimension of the array must be at least max (1,ncolb)$\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{ncolb}}\right)$
With ncolb > 0${\mathbf{ncolb}}>0$, the leading n$n$ by ncolb$\mathit{ncolb}$ part of the array b must contain the matrix to be transformed.
3:     wantq – logical scalar
Must be true if the matrix Q$Q$ is required.
If wantq = false${\mathbf{wantq}}=\mathbf{false}$, the array q is not referenced.
4:     wantp – logical scalar
Must be true if the matrix PT${P}^{\mathrm{T}}$ is required, in which case PT${P}^{\mathrm{T}}$ is overwritten on the array a, otherwise wantp must be false.

### Optional Input Parameters

1:     n – int64int32nag_int scalar
Default: The first dimension of the array a The second dimension of the array a.
n$n$, the order of the matrix R$R$.
If n = 0${\mathbf{n}}=0$, an immediate return is effected.
Constraint: n0${\mathbf{n}}\ge 0$.
2:     ncolb – int64int32nag_int scalar
Default: The second dimension of the array b.
ncolb$\mathit{ncolb}$, the number of columns of the matrix B$B$.
If ncolb = 0${\mathbf{ncolb}}=0$, the array b is not referenced.
Constraint: ncolb0${\mathbf{ncolb}}\ge 0$.

lda ldb ldq

### Output Parameters

1:     a(lda, : $:$) – double array
The first dimension of the array a will be max (1,n)$\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{n}}\right)$
The second dimension of the array will be max (1,n)$\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{n}}\right)$
ldamax (1,n)$\mathit{lda}\ge \mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{n}}\right)$.
If wantp = true${\mathbf{wantp}}=\mathbf{true}$, the n$n$ by n$n$ part of a will contain the n$n$ by n$n$ orthogonal matrix PT${P}^{\mathrm{T}}$, otherwise the n$n$ by n$n$ upper triangular part of a is used as internal workspace, but the strictly lower triangular part of a is not referenced.
2:     b(ldb, : $:$) – double array
The first dimension, ldb, of the array b will be
• if ncolb > 0${\mathbf{ncolb}}>0$, ldbmax (1,n)$\mathit{ldb}\ge \mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{n}}\right)$;
• otherwise ldb1$\mathit{ldb}\ge 1$.
The second dimension of the array will be max (1,ncolb)$\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{ncolb}}\right)$
The leading n$n$ by ncolb$\mathit{ncolb}$ part of the array b stores the matrix QTB${Q}^{\mathrm{T}}B$.
3:     q(ldq, : $:$) – double array
The first dimension, ldq, of the array q will be
• if wantq = true${\mathbf{wantq}}=\mathbf{true}$, ldqmax (1,n)$\mathit{ldq}\ge \mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{n}}\right)$;
• otherwise ldq1$\mathit{ldq}\ge 1$.
The second dimension of the array will be max (1,n)$\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{n}}\right)$ if wantq = true${\mathbf{wantq}}=\mathbf{true}$, and at least 1$1$ otherwise
With wantq = true${\mathbf{wantq}}=\mathbf{true}$, the leading n$n$ by n$n$ part of the array q will contain the orthogonal matrix Q$Q$. Otherwise the array q is not referenced.
4:     sv(n) – double array
The array sv will contain the n$n$ diagonal elements of the matrix S$S$.
5:     work( : $:$) – double array
Note: the dimension of the array work must be at least max (1,2 × (n1))$\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,2×\left({\mathbf{n}}-1\right)\right)$ if ncolb = 0${\mathbf{ncolb}}=0$ and wantq = false${\mathbf{wantq}}=\mathbf{false}$ and wantp = false${\mathbf{wantp}}=\mathbf{false}$, max (1,3 × (n1))$\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,3×\left({\mathbf{n}}-1\right)\right)$ if (ncolb = 0${\mathbf{ncolb}}=0$ and wantq = false${\mathbf{wantq}}=\mathbf{false}$ and wantp = true${\mathbf{wantp}}=\mathbf{true}$) or (wantp = false${\mathbf{wantp}}=\mathbf{false}$ and (ncolb > 0${\mathbf{ncolb}}>0$ or wantq = true${\mathbf{wantq}}=\mathbf{true}$)), and at least max (1,5 × (n1))$\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,5×\left({\mathbf{n}}-1\right)\right)$ otherwise.
work(n)${\mathbf{work}}\left({\mathbf{n}}\right)$ contains the total number of iterations taken by the QR$QR$ algorithm.
The rest of the array is used as internal workspace.
6:     ifail – int64int32nag_int scalar
${\mathrm{ifail}}={\mathbf{0}}$ unless the function detects an error (see [Error Indicators and Warnings]).

## Error Indicators and Warnings

Errors or warnings detected by the function:

Cases prefixed with W are classified as warnings and do not generate an error of type NAG:error_n. See nag_issue_warnings.

ifail = 1${\mathbf{ifail}}=-1$
 On entry, n < 0${\mathbf{n}}<0$, or lda < n$\mathit{lda}<{\mathbf{n}}$, or ncolb < 0${\mathbf{ncolb}}<0$, or ldb < n$\mathit{ldb}<{\mathbf{n}}$ and ncolb > 0${\mathbf{ncolb}}>0$, or ldq < n$\mathit{ldq}<{\mathbf{n}}$ and wantq = true${\mathbf{wantq}}=\mathbf{true}$.
W ifail > 0${\mathbf{ifail}}>0$
The QR$QR$ algorithm has failed to converge in 50 × n$50×{\mathbf{n}}$ iterations. In this case sv(1),sv(2),,sv(ifail)${\mathbf{sv}}\left(1\right),{\mathbf{sv}}\left(2\right),\dots ,{\mathbf{sv}}\left({\mathbf{ifail}}\right)$ may not have been found correctly and the remaining singular values may not be the smallest. The matrix R$R$ will nevertheless have been factorized as R = QEPT$R=QE{P}^{\mathrm{T}}$, where E$E$ is a bidiagonal matrix with sv(1),sv(2),,sv(n)${\mathbf{sv}}\left(1\right),{\mathbf{sv}}\left(2\right),\dots ,{\mathbf{sv}}\left(n\right)$ as the diagonal elements and work(1),work(2),,work(n1)${\mathbf{work}}\left(1\right),{\mathbf{work}}\left(2\right),\dots ,{\mathbf{work}}\left(n-1\right)$ as the superdiagonal elements.
This failure is not likely to occur.

## Accuracy

The computed factors Q$Q$, S$S$ and P$P$ satisfy the relation
 QSPT = R + E, $QSPT=R+E,$
where
 ‖E‖ ≤ cε ‖A‖ , $‖E‖ ≤cε ‖A‖ ,$
ε$\epsilon$ is the machine precision, c$c$ is a modest function of n$n$ and . $‖.‖$ denotes the spectral (two) norm. Note that A = sv(1)$‖A‖={\mathbf{sv}}\left(1\right)$.
A similar result holds for the computed matrix QTB${Q}^{\mathrm{T}}B$.
The computed matrix Q$Q$ satisfies the relation
 Q = T + F, $Q=T+F,$
where T$T$ is exactly orthogonal and
 ‖F‖ ≤ dε, $‖F‖ ≤dε,$
where d$d$ is a modest function of n$n$. A similar result holds for P$P$.

For given values of ncolb, wantq and wantp, the number of floating point operations required is approximately proportional to n3${n}^{3}$.
>Following the use of nag_eigen_real_triang_svd (f02wu) the rank of R$R$ may be estimated as follows:
```tol = eps;
irank = 1;
while (irank <= numel(sv) && sv(irank) >= tol*sv(1) )
irank = irank + 1;
end
```
returns the value k$k$ in irank, where k$k$ is the smallest integer for which sv(k) < tol × sv(1)${\mathbf{sv}}\left(k\right)<\mathit{tol}×{\mathbf{sv}}\left(1\right)$, where tol$\mathit{tol}$ is typically the machine precision, so that irank is an estimate of the rank of S$S$ and thus also of R$R$.

## Example

```function nag_eigen_real_triang_svd_example
a = [-4, -2, -3;
0, -3, -2;
0, 0, -4];
b = [-1; -1; -1];
wantq = true;
wantp = true;
[aOut, bOut, q, sv, work, ifail] = nag_eigen_real_triang_svd(a, b, wantq, wantp)
```
```

aOut =

-0.4694   -0.4324   -0.7699
0.7845    0.1961   -0.5883
-0.4054    0.8801   -0.2471

bOut =

-1.6716
-0.3922
0.2276

q =

0.7699   -0.5883    0.2471
0.4324    0.1961   -0.8801
0.4694    0.7845    0.4054

sv =

6.5616
3.0000
2.4384

work =

0
0
1.0000
0.5547
0
-0.8321
0
0
0
0

ifail =

0

```
```function f02wu_example
a = [-4, -2, -3;
0, -3, -2;
0, 0, -4];
b = [-1; -1; -1];
wantq = true;
wantp = true;
[aOut, bOut, q, sv, work, ifail] = f02wu(a, b, wantq, wantp)
```
```

aOut =

-0.4694   -0.4324   -0.7699
0.7845    0.1961   -0.5883
-0.4054    0.8801   -0.2471

bOut =

-1.6716
-0.3922
0.2276

q =

0.7699   -0.5883    0.2471
0.4324    0.1961   -0.8801
0.4694    0.7845    0.4054

sv =

6.5616
3.0000
2.4384

work =

0
0
1.0000
0.5547
0
-0.8321
0
0
0
0

ifail =

0

```