![]() |
Fast Auxiliary Space Preconditioning 2.7.7 Aug/28/2022
|
Incomplete LU decomposition: ILUk, ILUt, ILUtp. More...
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... | |
Incomplete LU decomposition: ILUk, ILUt, ILUtp.
Translated from SparseKit (Fortran code) by Chunsheng Feng, 09/03/2016
Copyright (C) 2016–Present by the FASP team. All rights reserved.
Definition in file BlaILU.c.
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.
n | row number of A |
a | nonzero entries of A |
ja | integer array of column for A |
ia | integer array of row pointers for A |
lfil | integer. 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. |
alu | matrix 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. |
jlu | integer array of length n containing the pointers to the beginning of each row of U in the matrix alu,jlu. |
iwk | integer. The minimum length of arrays alu, jlu, and levs. |
ierr | integer 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. |
nzlu | integer pointer. Return number of nonzero entries for alu and jlu |
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.
n | row number of A |
a | nonzero entries of A |
ja | integer array of column for A |
ia | integer array of row pointers for A |
lfil | integer. 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. |
droptol | real*8. Sets the threshold for dropping small terms in the factorization. See below for details on dropping strategy. |
alu | matrix 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. |
jlu | integer array of length n containing the pointers to the beginning of each row of U in the matrix alu,jlu. |
iwk | integer. 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. |
ierr | integer 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. |
nz | integer pointer. Return number of nonzero entries for alu and jlu |
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.
n | row number of A |
a | nonzero entries of A |
ja | integer array of column for A |
ia | integer array of row pointers for A |
lfil | integer. 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. |
droptol | real*8. Sets the threshold for dropping small terms in the factorization. See below for details on dropping strategy. |
permtol | tolerance 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] |
mbloc | integer.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. |
alu | matrix 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. |
jlu | integer array of length n containing the pointers to the beginning of each row of U in the matrix alu,jlu. |
iperm | permutation arrays |
iwk | integer. 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. |
ierr | integer 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. |
nz | integer pointer. Return number of nonzero entries for alu and jlu |
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.
n | row number of A |
colind | integer array of column for A |
rwptr | integer array of row pointers for A |
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. |
nzlu | integer pointer. Return number of nonzero entries for alu and jlu |
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 | integer pointer. Return error message with the following meaning. 0 --> successful return. 1 --> not enough storage; check mneed. |
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)". *
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.
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
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.
ifix, float, min0, srtr