{'fraction_non_alphanumeric': 0.05382607708892839, 'fraction_numerical': 0.018458977947573476, 'mean_word_length': 3.782561505065123, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 8, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 22, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Multilevel Preconditioners for Reaction-Diffusion Problems with Discontinuous Coefficients

In this paper, we extend some of the multilevel convergence results obtained by Xu and Zhu in [Xu and Zhu, M3AS 2008], to the case of second order linear reaction-diffusion equations. Specifically, we consider the multilevel preconditioners for solving the linear systems arising from the linear finite element approximation of the problem, where both diffusion and reaction coefficients are piecewise-constant functions. We discuss in detail the influence of both the discontinuous reaction and diffusion coefficients to the performance of the classical BPX and multigrid V-cycle preconditioners.


INTRODUCTION
In this paper, we will discuss the convergence of multilevel preconditioners for the linear finite element approximation of the second order elliptic boundary value problem with discontinuous coefficients: where Ω ∈ R d (d = 2 or 3) is a polygonal or polyhedral domain with Dirichlet boundary Γ D and Neumann boundary Γ N . While such problems arise in a wide variety of practical applications, our interest in (1.1) is motivated by the subspace problems in auxiliaryspace preconditioners for the definite Maxwell equations [17,18]. Multigrid algorithms are a family of powerful solution techniques which are frequently applied to the finite element discretizations of (1.1). When the coefficients ω > 0 and ρ ≥ 0 are constant, it is well known that Multigrid is an efficient optimal solver; while its additive version, the BPX algorithm, is an optimal preconditioner (see for example [5,15].) In many practical applications, however, the coefficients ω and ρ of (1.1) describe material properties, which can be considered constant in the material subdomains, but may have large jumps on the material interfaces.
There have been a lot of works devoted to developing efficient iterative solvers for solving the finite element discretization of (1.1), which are robust with respect to the jumps in the diffusion coefficient ω (when ρ ≡ 0), (see [6,28,29,10,22,14] for examples). For general cases, one usually need some special techniques to obtain robust iterative methods, (cf. [7,24,13,1]). Recently, Xu and Zhu addressed in [33,35] the performance of the BPX and Multigrid V-cycle preconditioners for (1.1) in the case of discontinuous ω and ρ = 0. It was shown that the jumps in ω affect only a small number of eigenvalues, and therefore the (asymptotic) convergence rate of the preconditioned conjugate gradient (PCG) method is uniform with respect to the jumps and the mesh size. See also [11,25,4,8,36] and the references cited therein for further developments in different directions.
All the analysis mentioned above focused on pure diffusion equation, and very little attention has been paid for the case when ρ is nonzero. In many applications, such as time discretization of heat conduction in composite materials, the equation involves a lower order term. In the case that ω i = ρ i , a robust overlapping domain decomposition method was developed for a two-dimensional model problem in [9]. Recently, the paper [19] discussed the performance of the algebraic multilevel iteration (AMLI) methods for the finite element discretizations for (1.1) in 2D, which are based on a multilevel block factorization and polynomial stabilization.
In this paper, we study the performance of the classical multilevel preconditioners (BPX and multigrid V-cycle) on the finite element discretization of equation (1.1), with emphasis on the discussion of the influence of both the discontinuous reaction and diffusion coefficients on the convergence of these multilevel preconditioners. We classify the coefficients in two different cases. In the first case, we require that both ω and ρ have the same coefficient distribution, namely, if ω i ≥ ω j then ρ i ≥ ρ j and vice versa. Note that this includes the case when ρ is a global constant. In this case, we recover the results from [33]. On the other hand, when ω and ρ have different distributions, it seems that the performance of the preconditioners deteriorate with the jumps (see the numerical examples in Section 5.3). In this case, we showed that the convergence rate depends on the minimal of the jumps in ω and ρ. As a special case, when ω is a global constant, or only varies moderately in the whole domain, we show that the multilevel preconditioners are robust with respect to both coefficients, and the mesh size.
The remainder of the paper is organized as follows: in the next section we investigate the Jacobi and Gauss-Seidel preconditioners. Then, in Section 3, we consider an interpolation operator which is needed in the analysis of the BPX and Multigrid algorithms, carried out in the following Section 4. The developed theory is illustrated by several numerical experiments collected in Section 5.
Throughout the paper we use the standard notation for Sobolev spaces and their norms. We will use the notation x 1 y 1 , and x 2 y 2 , whenever there exist constants C 1 , C 2 independent of the mesh size h and the coefficients ω and ρ, and such that x 1 ≤ C 1 y 1 and x 2 ≥ C 2 y 2 , respectively. We also use the notation x y for C 1 x ≤ y ≤ C 2 x.

NOTATION AND PRELIMINARIES
In this section, we establish the notations and review a few preliminary tools that will be needed for the subsequent analysis, following those in [33]. We consider solving the model equation (1.1) in a polyhedral domain Ω ⊂ R d for d = 2 or 3, and assume that there is a set of non-overlapping subdomains {Ω m } M m=1 such that Ω = ∪ M m=1 Ω m , on which the diffusion coefficient ω(x) and the reaction coefficient ρ(x) are constants, denoted by ω m := ω(x)| Ωm and ρ m := ρ(x)| Ωm for each m = 1, · · · , M respectively.
Let V = H 1 D (Ω) be the space that consists of H 1 (Ω) functions with vanishing trace on the Dirichlet boundary Γ D ⊂ ∂Ω. The variational problem of (1.1) reads: Given f ∈ L 2 (Ω) and g N ∈ H 1/2 (Γ N ), find u ∈ V such that where the bilinear form a(·, ·) is given by Ωm For the analysis of this paper, we will need the following weighted semi-norms and norms. Given any piecewise-constant coefficient τ = {τ 1 , · · · , τ M } > 0, we define weighted L 2 norm and H 1 semi-norm by In this notation, the bilinear form of interest is a(u, u) = |u| 2 1,ω + u 2 0,ρ . With a little abuse of the notation, we will use the same notation for the case when ρ = 0 in some subdomains of Ω, although in this case · 0,ρ is not a norm.
Let T h be a quasi-uniform triangulation of Ω. We assume that all Ω m are of unit size, and that their geometries are resolved exactly by the triangulation. Let V h ⊂ V be the corresponding linear Lagrangian finite element spaces. Then the finite element discretization of (1.1) reads: We define a linear symmetric positive definite (SPD) operator A : and use the notation · A = a(·, ·) to denote the energy norm. We need to solve the following operator equation, Since A is SPD, by the classical PCG theory we know that the convergence rate of the iterative method for A with a preconditioner, say B, is determined, by the (generalized) condition number of the preconditioned system: κ(BA) := λ N (BA)/λ 1 (BA), where λ i (BA) for i = 1, · · · , N are the eigenvalues of BA satisfying λ 1 ≤ λ 2 ≤ · · · ≤ λ N . However, if we know a priori that the spectrum σ(BA) of BA satisfies σ(BA) = σ 0 (BA) ∪ σ 1 (BA), where σ 0 (BA) = {λ 1 , . . . , λ m } contains all extreme ("bad") eigenvalues, and the remaining eigenvalues contained in σ 1 (BA) = {λ m+1 , . . . , λ N } are bounded from above and below, i.e., λ j ∈ [α, β] for j = m + 1, . . . , N , then the error at the k-th iteration of the PCG algorithm can be bounded by (cf. e.g. [2,16,3]): Specifically, if the number of extreme eigenvalues m is small, then the asymptotic convergence rate of the resulting PCG method will be determined by the ratio (β/α), which is the so-called effective condition number (cf. [21,33]).
Definition 2.1. Let T : V h → V h be a symmetric positive definite linear operator, with eigenvalues 0 < λ 1 ≤ · · · ≤ λ N . For m = 0, 1, · · · N − 1, the m-th effective condition number of T is defined by Remark 2.2. To estimate the effective condition number, and specifically λ m+1 (T ), a standard tool is the min-max principle (see e.g. [12, Theorem 8.1.2]): In particular, for any subspace As in [33], we introduce a subspace H 1 D (Ω) ⊂ H 1 D (Ω) as: We shall emphasize that m 0 is a fixed number which depends only on the distribution of the coefficient ω on the domain. We conclude this section by a discussion on the simple (but commonly used) Jacobi and Gauss-Seidel preconditioners. Note that the jumps in ρ do not influence the condition number estimates. Theorem 2.3 (cf. Theorem 2.2 in [33]). Let A be the stiffness matrix corresponding to a(·, ·) in V h , and let D be its diagonal. The condition number of D −1 A (Jacobi preconditioning) depends on the mesh size and the coefficient ω: On the other hand, the m 0 -th effective condition number is independent of the coefficients ω and ρ: Here m 0 = |I|, is the number of interior subdomains Ω m .
Proof. For ease of presentation, we introduce a mesh dependent coefficient defined as Note that when ρ = 0, we have ω h = ω as in [33]. Given any v h ∈ V h , let v be its vector representation in the nodal basis of V h .

First of all, by inverse inequality
This inequality implies that λ max (D −1 A) 1. On the other hand, by Poincaré-Friedrichs inequality, we have Since h 2 1, we obtain Then by the min-max principle (cf. Remark 2.2), we obtain that Thus we obtain the desired estimate for the m 0 -th effective condition number κ m 0 (D −1 A).
Remark 2.4. Analogous results hold for symmetric Gauss-Seidel preconditioner based on certain spectral equivalence between Jacobi and Gauss-Seidel iterations for SPD matrices (see [27] for more details).

INTERPOLATION OPERATOR
The analysis of the multilevel preconditioner relies on the approximation and stability of certain interpolation operator. In this section we describe the dual basis-based interpolation operator from [26], and show how it can be used to derive simultaneous estimates in two different weighted inner products.
Let T ∈ T h be a fixed mesh element and {λ T,i } be the set of its linear finite element shape functions. The local mass matrix on T has entries and it is easy to check that M T is spectrally equivalent to a diagonal matrix: Given v ∈ L 2 (Ω) we define Π h v ∈ V h by specifying its values in the vertices of T h . Specifically, the value at a vertex x is determined using an associated element T x ∈ T h : and for any w ∈ span {λ Tx,i }. Then Π h can be equivalently defined by We remark that the choice of T x is not unique. Given a particular ordering of the subdomains, say Ω 1 , · · · , Ω M , we choose T x ⊂ Ω k where k is the minimal index of all the subdomains that contain x. Note that this ordering has nothing to do with the actual geometry distribution of the coefficients. In order to make Π h satisfy certain stability property in the weighted norms, we may label the subdomains such that ω 1 ≥ ω 2 ≥ · · · ≥ ω M , or ρ 1 ≥ ρ 2 ≥ · · · ≥ ρ M . In this case, the choice of T x guarantees that the coefficient in T x is the maximum of all the coefficients in the neighborhood of x. By a standard argument, we have the following result on Π h .
In the last step, we used the property (3.1) on µ x i . Notice that T x i ⊂ Ω k where k is the minimal index of all the subdomains that intersect at x i . Therefore, the L 2 stability of Π h (3.3) follows by summing up all the elements in Ω m on both sides.
Based on this lemma, we have the following corollary on the stability of Π h in the weighted L 2 norms.
where k is the minimal index of all the subdomains that intersect at x.
(1) Π h is stable in the standard L 2 norm: (3) In the worst scenario, for any piecewise-constant coefficient τ > 0, we have where J (τ ) is the measure of the variation of τ defined by (2.4).
Remark 3.4. Other projection operators that are stable in both the L 2 and weighted L 2 norms are also available. For example, let The fact that it is also stable in the ρ-weighted L 2 norm follows from Now we turn to study the approximation and stability properties in the ω-weighted norms. It is standard (cf. [26]) that the Π h : H 1 D (Ω) → V h has the following classical approximation and stability estimates: In the ω-weighted norms, we have the following approximation and stability estimates.
where in the last step we used Lemma 3.2 for the stability of Π h on the subdomain Ω m .
Similarly, to show the weighted H 1 stability (3.10), we have on each subdomain Ω m : Then, the inequality (3.10) follows from the approximation and stability estimates of Q ω h (see for example [33,Lemma 3.3]). This completes the proof.

MULTILEVEL PRECONDITIONERS
In this section, we present the BPX and multigrid V-cycle preconditioners based on the subspace correction methods [31,34]. We present the main results of the robustness of these preconditioners with respect to the jump in the coefficients.
Let T 0 be an initial conforming mesh which resolves the jump interfaces. We obtain a nested sequence of triangulation {T k } L k=0 by a uniform refinement. Let h k be the mesh size of T k for k = 0, · · · , L, then we have h k γ k h 0 for some γ ∈ (0, 1), and L | log h L |. For simplicity, we denote h L = h. On each triangulation T k , let V k be the corresponding finite element space over T k . Then we obtain a sequence of nested spaces: and simply denote A = A L . A key ingredient in analyzing the multilevel preconditioners is the stable decomposition derived below.

Stable Decomposition.
With the help of the properties of the interpolation operator Π h , we now show several stable results of the subspace decomposition described above. In the multilevel context, we will use the notation Π k := Π h k . Also, we notice that Π L | V h = Id, i.e., the restriction of Π L on the finite element space V is identity. In particular, for any where Below, we discuss the stability of this decomposition in terms of the energy norm a(·, ·) = (·, ·) A , which involves both the ω-weighted H 1 semi-norm and the ρ-weighted L 2 norm. First, we consider the stability in terms of the ω-weighted H 1 semi-norm.
Lemma 4.1. The decomposition (4.1) satisfies the following properties: Proof. Given any v ∈ V h , to show (4.2), we notice that We have |Π 0 v| H 1 (Ω) ≤ |v| H 1 (Ω) . To estimate the second term on the right hand side of above inequality, we use the fact that Π k is stable in L 2 and Π k Q k = Q k , where Q k is the standard L 2 -projection on V k for k = 0, 1, · · · , L. Specifically, we have The estimate of the last sum is classical in the BPX theory, see [5,22]. Therefore, we have To prove (4.3) for any v ∈ V h , by the approximation and stability estimates (3.9)-(3.10) of Π k (k = 0, 1, · · · , L) in Lemma 3.5 and triangle inequality, we obtain This proves the inequality (4.3).
Now we establish the stable decomposition in the ρ-weighted L 2 norms. We have the following result. (1) For any u ∈ V h , the decomposition u 0 = Π 0 u and u k = (Π k − Π k−1 )u ∈ V k for k = 1, · · · , L satisfies: (2) If the reaction coefficients satisfy ρ 1 ≥ ρ 2 ≥ · · · ≥ ρ M , then for any u ∈ V h , the decomposition u 0 = Π 0 u and u k = (Π k − Π k−1 )u ∈ V k for k = 1, · · · , L is stable: Proof. Below, we only give the detailed proof of (4.5). The proof of (4.4) can be reduced to show the estimate which is a special case of (4.5). Given any u ∈ V h , by the stability (3.5) of Π 0 in the ρ-weighted L 2 norm, we have To estimate the summation term, we let Q ρ k be the ρ-weighted L 2 -projection on V k . Note that By the stability of (3.5) of Π k (k = 0, 1, · · · , L) in the ρ-weighted L 2 norm, Π k Q ρ k = Q ρ k and triangle inequality, we obtain To bound the sum on the right, we need to introduce some additional notation. Let V k be the space of discontinuous piecewise linear polynomials, associated with the same mesh as V k , and let Q k be the piecewise local L 2 -projection onto V k . Note that by defi- Here we used the (local) approximation property of Q ρ k , see [6,Lemma 4.5]:

This implies
Thus, we obtain This completes the proof.

4.2.
BPX Preconditioning. Now we are in position to discuss the performance of the multilevel preconditioners. For simplicity, we introduce the mesh dependent coefficient on each level k = 1, · · · , L ω k = ω + h 2 k ρ . (4.6) On each level k = 1, · · · , L, let R k be a smoother based on the SPD operator A k , and let R 0 = A −1 0 . We assume that the smoothers satisfy (R k u k , u k ) h 2 k u k 2 0,ω k ∀u k ∈ V k , k = 1, · · · , L which holds for the Jacobi and Gauss-Seidel smoothers, as shown in Section 2. Then the BPX preconditioner B : V → V is defined by where Q k : V h → V k be the L 2 projection on V k for k = 0, 1, · · · , L. The BPX preconditioner B satisfies the following well-known identity ( [30,32,34]): Based on the assumption on R k , it satisfies To analyze the BPX preconditioner, we make use of the following strengthened Cauchy Schwarz inequality. . For j, k = 0, · · · , L and j ≤ k, we have When ρ = 0, the largest eigenvalue of the matrix BA is bounded by a constant, as shown for example in [33]. For general ρ we only get a sub-optimal estimate. Lemma 4.4. The largest eigenvalue of BA is independent of the coefficients ρ and ω, but depends logarithmically on the mesh size: Proof. Given any u ∈ V h , let u = L k=0 u k with u k ∈ V k for (k = 0, · · · , L) be an arbitrary decomposition of u. Then by the Strengthened Cauchy Schwarz inequality (4.8), we obtain On the other hand, by the Schwarz inequality, we obtain Since the decomposition is arbitrary, we have (Au, u) |log h| (B −1 u, u), which completes the proof.
Remark 4.5. The estimate (4.9) can not be improved in general. This can be seen by taking u ∈ V 1 and decomposing it using u k = 1 L u for 1 ≤ k ≤ L. To estimate the smallest eigenvalue, we classify the coefficients in two different cases: (C1) The coefficients ω and ρ have the same distribution. Namely, if ω i ≥ ω j then ρ i ≥ ρ j and vice versa. (C2) The coefficients ω and ρ have different distribution.
In the case of (C1), we may label the subdomains based on the ordering ω 1 ≥ ω 2 ≥ · · · ≥ ω M . By the definition of Π h , it satisfies simultaneously the stable decomposition (4.5) in the ρ-weighted L 2 norm, and the stable decomposition (4.3) in the ω-weighted H 1 semi-norm. Based on these properties, we have the following result. Proof. For any u ∈ V h , we consider the decomposition Now we turn to discuss the case (C2) when ω and ρ have different distribution, e.g., there exists at least a pair of (neighboring) subdomains Ω i and Ω j on which ω i > ω j but ρ i < ρ j . In this case, the interpolation operator Π h does not satisfy the simultaneous stability in both the ρ-weighted L 2 norm and ω-weighted H 1 semi norm. Therefore, in this case, we only get some pessimistic estimates which depend on the jumps in the coefficients.
When J (ω) < J (ρ), we should label the subdomains based on the order of ρ to guarantee the ρ-weighted L 2 stability (4.5). Note that this includes the case when ρ = 0 in some subdomains, but not globally 0. While for the ω-weighted approximation and stability estimate, we apply the decomposition (4.2) instead of (4.3).
On the other hand, if J (ω) > J (ρ), then we should label the subdomains based on the ordering of ω to guarantee the stable decomposition (4.3) in the ω-weighted H 1 seminorm. For the stable decomposition in term of ρ-weighted L 2 norm, we can not apply (4.5) directly, but we may use the estimate (4.4). So in this case, we have the following result.
In summary, we have the following results for the BPX preconditioner.  In particular, if ω is a global constant, then the condition number of BA is independent of the jumps in both of ω and ρ.

4.3.
Multigrid V-cycle. Now we consider the Multigrid V-cycle as a solution algorithm and as a preconditioner to our original elliptic problem (1.1). We first introduce some standard notation. For each level k = 0, 1, . . . , L, we define the projections P k : V h → V k by a(P k u, v k ) = a(u, v k ) ∀v k ∈ V k . At each level, let R k : V k → V k be the smoothing operator. Here we use point Gauss-Seidel as the smoother. Then that standard multigrid V-cycle algorithm solves (2.1) by the iterative method , where the operator B k : V k → V k is defined recursively as follows: 0 , for k > 0 and g ∈ V k , define (1) Presmoothing : w 1 = R k g; (2) Correction: . We denote B L = B for simplicity. Following the same analysis in [33], it is clear that λ max (BA) ≤ 1. To estimate the smallest eigenvalue of BA, we consider the error propagation operator I − BA. By the XZ-identity (cf. [34]), we can get the following estimate, which is a straightforward generalization of [33,Lemma 5.2].
If we restrict to the subspace V h , we have a similar estimate: From Lemma 4.10, as in [33] we can deduce by min-max principle (cf. Remark 2.2): where m 0 = |I| is the number of floating subdomains. According to the above result, the convergence of the multigrid V-cycle method, and the condition number estimate of the multigrid preconditioner rely on the estimate on the constant c 0 ; while the estimate on the effective condition number relies on the estimate onc 0 . Both of these estimates follow from the stable decompositions (4.10) and (4.11). Now, based on the discussion for the BPX preconditioner case, we can obtain similar results for the multigrid V-cycle. Here m 0 = |I|, is the number of floating subdomains. (3) If the coefficients ω and ρ satisfy (C2), with J (ω) ≤ J (ρ), the condition number of BA is independent of the jump in ρ: In particular, if ω is a global constant, then the condition number of BA is independent of the jumps in both of ω and ρ.

NUMERICAL EXPERIMENTS
This section contains a set of numerical experiments performed with a version of the finite element library MFEM [20], which illustrate the convergence theory developed in the preceding sections. We focus on the commonly used V (1, 1)-cycle Multigrid method, and use a symmetric Gauss-Seidel iteration as a smoother. The same smoother was also used in the BPX algorithm, whose optimal implementation can be found in [5].
The jump-independent estimates of the effective condition number in Section 4.2 and Section 4.3 imply that a preconditioned conjugate gradient (PCG) acceleration will result in a solver which is optimal with respect top the mesh size. To investigate this, we report the number of PCG iterations needed to reduce the relative residual by a factor of 10 −12 .  We use the abbreviations GS-CG, BPX-CG and MG-CG to denote the symmetric Gauss-Seidel, BPX and Multigrid preconditioners respectively.
We run a simple test problem on the unit cube, which is a model of a soft/hard material enclosure. As in [33], we only consider the two material subdomains case pictured in Figure 1, and we let Ω 2 be the union of the two internal cubes, while Ω 1 denotes the rest of the domain. The problem was discretized with linear finite elements on regular tetrahedral mesh, using zero Dirichlet boundary conditions on the boundary. The righthand side in all the tests, was chosen to correspond to the unit constant function, and the initial guess was a vector of zeros. Some of the computed numerical solutions are plotted in Figure 2.
Since we can always rescale the original equation, we can assume, without a loss of generality, that ω 2 = 1. In particular, when ω is a constant we will set it equal to one. This is the case that we set to explore first.  Table 1. Here denotes the refinement level corresponding to problem size N and mesh size h. We also use the "scientific" notation 1e+p to denote the number 10 p .
Several things are apparent from Table 1. First, when ρh 2 ω (the lower right corner in the tables) the problem is well conditioned and GS-CG is an efficient solver. Second, the convergence is largely independent of the jumps in ρ and the number of iterations is proportional to h −1 , as expected by Theorem 2.3. Finally, it is clear that the problem of hard enclosure, when ρ 2 > ρ 1 , is more difficult than the one of soft enclosure (ρ 1 > ρ 2 ).

5.2.
The case of constant ρ. Next, we consider the case when the mass term coefficient is a constant. As in the previous section, we first perform a parameter study to determine an appropriate scaling of ρ when ω 2 is fixed to be one. The results are presented in Table  6, and in many respects are similar to those from Table 1. For example, the number of GS-GC iterations doubles from one level to the next, though the actual numbers are several times larger than those in Table 1.
The results for GS-CG are shown in Table 7. Clearly, the problem of hard enclosure, when ω 1 is small, is much more challenging than the case of large ω 1 . In contrast to Table 2, the number of iterations increases significantly with the magnitude of the jump. This is due to the fact that the condition number is proportional to J (ω), see Theorem 2.3 and the discussion after Theorem 2.1 in [33].
Using Multigrid as a preconditioner resolves this problem, since there are only finite number of small eigenvalues corresponding to the jump in ω. The results in Table 10 demonstrate a nearly optimal convergence with respect to the mesh size.

5.3.
The case of discontinuous ω and ρ. In this section we present a numerical investigation of the general case when both ω and ρ are discontinuous. Note that the theory developed in this paper can be applied only if we can construct an interpolation operator which is stable in both the ρ-weighted and the ω-weighted L 2 -inner products. This is the case, for example if ω 1 ≤ ω 2 and ρ 1 ≤ ρ 2 . TABLE 11. Number of GS-CG iterations when ω 2 = 1, while ω 1 , ρ 1 and ρ 2 are allowed to vary. Each cell in the table represents a maximum over a range of values for ρ. Problem size N = 274, 625.
In Table 11 we show the results of a parameter study based on Gauss-Seidel preconditioning. We emphasize that each cell in this table represents a maximum over several possible values for ρ, which result in a jump of the same magnitude ρ 1 /ρ 2 . Clearly, the difficulty of the problem is determined mostly by the jump in ω, so we choose to concentrate on the most challenging case ω 1 = 10 −8 .