NAG Library Routine Document

-norm of a complex 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 f04zdf (

irevcm, m, n, x, ldx, y, ldy, estnrm, t, seed, work, rwork, iwork, ifail)

Integer, Intent (In)	::	m, n, ldx, ldy, t, seed
Integer, Intent (Inout)	::	irevcm, iwork(2n+5t+20), ifail
Real (Kind=nag_wp), Intent (Inout)	::	estnrm, rwork(2*n)
Complex (Kind=nag_wp), Intent (Inout)	::	x(ldx,), y(ldy,), work(m*t)

C Header Interface

#include <nagmk26.h>

void	f04zdf_ (Integer irevcm, const Integer m, const Integer n, Complex x[], const Integer ldx, Complex y[], const Integer ldy, double estnrm, const Integer t, const Integer seed, Complex work[], double rwork[], Integer iwork[], Integer *ifail)

3

Description

f04zdf 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 routine 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 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 the formula (1) above.

The main use of the routine 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 f04zdf may be computed by back- and forward-substitutions, without computing

B^{- 1}

The routine 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}

, f04zdf 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 routine 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$ – IntegerInput/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 f04zdf once again, with all the arguments unchanged.

On intermediate re-entry: irevcm must be unchanged.

On final exit:

irevcm = 0

Note: any values you return to f04zdf 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 f04zdf. If your code does inadvertently return any NaNs or infinities, f04zdf is likely to produce unexpected results.

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 (ldx *)$ – Complex (Kind=nag_wp) arrayInput/Output

Note: the second dimension of the array x must be at least

t

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: $ldx$ – IntegerInput

On initial entry: the leading dimension of the array x as declared in the (sub)program from which f04zdf is called.

Constraint:

ldx \geq n

6: $y (ldy *)$ – Complex (Kind=nag_wp) arrayInput/Output

Note: the second dimension of the array y must be at least

t

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: $ldy$ – IntegerInput

On initial entry: the leading dimension of the array y as declared in the (sub)program from which f04zdf is called.

Constraint:

ldy \geq m

8: $estnrm$ – Real (Kind=nag_wp)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)$ – Complex (Kind=nag_wp) arrayCommunication Array

12: $rwork (2 \times n)$ – Real (Kind=nag_wp) arrayCommunication Array

13: $iwork (2 \times n + 5 \times t + 20)$ – Integer arrayCommunication Array

On initial entry: need not be set.

On intermediate re-entry: must not be changed.

14: $ifail$ – IntegerInput/Output

On entry: ifail must be set to

0

- 1 or 1

. If you are unfamiliar with this argument you should refer to Section 3.4 in How to Use the NAG Library and its Documentation for details.

For environments where it might be inappropriate to halt program execution when an error is detected, the value

- 1 or 1

is recommended. If the output of error messages is undesirable, then the value

1

is recommended. Otherwise, if you are not familiar with this argument, the recommended value is

0

. When the value $- 1 or 1$ is used it is essential to test the value of ifail on exit.

On exit:

ifail = 0

unless the routine detects an error or a warning has been flagged (see Section 6).

6

Error Indicators and Warnings

If on entry

ifail = 0

- 1

, explanatory error messages are output on the current error message unit (as defined by x04aaf).

Errors or warnings detected by the routine:

$ifail = 1$: Internal error; please contact NAG.

$ifail = - 1$: On entry, $irevcm = 〈value〉$ .
Constraint: $irevcm = 0$ , $1$ or $2$ .

On initial entry, $irevcm = 〈value〉$ .
Constraint: $irevcm = 0$ .

$ifail = - 2$: On entry, $m = 〈value〉$ .
Constraint: $m \geq 0$ .

$ifail = - 3$: On entry, $n = 〈value〉$ .
Constraint: $n \geq 0$ .

$ifail = - 5$: On entry, $ldx = 〈value〉$ and $n = 〈value〉$ .
Constraint: $ldx \geq n$ .

$ifail = - 7$: On entry, $ldy = 〈value〉$ and $m = 〈value〉$ .
Constraint: $ldy \geq m$ .

$ifail = - 9$: On entry, $m = 〈value〉$ and $t = 〈value〉$ .
Constraint: $1 \leq t \leq m$ .

$ifail = - 10$: On entry, $t = 〈value〉$ and $seed = 〈value〉$ .
Constraint: if $t > 1$ , $seed \geq 1$ .

$ifail = - 99$: An unexpected error has been triggered by this routine. Please contact NAG.
See Section 3.9 in How to Use the NAG Library and its Documentation for further information.

$ifail = - 399$: Your licence key may have expired or may not have been installed correctly.
See Section 3.8 in How to Use the NAG Library and its Documentation for further information.

$ifail = - 999$: Dynamic memory allocation failed.
See Section 3.7 in How to Use the NAG Library and its Documentation 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

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

f04zdf 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 implementation-specific information.

9

Further Comments

9.1

Timing

For most problems the time taken during calls to f04zdf will be negligible compared with the time spent evaluating matrix products between calls to f04zdf.

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

, f04zdf 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:

      ...  code to factorize A ...
      If (A is not singular) Then
         irevcm = 0
10       Call f04zdf (irevcm,m,n,x,ldx,y,ldy,estnrm,t,seed,work,  &
                      rwork,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(herm(A))y ...
            End If
            Go To 10
         End If
      End If

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 routine used to factorize

A

The example program in Section 10 illustrates how f04zdf can be used in conjunction with NAG Library routine for

L U

factorization of complex matrices f07arf (zgetrf)).

It is also straightforward to use f04zdf for Hermitian positive definite matrices, using f06tff, f07frf (zpotrf) and f07fsf (zpotrs) for factorization and solution.

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

Y = A^{- 1} X

and

X = A^{- H} Y

may be computed by calls to f06sjf (ztrsv) (or f06skf (ztbsv) if the matrix is banded, or f06slf (ztpsv) 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 Routine Document

f04zdf (complex_gen_norm_rcomm)

▸▿ 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

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