NAG FL Interface
f04ydf (real_gen_norm_rcomm)
1
Purpose
f04ydf estimates the $1$norm of a real rectangular matrix without accessing the matrix explicitly. It uses reverse communication for evaluating matrix products. The routine may be used for estimating condition numbers of square matrices.
2
Specification
Fortran Interface
Subroutine f04ydf ( 
irevcm, m, n, x, ldx, y, ldy, estnrm, t, seed, work, iwork, ifail) 
Integer, Intent (In) 
:: 
m, n, ldx, ldy, t, seed 
Integer, Intent (Inout) 
:: 
irevcm, iwork(2*n+5*t+20), ifail 
Real (Kind=nag_wp), Intent (Inout) 
:: 
x(ldx,*), y(ldy,*), estnrm, work(m*t) 

C Header Interface
#include <nag.h>
void 
f04ydf_ (Integer *irevcm, const Integer *m, const Integer *n, double x[], const Integer *ldx, double y[], const Integer *ldy, double *estnrm, const Integer *t, const Integer *seed, double work[], Integer iwork[], Integer *ifail) 

C++ Header Interface
#include <nag.h> extern "C" {
void 
f04ydf_ (Integer &irevcm, const Integer &m, const Integer &n, double x[], const Integer &ldx, double y[], const Integer &ldy, double &estnrm, const Integer &t, const Integer &seed, double work[], Integer iwork[], Integer &ifail) 
}

The routine may be called by the names f04ydf or nagf_linsys_real_gen_norm_rcomm.
3
Description
f04ydf computes an estimate (a lower bound) for the
$1$norm
of an
$m$ by
$n$ real matrix
$A=\left({a}_{ij}\right)$. The routine 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{T}}Y$. A reverse communication interface is used; thus control is returned to the calling program whenever a matrix product is required.
Note: this routine 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 routine 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{T}}Y$ required by f04ydf may be computed by back and forwardsubstitutions, without computing ${B}^{1}$.
The routine 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 ${\Vert A\Vert}_{\infty}={\Vert {A}^{\mathrm{T}}\Vert}_{1}$, f04ydf can be used to estimate the $\infty $norm of $A$ by working with ${A}^{\mathrm{T}}$ 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 routine 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{T}}Y$ and store the result in x,

(b)call f04ydf once again, with all the other arguments unchanged.
On intermediate reentry:
irevcm must be unchanged.
On final exit: ${\mathbf{irevcm}}=0$.
Note: any values you return to f04ydf as part of the reverse communication procedure should not include floatingpoint NaN (Not a Number) or infinity values, since these are not handled by f04ydf. If your code does inadvertently return any NaNs or infinities, f04ydf 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 entry: $n$, the number of columns of the matrix $A$.
Constraint:
${\mathbf{n}}\ge 0$.

4:
$\mathbf{x}\left({\mathbf{ldx}},*\right)$ – Real (Kind=nag_wp) array
Input/Output

Note: the second dimension of the array
x
must be at least
${\mathbf{t}}$.
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{T}}Y$.
On final exit: the array is undefined.

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

On initial entry: the leading dimension of the array
x as declared in the (sub)program from which
f04ydf is called.
Constraint:
${\mathbf{ldx}}\ge {\mathbf{n}}$.

6:
$\mathbf{y}\left({\mathbf{ldy}},*\right)$ – Real (Kind=nag_wp) array
Input/Output

Note: the second dimension of the array
y
must be at least
${\mathbf{t}}$.
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{ldy}$ – Integer
Input

On initial entry: the leading dimension of the array
y as declared in the (sub)program from which
f04ydf is called.
Constraint:
${\mathbf{ldy}}\ge {\mathbf{m}}$.

8:
$\mathbf{estnrm}$ – Real (Kind=nag_wp)
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)$ – Real (Kind=nag_wp) array
Communication Array

12:
$\mathbf{iwork}\left(2\times {\mathbf{n}}+5\times {\mathbf{t}}+20\right)$ – Integer array
Communication Array

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

13:
$\mathbf{ifail}$ – Integer
Input/Output

On initial entry:
ifail must be set to
$0$,
$1$ or
$1$ to set behaviour on detection of an error; these values have no effect when no error is detected.
A value of $0$ causes the printing of an error message and program execution will be halted; otherwise program execution continues. A value of $1$ means that an error message is printed while a value of $1$ means that it is not.
If halting is not appropriate, the value
$1$ or
$1$ is recommended. If message printing is undesirable, then the value
$1$ is recommended. Otherwise, the value
$1$ is recommended since useful values can be provided in some output arguments even when
${\mathbf{ifail}}\ne {\mathbf{0}}$ on exit.
When the value $\mathbf{1}$ or $\mathbf{1}$ is used it is essential to test the value of ifail on exit.
On final exit:
${\mathbf{ifail}}={\mathbf{0}}$ unless the routine detects an error or a warning has been flagged (see
Section 6).
6
Error Indicators and Warnings
If on entry
${\mathbf{ifail}}=0$ or
$1$, explanatory error messages are output on the current error message unit (as defined by
x04aaf).
Errors or warnings detected by the routine:
 ${\mathbf{ifail}}=1$

Internal error; please contact
NAG.
 ${\mathbf{ifail}}=1$

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

On entry, ${\mathbf{m}}=\u2329\mathit{\text{value}}\u232a$.
Constraint: ${\mathbf{m}}\ge 0$.
 ${\mathbf{ifail}}=3$

On entry, ${\mathbf{n}}=\u2329\mathit{\text{value}}\u232a$.
Constraint: ${\mathbf{n}}\ge 0$.
 ${\mathbf{ifail}}=5$

On entry, ${\mathbf{ldx}}=\u2329\mathit{\text{value}}\u232a$ and ${\mathbf{n}}=\u2329\mathit{\text{value}}\u232a$.
Constraint: ${\mathbf{ldx}}\ge {\mathbf{n}}$.
 ${\mathbf{ifail}}=7$

On entry, ${\mathbf{ldy}}=\u2329\mathit{\text{value}}\u232a$ and ${\mathbf{m}}=\u2329\mathit{\text{value}}\u232a$.
Constraint: ${\mathbf{ldy}}\ge {\mathbf{m}}$.
 ${\mathbf{ifail}}=9$

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}}$.
 ${\mathbf{ifail}}=10$

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$.
 ${\mathbf{ifail}}=99$
An unexpected error has been triggered by this routine. Please
contact
NAG.
See
Section 7 in the Introduction to the NAG Library FL Interface for further information.
 ${\mathbf{ifail}}=399$
Your licence key may have expired or may not have been installed correctly.
See
Section 8 in the Introduction to the NAG Library FL Interface for further information.
 ${\mathbf{ifail}}=999$
Dynamic memory allocation failed.
See
Section 9 in the Introduction to the NAG Library FL 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
f04ydf 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 routine. Please also consult the
Users' Note for your implementation for any additional implementationspecific information.
9.1
Timing
For most problems the time taken during calls to f04ydf will be negligible compared with the time spent evaluating matrix products between calls to f04ydf.
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{T}}Y$ will be required. The number of iterations is independent of the choice of $t$.
9.2
Overflow
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$, f04ydf 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.
9.4
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) Then
irevcm = 0
10 Call f04ydf (irevcm,m,n,x,ldx,y,ldy,estnrm,t,seed,work, &
iwork,ifail)
If (irevcm.ne.0) Then
If (irevcm.eq.1) Then
... code to compute y=inv(A)x ...
Else
... code to compute x=inv(transpose(A))y ...
End If
Go To 10
End If
End If
To compute
${A}^{1}X$ or
${A}^{\mathrm{T}}Y$, solve the equation
$AY=X$ or
${A}^{\mathrm{T}}X=Y$, storing the result in
y or
x respectively.
The factorization will normally have been performed by a suitable routine from
Chapters F01,
F03 or
F07. Note also that many of the ‘Black Box’ routines in
Chapter F04 for solving systems of equations also return a factorization of the matrix.
It is straightforward to use
f04ydf for the following other types of matrix, using the named routines for factorization and solution:
For upper or lower triangular matrices, no factorization routine is needed:
$Y={A}^{1}X$ and
$X={A}^{\mathrm{T}}Y$ may be computed by calls to
f06pjf (or
f06pkf if the matrix is banded, or
f06plf if the matrix is stored in packed form).
10
Example
For this routine two examples are provided. There is a single example program for f04ydf, 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)
This example estimates the condition number
${\Vert A\Vert}_{1}{\Vert {A}^{1}\Vert}_{1}$ of the matrix
$A$ given by
Example 2 (EX2)
This example estimates the condition number of the sparse matrix
$A$ (stored in symmetric coordinate storage format) given by
10.1
Program Text
10.2
Program Data
10.3
Program Results