NAG Library Function Document

nag_linsys_complex_gen_norm_rcomm (f04zdc)

-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

#include <nag.h>

#include <nagf04.h>

void	nag_linsys_complex_gen_norm_rcomm (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)

3 Description

nag_linsys_complex_gen_norm_rcomm (f04zdc) 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

complex 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^{H} Y

, where

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

κ_{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^{- H} Y

required by nag_linsys_complex_gen_norm_rcomm (f04zdc) 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 in Higham (1988).

Since

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

, nag_linsys_complex_gen_norm_rcomm (f04zdc) can be used to estimate the

\infty

-norm of

A

by working with

A^{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 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^{H} Y$ and store the result in x, where $A^{H}$ is the complex conjugate transpose;
(b)	call nag_linsys_complex_gen_norm_rcomm (f04zdc) once again, with all the 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 initial entry:

n

, the number of columns of the matrix

A

Constraint:

n \geq 0

4: x[ $\dim$ ] – ComplexInput/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^{H} 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$ ] – ComplexInput/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$ ] – ComplexCommunication Array

12: rwork[ $2 \times n$ ] – doubleCommunication Array

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

14: 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_complex_gen_norm_rcomm (f04zdc) will be negligible compared with the time spent evaluating matrix products between calls to nag_linsys_complex_gen_norm_rcomm (f04zdc).

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^{H} 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_complex_gen_norm_rcomm (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 \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 {
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

A^{- H} Y

, solve the equation

A Y = X

A^{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 nag_linsys_complex_gen_norm_rcomm (f04zdc) can be used in conjunction with NAG C Library function for

L U

factorization of complex matrices nag_zgetrf (f07arc)).

It is also straightforward to use nag_linsys_complex_gen_norm_rcomm (f04zdc) for Hermitian positive definite matrices, using nag_zge_copy (f16tfc), nag_zpotrf (f07frc) and nag_zpotrs (f07fsc) for factorization and solution.

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

Y = A^{- 1} X

and

X = A^{- H} Y

may be computed by calls to nag_ztrsv (f16sjc) (or nag_ztbsv (f16skc) if the matrix is banded, or nag_ztpsv (f16slc) 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{array}{r} 0.7 + 0.1 i & - 0.2 + 0.0 i & 1.0 + 0.0 i & 0.0 + 0.0 i & 0.0 + 0.0 i & 0.1 + 0.0 i \\ 0.3 + 0.0 i & 0.7 + 0.0 i & 0.0 + 0.0 i & 1.0 + 0.2 i & 0.9 + 0.0 i & 0.2 + 0.0 i \\ 0.0 + 5.9 i & 0.0 + 0.0 i & 0.2 + 0.0 i & 0.7 + 0.0 i & 0.4 + 6.1 i & 1.1 + 0.4 i \\ 0.0 + 0.1 i & 0.0 + 0.1 i & - 0.7 + 0.0 i & 0.2 + 0.0 i & 0.1 + 0.0 i & 0.1 + 0.0 i \\ 0.0 + 0.0 i & 4.0 + 0.0 i & 0.0 + 0.0 i & 1.0 + 0.0 i & 9.0 + 0.0 i & 0.0 + 0.1 i \\ 4.5 + 6.7 i & 0.1 + 0.4 i & 0.0 + 3.2 i & 1.2 + 0.0 i & 0.0 + 0.0 i & 7.8 + 0.2 i \end{array}) .

NAG Library Function Documentnag_linsys_complex_gen_norm_rcomm (f04zdc)

+− 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_complex_gen_norm_rcomm (f04zdc)

9.3 Choice of $t$