This preview shows page 1. Sign up to view the full content.
Unformatted text preview: oting is performed
and sparsity is not included. 52 Solution Techniques for Elliptic Problems for i = 2 to N do
for k = 1 to i ; 1 do aik = aik =akk
for j = k + 1 to N do
aij = aij ; aik akj end for
end for
end for Figure 9.3.9: Gaussian elimination loop.
The Gaussian elimination update in Figure 9.3.9 is aij := aij ; aik akj :
Let the current level of ll of aij be denoted as levij . Then the size of updated element
aij should be
levij ; levik levkj = levij ; (levik +levkj ) :
Roughly speaking, the size of the updated element aij will be the maximum of the two
sizes levij and (levik +levkj ) . Thus, it makes sense to de ne the new level of ll as levij := min(levij levik + levkj ):
It is common (cf. 4], Section 10.3) to shift all levels by 1 from this de nition. Thus, the
initial level of ll of aij is set as
aij 6= 0
levij := 0 ioftherwise or i = j :
(9.3.33a)
1
The updated level of ll is levij := min(levij levik + levkj + 1): (9.3.33b) We see that levij cannot increase during the factorization. Thus, if aij 6= 0 in the original
matrix A, levij never becomes nonzero. Equations (9.3.33) suggest a natural way of
discarding elements during the incomplete factorization. With an ILU(p) algorithm,
we'll keep all elements whose level of ll does not exceed p. The zero pattern of A for
the ILU(p) procedure may be de ned as the set Zp(A) := f(i j )jlevij > pg: (9.3.33c) 9.3. Conjugate Gradient Methods 53 This de nition is consistent with the ILU(0) preconditioning described in Example 9.3.5.
De ne levij according to (9.3.33a)
for i = 2 to N do
for k = 1 to i ; 1 do
if levik p then
aik = aik =akk
for j = k + 1 to N do
if levij p then
aij = aij ; aik akj
Update levij according to (9.3.33b) end if
end for
for j = k + 1 to N do
if levij > p then
aij = 0
levij = 1 end if
end for
end if
end for
end for Figure 9.3.10: Incomplete ILU(p) factorization.
A pseudocode ILU(p) algorithm appears in Figure 9.3.10. It is intended to give the
general idea of the procedure and, as stated, is not ve...
View
Full
Document
This document was uploaded on 03/16/2014 for the course CSCI 6840 at Rensselaer Polytechnic Institute.
 Spring '14
 JosephE.Flaherty

Click to edit the document details