SIAM J. SCI. COMPUT.
c
2007 Society for Industrial and Applied Mathematics
Vol. 29, No. 4, pp. 1476–1494
INCOMPLETE LU PRECONDITIONING WITH THE MULTILEVEL
FAST MULTIPOLE ALGORITHM FOR ELECTROMAGNETIC
SCATTERING
∗
TAH
˙
IR MALAS
†
AND LEVENT G
¨
UREL
†‡
Abstract. Iterative solution of large-scale scattering problems in computational electromagnet-
ics with the multilevel fast multipole algorithm (MLFMA) requires strong preconditioners, especially
for the electric-field integral equation (EFIE) formulation. Incomplete LU (ILU) preconditioners are
widely used and available in several solver packages. However, they lack robustness due to potential
instability problems. In this study, we consider various ILU-class preconditioners and investigate
the parameters that render them safely applicable to common surface integral formulations without
increasing the O(n log n) complexity of MLFMA. We conclude that the no-fill ILU(0) preconditioner
is an optimal choice for the combined-field integral equation (CFIE). For EFIE, we establish the need
to resort to methods depending on drop tolerance and apply pivoting for problems with high condi-
tion estimate. We propose a strategy for the selection of the parameters so that the preconditioner
can be used as a black-box method. Robustness and efficiency of the employed preconditioners are
demonstrated over several test problems.
Key words. preconditioning, incomplete LU preconditioners, multilevel fast multipole algo-
rithm, electromagnetic scattering
AMS subject classifications. 31A10, 65F10, 78A45, 78M05
DOI. 10.1137/060659107
1. Introduction. A popular approach in studying wave scattering phenomena
in computational electromagnetics (CEM) is to solve discretized surface integral equa-
tions, which give rise to large, dense, complex systems in the form of
A · x = b.Two
kinds of surface integral equations are commonly used. The electric-field integral
equation (EFIE) can be used for both open and closed geometries, but it results
in poorly conditioned systems, especially when the geometry is large in terms of the
wavelength. On the other hand, the combined-field integral equation (CFIE) produces
well-conditioned systems, but it is applicable to closed geometries only [21].
For the solution of such dense systems, direct methods based on Gaussian elimina-
tion are still widely used due to their robustness [34]. However, the large problem sizes
confronted in computational electromagnetics prohibit the use of these methods which
have O(n
2
) memory and O(n
3
) computational complexity for n unknowns. On the
other hand, by making use of the multilevel fast multipole algorithm (MLFMA) [12],
the dense matrix-vector products required at least once in each step of the iterative
solvers can be performed in O(n log n) time and by storing only the sparse near-field
matrix elements, rendering these solvers very attractive for large problems.
However, the iterative solver may not converge, or convergence may require too
many iterations. We need to have a suitable preconditioner to reach convergence in
∗
Received by the editors May 3, 2006; accepted for publication (in revised form) January 19, 2007;
published electronically June 21, 2007. This work was supported by the Scientific and Technical
Research Council of Turkey (TUBITAK) under research grant 105E172, by the Turkish Academy of
Sciences in the framework of the Young Scientist Award Program (LG/TUBA-GEBIP/2002-1-12),
and by contracts from ASELSAN and SSM.
http://www.siam.org/journals/sisc/29-4/65910.html
†
Department of Electrical and Electronics Engineering, Bilkent University, TR-06800, Bilkent,
Ankara, Turkey (tmalas@ee.bilkent.edu.tr, lgurel@bilkent.edu.tr).
‡
Computational Electromagnetics Research Center (BiLCEM), Bilkent University, TR-06800,
Bilkent, Ankara, Turkey.
1476
Downloaded 09/28/17 to 139.179.72.198. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php
INCOMPLETE LU PRECONDITIONING AND FAST MULTIPOLE 1477
a reasonable number of iterations and retain the O(n log n) complexity of MLFMA.
We can use the sparse near-field matrix
N to construct a preconditioner. Since this
matrix is the best available approximation to the coefficient matrix
A, it makes sense
to use the near-field matrix as a preconditioner and solve (for example) the left-
preconditioned system
(1.1)
N
−1
· A · x = N
−1
· b.
The inversion of the near-field matrix
N can be accomplished using direct meth-
ods which decompose the matrix into a product of a unit lower-triangular matrix
L and an upper-triangular matrix U. However, during the factorization of sparse
matrices, in general, fill-in occurs and the resulting factors lose their sparsity [28].
This may make it difficult to preserve the O(n log n) complexity of MLFMA. Nev-
ertheless, we can discard part of the fill-in and partially incorporate the robustness
of the LU factorization into the iterative method by using the incomplete factors of
N as a preconditioner. This is the general idea behind the incomplete LU (ILU)
preconditioners.
In a general setting, depending on the dropping strategy, we can talk about two
kinds of ILU-class preconditioners. The first one depends on the matrix structure and
the entries are dropped by their position. A “levels of fill-in” concept is introduced and
stronger preconditioners can be constructed by increasing the level of fill-in [31]. Since
this technique does not consider numerical values, it becomes ineffective in predicting
the locations of the largest entries, particularly for matrices that are far from being
diagonally dominant and indefinite [13]. This is the case for the matrices arising from
the EFIE formulation. Alternatively, one can drop the matrix elements depending on
their magnitudes, and the zero pattern is generated dynamically during the factor-
ization. Among such methods, ILUT(τ,p) proposed by Saad has been successful for
many general systems [3]. During the factorization, ILUT drops matrix elements that
are smaller than τ times the 2-norm of the current row; and of all the remaining en-
tries no more than the p largest ones are kept. ILUT is known to yield more accurate
factorizations than the level-of-fill methods with the same amount of fill-in [13].
Although ILUT is more robust than its counterparts depending on the level of
fill-in, it may occasionally encounter problems of instability for real-life problems.
Even when factorization terminates normally, the resulting incomplete factors may
sometimes be unstable. The common reasons of instability are in general excessive
dropping and small pivots [13]. If the problem is related to the small pivots, one
can significantly increase the quality of the ILUT preconditioner by using partial
pivoting as in the complete factorization case. The resulting preconditioner is called
ILUTP [31].
In order to understand the quality of the preconditioner, or to understand the
reason for failure when it occurs, we can use (
L · U)
−1
· e
∞
, where e is the vector
of ones. This statistic is called condest (for condition estimate) and it provides an
upper bound for (
L · U)
−1
∞
[13]. If the condest value is not very high, but the
preconditioner still does not work, one can deduce that fill-in should be increased to
achieve a successful preconditioner. On the other hand, if the condest value is high,
one can first try pivoting to remedy the situation instead of including more elements
in the incomplete factors.
Considering the remarkable success of ILU-class preconditioners for general non-
symmetric and indefinite systems [13] and the wide availability of ILU-class precondi-
tioners in various packages [2, 24, 6], the present study aims to develop a strategy for
Downloaded 09/28/17 to 139.179.72.198. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php
1478 TAH
˙
IR MALAS AND LEVENT G
¨
UREL
both selecting the most appropriate ILU-class preconditioner and determining their
parameters to use them as black-box preconditioners for CFIE and EFIE formula-
tions. We perform tests on canonical, quasi-canonical, and real-life problems with
increasing number of unknowns and show that when these preconditioners fail for the
reasons stated earlier, the failure can be circumvented using pivoting strategies with-
out increasing the memory cost. We also show that the condest value is very useful
for determining the quality of the resulting factorization before starting the iterative
solution.
ILU-class preconditioners have been tested for electromagnetic problems in [8,
32, 26]. In [8], ILU(0) was tried on systems resulting from EFIE formulation with
discouraging results in all test problems. Sertel and Volakis [32] tried ILU(0) on two
model problems. For the very small problem of 480 unknowns, ILU(0) was successful
with the EFIE and CFIE formulations, but clearly such a small problem is not repre-
sentative of large-scale CEM simulations. They presented only CFIE results for the
50,000-unknown problem; in this case, ILU(0) was quite successful in reducing the
number of iterations. Probably the most impressive results are those of Lee, Zhang,
and Lu [26], who tried the ILUT preconditioner on hybrid surface-volume integral
equations and showed it to be successful on many test problems. However, they nei-
ther tried commonly used EFIE or CFIE formulations nor applied pivoting or any
other techniques to increase the effectiveness of the preconditioner.
The rest of the paper is organized as follows. In the next section, we outline the
surface integral equations of CEM, their discretizations, resulting matrix equations,
and MLFMA. Spectral properties of such matrices are analyzed in section 3. Then,
in section 4, we comment on the preconditioning of the systems arising from integral-
equation formulations. In section 5, we present experimental results for the CFIE and
EFIE formulations. In the last section, we make some suggestions for the effective
and safe usage of ILU-class preconditioners with EFIE and CFIE.
2. Surface integral equations and fast solvers. EFIE and CFIE are surface
integral equations that are formed from the application of physical boundary condi-
tions; they are formulated to solve the radiation and scattering problems of arbitrarily
shaped geometries. This leads to the reduction of a three-dimensional problem into a
two-dimensional problem, but the resulting matrix equation becomes fully populated.
The boundary condition stating that the total tangential electric field should
vanish on a conducting surface can be mathematically expressed to obtain EFIE as
(2.1)
ˆ
t ·
S
dr
G(r, r
) · J(r
)=
i
kη
ˆ
t · E
inc
(r),
where E
inc
represents the incident electric field, S
is the surface of the object,
ˆ
t is
any tangential unit vector on S
, J(r
) is the unknown induced current residing on
the surface,
(2.2)
G(r, r
)=
I +
∇∇
k
2
g(r, r
)
is the dyadic Green’s function, and
(2.3) g(r, r
)=
e
ik|r−r
|
4π|r − r
|
is the scalar Green’s function for the three-dimensional scalar Helmholtz equation.
Green’s function represents the response at r due to a point source located at r
.
Downloaded 09/28/17 to 139.179.72.198. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php
INCOMPLETE LU PRECONDITIONING AND FAST MULTIPOLE 1479
Similarly, the boundary condition for the tangential magnetic field on a conduct-
ing surface is used to derive the magnetic-field integral equation (MFIE) [22, 14] as
(2.4) J(r) −
ˆ
n ×
S
dr
J(r
) ×∇
g(r, r
)=
ˆ
n × H
inc
(r),
where
ˆ
n is any unit normal on S
and H
inc
(r) is the incident magnetic field. It should
be noted that EFIE is valid for both open and closed geometries, whereas MFIE is
valid only for closed surfaces [18, 16].
Combining EFIE and MFIE, we obtain CFIE, i.e.,
(2.5) CFIE = αEFIE + (1 − α)MFIE,
where α is a parameter between 0 and 1. CFIE is free from the internal resonance
problems of both EFIE and MFIE and leads to well-conditioned systems [19]. On
the other hand, CFIE is not applicable to open geometries since it contains MFIE.
Therefore, CFIE is preferred over MFIE for closed geometries, but EFIE is the only
choice for open geometries. CFIE can also be extended for geometries containing both
closed and open surfaces [23].
Surface integral equations can be converted to linear systems of equations using
the method of moments (MOM). In fact, MOM can be used to reduce any linear-
operator equation to a matrix equation. Using a linear operator L, we can represent
EFIE and CFIE as
(2.6) L{f} = g,
where f is the unknown vector function, and g is the excitation. To discretize this
equation, f is first approximated by a set of vector basis functions via the expansion
(2.7) f ≈
n
j=1
x
j
f
j
,
where x
j
is the unknown coefficient of the jth basis function. Then this discretized
equation is tested at n points and a linear system is obtained, i.e.,
(2.8)
n
j=1
x
j
t
i
,L{f
j
} = t
i
, g,i=1,...,n.
In (2.8), t
i
is the vector testing function and the inner product is defined as
(2.9) t, f =
S
drt(r) · f (r).
By defining the elements of the coefficient matrix (or the so-called impedance
matrix) as
(2.10) A
ij
= t
i
,L{f
j
}
and the elements of the right-hand-side vector (or the excitation vector) as
(2.11) b
i
= t
i
, g,
Downloaded 09/28/17 to 139.179.72.198. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php
1480 TAH
˙
IR MALAS AND LEVENT G
¨
UREL
the linear system of equations can be expressed as
(2.12)
A · x = b.
In order to implement a Galerkin approach for discretization, we use Rao–Wilton–
Glisson (RWG) functions [29] for both basis and testing functions. RWG functions
are linearly varying vector functions defined on planar triangular domains. They are
widely used in MOM applications to discretize surface current distributions.
Surfaces of geometries are meshed with planar triangles in accordance with the
RWG functions. Each basis function is associated with an edge; hence the number
of unknowns n is equal to the total number of edges in a triangulation. For high
accuracy, triangle sizes have to be small, but this leads to problems with large numbers
of unknowns. As a rule of thumb, we choose the average size of the mesh about one-
tenth of the wavelength.
Since the integrodifferential operator L for both EFIE and CFIE involves long-
distance interactions, the resulting coefficient matrix is dense, which is expensive to
store and to solve. However, when we use iterative solvers, direct solution is replaced
with a series of matrix-vector multiplications, and we have the opportunity to perform
the matrix-vector product with O(n log n) complexity using MLFMA.
MLFMA can be described as the multilevel application of the fast multipole
method (FMM). An important concept, on which FMM relies, is the diagonalized
factorization of the Green’s function, i.e.,
(2.13) g(r, r
)=
e
ik|r−r
|
4π|r − r
|
=
e
ik|D+d|
4π|D + d|
≈
1
4π
d
2
ˆ
ke
i
ˆ
k·d
T
L
(k, D, ψ),
where D = |D| <d= |d|, the integration is on the unit sphere,
ˆ
k is the unit vector
normal to the unit sphere, and T
L
is the truncated sum
(2.14) T
L
(k, D, ψ)=
ik
4π
L
l=0
i
l
(2l +1)h
(1)
l
(kD)P
l
cos(ψ),
which is called the translation function [17]. In (2.14), h
(1)
l
(x) is the spherical Hankel
function of the first kind, P
l
is the Legendre polynomial, and ψ is the angle between
unit vectors
ˆ
k and
ˆ
D. The truncation number L is determined by the formula given
in [12],
(2.15) L ≈ kd +1.8d
2/3
0
(kd)
1/3
,
where d
0
is the number of accurate digits required in the matrix-vector multiplication.
In a physical setting, we can think of a matrix-vector multiplication as a set of
electromagnetic interactions between the basis and testing functions. When the basis
and testing functions are clustered according to the proximity of their locations in
space, the same translation function can be used for all interactions between pairs
of functions in any two clusters, instead of calculating the interactions separately.
For this purpose, radiation patterns of the basis functions f
j
weighted with the cor-
responding coefficients x
j
are evaluated on the unit sphere at K directions. Next,
in each basis group, the radiation patterns are added at each direction, a process
called aggregation. For each pair of basis and testing groups, a translation function
is defined. The overall radiation pattern of each basis group is multiplied by a trans-
lation function and thereby translated to the center of a testing cluster. Following
Downloaded 09/28/17 to 139.179.72.198. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php