factorization of a complex sparse non-Hermitian 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 nag_sparse_nherm_basic_solver (f11bsc) or nag_sparse_nherm_precon_bdilu_solve (f11duc).

2

Specification

#include <nag.h>

#include <nagf11.h>

void

nag_sparse_nherm_precon_bdilu (Integer n, Integer nnz, Complex a[], Integer la, Integer irow[], Integer icol[], Integer nb, const Integer istb[], const Integer indb[], Integer lindb, const Integer lfill[], const double dtol[], const Nag_SparseNsym_Piv pstrat[], const Nag_SparseNsym_Fact milu[], Integer ipivp[], Integer ipivq[], Integer istr[], Integer idiag[], Integer *nnzc, Integer npivm[], NagError *fail)

3

Description

nag_sparse_nherm_precon_bdilu (f11dtc) computes an incomplete

L U

factorization (see Meijerink and Van der Vorst (1977) and Meijerink and Van der Vorst (1981)) of the (possibly overlapping) diagonal blocks

A_{b}

b = 1, 2, \dots, nb

, of a complex sparse non-Hermitian

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 nag_sparse_nherm_basic_solver (f11bsc) and nag_sparse_nherm_precon_bdilu_solve (f11duc).

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

A_{b}

b = 1, 2, \dots, nb

, is factorized as:

A_{b} = M_{b} + R_{b}

where

M_{b} = P_{b} L_{b} D_{b} U_{b} Q_{b}

and

L_{b}

is lower triangular with unit diagonal elements,

D_{b}

is diagonal,

U_{b}

is upper triangular with unit diagonals,

P_{b}

and

Q_{b}

are permutation matrices, and

R_{b}

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 - 1]

, or the drop tolerance

dtol [b - 1]

The parameter

pstrat [b - 1]

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

M_{b}

b = 1, 2, \dots, nb

, are returned in terms of the CS representations of the matrices

C_{b} = L_{b} + {D^{- 1}}_{b} + U_{b} - 2 I .

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$ – IntegerInput

On entry:

n

, the order of the matrix

A

Constraint:

n \geq 1

2: $nnz$ – IntegerInput

On entry: the number of nonzero elements in the matrix

A

Constraint:

1 \leq nnz \leq n^{2}

3: $a [la]$ – ComplexInput/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 function nag_sparse_nherm_sort (f11znc) 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

C_{b}

, for

b = 1, 2, \dots, 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$ – IntegerInput

On entry: the dimension of the arrays a, irow and icol. 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:

la \geq 2 \times nnz

5: $irow [la]$ – IntegerInput/Output

6: $icol [la]$ – IntegerInput/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 nag_sparse_nherm_sort (f11znc)):

$1 \leq irow [i - 1] \leq n$ and $1 \leq icol [i - 1] \leq n$ , for $i = 1, 2, \dots, nnz$ ;
either $irow [i - 1] < irow [i]$ or both $irow [i - 1] = irow [i]$ and $icol [i - 1] < icol [i]$ , for $i = 1, 2, \dots, nnz$ .

On exit: the row and column indices of the nonzero elements returned in a.

7: $nb$ – IntegerInput

On entry: the number of diagonal blocks to factorize.

Constraint:

1 \leq nb \leq n

8: $istb [nb + 1]$ – const IntegerInput

On entry:

istb [b - 1]

, for

b = 1, 2, \dots, nb

, holds the indices in arrays indb, ipivp, ipivq and idiag that, on successful exit from this function, define block

b

. Let

r_{b}

denote the number of rows in block

b

; then

istb [b] = istb [b - 1] + r_{b}

, for

b = 1, 2, \dots, nb

. Thus,

istb [nb]

holds the sum of the number of rows in all blocks plus

istb [0]

Constraint:

istb [0] \geq 1, istb [b - 1] < istb [b]

, for

b = 1, 2, \dots, nb

9: $indb [lindb]$ – const IntegerInput

On entry: indb must hold the row indices appearing in each diagonal block, stored consecutively. Thus the elements

indb [k_{b} - 1]

, for

k_{b} = istb [b - 1], istb [b - 1] + 1, \dots, istb [b] - 2, istb [b] - 1

, are the row indices in the

b

th block, for

b = 1, 2, \dots, nb

Constraint:

1 \leq indb [m - 1] \leq n

, for

m = istb [0], istb [0] + 1, \dots, istb [nb] - 1

10: $lindb$ – IntegerInput

On entry: the dimension of the arrays indb, ipivp, ipivq and idiag.

Constraint:

lindb \geq istb [nb] - 1

11: $lfill [nb]$ – const IntegerInput

On entry: if

lfill [b - 1] \geq 0

its value is the maximum level of fill allowed in the decomposition of the block

b

(see Section 9.2 in nag_sparse_nherm_fac (f11dnc)). A negative value of

lfill [b - 1]

indicates that

dtol [b - 1]

will be used to control the fill in block

b

instead.

12: $dtol [nb]$ – const doubleInput

On entry: if

lfill [b - 1] < 0

then

dtol [b - 1]

is used as a drop tolerance in block

b

to control the fill-in (see Section 9.2 in nag_sparse_nherm_fac (f11dnc)); otherwise

dtol [b - 1]

is not referenced.

Constraint: if

lfill [b - 1] < 0

dtol [b - 1] \geq 0.0

, for

b = 1, 2, \dots, nb

13: $pstrat [nb]$ – const Nag_SparseNsym_PivInput

On entry:

pstrat [b - 1]

, for

b = 1, 2, \dots, nb

, specifies the pivoting strategy to be adopted in the block as follows:

$pstrat [b - 1] = Nag_SparseNsym_NoPiv$: No pivoting is carried out.
$pstrat [b - 1] = Nag_SparseNsym_UserPiv$: Pivoting is carried out according to the user-defined input values of ipivp and ipivq.
$pstrat [b - 1] = Nag_SparseNsym_PartialPiv$: Partial pivoting by columns for stability is carried out.
$pstrat [b - 1] = Nag_SparseNsym_CompletePiv$: Complete pivoting by rows for sparsity, and by columns for stability, is carried out.

Suggested value:

pstrat [b - 1] = Nag_SparseNsym_CompletePiv

, for

b = 1, 2, \dots, nb

Constraint:

pstrat [b - 1] = Nag_SparseNsym_NoPiv

Nag_SparseNsym_UserPiv

Nag_SparseNsym_PartialPiv

Nag_SparseNsym_CompletePiv

, for

b = 1, 2, \dots, nb

14: $milu [nb]$ – const Nag_SparseNsym_FactInput

On entry:

milu [b - 1]

, for

b = 1, 2, \dots, nb

, indicates whether or not the factorization in block

b

should be modified to preserve row-sums (see Section 9.4 in nag_sparse_nherm_fac (f11dnc)).

$milu [b - 1] = Nag_SparseNsym_ModFact$: The factorization is modified.
$milu [b - 1] = Nag_SparseNsym_UnModFact$: The factorization is not modified.

Constraint:

milu [b - 1] = Nag_SparseNsym_ModFact

Nag_SparseNsym_UnModFact

, for

b = 1, 2, \dots, nb

15: $ipivp [lindb]$ – IntegerInput/Output

16: $ipivq [lindb]$ – IntegerInput/Output

On entry: if

pstrat [b - 1] = Nag_SparseNsym_UserPiv

ipivp [istb [b - 1] + k - 2]

and

ipivq [istb [b - 1] + k - 2]

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 - 1] = Nag_SparseNsym_UserPiv

, the elements

istb [b - 1] - 1

istb [b] - 2

of ipivp and ipivq must both hold valid permutations of the integers on

[1, istb [b] - istb [b - 1]]

On exit: the row and column indices of the pivot elements, arranged consecutively for each block, as for indb. If

ipivp [istb [b - 1] + k - 2] = i

and

ipivq [istb [b - 1] + k - 2] = j

, the element in row

i

and column

j

A_{b}

was used as the pivot at elimination stage

k

17: $istr [lindb + 1]$ – IntegerOutput

On exit:

istr [istb [b - 1] + k - 2]

, gives the index in the arrays a, irow and icol of row

k

of the matrix

C_{b}

, for

b = 1, 2, \dots, nb

and

k = 1, 2, \dots, istb [b] - istb [b - 1]

istr [istb [nb] - 1]

contains

nnz + nnzc + 1

18: $idiag [lindb]$ – IntegerOutput

On exit:

idiag [istb [b - 1] + k - 2]

, gives the index in the arrays a, irow and icol of the diagonal element in row

k

of the matrix

C_{b}

, for

b = 1, 2, \dots, nb

and

k = 1, 2, \dots, istb [b] - istb [b - 1]

19: $nnzc$ – Integer *Output

On exit: the sum total number of nonzero elements in the matrices

C_{b}

, for

b = 1, 2, \dots, nb

20: $npivm [nb]$ – IntegerOutput

On exit: if

npivm [b - 1] > 0

it gives the number of pivots which were modified during the factorization to ensure that

M_{b}

exists.

npivm [b - 1] = - 1

no pivot modifications were required, but a local restart occurred (see Section 9.3 in nag_sparse_nherm_fac (f11dnc)). The quality of the preconditioner will generally depend on the returned values of

npivm [b - 1]

, for

b = 1, 2, \dots, nb

npivm [b - 1]

is large, for some

b

, the preconditioner may not be satisfactory. In this case it may be advantageous to call nag_sparse_nherm_precon_bdilu (f11dtc) again with an increased value of

lfill [b - 1]

, a reduced value of

dtol [b - 1]

, or

pstrat [b - 1] = Nag_SparseNsym_CompletePiv

21: $fail$ – NagError *Input/Output

The NAG error argument (see Section 3.7 in How to Use the NAG Library and its Documentation).

6

Error Indicators and Warnings

NE_ALLOC_FAIL: Dynamic memory allocation failed.
See Section 2.3.1.2 in How to Use the NAG Library and its Documentation for further information.
NE_BAD_PARAM: On entry, argument $〈value〉$ had an illegal value.
NE_INT: On entry, $istb [0] = 〈value〉$ .
Constraint: $istb [0] \geq 1$ .

On entry, $n = 〈value〉$ .
Constraint: $n \geq 1$ .

On entry, $nnz = 〈value〉$ .
Constraint: $nnz \geq 1$ .
NE_INT_2: On entry, $la = 〈value〉$ and $nnz = 〈value〉$ .
Constraint: $la \geq 2 \times nnz$ .

On entry, $nb = 〈value〉$ and $n = 〈value〉$ .
Constraint: $1 \leq nb \leq n$ .

On entry, $nnz = 〈value〉$ and $n = 〈value〉$ .
Constraint: $nnz \leq n^{2}$ .
NE_INT_3: On entry, $lindb = 〈value〉$ , $istb [nb] - 1 = 〈value〉$ and $nb = 〈value〉$ .
Constraint: $lindb \geq istb [nb] - 1$ .
NE_INT_ARRAY: On entry, $indb [〈value〉] = 〈value〉$ and $n = 〈value〉$ .
Constraint: $1 \leq indb [m - 1] \leq n$ , for $m = istb [0], istb [0] + 1, \dots, istb [nb] - 1$ .
NE_INTERNAL_ERROR: An internal error has occurred in this function. Check the function call and any array sizes. If the call is correct then please contact NAG for assistance.
See Section 2.7.6 in How to Use the NAG Library and its Documentation for further information.
NE_INVALID_CS: On entry, $icol [〈value〉] = 〈value〉$ and $n = 〈value〉$ .
Constraint: $1 \leq icol [j - 1] \leq n$ , for $j = 1, 2, \dots, nnz$ .

On entry, $irow [〈value〉] = 〈value〉$ and $n = 〈value〉$ .
Constraint: $1 \leq irow [i - 1] \leq n$ , for $i = 1, 2, \dots, nnz$ .
NE_INVALID_ROWCOL_PIVOT: 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.
NE_NO_LICENCE: Your licence key may have expired or may not have been installed correctly.
See Section 2.7.5 in How to Use the NAG Library and its Documentation for further information.
NE_NOT_STRICTLY_INCREASING: On entry, element $〈value〉$ of a was out of order.

On entry, for $b = 〈value〉$ , $istb [b] = 〈value〉$ and $istb [b - 1] = 〈value〉$ .
Constraint: $istb [b] > istb [b - 1]$ , for $b = 1, 2, \dots, nb$ .

On entry, location $〈value〉$ of $(irow, icol)$ was a duplicate.
NE_REAL_ARRAY: On entry, $dtol [〈value〉] = 〈value〉$ .
Constraint: $dtol [b - 1] \geq 0.0$ , for $b = 1, 2, \dots, nb$ .
NE_TOO_SMALL: 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.

7

Accuracy

The accuracy of the factorization of each block

A_{b}

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

A_{b}

. The factorization can generally be made more accurate by increasing the level of fill

lfill [b - 1]

, or by reducing the drop tolerance

dtol [b - 1]

with

lfill [b - 1] < 0

If nag_sparse_nherm_precon_bdilu (f11dtc) is used in combination with nag_sparse_nherm_basic_solver (f11bsc) or nag_sparse_nherm_precon_bdilu_solve (f11duc), 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

nag_sparse_nherm_precon_bdilu (f11dtc) 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 function. Please also consult the Users' Note for your implementation for any additional implementation-specific information.

9

Further Comments

nag_sparse_nherm_precon_bdilu (f11dtc) calls nag_sparse_nherm_fac (f11dnc) internally for each block

A_{b}

. The comments and advice provided in Section 9 in nag_sparse_nherm_fac (f11dnc) on timing, control of fill, algorithmic details, and choice of parameters, are all therefore relevant to nag_sparse_nherm_precon_bdilu (f11dtc), 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

L U

factorization suitable for use as an additive Schwarz preconditioner. Such a factorization is used for this purpose in the example program of nag_sparse_nherm_precon_bdilu_solve (f11duc).

NAG Library Function Document

nag_sparse_nherm_precon_bdilu (f11dtc)

▸▿ Contents

1 Purpose

2 Specification

3 Description

4 References

5 Arguments

6 Error Indicators and Warnings

7 Accuracy

8 Parallelism and Performance

9 Further Comments

10 Example

10.1 Program Text

10.2 Program Data

10.3 Program Results

1

Purpose

2

Specification

3

Description

4

References

5

Arguments

6

Error Indicators and Warnings

7

Accuracy

8

Parallelism and Performance

9

Further Comments

10

Example

10.1

Program Text

10.2

Program Data

10.3

Program Results