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_complex_norm_rcomm (f04zc)

## Purpose

nag_linsys_complex_norm_rcomm (f04zc) estimates the $1$-norm of a complex 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 f04zc in Advice on Replacement Calls for Withdrawn/Superseded Routines..

## Syntax

[icase, x, estnrm, work, ifail] = f04zc(icase, x, estnrm, work, 'n', n)
[icase, x, estnrm, work, ifail] = nag_linsys_withdraw_complex_norm_rcomm(icase, x, estnrm, work, 'n', n)

## Description

nag_linsys_complex_norm_rcomm (f04zc) computes an estimate (a lower bound) for the $1$-norm
 $A1 = max 1≤j≤n ∑ i=1 n aij$ (1)
of an $n$ by $n$ complex 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{H}}x$, where ${A}^{\mathrm{H}}$ is the complex conjugate transpose. 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 the 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{H}}x$ required by nag_linsys_complex_norm_rcomm (f04zc) 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 in Higham (1988).
Since ${‖A‖}_{\infty }={‖{A}^{\mathrm{H}}‖}_{1}$, nag_linsys_complex_norm_rcomm (f04zc) can be used to estimate the $\infty$-norm of $A$ by working with ${A}^{\mathrm{H}}$ instead of $A$.
The algorithm used is based on a method given in Hager (1984) and is described in Higham (1988). A comparison of several techniques for condition number estimation is given in Higham (1987).
Note: nag_linsys_complex_gen_norm_rcomm (f04zd) can also be used to estimate the norm of a real matrix. nag_linsys_complex_gen_norm_rcomm (f04zd) uses a more recent algorithm than nag_linsys_complex_norm_rcomm (f04zc) and it is recommended that nag_linsys_complex_gen_norm_rcomm (f04zd) be used in place of nag_linsys_complex_norm_rcomm (f04zc).

## 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)$ – complex array
On initial entry: need not be set.
On intermediate re-entry: must contain $Ax$ (if ${\mathbf{icase}}=1$) or ${A}^{\mathrm{H}}x$ (if ${\mathbf{icase}}=2$).
3:     $\mathrm{estnrm}$ – double scalar
On initial entry: need not be set.
4:     $\mathrm{work}\left({\mathbf{n}}\right)$ – complex array
On initial entry: need not be set.

### Optional Input Parameters

1:     $\mathrm{n}$int64int32nag_int scalar
Default: the dimension of the arrays x, work. (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{H}}x$ (if ${\mathbf{icase}}=2$), where ${A}^{\mathrm{H}}$ is the complex conjugate transpose; (b) place the result in x; and, (c) call nag_linsys_complex_norm_rcomm (f04zc) once again, with all the other arguments unchanged.
On final exit: ${\mathbf{icase}}=0$.
2:     $\mathrm{x}\left({\mathbf{n}}\right)$ – complex 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)$ – complex 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{ifail}$int64int32nag_int scalar
${\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 by nag_linsys_complex_norm_rcomm (f04zc) is proportional to $n$. For most problems the time taken during calls to nag_linsys_complex_norm_rcomm (f04zc) will be negligible compared with the time spent evaluating matrix-vector products between calls to nag_linsys_complex_norm_rcomm (f04zc).
The number of matrix-vector products required varies from $5$ to $11$ (or is $1$ if $n=1$). In most cases $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_complex_norm_rcomm (f04zc) 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, ifail] = f04zc(icase, x, estnrm, work);
while (icase ~= 0)
if (icase == 1)
...  code to compute A(-1)x ...
else
...  code to compute (A(-1)(H)) x ...
end
[icase, x, estnrm, work, ifail] = f04zc(icase, x, estnrm, work);
end
end```
To compute ${A}^{-1}x$ or ${A}^{-\mathrm{H}}x$, solve the equation $Ay=x$ or ${A}^{\mathrm{H}}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 Hermitian matrix, then $A={A}^{\mathrm{H}}$, and the if statement after the while can be reduced to:
```    ...  code to compute A(-1)x ...
```
The example program in Example illustrates how nag_linsys_complex_norm_rcomm (f04zc) can be used in conjunction with NAG Toolbox functions for complex band matrices (factorized by nag_lapack_zgbtrf (f07br)).
It is also straightforward to use nag_linsys_complex_norm_rcomm (f04zc) for Hermitian positive definite matrices, using nag_lapack_zpotrf (f07fr) and nag_lapack_zpotrs (f07fs) for factorization and solution.

## Example

This example estimates the condition number ${‖A‖}_{1}{‖{A}^{-1}‖}_{1}$ of the order $5$ matrix
 $A = 1+0i 2+0i 1+2i 0i+0 0i+0 2i 3+5i 1+3i 2+0i 0i+0 0i+0 -2+6i 5+7i 6i 1-0i 0i+0 0i+0 3+9i 4i 4-3i 0i+0 0i+0 0i+0 -1+8i 10-3i$
where $A$ is a band matrix stored in the packed format required by nag_lapack_zgbtrf (f07br) and nag_lapack_zgbtrs (f07bs).
Further examples of the technique for condition number estimation in the case of double matrices can be seen in the example program section of nag_linsys_real_norm_rcomm (f04yc).
```function f04zc_example

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

a = [ 1.0 + 1.0i,  2.0 + 1.0i,  1.0 + 2.0i,  0.0 + 0.0i,  0.0 + 0.0i;
0.0 + 2.0i,  3.0 + 5.0i,  1.0 + 3.0i,  2.0 + 1.0i,  0.0 + 0.0i,
0.0 + 0.0i, -2.0 + 6.0i,  5.0 + 7.0i,  0.0 + 6.0i,  1.0 - 1.0i;
0.0 + 0.0i,  0.0 + 0.0i,  3.0 + 9.0i,  0.0 + 4.0i,  4.0 - 3.0i;
0.0 + 0.0i,  0.0 + 0.0i,  0.0 + 0.0i, -1.0 + 8.0i, 10.0 - 3.0i];

x     = complex(zeros(5, 1));
work  = complex(zeros(5,1));
anorm = norm(a,1);
icase = int64(0);
estnrm = 0;

done = false;
while (~done)
[icase, x, estnrm, work, ifail] = ...
f04zc(icase, x, estnrm, work);
if (icase == 0)
done = true;
elseif (icase == 1)
x = inv(a)*x;
else
x = conj(transpose(inv(a)))*x;
end
end
fprintf('Computed norm of a              = %6.4g\n', anorm);
fprintf('Estimated norm of inverse(A)    = %6.4g\n', estnrm);
fprintf('Estimated condition number of A = %6.1f\n', estnrm*anorm);

```
```f04zc example results

Computed norm of a              =  23.49
Estimated norm of inverse(A)    =  37.04
Estimated condition number of A =  870.0
```