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_complex_gen_norm_rcomm (f04zd)

▸▿ 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_complex_gen_norm_rcomm (f04zd) estimates the 1-norm of a complex 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, rwork, iwork, ifail] = f04zd(irevcm, x, y, estnrm, seed, work, rwork, iwork, 'm', m, 'n', n, 't', t)
[irevcm, x, y, estnrm, work, rwork, iwork, ifail] = nag_linsys_complex_gen_norm_rcomm(irevcm, x, y, estnrm, seed, work, rwork, iwork, 'm', m, 'n', n, 't', t)

Description

nag_linsys_complex_gen_norm_rcomm (f04zd) computes an estimate (a lower bound) for the 1-norm
A1 = max 1≤j≤n ∑ i=1 m aij (1)
of an m by n complex matrix A=aij. The function regards the matrix A as being defined by a user-supplied ‘Black Box’ which, given an n×t matrix X (with t≪n) or an m×t matrix Y, can return AX or AHY, where AH 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 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 the formula (1) above.
The main use of the function is for estimating B-11 for a square matrix B, 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 products B-1X and B-HY required by nag_linsys_complex_gen_norm_rcomm (f04zd) 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 in Higham (1988).
Since A∞=AH1, nag_linsys_complex_gen_norm_rcomm (f04zd) can be used to estimate the ∞-norm of A by working with AH instead of A.
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 1-norm estimation, with an application to 1-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 argument irevcm. Between intermediate exits and re-entries, all arguments other than x and y must remain unchanged.

Compulsory Input Parameters

1:     irevcm – int64int32nag_int scalar
On initial entry: must be set to 0.
On intermediate re-entry: irevcm must be unchanged.
2:     xldx: – complex array
The first dimension of the array x must be at least n.
The second dimension of the array x must be at least max1,t.
On initial entry: need not be set.
On intermediate re-entry: if irevcm=2, must contain AHY.
3:     yldy: – complex array
The first dimension of the array y must be at least m.
The second dimension of the array y must be at least max1,t.
On initial entry: need not be set.
On intermediate re-entry: if irevcm=1, must contain AX.
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=1, seed is not used.
Constraint: if t>1, seed≥1.
6:     workm×t – complex array
7:     rwork2×n – double array
8:     iwork2×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 A.
Constraint: m≥0.
2:     n – int64int32nag_int scalar
Default: the first dimension of the array x.
On initial entry: n, the number of columns of the matrix A.
Constraint: n≥0.
3:     t – int64int32nag_int scalar
Suggested value: t=2.
Default: the second dimension of the arrays x, y. (An error is raised if these dimensions are not equal.)
The number of columns t of the matrices X and Y. This is a 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.
If t≥2 then a partly random starting matrix is used in the algorithm.
Constraint: 1≤t≤m.

Output Parameters

1:     irevcm – int64int32nag_int scalar
On intermediate exit: irevcm=1 or 2, and x contains the n×t matrix X and y contains the m×t matrix Y. The calling program must
(a) if irevcm=1, evaluate AX and store the result in y
or
if irevcm=2, evaluate AHY and store the result in x, where AH is the complex conjugate transpose;
(b) call nag_linsys_complex_gen_norm_rcomm (f04zd) once again, with all the arguments unchanged.
On final exit: irevcm=0.
2:     xldx: – complex array
The first dimension of the array x will be n.
The second dimension of the array x will be max1,t.
On intermediate exit: if irevcm=1, contains the current matrix X.
On final exit: the array is undefined.
3:     yldy: – complex array
The first dimension of the array y will be m.
The second dimension of the array y will be max1,t.
On intermediate exit: if irevcm=2, contains the current matrix Y.
On final exit: the array is undefined.
4:     estnrm – double scalar
On final exit: an estimate (a lower bound) for A1.
5:     workm×t – complex array
6:     rwork2×n – double array
7:     iwork2×n+5×t+20 – int64int32nag_int array
8:     ifail – int64int32nag_int scalar
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
Internal error; please contact NAG.
   ifail=-1
Constraint: irevcm=0, 1 or 2.
On initial entry, irevcm=_.
Constraint: irevcm=0.
   ifail=-2
Constraint: m≥0.
   ifail=-3
Constraint: n≥0.
   ifail=-5
Constraint: ldx≥n.
   ifail=-7
Constraint: ldy≥m.
   ifail=-9
Constraint: 1≤t≤m.
   ifail=-10
Constraint: if t>1, seed≥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 m=n=450 the estimate estnrm has been found always to be within a factor two 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 and Tisseur (2000) for further details.

Further Comments

Timing

For most problems the time taken during calls to nag_linsys_complex_gen_norm_rcomm (f04zd) will be negligible compared with the time spent evaluating matrix products between calls to nag_linsys_complex_gen_norm_rcomm (f04zd).
The number of matrix products required depends on the matrix A. At most six products of the form Y=AX and five products of the form X=AHY will be required. The number of iterations is independent of the choice of t.

Overflow

It is your responsibility to guard against potential overflows during evaluation of the matrix products. In particular, when estimating B-11 using a triangular factorization of B, nag_linsys_complex_gen_norm_rcomm (f04zd) 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 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≥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.

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, ifail] = f04zd(icase, x, estnrm, work);
  while (icase ~= 0)
    if (icase == 1) 
      ...  code to compute A(-1)x ...
    else
      ...  code to compute (A(-1)(H)) x ...
    end
    [icase, x, estnrm, work, ifail] = f04zd(icase, x, estnrm, work);
  end
end
To compute A-1X or A-HY, solve the equation AY=X or AHX=Y storing the result in y or x respectively. The code will vary, depending on the type of the matrix A, and the NAG function used to factorize A.
The example program in Example illustrates how nag_linsys_complex_gen_norm_rcomm (f04zd) can be used in conjunction with NAG Toolbox function for LU factorization of complex matrices nag_lapack_zgetrf (f07ar)).
It is also straightforward to use nag_linsys_complex_gen_norm_rcomm (f04zd) for Hermitian positive definite matrices, using nag_lapack_zpotrf (f07fr) and nag_lapack_zpotrs (f07fs) for factorization and solution.

Example

This example estimates the condition number A1A-11 of the matrix A given by
A = 0.7+0.1i -0.2+0.0i 1.0+0.0i 0.0+0.0i 0.0+0.0i 0.1+0.0i 0.3+0.0i 0.7+0.0i 0.0+0.0i 1.0+0.2i 0.9+0.0i 0.2+0.0i 0.0+5.9i 0.0+0.0i 0.2+0.0i 0.7+0.0i 0.4+6.1i 1.1+0.4i 0.0+0.1i 0.0+0.1i -0.7+0.0i 0.2+0.0i 0.1+0.0i 0.1+0.0i 0.0+0.0i 4.0+0.0i 0.0+0.0i 1.0+0.0i 9.0+0.0i 0.0+0.1i 4.5+6.7i 0.1+0.4i 0.0+3.2i 1.2+0.0i 0.0+0.0i 7.8+0.2i .  

Open in the MATLAB editor: f04zd_example

function f04zd_example


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

a = [0.7+0.1i, -0.2+0.0i,  1.0+0.0i, 0.0+0.0i, 0.0+0.0i, 0.1+0.0i;
     0.3+0.0i,  0.7+0.0i,  0.0+0.0i, 1.0+0.2i, 0.9+0.0i, 0.2+0.0i;
     0.0+5.9i,  0.0+0.0i,  0.2+0.0i, 0.7+0.0i, 0.4+6.1i, 1.1+0.4i;
     0.0+0.1i,  0.0+0.1i, -0.7+0.0i, 0.2+0.0i, 0.1+0.0i, 0.1+0.0i;
     0.0+0.0i,  4.0+0.0i,  0.0+0.0i, 1.0+0.0i, 9.0+0.0i, 0.0+0.1i;
     4.5+6.7i,  0.1+0.4i,  0.0+3.2i, 1.2+0.0i, 0.0+0.0i, 7.8+0.2i];

t = int64(2);
m = 6;
n = 6;
x = complex(zeros(n, t));
y = complex(zeros(m, t));
estnrm = 0;
seed = int64(652);
work  = complex(zeros(m*t, 1));
rwork = zeros(2*n, 1);
iwork = zeros(2*n+5*t+20, 1, 'int64');

nrma =  norm(a, 1);

% 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] = f07ar(a);

done = false;
irevcm = int64(0);

while ~done

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

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

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



f04zd example results

The norm of a                    =  16.11
The estimated norm of inverse(a) =  24.02
Estimated condition number of a  = 387.08
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