Frequently, in practical problems the data  and  are not known exactly, and it is then important to understand how uncertainties or perturbations in the data can affect the solution.

If  is the exact solution to , and  is the exact solution to a perturbed problem , then

where  is the condition number of  defined by

(3)

In other words, relative errors in  or  may be amplified in  by a factor . Section [Estimating condition numbers] discusses how to compute or estimate .

Similar considerations apply when we study the effects of rounding errors introduced by computation in finite precision. The effects of rounding errors can be shown to be equivalent to perturbations in the original data, such that  and  are usually at most , where  is the machine precision and  is an increasing function of  which is seldom larger than  (although in theory it can be as large as ).

In other words, the computed solution  is the exact solution of a linear system  which is close to the original system in a normwise sense.

The previous section has emphasised the usefulness of the quantity  in understanding the sensitivity of the solution of . To compute the value of  from equation (3) is more expensive than solving  in the first place. Hence it is standard practice to estimate , in either the 1-norm or the -norm, by a method which only requires  additional operations, assuming that a suitable factorization of  is available.

The method used in this chapter is Higham's modification of Hager's method (see Higham (1988)). It yields an estimate which is never larger than the true value, but which seldom falls short by more than a factor of 3 (although artificial examples can be constructed where it is much smaller). This is acceptable since it is the order of magnitude of  which is important rather than its precise value.

Because  is infinite if  is singular, the functions in this chapter actually return the reciprocal of .

The condition of a matrix and hence the accuracy of the computed solution, may be improved by scaling; thus if  and  are diagonal matrices with positive diagonal elements, then

is the scaled matrix. A general matrix is said to be equilibrated if it is scaled so that the lengths of its rows and columns have approximately equal magnitude. Similarly a general matrix is said to be row-equilibrated (column-equilibrated) if it is scaled so that the lengths of its rows (columns) have approximately equal magnitude. Note that row scaling can affect the choice of pivot when partial pivoting is used in the factorization of .

A symmetric or Hermitian positive-definite matrix is said to be equilibrated if the diagonal elements are all approximately equal to unity.

For further information on scaling and equilibration see Section 3.5.2 of Golub and Van Loan (1996), Section 7.2, 7.3 and 9.8 of Higham (1988) and Section 5 of Chapter 4 of Wilkinson (1965).

Functions are provided to return the scaling factors that equilibrate a matrix for general, general band, symmetric and Hermitian positive-definite and symmetric and Hermitian positive-definite band matrices.

A disadvantage of normwise error bounds is that they do not reflect any special structure in the data  and  – that is, a pattern of elements which are known to be zero – and the bounds are dominated by the largest elements in the data.

Componentwise error bounds overcome these limitations. Instead of the normwise relative error, we can bound the relative error in each component of  and :

where the component-wise backward error bound  is given by

Functions are provided in this chapter which compute , and also compute a forward error bound which is sometimes much sharper than the normwise bound given earlier:

Care is taken when computing this bound to allow for rounding errors in computing . The norm  is estimated cheaply (without computing ) by a modification of the method used to estimate .

If  is an approximate computed solution to , and  is the corresponding residual, then a procedure for iterative refinement of  can be defined as follows, starting with :

In Chapter f04, functions are provided which perform this procedure using additional precision to compute , and are thus able to reduce the forward error to the level of machine precision.

The functions in this chapter do not use additional precision to compute , and cannot guarantee a small forward error, but can guarantee a small backward error (except in rare cases when  is very ill-conditioned, or when  and  are sparse in such a way that  has a zero or very small component). The iterations continue until the backward error has been reduced as much as possible; usually only one iteration is needed, and at most five iterations are allowed.

Matrices may be stored column-wise or row-wise as described in Section [ ] of the NAG Library Introduction: a matrix  is stored in a one-dimensional array a, with matrix element  stored column-wise in array element  or row-wise in array element  where pda is the principle dimension of the array (i.e., the stride separating row or column elements of the matrix respectively). Most functions in this chapter contain the storage order argument which can be set to Nag_ColMajor for column-wise storage or Nag_RowMajor for row-wise storage of matrices. Where groups of functions are intended to be used together, the value of the storage order argument passed must be consistent throughout.

If a matrix is triangular (upper or lower, as specified by the argument uplo), only the elements of the relevant triangle are stored; the remaining elements of the array need not be set. Such elements are indicated by * in the examples below.

For example, when :

Functions which handle symmetric or Hermitian matrices allow for either the upper or lower triangle of the matrix (as specified by uplo) to be stored in the corresponding elements of the array; the remaining elements of the array need not be set.

For example, when :

Symmetric, Hermitian or triangular matrices may be stored more compactly, if the relevant triangle (again as specified by uplo) is packed by columns or rows in a one-dimensional array. In Chapters f07 and f08, arrays which hold matrices in packed storage have names ending in . The storage of matrix elements  in the packed array ap is as follows:

For example:

Note that for real symmetric matrices, packing the upper triangle by columns is equivalent to packing the lower triangle by rows; packing the lower triangle by columns is equivalent to packing the upper triangle by rows. (For complex Hermitian matrices, the only difference is that the off-diagonal elements are conjugated.)

A band matrix with  subdiagonals and  superdiagonals may be stored compactly in a notional two-dimensional array with  rows and  columns if stored column-wise or  rows and  columns if stored row-wise. In column-major order, elements of a column of the matrix are stored contiguously in the array, and elements of the diagonals of the matrix are stored with constant stride (i.e., in a row of the two-dimensional array). In row-major order, elements of a row of the matrix are stored contiguously in the array, and elements of a diagonal of the matrix are stored with constant stride (i.e., in a column of the two-dimensional array). These storage schemes should only be used in practice if , , although the functions in Chapters f07 and f08 work correctly for all values of  and . In Chapters f07 and f08 arrays which hold matrices in band storage have names ending in .

To be precise, elements of matrix elements  are stored as follows:

where  is the stride between diagonal elements and where .

For example, when ,  and :

 The elements marked  in the upper left and lower right corners of the array ab need not be set, and are not referenced by the functions. In this example, if  and pdab takes the minimum value of 4, then  need not be set, . On the other hand, if  (), then  and  need not be set, .

Note: when a general band matrix is supplied for  factorization, space must be allowed to store an additional  superdiagonals, generated by fill-in as a result of row interchanges. This means that the matrix is stored according to the above scheme, but with  superdiagonals; it also means that the principal dimension has the constraint .

Triangular band matrices are stored in the same format, with either  if upper triangular, or  if lower triangular.

For symmetric or Hermitian band matrices with  subdiagonals or superdiagonals, only the upper or lower triangle (as specified by uplo) need be stored:

where  is the stride separating diagonal matrix elements in the array ab.

For example, when  and :

 Note that different storage schemes for band matrices are used by some functions in Chapters f01, f02, f03 and f04. In the above example, if  and , then for , ; while for , . If  (), then for , ; while for , .

Some functions in this chapter have an option to handle unit triangular matrices (that is, triangular matrices with diagonal elements ). This option is specified by an argument diag. If  (Unit triangular), the diagonal elements of the matrix need not be stored, and the corresponding array elements are not referenced by the functions. The storage scheme for the rest of the matrix (whether conventional, packed or band) remains unchanged.

Complex Hermitian matrices have diagonal elements that are by definition purely real. In addition, complex triangular matrices which arise in Cholesky factorization are defined by the algorithm to have real diagonal elements.

If such matrices are supplied as input to functions in this chapter, the imaginary parts of the diagonal elements are not referenced, but are assumed to be zero. If such matrices are returned as output by the functions, the computed imaginary parts are explicitly set to zero.

In addition to the storage order argument of type Nag_OrderType, most functions in this Chapter have one or more option arguments of various types; only options of the correct type may be supplied.

For example,

 f07fdc("Nag_Upper",...)

It is permissible for the problem dimensions (for example, m in f07adc (nag_dgetrf), n or nrhs in f07aec (nag_dgetrs)) to be passed as zero, in which case the computation (or part of it) is skipped. Negative dimensions are regarded as an error.

 Each entry gives:

Table 1

Each entry gives:

Table 2