factorization is output in the form of four one-dimensional arrays: integer arrays il and iu and real-valued arrays lval and uval. These describe the sparsity pattern and numerical values in the

L

and

U

matrices. The minimum required dimensions of these arrays cannot be given as a simple function of the size arguments (order and number of nonzero values) of the matrix

A

. This is due to unpredictable fill-in created by partial pivoting. nag_sparse_direct_real_gen_lu (f11me) will, on return, indicate which dimensions of these arrays were not adequate for the computation or (in the case of one of them) give a firm bound. You should then allocate more storage and try again.

References

Demmel J W, Eisenstat S C, Gilbert J R, Li X S and Li J W H (1999) A supernodal approach to sparse partial pivoting SIAM J. Matrix Anal. Appl. 20 720–755

Demmel J W, Gilbert J R and Li X S (1999) An asynchronous parallel supernodal algorithm for sparse gaussian elimination SIAM J. Matrix Anal. Appl. 20 915–952

Parameters

Compulsory Input Parameters

1: $n$ – int64int32nag_int scalar: $n$ , the order of the matrix $A$ .

Constraint: $n \geq 0$ .
2: $irowix (:)$ – int64int32nag_int array: The dimension of the array irowix must be at least $nnz$ , the number of nonzeros of the sparse matrix $A$

The row index array of sparse matrix $A$ .
3: $a (:)$ – double array: The dimension of the array a must be at least $nnz$ , the number of nonzeros of the sparse matrix $A$

The array of nonzero values in the sparse matrix $A$ .
4: $iprm (7 \times n)$ – int64int32nag_int array: Contains the column permutation which defines the permutation $P_{c}$ and associated data structures as computed by function nag_sparse_direct_real_gen_setup (f11md).
5: $thresh$ – double scalar: Suggested value: $thresh = 1.0$ . Smaller values may result in a faster factorization, but the benefits are likely to be small in most cases. It might be possible to use $thresh = 0.0$ if you are confident about the stability of the factorization, for example, if $A$ is diagonally dominant.
The diagonal pivoting threshold, $t$ . At step $j$ of the Gaussian elimination, if $|A_{j j}| \geq t (\max_{i \geq j} |A_{i j}|)$ , use $A_{j j}$ as a pivot, otherwise use $\max_{i \geq j} |A_{i j}|$ . A value of $t = 1$ corresponds to partial pivoting, a value of $t = 0$ corresponds to always choosing the pivot on the diagonal (unless it is zero).

Constraint: $0.0 \leq thresh \leq 1.0$ .
6: $nzlmx$ – int64int32nag_int scalar: Indicates the available size of array il. The dimension of il should be at least $7 \times n + nzlmx + 4$ . A good range for nzlmx that works for many problems is $nnz$ to $8 \times nnz$ , where $nnz$ is the number of nonzeros in the sparse matrix $A$ . If, on exit, $ifail = 2$ , the given nzlmx was too small and you should attempt to provide more storage and call the function again.

Constraint: $nzlmx \geq 1$ .
7: $nzlumx$ – int64int32nag_int scalar: Indicates the available size of array lval. The dimension of lval should be at least nzlumx.

Constraint: $nzlumx \geq 1$ .
8: $nzumx$ – int64int32nag_int scalar: Indicates the available sizes of arrays iu and uval. The dimension of iu should be at least $2 \times n + nzumx + 1$ and the dimension of uval should be at least nzumx. A good range for nzumx that works for many problems is $nnz$ to $8 \times nnz$ , where $nnz$ is the number of nonzeros in the sparse matrix $A$ . If, on exit, $ifail = 3$ , the given nzumx was too small and you should attempt to provide more storage and call the function again.

Constraint: $nzumx \geq 1$ .

Optional Input Parameters

None.

Output Parameters

1: $iprm (7 \times n)$ – int64int32nag_int array: Part of the array is modified to record the row permutation $P_{r}$ determined by pivoting.
2: $nzlumx$ – int64int32nag_int scalar: If $ifail = 4$ , the given nzlumx was too small and is reset to a value that will be sufficient. You should then provide the indicated storage and call the function again.
3: $il (7 \times n + nzlmx + 4)$ – int64int32nag_int array: Encapsulates the sparsity pattern of matrix $L$ .
4: $lval (nzlumx)$ – double array: Records the nonzero values of matrix $L$ and some of the nonzero values of matrix $U$ .
5: $iu (2 \times n + nzumx + 1)$ – int64int32nag_int array: Encapsulates the sparsity pattern of matrix $U$ .
6: $uval (nzumx)$ – double array: Records some of the nonzero values of matrix $U$ .
7: $nnzl$ – int64int32nag_int scalar: The number of nonzero values in the matrix $L$ .
8: $nnzu$ – int64int32nag_int scalar: The number of nonzero values in the matrix $U$ .
9: $flop$ – double scalar: The number of floating-point operations performed.
10: $ifail$ – int64int32nag_int scalar: $ifail = 0$ unless the function detects an error (see Error Indicators and Warnings).

Error Indicators and Warnings

Errors or warnings detected by the function:

$ifail = 1$: Constraint: $0.0 \leq thresh \leq 1.0$ .

Constraint: $n \geq 0$ .

Constraint: $nzlmx \geq 1$ .

Constraint: $nzlumx \geq 1$ .

Constraint: $nzumx \geq 1$ .

$ifail = 2$: Insufficient nzlmx.

$ifail = 3$: Insufficient nzumx.

$ifail = 4$: Insufficient nzlumx.

$ifail = 5$: The matrix is singular – no factorization possible.

$ifail = - 99$: An unexpected error has been triggered by this routine. Please contact NAG.

$ifail = - 399$: Your licence key may have expired or may not have been installed correctly.

$ifail = - 999$: Dynamic memory allocation failed.

Accuracy

The computed factors

L

and

U

are the exact factors of a perturbed matrix

A + E

, where

|E| \leq c (n) ε |L| |U|,

c (n)

is a modest linear function of

n

, and

ε

is the machine precision, when partial pivoting is used. If no partial pivoting is used, the factorization accuracy can be considerably worse. A call to nag_sparse_direct_real_gen_diag (f11mm) after nag_sparse_direct_real_gen_lu (f11me) can help determine the quality of the factorization.

Further Comments

The total number of floating-point operations depends on the sparsity pattern of the matrix

A

A call to nag_sparse_direct_real_gen_lu (f11me) may be followed by calls to the functions:

nag_sparse_direct_real_gen_solve (f11mf) to solve $A X = B$ or $A^{T} X = B$ ;
nag_sparse_direct_real_gen_cond (f11mg) to estimate the condition number of $A$ ;
nag_sparse_direct_real_gen_diag (f11mm) to estimate the reciprocal pivot growth of the $L U$ factorization.

Example

This example computes the

L U

factorization of the matrix

A

, where

A = (\begin{matrix} 2.00 & 1.00 & 0 & 0 & 0 \\ 0 & 0 & 1.00 & - 1.00 & 0 \\ 4.00 & 0 & 1.00 & 0 & 1.00 \\ 0 & 0 & 0 & 1.00 & 2.00 \\ 0 & - 2.00 & 0 & 0 & 3.00 \end{matrix}) .

Open in the MATLAB editor: f11me_example

function f11me_example


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

% LU factorize sparse A and solve Ax = b
n  = int64(5);
nz = int64(11);
icolzp = int64([1; 3; 5; 7; 9; 12]);
irowix = int64([1; 3;
                  1; 5;
                  2; 3;
                  2; 4;
                  3; 4; 5]);
a = [ 2;  4;
      1; -2;
      1;  1;
     -1;  1;
      1;  2; 3];

% Permute matrix
iprm   = zeros(1, 7*n, 'int64');
spec   = 'M';
% Calculate COLAMD permutation
[iprm, ifail] = f11md( ...
                       spec, n, icolzp, irowix, iprm);

% Factorize
thresh = 1;
nzlmx  = int64(8*nz);
nzlumx = nzlmx;
nzumx  = nzlmx;
[iprm, nzlumx, il, lval, iu, uval, nnzl, nnzu, flop, ifail] = ...
    f11me(...
          n, irowix, a, iprm, thresh, nzlmx, nzlumx, nzumx);

fprintf('\nNumber of nonzeros in factors (excluding unit diagonal)');
fprintf('%8d\n',nnzl + nnzu - n);

disp('Factor elements in lval');
fprintf('%7.2f%7.2f%7.2f%7.2f%7.2f\n',lval(1:10)');
disp('Factor elements in uval');
fprintf('%7.2f%7.2f%7.2f%7.2f\n',uval(1:4)');

f11me example results


Number of nonzeros in factors (excluding unit diagonal)      14
Factor elements in lval
  -2.00  -0.50   4.00   0.50   2.00
   0.50  -1.00   0.50   1.00  -1.00
Factor elements in uval
   1.00   3.00   1.00   1.00

PDF version (NAG web site, 64-bit version, 64-bit version)

Chapter Contents

Chapter Introduction

NAG Toolbox