PDF version (NAG web site
, 64bit version, 64bit version)
NAG Toolbox: nag_linsys_withdraw_complex_norm_rcomm (f04zc)
Purpose
nag_linsys_complex_norm_rcomm (f04zc) estimates the $1$norm of a complex matrix without accessing the matrix explicitly. It uses reverse communication for evaluating matrixvector products. The function may be used for estimating matrix condition numbers.
Note: this function is scheduled to be withdrawn, please see
f04zc in
Advice on Replacement Calls for Withdrawn/Superseded Routines..
Syntax
Description
nag_linsys_complex_norm_rcomm (f04zc) computes an estimate (a lower bound) for the
$1$norm
of an
$n$ by
$n$ complex matrix
$A=\left({a}_{ij}\right)$. The function regards the matrix
$A$ as being defined by a usersupplied ‘Black Box’ which, given an input vector
$x$, can return either of the matrixvector products
$Ax$ or
${A}^{\mathrm{H}}x$, where
${A}^{\mathrm{H}}$ is the complex conjugate transpose. A reverse communication interface is used; thus control is returned to the calling program whenever a matrixvector 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 ${\Vert {B}^{1}\Vert}_{1}$, and hence the condition number
${\kappa}_{1}\left(B\right)={\Vert B\Vert}_{1}{\Vert {B}^{1}\Vert}_{1}$, without forming ${B}^{1}$ explicitly ($A={B}^{1}$ above).
If, for example, an $LU$ factorization of $B$ is available, the matrixvector products ${B}^{1}x$ and ${B}^{\mathrm{H}}x$ required by nag_linsys_complex_norm_rcomm (f04zc) may be computed by back and forwardsubstitutions, without computing ${B}^{1}$.
The function can also be used to estimate
$1$norms of matrix products such as
${A}^{1}B$ and
$ABC$, without forming the products explicitly. Further applications are described in
Higham (1988).
Since ${\Vert A\Vert}_{\infty}={\Vert {A}^{\mathrm{H}}\Vert}_{1}$, nag_linsys_complex_norm_rcomm (f04zc) can be used to estimate the $\infty $norm of $A$ by working with ${A}^{\mathrm{H}}$ instead of $A$.
The algorithm used is based on a method given in
Hager (1984) and is described in
Higham (1988). A comparison of several techniques for condition number estimation is given in
Higham (1987).
Note: nag_linsys_complex_gen_norm_rcomm (f04zd) can also be used to estimate the norm of a real matrix.
nag_linsys_complex_gen_norm_rcomm (f04zd) uses a more recent algorithm than
nag_linsys_complex_norm_rcomm (f04zc) and it is recommended that
nag_linsys_complex_gen_norm_rcomm (f04zd) be used in place of
nag_linsys_complex_norm_rcomm (f04zc).
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 onenorm 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 reentries, and a final exit, as indicated by the
argument icase. Between intermediate exits and reentries,
all arguments other than x must remain unchanged.
Compulsory Input Parameters
 1:
$\mathrm{icase}$ – int64int32nag_int scalar

On initial entry: must be set to $0$.
 2:
$\mathrm{x}\left({\mathbf{n}}\right)$ – complex array

On initial entry: need not be set.
On intermediate reentry: must contain $Ax$ (if ${\mathbf{icase}}=1$) or ${A}^{\mathrm{H}}x$ (if ${\mathbf{icase}}=2$).
 3:
$\mathrm{estnrm}$ – double scalar

On initial entry: need not be set.
 4:
$\mathrm{work}\left({\mathbf{n}}\right)$ – complex array

On initial entry: need not be set.
Optional Input Parameters
 1:
$\mathrm{n}$ – int64int32nag_int scalar

Default:
the dimension of the arrays
x,
work. (An error is raised if these dimensions are not equal.)
On initial entry: $n$, the order of the matrix $A$.
Constraint:
${\mathbf{n}}\ge 1$.
Output Parameters
 1:
$\mathrm{icase}$ – int64int32nag_int scalar

On intermediate exit:
${\mathbf{icase}}=1$ or
$2$, and
${\mathbf{x}}\left(\mathit{i}\right)$, for
$\mathit{i}=1,2,\dots ,n$, contain the elements of a vector
$x$. The calling program must
(a) 
evaluate $Ax$ (if ${\mathbf{icase}}=1$) or ${A}^{\mathrm{H}}x$ (if ${\mathbf{icase}}=2$), where ${A}^{\mathrm{H}}$ is the complex conjugate transpose; 
(b) 
place the result in x; and, 
(c) 
call nag_linsys_complex_norm_rcomm (f04zc) once again, with all the other arguments unchanged. 
On final exit: ${\mathbf{icase}}=0$.
 2:
$\mathrm{x}\left({\mathbf{n}}\right)$ – complex array

On intermediate exit:
contains the current vector $x$.
On final exit: the array is undefined.
 3:
$\mathrm{estnrm}$ – double scalar

On intermediate exit:
should not be changed.
On final exit: an estimate (a lower bound) for ${\Vert A\Vert}_{1}$.
 4:
$\mathrm{work}\left({\mathbf{n}}\right)$ – complex array

On final exit: contains a vector
$v$ such that
$v=Aw$ where
${\mathbf{estnrm}}={\Vert v\Vert}_{1}/{\Vert w\Vert}_{1}$ (
$w$ is not returned). If
$A={B}^{1}$ and
estnrm is large, then
$v$ is an approximate null vector for
$B$.
 5:
$\mathrm{ifail}$ – int64int32nag_int scalar
${\mathbf{ifail}}={\mathbf{0}}$ unless the function detects an error (see
Error Indicators and Warnings).
Error Indicators and Warnings
Errors or warnings detected by the function:
 ${\mathbf{ifail}}=1$

On entry,  ${\mathbf{n}}<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
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
${\Vert A\Vert}_{1}$; often the estimate has many correct figures. However, matrices exist for which the estimate is smaller than
${\Vert A\Vert}_{1}$ 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 by nag_linsys_complex_norm_rcomm (f04zc) is proportional to $n$. For most problems the time taken during calls to nag_linsys_complex_norm_rcomm (f04zc) will be negligible compared with the time spent evaluating matrixvector products between calls to nag_linsys_complex_norm_rcomm (f04zc).
The number of matrixvector products required varies from $5$ to $11$ (or is $1$ if $n=1$). In most cases $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 matrixvector products. In particular, when estimating ${\Vert {B}^{1}\Vert}_{1}$ using a triangular factorization of $B$, nag_linsys_complex_norm_rcomm (f04zc) 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, ifail] = f04zc(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] = f04zc(icase, x, estnrm, work);
end
end
To compute ${A}^{1}x$ or ${A}^{\mathrm{H}}x$, solve the equation $Ay=x$ or ${A}^{\mathrm{H}}y=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
Hermitian matrix, then
$A={A}^{\mathrm{H}}$, and the
if statement after the
while
can be reduced to:
... code to compute A(1)x ...
The example program in
Example illustrates how
nag_linsys_complex_norm_rcomm (f04zc) can be used in conjunction with NAG Toolbox functions for complex band matrices (factorized by
nag_lapack_zgbtrf (f07br)).
It is also straightforward to use
nag_linsys_complex_norm_rcomm (f04zc) 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
${\Vert A\Vert}_{1}{\Vert {A}^{1}\Vert}_{1}$ of the order
$5$ matrix
where
$A$ is a band matrix stored in the packed format required by
nag_lapack_zgbtrf (f07br) and
nag_lapack_zgbtrs (f07bs).
Further examples of the technique for condition number estimation in the case of double matrices can be seen in the example program section of
nag_linsys_real_norm_rcomm (f04yc).
Open in the MATLAB editor:
f04zc_example
function f04zc_example
fprintf('f04zc example results\n\n');
a = [ 1.0 + 1.0i, 2.0 + 1.0i, 1.0 + 2.0i, 0.0 + 0.0i, 0.0 + 0.0i;
0.0 + 2.0i, 3.0 + 5.0i, 1.0 + 3.0i, 2.0 + 1.0i, 0.0 + 0.0i,
0.0 + 0.0i, 2.0 + 6.0i, 5.0 + 7.0i, 0.0 + 6.0i, 1.0  1.0i;
0.0 + 0.0i, 0.0 + 0.0i, 3.0 + 9.0i, 0.0 + 4.0i, 4.0  3.0i;
0.0 + 0.0i, 0.0 + 0.0i, 0.0 + 0.0i, 1.0 + 8.0i, 10.0  3.0i];
x = complex(zeros(5, 1));
work = complex(zeros(5,1));
anorm = norm(a,1);
icase = int64(0);
estnrm = 0;
done = false;
while (~done)
[icase, x, estnrm, work, ifail] = ...
f04zc(icase, x, estnrm, work);
if (icase == 0)
done = true;
elseif (icase == 1)
x = inv(a)*x;
else
x = conj(transpose(inv(a)))*x;
end
end
fprintf('Computed norm of a = %6.4g\n', anorm);
fprintf('Estimated norm of inverse(A) = %6.4g\n', estnrm);
fprintf('Estimated condition number of A = %6.1f\n', estnrm*anorm);
f04zc example results
Computed norm of a = 23.49
Estimated norm of inverse(A) = 37.04
Estimated condition number of A = 870.0
PDF version (NAG web site
, 64bit version, 64bit version)
© The Numerical Algorithms Group Ltd, Oxford, UK. 2009–2015