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_gen_norm_rcomm (f04yd)

Purpose

nag_linsys_real_gen_norm_rcomm (f04yd) estimates the 11-norm of a real 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.

Syntax

[irevcm, x, y, estnrm, work, iwork, ifail] = f04yd(irevcm, x, y, estnrm, seed, work, iwork, 'm', m, 'n', n, 't', t)
[irevcm, x, y, estnrm, work, iwork, ifail] = nag_linsys_real_gen_norm_rcomm(irevcm, x, y, estnrm, seed, work, iwork, 'm', m, 'n', n, 't', t)

Description

nag_linsys_real_gen_norm_rcomm (f04yd) computes an estimate (a lower bound) for the 11-norm
m
A1 = max |aij|
1jn i = 1
A1 = max 1jn i=1 m |aij|
(1)
of an mm 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 n × tn×t matrix XX (with tntn) or an m × tm×t matrix YY, can return AXAX or ATYATY. 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 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 for a square matrix, BB, 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 products B1XB-1X and BTYB-TY required by nag_linsys_real_gen_norm_rcomm (f04yd) 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_gen_norm_rcomm (f04yd) can be used to estimate the -norm of AA by working with ATAT instead of AA.
The algorithm used is described in Higham and Tisseur (2000).

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 11-norm estimation, with an application to 11-norm pseudospectra SIAM J. Matrix. Anal. Appl. 21 1185–1201

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

Compulsory Input Parameters

1:     irevcm – int64int32nag_int scalar
On initial entry: must be set to 00.
On intermediate re-entry: irevcm must be unchanged.
2:     x(ldx, : :) – double array
The first dimension of the array x must be at least nn
The second dimension of the array must be at least max (1,t)max(1,t)
On initial entry: need not be set.
On intermediate re-entry: if irevcm = 2irevcm=2, must contain ATYATY.
3:     y(ldy, : :) – double array
The first dimension of the array y must be at least mm
The second dimension of the array must be at least max (1,t)max(1,t)
On initial entry: need not be set.
On intermediate re-entry: if irevcm = 1irevcm=1, must contain AXAX.
4:     estnrm – double scalar
On initial entry: need not be set.
On intermediate re-entry: must not be changed.
5:     seed – int64int32nag_int scalar
The seed used for random number generation.
If t = 1t=1, seed is not used.
6:     work(m × tm×t) – double array
7:     iwork(2 × n + 5 × t + 202×n+5×t+20) – int64int32nag_int array
On initial entry: need not be set.
On intermediate re-entry: must not be changed.

Optional Input Parameters

1:     m – int64int32nag_int scalar
Default: The first dimension of the arrays y, work. (An error is raised if these dimensions are not equal.)
The number of rows of the matrix AA.
Constraint: m0m0.
2:     n – int64int32nag_int scalar
Default: The first dimension of the array x and the dimension of the array iwork. (An error is raised if these dimensions are not equal.)
nn, the number of columns of the matrix AA.
Constraint: n0n0.
3:     t – int64int32nag_int scalar
Default: The second dimension of the arrays x, y.
The number of columns tt of the matrices XX and YY. This is a parameter that can be used to control the accuracy and reliability of the estimate and corresponds roughly to the number of columns of AA that are visited during each iteration of the algorithm.
If t2t2 then a partly random starting matrix is used in the algorithm.
Default: t = 2t=2
Constraint: 1tm1tm.

Input Parameters Omitted from the MATLAB Interface

ldx ldy

Output Parameters

1:     irevcm – int64int32nag_int scalar
On intermediate exit: irevcm = 1irevcm=1 or 22, and x and y contain the elements of n × tn×t and m × tm×t matrices respectively. The calling program must
(a) if irevcm = 1irevcm=1, evaluate AXAX and store the result in y
or
if irevcm = 2irevcm=2, evaluate ATYATY and store the result in x,
(b) call nag_linsys_real_gen_norm_rcomm (f04yd) once again, with all the other parameters unchanged.
On final exit: irevcm = 0irevcm=0.
2:     x(ldx, : :) – double array
The first dimension of the array x will be nn
The second dimension of the array will be max (1,t)max(1,t)
ldxnldxn.
On intermediate exit: if irevcm = 1irevcm=1, contains the current matrix XX.
ldxnldxn.
On final exit: the array is undefined.
3:     y(ldy, : :) – double array
The first dimension of the array y will be mm
The second dimension of the array will be max (1,t)max(1,t)
ldymldym.
On intermediate exit: if irevcm = 2irevcm=2, contains the current matrix YY.
ldymldym.
On final exit: the array is undefined.
4:     estnrm – double scalar
On final exit: an estimate (a lower bound) for A1A1.
5:     work(m × tm×t) – double array
6:     iwork(2 × n + 5 × t + 202×n+5×t+20) – int64int32nag_int array
7:     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
Internal error; please contact NAG.
  ifail = 1ifail=-1
Constraint: irevcm = 0irevcm=0, 11 or 22.
On initial entry.
Constraint: irevcm = 0irevcm=0.
  ifail = 2ifail=-2
Constraint: m0m0.
  ifail = 3ifail=-3
Constraint: n0n0.
  ifail = 5ifail=-5
Constraint: ldxnldxn.
  ifail = 7ifail=-7
Constraint: ldymldym.
  ifail = 9ifail=-9
Constraint: 1tm1tm.

Accuracy

In extensive tests on random matrices of size up to m = n = 450m=n=450 the estimate estnrm has been found always to be within a factor two 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 and Tisseur (2000) for further details.

Further Comments

Timing

For most problems the time taken during calls to nag_linsys_real_gen_norm_rcomm (f04yd) will be negligible compared with the time spent evaluating matrix products between calls to nag_linsys_real_gen_norm_rcomm (f04yd).
The number of matrix products required depends on the matrix AA. At most six products of the form Y = AXY=AX and five products of the form X = ATYX=ATY will be required. The number of iterations is independent of the choice of tt.

Overflow

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

Choice of t

The parameter tt controls the accuracy and reliability of the estimate. For t = 1t=1, the algorithm behaves similarly to the LAPACK estimator xLACON. Increasing tt typically improves the estimate, without increasing the number of iterations required.
For t2t2, random matrices are used in the algorithm, so for repeatable results the value of seed should be kept constant.
A value of t = 2t=2 is recommended for new users.

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] = f04yd(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] = f04yd(icase, x, estnrm, work, iwork);
  end 
end
To compute A1XA-1X or ATYA-TY, solve the equation AY = XAY=X or ATX = YATX=Y, storing the result in y or x respectively. The code will vary, depending on the type of the matrix AA, and the NAG function used to factorize AA.
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_gen_norm_rcomm (f04yd) can be used in conjunction with NAG Toolbox functions for LULU factorization of a real matrix nag_lapack_dgetrf (f07ad).
It is straightforward to use nag_linsys_real_gen_norm_rcomm (f04yd) for the following other types of matrix, using the named functions for factorization and solution:

Example

function nag_linsys_real_gen_norm_rcomm_example

% Example 1: Compute the condition number of a dense matrix

a = [0.7,  -0.2,   1.0,   0.0,   2.0,   0.1;
     0.3,   0.7,   0.0,   1.0,   0.9,   0.2;
     0.0,   0.0,   0.2,   0.7,   0.0,  -1.1;
     0.0,   3.4,  -0.7,   0.2,   0.1,   0.1;
     0.0,  -4.0,   0.0,   1.0,   9.0,   0.0;
     0.4,   1.2,   4.3,   0.0,   6.2,   5.9];
t = int64(2);
m = 6;
n = 6;
x = zeros(n, t);
y = zeros(m, t);
estnrm = 0;
seed = int64(354);
irevcm = int64(0);
work = zeros(n*t, 1);
iwork = zeros(2*n+5*t+20, 1, 'int64');

nrma =  norm(a, 1);
fprintf('\nExample 1\n');
fprintf('\nThe norm of a is %6.2f\n', nrma);

% Estimate the norm of a^(-1) without explicitly forming a^(-1)

% Perform an LU factorization so that A=LU where L and U are lower and upper
% triangular.
[a, ipiv, info] = nag_lapack_dgetrf(a);

first = true;

while first || (irevcm ~= 0)
  first = false;

  [irevcm, x, y, estnrm, work, iwork, ifail] = ...
      nag_linsys_real_gen_norm_rcomm(irevcm, x, y, estnrm, seed, work, iwork);

  switch irevcm
    case 1
      % Compute y = inv(a)*x
      [y, info] = nag_lapack_dgetrs('n', a, ipiv, x);
    case 2
      % Compute x = transpose(inv(a))*y
      [x, info] = nag_lapack_dgetrs('t', a, ipiv, y);
    otherwise
  end
end

fprintf('The estimated norm of inverse(a) is: %6.2f\n', estnrm);
fprintf('\nEstimated condition number of a: %6.2f\n', estnrm*nrma);

% Example 2: Compute the condition number of a sparse matrix
% (stored in symmetric coordinate storage format)

t = int64(2);
n = int64(5);
nz = int64(10);
a = zeros(4*nz, 1);
icn = zeros(4*nz, 1, 'int64');
irn = zeros(4*nz, 1, 'int64');
a(1:nz) = [3; 2; 1; 2; 1; 4; 2; 1; 2; 5];
irn(1:nz) = [2; 4; 2; 3; 5; 4; 5; 1; 3; 4];
icn(1:nz) = [1; 1; 2; 2; 2; 3; 3; 4; 4; 5];

x = zeros(n, t);
y = zeros(n, t);
estnrm = 0;
seed = int64(412);
irevcm = int64(0);
work = zeros(n*t, 1);
iwork = zeros(2*n+5*t+20, 1, 'int64');

% Compute 1-norm of a
nrma =  0;
for i = 1:n
  asum = 0;
  for j = 1:nz
    if (icn(j)==i)
      asum = asum + abs(a(j));
    end
  end
  nrma = max(nrma,asum);
end
fprintf('\nExample 2\n');
fprintf('\nThe norm of a is %6.2f\n', nrma);

% Perform an LU factorization so that A=LU where L and U are lower and upper
% triangular.
abort = [true; true; false; false];
[a, irn, icn, ikeep, w, idisp, ifail] = nag_matop_real_gen_sparse_lu(n, nz, a, irn, icn, abort);


% Compute an estimate of the 1-norm of inv(a)

first = true;

while first || (irevcm ~= 0)
  first = false;

  [irevcm, x, y, estnrm, work, iwork, ifail] = ...
      nag_linsys_real_gen_norm_rcomm(irevcm, x, y, estnrm, seed, work, iwork);

  switch irevcm
    case 1
      % Compute y = inv(a)*x
      for i=1:2
        [y(:, i), resid] = nag_linsys_real_sparse_fac_solve(a, icn, ikeep, x(:, i), irevcm, idisp);
      end
    case 2
      % Compute x = transpose(inv(a))*y
      for i=1:2
        [x(:, i), resid] = nag_linsys_real_sparse_fac_solve(a, icn, ikeep, y(:, i), irevcm, idisp);
      end
    otherwise
  end
end

fprintf('The estimated norm of inverse(a) is: %6.2f\n', estnrm);
fprintf('\nEstimated condition number of a: %6.2f\n', estnrm*nrma);
 

Example 1

The norm of a is  18.20
The estimated norm of inverse(a) is:   2.97

Estimated condition number of a:  54.14

Example 2

The norm of a is   6.00
The estimated norm of inverse(a) is:   3.37

Estimated condition number of a:  20.20

function f04yd_example

% Example 1: Compute the condition number of a dense matrix

a = [0.7,  -0.2,   1.0,   0.0,   2.0,   0.1;
     0.3,   0.7,   0.0,   1.0,   0.9,   0.2;
     0.0,   0.0,   0.2,   0.7,   0.0,  -1.1;
     0.0,   3.4,  -0.7,   0.2,   0.1,   0.1;
     0.0,  -4.0,   0.0,   1.0,   9.0,   0.0;
     0.4,   1.2,   4.3,   0.0,   6.2,   5.9];
t = int64(2);
m = 6;
n = 6;
x = zeros(n, t);
y = zeros(m, t);
estnrm = 0;
seed = int64(354);
irevcm = int64(0);
work = zeros(n*t, 1);
iwork = zeros(2*n+5*t+20, 1, 'int64');

nrma =  norm(a, 1);
fprintf('\nExample 1\n');
fprintf('\nThe norm of a is %6.2f\n', nrma);

% Estimate the norm of a^(-1) without explicitly forming a^(-1)

% Perform an LU factorization so that A=LU where L and U are lower and upper
% triangular.
[a, ipiv, info] = f07ad(a);

first = true;

while first || (irevcm ~= 0)
  first = false;

  [irevcm, x, y, estnrm, work, iwork, ifail] = ...
      f04yd(irevcm, x, y, estnrm, seed, work, iwork);

  switch irevcm
    case 1
      % Compute y = inv(a)*x
      [y, info] = f07ae('n', a, ipiv, x);
    case 2
      % Compute x = transpose(inv(a))*y
      [x, info] = f07ae('t', a, ipiv, y);
    otherwise
  end
end

fprintf('The estimated norm of inverse(a) is: %6.2f\n', estnrm);
fprintf('\nEstimated condition number of a: %6.2f\n', estnrm*nrma);

% Example 2: Compute the condition number of a sparse matrix
% (stored in symmetric coordinate storage format)

t = int64(2);
n = int64(5);
nz = int64(10);
a = zeros(4*nz, 1);
icn = zeros(4*nz, 1, 'int64');
irn = zeros(4*nz, 1, 'int64');
a(1:nz) = [3; 2; 1; 2; 1; 4; 2; 1; 2; 5];
irn(1:nz) = [2; 4; 2; 3; 5; 4; 5; 1; 3; 4];
icn(1:nz) = [1; 1; 2; 2; 2; 3; 3; 4; 4; 5];

x = zeros(n, t);
y = zeros(n, t);
estnrm = 0;
seed = int64(412);
irevcm = int64(0);
work = zeros(n*t, 1);
iwork = zeros(2*n+5*t+20, 1, 'int64');

% Compute 1-norm of a
nrma =  0;
for i = 1:n
  asum = 0;
  for j = 1:nz
    if (icn(j)==i)
      asum = asum + abs(a(j));
    end
  end
  nrma = max(nrma,asum);
end
fprintf('\nExample 2\n');
fprintf('\nThe norm of a is %6.2f\n', nrma);

% Perform an LU factorization so that A=LU where L and U are lower and upper
% triangular.
abort = [true; true; false; false];
[a, irn, icn, ikeep, w, idisp, ifail] = f01br(n, nz, a, irn, icn, abort);


% Compute an estimate of the 1-norm of inv(a)

first = true;

while first || (irevcm ~= 0)
  first = false;

  [irevcm, x, y, estnrm, work, iwork, ifail] = ...
      f04yd(irevcm, x, y, estnrm, seed, work, iwork);

  switch irevcm
    case 1
      % Compute y = inv(a)*x
      for i=1:2
        [y(:, i), resid] = f04ax(a, icn, ikeep, x(:, i), irevcm, idisp);
      end
    case 2
      % Compute x = transpose(inv(a))*y
      for i=1:2
        [x(:, i), resid] = f04ax(a, icn, ikeep, y(:, i), irevcm, idisp);
      end
    otherwise
  end
end

fprintf('The estimated norm of inverse(a) is: %6.2f\n', estnrm);
fprintf('\nEstimated condition number of a: %6.2f\n', estnrm*nrma);
 

Example 1

The norm of a is  18.20
The estimated norm of inverse(a) is:   2.97

Estimated condition number of a:  54.14

Example 2

The norm of a is   6.00
The estimated norm of inverse(a) is:   3.37

Estimated condition number of a:  20.20


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