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_withdraw_real_norm_rcomm (f04yc)

▸▿ Contents

    1  Purpose
    2  Syntax
    3  Description
    4  References
▸▿ 5  Parameters
    6  Error Indicators and Warnings
    7  Accuracy
▸▿ 8  Further Comments
    9  Example

Purpose

nag_linsys_real_norm_rcomm (f04yc) estimates the 1-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_withdraw_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 1-norm
A1 = max 1≤j≤n ∑i=1n aij (1)
of an n by n real matrix A=aij. The function 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 Ax or ATx. 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 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-11, and hence the condition number κ1B=B1B-11, without forming B-1 explicitly (A=B-1 above).
If, for example, an LU factorization of B is available, the matrix-vector products B-1x and B-Tx required by nag_linsys_real_norm_rcomm (f04yc) 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-1B and ABC, without forming the products explicitly. Further applications are described by Higham (1988).
Since A∞=AT1, nag_linsys_real_norm_rcomm (f04yc) can be used to estimate the ∞-norm of A by working with AT 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: 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 argument icase. Between intermediate exits and re-entries, all arguments other than x must remain unchanged.

Compulsory Input Parameters

1:     icase – int64int32nag_int scalar
On initial entry: must be set to 0.
2:     xn – double array
On initial entry: need not be set.
On intermediate re-entry: must contain Ax (if icase=1) or ATx (if icase=2).
3:     estnrm – double scalar
On initial entry: need not be set.
4:     workn – double array
On initial entry: need not be set.
5:     iworkn – 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: n, the order of the matrix A.
Constraint: n≥1.

Output Parameters

1:     icase – int64int32nag_int scalar
On intermediate exit: icase=1 or 2, and xi, for i=1,2,…,n, contain the elements of a vector x. The calling program must
(a) evaluate Ax (if icase=1) or ATx (if icase=2),
(b) place the result in x, and
(c) call nag_linsys_real_norm_rcomm (f04yc) once again, with all the other arguments unchanged.
On final exit: icase=0.
2:     xn – double array
On intermediate exit: contains the current vector x.
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 A1.
4:     workn – double array
On final exit: contains a vector v such that v=Aw where estnrm=v1/w1 (w is not returned). If A=B-1 and estnrm is large, then v is an approximate null vector for B.
5:     iworkn – int64int32nag_int array
6:     ifail – int64int32nag_int scalar
On final exit: ifail=0 unless the function detects an error (see Error Indicators and Warnings).

Error Indicators and Warnings

Errors or warnings detected by the function:
   ifail=1
On entry, n<1.
   ifail=-99
An unexpected error has been triggered by this routine. Please contact NAG.
   ifail=-399
Your licence key may have expired or may not have been installed correctly.
   ifail=-999
Dynamic memory allocation failed.

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 A1; often the estimate has many correct figures. However, matrices exist for which the estimate is smaller than A1 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 n. 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 4 to 11 (or is 1 if n=1). In most cases 4 or 5 products are required; it is rare for more than 7 to be needed.

Overflow

It is your responsibility to guard against potential overflows during evaluation of the matrix-vector products. In particular, when estimating B-11 using a triangular factorization of B, 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 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)
  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 A-1x or A-Tx, solve the equation Ay=x or ATy=x for y, overwriting y on x. The code will vary, depending on the type of the matrix A, and the NAG function used to factorize A.
Note that if A is any type of symmetric matrix, then A=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 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

For this function two examples are presented. There is a single example program for nag_linsys_real_norm_rcomm (f04yc), 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 A1A-11 of the matrix A given by
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 .  
Example 2 (EX2)
To estimate the condition number A1A-11 of the matrix A given by
A= 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 .  

Open in the MATLAB editor: f04yc_example

function f04yc_example


fprintf('f04yc example results\n\n');

ex1;
ex2;



function ex1
  fprintf('Example 1:\n\n');
  % Estimate condition number of A
  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];
  anorm = norm(a,1);
  ainv = inv(a);

  % reverse communication initializations
  icase = int64(0);
  estnrm = 0;
  x     = zeros(5, 1);
  work  = zeros(5, 1);
  iwork = zeros(5, 1, 'int64');

  % reverse communication loop
  done = false;
  while (~done)
    [icase, x, estnrm, work, iwork, ifail] = ...
    f04yc(icase, x, estnrm, work, iwork);
    if (icase == 0)
      done = true;
    elseif (icase == 1)
      x = ainv*x;
    elseif (icase == 2)
      x = ainv'*x;
    end
  end

  fprintf('Computed norm of A              = %8.4f\n', anorm);
  fprintf('Estimated norm of inverse(A)    = %8.4f\n', estnrm);
  cond = anorm*estnrm;
  fprintf('Estimated condition number of A = %8.1f\n', cond);

function ex2
  fprintf('\nExample 2:\n\n');
  % Estimate condition number of sparse A
  n  = int64(6);
  nz = int64(15);
  a   = zeros(150,1);
  irn = zeros(75,1,'int64');
  icn = zeros(150,1,'int64');
  a(1:nz) = [ 5                       ...
                  2  -1   2           ...
                      3               ...
             -2           1   1       ...
             -1          -1   2  -3   ...
             -1  -1               6];
  anorm = 9;
  irn(1:15) = [int64(1); 2;2;2; 3; 4;4;4; 5;5;5;5; 6;6;6];
  icn(1:15) = [int64(1); 2;3;4; 3; 1;4;5; 1;4;5;6; 1;2;6];

  % Factorize A
  abort     = [true;     true;     false;     true];
  [a, irn, icn, ikeep, w, idisp, ifail] = ...
  f01br( ...
         n, nz, a, irn, icn, abort);

  % reverse communication initializations
  icase = int64(0);
  estnrm = 0;
  x     = zeros(n, 1);
  work  = zeros(n, 1);
  iwork = zeros(n, 1, 'int64');
  
  % reverse communication loop
  done = false;
  while (~done)
    [icase, x, estnrm, work, iwork, ifail] = ...
    f04yc(icase, x, estnrm, work, iwork);
    if (icase == 0)
      done = true;
    else
      % Solve Ax = b or A'x = b
      [x, resid] = f04ax( ...
                          a, icn, ikeep, x, icase, idisp);
    end
  end
  
  fprintf('Computed norm of A              = %8.4f\n', anorm);
  fprintf('Estimated norm of inverse(A)    = %8.4f\n', estnrm);
  cond = anorm*estnrm;
  fprintf('Estimated condition number of A = %8.1f\n', cond);

f04yc example results

Example 1:

Computed norm of A              =  15.9000
Estimated norm of inverse(A)    =   1.7635
Estimated condition number of A =     28.0

Example 2:

Computed norm of A              =   9.0000
Estimated norm of inverse(A)    =   1.9333
Estimated condition number of A =     17.4
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–2015