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_rand_matrix_orthog (g05px)

## Purpose

nag_rand_matrix_orthog (g05px) generates a random orthogonal matrix.

## Syntax

[state, a, ifail] = g05px(side, init, state, a, 'm', m, 'n', n)
[state, a, ifail] = nag_rand_matrix_orthog(side, init, state, a, 'm', m, 'n', n)

## Description

nag_rand_matrix_orthog (g05px) pre- or post-multiplies an m$m$ by n$n$ matrix A$A$ by a random orthogonal matrix U$U$, overwriting A$A$. The matrix A$A$ may optionally be initialized to the identity matrix before multiplying by U$U$, hence U$U$ is returned. U$U$ is generated using the method of Stewart (1980). The algorithm can be summarised as follows.
Let x1,x2,,xn1${x}_{1},{x}_{2},\dots ,{x}_{n-1}$ follow independent multinormal distributions with zero mean and variance Iσ2$I{\sigma }^{2}$ and dimensions n,n1,,2$n,n-1,\dots ,2$; let Hj = diag(Ij1,Hj * )${H}_{j}=\mathrm{diag}\left({I}_{j-1},{H}_{j}^{*}\right)$, where Ij1${I}_{j-1}$ is the identity matrix and Hj * ${H}_{j}^{*}$ is the Householder transformation that reduces xj${x}_{j}$ to rjje1${r}_{jj}{e}_{1}$, e1${e}_{1}$ being the vector with first element one and the remaining elements zero and rjj${r}_{jj}$ being a scalar, and let D = diag(sign(r11),sign(r22),,sign(rnn))$D=\mathrm{diag}\left(\mathrm{sign}\left({r}_{11}\right),\mathrm{sign}\left({r}_{22}\right),\dots ,\mathrm{sign}\left({r}_{nn}\right)\right)$. Then the product U = DH1H2Hn1$U=D{H}_{1}{H}_{2}\dots {H}_{n-1}$ is a random orthogonal matrix distributed according to the Haar measure over the set of orthogonal matrices of n$n$. See Theorem 3.3 in Stewart (1980).
One of the initialization functions nag_rand_init_repeat (g05kf) (for a repeatable sequence if computed sequentially) or nag_rand_init_nonrepeat (g05kg) (for a non-repeatable sequence) must be called prior to the first call to nag_rand_matrix_orthog (g05px).

## References

Stewart G W (1980) The efficient generation of random orthogonal matrices with an application to condition estimates SIAM J. Numer. Anal. 17 403–409

## Parameters

### Compulsory Input Parameters

1:     side – string (length ≥ 1)
Indicates whether the matrix A$A$ is multiplied on the left or right by the random orthogonal matrix U$U$.
side = 'L'${\mathbf{side}}=\text{'L'}$
The matrix A$A$ is multiplied on the left, i.e., premultiplied.
side = 'R'${\mathbf{side}}=\text{'R'}$
The matrix A$A$ is multiplied on the right, i.e., post-multiplied.
Constraint: side = 'L'${\mathbf{side}}=\text{'L'}$ or 'R'$\text{'R'}$.
2:     init – string (length ≥ 1)
Indicates whether or not a should be initialized to the identity matrix.
init = 'I'${\mathbf{init}}=\text{'I'}$
a is initialized to the identity matrix.
init = 'N'${\mathbf{init}}=\text{'N'}$
a is not initialized and the matrix A$A$ must be supplied in a.
Constraint: init = 'I'${\mathbf{init}}=\text{'I'}$ or 'N'$\text{'N'}$.
3:     state( : $:$) – int64int32nag_int array
Note: the actual argument supplied must be the array state supplied to the initialization routines nag_rand_init_repeat (g05kf) or nag_rand_init_nonrepeat (g05kg).
Contains information on the selected base generator and its current state.
4:     a(lda,n) – double array
lda, the first dimension of the array, must satisfy the constraint ldam$\mathit{lda}\ge {\mathbf{m}}$.
If init = 'N'${\mathbf{init}}=\text{'N'}$, a must contain the matrix A$A$.

### Optional Input Parameters

1:     m – int64int32nag_int scalar
Default: The first dimension of the array a.
m$m$, the number of rows of the matrix A$A$.
Constraints:
• if side = 'L'${\mathbf{side}}=\text{'L'}$, m > 1${\mathbf{m}}>1$;
• otherwise m1${\mathbf{m}}\ge 1$.
2:     n – int64int32nag_int scalar
Default: The second dimension of the array a.
n$n$, the number of columns of the matrix A$A$.
Constraints:
• if side = 'R'${\mathbf{side}}=\text{'R'}$, n > 1${\mathbf{n}}>1$;
• otherwise n1${\mathbf{n}}\ge 1$.

lda

### Output Parameters

1:     state( : $:$) – int64int32nag_int array
Note: the actual argument supplied must be the array state supplied to the initialization routines nag_rand_init_repeat (g05kf) or nag_rand_init_nonrepeat (g05kg).
Contains updated information on the state of the generator.
2:     a(lda,n) – double array
ldam$\mathit{lda}\ge {\mathbf{m}}$.
The matrix UA$UA$ when side = 'L'${\mathbf{side}}=\text{'L'}$ or the matrix A U$AU$ when side = 'R'${\mathbf{side}}=\text{'R'}$.
3:     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:
ifail = 1${\mathbf{ifail}}=1$
 On entry, side ≠ 'L'${\mathbf{side}}\ne \text{'L'}$ or 'R'$\text{'R'}$.
ifail = 2${\mathbf{ifail}}=2$
 On entry, init ≠ 'I'${\mathbf{init}}\ne \text{'I'}$ or 'N'$\text{'N'}$.
ifail = 3${\mathbf{ifail}}=3$
 On entry, m < 1${\mathbf{m}}<1$.
ifail = 4${\mathbf{ifail}}=4$
 On entry, n < 1${\mathbf{n}}<1$.
ifail = 5${\mathbf{ifail}}=5$
 On entry, state vector was not initialized or has been corrupted.
ifail = 7${\mathbf{ifail}}=7$
 On entry, lda < m$\mathit{lda}<{\mathbf{m}}$.
ifail = 8${\mathbf{ifail}}=8$
On entry, an orthogonal matrix of dimension 1$1$ has been requested.

## Accuracy

The maximum error in UT U${U}^{\mathrm{T}}U$ should be a modest multiple of machine precision (see Chapter X02).

None.

## Example

```function nag_rand_matrix_orthog_example
% Initialize the seed
seed = [int64(1762543)];
% genid and subid identify the base generator
genid = int64(1);
subid =  int64(1);

side = 'Right';
init = 'Initialize';
a = zeros(4, 4);
% Initialize the generator to a repeatable sequence
[state, ifail] = nag_rand_init_repeat(genid, subid, seed);

[state, a, ifail] = nag_rand_matrix_orthog(side, init, state, a)
```
```

state =

17
1234
1
0
17510
6613
28263
24312
17917
13895
19930
8
0
1234
1
1
1234

a =

0.1756    0.7401   -0.3067   -0.5722
0.6593   -0.5781   -0.2191   -0.4279
0.6680    0.3172    0.6077    0.2895
-0.2971   -0.1323    0.6990   -0.6369

ifail =

0

```
```function g05px_example
% Initialize the seed
seed = [int64(1762543)];
% genid and subid identify the base generator
genid = int64(1);
subid =  int64(1);

side = 'Right';
init = 'Initialize';
a = zeros(4, 4);
% Initialize the generator to a repeatable sequence
[state, ifail] = g05kf(genid, subid, seed);

[state, a, ifail] = g05px(side, init, state, a)
```
```

state =

17
1234
1
0
17510
6613
28263
24312
17917
13895
19930
8
0
1234
1
1
1234

a =

0.1756    0.7401   -0.3067   -0.5722
0.6593   -0.5781   -0.2191   -0.4279
0.6680    0.3172    0.6077    0.2895
-0.2971   -0.1323    0.6990   -0.6369

ifail =

0

```