A MULTIGRID-BASED SHIFTED-LAPLACIAN PRECONDITIONER
FOR A FOURTH-ORDER HELMHOLTZ DISCRETIZATION.
N.UMETANI∗, S.P.MACLACHLAN†, C.W. OOSTERLEE‡ §Abstract. In this paper, an iterative solution method for a fourth-order accurate discretizationof the Helmholtz equation is presented. The method is a generalization of that presented in [10],where multigrid was employed as a preconditioner for a Krylov subspace iterative method. This
multigrid preconditioner is based on the solution of a second Helmholtz operator with a complexvalued shift. In particular, we compare preconditioners based on a point-wise Jacobi smootherwith those using an ILU(0) smoother, we compare using the prolongation operator developed by deZeeuw in [37] with interpolation operators based on algebraic multigrid principles, and we compare
the performance of the Krylov subspace method Bi-CGSTAB with the recently introduced induced
dimension reduction method, IDR(s). These three improvements are combined to yield an efficientsolver for heterogeneous high-wavenumber problems.
Key words. Helmholtz equation, non-constant high wavenumber, complex-valued multigrid
preconditioner, algebraic multigrid, ILU smoother
AMS subject classifications. 65N55, 65F10, 78A45, 35J05

1. Introduction. Many authors, e.g. [5, 7, 13, 19], have contributed to the development of appropriate multigrid methods for the Helmholtz equation, but an efficient multigrid treatment of heterogeneous problems with high wavenumers arisingin engineering settings has not yet been proposed in the literature. The multigrid
method [4, 14] is known to be a highly efficient iterative method, for example, for discrete Poisson-type equations, even with fourth-order accurate discretizations [6, 33].
The Helmholtz equation, however, does not belong to the class of PDEs for which
off-the-shelf multigrid methods perform efficiently. Convergence degradation and,
consequently, loss of O(N ) complexity are caused by difficulties encountered in thesmoothing and coarse-grid correction components; see [7, 33] for a discussion.
We present an efficient numerical solution technique for the heterogeneous highwavenumber Helmholtz equation, discretized by fourth-order finite differences. Recently, in [10], a robust preconditioned Bi-CGSTAB method has been proposed forsolving these problems, in which the preconditioner is based on a second Helmholtz
equation with an imaginary shift. This preconditioner is a member of the familyof shifted Laplacian operators, introduced in [20], and its inverse can be efficiently
approximated by means of a multigrid iteration. Two-dimensional results, representative for geophysical applications, generated by second-order finite differences, havebeen presented in [26] and 3D results in [27].
In this paper, we generalize this solver and include, in particular, a fourth-order
discretization of the Helmholtz operator in our discussion. The multigrid preconditioner is enhanced, in the sense that we replace the point-wise Jacobi smoother inthe multigrid preconditioner by a variant of the incomplete lower-upper factorization smoother, ILU(0). Furthermore, we evaluate the performance of a prolongation∗ Delft University of Technology, Delft, the Netherlands. Exchange student from Tokyo University† Delft University of Technology, the Netherlands. Now at Tufts University, Medford, MA, US,
The work of this author was supported by the European Community’s Sixth Framework Programme,
through a Marie Curie International Incoming Fellowship, MIF1-CT-2006-021927.
‡ CWI, Center for Mathematics and Computer Science, Amsterdam, and Delft University ofTechnology, Delft, the Netherlands.
§ Corresponding author: c.w.oosterlee@cwi.nl

1

scheme that originates from algebraic multigrid (AMG) [28]. Next to this, we also
compare a recently introduced Krylov subspace method, called Induced Dimension
Reduction [32] (IDR), with the commonly used Bi-CGSTAB [34] method. It is expected that these enhancements to the iterative solver can reduce both the number of
iterations and the total CPU time needed for convergence. Moreover, we aim to reduce the size of the imaginary shift parameter in the shifted Laplacian preconditioner,so that an even faster solution method is obtained. Furthermore, the overall solution
method with these algorithmic improvements is not strictly restricted to structured
Cartesian grids, as it can be set up fully algebraically (a similar goal has been pursuedin [1]). Although our method extends to solving problems on unstructured grids, we
focus here on heterogeneous Helmholtz problems on Cartesian grids. We focus on
the two-dimensional case; however, all of the method’s ingredients can be directly
generalized to three dimensions.
This paper is set up as follows. In Section 2, we briefly discuss the 2D Helmholtz
equation, its field of application, and the discrete finite-difference formulations ofsecond and fourth order. The iterative solution method, including the preconditioner
and its components are presented in Section 3. Numerical results are presented in
Section 4, where we evaluate the iterative solver.
2. The Helmholtz Equation and Its Applications. Accurately imaging the
Earth is one of the major challenges in the hydrocarbon industry. Subsurface formations are mapped by measuring the time required for a seismic pulse to return tothe surface after reflection from interfaces between formations with different physicalproperties. Variations in these reflection times, as recorded on the Earth’s surface,usually indicate structural features in the strata below. Depths to reflecting interfaces
can be determined from the times, using velocity information that can be obtained
from the reflected signals themselves.
In geophysics, numerical methods can be employed in acoustic imaging techniques
to gain insight into the geological structures deep within the Earth’s subsurface [23,24].
Traditional mathematical imaging techniques based on the wave equation have beensuccessfully developed based on a high-frequency approximation. When the medium
is very complex (containing heterogeneities that result in strong lateral velocity variations), however, these so-called ray-based high-frequency migration techniques reach
their limits. The industry is gradually moving to finite-difference-wave-equation migration [23], as the linearized acoustic wave equation can be used to accurately image
steep dips. The numerical solution of the wave equation has greatly helped both forward modeling and migration of seismic wavefields in complex Earth structures, andit also serves as a starting point for solving the full inverse problem [25].
The wave equation can be solved in either the time or the frequency domain.
When applying the Fourier transformation with respect to time to the acoustic wave
equation, the Helmholtz equation is obtained. One advantage of this transformation
to the frequency domain is that the problem becomes easier to solve, because the time
variable, t, decouples. Instead of a single (coupled) equation in t, we get a series of
decoupled equations for each wavenumber, k, that can be solved in parallel. Oncea pressure solution in the frequency domain is obtained, the solution of the wave
equation can be computed efficiently by using the forward Fourier transformation. To
represent a solution of the wave equation for a given time interval accurately, however,a wide range of frequencies has to be resolved and the Helmholtz equation has to besolved independently for each frequency.
From the exploration-seismology point of view, the Earth is a heterogeneous semi2

infinite medium. The wavenumber can be large, which implies that the discretizedHelmholtz operator gives rise to both positive and negative eigenvalues and, therefore,
the discretization matrix, Ah , is indefinite. For 2D problems, however, the computation can be performed efficiently by using, for example, direct methods combined with
nested-dissection reordering [12]. Only one LU decomposition is needed to calculate
the solutions at multiple source locations. The result can be used for the computation of all of the wavefields, for all shots and, also, for the back-propagated receiverwavefields [24]. However, for 3D problems, the matrix sizes and bandwidths rapidly
become too large and one has to fall back on iterative methods. In that case, one nolonger has the advantages in the frequency domain related to the LU decomposition.
For the Helmholtz equation, unfortunately, many iterative methods suffer fromslow convergence, especially if high frequencies need to be resolved, due to the indefiniteness. The development of fast iterative methods for high-frequency Helmholtz
problems remains a subject of active research. One approach to iteratively solving
this equation is presented below.
2.1. Mathematical Problem Definition. We start with the description of the
2D Helmholtz problem which we would like to solve,−∇2 u(x, ω) − k(x)2 (1 − αi)u(x, ω) = g(x, ω), x ∈ Ω.

(2.1)

Unknown u(x, ω) represents the pressure field in the frequency domain, ∇2 is the
Laplacian operator, k(x) = ω/c(x) is the wavenumber, with c(x), the acoustic-wave
velocity, which varies with position, and ω = 2πf denotes angular frequency, a scalarmeasure of rotation rate (f is the frequency in Hertz). Wavenumber k depends on x
because of a spatially dependent speed of sound, c(x). The source term is denoted
by g. The medium is called barely attenuative√ if 0 ≤ α  1, with α indicating thefraction of damping in the medium (and i = −1, the imaginary unit). In geophysicalapplications, which are of our main interest, this damping can be up to 5% (α = 0.05).
While Equation (2.1) arises through the Fourier transform of a wave equation with avery simple model of damping, (−∇2 + (1 − αi)∂t2 )u = g, it is closely related to the
Fourier transform of the strongly damped wave equation, (−∇2 + τ ∂t ∇2 + ∂t2 )u = gk2that yields, after scaling, (−∇2 − 1+τ
ik )u. For small values of α = τ k, Equation (2.1)
is an accurate approximation of the Fourier-domain strongly damped wave equation.
The semi-infinite physical domain needs to be truncated for a numerical treatment. A popular approach in geophysics in order to obtain a satisfactory nearboundary solution, without many artificial reflections, is to use the absorbing boundary layer (ABL) approach; see, for example, [16] or [21]. This unphysical boundary
layer is used to gradually damp out the outgoing waves by adding dissipation in theequation outside the domain of interest. An efficient numerical technique should berobust with respect to this kind of feature. The absorption layers (denoted by Ωe )are attached to the physical domain, Ω, (see Figure 2.1). In Ωe , a damped Helmholtz
equation (2.1) should be satisfied [30], with
α = 0.25

k x − xd k2, x ∈ Ωe ,k xe − xd k2(2.2)

where point xd is a point at the boundary, Γ, and xe a point at Γe (see Figure 2.1).
At Γe , this equation is supplemented by second-order absorbing boundary conditions,as described in [2, 8].
3

Absorption layer
Γ
Ω

e

e

1000m/s

2000 m/sΩ1000 m/s

ΓPhysical domain

Fig. 2.1. A 2D domain with ABL in the case of a regular heterogeneous wedge medium.

2.2. Discretization. The equations are discretized either by a second- or afourth-order finite-difference scheme, resulting in the linear system:
Ah φh = bh ,(2.3)

where φh and bh represent the discrete frequency-domain pressure field and the source,
respectively.
In a heterogeneous medium, the smallest velocity is usually selected based on
the representative wavelength, λf . The number of wavelengths in a domain of size Lequals L/λf . A dimensionless wavenumber, k, on a non-dimensional [0, 1]2 domain isdefined by k = 2πf L/c, and a corresponding mesh size by h = λf /(nw L), with nw thenumber of points per wavelength. With domain size L = 1, an accuracy requirement,
for second-order discretizations, is that kh ≤ π/5(≈ 0.63) for nw = 10 points per
wavelength, and kh ≤ 0.53 with nw = 12 points per wavelength. The number ofgrid points used assumes a linear connection between k and h. In order to avoid areduction of accuracy due to pollution of the solution, however, k 2 h3 should be chosen
constant, as stated in [3, 15]. For an iterative solution method, the requirement thatkh should be constant is more severe and, so, this is the constraint that we consider
here.
For second-order finite differences, the 5-point discretization stencil used to set up
matrix Ah is well-known. For the absorbing boundary conditions at Γe , we also apply
central differences. The discretization that is particularly of interest in this work is
the O(h4 ) accurate discretization based on the Padé approximation. It is called the
HO discretization in [29], with stencil,

∧

AHO
=
h

− 16

1 
 2 (kh)2 (1 − αi)2h  −3 −121−6(kh)2 (1 − αi)1210 − 2(kh)2 (1 − αi)33(kh)2 (1 − αi)
2
−3 −
12
− 32 −− 16(kh)2 (1 − αi) 
2
.
−3 −
12

1
−6
(2.4)

An important reason for choosing a higher-order discretization method is that the
number of grid points per wavelength can be reduced compared to a second-order
discretization. This results in smaller matrices for the same level of accuracy and,4

thus, may lead to an algorithm that is more efficient overall, if the matrices associatedwith the higher-order discretization can be solved efficiently.
These matrices remain positive definite as long as k 2 is smaller than the first
eigenvalue of the discrete Laplacian. The wavenumber in geophysical applications
can, however, be large, which implies that the discretized Helmholtz equation gives
rise to both positive and negative eigenvalues and, therefore, the discretization matrix,
Ah , is indefinite. The size of the system of linear equations (2.3) gets very large forhigh frequencies. So, Ah in (2.3) is a large but sparse matrix, with complex-valued
entries, because of the absorbing boundary conditions and the attenuative medium.
It is symmetric but non-Hermitian.
2.2.1. Validation of the discretization. In order to validate the choice ofboundary condition, ABL, and discretization, we compute the solution for a constantwavenumber problem in a homogeneous medium with the source function, representative for a seismic pulse, chosen as
gh =

1
1
1
δ(x1 − , x2 − ).
h2232Here, δ(·, ·) represents the Dirac delta function, which is 1 when its argument is(0, 0), and 0 elsewhere. The scaling by h2 guarantees that the solutions on fine andcoarse grids are of the same amplitude, giving a discrete approximation of a δ-function
distribution.
Two formulations of the boundary discretization are compared here. In the first,
we prescribe the second-order absorbing boundary conditions directly at the physicalboundaries, whereas, in the second formulation, the boundary discretization is basedon an extra absorbing boundary layer, placed along all physical domain boundaries.
An ABL of n/4 points is added to each side.
In this first numerical experiment, we fix the wavenumber, k = 40, and use amodel domain, (0, 1)2 , covered by a fine grid consisting of 2562 points (kh = 0.156).
Figure 2.2 presents the two corresponding solutions with the second-order discretization. An unphysical damping of the solution without the ABL can be observed near
the domain boundaries.
Figure 2.3 presents solution profiles along the line x = 0.125, for the secondand fourth-order discretizations, with and without the ABL, on three meshes with

Fig. 2.2. Numerical solutions for k = 40 and h = 1/256, without (left-side) and with (rightside) ABL.
5

0.05

0.06
256 POINTS
128 POINTS
64 POINTS

0.04

256 POINTS
128 POINTS
64 POINTS
0.04

0.03
0.02

0.02

0.01
0

0

−0.01
−0.02

−0.02
−0.03

−0.04

−0.04
−0.05

−0.06
0

0.2

0.4

0.6

0.8

1

0

(a) No ABL, 2nd order

0.2

0.4

0.6

0.8

1

(b) With ABL, 2nd order

0.05

0.06
256 POINTS
128 POINTS
64 POINTS

0.04

256 POINTS
128 POINTS
64 POINTS
0.04

0.03
0.02

0.02

0.01
0

0

−0.01
−0.02

−0.02
−0.03

−0.04

−0.04
−0.05

−0.06
0

0.2

0.4

0.6

0.8

1

0

(c) No ABL, 4th order

0.2

0.4

0.6

0.8

1

(d) With ABL, 4th orderFig. 2.3. Comparison of the vertical line solutions, for three consecutive grid sizes, at x = 0.125,
k = 40, with second- and fourth-order discretizations; Left side: No ABL, right side: with ABL.
0.06

0.06
48 POINTS
64 POINTS

32 POINTS
48 POINTS
64 POINTS

0.04

0.04

0.02

0.02

0

0

−0.02

−0.02

−0.04

−0.04

−0.06

0

0.2

0.4

0.6

0.8

1

−0.06

0

0.2

0.4

0.6

0.8

1

Fig. 2.4. Comparison of the vertical line solutions on coarse grids, at x = 0.125, k = 40, left
side: second-order discretization, right side: fourth-order discretization.

642 , 1282 and 2562 points. Grid convergence for the three mesh sizes is achievedin both cases, but the solution profiles, with and without ABL, differ significantly,
especially near the domain boundaries. We also compare, for the same problem, thesolution profiles on coarse grids, with the second- and fourth-order discretizations.
The ABL is now attached to the domain. In Figure 2.4a, the profiles for the secondorder discretization with 642 and 482 interior points are presented. A non-physicalphase shift is observed on the grid that is too coarse (kh > 0.8) with this second-order
discretization. In comparison, Figure 2.4b presents the solution profiles with the HO
discretization on 322 , 482 and 642 grids, which converge nicely towards the physicalsolution. Thus, in this case, the fourth-order discretization with the ABL leads to an
accurate numerical solution, already on relatively coarse grids.
63. Iterative Solution Method. Iterative solution methods for complex-valued
indefinite systems based on Krylov subspace methods [31] are typically generalizationsof the conjugate-gradient (CG) method. The Bi-conjugate gradient stabilized (BiCGSTAB) algorithm [34] is one of the better known Krylov subspace algorithms for
non-Hermitian problems, which has been used for Helmholtz problems, for example,in [1, 10]. One of the advantages of Bi-CGSTAB, compared to full GMRES, is its
limited memory requirements.
Without a preconditioner, however, the Krylov subspace methods converge very
slowly, or not at all, for the problems of interest [9]. By preconditioning with a matrix,
Mh−1 , we solve an equivalent linear system,Ah Mh−1 φ̃h = bh ,φ̃h = Mh φh .

(3.1)

The challenge, then, is to find a matrix, Mh , such that Ah Mh−1 has a spectrum thatis favorable for iterative solution with Krylov subspace methods, and whose inverse,
Mh−1 , can be efficiently approximated.
In [10], a shifted-Laplacian operator was proposed as a preconditioner for theHelmholtz equation, with Mh defined as a discretization of
M = −∇2 − k 2 (x)(β1 − β2 i).

(3.2)

Boundary conditions were set identically to those for the original Helmholtz equation.
The influence of parameters β1 and β2 was evaluated in [10], and the optimal values forthe solver proposed there were (β1 , β2 ) = (1, 0.5). Here, we will also consider β2 = 0.4(smaller values of β2 do not lead to a converging algorithm with the components to
be introduced below). The matrix after discretization of (3.2), Mh , is obtained fromeither the 5−point, O(h2 ), or the 9-point, O(h4 ), finite-difference discretization.
3.1. Fourier Analysis. The discrete Helmholtz matrix, Ah , as well as the preconditioner, Mh , allow us, assuming a constant wavenumber and Dirichlet boundary
conditions, to apply Fourier analysis on the basis of discrete sine-eigenfunctions,
vhp,q = sin(pπx) sin(qπy),(3.3)

to gain insight into the spectrum of Ah Mh−1 . With these discrete sine functions,Ah Mh−1 is diagonalized, and the eigenvalues can easily be determined. As long as
k 2 is not equal to any of the eigenvalues of the discrete Laplace operator, Ah Mh−1 isnonsingular. Otherwise, the matrix is singular and its nullspace is spanned by thecorresponding eigenfunctions (3.3).
We perform Fourier analysis here to visualize the effect of the choice of the parameter, β2 , as well as the choice of discretization on the clustering of the eigenvalues ofthe preconditioned system. This analysis gives a first indication of what we can expectfrom the solver. For both Ah and Mh , we use either the second-order discretizationsor the fourth-order, HO stencils. We do not include damping in Ah in the analysis (wetake α = 0 in (2.1)), and concentrate here on the case k = 100 (k 2 = 104 ), h = 1/160.
Preconditioners with complex Helmholtz terms give rise to a curved spectrum, seealso [11]. Figure 3.1 presents the curved spectrum of Ah Mh−1 for (β1 , β2 ) = (1, 0.4) in
Mh , where both operators, Ah and Mh , are discretized by the fourth-order stencil. Avery similar eigenvalue distribution is obtained with the second-order discretization.
For β2 = 0.5, we find essentially the same spectrum; however, more eigenvalues arein the vicinity of the origin.
7

0.6

0.4

0.2

0

−0.2

−0.4

−0.6
−0.2

0

0.2

0.4

0.6

0.8

1

1.2

Fig. 3.1. Spectral picture of Mh−1 Ah with α = 0, (β1 , β2 ) = (1, 0.4), k = 100, h = 1/160. Bothoperators are discretized by the fourth-order, HO discretization.
0.04

0.04

2nd
4th

0.03
0.02

0.02

0.01

0.01

0

0

−0.01

−0.01

−0.02

−0.02

−0.03

−0.03

−0.04
−0.000500.00050.001

0.00150.002

0.0025β 2 = 0.5
β 2 = 0.4

0.03

0.003

−0.04
−0.0005(a) 2nd versus 4th order, β2 = 0.5

00.00050.001

0.00150.002

0.0025

0.003(b) 4th order, β2 = 0.4 vs. β2 = 0.5

Fig. 3.2. Comparison of zoomed spectral pictures of Mh−1 Ah with α = 0, k = 100, h = 1/160.
(a): Comparison of eigenvalues near the origin discretized with second-order and with fourth-order
discretizations; (b) Comparison of eigenvalues for β2 = 0.4 and β2 = 0.5, fourth-order discretizations.

These near-zero eigenvalues are problematic for the convergence of the Krylov
subspace method. Figure 3.2a shows a zoom of the spectrum near the origin, comparing (for β2 = 0.4) the location of the eigenvalues near the origin with the secondand the fourth-order discretizations. In Figure 3.2b, we also compare the locationof the eigenvalues near the origin for β2 = 0.4 and β2 = 0.5, with the fourth-order
discretization, keeping k = 100, h = 1/160. With the fourth-order discretization, theeigenvalues stay further from the origin as compared to the second-order discretization. This has a positive effect on the convergence of the Krylov subspace method.
The same is true when comparing the eigenvalues with β = 0.4 and β2 = 0.5 where,
as expected, the clustering with β2 = 0.4 is more favorable for iterative solution.
When discretized with second-order finite differences, Mh , with (β1 , β2 ) = (1, 0.5),
can be relatively easily handled by a multigrid method, which is confirmed by Local
Fourier Analysis, a quantitative multigrid analysis tool [33]. This is due to the imaginary term, β2 i in the shifted Laplacian. LFA also indicates that Mh based on the
fourth-order discretization can be dealt with in multigrid as efficiently as the secondorder discretization. With the robust multigrid components presented below, includ8

ing the ILU(0) smoother used within the preconditioner, we aim to decrease the valueof β2 , to β2 = 0.4, and obtain an efficient preconditioned Krylov subspace method.
3.2. Multigrid Preconditioner. One multigrid cycle, based on standard gridcoarsening and point-wise smoothing, can be used as an approximation to Mh−1 .
In [10], an F(1,1)-cycle [33], with one pre- and one post-smoothing iteration, with a
Jacobi smoother with under-relaxation parameter ω = 0.5 was chosen for the highwavenumber problems. The other multigrid components were:(i) Restriction operators, IhH , based on 2D full weighting, whose stencil [33] reads:

1
1
∧
 2
IhH =
16
1

2
4
2

H
1
2  ,
1 h

(3.4)

with h denoting the fine-, and H denoting the coarse-mesh size.
h(ii) Prolongation operators, IH, were 2D matrix-dependent interpolation, basedon de Zeeuw’s interpolation weights [37].
(iii) Coarse-grid matrices were based on Galerkin coarse-grid discretizations, deh.
fined as MH = IhH Mh IH
In [10], it was shown that the full-weighting restriction combined with the matrixdependent prolongation resulted in robust convergence for a variety of problems withirregular heterogeneities and strong contrasts. The inclusion of an ABL in the discretization does not lead to any multigrid convergence difficulties, as the components
chosen are especially designed for problems with varying coefficients.
With a more powerful smoother, we expect that a robust multigrid method can be
developed for approximately inverting matrices Mh that originate from a fourth-order
discretization. As the smoother in the multigrid preconditioner, we replace the pointwise Jacobi smoother by an ILU smoother, well-known in the multigrid literature [17,18, 35, 36, 38]. We choose here the ILU(0) variant, meaning that we do not allow
any additional fill-in in the lower- and upper-triangular factors outside of the nonzero
pattern of matrix Mh . An ILU(0) smoother is known to be more powerful than apoint-wise Jacobi smoother for a number of test problems [35]. Strictly speaking, ILU
methods do not only have a smoothing effect on the errors. A lexicographical version
may also reduce low-frequency errors, especially when the entries of the remainder
matrix, Rh , in
Mh = L̂h Ûh − Rh ,are relatively small.
3.3. AMG Type Interpolation. An efficient multigrid scheme relies on theeffective complementarity of the chosen relaxation and interpolation procedures in reducing the error components in an approximate solution. The coarse-grid correction
operator is designed to reduce errors that the chosen smoother is slow to attenuate.
Such errors should lie in the range of interpolation, so that the coarse-grid correction
may be effective. Here, we consider a fixed choice of coarse grid, using the same Cartesian coarse grids (doubling the mesh size in each direction) as in geometric multigrid
and the multigrid method examined in [10], but evaluate an interpolation operator
that is chosen based on algebraic multigrid (AMG) principles. The interpolation developed is largely based on the real-valued AMG interpolation from [28], and discussed
for complex-valued equations in [22].
9Consider, then, an error, eh , that is not quickly reduced by relaxation. For many
standard problems and smoothers, these errors coincide with those vectors that yieldsmall residuals. For the purpose of interpolation, AMG assumes that the error, eh , ismuch larger than its residual when measured point-wise, (Ah eh )j  (eh )j , for eachfine-grid index j. Based on this property, we haveX
(Ah eh )j ≈ 0 ⇒ ajj (eh )j ≈ −ajk (eh )k ,
(3.5)
k6=j

meaning that the value of the error at a fine-grid node, j, can be accurately approximated by the values from its neighboring nodes. If all neighboring nodes are alsocoarse-grid nodes, then (3.5) is easily turned into an interpolation formula.
With either the fixed coarsening considered here, or with any other sufficiently
rapid coarsening procedure, fine-grid node j will have both fine-grid and coarse-gridnodes as neighbors. Designing an interpolation procedure can, then, be thought of as
modifying the balance in (3.5) in such a way as to remove connections to other fine-gridneighbors of j while preserving the overall balance. This is typically done by applying
a partition to the neighboring nodes of j that identifies some nodes as important, or
strong, connections and other nodes as unimportant, or weak connections. That is,
we write the set, {k 6= j} = Cj ∪ Fjs ∪ Fjw , where Cj is the set of strongly connectedcoarse-grid neighbors of j, and the disjoint sets, Fjs and Fjw , denote the strong fine-gridand weak connections, respectively.
The matrix arising from the Helmholtz equation is complex and, typically, thesum of the moduli of the off-diagonal elements is larger than that of the diagonal
element in each row. In this case, a different criterion should be considered as ameasure of the strong connections. Here, we give two common criteria for defining
the set, Sj , of strong connections for node j, defining
Sj = k : |ajk | ≥ θ max |ajl | ,
l6=j

or

Sj =


k : Re(ajk ) ≥ θ max Re(ajl ) .
l6=jParameter θ allows some adjustment of the number of connections chosen as strong(relative to the strongest connection); for many problems, θ = 0.25 is considered to bea standard choice. Numerical experiments with the discrete complex-valued shifted
Laplacian have revealed that sometimes divergence is observed for high wavenumber
problems if we use the measure based on the norm. The measure based on the realpart of the matrix elements gave a robust multigrid performance over a large rangeof wavenumbers and, thus, is used in the numerical results that follow.
It is expected that the weak connections of fine-grid node j can be discarded
from the balance in (3.5). To remove these terms (in particular, the appearance of(eh )k for k ∈/ Sj ) without upsetting the balance, these terms are “lumped to the
diagonal”. In effect, this means that we make the approximation that (eh )k ≈ (eh )jfor k ∈ Fjw ; while this approximation may not be very accurate, it is not harmful
to make such a choice, since the connections involved are not important. Treating
the strongly connected fine-grid neighbors of j is much more important, as these areconnections that (by definition) cannot be easily dropped. In classical AMG methods,
one assumes that these connections are well-represented on the coarse grid, by their
10values at neighboring points. Then, an approximation may be made by considering
the weighted average of the values at common coarse-grid neighbors of node j and its
fine-grid neighbor, node k, resulting in the expression,P

akl (eh )l

l∈C

Pj

(eh )k ≈l∈Cj akl.

P
If there is no point in Cj such that akl 6= 0 (or ifl∈Cj akl = 0), then node k
is neglected in the interpolation formula. Making these substitutions in (3.5) andchoosing for equality, we then haveP
ajj (eh )j = −Xajk (eh )k −

l∈C

Pj

ajk

akl (eh )ll∈Cj akl

k∈Fis

k∈Cj

or (eh )j =

X

−

X

ajk (eh )j ,

k∈Fiw

P

k∈Cj wjk (eh )k , forajk +wjk = −

P

m∈Fis

ajj +

Pajm amkl∈C amljPm∈Fiw ajm.

With these weights, we can form the coarse-to-fine transfer matrix, W , from whichh, as
we can express the overall prolongation matrix, IH
h
IH
=



W
I


.

We stress that while we only investigate the use of this interpolation for structured
grids in this work, the use of these multigrid components enable the solution of
unstructured-grid Helmholtz problems, which will be the subject of future work.
3.4. Induced Dimension Reduction, IDR(s). Finally, we discuss an alternative to the Bi-CGSTAB iterative method. Recently, an efficient alternative to BiCGSTAB has been proposed in [32], also with limited memory requirements, called
Induced Dimension Reduction, IDR(s). The parameter, s, in IDR(s) determines thenumber of pre-defined vectors used to enhance the method’s convergence. In [32], it
has been shown that IDR(1) has similar cost in terms of the memory requirements andcomputational complexity as Bi-CGSTAB. With higher values of s, the algorithm’s
storage requirements increase, but, typically, its performance improves accordingly.
With s = 4, for example, 17 vectors need to be stored, whereas Bi-CGSTAB requires
only 7. Bi-CGSTAB is based on the idea of computing two mutually bi-orthogonalbases for the Krylov subspaces based on matrix, Ah , and its conjugate transpose, AHh ;IDR(s), on the other hand, is based on the generation of residuals that are forced to bein subspaces of decreasing dimension. The IDR(s) algorithm is as easy to implement
as Bi-CGSTAB and other Krylov methods. In pseudo-code, the algorithm reads:
11while ||rl || > TOL or l < MAXIT dofor k = 0 to s do
Solve c from P H dRl c = P H rlv = rl − dRl c; t = Ah v;
if k = 0 thenω = (tH v)/(tH t) ;endif
drl = −dRl c − ωt; dxl = −dXl c + ωv;
rl+1 = rl + drl ; xl+1 = xl + dxl ;l = l + 1;
dRl = (drl−1 , . . . , drl−s ); dXl = (dxl−1 , . . . , dxl−s );end forend whileHere, we make a comparison between Bi-CGSTAB and the IDR(s) method for theheterogeneous high-wavenumber problems with realistic absorbing boundary conditions.
4. Numerical Experiments. In this section, we perform several 2D numericalexperiments of increasing complexity. We start with the constant wavenumber problem, which serves as a benchmark for the algorithmic choices, after which we evaluate
the method’s performance for a Helmholtz problem with a wedge heterogeneity ,anda model of the Sigsbee field.
4.1. Homogeneous problem. The first numerical experiments are based on
the homogeneous Helmholtz problem on the square domain, (0, 1)2 , to gain insight
into the overall performance of the solvers. The pulse source, g, is located near the1
surface, at ( 12 , 32
) and is represented by the scaled delta function. We will evaluate
the following three multigrid preconditioners:
1 A multigrid V(1,1)-cycle with de Zeeuw’s prolongation operator, FW restriction and Jacobi smoothing with relaxation parameter ω = 0.5. This is thesolver from [10].
2 A multigrid V(0,1)-cycle with de Zeeuw’s prolongation operator, FW restriction and ILU post-smoothing.
3 A multigrid V(0,1)-cycle with AMG’s prolongation operator, FW restriction
and ILU post-smoothing.
These preconditioners are combined with the Bi-CGSTAB and IDR(s) Krylov subspace solvers. The value of β1 in the shifted Helmholtz preconditioner equals 1, β2 is
set to either 0.4 or 0.5.
4.1.1. Systematic comparison of the multigrid schemes. We start by evaluating the choices made for the multigrid components in the preconditioner. Thus, we
compare ILU(0) relaxation with the point-wise Jacobi smoother, and the prolongation operator of de Zeeuw [37] with AMG interpolation (based on the real part of theoperator elements). We first consider the fourth-order discrete preconditioner in (3.2)
and evaluate its solution using multigrid. Here, we choose a 2D homogeneous model,
without an ABL. Starting with a random initial guess, zero boundary conditions andright-hand side, we present the asymptotic multigrid convergence factor, ρ, measured
after 300 iterations.
For a given h, we fix the wavenumber, k, by requiring that kh = 0.625, meaning
there are approximately ten grid points per wavelength. Table 4.1 presents the multi12

multigrid
preconditioner
ω-Jacobi
Zeeuw-V(1,1)
ILU(0)
Zeeuw-V(0,1)
ILU(0)AMG-V(0,1)

β2
0.4
0.5
0.4
0.5
0.4
0.5

1/64
0.66
0.65
0.39
0.27
0.38
0.27

1/128
0.87
0.65
0.44
0.27
0.46
0.27

h:
1/256
0.90
0.65
0.60
0.30
0.56
0.30

1/512
0.91
0.65
0.77
0.33
0.68
0.34

Table 4.1
Multigrid convergence factors for solving Mh , the discrete (fourth-order discretization) versionof Equation (3.2).

grid convergence factors for two values of β2 (with β1 = 1) and varying mesh size h.
It shows that for β2 = 0.5, as in [10], a stable multigrid convergence for decreasing
mesh sizes (and increasing wavenumbers) is obtained. For β2 = 0.4, the multigrid
convergence rates increase; however, the combination of the V(0,1)-cycle, the ILU(0)smoother and the AMG prolongation operator still performs well for β2 = 0.4 andh = 1/512. The ILU(0) smoother shows a significant improvement in the multigrid
convergence, compared with those using a V(1,1)-cycle with an ω-Jacobi smoother.
Moreover, the multigrid convergence for the fourth-order discretization is similar tothat of the second-order discretization (not shown here). In fact, for the second-order
discretization and β2 = 0.4, we observed a multigrid divergence for the combinationof ILU(0) smoothing and AMG interpolation, whereas Table 4.1 presents a highly
satisfactory convergence (even better than that with the weights based on de Zeeuw’s
interpolation).
4.1.2. Second- and fourth-order discretizations. We compare the convergence of the Krylov subspace solvers for the second- and fourth-order discretizationsof both the original operator and the preconditioner. The ABL is also not included inthis experiment. We test each setting with a random initial guess and a point source
as a right-hand side. The iteration is terminated as soon as the relative residual isreduced to a prescribed tolerance of 10−6 ,||ri ||≤ 10−6 .
||r0 ||(4.1)

Tables 4.2 and 4.3 present, for fixed kh = 0.625, the Bi-CGSTAB performance onfour meshes, with the three multigrid preconditioners, for the second- and fourthorder discretizations, respectively.
For all of the solvers, we observe a linear increase in the number of iterations
for increasing wavenumbers. The performance of the V(0,1) multigrid preconditioner
with de Zeeuw’s prolongation weights and ILU(0) smoothing appears to be the mostrobust among these choices. Its convergence for both values of β2 , and for both
discretizations, is very satisfactory. An interesting observation, however, is that theperformance of the AMG interpolation (as in Table 4.3) is, especially on the finest
grid, significantly better for the fourth-order discretization. The CPU times per iteration reported for the fourth-order discretization are always somewhat higher thanfor the second-order problem, as we deal with 9-point discretization stencils on all
grids. However, as stated earlier, we do not need the extremely fine grids for highwavenumber problems in combination with the fourth-order schemes.
13h:preconditioner
ω-Jacobi
Zeeuw-V(1,1)
ILU(0)
Zeeuw-V(0,1)
ILU(0)AMG-V(0,1)

β2
0.4
0.5
0.4
0.5
0.4
0.5

1/64
39 (0.26)
36 (0.22)
17 (0.13)
19 (0.14)
26 (0.17)
27 (0.19)

1/128
75 (1.7)
69 (1.5)28 (0.63)
32 (0.70)
43 (0.94)
42 (0.92)

1/256
139 (12.6)
125 (11.3)
48 (4.3)
52 (4.7)
88 (7.9)
83 (7.4)

1/512
266 (99)
236 (88)
94 (36)98 (37)218 (83)
162 (61)

Table 4.2
Bi-CGSTAB performance for the homogeneous model (second-order discretization) in terms ofnumber of iterations and CPU time in seconds (in brackets).

multigridpreconditioner
ω-Jacobi
Zeeuw-V(1,1)
ILU(0)
Zeeuw-V(0,1)
ILU(0)AMG-V(0,1)

h:
β2
0.4
0.5
0.4
0.5
0.4
0.5

1/64
35 (0.30)
32 (0.25)
16 (0.17)
20 (0.19)
16 (0.13)
19 (0.16)

1/128
70 (2.0)62 (1.8)
26 (0.78)
30 (0.89)
26 (0.78)
30 (0.91)

1/256
122 (14.1)
110 (13.0)
46 (5.9)
51 (6.5)
45 (5.7)
52 (6.6)

1/512
215 (103)
201 (96)
85 (46)
96 (51)84 (45)95 (51)

Table 4.3
Bi-CGSTAB performance for the homogeneous model (fourth-order discretization) in terms ofnumber of iterations and CPU time in seconds (in brackets).

4.1.3. Comparison with IDR(s). In this subsection, we fix our discretization
to be fourth order and compare the convergence with and without ABL. The preconditioner chosen in all tests is the multigrid V(0,1)-cycle with ILU(0) smoothing
and AMG interpolation, that performed the best for the fourth-order problems tested
above.
Table 4.4 presents the number of matrix-vector products for preconditioned BiCGSTAB, IDR(2) and IDR(4), with the CPU times given in brackets. The presentation of the number of matrix-vector products, rather than the number of iterations,enables us to compare fairly among the different Krylov subspace methods. We choose
h = 1/512 and k = 320, so that kh = 0.625. For the computations with the ABL, we
add n/4 points to all sides of the domain. Thus, there is an increase in CPU time periteration, as seen in the table, comparing the results with and without ABL, whereas
the number of iterations is reduced with the ABL. Adding n/4 points in the ABL
increases the overall problem size by a factor of 2.25.
From the results in Table 4.4, we conclude that the performance of Bi-CGSTAB
in combination with the multigrid preconditioner chosen is optimal in terms of performance relative to the cost in memory. This may be somewhat surprising as the use of
the IDR(s) method is especially preferred when there is a significant gap in the number of iterations between Bi-CGSTAB and the full GMRES performance. IDR(s) isable to close this gap; i.e., with limited extra storage and work, IDR(s) may converge
like GMRES even when Bi-CGSTAB performs poorly. For the Helmholtz problems
evaluated here, however, the difference in the number of iterations between full GMRES and Bi-CGSTAB is less than 10 %, confirming the conclusion that Bi-CGSTAB
is favorable in these cases. This may be a result of the powerful preconditioner.
14

discretization

no
ABL

with
ABL

multigrid
prec.

β2

Krylov subspace method
Bi-CGSTAB
IDR(2)
IDR(4)

ILU
Zeeuw V(0,1)
ILU
AMG V(0.1)
ILU
Zeeuw V(0,1)
ILU
AMG V(0,1)

0.4
0.5
0.4
0.5
0.4
0.5
0.4
0.5

170 (46)
192 (51)
168(48)190(51)136 (84)
168 (104)134 (82)
166 (103)145 (45)
177 (54)
145(46)
176(56)122 (83)
149 (103)125 (85)
143 (98)

149 (50)
174 (59)
149(50)172(57)129 (93)154 (112)127 (92)
149 (108)

Table 4.4
Solver performance with Bi-CGSTAB and IDR(s) for the homogeneous problem, fourth-order
discretizations, with and without the ABL, in terms of the number of matrix-vector products andCPU time in seconds (in brackets).

number of points
n
32
64
128
256

β2
0.4
0.5
0.4
0.5
0.4
0.5
0.4
0.5

discretization
2nd order 4th order
32 (0.16) 22 (0.11)
31 (0.11) 27 (0.13)
20 (0.27) 17 (0.30)
20 (0.27) 20 (0.34)
17 (0.84) 17 (1.20)
21 (1.08) 19 (1.36)
18 (3.94) 17 (5.20)
21 (4.66) 19 (5.88)

Table 4.5
Bi-CGSTAB performance for second- and fourth-order discretizations for the homogeneousmodel in terms of number of iterations and CPU time in seconds (in brackets).

4.1.4. Fixed wavenumber, increasing mesh sizes. Finally, we reconsider
the 2D homogeneous model from Subsection 2.2.1, discretized on the unit square with
the ABL. Wavenumber k = 40 is set in this experiment. The number of grid points
increases in order to confirm grid independent convergence of the preconditioned BiCGSTAB solver for a fixed continuum problem. We terminate the iterations as soon
as the relative residual is less than 10−6 .
The iterative solver in this experiment is based on the Bi-CGSTAB method witha V(0,1) multigrid preconditioner, with β1 = 1, in which the ILU(0) smoother and theAMG prolongation are incorporated. Table 4.5 presents the number of Bi-CGSTAB
iterations, plus the CPU time to reach the termination criterion, for two values of β2 ,β2 = 0.4 and β2 = 0.5, with the second- and fourth-order discretizations. We observe
the h-independent convergence rate for the iterative solver; with fixed wavenumber
and h decreasing, approximately the same number of iterations is needed to satisfy
the termination criterion, for both discretizations. The convergence for the two discretizations is very similar.
4.2. The Wedge Problem. In this section, we present numerical results forthe so-called wedge problem. The domain, as in Figure 2.1, is a box, (0, 1)2 , in which
a wedge-shaped heterogeneity is placed, and the location of the source is (1/2, 1/32).
15

(a) k = 40(b) k = 80

Fig. 4.1. Numerical solutions for the wedge problem with k = 40 and k = 80 with ABL.

The wave number inside the wedge region is k, and outside the wedge it is set to k/2.
We employ the fourth-order discretization with the ABL here, with n/4 points on
both sides in the ABL.
Figure 4.1 presents the solutions of the wedge problem for k = 40 and k = 80.
Furthermore, for k = 80, we present in Figure 4.2 the vertical centerline solution profile
obtained on three different meshes. These profiles converge as expected for increasing
mesh sizes. With the fourth-order discretization, 64 points already seem sufficient foran accurate representation, whereas 128 points should be chosen according to the rulekh = 0.625.

Fig. 4.2. Centerline solution profile for the wedge problem with k = 80 with ABL, for different
mesh sizes.

Next, we examine the convergence of the preconditioned Bi-CGSTAB method,with the shifted Laplacian V(0,1)-multigrid preconditioner with the ILU(0) smootherand the AMG prolongation. Parameter kh is set to 0.625. Here, we also present
16

β2
0.4

0.5

damping
α
0.0%
1.0%
2.5%
5.0%
0.0%
1.0%
2.5%
5.0%

number of points (kh constant)
64
128
256
512
18 (0.31) 31 (2.16) 55 (16.4) 96 (117)
17 (0.34) 26 (1.84) 39 (11.8)
58 (72)15 (0.28) 21 (1.52) 28 (8.5)
44 (54)13 (0.28) 16 (1.22) 21 (6.5)
25 (32)23 (0.50) 38 (2.66) 68 (20.2) 121 (150)
21 (0.34) 31 (2.22) 46 (13.9)
71 (88)18 (0.34) 27 (1.92) 37 (11.6)
49 (61)15 (0.28) 19 (1.41) 26 (7.9)
28 (35)

Table 4.6
Bi-CGSTAB performance for the fourth-order discretization of the wedge model, with α %damping, in terms of number of iterations and CPU time in seconds (in brackets).

results with some damping included in the original Helmholtz equation. Parameterα in (2.1) varies between 0 and 0.05. The number of iterations and CPU time (inseconds) are presented in Table 4.6. We notice a significant improvement of themethod’s convergence, already when 1% damping is included in the original problem.
With 5% damping, we even observe a constant number of iterations, for varying
k. Compared to the performance of the solution method with the damped Jacobi
smoother (not shown here), the results in Table 4.6 are significantly improved, bothin terms of the number of iterations and in terms of the CPU time. Especially on thefine meshes, we see in Table 4.6 that it is beneficial to choose β2 = 0.4.
4.3. The Sigsbee Problem. The Sigsbee2A synthetic dataset models the geologic setting found on the Sigsbee escarpment in the deep-water Gulf of Mexico.
There is a substantial uniform layer of water at the top of the model. Here, we use ascaled version of the original Sigsbee model to test our iterative Helmholtz solver, see
Figure 4.3. The size of the domain is 150002 m and a source is placed at (7500, 117),
near the top wall. The frequency chosen for this computation is 5Hz.
The grid size consists of 5122 points with an ABL of 128 points on each edge. The
largest value of kh is 0.6135. Figure 4.4 presents the solution of this Helmholtz problem, where the fourth-order discretization is used. With the linear solver based on
preconditioned Bi-CGSTAB with the V(0,1) multigrid preconditioner for the shifted
Laplacian, using the ILU(0) smoother and AMG-based interpolation as essential components, we solve this problem in 61 iterations and 74.2 CPU seconds with β2 = 0.4,
and in 68 iterations and 85.5 CPU seconds for β2 = 0.5. This convergence is highly
satisfactory for this real-life setting. As a comparison, the original solver in [10] with
β2 = 0.5, the multigrid V(1,1)-cycle preconditioner, point-wise Jacobi smoothing andde Zeeuw’s interpolation needed 216 iterations and 237 seconds CPU time.
5. Conclusion. In this paper, we have presented a fast iterative solver for heterogeneous high-wavenumber Helmholtz problems. The solver is an improvement of
the work in [10]. A preconditioned Bi-CGSTAB solver has been developed where
the preconditioner is based on a shifted Laplacian with a complex-valued shift. We
have shown that it is possible to work with fourth-order finite differences, both in
the discrete original problem, as well as in the preconditioner without any problem.
An absorbing boundary layer improves the quality of the solution significantly anddoes not pose difficulties to the solution method proposed. We have focussed on
discretizations obeying a linear relation between the wavenumber and the mesh size,17

(a) Speed of sound, c(x) for Sigsbee
scaled domain(b) Distribution of k

Fig. 4.3. Domain for the scaled Sigsbee problem and the distribution of velocity, c(x), andwavenumber, k.

Fig. 4.4. Solution of the scaled Sigsbee problem with ABL, frequency 5Hz.

kh = 0.625, here.
The fourth-order accurate shifted Laplacian preconditioner can easily be approximated by one V(0,1)-cycle of multigrid. In the multigrid preconditioner, we haveincluded a powerful ILU(0) smoother and an AMG-based prolongation scheme. This
enables us to choose a somewhat smaller imaginary shift parameter (β2 = 0.4) in thepreconditioner, which improves the solver’s convergence (especially for high wavenumbers on fine meshes). The performance of the solver is significantly improved compared to the convergence results in [10]. With these additions, Helmholtz problems
on unstructured grids can, in principle, be solved in the same manner.
18

REFERENCES[1] T. Airaksinen, E. Heikkola, A. Pennanen, J. Toivanen, An algebraic multigrid based shiftedLaplacian preconditioner for the Helmholtz equation. J. Comp. Physics, 226: 1196–1210,
2007.
[2] A. Bamberger, P. Joly, J.E. Roberts, Second-order absorbing boundary conditions for the wave
equations: A solution for the corner problem. SIAM J. Numer. Anal., 27: 323–352, 1990.
[3] A. Bayliss, C.I. Goldstein and E. Turkel, On accuracy conditions for the numerical computationof waves. J. Comput. Phys., 59: 396–404, 1985.
[4] A. Brandt, Multi-level adaptive solutions to boundary-value problems. Math. Comput. 31: 333–
390, 1977.
[5] A. Brandt, I. Livshits, Wave-ray multigrid methods for standing wave equations. Elect. Trans.
Numer. Anal. 6: 162–181, 1997.
[6] W.L. Briggs, V.E. Henson, and S.F. McCormick, A multigrid tutorial: SIAM, Philadelphia,
USA, 2000.
[7] H. R. Elman, O. G. Ernst and D. P. O’Leary, A multigrid method enhanced by Krylov subspace
iteration for discrete Helmholtz equations, SIAM J. Sci. Comput., 23: 1291–1315, 2001.
[8] B. Engquist, A. Majda, Absorbing boundary conditions for the numerical simulation of waves,
Math. Comput., 31: 629–651, 1977.
[9] Y.A. Erlangga, C. Vuik, C.W. Oosterlee, Comparison of multigrid and incomplete LU shiftedLaplace preconditioners for the inhomogeneous Helmholtz equation. Applied Num. Math.
56: 648-666, 2006.
[10] Y.A. Erlangga, C.W. Oosterlee, C. Vuik, A novel multigrid based preconditioner for heterogeneous Helmholtz problems. SIAM J. Sci. Comput. 27: 1471-1492, 2006.
[11] M. B. van Gijzen, Y. A. Erlangga, C. Vuik, Spectral Analysis of the Discrete Helmholtz Operator Preconditioned with a Shifted Laplacian SIAM J. Sci. Comput., 29: 1942-1958,2007.
[12] A. George, J.W. Liu, Computer solution of large sparse positive definite systems, Prentice-Hall,
New Jersey, 1981.
[13] J. Gozani, A. Nachshon, E. Turkel, Conjugate gradient coupled with multigrid for an indefinite
problem, in Advances in Comput. Methods for PDEs V, 425–427, 1984.
[14] W. Hackbusch, Multi-grid methods and applications. Springer, Berlin, 1985.
[15] F. Ihlenburg and I. Babuska, Finite element solution to the Helmholtz equation with high wave
numbers. Comput. Math. Appl., 30: 9-37, 1995.
[16] C.-H. Jo, C. Shin, J.H. Suh, An optimal 9-point, finite-difference, frequency space, 2-D scalar
wave extrapolator, Geophysics 61(2): 529–537, 1996.
[17] R. Kettler, Analysis and comparison of relaxation schemes in robust multigrid and preconditioned conjugate gradient methods. In:W. Hackbusch and U. Trottenberg (eds.), Multigrid
methods, Lecture Notes in Mathematics 960: 502–534, Springer, Berlin, 1982.
[18] M. Khalil, Analysis of linear multigrid methods for elliptic differential equations with discontinuous and anisotropic coefficients. Ph.D. Thesis, Delft University of Technology, Delft,
Netherlands, 1989.
[19] S. Kim, S. Kim, Multigrid simulation for high-frequency solutions of the Helmholtz problem inheterogeneous media, SIAM J. Sci. Comput. 24: 684–701, 2002.
[20] A. L. Laird, M. B. Giles, Preconditioned iterative solution of the 2D Helmholtz equation. Report
NA 02-12, Comp. Lab., Oxford Univ., 2002.
[21] Q. Liao, G.A. McMechan, Multifrequency viscoacoustic modeling and inversion. Geophysics
61(5): 1371–1378, 1996.
[22] S.P. MacLachlan, C.W.Oosterlee, Algebraic multigrid solvers for complex-valued matrices. To
appear in SIAM J. Sci. Comput., 2008.
[23] W.A. Mulder, R.-E. Plessix, One-way and two-way wave-equation migration, Geophysics, 69:
1491–1504, 2004.
[24] W.A. Mulder, R.-R. Plessix, How to choose a subset of frequencies in frequency-domain finitedifference migration, Geophys. J. Int. 158: 801–812, 2004.
[25] R.G. Pratt, Seismic waveform inversion in frequency domain. Part I: theory and verification ina physical scale domain, Geophysics, 64: 888–901, 1999.
[26] C.D. Riyanti, Y.A. Erlangga, R.-E. Plessix, W.A. Mulder, C.W. Oosterlee, C. Vuik, A new
iterative solver for the time-harmonic wave equation, Geophysics 71: 57-63, 2006.
[27] C.D. Riyanti, A. Kononov, Y.A. Erlangga, C. Vuik, C.W. Oosterlee, R-E Plessix and W.A.
Mulder, A parallel multigrid-based preconditioner for the 3D heterogeneous high-frequency
Helmholtz equation. J. Comp. Physics 224: 431-448, 2007.
[28] J.W. Ruge and K. Stüben, Algebraic Multigrid (AMG). In: S.F. McCormick (ed.), Multigrid
19

Methods, Frontiers in Appl. Math., SIAM Philadelphia, 5: 73–130, 1987.
[29] I. Singer and E. Turkel, High Order Finite Difference Methods for the Helmholtz Equation,
Comp. Meth. Appl. Mech. Eng. 163:343-358, 1998.
[30] I. Singer, E. Turkel, A Perfectly Matched Layer for the Helmholtz Equation in a Semi-infinite
Strip. J. Comput. Phys. 201 (2004) 439-465.
[31] Y. Saad, Iterative Methods for Sparse Linear Systems. SIAM, Philadelphia, 2003[32] P. Sonneveld, M.B. van Gijzen, IDR(s): a new family of efficient algorithms for solving large
nonsymmetric linear systems. Delft Univ. Techn. Report 07-07, 2007.
[33] U. Trottenberg, C.W. Oosterlee and A. Schüller, Multigrid, Academic Press, London, 2001.
[34] H.A. van der Vorst, Bi-CGSTAB: a fast and smoothly converging variant of Bi-CG for thesolution of nonsymmetric linear systems, SIAM J. Sci. Comput., 13: 631-645, 1992.
[35] P. Wesseling, An introduction to multigrid methods. John Wiley, Chichester, 1992.
[36] G. Wittum, On the robustness of ILU smoothing. SIAM J. Sci. Comput. 10: 699–717, 1989.
[37] P. M. de Zeeuw, Matrix-dependent prolongations and restrictions in a blackbox multigrid solver,
J. Comput. Appl. Math., 33: 1–27, 1990.
[38] P.M. de Zeeuw, Incomplete line LU as smoother and as preconditioner. In: W. Hackbusch,
G. Wittum (eds.), Incompete decompositions (ILU) - algorithms, theory, and applications.
Vieweg, Braunschweig, 215–224, 1993. Int. J. Num. Methods in Fluids 20: 59–74, 1995.

20