# NAG Library Function Document

## 1Purpose

nag_linsys_real_gen_norm_rcomm (f04ydc) estimates the $1$-norm of a real 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.

## 2Specification

 #include #include
 void nag_linsys_real_gen_norm_rcomm (Integer *irevcm, Integer m, Integer n, double x[], Integer pdx, double y[], Integer pdy, double *estnrm, Integer t, Integer seed, double work[], Integer iwork[], NagError *fail)

## 3Description

nag_linsys_real_gen_norm_rcomm (f04ydc) 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$ 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 $n×t$ matrix $X$ (with $t\ll n$) or an $m×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 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}$ 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{T}}Y$ required by nag_linsys_real_gen_norm_rcomm (f04ydc) 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_gen_norm_rcomm (f04ydc) 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).

## 4References

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

## 5Arguments

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.
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×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{T}}Y$ and store the result in x, (b) call nag_linsys_real_gen_norm_rcomm (f04ydc) once again, with all the other arguments unchanged.
On intermediate re-entry: irevcm must be unchanged.
On final exit: ${\mathbf{irevcm}}=0$.
Note: any values you return to nag_linsys_real_gen_norm_rcomm (f04ydc) as part of the reverse communication procedure should not include floating-point NaN (Not a Number) or infinity values, since these are not handled by nag_linsys_real_gen_norm_rcomm (f04ydc). If your code inadvertently does return any NaNs or infinities, nag_linsys_real_gen_norm_rcomm (f04ydc) is likely to produce unexpected results.
2:    $\mathbf{m}$IntegerInput
On entry: the number of rows of the matrix $A$.
Constraint: ${\mathbf{m}}\ge 0$.
3:    $\mathbf{n}$IntegerInput
On entry: $n$, the number of columns of the matrix $A$.
Constraint: ${\mathbf{n}}\ge 0$.
4:    $\mathbf{x}\left[\mathit{dim}\right]$doubleInput/Output
Note: the dimension, dim, of the array x must be at least ${\mathbf{pdx}}×{\mathbf{t}}$.
The $\left(i,j\right)$th element of the matrix $X$ is stored in ${\mathbf{x}}\left[\left(j-1\right)×{\mathbf{pdx}}+i-1\right]$.
On initial entry: need not be set.
On intermediate exit: if ${\mathbf{irevcm}}=1$, contains the current matrix $X$.
On intermediate re-entry: if ${\mathbf{irevcm}}=2$, must contain ${A}^{\mathrm{T}}Y$.
On final exit: the array is undefined.
5:    $\mathbf{pdx}$IntegerInput
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]$doubleInput/Output
Note: the dimension, dim, of the array y must be at least ${\mathbf{pdy}}×{\mathbf{t}}$.
The $\left(i,j\right)$th element of the matrix $Y$ is stored in ${\mathbf{y}}\left[\left(j-1\right)×{\mathbf{pdy}}+i-1\right]$.
On initial entry: need not be set.
On intermediate exit: if ${\mathbf{irevcm}}=2$, contains the current matrix $Y$.
On intermediate re-entry: if ${\mathbf{irevcm}}=1$, must contain $AX$.
On final exit: the array is undefined.
7:    $\mathbf{pdy}$IntegerInput
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 re-entry: must not be changed.
On final exit: an estimate (a lower bound) for ${‖A‖}_{1}$.
9:    $\mathbf{t}$IntegerInput
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}$IntegerInput
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}}×{\mathbf{t}}\right]$doubleCommunication Array
12:  $\mathbf{iwork}\left[2×{\mathbf{n}}+5×{\mathbf{t}}+20\right]$IntegerCommunication Array
On initial entry: need not be set.
On intermediate re-entry: must not be changed.
13:  $\mathbf{fail}$NagError *Input/Output
The NAG error argument (see Section 3.7 in How to Use the NAG Library and its Documentation).

## 6Error Indicators and Warnings

NE_ALLOC_FAIL
Dynamic memory allocation failed.
See Section 2.3.1.2 in How to Use the NAG Library and its Documentation for further information.
On entry, argument $〈\mathit{\text{value}}〉$ had an illegal value.
NE_INT
On entry, ${\mathbf{irevcm}}=〈\mathit{\text{value}}〉$.
Constraint: ${\mathbf{irevcm}}=0$, $1$ or $2$.
On entry, ${\mathbf{m}}=〈\mathit{\text{value}}〉$.
Constraint: ${\mathbf{m}}\ge 0$.
On entry, ${\mathbf{n}}=〈\mathit{\text{value}}〉$.
Constraint: ${\mathbf{n}}\ge 0$.
On initial entry, ${\mathbf{irevcm}}=〈\mathit{\text{value}}〉$.
Constraint: ${\mathbf{irevcm}}=0$.
NE_INT_2
On entry, ${\mathbf{m}}=〈\mathit{\text{value}}〉$ and ${\mathbf{t}}=〈\mathit{\text{value}}〉$.
Constraint: $1\le {\mathbf{t}}\le {\mathbf{m}}$.
On entry, ${\mathbf{pdx}}=〈\mathit{\text{value}}〉$ and ${\mathbf{n}}=〈\mathit{\text{value}}〉$.
Constraint: ${\mathbf{pdx}}\ge {\mathbf{n}}$.
On entry, ${\mathbf{pdy}}=〈\mathit{\text{value}}〉$ and ${\mathbf{m}}=〈\mathit{\text{value}}〉$.
Constraint: ${\mathbf{pdy}}\ge {\mathbf{m}}$.
On entry, ${\mathbf{t}}=〈\mathit{\text{value}}〉$ and ${\mathbf{seed}}=〈\mathit{\text{value}}〉$.
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 2.7.6 in How to Use the NAG Library and its Documentation for further information.
NE_NO_LICENCE
Your licence key may have expired or may not have been installed correctly.
See Section 2.7.5 in How to Use the NAG Library and its Documentation for further information.

## 7Accuracy

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.

## 8Parallelism and Performance

nag_linsys_real_gen_norm_rcomm (f04ydc) 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 implementation-specific information.

### 9.1Timing

For most problems the time taken during calls to nag_linsys_real_gen_norm_rcomm (f04ydc) will be negligible compared with the time spent evaluating matrix products between calls to nag_linsys_real_gen_norm_rcomm (f04ydc).
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.2Overflow

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_real_gen_norm_rcomm (f04ydc) should not be called if one of the factors is exactly singular – otherwise division by zero may occur in the substitutions.

### 9.3Choice 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.4Use 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:
```do {
f04ydc(&irevcm,m,n,x,pdx,y,pdy,&estnrm,t,seed,work,iwork,&fail);
if (irevcm == 1){
.. Code to compute y = A^(-1) x ..
}
else if  (irevcm == 2){
.. Code to compute x = A^(-T) y ..
}
} (while irevcm != 0)```
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 code will vary, depending on the type of the matrix $A$, and the NAG function used to factorize $A$.
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 Section 10 illustrates how nag_linsys_real_gen_norm_rcomm (f04ydc) can be used in conjunction with NAG C Library functions for $LU$ factorization of a real matrix nag_dgetrf (f07adc).
It is straightforward to use nag_linsys_real_gen_norm_rcomm (f04ydc) for the following other types of matrix, using the named functions for factorization and solution:
For upper or lower triangular matrices, no factorization function is needed: $Y={A}^{-1}X$ and $X={A}^{-\mathrm{T}}Y$ may be computed by calls to nag_dtrsv (f16pjc) (or nag_dtbsv (f16pkc) if the matrix is banded, or nag_dtpsv (f16plc) if the matrix is stored in packed form).

## 10Example

This example estimates the condition number ${‖A‖}_{1}{‖{A}^{-1}‖}_{1}$ of the matrix $A$ given by
 $A= 0.7 -0.2 1.0 0.0 2.0 0.1 0.3 0.7 0.0 1.0 0.9 0.2 0.0 0.0 0.2 0.7 0.0 -1.1 0.0 3.4 -0.7 0.2 0.1 0.1 0.0 -4.0 0.0 1.0 9.0 0.0 0.4 1.2 4.3 0.0 6.2 5.9 .$

### 10.1Program Text

Program Text (f04ydce.c)

### 10.2Program Data

Program Data (f04ydce.d)

### 10.3Program Results

Program Results (f04ydce.r)

© The Numerical Algorithms Group Ltd, Oxford, UK. 2017