Fast Auxiliary Space Preconditioning 2.7.7 Aug/28/2022
Loading...
Searching...
No Matches
BlaILU.c File Reference

Incomplete LU decomposition: ILUk, ILUt, ILUtp. More...

#include <math.h>
#include <time.h>
#include "fasp.h"
#include "fasp_functs.h"

Go to the source code of this file.

Functions

void fasp_iluk (INT n, REAL *a, INT *ja, INT *ia, INT lfil, REAL *alu, INT *jlu, INT iwk, INT *ierr, INT *nzlu)
 Get ILU factorization with level of fill-in k (ilu(k)) for a CSR matrix A. More...
 
void fasp_ilut (INT n, REAL *a, INT *ja, INT *ia, INT lfil, REAL droptol, REAL *alu, INT *jlu, INT iwk, INT *ierr, INT *nz)
 Get incomplete LU factorization with dual truncations of a CSR matrix A. More...
 
void fasp_ilutp (INT n, REAL *a, INT *ja, INT *ia, INT lfil, REAL droptol, REAL permtol, INT mbloc, REAL *alu, INT *jlu, INT *iperm, INT iwk, INT *ierr, INT *nz)
 Get incomplete LU factorization with pivoting dual truncations of a CSR matrix A. More...
 
void fasp_symbfactor (INT n, INT *colind, INT *rwptr, INT levfill, INT nzmax, INT *nzlu, INT *ijlu, INT *uptr, INT *ierr)
 Symbolic factorization of a CSR matrix A in compressed sparse row format, with resulting factors stored in a single MSR data structure. More...
 

Detailed Description

Incomplete LU decomposition: ILUk, ILUt, ILUtp.

Note
This file contains Level-1 (Bla) functions. It requires: AuxMemory.c

Translated from SparseKit (Fortran code) by Chunsheng Feng, 09/03/2016


Copyright (C) 2016–Present by the FASP team. All rights reserved.

Released under the terms of the GNU Lesser General Public License 3.0 or later.

Definition in file BlaILU.c.

Function Documentation

◆ fasp_iluk()

void fasp_iluk ( INT  n,
REAL a,
INT ja,
INT ia,
INT  lfil,
REAL alu,
INT jlu,
INT  iwk,
INT ierr,
INT nzlu 
)

Get ILU factorization with level of fill-in k (ilu(k)) for a CSR matrix A.

Parameters
nrow number of A
anonzero entries of A
jainteger array of column for A
iainteger array of row pointers for A
lfilinteger. The fill-in parameter. Each row of L and each row of U will have a maximum of lfil elements (excluding the diagonal element). lfil must be .ge. 0.
alumatrix stored in Modified Sparse Row (MSR) format containing the L and U factors together. The diagonal (stored in alu(1:n) ) is inverted. Each i-th row of the alu,jlu matrix contains the i-th row of L (excluding the diagonal entry=1) followed by the i-th row of U.
jluinteger array of length n containing the pointers to the beginning of each row of U in the matrix alu,jlu.
iwkinteger. The minimum length of arrays alu, jlu, and levs.
ierrinteger pointer. Return error message with the following meaning. 0 --> successful return. >0 --> zero pivot encountered at step number ierr. -1 --> Error. input matrix may be wrong. (The elimination process has generated a row in L or U whose length is .gt. n.) -2 --> The matrix L overflows the array al. -3 --> The matrix U overflows the array alu. -4 --> Illegal value for lfil. -5 --> zero row encountered.
nzluinteger pointer. Return number of nonzero entries for alu and jlu
Note
: All the diagonal elements of the input matrix must be nonzero.
Author
Chunsheng Feng
Date
09/06/2016

Definition at line 72 of file BlaILU.c.

◆ fasp_ilut()

void fasp_ilut ( INT  n,
REAL a,
INT ja,
INT ia,
INT  lfil,
REAL  droptol,
REAL alu,
INT jlu,
INT  iwk,
INT ierr,
INT nz 
)

Get incomplete LU factorization with dual truncations of a CSR matrix A.

Parameters
nrow number of A
anonzero entries of A
jainteger array of column for A
iainteger array of row pointers for A
lfilinteger. The fill-in parameter. Each row of L and each row of U will have a maximum of lfil elements (excluding the diagonal element). lfil must be .ge. 0.
droptolreal*8. Sets the threshold for dropping small terms in the factorization. See below for details on dropping strategy.
alumatrix stored in Modified Sparse Row (MSR) format containing the L and U factors together. The diagonal (stored in alu(1:n) ) is inverted. Each i-th row of the alu,jlu matrix contains the i-th row of L (excluding the diagonal entry=1) followed by the i-th row of U.
jluinteger array of length n containing the pointers to the beginning of each row of U in the matrix alu,jlu.
iwkinteger. The lengths of arrays alu and jlu. If the arrays are not big enough to store the ILU factorizations, ilut will stop with an error message.
ierrinteger pointer. Return error message with the following meaning. 0 --> successful return. >0 --> zero pivot encountered at step number ierr. -1 --> Error. input matrix may be wrong. (The elimination process has generated a row in L or U whose length is .gt. n.) -2 --> The matrix L overflows the array al. -3 --> The matrix U overflows the array alu. -4 --> Illegal value for lfil. -5 --> zero row encountered.
nzinteger pointer. Return number of nonzero entries for alu and jlu
Note
All the diagonal elements of the input matrix must be nonzero.
Author
Chunsheng Feng
Date
09/06/2016

Definition at line 467 of file BlaILU.c.

◆ fasp_ilutp()

void fasp_ilutp ( INT  n,
REAL a,
INT ja,
INT ia,
INT  lfil,
REAL  droptol,
REAL  permtol,
INT  mbloc,
REAL alu,
INT jlu,
INT iperm,
INT  iwk,
INT ierr,
INT nz 
)

Get incomplete LU factorization with pivoting dual truncations of a CSR matrix A.

Parameters
nrow number of A
anonzero entries of A
jainteger array of column for A
iainteger array of row pointers for A
lfilinteger. The fill-in parameter. Each row of L and each row of U will have a maximum of lfil elements (excluding the diagonal element). lfil must be .ge. 0.
droptolreal*8. Sets the threshold for dropping small terms in the factorization. See below for details on dropping strategy.
permtoltolerance ratio used to determne whether or not to permute two columns. At step i columns i and j are permuted when abs(a(i,j))*permtol .gt. abs(a(i,i)) [0 --> never permute; good values 0.1 to 0.01]
mblocinteger.If desired, permuting can be done only within the diagonal blocks of size mbloc. Useful for PDE problems with several degrees of freedom.. If feature not wanted take mbloc=n.
alumatrix stored in Modified Sparse Row (MSR) format containing the L and U factors together. The diagonal (stored in alu(1:n) ) is inverted. Each i-th row of the alu,jlu matrix contains the i-th row of L (excluding the diagonal entry=1) followed by the i-th row of U.
jluinteger array of length n containing the pointers to the beginning of each row of U in the matrix alu,jlu.
ipermpermutation arrays
iwkinteger. The lengths of arrays alu and jlu. If the arrays are not big enough to store the ILU factorizations, ilut will stop with an error message.
ierrinteger pointer. Return error message with the following meaning. 0 --> successful return. >0 --> zero pivot encountered at step number ierr. -1 --> Error. input matrix may be wrong. (The elimination process has generated a row in L or U whose length is .gt. n.) -2 --> The matrix L overflows the array al. -3 --> The matrix U overflows the array alu. -4 --> Illegal value for lfil. -5 --> zero row encountered.
nzinteger pointer. Return number of nonzero entries for alu and jlu
Note
: All the diagonal elements of the input matrix must be nonzero.
Author
Chunsheng Feng
Date
09/06/2016

Definition at line 906 of file BlaILU.c.

◆ fasp_symbfactor()

void fasp_symbfactor ( INT  n,
INT colind,
INT rwptr,
INT  levfill,
INT  nzmax,
INT nzlu,
INT ijlu,
INT uptr,
INT ierr 
)

Symbolic factorization of a CSR matrix A in compressed sparse row format, with resulting factors stored in a single MSR data structure.

Parameters
nrow number of A
colindinteger array of column for A
rwptrinteger array of row pointers for A
levfillinteger. Level of fill-in allowed
nzmaxinteger. The maximum number of nonzero entries in the approximate factorization of a. This is the amount of storage allocated for ijlu.
nzluinteger pointer. Return number of nonzero entries for alu and jlu
ijluinteger array of length nzlu containing pointers to delimit rows and specify column number for stored elements of the approximate factors of A. the L and U factors are stored as one matrix.
uptrinteger array of length n containing the pointers to upper trig matrix
ierrinteger pointer. Return error message with the following meaning. 0 --> successful return. 1 --> not enough storage; check mneed.
Author
Chunsheng Feng
Date
09/06/2016

Symbolic factorization of a matrix in compressed sparse row format, * with resulting factors stored in a single MSR data structure. *

This routine uses the CSR data structure of A in two integer vectors * colind, rwptr to set up the data structure for the ILU(levfill) * factorization of A in the integer vectors ijlu and uptr. Both L * and U are stored in the same structure, and uptr(i) is the pointer * to the beginning of the i-th row of U in ijlu. *


Method Used * =========== *

The implementation assumes that the diagonal entries are * nonzero, and remain nonzero throughout the elimination * process. The algorithm proceeds row by row. When computing * the sparsity pattern of the i-th row, the effect of row * operations from previous rows is considered. Only those * preceding rows j for which (i,j) is nonzero need be considered, * since otherwise we would not have formed a linear combination * of rows i and j. *

The method used has some variations possible. The definition * of ILU(s) is not well specified enough to get a factorization * that is uniquely defined, even in the sparsity pattern that * results. For s = 0 or 1, there is not much variation, but for * higher levels of fill the problem is as follows: Suppose * during the decomposition while computing the nonzero pattern * for row i the following principal submatrix is obtained: * _______________________ * | | | * | | | * | j,j | j,k | * | | | * |__________|___________| * | | | * | | | * | i,j | i,k | * | | | * |__________|___________| *

Furthermore, suppose that entry (i,j) resulted from an earlier * fill-in and has level s1, and (j,k) resulted from an earlier * fill-in and has level s2: * _______________________ * | | | * | | | * | level 0 | level s2 | * | | | * |__________|___________| * | | | * | | | * | level s1 | | * | | | * |__________|___________| *

When using A(j,j) to annihilate A(i,j), fill-in will be incurred * in A(i,k). How should its level be defined? It would not be * operated on if A(i,j) or A(j,m) had not been filled in. The * version used here is to define its level as s1 + s2 + 1. However, * other reasonable choices would have been min(s1,s2) or max(s1,s2). * Using the sum gives a more conservative strategy in terms of the * growth of the number of nonzeros as s increases. *

levels(n+2:nzlu ) stores the levels from previous rows, * that is, the s2's above. levels(1:n) stores the fill-levels * of the current row (row i), which are the s1's above. * levels(n+1) is not used, so levels is conformant with MSR format. *

Vectors used: * ============= *

lastcol(n): * The integer lastcol(k) is the row index of the last row * to have a nonzero in column k, including the current * row, and fill-in up to this point. So for the matrix *

|-----------------------—| * | 11 12 15 | * | 21 22 26| * | 32 33 34 | * | 41 43 44 | * | 52 54 55 56| * | 62 66| * ------------------------— *

after step 1, lastcol() = [1 0 0 0 1 0] * after step 2, lastcol() = [2 2 0 0 2 2] * after step 3, lastcol() = [2 3 3 3 2 3] * after step 4, lastcol() = [4 3 4 4 4 3] * after step 5, lastcol() = [4 5 4 5 5 5] * after step 6, lastcol() = [4 6 4 5 5 6] *

Note that on step 2, lastcol(5) = 2 because there is a * fillin position (2,5) in the matrix. lastcol() is used * to determine if a nonzero occurs in column j because * it is a nonzero in the original matrix, or was a fill. *

rowll(n): * The integer vector rowll is used to keep a linked list of * the nonzeros in the current row, allowing fill-in to be * introduced sensibly. rowll is initialized with the * original nonzeros of the current row, and then sorted * using a shell sort. A pointer called head * (what ingenuity) is initialized. Note that at any * point rowll may contain garbage left over from previous * rows, which the linked list structure skips over. * For row 4 of the matrix above, first rowll is set to * rowll() = [3 1 2 5 - -], where - indicates any integer. * Then the vector is sorted, which yields * rowll() = [1 2 3 5 - -]. The vector is then expanded * to linked list form by setting head = 1 and * rowll() = [2 3 5 - 7 -], where 7 indicates termination. *

ijlu(nzlu): * The returned nonzero structure for the LU factors. * This is built up row by row in MSR format, with both L * and U stored in the data structure. Another vector, uptr(n), * is used to give pointers to the beginning of the upper * triangular part of the LU factors in ijlu. *

levels(n+2:nzlu): * This vector stores the fill level for each entry from * all the previous rows, used to compute if the current entry * will exceed the allowed levels of fill. The value in * levels(m) is added to the level of fill for the element in * the current row that is being reduced, to figure if * a column entry is to be accepted as fill, or rejected. * See the method explanation above. *

levels(1:n): * This vector stores the fill level number for the current * row's entries. If they were created as fill elements * themselves, this number is added to the corresponding * entry in levels(n+2:nzlu) to see if a particular column * entry will * be created as new fill or not. NOTE: in practice, the * value in levels(1:n) is one larger than the "fill" level of * the corresponding row entry, except for the diagonal * entry. That is why the accept/reject test in the code * is "if (levels(j) + levels(m) .le. levfill + 1)". *


on entry:

n = The order of the matrix A. ija = Integer array. Matrix A stored in modified sparse row format. levfill = Integer. Level of fill-in allowed. nzmax = Integer. The maximum number of nonzero entries in the approximate factorization of a. This is the amount of storage allocated for ijlu.

on return:

nzlu = The actual number of entries in the approximate factors, plus one. ijlu = Integer array of length nzlu containing pointers to delimit rows and specify column number for stored elements of the approximate factors of a. the l and u factors are stored as one matrix. uptr = Integer array of length n containing the pointers to upper trig matrix

ierr is an error flag: ierr = -i --> near zero pivot in step i ierr = 0 --> all's OK ierr = 1 --> not enough storage; check mneed. ierr = 2 --> illegal parameter

mneed = contains the actual number of elements in ldu, or the amount of additional storage needed for ldu

work arrays:

lastcol = integer array of length n containing last update of the corresponding column. levels = integer array of length n containing the level of fill-in in current row in its first n entries, and level of fill of previous rows of U in remaining part. rowll = integer array of length n containing pointers to implement a linked list for the fill-in elements.

external functions:

ifix, float, min0, srtr


Definition at line 1372 of file BlaILU.c.