NAG Library Routine Document

F04YCF

-norm of a real matrix without accessing the matrix explicitly. It uses reverse communication for evaluating matrix-vector products. The routine may be used for estimating matrix condition numbers.

2 Specification

SUBROUTINE F04YCF (

ICASE, N, X, ESTNRM, WORK, IWORK, IFAIL)

INTEGER	ICASE, N, IWORK(N), IFAIL
REAL (KIND=nag_wp)	X(N), ESTNRM, WORK(N)

3 Description

F04YCF computes an estimate (a lower bound) for the

1

-norm

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

(1)

of an

n

n

real matrix

A = (a_{i j})

. The routine regards the matrix

A

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

x

, can return either of the matrix-vector products

A x

A^{T} x

. A reverse communication interface is used; thus control is returned to the calling program whenever a matrix-vector 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

{‖B^{- 1}‖}_{1}

, 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-vector products

B^{- 1} x

and

B^{- T} x

required by F04YCF 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 by Higham (1988).

Since

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

, F04YCF can be used to estimate the

\infty

-norm of

A

by working with

A^{T}

instead of

A

The algorithm used is based on a method given by Hager (1984) and is described by Higham (1988). A comparison of several techniques for condition number estimation is given by Higham (1987).

Note: F04YDF can also be used to estimate the norm of a real matrix. F04YDF uses a more recent algorithm than F04YCF and it is recommended that F04YDF be used in place of F04YCF.

4 References

Hager W W (1984) Condition estimates SIAM J. Sci. Statist. Comput. 5 311–316

Higham N J (1987) A survey of condition number estimation for triangular matrices SIAM Rev. 29 575–596

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

5 Parameters

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 parameter ICASE. Between intermediate exits and re-entries, all parameters other than X must remain unchanged.

1: $ICASE$ – INTEGERInput/Output

On initial entry: must be set to

0

On intermediate exit:

ICASE = 1

2

, and

X (i)

, for

i = 1, 2, \dots, n

, contain the elements of a vector

x

. The calling program must

(a)	evaluate $A x$ (if $ICASE = 1$ ) or $A^{T} x$ (if $ICASE = 2$ ),
(b)	place the result in X, and
(c)	call F04YCF once again, with all the other parameters unchanged.

On final exit:

ICASE = 0

2: $N$ – INTEGERInput

On initial entry:

n

, the order of the matrix

A

Constraint:

N \geq 1

3: $X (N)$ – REAL (KIND=nag_wp) arrayInput/Output

On initial entry: need not be set.

On intermediate exit: contains the current vector

x

On intermediate re-entry: must contain

A x

(if

ICASE = 1

) or

A^{T} x

(if

ICASE = 2

On final exit: the array is undefined.

4: $ESTNRM$ – REAL (KIND=nag_wp)Input/Output

On initial entry: need not be set.

On intermediate exit: should not be changed.

On final exit: an estimate (a lower bound) for

{‖A‖}_{1}

5: $WORK (N)$ – REAL (KIND=nag_wp) arrayInput/Output

On initial entry: need not be set.

On final exit: contains a vector

v

such that

v = A w

where

ESTNRM = {‖v‖}_{1} / {‖w‖}_{1}

(

w

is not returned). If

A = B^{- 1}

and ESTNRM is large, then

v

is an approximate null vector for

B

6: $IWORK (N)$ – INTEGER arrayCommunication Array

7: $IFAIL$ – INTEGERInput/Output

On initial entry: IFAIL must be set to

0

- 1 ​ or ​ 1

. If you are unfamiliar with this parameter you should refer to Section 3.3 in the Essential Introduction 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, because for this routine the values of the output parameters may be useful even if

IFAIL \neq 0

on exit, the recommended value is

- 1

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

On final 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$: On entry, $N < 1$ .

$IFAIL = - 99$: An unexpected error has been triggered by this routine. Please contact NAG.
See Section 3.8 in the Essential Introduction for further information.

$IFAIL = - 399$: Your licence key may have expired or may not have been installed correctly.
See Section 3.7 in the Essential Introduction for further information.

$IFAIL = - 999$: Dynamic memory allocation failed.
See Section 3.6 in the Essential Introduction for further information.

7 Accuracy

In extensive tests on random matrices of size up to

n = 100

the estimate ESTNRM has been found always to be within a factor eleven 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 (1988) for further details.

8 Parallelism and Performance

F04YCF is not threaded by NAG in any implementation.

F04YCF makes calls to BLAS and/or LAPACK routines, which may be threaded within the vendor library used by this implementation. Consult the documentation for the vendor library for further information.

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

The total time taken within F04YCF is proportional to

n

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

The number of matrix-vector products required varies from

4

11

(or is

1

n = 1

). In most cases

4

5

products are required; it is rare for more than

7

to be needed.

9.2 Overflow

It is your responsibility to guard against potential overflows during evaluation of the matrix-vector products. In particular, when estimating

{‖B^{- 1}‖}_{1}

using a triangular factorization of

B

, F04YCF should not be called if one of the factors is exactly singular – otherwise division by zero may occur in the substitutions.

9.3 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 
         ICASE = 0 
10       CALL F04YCF (ICASE,N,X,ESTNRM,WORK,IWORK,IFAIL) 
         IF (ICASE.NE.0) THEN 
            IF (ICASE.EQ.1) THEN 
               ...  code to compute inv(A)*x ...  
            ELSE 
               ...  code to compute inv(transpose(A))*x ...  
            END IF 
            GO TO 10 
         END IF 
      END IF

To compute

A^{- 1} x

A^{- T} x

, solve the equation

A y = x

A^{T} y = x

for

y

, overwriting

y

x

. The code will vary, depending on the type of the matrix

A

, and the NAG routine used to factorize

A

Note that if

A

is any type of symmetric matrix, then

A = A^{T}

, and the code following the call of F04YCF can be reduced to:

 IF (ICASE.NE.0) THEN 
    ...  code to compute inv(A)*x ...
    GO TO 10
 END IF

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. The example program in Section 10 illustrates how F04YCF can be used in conjunction with NAG Library routines for two important types of matrix: full nonsymmetric matrices (factorized by F07ADF (DGETRF)) and sparse nonsymmetric matrices (factorized by F01BRF).

It is straightforward to use F04YCF for the following other types of matrix, using the named routines for factorization and solution:

nonsymmetric tridiagonal (F01LEF and F04LEF);
nonsymmetric almost block-diagonal (F01LHF and F04LHF);
nonsymmetric band (F07BDF (DGBTRF) and F07BEF (DGBTRS));
symmetric positive definite (F03BFF, or F07FDF (DPOTRF) and F07FEF (DPOTRS));
symmetric positive definite band (F07HDF (DPBTRF) and F07HEF (DPBTRS));
symmetric positive definite tridiagonal (F07JAF (DPTSV), F07JDF (DPTTRF) and F07JEF (DPTTRS));
symmetric positive definite variable bandwidth (F01MCF and F04MCF);
symmetric positive definite sparse (F11JAF and F11JBF);
symmetric indefinite (F07PDF (DSPTRF) and F07PEF (DSPTRS)).

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

A^{- 1} x

and

A^{- T} x

may be computed by calls to F06PJF (DTRSV) (or F06PKF (DTBSV) if the matrix is banded, or F06PLF (DTPSV) if the matrix is stored in packed form).

10 Example

For this routine two examples are presented. There is a single example program for F04YCF, 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)

To estimate the condition number

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

of the matrix

A

given by

A = (\begin{matrix} 1.5 & 2.0 & 3.0 & - 2.1 & 0.3 \\ 2.5 & 3.0 & - 4.0 & 2.3 & - 1.1 \\ 3.5 & 4.0 & 0.5 & - 3.1 & - 1.4 \\ - 0.4 & - 3.2 & - 2.1 & 3.1 & 2.1 \\ 1.7 & 3.7 & 1.9 & - 2.2 & - 3.3 \end{matrix}) .

Example 2 (EX2)

To estimate the condition number

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

of the matrix

A

given by

A = (\begin{matrix} 5.0 & 0.0 & 0.0 & 0.0 & 0.0 & 0.0 \\ 0.0 & 0.0 & - 1.0 & 2.0 & 0.0 & 0.0 \\ 0.0 & 2.0 & 3.0 & 0.0 & 0.0 & 0.0 \\ - 2.0 & 0.0 & 0.0 & 1.0 & 1.0 & 0.0 \\ - 1.0 & 0.0 & 0.0 & - 1.0 & 2.0 & - 3.0 \\ - 1.0 & - 1.0 & 0.0 & 0.0 & 0.0 & 6.0 \end{matrix}) .

NAG Library Routine DocumentF04YCF

▸▿ Contents

1 Purpose

2 Specification

3 Description

4 References

5 Parameters

6 Error Indicators and Warnings

7 Accuracy

8 Parallelism and Performance

9 Further Comments

9.1 Timing

9.2 Overflow

9.3 Use in Conjunction with NAG Library Routines

10 Example

10.1 Program Text

10.2 Program Data

10.3 Program Results

NAG Library Routine Document

F04YCF