NAG Toolbox: nag_linsys_real_gen_norm_rcomm (f04yd)

Purpose

Syntax

Description

References

Parameters

Compulsory Input Parameters

1:     $\mathrm{irevcm}$ – int64int32nag_int scalar

On initial entry: must be set to $0$.

On intermediate re-entry: irevcm must be unchanged.

2:     $\mathrm{x}\left(\mathit{ldx},:\right)$ – double array

The first dimension of the array x must be at least $\mathbf{n}$.

The second dimension of the array x must be at least $\max \left(1,\mathbf{t}\right)$.

On initial entry: need not be set.

On intermediate re-entry: if $\mathbf{irevcm}=2$, must contain $A^{\mathrm{T}}Y$.

3:     $\mathrm{y}\left(\mathit{ldy},:\right)$ – double array

The first dimension of the array y must be at least $\mathbf{m}$.

The second dimension of the array y must be at least $\max \left(1,\mathbf{t}\right)$.

On initial entry: need not be set.

On intermediate re-entry: if $\mathbf{irevcm}=1$, must contain $AX$.

4:     $\mathrm{estnrm}$ – double scalar

On initial entry: need not be set.

On intermediate re-entry: must not be changed.

5:     $\mathrm{seed}$ – int64int32nag_int scalar

The seed used for random number generation.

If $\mathbf{t}=1$, seed is not used.

Constraint: if $\mathbf{t}>1$, $\mathbf{seed}\ge 1$.

6:     $\mathrm{work}\left(\mathbf{m}\times \mathbf{t}\right)$ – double array

7:     $\mathrm{iwork}\left(2\times \mathbf{n}+5\times \mathbf{t}+20\right)$ – int64int32nag_int array

On initial entry: need not be set.

On intermediate re-entry: must not be changed.

Optional Input Parameters

1:     $\mathrm{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: $\mathbf{m}\ge 0$.

2:     $\mathrm{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.)

$n$, the number of columns of the matrix $A$.

Constraint: $\mathbf{n}\ge 0$.

3:     $\mathrm{t}$ – int64int32nag_int scalar

Suggested value: $\mathbf{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 $\mathbf{t}\ge 2$ then a partly random starting matrix is used in the algorithm.

Constraint: $1\le \mathbf{t}\le \mathbf{m}$.

Output Parameters

1:     $\mathrm{irevcm}$ – int64int32nag_int scalar

On intermediate exit: $\mathbf{irevcm}=1$ or $2$, and x contains the $n\times t$ matrix $X$ and y contains the $m\times t$ matrix $Y$. The calling program must

(a)	if $\mathbf{irevcm}=1$, evaluate $AX$ and store the result in y or if $\mathbf{irevcm}=2$, evaluate $A^{\mathrm{T}}Y$ and store the result in x,
(b)	call nag_linsys_real_gen_norm_rcomm (f04yd) once again, with all the other arguments unchanged.

On final exit: $\mathbf{irevcm}=0$.

2:     $\mathrm{x}\left(\mathit{ldx},:\right)$ – double array

The first dimension of the array x will be $\mathbf{n}$.

The second dimension of the array x will be $\max \left(1,\mathbf{t}\right)$.

On intermediate exit: if $\mathbf{irevcm}=1$, contains the current matrix $X$.

On final exit: the array is undefined.

3:     $\mathrm{y}\left(\mathit{ldy},:\right)$ – double array

The first dimension of the array y will be $\mathbf{m}$.

The second dimension of the array y will be $\max \left(1,\mathbf{t}\right)$.

On intermediate exit: if $\mathbf{irevcm}=2$, contains the current matrix $Y$.

On final exit: the array is undefined.

4:     $\mathrm{estnrm}$ – double scalar

On final exit: an estimate (a lower bound) for ${‖A‖}_1$.

5:     $\mathrm{work}\left(\mathbf{m}\times \mathbf{t}\right)$ – double array

6:     $\mathrm{iwork}\left(2\times \mathbf{n}+5\times \mathbf{t}+20\right)$ – int64int32nag_int array

7:     $\mathrm{ifail}$ – int64int32nag_int scalar

On final exit: $\mathbf{ifail}=\mathbf{0}$ unless the function detects an error (see Error Indicators and Warnings).

Error Indicators and Warnings

   $\mathbf{ifail}=1$

Internal error; please contact NAG.

   $\mathbf{ifail}=-1$

Constraint: $\mathbf{irevcm}=0$, $1$ or $2$.

On initial entry, $\mathbf{irevcm}=\_$.
Constraint: $\mathbf{irevcm}=0$.

   $\mathbf{ifail}=-2$

Constraint: $\mathbf{m}\ge 0$.

   $\mathbf{ifail}=-3$

Constraint: $\mathbf{n}\ge 0$.

   $\mathbf{ifail}=-5$

Constraint: $\mathit{ldx}\ge \mathbf{n}$.

   $\mathbf{ifail}=-7$

Constraint: $\mathit{ldy}\ge \mathbf{m}$.

   $\mathbf{ifail}=-9$

Constraint: $1\le \mathbf{t}\le \mathbf{m}$.

   $\mathbf{ifail}=-10$

Constraint: if $\mathbf{t}>1$, $\mathbf{seed}\ge 1$.

   $\mathbf{ifail}=-99$

An unexpected error has been triggered by this routine. Please contact NAG.

   $\mathbf{ifail}=-399$

Your licence key may have expired or may not have been installed correctly.

   $\mathbf{ifail}=-999$

Dynamic memory allocation failed.

Accuracy

Further Comments

Timing

Overflow

Choice of t

Use in Conjunction with NAG Library Routines

...  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

Example

Open in the MATLAB editor: f04yd_example

function f04yd_example


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


% Example 1: Compute the condition number of a dense matrix
fprintf('\nExample 1\n');

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('\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)
fprintf('\nExample 2\n');

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('\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);

f04yd example results


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