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_linsys_real_gen_norm_rcomm (f04yd)

Purpose

nag_linsys_real_gen_norm_rcomm (f04yd) estimates the 1$1$-norm of a real rectangular matrix without accessing the matrix explicitly. It uses reverse communication for evaluating matrix products. The function may be used for estimating condition numbers of square matrices.

Syntax

[irevcm, x, y, estnrm, work, iwork, ifail] = f04yd(irevcm, x, y, estnrm, seed, work, iwork, 'm', m, 'n', n, 't', t)
[irevcm, x, y, estnrm, work, iwork, ifail] = nag_linsys_real_gen_norm_rcomm(irevcm, x, y, estnrm, seed, work, iwork, 'm', m, 'n', n, 't', t)

Description

nag_linsys_real_gen_norm_rcomm (f04yd) computes an estimate (a lower bound) for the 1$1$-norm
 m ‖A‖1 = max ∑ |aij| 1 ≤ j ≤ n i = 1
$‖A‖1 = max 1≤j≤n ∑ i=1 m |aij|$
(1)
of an m$m$ by n$n$ real matrix A = (aij)$A=\left({a}_{ij}\right)$. The function regards the matrix A$A$ as being defined by a user-supplied ‘Black Box’ which, given an n × t$n×t$ matrix X$X$ (with tn$t\ll n$) or an m × t$m×t$ matrix Y$Y$, can return AX$AX$ or ATY${A}^{\mathrm{T}}Y$. A reverse communication interface is used; thus control is returned to the calling program whenever a matrix product is required.
Note:  this function is not recommended for use when the elements of A$A$ are known explicitly; it is then more efficient to compute the 1$1$-norm directly from formula (1) above.
The main use of the function is for estimating B11${‖{B}^{-1}‖}_{1}$ for a square matrix, B$B$, and hence the condition number κ1(B) = B1B11${\kappa }_{1}\left(B\right)={‖B‖}_{1}{‖{B}^{-1}‖}_{1}$, without forming B1${B}^{-1}$ explicitly (A = B1$A={B}^{-1}$ above).
If, for example, an LU$LU$ factorization of B$B$ is available, the matrix products B1X${B}^{-1}X$ and BTY${B}^{-\mathrm{T}}Y$ required by nag_linsys_real_gen_norm_rcomm (f04yd) may be computed by back- and forward-substitutions, without computing B1${B}^{-1}$.
The function can also be used to estimate 1$1$-norms of matrix products such as A1B${A}^{-1}B$ and ABC$ABC$, without forming the products explicitly. Further applications are described by Higham (1988).
Since A = AT1${‖A‖}_{\infty }={‖{A}^{\mathrm{T}}‖}_{1}$, nag_linsys_real_gen_norm_rcomm (f04yd) can be used to estimate the $\infty$-norm of A$A$ by working with AT${A}^{\mathrm{T}}$ instead of A$A$.
The algorithm used is described in Higham and Tisseur (2000).

References

Higham N J (1988) FORTRAN codes for estimating the one-norm of a real or complex matrix, with applications to condition estimation ACM Trans. Math. Software 14 381–396
Higham N J and Tisseur F (2000) A block algorithm for matrix 1$1$-norm estimation, with an application to 1$1$-norm pseudospectra SIAM J. Matrix. Anal. Appl. 21 1185–1201

Parameters

Note: this function uses reverse communication. Its use involves an initial entry, intermediate exits and re-entries, and a final exit, as indicated by the parameter irevcm. Between intermediate exits and re-entries, all parameters other than x and y must remain unchanged.

Compulsory Input Parameters

1:     irevcm – int64int32nag_int scalar
On initial entry: must be set to 0$0$.
On intermediate re-entry: irevcm must be unchanged.
2:     x(ldx, : $:$) – double array
The first dimension of the array x must be at least n${\mathbf{n}}$
The second dimension of the array must be at least max (1,t)$\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{t}}\right)$
On initial entry: need not be set.
On intermediate re-entry: if irevcm = 2${\mathbf{irevcm}}=2$, must contain ATY${A}^{\mathrm{T}}Y$.
3:     y(ldy, : $:$) – double array
The first dimension of the array y must be at least m${\mathbf{m}}$
The second dimension of the array must be at least max (1,t)$\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{t}}\right)$
On initial entry: need not be set.
On intermediate re-entry: if irevcm = 1${\mathbf{irevcm}}=1$, must contain AX$AX$.
4:     estnrm – double scalar
On initial entry: need not be set.
On intermediate re-entry: must not be changed.
5:     seed – int64int32nag_int scalar
The seed used for random number generation.
If t = 1${\mathbf{t}}=1$, seed is not used.
6:     work(m × t${\mathbf{m}}×{\mathbf{t}}$) – double array
7:     iwork(2 × n + 5 × t + 20$2×{\mathbf{n}}+5×{\mathbf{t}}+20$) – int64int32nag_int array
On initial entry: need not be set.
On intermediate re-entry: must not be changed.

Optional Input Parameters

1:     m – int64int32nag_int scalar
Default: The first dimension of the arrays y, work. (An error is raised if these dimensions are not equal.)
The number of rows of the matrix A$A$.
Constraint: m0${\mathbf{m}}\ge 0$.
2:     n – int64int32nag_int scalar
Default: The first dimension of the array x and the dimension of the array iwork. (An error is raised if these dimensions are not equal.)
n$n$, the number of columns of the matrix A$A$.
Constraint: n0${\mathbf{n}}\ge 0$.
3:     t – int64int32nag_int scalar
Default: The second dimension of the arrays x, y.
The number of columns t$t$ of the matrices X$X$ and Y$Y$. This is a parameter that can be used to control the accuracy and reliability of the estimate and corresponds roughly to the number of columns of A$A$ that are visited during each iteration of the algorithm.
If t2${\mathbf{t}}\ge 2$ then a partly random starting matrix is used in the algorithm.
Default: t = 2${\mathbf{t}}=2$
Constraint: 1tm$1\le {\mathbf{t}}\le {\mathbf{m}}$.

ldx ldy

Output Parameters

1:     irevcm – int64int32nag_int scalar
On intermediate exit: irevcm = 1${\mathbf{irevcm}}=1$ or 2$2$, and x and y contain the elements of n × t$n×t$ and m × t$m×t$ matrices respectively. The calling program must
 (a) if irevcm = 1${\mathbf{irevcm}}=1$, evaluate AX$AX$ and store the result in y or if irevcm = 2${\mathbf{irevcm}}=2$, evaluate ATY${A}^{\mathrm{T}}Y$ and store the result in x, (b) call nag_linsys_real_gen_norm_rcomm (f04yd) once again, with all the other parameters unchanged.
On final exit: irevcm = 0${\mathbf{irevcm}}=0$.
2:     x(ldx, : $:$) – double array
The first dimension of the array x will be n${\mathbf{n}}$
The second dimension of the array will be max (1,t)$\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{t}}\right)$
ldxn$\mathit{ldx}\ge {\mathbf{n}}$.
On intermediate exit: if irevcm = 1${\mathbf{irevcm}}=1$, contains the current matrix X$X$.
ldxn$\mathit{ldx}\ge {\mathbf{n}}$.
On final exit: the array is undefined.
3:     y(ldy, : $:$) – double array
The first dimension of the array y will be m${\mathbf{m}}$
The second dimension of the array will be max (1,t)$\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{t}}\right)$
ldym$\mathit{ldy}\ge {\mathbf{m}}$.
On intermediate exit: if irevcm = 2${\mathbf{irevcm}}=2$, contains the current matrix Y$Y$.
ldym$\mathit{ldy}\ge {\mathbf{m}}$.
On final exit: the array is undefined.
4:     estnrm – double scalar
On final exit: an estimate (a lower bound) for A1${‖A‖}_{1}$.
5:     work(m × t${\mathbf{m}}×{\mathbf{t}}$) – double array
6:     iwork(2 × n + 5 × t + 20$2×{\mathbf{n}}+5×{\mathbf{t}}+20$) – int64int32nag_int array
7:     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$
ifail = 1${\mathbf{ifail}}=-1$
Constraint: irevcm = 0${\mathbf{irevcm}}=0$, 1$1$ or 2$2$.
On initial entry.
Constraint: irevcm = 0${\mathbf{irevcm}}=0$.
ifail = 2${\mathbf{ifail}}=-2$
Constraint: m0${\mathbf{m}}\ge 0$.
ifail = 3${\mathbf{ifail}}=-3$
Constraint: n0${\mathbf{n}}\ge 0$.
ifail = 5${\mathbf{ifail}}=-5$
Constraint: ldxn$\mathit{ldx}\ge {\mathbf{n}}$.
ifail = 7${\mathbf{ifail}}=-7$
Constraint: ldym$\mathit{ldy}\ge {\mathbf{m}}$.
ifail = 9${\mathbf{ifail}}=-9$
Constraint: 1tm$1\le {\mathbf{t}}\le {\mathbf{m}}$.

Accuracy

In extensive tests on random matrices of size up to m = n = 450$m=n=450$ the estimate estnrm has been found always to be within a factor two of A1${‖A‖}_{1}$; often the estimate has many correct figures. However, matrices exist for which the estimate is smaller than A1${‖A‖}_{1}$ by an arbitrary factor; such matrices are very unlikely to arise in practice. See Higham and Tisseur (2000) for further details.

Timing

For most problems the time taken during calls to nag_linsys_real_gen_norm_rcomm (f04yd) will be negligible compared with the time spent evaluating matrix products between calls to nag_linsys_real_gen_norm_rcomm (f04yd).
The number of matrix products required depends on the matrix A$A$. At most six products of the form Y = AX$Y=AX$ and five products of the form X = ATY$X={A}^{\mathrm{T}}Y$ will be required. The number of iterations is independent of the choice of t$t$.

Overflow

It is your responsibility to guard against potential overflows during evaluation of the matrix products. In particular, when estimating B11${‖{B}^{-1}‖}_{1}$ using a triangular factorization of B$B$, nag_linsys_real_gen_norm_rcomm (f04yd) should not be called if one of the factors is exactly singular – otherwise division by zero may occur in the substitutions.

Choice of t

The parameter t$t$ controls the accuracy and reliability of the estimate. For t = 1$t=1$, the algorithm behaves similarly to the LAPACK estimator xLACON. Increasing t$t$ typically improves the estimate, without increasing the number of iterations required.
For t2$t\ge 2$, random matrices are used in the algorithm, so for repeatable results the value of seed should be kept constant.
A value of t = 2$t=2$ is recommended for new users.

Use in Conjunction with NAG Library Routines

To estimate the 1$1$-norm of the inverse of a matrix A$A$, the following skeleton code can normally be used:
```...  code to factorize A ...
if (A is not singular)
icase = 0;
[icase, x, estnrm, work, iwork, ifail] = f04yd(icase, x, estnrm, work, iwork);
while (icase ~= 0)
if (icase == 1)
...  code to compute inv(A)*x ...
else
...  code to compute inv(transpose(A))*x ...
end
[icase, x, estnrm, work, iwork, ifail] = f04yd(icase, x, estnrm, work, iwork);
end
end
```
To compute A1X${A}^{-1}X$ or ATY${A}^{-\mathrm{T}}Y$, solve the equation AY = X$AY=X$ or ATX = Y${A}^{\mathrm{T}}X=Y$, storing the result in y or x respectively. The code will vary, depending on the type of the matrix A$A$, and the NAG function used to factorize A$A$.
The factorization will normally have been performed by a suitable function from Chapters F01, F03 or F07. Note also that many of the ‘Black Box’ functions in Chapter F04 for solving systems of equations also return a factorization of the matrix. The example program in Section [Example] illustrates how nag_linsys_real_gen_norm_rcomm (f04yd) can be used in conjunction with NAG Toolbox functions for LU$LU$ factorization of a real matrix nag_lapack_dgetrf (f07ad).
It is straightforward to use nag_linsys_real_gen_norm_rcomm (f04yd) for the following other types of matrix, using the named functions for factorization and solution:

Example

```function nag_linsys_real_gen_norm_rcomm_example

% Example 1: Compute the condition number of a dense matrix

a = [0.7,  -0.2,   1.0,   0.0,   2.0,   0.1;
0.3,   0.7,   0.0,   1.0,   0.9,   0.2;
0.0,   0.0,   0.2,   0.7,   0.0,  -1.1;
0.0,   3.4,  -0.7,   0.2,   0.1,   0.1;
0.0,  -4.0,   0.0,   1.0,   9.0,   0.0;
0.4,   1.2,   4.3,   0.0,   6.2,   5.9];
t = int64(2);
m = 6;
n = 6;
x = zeros(n, t);
y = zeros(m, t);
estnrm = 0;
seed = int64(354);
irevcm = int64(0);
work = zeros(n*t, 1);
iwork = zeros(2*n+5*t+20, 1, 'int64');

nrma =  norm(a, 1);
fprintf('\nExample 1\n');
fprintf('\nThe norm of a is %6.2f\n', nrma);

% Estimate the norm of a^(-1) without explicitly forming a^(-1)

% Perform an LU factorization so that A=LU where L and U are lower and upper
% triangular.
[a, ipiv, info] = nag_lapack_dgetrf(a);

first = true;

while first || (irevcm ~= 0)
first = false;

[irevcm, x, y, estnrm, work, iwork, ifail] = ...
nag_linsys_real_gen_norm_rcomm(irevcm, x, y, estnrm, seed, work, iwork);

switch irevcm
case 1
% Compute y = inv(a)*x
[y, info] = nag_lapack_dgetrs('n', a, ipiv, x);
case 2
% Compute x = transpose(inv(a))*y
[x, info] = nag_lapack_dgetrs('t', a, ipiv, y);
otherwise
end
end

fprintf('The estimated norm of inverse(a) is: %6.2f\n', estnrm);
fprintf('\nEstimated condition number of a: %6.2f\n', estnrm*nrma);

% Example 2: Compute the condition number of a sparse matrix
% (stored in symmetric coordinate storage format)

t = int64(2);
n = int64(5);
nz = int64(10);
a = zeros(4*nz, 1);
icn = zeros(4*nz, 1, 'int64');
irn = zeros(4*nz, 1, 'int64');
a(1:nz) = [3; 2; 1; 2; 1; 4; 2; 1; 2; 5];
irn(1:nz) = [2; 4; 2; 3; 5; 4; 5; 1; 3; 4];
icn(1:nz) = [1; 1; 2; 2; 2; 3; 3; 4; 4; 5];

x = zeros(n, t);
y = zeros(n, t);
estnrm = 0;
seed = int64(412);
irevcm = int64(0);
work = zeros(n*t, 1);
iwork = zeros(2*n+5*t+20, 1, 'int64');

% Compute 1-norm of a
nrma =  0;
for i = 1:n
asum = 0;
for j = 1:nz
if (icn(j)==i)
asum = asum + abs(a(j));
end
end
nrma = max(nrma,asum);
end
fprintf('\nExample 2\n');
fprintf('\nThe norm of a is %6.2f\n', nrma);

% Perform an LU factorization so that A=LU where L and U are lower and upper
% triangular.
abort = [true; true; false; false];
[a, irn, icn, ikeep, w, idisp, ifail] = nag_matop_real_gen_sparse_lu(n, nz, a, irn, icn, abort);

% Compute an estimate of the 1-norm of inv(a)

first = true;

while first || (irevcm ~= 0)
first = false;

[irevcm, x, y, estnrm, work, iwork, ifail] = ...
nag_linsys_real_gen_norm_rcomm(irevcm, x, y, estnrm, seed, work, iwork);

switch irevcm
case 1
% Compute y = inv(a)*x
for i=1:2
[y(:, i), resid] = nag_linsys_real_sparse_fac_solve(a, icn, ikeep, x(:, i), irevcm, idisp);
end
case 2
% Compute x = transpose(inv(a))*y
for i=1:2
[x(:, i), resid] = nag_linsys_real_sparse_fac_solve(a, icn, ikeep, y(:, i), irevcm, idisp);
end
otherwise
end
end

fprintf('The estimated norm of inverse(a) is: %6.2f\n', estnrm);
fprintf('\nEstimated condition number of a: %6.2f\n', estnrm*nrma);
```
```

Example 1

The norm of a is  18.20
The estimated norm of inverse(a) is:   2.97

Estimated condition number of a:  54.14

Example 2

The norm of a is   6.00
The estimated norm of inverse(a) is:   3.37

Estimated condition number of a:  20.20

```
```function f04yd_example

% Example 1: Compute the condition number of a dense matrix

a = [0.7,  -0.2,   1.0,   0.0,   2.0,   0.1;
0.3,   0.7,   0.0,   1.0,   0.9,   0.2;
0.0,   0.0,   0.2,   0.7,   0.0,  -1.1;
0.0,   3.4,  -0.7,   0.2,   0.1,   0.1;
0.0,  -4.0,   0.0,   1.0,   9.0,   0.0;
0.4,   1.2,   4.3,   0.0,   6.2,   5.9];
t = int64(2);
m = 6;
n = 6;
x = zeros(n, t);
y = zeros(m, t);
estnrm = 0;
seed = int64(354);
irevcm = int64(0);
work = zeros(n*t, 1);
iwork = zeros(2*n+5*t+20, 1, 'int64');

nrma =  norm(a, 1);
fprintf('\nExample 1\n');
fprintf('\nThe norm of a is %6.2f\n', nrma);

% Estimate the norm of a^(-1) without explicitly forming a^(-1)

% Perform an LU factorization so that A=LU where L and U are lower and upper
% triangular.

first = true;

while first || (irevcm ~= 0)
first = false;

[irevcm, x, y, estnrm, work, iwork, ifail] = ...
f04yd(irevcm, x, y, estnrm, seed, work, iwork);

switch irevcm
case 1
% Compute y = inv(a)*x
[y, info] = f07ae('n', a, ipiv, x);
case 2
% Compute x = transpose(inv(a))*y
[x, info] = f07ae('t', a, ipiv, y);
otherwise
end
end

fprintf('The estimated norm of inverse(a) is: %6.2f\n', estnrm);
fprintf('\nEstimated condition number of a: %6.2f\n', estnrm*nrma);

% Example 2: Compute the condition number of a sparse matrix
% (stored in symmetric coordinate storage format)

t = int64(2);
n = int64(5);
nz = int64(10);
a = zeros(4*nz, 1);
icn = zeros(4*nz, 1, 'int64');
irn = zeros(4*nz, 1, 'int64');
a(1:nz) = [3; 2; 1; 2; 1; 4; 2; 1; 2; 5];
irn(1:nz) = [2; 4; 2; 3; 5; 4; 5; 1; 3; 4];
icn(1:nz) = [1; 1; 2; 2; 2; 3; 3; 4; 4; 5];

x = zeros(n, t);
y = zeros(n, t);
estnrm = 0;
seed = int64(412);
irevcm = int64(0);
work = zeros(n*t, 1);
iwork = zeros(2*n+5*t+20, 1, 'int64');

% Compute 1-norm of a
nrma =  0;
for i = 1:n
asum = 0;
for j = 1:nz
if (icn(j)==i)
asum = asum + abs(a(j));
end
end
nrma = max(nrma,asum);
end
fprintf('\nExample 2\n');
fprintf('\nThe norm of a is %6.2f\n', nrma);

% Perform an LU factorization so that A=LU where L and U are lower and upper
% triangular.
abort = [true; true; false; false];
[a, irn, icn, ikeep, w, idisp, ifail] = f01br(n, nz, a, irn, icn, abort);

% Compute an estimate of the 1-norm of inv(a)

first = true;

while first || (irevcm ~= 0)
first = false;

[irevcm, x, y, estnrm, work, iwork, ifail] = ...
f04yd(irevcm, x, y, estnrm, seed, work, iwork);

switch irevcm
case 1
% Compute y = inv(a)*x
for i=1:2
[y(:, i), resid] = f04ax(a, icn, ikeep, x(:, i), irevcm, idisp);
end
case 2
% Compute x = transpose(inv(a))*y
for i=1:2
[x(:, i), resid] = f04ax(a, icn, ikeep, y(:, i), irevcm, idisp);
end
otherwise
end
end

fprintf('The estimated norm of inverse(a) is: %6.2f\n', estnrm);
fprintf('\nEstimated condition number of a: %6.2f\n', estnrm*nrma);
```
```

Example 1

The norm of a is  18.20
The estimated norm of inverse(a) is:   2.97

Estimated condition number of a:  54.14

Example 2

The norm of a is   6.00
The estimated norm of inverse(a) is:   3.37

Estimated condition number of a:  20.20

```