NAG Library Function Document

nag_linsys_real_gen_norm_rcomm (f04ydc)

-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.

2 Specification

#include <nag.h>

#include <nagf04.h>

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)

3 Description

nag_linsys_real_gen_norm_rcomm (f04ydc) computes an estimate (a lower bound) for the

1

-norm

{‖A‖}_{1} = \max_{1 \leq j \leq n} \sum_{i = 1}^{m} |a_{i j}|

(1)

of an

m

n

real matrix

A = (a_{i j})

. The function regards the matrix

A

as being defined by a user-supplied ‘Black Box’ which, given an

n \times t

matrix

X

(with

t ≪ n

) or an

m \times t

matrix

Y

, can return

A X

A^{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

κ_{1} (B) = {‖B‖}_{1} {‖B^{- 1}‖}_{1}

, without forming

B^{- 1}

explicitly (

A = B^{- 1}

above).

If, for example, an

L U

factorization of

B

is available, the matrix products

B^{- 1} X

and

B^{- 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

A B C

, without forming the products explicitly. Further applications are described by Higham (1988).

Since

{‖A‖}_{\infty} = {‖A^{T}‖}_{1}

, nag_linsys_real_gen_norm_rcomm (f04ydc) can be used to estimate the

\infty

-norm of

A

by working with

A^{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 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

5 Arguments

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: irevcm – Integer *Input/Output

On initial entry: must be set to

0

On intermediate exit:

irevcm = 1

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 $irevcm = 1$ , evaluate $A X$ and store the result in y or if $irevcm = 2$ , evaluate $A^{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:

irevcm = 0

2: m – IntegerInput

On entry: the number of rows of the matrix

A

Constraint:

m \geq 0

3: n – IntegerInput

On entry:

n

, the number of columns of the matrix

A

Constraint:

n \geq 0

4: x[ $\dim$ ] – doubleInput/Output

Note: the dimension, dim, of the array x must be at least

pdx \times t

The

(i, j)

th element of the matrix

X

is stored in

x [(j - 1) \times pdx + i - 1]

On initial entry: need not be set.

On intermediate exit: if

irevcm = 1

, contains the current matrix

X

On intermediate re-entry: if

irevcm = 2

, must contain

A^{T} Y

On final exit: the array is undefined.

5: pdx – IntegerInput

On entry: the stride separating matrix row elements in the array x.

Constraint:

pdx \geq n

6: y[ $\dim$ ] – doubleInput/Output

Note: the dimension, dim, of the array y must be at least

pdy \times t

The

(i, j)

th element of the matrix

Y

is stored in

y [(j - 1) \times pdy + i - 1]

On initial entry: need not be set.

On intermediate exit: if

irevcm = 2

, contains the current matrix

Y

On intermediate re-entry: if

irevcm = 1

, must contain

A X

On final exit: the array is undefined.

7: pdy – IntegerInput

On entry: the stride separating matrix row elements in the array y.

Constraint:

pdy \geq m

8: 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: 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.

t \geq 2

then a partly random starting matrix is used in the algorithm.

Suggested value:

t = 2

Constraint:

1 \leq t \leq m

10: seed – IntegerInput

On entry: the seed used for random number generation.

t = 1

, seed is not used.

Constraint: if

t > 1

seed \geq 1

11: work[ $m \times t$ ] – doubleCommunication Array

12: iwork[ $2 \times n + 5 \times t + 20$ ] – IntegerCommunication Array

On initial entry: need not be set.

On intermediate re-entry: must not be changed.

13: fail – NagError *Input/Output

The NAG error argument (see Section 3.6 in the Essential Introduction).

6 Error Indicators and Warnings

NE_BAD_PARAM: On entry, argument $⟨value⟩$ had an illegal value.
NE_INT: On entry, $irevcm = ⟨value⟩$ .
Constraint: $irevcm = 0$ , $1$ or $2$ .
On entry, $m = ⟨value⟩$ .
Constraint: $m \geq 0$ .
On entry, $n = ⟨value⟩$ .
Constraint: $n \geq 0$ .
On initial entry, $irevcm = ⟨value⟩$ .
Constraint: $irevcm = 0$ .
NE_INT_2: On entry, $m = ⟨value⟩$ and $t = ⟨value⟩$ .
Constraint: $1 \leq t \leq m$ .
On entry, $pdx = ⟨value⟩$ and $n = ⟨value⟩$ .
Constraint: $pdx \geq n$ .
On entry, $pdy = ⟨value⟩$ and $m = ⟨value⟩$ .
Constraint: $pdy \geq m$ .
On entry, $t = ⟨value⟩$ and $seed = ⟨value⟩$ .
Constraint: if $t > 1$ , $seed \geq 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.

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

{‖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.

8 Parallelism and Performance

Not applicable.

9 Further Comments

9.1 Timing

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 = A X

and five products of the form

X = A^{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

{‖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.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 \geq 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:

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

A^{- T} Y

, solve the equation

A Y = X

A^{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

L U

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:

nonsymmetric band (nag_dgbtrf (f07bdc) and nag_dgbtrs (f07bec));
symmetric positive definite (nag_dpotrf (f07fdc) and nag_dpotrs (f07fec));
symmetric positive definite band (nag_dpbtrf (f07hdc) and nag_dpbtrs (f07hec));
symmetric positive definite tridiagonal (nag_dptsv (f07jac), nag_dpttrf (f07jdc) and nag_dpttrs (f07jec));
symmetric positive definite variable bandwidth (nag_real_cholesky_skyline (f01mcc) and nag_real_cholesky_skyline_solve (f04mcc));
symmetric positive definite sparse (nag_sparse_sym_chol_fac (f11jac) and nag_sparse_sym_precon_ichol_solve (f11jbc));
symmetric indefinite (nag_dsptrf (f07pdc) and nag_dsptrs (f07pec));
nonsymmetric sparse (nag_superlu_lu_factorize (f11mec) and nag_superlu_solve_lu (f11mfc); note that nag_superlu_condition_number_lu (f11mgc) can also be used here).

For upper or lower triangular matrices, no factorization function is needed:

Y = A^{- 1} X

and

X = A^{- 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).

10 Example

This example estimates the condition number

{‖A‖}_{1} {‖A^{- 1}‖}_{1}

of the matrix

A

given by

A = (\begin{matrix} 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 \end{matrix}) .

NAG Library Function Documentnag_linsys_real_gen_norm_rcomm (f04ydc)

+− Contents

1 Purpose

2 Specification

3 Description

4 References

5 Arguments

6 Error Indicators and Warnings

7 Accuracy

8 Parallelism and Performance

9 Further Comments

9.1 Timing

9.2 Overflow

9.3 Choice of t

9.4 Use in Conjunction with NAG Library Routines

10 Example

10.1 Program Text

10.2 Program Data

10.3 Program Results

NAG Library Function Document

nag_linsys_real_gen_norm_rcomm (f04ydc)

9.3 Choice of $t$