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_complex_gen_norm_rcomm (f04zd)

## Purpose

nag_linsys_complex_gen_norm_rcomm (f04zd) estimates the $1$-norm of a complex 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, rwork, iwork, ifail] = f04zd(irevcm, x, y, estnrm, seed, work, rwork, iwork, 'm', m, 'n', n, 't', t)
[irevcm, x, y, estnrm, work, rwork, iwork, ifail] = nag_linsys_complex_gen_norm_rcomm(irevcm, x, y, estnrm, seed, work, rwork, iwork, 'm', m, 'n', n, 't', t)

## Description

nag_linsys_complex_gen_norm_rcomm (f04zd) computes an estimate (a lower bound) for the $1$-norm
 $A1 = max 1≤j≤n ∑ i=1 m aij$ (1)
of an $m$ 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 $n×t$ matrix $X$ (with $t\ll n$) or an $m×t$ matrix $Y$, can return $AX$ or ${A}^{\mathrm{H}}Y$, 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 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}$ for a square matrix $B$, 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 products ${B}^{-1}X$ and ${B}^{-\mathrm{H}}Y$ required by nag_linsys_complex_gen_norm_rcomm (f04zd) 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_gen_norm_rcomm (f04zd) can be used to estimate the $\infty$-norm of $A$ by working with ${A}^{\mathrm{H}}$ instead of $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$-norm estimation, with an application to $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 argument irevcm. Between intermediate exits and re-entries, all arguments other than x and y must remain unchanged.

### Compulsory Input Parameters

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

### Optional Input Parameters

1:     $\mathrm{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$.
Constraint: ${\mathbf{m}}\ge 0$.
2:     $\mathrm{n}$int64int32nag_int scalar
Default: the first dimension of the array x.
On initial entry: $n$, the number of columns of the matrix $A$.
Constraint: ${\mathbf{n}}\ge 0$.
3:     $\mathrm{t}$int64int32nag_int scalar
Suggested value: ${\mathbf{t}}=2$.
Default: the second dimension of the arrays x, y. (An error is raised if these dimensions are not equal.)
The number of columns $t$ of the matrices $X$ and $Y$. This is a argument that can be used to control the accuracy and reliability of the estimate and corresponds roughly to the number of columns of $A$ that are visited during each iteration of the algorithm.
If ${\mathbf{t}}\ge 2$ then a partly random starting matrix is used in the algorithm.
Constraint: $1\le {\mathbf{t}}\le {\mathbf{m}}$.

### Output Parameters

1:     $\mathrm{irevcm}$int64int32nag_int scalar
On intermediate exit: ${\mathbf{irevcm}}=1$ or $2$, and x contains the $n×t$ matrix $X$ and y contains the $m×t$ matrix $Y$. The calling program must
 (a) if ${\mathbf{irevcm}}=1$, evaluate $AX$ and store the result in y or if ${\mathbf{irevcm}}=2$, evaluate ${A}^{\mathrm{H}}Y$ and store the result in x, where ${A}^{\mathrm{H}}$ is the complex conjugate transpose; (b) call nag_linsys_complex_gen_norm_rcomm (f04zd) once again, with all the arguments unchanged.
On final exit: ${\mathbf{irevcm}}=0$.
2:     $\mathrm{x}\left(\mathit{ldx},:\right)$ – complex array
The first dimension of the array x will be ${\mathbf{n}}$.
The second dimension of the array x will be $\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{t}}\right)$.
On intermediate exit: if ${\mathbf{irevcm}}=1$, contains the current matrix $X$.
On final exit: the array is undefined.
3:     $\mathrm{y}\left(\mathit{ldy},:\right)$ – complex array
The first dimension of the array y will be ${\mathbf{m}}$.
The second dimension of the array y will be $\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{t}}\right)$.
On intermediate exit: if ${\mathbf{irevcm}}=2$, contains the current matrix $Y$.
On final exit: the array is undefined.
4:     $\mathrm{estnrm}$ – double scalar
On final exit: an estimate (a lower bound) for ${‖A‖}_{1}$.
5:     $\mathrm{work}\left({\mathbf{m}}×{\mathbf{t}}\right)$ – complex array
6:     $\mathrm{rwork}\left(2×{\mathbf{n}}\right)$ – double array
7:     $\mathrm{iwork}\left(2×{\mathbf{n}}+5×{\mathbf{t}}+20\right)$int64int32nag_int array
8:     $\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$
${\mathbf{ifail}}=-1$
Constraint: ${\mathbf{irevcm}}=0$, $1$ or $2$.
On initial entry, ${\mathbf{irevcm}}=_$.
Constraint: ${\mathbf{irevcm}}=0$.
${\mathbf{ifail}}=-2$
Constraint: ${\mathbf{m}}\ge 0$.
${\mathbf{ifail}}=-3$
Constraint: ${\mathbf{n}}\ge 0$.
${\mathbf{ifail}}=-5$
Constraint: $\mathit{ldx}\ge {\mathbf{n}}$.
${\mathbf{ifail}}=-7$
Constraint: $\mathit{ldy}\ge {\mathbf{m}}$.
${\mathbf{ifail}}=-9$
Constraint: $1\le {\mathbf{t}}\le {\mathbf{m}}$.
${\mathbf{ifail}}=-10$
Constraint: if ${\mathbf{t}}>1$, ${\mathbf{seed}}\ge 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 $m=n=450$ the estimate estnrm has been found always to be within a factor two 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 and Tisseur (2000) for further details.

### Timing

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

### Overflow

It is your responsibility to guard against potential overflows during evaluation of the matrix products. In particular, when estimating ${‖{B}^{-1}‖}_{1}$ using a triangular factorization of $B$, nag_linsys_complex_gen_norm_rcomm (f04zd) 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 argument $t$ controls the accuracy and reliability of the estimate. For $t=1$, the algorithm behaves similarly to the LAPACK estimator xLACON. Increasing $t$ typically improves the estimate, without increasing the number of iterations required.
For $t\ge 2$, random matrices are used in the algorithm, so for repeatable results the same value of seed should be used each time.
A value of $t=2$ is recommended for new users.

### 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] = f04zd(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] = f04zd(icase, x, estnrm, work);
end
end```
To compute ${A}^{-1}X$ or ${A}^{-\mathrm{H}}Y$, solve the equation $AY=X$ or ${A}^{\mathrm{H}}X=Y$ storing the result in y or x respectively. The code will vary, depending on the type of the matrix $A$, and the NAG function used to factorize $A$.
The example program in Example illustrates how nag_linsys_complex_gen_norm_rcomm (f04zd) can be used in conjunction with NAG Toolbox function for $LU$ factorization of complex matrices nag_lapack_zgetrf (f07ar)).
It is also straightforward to use nag_linsys_complex_gen_norm_rcomm (f04zd) 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 matrix $A$ given by
 $A = 0.7+0.1i -0.2+0.0i 1.0+0.0i 0.0+0.0i 0.0+0.0i 0.1+0.0i 0.3+0.0i 0.7+0.0i 0.0+0.0i 1.0+0.2i 0.9+0.0i 0.2+0.0i 0.0+5.9i 0.0+0.0i 0.2+0.0i 0.7+0.0i 0.4+6.1i 1.1+0.4i 0.0+0.1i 0.0+0.1i -0.7+0.0i 0.2+0.0i 0.1+0.0i 0.1+0.0i 0.0+0.0i 4.0+0.0i 0.0+0.0i 1.0+0.0i 9.0+0.0i 0.0+0.1i 4.5+6.7i 0.1+0.4i 0.0+3.2i 1.2+0.0i 0.0+0.0i 7.8+0.2i .$
```function f04zd_example

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

a = [0.7+0.1i, -0.2+0.0i,  1.0+0.0i, 0.0+0.0i, 0.0+0.0i, 0.1+0.0i;
0.3+0.0i,  0.7+0.0i,  0.0+0.0i, 1.0+0.2i, 0.9+0.0i, 0.2+0.0i;
0.0+5.9i,  0.0+0.0i,  0.2+0.0i, 0.7+0.0i, 0.4+6.1i, 1.1+0.4i;
0.0+0.1i,  0.0+0.1i, -0.7+0.0i, 0.2+0.0i, 0.1+0.0i, 0.1+0.0i;
0.0+0.0i,  4.0+0.0i,  0.0+0.0i, 1.0+0.0i, 9.0+0.0i, 0.0+0.1i;
4.5+6.7i,  0.1+0.4i,  0.0+3.2i, 1.2+0.0i, 0.0+0.0i, 7.8+0.2i];

t = int64(2);
m = 6;
n = 6;
x = complex(zeros(n, t));
y = complex(zeros(m, t));
estnrm = 0;
seed = int64(652);
work  = complex(zeros(m*t, 1));
rwork = zeros(2*n, 1);
iwork = zeros(2*n+5*t+20, 1, 'int64');

nrma =  norm(a, 1);

% 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] = f07ar(a);

done = false;
irevcm = int64(0);

while ~done

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

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

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

```
```f04zd example results

The norm of a                    =  16.11
The estimated norm of inverse(a) =  24.02
Estimated condition number of a  = 387.08
```