Integer type:** int32**** int64**** nag_int** show int32 show int32 show int64 show int64 show nag_int show nag_int

nag_linsys_complex_norm_rcomm (f04zc) estimates the 1$1$-norm of a complex 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.

nag_linsys_complex_norm_rcomm (f04zc) computes an estimate (a lower bound) for the 1$1$-norm

of an n$n$ by n$n$ complex matrix A = (a_{ij})$A=\left({a}_{ij}\right)$. The function regards the matrix A$A$ as being defined by a user-supplied ‘Black Box’ which, given an input vector x$x$, can return either of the matrix-vector products Ax$Ax$ or A^{H}x${A}^{\mathrm{H}}x$, where A^{H}${A}^{\mathrm{H}}$ is the complex conjugate transpose. A reverse communication interface is used; thus control is returned to the calling program whenever a matrix-vector product is required.

$${\Vert A\Vert}_{1}=\underset{1\le j\le n}{\mathrm{max}}\phantom{\rule{0.25em}{0ex}}\sum _{i=1}^{n}\left|{a}_{ij}\right|$$ | (1) |

The **main**
**use** of the function is for estimating ‖B^{ − 1}‖_{1}${\Vert {B}^{-1}\Vert}_{1}$, and hence the **condition number**
κ_{1}(B) = ‖B‖_{1}‖B^{ − 1}‖_{1}${\kappa}_{1}\left(B\right)={\Vert B\Vert}_{1}{\Vert {B}^{-1}\Vert}_{1}$, without forming B^{ − 1}${B}^{-1}$ explicitly (A = B^{ − 1}$A={B}^{-1}$ above).

If, for example, an LU$LU$ factorization of B$B$ is available, the matrix-vector products B^{ − 1}x${B}^{-1}x$ and B^{ − H}x${B}^{-\mathrm{H}}x$ required by nag_linsys_complex_norm_rcomm (f04zc) may be computed by back- and forward-substitutions, without computing B^{ − 1}${B}^{-1}$.

The function can also be used to estimate 1$1$-norms of matrix products such as A^{ − 1}B${A}^{-1}B$ and ABC$ABC$, without forming the products explicitly. Further applications are described in Higham (1988).

Since ‖A‖_{∞} = ‖A^{H}‖_{1}${\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$A$ by working with A^{H}${A}^{\mathrm{H}}$ instead of A$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).

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

- 1: icase – int64int32nag_int scalar
*On initial entry*: must be set to 0$0$.- 2: x(n) – complex array
*On initial entry*: need not be set.- 3: estnrm – double scalar
*On initial entry*: need not be set.- 4: work(n) – complex array
*On initial entry*: need not be set.

- 1: 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$n$, the order of the matrix A$A$.

None.

- 1: icase – int64int32nag_int scalar
*On intermediate exit*: icase = 1${\mathbf{icase}}=1$ or 2$2$, and x(i)${\mathbf{x}}\left(\mathit{i}\right)$, for i = 1,2, … ,n$\mathit{i}=1,2,\dots ,n$, contain the elements of a vector x$x$. The calling program must(a) evaluate Ax$Ax$ (if icase = 1${\mathbf{icase}}=1$) or A ^{H}x${A}^{\mathrm{H}}x$ (if icase = 2${\mathbf{icase}}=2$), where A^{H}${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 parameters unchanged. - 2: x(n) – complex array
*On intermediate exit*: contains the current vector x$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 ‖A‖_{1}${\Vert A\Vert}_{1}$.- 4: work(n) – complex array
- 5: ifail – int64int32nag_int scalar
- ifail = 0${\mathrm{ifail}}={\mathbf{0}}$ unless the function detects an error (see [Error Indicators and Warnings]).

Errors or warnings detected by the function:

On entry, n < 1${\mathbf{n}}<1$.

In extensive tests on **random** matrices of size up to n = 100$n=100$ the estimate estnrm has been found always to be within a factor eleven of ‖A‖_{1}${\Vert A\Vert}_{1}$; often the estimate has many correct figures. However, matrices exist for which the estimate is smaller than ‖A‖_{1}${\Vert A\Vert}_{1}$ by an arbitrary factor; such matrices are very unlikely to arise in practice. See Higham (1988) for further details.

The total time taken by nag_linsys_complex_norm_rcomm (f04zc) is proportional to n$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 matrix-vector products between calls to nag_linsys_complex_norm_rcomm (f04zc).

The number of matrix-vector products required varies from 5$5$ to 11$11$ (or is 1$1$ if n = 1$n=1$). In most cases 5$5$ products are required; it is rare for more than 7$7$ to be needed.

It is your responsibility to guard against potential overflows during evaluation of the matrix-vector products. In particular, when estimating ‖B^{ − 1}‖_{1}${\Vert {B}^{-1}\Vert}_{1}$ using a triangular factorization of B$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.

To estimate the 1$1$-norm of the inverse of a matrix A$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${A}^{-1}x$ or A^{ − H}x${A}^{-\mathrm{H}}x$, solve the equation Ay = x$Ay=x$ or A^{H}y = x${A}^{\mathrm{H}}y=x$ for y$y$, overwriting y$y$ on x$x$. The code will vary, depending on the type of the matrix A$A$, and the NAG function used to factorize A$A$.

Note that if A$A$ is any type of **Hermitian** matrix, then A = A^{H}$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 Section [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.

Open in the MATLAB editor: nag_linsys_complex_norm_rcomm_example

function nag_linsys_complex_norm_rcomm_examplea = [ 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]; icase = int64(0); x = complex(zeros(5, 1)); estnrm = 0; work = complex(zeros(5,1)); anorm = norm(a,1); fprintf('Computed norm of a = %6.4g\n', anorm); [icase, x, estnrm, work, ifail] = nag_linsys_complex_norm_rcomm(icase, x, estnrm, work); while (icase ~= 0) if (icase == 1) x = inv(a)*x; else x = conj(transpose(inv(a)))*x; end [icase, x, estnrm, work, ifail] = nag_linsys_complex_norm_rcomm(icase, x, estnrm, work); end fprintf('Estimated norm of inverse(a) = %6.4g\n', estnrm); fprintf('Estimated condition number of A = %6.1f\n', estnrm*anorm);

Computed norm of a = 23.49 Estimated norm of inverse(a) = 37.04 Estimated condition number of A = 870.0

Open in the MATLAB editor: f04zc_example

function f04zc_examplea = [ 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]; icase = int64(0); x = complex(zeros(5, 1)); estnrm = 0; work = complex(zeros(5,1)); anorm = norm(a,1); fprintf('Computed norm of a = %6.4g\n', anorm); [icase, x, estnrm, work, ifail] = f04zc(icase, x, estnrm, work); while (icase ~= 0) if (icase == 1) x = inv(a)*x; else x = conj(transpose(inv(a)))*x; end [icase, x, estnrm, work, ifail] = f04zc(icase, x, estnrm, work); end fprintf('Estimated norm of inverse(a) = %6.4g\n', estnrm); fprintf('Estimated condition number of A = %6.1f\n', estnrm*anorm);

Computed norm of a = 23.49 Estimated norm of inverse(a) = 37.04 Estimated condition number of A = 870.0

© The Numerical Algorithms Group Ltd, Oxford, UK. 2009–2013