NAG FL Interface
f11dff (real_​gen_​precon_​bdilu)

Settings help

FL Name Style:


FL Specification Language:


1 Purpose

f11dff computes a block diagonal incomplete LU factorization of a real sparse nonsymmetric matrix, represented in coordinate storage format. The diagonal blocks may be composed of arbitrary rows and the corresponding columns, and may overlap. This factorization can be used to provide a block Jacobi or additive Schwarz preconditioner, for use in combination with f11bef or f11dgf.

2 Specification

Fortran Interface
Subroutine f11dff ( n, nnz, a, la, irow, icol, nb, istb, indb, lindb, lfill, dtol, pstrat, milu, ipivp, ipivq, istr, idiag, nnzc, npivm, iwork, liwork, ifail)
Integer, Intent (In) :: n, nnz, la, nb, istb(nb+1), indb(lindb), lindb, lfill(nb), liwork
Integer, Intent (Inout) :: irow(la), icol(la), ipivp(lindb), ipivq(lindb), ifail
Integer, Intent (Out) :: istr(lindb+1), idiag(lindb), nnzc, npivm(nb), iwork(liwork)
Real (Kind=nag_wp), Intent (In) :: dtol(nb)
Real (Kind=nag_wp), Intent (Inout) :: a(la)
Character (1), Intent (In) :: pstrat(nb), milu(nb)
C Header Interface
#include <nag.h>
void  f11dff_ (const Integer *n, const Integer *nnz, double a[], const Integer *la, Integer irow[], Integer icol[], const Integer *nb, const Integer istb[], const Integer indb[], const Integer *lindb, const Integer lfill[], const double dtol[], const char pstrat[], const char milu[], Integer ipivp[], Integer ipivq[], Integer istr[], Integer idiag[], Integer *nnzc, Integer npivm[], Integer iwork[], const Integer *liwork, Integer *ifail, const Charlen length_pstrat, const Charlen length_milu)
The routine may be called by the names f11dff or nagf_sparse_real_gen_precon_bdilu.

3 Description

f11dff computes an incomplete LU factorization (see Meijerink and Van der Vorst (1977) and Meijerink and Van der Vorst (1981)) of the (possibly overlapping) diagonal blocks Ab, for b=1,2,,nb, of a real sparse nonsymmetric n×n matrix A. The factorization is intended primarily for use as a block Jacobi or additive Schwarz preconditioner (see Saad (1996)), with one of the iterative solvers f11bef and f11dgf.
The nb diagonal blocks need not consist of consecutive rows and columns of A, but may be composed of arbitrarily indexed rows, and the corresponding columns, as defined in the arguments indb and istb. Any given row or column index may appear in more than one diagonal block, resulting in overlap. Each diagonal block Ab, for b=1,2,,nb, is factorized as:
Ab = Mb+Rb  
where
Mb = Pb Lb Db Ub Qb  
and Lb is lower triangular with unit diagonal elements, Db is diagonal, Ub is upper triangular with unit diagonals, Pb and Qb are permutation matrices, and Rb is a remainder matrix.
The amount of fill-in occurring in the factorization of block b can vary from zero to complete fill, and can be controlled by specifying either the maximum level of fill lfill(b), or the drop tolerance dtol(b).
The parameter pstrat(b) defines the pivoting strategy to be used in block b. The options currently available are no pivoting, user-defined pivoting, partial pivoting by columns for stability, and complete pivoting by rows for sparsity and by columns for stability. The factorization may optionally be modified to preserve the row-sums of the original block matrix.
The sparse matrix A is represented in coordinate storage (CS) format (see Section 2.1.1 in the F11 Chapter Introduction). The array a stores all the nonzero elements of the matrix A, while arrays irow and icol store the corresponding row and column indices respectively. Multiple nonzero elements may not be specified for the same row and column index.
The preconditioning matrices Mb, for b=1,2,,nb, are returned in terms of the CS representations of the matrices
Cb = Lb + D-1b + Ub -2I .  

4 References

Meijerink J and Van der Vorst H (1977) An iterative solution method for linear systems of which the coefficient matrix is a symmetric M-matrix Math. Comput. 31 148–162
Meijerink J and Van der Vorst H (1981) Guidelines for the usage of incomplete decompositions in solving sets of linear equations as they occur in practical problems J. Comput. Phys. 44 134–155
Saad Y (1996) Iterative Methods for Sparse Linear Systems PWS Publishing Company, Boston, MA

5 Arguments

1: n Integer Input
On entry: n, the order of the matrix A.
Constraint: n1.
2: nnz Integer Input
On entry: the number of nonzero elements in the matrix A.
Constraint: 1nnzn2.
3: a(la) Real (Kind=nag_wp) array Input/Output
On entry: the nonzero elements in the matrix A, ordered by increasing row index, and by increasing column index within each row. Multiple entries for the same row and column indices are not permitted. The routine f11zaf may be used to order the elements in this way.
On exit: the first nnz entries of a contain the nonzero elements of A and the next nnzc entries contain the elements of the matrices Cb, for b=1,2,,nb stored consecutively. Within each block the matrix elements are ordered by increasing row index, and by increasing column index within each row.
4: la Integer Input
On entry: the dimension of the arrays a, irow and icol as declared in the (sub)program from which f11dff is called. These arrays must be of sufficient size to store both A (nnz elements) and C (nnzc elements).
Note: the minimum value for la is only appropriate if lfill and dtol are set such that minimal fill-in occurs. If this is not the case then we recommend that la is set much larger than the minimum value indicated in the constraint.
Constraint: la2×nnz.
5: irow(la) Integer array Input/Output
6: icol(la) Integer array Input/Output
On entry: the row and column indices of the nonzero elements supplied in a.
Constraints:
irow and icol must satisfy these constraints (which may be imposed by a call to f11zaf):
  • 1irow(i)n and 1icol(i)n, for i=1,2,,nnz;
  • either irow(i-1)<irow(i) or both irow(i-1)=irow(i) and icol(i-1)<icol(i), for i=2,3,,nnz.
On exit: the row and column indices of the nonzero elements returned in a.
7: nb Integer Input
On entry: the number of diagonal blocks to factorize.
Constraint: 1nbn.
8: istb(nb+1) Integer array Input
On entry: istb(b), for b=1,2,,nb, holds the indices in arrays indb, ipivp, ipivq and idiag that, on successful exit from this function, define block b. istb(nb+1) holds the sum of the number of rows in all blocks plus istb(1).
Constraint: istb(1)1, istb(b)< istb(b+1) , for b=1,2,,nb.
9: indb(lindb) Integer array Input
On entry: indb must hold the row indices appearing in each diagonal block, stored consecutively. Thus the elements indb(istb(b)) to indb(istb(b+1)-1) are the row indices in the bth block, for b=1,2,,nb.
Constraint: 1indb(m)n, for m=1,2,,istb(nb+1)-1.
10: lindb Integer Input
On entry: the dimension of the arrays indb, ipivp, ipivq and idiag as declared in the (sub)program from which f11dff is called.
Constraint: lindbistb(nb+1)-1.
11: lfill(nb) Integer array Input
On entry: if lfill(b)0 its value is the maximum level of fill allowed in the decomposition of the block (see Section 9.2 in f11daf). A negative value of lfill(b) indicates that dtol(b) will be used to control the fill in the block instead.
12: dtol(nb) Real (Kind=nag_wp) array Input
On entry: if lfill(b)<0 then dtol(b) is used as a drop tolerance in the block to control the fill-in (see Section 9.2 in f11daf); otherwise dtol(b) is not referenced.
Constraint: if lfill(b)<0, dtol(b)0.0, for b=1,2,,nb.
13: pstrat(nb) Character(1) array Input
On entry: pstrat(b), for b=1,2,,nb, specifies the pivoting strategy to be adopted in the block as follows:
pstrat(b)='N'
No pivoting is carried out.
pstrat(b)='U'
Pivoting is carried out according to the user-defined input values of ipivp and ipivq.
pstrat(b)='P'
Partial pivoting by columns for stability is carried out.
pstrat(b)='C'
Complete pivoting by rows for sparsity, and by columns for stability, is carried out.
Suggested value: pstrat(b)='C', for b=1,2,,nb.
Constraint: pstrat(b)='N', 'U', 'P' or 'C', for b=1,2,,nb.
14: milu(nb) Character(1) array Input
On entry: milu(b), for b=1,2,,nb, indicates whether or not the factorization in the block should be modified to preserve row-sums (see Section 9.4 in f11daf).
milu(b)='M'
The factorization is modified.
milu(b)='N'
The factorization is not modified.
Constraint: milu(b)='M' or 'N', for b=1,2,,nb.
15: ipivp(lindb) Integer array Input/Output
16: ipivq(lindb) Integer array Input/Output
On entry: if pstrat(b)='U', ipivp(istb(b)+k-1) and ipivq(istb(b)+k-1) must specify the row and column indices of the element used as a pivot at elimination stage k of the factorization of the block. Otherwise ipivp and ipivq need not be initialized.
Constraint: if pstrat(b)='U', the elements istb(b) to istb(b+1)-1 of ipivp and ipivq must both hold valid permutations of the integers on [1,istb(b+1)-istb(b)].
On exit: the row and column indices of the pivot elements, arranged consecutively for each block, as for indb. If ipivp(istb(b)+k-1)=i and ipivq(istb(b)+k-1)=j, the element in row i and column j of Ab was used as the pivot at elimination stage k.
17: istr(lindb+1) Integer array Output
On exit: istr(istb(b)+k-1), gives the index in the arrays a, irow and icol of row k of the matrix Cb, for b=1,2,,nb and k=1,2,,istb(b+1)-istb(b).
istr(istb(nb+1)) contains nnz+nnzc+1.
18: idiag(lindb) Integer array Output
On exit: idiag(istb(b)+k-1), gives the index in the arrays a, irow and icol of the diagonal element in row k of the matrix Cb, for b=1,2,,nb and k=1,2,,istb(b+1)-istb(b).
19: nnzc Integer Output
On exit: the sum total number of nonzero elements in the matrices Cb, for b=1,2,,nb.
20: npivm(nb) Integer array Output
On exit: if npivm(b)>0 it gives the number of pivots which were modified during the factorization to ensure that Mb exists.
If npivm(b)=−1 no pivot modifications were required, but a local restart occurred (see Section 9.3 in f11daf). The quality of the preconditioner will generally depend on the returned values of npivm(b), for b=1,2,,nb.
If npivm(b) is large, for some block, the preconditioner may not be satisfactory. In this case it may be advantageous to call f11dff again with an increased value of lfill(b), a reduced value of dtol(b), or pstrat(b)='C'.
21: iwork(liwork) Integer array Workspace
22: liwork Integer Input
On entry: the dimension of the array iwork as declared in the (sub)program from which f11dff is called.
Constraint: liwork9×n+3.
23: ifail Integer Input/Output
On entry: ifail must be set to 0, −1 or 1 to set behaviour on detection of an error; these values have no effect when no error is detected.
A value of 0 causes the printing of an error message and program execution will be halted; otherwise program execution continues. A value of −1 means that an error message is printed while a value of 1 means that it is not.
If halting is not appropriate, the value −1 or 1 is recommended. If message printing is undesirable, then the value 1 is recommended. Otherwise, the value 0 is recommended. When the value -1 or 1 is used it is essential to test the value of ifail on exit.
On exit: ifail=0 unless the routine detects an error or a warning has been flagged (see Section 6).

6 Error Indicators and Warnings

If on entry ifail=0 or −1, explanatory error messages are output on the current error message unit (as defined by x04aaf).
Errors or warnings detected by the routine:
ifail=1
On entry, dtol(value)=value.
Constraint: dtol(b)0.0, for b=1,2,,nb.
On entry, for b=value, istb(b+1)=value and istb(b)=value.
Constraint: istb(b+1)>istb(b), for b=1,2,,nb.
On entry, indb(value)=value and n=value.
Constraint: 1indb(m)n, for m=1,2,,istb(nb+1)-1
On entry, istb(1)=value.
Constraint: istb(1)1.
On entry, la=value and nnz=value.
Constraint: la2×nnz.
On entry, lindb=value, istb(nb+1)-1=value and nb=value.
Constraint: lindbistb(nb+1)-1.
On entry, liwork=value.
Constraint: liworkvalue.
On entry, milu(value)=value.
Constraint: milu(b)='M' or 'N' for all b.
On entry, n=value.
Constraint: n1.
On entry, nb=value and n=value.
Constraint: 1nbn.
On entry, nnz=value.
Constraint: nnz1.
On entry, nnz=value and n=value.
Constraint: nnzn2.
On entry, pstrat(value)=value.
Constraint: pstrat(b)='N', 'U', 'P' or 'C' for all b.
ifail=2
On entry, element value of a was out of order.
On entry, icol(value)=value and n=value.
Constraint: 1icol(j)n, for j=1,2,,nnz.
On entry, irow(value)=value and n=value.
Constraint: 1irow(i)n, for i=1,2,,nnz.
On entry, location value of (irow,icol) was a duplicate.
ifail=3
On entry, the user-supplied value of ipivp for block value lies outside its range.
On entry, the user-supplied value of ipivp for block value was repeated.
On entry, the user-supplied value of ipivq for block value lies outside its range.
On entry, the user-supplied value of ipivq for block value was repeated.
ifail=4
The number of nonzero entries in the decomposition is too large.
The decomposition has been terminated before completion.
Either increase la, or reduce the fill by reducing lfill, or increasing dtol.
ifail=-99
An unexpected error has been triggered by this routine. Please contact NAG.
See Section 7 in the Introduction to the NAG Library FL Interface for further information.
ifail=-399
Your licence key may have expired or may not have been installed correctly.
See Section 8 in the Introduction to the NAG Library FL Interface for further information.
ifail=-999
Dynamic memory allocation failed.
See Section 9 in the Introduction to the NAG Library FL Interface for further information.

7 Accuracy

The accuracy of the factorization of each block Ab will be determined by the size of the elements that are dropped and the size of any modifications made to the pivot elements. If these sizes are small then the computed factors will correspond to a matrix close to Ab. The factorization can generally be made more accurate by increasing the level of fill lfill(b), or by reducing the drop tolerance dtol(b) with lfill(b)<0.
If f11dff is used in combination with f11bef or f11dgf, the more accurate the factorization the fewer iterations will be required. However, the cost of the decomposition will also generally increase.

8 Parallelism and Performance

Background information to multithreading can be found in the Multithreading documentation.
f11dff is threaded by NAG for parallel execution in multithreaded implementations of the NAG Library.
Please consult the X06 Chapter Introduction for information on how to control and interrogate the OpenMP environment used within this routine. Please also consult the Users' Note for your implementation for any additional implementation-specific information.

9 Further Comments

f11dff calls f11daf internally for each block Ab. The comments and advice provided in Section 9 in f11daf on timing, control of fill, algorithmic details, and choice of parameters, are all therefore, relevant to f11dff, if interpreted blockwise.

10 Example

This example program reads in a sparse matrix A and then defines a block partitioning of the row indices with a user-supplied overlap and computes an overlapping incomplete LU factorization suitable for use as an additive Schwarz preconditioner. Such a factorization is used for this purpose in the example program of f11dgf.

10.1 Program Text

Program Text (f11dffe.f90)

10.2 Program Data

Program Data (f11dffe.d)

10.3 Program Results

Program Results (f11dffe.r)