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_norm_rcomm (f04yc)

Purpose

nag_linsys_real_norm_rcomm (f04yc) estimates the 1$1$-norm of a real matrix without accessing the matrix explicitly. It uses reverse communication for evaluating matrix-vector products. The function may be used for estimating matrix condition numbers.
Note: this function is scheduled to be withdrawn, please see f04yc in Advice on Replacement Calls for Withdrawn/Superseded Routines..

Syntax

[icase, x, estnrm, work, iwork, ifail] = f04yc(icase, x, estnrm, work, iwork, 'n', n)
[icase, x, estnrm, work, iwork, ifail] = nag_linsys_real_norm_rcomm(icase, x, estnrm, work, iwork, 'n', n)

Description

nag_linsys_real_norm_rcomm (f04yc) computes an estimate (a lower bound) for the 1$1$-norm
 n ‖A‖1 = max ∑ |aij| 1 ≤ j ≤ n i = 1
$‖A‖1 = max 1≤j≤n ∑i=1n |aij|$
(1)
of an n$n$ 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 input vector x$x$, can return either of the matrix-vector products Ax$Ax$ or ATx${A}^{\mathrm{T}}x$. A reverse communication interface is used; thus control is returned to the calling program whenever a matrix-vector 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}$, 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-vector products B1x${B}^{-1}x$ and BTx${B}^{-\mathrm{T}}x$ required by nag_linsys_real_norm_rcomm (f04yc) 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_norm_rcomm (f04yc) 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 based on a method given by Hager (1984) and is described by Higham (1988). A comparison of several techniques for condition number estimation is given by Higham (1987).
Note: nag_linsys_real_gen_norm_rcomm (f04yd) can also be used to estimate the norm of a real matrix. nag_linsys_real_gen_norm_rcomm (f04yd) uses a more recent algorithm than nag_linsys_real_norm_rcomm (f04yc) and it is recommended that nag_linsys_real_gen_norm_rcomm (f04yd) be used in place of nag_linsys_real_norm_rcomm (f04yc).

References

Hager W W (1984) Condition estimates SIAM J. Sci. Statist. Comput. 5 311–316
Higham N J (1987) A survey of condition number estimation for triangular matrices SIAM Rev. 29 575–596
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

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 icase. Between intermediate exits and re-entries, all parameters other than x must remain unchanged.

Compulsory Input Parameters

1:     icase – int64int32nag_int scalar
On initial entry: must be set to 0$0$.
2:     x(n) – double array
n, the dimension of the array, must satisfy the constraint n1${\mathbf{n}}\ge 1$.
On initial entry: need not be set.
n, the dimension of the array, must satisfy the constraint n1${\mathbf{n}}\ge 1$.
On intermediate re-entry: must contain Ax$Ax$ (if icase = 1${\mathbf{icase}}=1$) or ATx${A}^{\mathrm{T}}x$ (if icase = 2${\mathbf{icase}}=2$).
3:     estnrm – double scalar
On initial entry: need not be set.
4:     work(n) – double array
n, the dimension of the array, must satisfy the constraint n1${\mathbf{n}}\ge 1$.
On initial entry: need not be set.
5:     iwork(n) – int64int32nag_int array

Optional Input Parameters

1:     n – int64int32nag_int scalar
Default: The dimension of the arrays x, work, iwork. (An error is raised if these dimensions are not equal.)
On initial entry: n$n$, the order of the matrix A$A$.
Constraint: n1${\mathbf{n}}\ge 1$.

None.

Output Parameters

1:     icase – int64int32nag_int scalar
On intermediate exit: icase = 1${\mathbf{icase}}=1$ or 2$2$, and x(i)${\mathbf{x}}\left(\mathit{i}\right)$, for i = 1,2,,n$\mathit{i}=1,2,\dots ,n$, contain the elements of a vector x$x$. The calling program must
 (a) evaluate Ax$Ax$ (if icase = 1${\mathbf{icase}}=1$) or ATx${A}^{\mathrm{T}}x$ (if icase = 2${\mathbf{icase}}=2$), (b) place the result in x, and (c) call nag_linsys_real_norm_rcomm (f04yc) once again, with all the other parameters unchanged.
On final exit: icase = 0${\mathbf{icase}}=0$.
2:     x(n) – double array
On intermediate exit: contains the current vector x$x$.
On final exit: the array is undefined.
3:     estnrm – double scalar
On intermediate exit: should not be changed.
On final exit: an estimate (a lower bound) for A1${‖A‖}_{1}$.
4:     work(n) – double array
On final exit: contains a vector v$v$ such that v = Aw$v=Aw$ where estnrm = v1 / w1${\mathbf{estnrm}}={‖v‖}_{1}/{‖w‖}_{1}$ (w$w$ is not returned). If A = B1$A={B}^{-1}$ and estnrm is large, then v$v$ is an approximate null vector for B$B$.
5:     iwork(n) – int64int32nag_int array
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:
ifail = 1${\mathbf{ifail}}=1$
On entry, n < 1${\mathbf{n}}<1$.

Accuracy

In extensive tests on random matrices of size up to n = 100$n=100$ the estimate estnrm has been found always to be within a factor eleven 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 (1988) for further details.

Timing

The total time taken within nag_linsys_real_norm_rcomm (f04yc) is proportional to n$n$. For most problems the time taken during calls to nag_linsys_real_norm_rcomm (f04yc) will be negligible compared with the time spent evaluating matrix-vector products between calls to nag_linsys_real_norm_rcomm (f04yc).
The number of matrix-vector products required varies from 4$4$ to 11$11$ (or is 1$1$ if n = 1$n=1$). In most cases 4$4$ or 5$5$ products are required; it is rare for more than 7$7$ to be needed.

Overflow

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

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] = f04yc(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] = f04yc(icase, x, estnrm, work, iwork);
end
end
```
To compute A1x${A}^{-1}x$ or ATx${A}^{-\mathrm{T}}x$, solve the equation Ay = x$Ay=x$ or ATy = x${A}^{\mathrm{T}}y=x$ for y$y$, overwriting y$y$ on x$x$. The code will vary, depending on the type of the matrix A$A$, and the NAG function used to factorize A$A$.
Note that if A$A$ is any type of symmetric matrix, then A = AT$A={A}^{\mathrm{T}}$, and the ifstatement after the while can be reduced to:
```       ...  code to compute inv(A)*x ...
```
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_norm_rcomm (f04yc) can be used in conjunction with NAG Toolbox functions for two important types of matrix: full nonsymmetric matrices (factorized by nag_lapack_dgetrf (f07ad)) and sparse nonsymmetric matrices (factorized by nag_matop_real_gen_sparse_lu (f01br)).
It is straightforward to use nag_linsys_real_norm_rcomm (f04yc) for the following other types of matrix, using the named functions for factorization and solution:

Example

```function nag_linsys_real_norm_rcomm_example
a = [ 1.5,  2.0,  3.0, -2.1,  0.3;
2.5,  3.0, -4.0,  2.3, -1.1;
3.5,  4.0,  0.5, -3.1, -1.4;
-0.4, -3.2, -2.1,  3.1,  2.1;
1.7,  3.7,  1.9, -2.2, -3.3];
icase = int64(0);
x = zeros(5, 1);
estnrm = 0;
work = zeros(5, 1);
iwork = zeros(5, 1, 'int64');
[icase, x, estnrm, work, iwork, ifail] = ...
nag_linsys_real_norm_rcomm(icase, x, estnrm, work, iwork);
while (icase > 0)
if (icase ==1)
x = a*x;
elseif (icase == 2)
x = transpose(a)*x;
end
[icase, x, estnrm, work, iwork, ifail] = ...
nag_linsys_real_norm_rcomm(icase, x, estnrm, work, iwork);
end
estnrm
```
```

estnrm =

15.9000

```
```function f04yc_example
a = [ 1.5,  2.0,  3.0, -2.1,  0.3;
2.5,  3.0, -4.0,  2.3, -1.1;
3.5,  4.0,  0.5, -3.1, -1.4;
-0.4, -3.2, -2.1,  3.1,  2.1;
1.7,  3.7,  1.9, -2.2, -3.3];
icase = int64(0);
x = zeros(5, 1);
estnrm = 0;
work = zeros(5, 1);
iwork = zeros(5, 1, 'int64');
[icase, x, estnrm, work, iwork, ifail] = f04yc(icase, x, estnrm, work, iwork);
while (icase > 0)
if (icase ==1)
x = a*x;
elseif (icase == 2)
x = transpose(a)*x;
end
[icase, x, estnrm, work, iwork, ifail] = f04yc(icase, x, estnrm, work, iwork);
end
estnrm
```
```

estnrm =

15.9000

```