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

nag_linsys_real_norm_rcomm (f04yc) estimates the 1$1$-norm of a real 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_real_norm_rcomm (f04yc) computes an estimate (a lower bound) for the 1$1$-norm

of an n$n$ by n$n$ real 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^{T}x${A}^{\mathrm{T}}x$. 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^{ − T}x${B}^{-\mathrm{T}}x$ required by nag_linsys_real_norm_rcomm (f04yc) 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 by Higham (1988).

Since ‖A‖_{∞} = ‖A^{T}‖_{1}${\Vert A\Vert}_{\infty}={\Vert {A}^{\mathrm{T}}\Vert}_{1}$, nag_linsys_real_norm_rcomm (f04yc) can be used to estimate the ∞$\infty $-norm of A$A$ by working with A^{T}${A}^{\mathrm{T}}$ instead of A$A$.

The algorithm used is based on a method given by Hager (1984) and is described by Higham (1988). A comparison of several techniques for condition number estimation is given by 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) – double array
*On initial entry*: need not be set.- 3: estnrm – double scalar
*On initial entry*: need not be set.- 4: work(n) – double array
*On initial entry*: need not be set.- 5: iwork(n) – int64int32nag_int array

- 1: n – int64int32nag_int scalar
*Default*: The dimension of the arrays x, work, iwork. (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 ^{T}x${A}^{\mathrm{T}}x$ (if icase = 2${\mathbf{icase}}=2$),(b) place the result in x, and (c) call nag_linsys_real_norm_rcomm (f04yc) once again, with all the other parameters unchanged. - 2: x(n) – double 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) – double array
- 5: iwork(n) – int64int32nag_int array
- 6: 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 within nag_linsys_real_norm_rcomm (f04yc) is proportional to n$n$. For most problems the time taken during calls to nag_linsys_real_norm_rcomm (f04yc) will be negligible compared with the time spent evaluating matrix-vector products between calls to nag_linsys_real_norm_rcomm (f04yc).

The number of matrix-vector products required varies from 4$4$ to 11$11$ (or is 1$1$ if n = 1$n=1$). In most cases 4$4$ or 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_real_norm_rcomm (f04yc) 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, iwork, ifail] = f04yc(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] = f04yc(icase, x, estnrm, work, iwork); end end

To compute A^{ − 1}x${A}^{-1}x$ or A^{ − T}x${A}^{-\mathrm{T}}x$, solve the equation Ay = x$Ay=x$ or A^{T}y = x${A}^{\mathrm{T}}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 **symmetric** matrix, then A = A^{T}$A={A}^{\mathrm{T}}$, and the ifstatement after the while can be reduced to:

... code to compute inv(A)*x ...

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_norm_rcomm (f04yc) can be used in conjunction with NAG Toolbox functions for two important types of matrix: full nonsymmetric matrices (factorized by nag_lapack_dgetrf (f07ad)) and sparse nonsymmetric matrices (factorized by nag_matop_real_gen_sparse_lu (f01br)).

It is straightforward to use nag_linsys_real_norm_rcomm (f04yc) for the following other types of matrix, using the named functions for factorization and solution:

- nonsymmetric tridiagonal (nag_matop_real_gen_tridiag_lu (f01le) and nag_linsys_real_tridiag_fac_solve (f04le));
- nonsymmetric almost block-diagonal (nag_matop_real_gen_blkdiag_lu (f01lh) and nag_linsys_real_blkdiag_fac_solve (f04lh));
- nonsymmetric band (nag_lapack_dgbtrf (f07bd) and nag_lapack_dgbtrs (f07be));
- symmetric positive definite (nag_det_real_sym (f03bf), or nag_lapack_dpotrf (f07fd) and nag_lapack_dpotrs (f07fe));
- symmetric positive definite band (nag_lapack_dpbtrf (f07hd) and nag_lapack_dpbtrs (f07he));
- symmetric positive definite tridiagonal (nag_lapack_dptsv (f07ja), nag_lapack_dpttrf (f07jd) and nag_lapack_dpttrs (f07je));
- symmetric positive definite variable bandwidth (nag_matop_real_vband_posdef_fac (f01mc) and nag_linsys_real_posdef_vband_solve (f04mc));
- symmetric positive definite sparse (nag_sparse_real_symm_precon_ichol (f11ja) and nag_sparse_real_symm_precon_ichol_solve (f11jb));
- symmetric indefinite (nag_lapack_dsptrf (f07pd) and nag_lapack_dsptrs (f07pe)).

Open in the MATLAB editor: nag_linsys_real_norm_rcomm_example

function nag_linsys_real_norm_rcomm_examplea = [ 1.5, 2.0, 3.0, -2.1, 0.3; 2.5, 3.0, -4.0, 2.3, -1.1; 3.5, 4.0, 0.5, -3.1, -1.4; -0.4, -3.2, -2.1, 3.1, 2.1; 1.7, 3.7, 1.9, -2.2, -3.3]; icase = int64(0); x = zeros(5, 1); estnrm = 0; work = zeros(5, 1); iwork = zeros(5, 1, 'int64'); [icase, x, estnrm, work, iwork, ifail] = ... nag_linsys_real_norm_rcomm(icase, x, estnrm, work, iwork); while (icase > 0) if (icase ==1) x = a*x; elseif (icase == 2) x = transpose(a)*x; end [icase, x, estnrm, work, iwork, ifail] = ... nag_linsys_real_norm_rcomm(icase, x, estnrm, work, iwork); end estnrm

estnrm = 15.9000

Open in the MATLAB editor: f04yc_example

function f04yc_examplea = [ 1.5, 2.0, 3.0, -2.1, 0.3; 2.5, 3.0, -4.0, 2.3, -1.1; 3.5, 4.0, 0.5, -3.1, -1.4; -0.4, -3.2, -2.1, 3.1, 2.1; 1.7, 3.7, 1.9, -2.2, -3.3]; icase = int64(0); x = zeros(5, 1); estnrm = 0; work = zeros(5, 1); iwork = zeros(5, 1, 'int64'); [icase, x, estnrm, work, iwork, ifail] = f04yc(icase, x, estnrm, work, iwork); while (icase > 0) if (icase ==1) x = a*x; elseif (icase == 2) x = transpose(a)*x; end [icase, x, estnrm, work, iwork, ifail] = f04yc(icase, x, estnrm, work, iwork); end estnrm

estnrm = 15.9000

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