hide long namesshow long names
hide short namesshow short names
Integer type:  int32  int64  nag_int  show int32  show int32  show int64  show int64  show nag_int  show nag_int

PDF version (NAG web site, 64-bit version, 64-bit version)
Chapter Contents
Chapter Introduction
NAG Toolbox

NAG Toolbox: nag_linsys_real_norm_rcomm (f04yc)

Purpose

nag_linsys_real_norm_rcomm (f04yc) estimates the 11-norm of a real matrix without accessing the matrix explicitly. It uses reverse communication for evaluating matrix-vector products. The function may be used for estimating matrix condition numbers.
Note: this function is scheduled to be withdrawn, please see f04yc in Advice on Replacement Calls for Withdrawn/Superseded Routines..

Syntax

[icase, x, estnrm, work, iwork, ifail] = f04yc(icase, x, estnrm, work, iwork, 'n', n)
[icase, x, estnrm, work, iwork, ifail] = nag_linsys_real_norm_rcomm(icase, x, estnrm, work, iwork, 'n', n)

Description

nag_linsys_real_norm_rcomm (f04yc) computes an estimate (a lower bound) for the 11-norm
n
A1 = max |aij|
1jn i = 1
A1 = max 1jn i=1n |aij|
(1)
of an nn by nn real matrix A = (aij)A=(aij). The function regards the matrix AA as being defined by a user-supplied ‘Black Box’ which, given an input vector xx, can return either of the matrix-vector products AxAx or ATxATx. A reverse communication interface is used; thus control is returned to the calling program whenever a matrix-vector product is required.
Note:  this function is not recommended for use when the elements of AA are known explicitly; it is then more efficient to compute the 11-norm directly from formula (1) above.
The main use of the function is for estimating B11B-11, and hence the condition number κ1(B) = B1B11κ1(B)=B1B-11, without forming B1B-1 explicitly (A = B1A=B-1 above).
If, for example, an LULU factorization of BB is available, the matrix-vector products B1xB-1x and BTxB-Tx required by nag_linsys_real_norm_rcomm (f04yc) may be computed by back- and forward-substitutions, without computing B1B-1.
The function can also be used to estimate 11-norms of matrix products such as A1BA-1B and ABCABC, without forming the products explicitly. Further applications are described by Higham (1988).
Since A = AT1A=AT1, nag_linsys_real_norm_rcomm (f04yc) can be used to estimate the -norm of AA by working with ATAT instead of AA.
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: nag_linsys_real_gen_norm_rcomm (f04yd) can also be used to estimate the norm of a real matrix. nag_linsys_real_gen_norm_rcomm (f04yd) uses a more recent algorithm than nag_linsys_real_norm_rcomm (f04yc) and it is recommended that nag_linsys_real_gen_norm_rcomm (f04yd) be used in place of nag_linsys_real_norm_rcomm (f04yc).

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

Parameters

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

Compulsory Input Parameters

1:     icase – int64int32nag_int scalar
On initial entry: must be set to 00.
2:     x(n) – double array
n, the dimension of the array, must satisfy the constraint n1n1.
On initial entry: need not be set.
n, the dimension of the array, must satisfy the constraint n1n1.
On intermediate re-entry: must contain AxAx (if icase = 1icase=1) or ATxATx (if icase = 2icase=2).
3:     estnrm – double scalar
On initial entry: need not be set.
4:     work(n) – double array
n, the dimension of the array, must satisfy the constraint n1n1.
On initial entry: need not be set.
5:     iwork(n) – int64int32nag_int array

Optional Input Parameters

1:     n – int64int32nag_int scalar
Default: The dimension of the arrays x, work, iwork. (An error is raised if these dimensions are not equal.)
On initial entry: nn, the order of the matrix AA.
Constraint: n1n1.

Input Parameters Omitted from the MATLAB Interface

None.

Output Parameters

1:     icase – int64int32nag_int scalar
On intermediate exit: icase = 1icase=1 or 22, and x(i)xi, for i = 1,2,,ni=1,2,,n, contain the elements of a vector xx. The calling program must
(a) evaluate AxAx (if icase = 1icase=1) or ATxATx (if icase = 2icase=2),
(b) place the result in x, and
(c) call nag_linsys_real_norm_rcomm (f04yc) once again, with all the other parameters unchanged.
On final exit: icase = 0icase=0.
2:     x(n) – double array
On intermediate exit: contains the current vector xx.
On final exit: the array is undefined.
3:     estnrm – double scalar
On intermediate exit: should not be changed.
On final exit: an estimate (a lower bound) for A1A1.
4:     work(n) – double array
On final exit: contains a vector vv such that v = Awv=Aw where estnrm = v1 / w1estnrm=v1/w1 (ww is not returned). If A = B1A=B-1 and estnrm is large, then vv is an approximate null vector for BB.
5:     iwork(n) – int64int32nag_int array
6:     ifail – int64int32nag_int scalar
ifail = 0ifail=0 unless the function detects an error (see [Error Indicators and Warnings]).

Error Indicators and Warnings

Errors or warnings detected by the function:
  ifail = 1ifail=1
On entry, n < 1n<1.

Accuracy

In extensive tests on random matrices of size up to n = 100n=100 the estimate estnrm has been found always to be within a factor eleven of A1A1; often the estimate has many correct figures. However, matrices exist for which the estimate is smaller than A1A1 by an arbitrary factor; such matrices are very unlikely to arise in practice. See Higham (1988) for further details.

Further Comments

Timing

The total time taken within nag_linsys_real_norm_rcomm (f04yc) is proportional to nn. For most problems the time taken during calls to nag_linsys_real_norm_rcomm (f04yc) will be negligible compared with the time spent evaluating matrix-vector products between calls to nag_linsys_real_norm_rcomm (f04yc).
The number of matrix-vector products required varies from 44 to 1111 (or is 11 if n = 1n=1). In most cases 44 or 55 products are required; it is rare for more than 77 to be needed.

Overflow

It is your responsibility to guard against potential overflows during evaluation of the matrix-vector products. In particular, when estimating B11B-11 using a triangular factorization of BB, nag_linsys_real_norm_rcomm (f04yc) should not be called if one of the factors is exactly singular – otherwise division by zero may occur in the substitutions.

Use in Conjunction with NAG Library Routines

To estimate the 11-norm of the inverse of a matrix AA, the following skeleton code can normally be used:
...  code to factorize A ...  
if (A is not singular)
  icase = 0;
  [icase, x, estnrm, work, iwork, ifail] = f04yc(icase, x, estnrm, work, iwork);
  while (icase ~= 0)
    if (icase == 1)
       ...  code to compute inv(A)*x ...  
    else 
       ...  code to compute inv(transpose(A))*x ...  
    end 
    [icase, x, estnrm, work, iwork, ifail] = f04yc(icase, x, estnrm, work, iwork);
  end 
end
To compute A1xA-1x or ATxA-Tx, solve the equation Ay = xAy=x or ATy = xATy=x for yy, overwriting yy on xx. The code will vary, depending on the type of the matrix AA, and the NAG function used to factorize AA.
Note that if AA is any type of symmetric matrix, then A = ATA=AT, and the ifstatement after the while can be reduced to:
       ...  code to compute inv(A)*x ...
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 [Example] illustrates how nag_linsys_real_norm_rcomm (f04yc) can be used in conjunction with NAG Toolbox functions for two important types of matrix: full nonsymmetric matrices (factorized by nag_lapack_dgetrf (f07ad)) and sparse nonsymmetric matrices (factorized by nag_matop_real_gen_sparse_lu (f01br)).
It is straightforward to use nag_linsys_real_norm_rcomm (f04yc) for the following other types of matrix, using the named functions for factorization and solution:

Example

function nag_linsys_real_norm_rcomm_example
a = [ 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];
icase = int64(0);
x = zeros(5, 1);
estnrm = 0;
work = zeros(5, 1);
iwork = zeros(5, 1, 'int64');
[icase, x, estnrm, work, iwork, ifail] = ...
    nag_linsys_real_norm_rcomm(icase, x, estnrm, work, iwork);
while (icase > 0)
  if (icase ==1)
    x = a*x;
  elseif (icase == 2)
    x = transpose(a)*x;
  end
  [icase, x, estnrm, work, iwork, ifail] = ...
    nag_linsys_real_norm_rcomm(icase, x, estnrm, work, iwork);
end
 estnrm
 

estnrm =

   15.9000


function f04yc_example
a = [ 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];
icase = int64(0);
x = zeros(5, 1);
estnrm = 0;
work = zeros(5, 1);
iwork = zeros(5, 1, 'int64');
[icase, x, estnrm, work, iwork, ifail] = f04yc(icase, x, estnrm, work, iwork);
while (icase > 0)
  if (icase ==1)
    x = a*x;
  elseif (icase == 2)
    x = transpose(a)*x;
  end
  [icase, x, estnrm, work, iwork, ifail] = f04yc(icase, x, estnrm, work, iwork);
end
 estnrm
 

estnrm =

   15.9000



PDF version (NAG web site, 64-bit version, 64-bit version)
Chapter Contents
Chapter Introduction
NAG Toolbox

© The Numerical Algorithms Group Ltd, Oxford, UK. 2009–2013