NAG CL Interface
f04zdc (complex_gen_norm_rcomm)
1
Purpose
f04zdc 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.
2
Specification
void 
f04zdc (Integer *irevcm,
Integer m,
Integer n,
Complex x[],
Integer pdx,
Complex y[],
Integer pdy,
double *estnrm,
Integer t,
Integer seed,
Complex work[],
double rwork[],
Integer iwork[],
NagError *fail) 

The function may be called by the names: f04zdc or nag_linsys_complex_gen_norm_rcomm.
3
Description
f04zdc computes an estimate (a lower bound) for the
$1$norm
of an
$m$ by
$n$ complex matrix
$A=\left({a}_{ij}\right)$. The function regards the matrix
$A$ as being defined by a usersupplied ‘Black Box’ which, given an
$n\times t$ matrix
$X$ (with
$t\ll n$) or an
$m\times 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 ${\Vert {B}^{1}\Vert}_{1}$ for a square matrix $B$, and hence the condition number
${\kappa}_{1}\left(B\right)={\Vert B\Vert}_{1}{\Vert {B}^{1}\Vert}_{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 f04zdc may be computed by back and forwardsubstitutions, 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 ${\Vert A\Vert}_{\infty}={\Vert {A}^{\mathrm{H}}\Vert}_{1}$, f04zdc 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).
4
References
Higham N J (1988) FORTRAN codes for estimating the onenorm 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
5
Arguments
Note: this function uses
reverse communication. Its use involves an initial entry, intermediate exits and reentries, and a final exit, as indicated by the argument
irevcm. Between intermediate exits and reentries,
all arguments other than x and y must remain unchanged.

1:
$\mathbf{irevcm}$ – Integer *
Input/Output

On initial entry: must be set to $0$.
On intermediate exit:
${\mathbf{irevcm}}=1$ or
$2$, and
x contains the
$n\times t$ matrix
$X$ and
y contains the
$m\times 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 f04zdc once again, with all the arguments unchanged.
On intermediate reentry:
irevcm must be unchanged.
On final exit: ${\mathbf{irevcm}}=0$.
Note: any values you return to f04zdc as part of the reverse communication procedure should not include floatingpoint NaN (Not a Number) or infinity values, since these are not handled by f04zdc. If your code inadvertently does return any NaNs or infinities, f04zdc is likely to produce unexpected results.

2:
$\mathbf{m}$ – Integer
Input

On entry: the number of rows of the matrix $A$.
Constraint:
${\mathbf{m}}\ge 0$.

3:
$\mathbf{n}$ – Integer
Input

On initial entry: $n$, the number of columns of the matrix $A$.
Constraint:
${\mathbf{n}}\ge 0$.

4:
$\mathbf{x}\left[\mathit{dim}\right]$ – Complex
Input/Output

Note: the dimension,
dim, of the array
x
must be at least
${\mathbf{pdx}}\times {\mathbf{t}}$.
The $\left(i,j\right)$th element of the matrix $X$ is stored in ${\mathbf{x}}\left[\left(j1\right)\times {\mathbf{pdx}}+i1\right]$.
On initial entry: need not be set.
On intermediate exit:
if ${\mathbf{irevcm}}=1$, contains the current matrix $X$.
On intermediate reentry: if ${\mathbf{irevcm}}=2$, must contain ${A}^{\mathrm{H}}Y$.
On final exit: the array is undefined.

5:
$\mathbf{pdx}$ – Integer
Input

On entry: the stride separating matrix row elements in the array
x.
Constraint:
${\mathbf{pdx}}\ge {\mathbf{n}}$.

6:
$\mathbf{y}\left[\mathit{dim}\right]$ – Complex
Input/Output

Note: the dimension,
dim, of the array
y
must be at least
${\mathbf{pdy}}\times {\mathbf{t}}$.
The $\left(i,j\right)$th element of the matrix $Y$ is stored in ${\mathbf{y}}\left[\left(j1\right)\times {\mathbf{pdy}}+i1\right]$.
On initial entry: need not be set.
On intermediate exit:
if ${\mathbf{irevcm}}=2$, contains the current matrix $Y$.
On intermediate reentry: if ${\mathbf{irevcm}}=1$, must contain $AX$.
On final exit: the array is undefined.

7:
$\mathbf{pdy}$ – Integer
Input

On entry: the stride separating matrix row elements in the array
y.
Constraint:
${\mathbf{pdy}}\ge {\mathbf{m}}$.

8:
$\mathbf{estnrm}$ – double *
Input/Output

On initial entry: need not be set.
On intermediate reentry: must not be changed.
On final exit: an estimate (a lower bound) for ${\Vert A\Vert}_{1}$.

9:
$\mathbf{t}$ – Integer
Input

On entry: the number of columns
$t$ of the matrices
$X$ and
$Y$. This is an 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.
Suggested value:
${\mathbf{t}}=2$.
Constraint:
$1\le {\mathbf{t}}\le {\mathbf{m}}$.

10:
$\mathbf{seed}$ – Integer
Input

On entry: the seed used for random number generation.
If
${\mathbf{t}}=1$,
seed is not used.
Constraint:
if ${\mathbf{t}}>1$, ${\mathbf{seed}}\ge 1$.

11:
$\mathbf{work}\left[{\mathbf{m}}\times {\mathbf{t}}\right]$ – Complex
Communication Array

12:
$\mathbf{rwork}\left[2\times {\mathbf{n}}\right]$ – double
Communication Array

13:
$\mathbf{iwork}\left[2\times {\mathbf{n}}+5\times {\mathbf{t}}+20\right]$ – Integer
Communication Array

On initial entry: need not be set.
On intermediate reentry: must not be changed.

14:
$\mathbf{fail}$ – NagError *
Input/Output

The NAG error argument (see
Section 7 in the Introduction to the NAG Library CL Interface).
6
Error Indicators and Warnings
 NE_ALLOC_FAIL

Dynamic memory allocation failed.
See
Section 3.1.2 in the Introduction to the NAG Library CL Interface for further information.
 NE_BAD_PARAM

On entry, argument $\u2329\mathit{\text{value}}\u232a$ had an illegal value.
 NE_INT

On entry, ${\mathbf{irevcm}}=\u2329\mathit{\text{value}}\u232a$.
Constraint: ${\mathbf{irevcm}}=0$, $1$ or $2$.
On entry, ${\mathbf{m}}=\u2329\mathit{\text{value}}\u232a$.
Constraint: ${\mathbf{m}}\ge 0$.
On entry, ${\mathbf{n}}=\u2329\mathit{\text{value}}\u232a$.
Constraint: ${\mathbf{n}}\ge 0$.
On initial entry, ${\mathbf{irevcm}}=\u2329\mathit{\text{value}}\u232a$.
Constraint: ${\mathbf{irevcm}}=0$.
 NE_INT_2

On entry, ${\mathbf{m}}=\u2329\mathit{\text{value}}\u232a$ and ${\mathbf{t}}=\u2329\mathit{\text{value}}\u232a$.
Constraint: $1\le {\mathbf{t}}\le {\mathbf{m}}$.
On entry, ${\mathbf{pdx}}=\u2329\mathit{\text{value}}\u232a$ and ${\mathbf{n}}=\u2329\mathit{\text{value}}\u232a$.
Constraint: ${\mathbf{pdx}}\ge {\mathbf{n}}$.
On entry, ${\mathbf{pdy}}=\u2329\mathit{\text{value}}\u232a$ and ${\mathbf{m}}=\u2329\mathit{\text{value}}\u232a$.
Constraint: ${\mathbf{pdy}}\ge {\mathbf{m}}$.
On entry, ${\mathbf{t}}=\u2329\mathit{\text{value}}\u232a$ and ${\mathbf{seed}}=\u2329\mathit{\text{value}}\u232a$.
Constraint: if ${\mathbf{t}}>1$, ${\mathbf{seed}}\ge 1$.
 NE_INTERNAL_ERROR

An internal error has occurred in this function. Check the function call and any array sizes. If the call is correct then please contact
NAG for assistance.
See
Section 7.5 in the Introduction to the NAG Library CL Interface for further information.
 NE_NO_LICENCE

Your licence key may have expired or may not have been installed correctly.
See
Section 8 in the Introduction to the NAG Library CL Interface for further information.
7
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
${\Vert A\Vert}_{1}$; often the estimate has many correct figures. However, matrices exist for which the estimate is smaller than
${\Vert A\Vert}_{1}$ by an arbitrary factor; such matrices are very unlikely to arise in practice. See
Higham and Tisseur (2000) for further details.
8
Parallelism and Performance
f04zdc is threaded by NAG for parallel execution in multithreaded implementations of the NAG Library.
Please consult the
X06 Chapter Introduction for information on how to control and interrogate the OpenMP environment used within this function. Please also consult the
Users' Note for your implementation for any additional implementationspecific information.
For most problems the time taken during calls to f04zdc will be negligible compared with the time spent evaluating matrix products between calls to f04zdc.
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$.
It is your responsibility to guard against potential overflows during evaluation of the matrix products. In particular, when estimating ${\Vert {B}^{1}\Vert}_{1}$ using a triangular factorization of $B$, f04zdc should not be called if one of the factors is exactly singular – otherwise division by zero may occur in the substitutions.
9.3
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.
To estimate the
$1$norm of the inverse of a matrix
$A$, the following skeleton code can normally be used:
do {
f04zdc(&irevcm,m,n,x,pdx,y,pdy,&estnrm,t,seed,work,rwork,iwork,&fail);
if (irevcm == 1){
.. Code to compute y = A^(1) x ..
}
else if (irevcm == 2){
.. Code to compute x = A^(H) y ..
}
} (while irevcm != 0)
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
Section 10 illustrates how
f04zdc can be used in conjunction with NAG Library function for
$LU$ factorization of complex matrices
f07arc).
It is also straightforward to use
f04zdc for Hermitian positive definite matrices, using
f16tfc,
f07frc and
f07fsc for factorization and solution.
For upper or lower triangular square matrices, no factorization function is needed:
$Y={A}^{1}X$ and
$X={A}^{\mathrm{H}}Y$ may be computed by calls to
f16sjc (or
f16skc if the matrix is banded, or
f16slc if the matrix is stored in packed form).
10
Example
This example estimates the condition number
${\Vert A\Vert}_{1}{\Vert {A}^{1}\Vert}_{1}$ of the matrix
$A$ given by
10.1
Program Text
10.2
Program Data
10.3
Program Results