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

## Purpose

nag_linsys_real_norm_rcomm (f04yc) estimates the $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_withdraw_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$-norm
 $A1 = max 1≤j≤n ∑i=1n aij$ (1)
of an $n$ by $n$ real matrix $A=\left({a}_{ij}\right)$. The function regards the matrix $A$ as being defined by a user-supplied ‘Black Box’ which, given an input vector $x$, can return either of the matrix-vector products $Ax$ or ${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$ are known explicitly; it is then more efficient to compute the $1$-norm directly from formula (1) above.
The main use of the function is for estimating ${‖{B}^{-1}‖}_{1}$, and hence the condition number ${\kappa }_{1}\left(B\right)={‖B‖}_{1}{‖{B}^{-1}‖}_{1}$, without forming ${B}^{-1}$ explicitly ($A={B}^{-1}$ above).
If, for example, an $LU$ factorization of $B$ is available, the matrix-vector products ${B}^{-1}x$ and ${B}^{-\mathrm{T}}x$ required by nag_linsys_real_norm_rcomm (f04yc) may be computed by back- and forward-substitutions, without computing ${B}^{-1}$.
The function can also be used to estimate $1$-norms of matrix products such as ${A}^{-1}B$ and $ABC$, without forming the products explicitly. Further applications are described by Higham (1988).
Since ${‖A‖}_{\infty }={‖{A}^{\mathrm{T}}‖}_{1}$, nag_linsys_real_norm_rcomm (f04yc) can be used to estimate the $\infty$-norm of $A$ by working with ${A}^{\mathrm{T}}$ instead of $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 argument icase. Between intermediate exits and re-entries, all arguments other than x must remain unchanged.

### Compulsory Input Parameters

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

### Optional Input Parameters

1:     $\mathrm{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$, the order of the matrix $A$.
Constraint: ${\mathbf{n}}\ge 1$.

### Output Parameters

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

## Accuracy

In extensive tests on random matrices of size up to $n=100$ the estimate estnrm has been found always to be within a factor eleven of ${‖A‖}_{1}$; often the estimate has many correct figures. However, matrices exist for which the estimate is smaller than ${‖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$. 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$ to $11$ (or is $1$ if $n=1$). In most cases $4$ or $5$ products are required; it is rare for more than $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 ${‖{B}^{-1}‖}_{1}$ using a triangular factorization of $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$-norm of the inverse of a matrix $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 ${A}^{-1}x$ or ${A}^{-\mathrm{T}}x$, solve the equation $Ay=x$ or ${A}^{\mathrm{T}}y=x$ for $y$, overwriting $y$ on $x$. The code will vary, depending on the type of the matrix $A$, and the NAG function used to factorize $A$.
Note that if $A$ is any type of symmetric matrix, then $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 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

For this function two examples are presented. There is a single example program for nag_linsys_real_norm_rcomm (f04yc), with a main program and the code to solve the two example problems is given in Example 1 (EX1) and Example 2 (EX2).
Example 1 (EX1)
To estimate the condition number ${‖A‖}_{1}{‖{A}^{-1}‖}_{1}$ of the matrix $A$ given by
 $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 .$
Example 2 (EX2)
To estimate the condition number ${‖A‖}_{1}{‖{A}^{-1}‖}_{1}$ of the matrix $A$ given by
 $A= 5.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -1.0 2.0 0.0 0.0 0.0 2.0 3.0 0.0 0.0 0.0 -2.0 0.0 0.0 1.0 1.0 0.0 -1.0 0.0 0.0 -1.0 2.0 -3.0 -1.0 -1.0 0.0 0.0 0.0 6.0 .$
```function f04yc_example

fprintf('f04yc example results\n\n');

ex1;
ex2;

function ex1
fprintf('Example 1:\n\n');
% Estimate condition number of A
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];
anorm = norm(a,1);
ainv = inv(a);

% reverse communication initializations
icase = int64(0);
estnrm = 0;
x     = zeros(5, 1);
work  = zeros(5, 1);
iwork = zeros(5, 1, 'int64');

% reverse communication loop
done = false;
while (~done)
[icase, x, estnrm, work, iwork, ifail] = ...
f04yc(icase, x, estnrm, work, iwork);
if (icase == 0)
done = true;
elseif (icase == 1)
x = ainv*x;
elseif (icase == 2)
x = ainv'*x;
end
end

fprintf('Computed norm of A              = %8.4f\n', anorm);
fprintf('Estimated norm of inverse(A)    = %8.4f\n', estnrm);
cond = anorm*estnrm;
fprintf('Estimated condition number of A = %8.1f\n', cond);

function ex2
fprintf('\nExample 2:\n\n');
% Estimate condition number of sparse A
n  = int64(6);
nz = int64(15);
a   = zeros(150,1);
irn = zeros(75,1,'int64');
icn = zeros(150,1,'int64');
a(1:nz) = [ 5                       ...
2  -1   2           ...
3               ...
-2           1   1       ...
-1          -1   2  -3   ...
-1  -1               6];
anorm = 9;
irn(1:15) = [int64(1); 2;2;2; 3; 4;4;4; 5;5;5;5; 6;6;6];
icn(1:15) = [int64(1); 2;3;4; 3; 1;4;5; 1;4;5;6; 1;2;6];

% Factorize A
abort     = [true;     true;     false;     true];
[a, irn, icn, ikeep, w, idisp, ifail] = ...
f01br( ...
n, nz, a, irn, icn, abort);

% reverse communication initializations
icase = int64(0);
estnrm = 0;
x     = zeros(n, 1);
work  = zeros(n, 1);
iwork = zeros(n, 1, 'int64');

% reverse communication loop
done = false;
while (~done)
[icase, x, estnrm, work, iwork, ifail] = ...
f04yc(icase, x, estnrm, work, iwork);
if (icase == 0)
done = true;
else
% Solve Ax = b or A'x = b
[x, resid] = f04ax( ...
a, icn, ikeep, x, icase, idisp);
end
end

fprintf('Computed norm of A              = %8.4f\n', anorm);
fprintf('Estimated norm of inverse(A)    = %8.4f\n', estnrm);
cond = anorm*estnrm;
fprintf('Estimated condition number of A = %8.1f\n', cond);
```
```f04yc example results

Example 1:

Computed norm of A              =  15.9000
Estimated norm of inverse(A)    =   1.7635
Estimated condition number of A =     28.0

Example 2:

Computed norm of A              =   9.0000
Estimated norm of inverse(A)    =   1.9333
Estimated condition number of A =     17.4
```