# A Modified SSOR Preconditioning Strategy for Helmholtz Equations

Authors: Shi-Liang Wu; Cui-Xia Li
Journal: Journal of Applied Mathematics
 (2012)

---

## Abstract

The finite difference method discretization of Helmholtz equations usually leads to
the large spare linear systems. Since the coefficient matrix is frequently indefinite, it
is difficult to solve iteratively. In this paper, a modified symmetric successive overrelaxation (MSSOR) preconditioning strategy is constructed based on the coefficient
matrix and employed to speed up the convergence rate of iterative methods. The idea
is to increase the values of diagonal elements of the coefficient matrix to obtain better
preconditioners for the original linear systems. Compared with SSOR preconditioner,
MSSOR preconditioner has no additional computational cost to improve the convergence rate of iterative methods. Numerical results demonstrate that this method can
reduce both the number of iterations and the computational time significantly with
low cost for construction and implementation of preconditioners.


---

## Body

## 1. Introduction

The finite difference method is one of the most effective and popular techniques in computational electromagnetics and seismology, such as time-harmonic wave propagations, scattering phenomena arising in acoustic and optical problems, and electromagnetics scattering from a large cavity. More information about applications of this method in electromagnetics can be found in [1–5].In this paper, we focus on the following form of the complex Helmholtz equation:(1.1)-Δu-pu+iqu=finΩ,u=gon⁡∂Ω.
Here Ω is a bounded region in ℝ2. p and q are real continuous coefficient functions on Ω¯, while f and g are given continuous functions on Ω¯ and ∂Ω, respectively.To conveniently find numerical solutions of (1.1), the Laplace operator is approximated by using the second-order accurate 5-point difference stencil:(1.2)Lh=̇1h2[11-411]h.
Making use of the above stencil, the following linear system is obtained:
(1.3)Ax=(A-h2D+ih2E)x=b,
where D(E) is a diagonal matrix whose diagonal elements are just the values of p(q) at the mesh points and A is the symmetric positive definite M-matrix arising from the discrete Laplace operator and is of the block tridiagonal form(1.4)A=[G1F2E2G2⋱⋱⋱FmEmGm],with(1.5)Gk=[4-1-14⋱⋱⋱-1-14]nk×nk,k=1,2,…,m,Ek=FkT={(Onk×pk-InkOnk×qk),pk,qk≥0,pk+qk=nk-1-nkifnk≤nk-1,(Onk-1×sk-Ink-1Onk-1×tk)T,sk,tk≥0,sk+tk=nk-nk-1ifnk>nk-1,k=2,3,…,m.Obviously, from the linear systems (1.3), it is not difficult to find that the matrix 𝒜 is a complex symmetric coefficient matrix. Matrix 𝒜 becomes highly indefinite and ill-conditioned as p is a sufficiently large positive number. So, large amount of computation times and memory are needed in order to solve the linear systems (1.3) efficiently.As is well known, direct methods and iterative methods can be employed to solve the linear systems (1.3). The former is widely employed when the order of the coefficient matrix 𝒜 is not too large and is usually regarded as robust methods. The memory and the computational requirements for solving the large sparse linear systems may seriously challenge the most efficient direct solution method available today. Currently, the latter employed to solve the large sparse linear systems is popular. The reason is that iterative methods are easier to implement efficiently on high performance computers than direct methods. In practice, a natural choice is that we make use of iterative methods instead of direct methods to solve the large sparse linear systems.At present, Krylov subspace methods are considered as one kind of the important and efficient iterative techniques for solving the large sparse linear systems because these methods are cheap to be implemented and are able to fully exploit the sparsity of the coefficient matrix. It is well known that the convergence speed of Krylov subspace methods depends on the distribution of the eigenvalues of the coefficient matrix [6]. When the coefficient matrix is typically extremely ill-conditioned and highly indefinite, the convergence of Krylov subspace methods can be unacceptably slow. In this case, Krylov subspace methods are not competitive without a good preconditioner. That is, preconditioning technique is a key ingredient for the success of Krylov subspace methods in applications. The idea of preconditioning technique is based on consideration of the linear system with the same solution as the original equation. The problem is that each preconditioning technique is suited for a different type of problem. Until current days, no robust preconditioning technique appears for all or at least much types of problems. Finding a good preconditioner to solve a given large sparse linear systems is often viewed as a combination of art and science.In recent years, a great deal of effort has been invested in solving indefinite linear systems from the discrete Helmholtz equations. Most of the work has been aimed at developing effective preconditioning techniques. In general, there exist two classes of preconditioners for Helmholtz equations: the “operator-based" preconditioning technique and the “matrix-based" preconditioning technique.The former is built based on an operator, such as the Laplace preconditioner [2, 3, 7–9], Analytic ILU [10], the Separation-of-variables [11]. The purpose of this class of preconditioners is that the spectrum of the corresponding preconditioned matrix is favorably clustered. Its advantage is that this operator does not have to be a representation of the inverse of the Helmholtz operator.The latter is established based on an approximation of the inverse of the coefficient matrix. For this class, one of the natural and simplest ways of structuring a preconditioner is to employ a diagonal or block diagonal of the coefficient matrix as a preconditioner [12]. The above two diagonal preconditioners have no remarkable reduce with respect to the iterative number and CPU time. Another one of the simplest preconditioners is to perform an incomplete factorization (ILU) of the coefficient matrix [1]. The main idea of ILU factorizations depends on the implementation of Gaussian elimination which is used, see the survey [13] and the related references therein.When the coefficient matrix of the linear systems (1.1) is complex symmetric and indefinite, it is difficult to solve iteratively. Using the symmetric successive overrelaxation (SSOR) as a preconditioner preserves the symmetry of the iterative matrix and also taking little initialization cost, which in some cases makes it preferable over other factorization methods such as ILU. So far, some variant SSOR preconditioning techniques have been proposed to improve the convergence rate of the corresponding iterative method for solving the linear systems. Mazzia and Alan [14] introduced a shift parameter to develop a shifted SSOR preconditioner for solving the linear systems from an electromagnetics application. Bai in [15] used a (block) diagonal matrix instead of the diagonal matrix of the coefficient matrix to establish a modified (block) SSOR preconditioner for the second-order selfadjoint elliptic boundary value problems. Chen et al. [16] used a diagonal matrix with a relaxation parameter instead of the diagonal matrix of the coefficient matrix to establish a modified block SSOR preconditioner for the Biot’s consolidation equations. Ran and Yuan in [17] also discussed a class of modified (block) SSOR preconditioners for linear systems from steady incompressible viscous flow problems. We refer the reader to [18, 19] for a general discussion.Although the SSOR preconditioner with Krylov subspace methods can improve convergence improvement, the disadvantage of the SSOR preconditioner is that the convergence rate may still remain unsatisfactorily slow in some cases, especially in indefinite linear systems. In this paper, a modified symmetric successive overrelaxation (MSSOR) preconditioning strategy is presented, which can significantly improve the convergence speed and CPU time. Our motivation for this method arises from the solution of complex symmetric and indefinite linear systems from the discrete Helmholtz equations. The idea is to increase the values of diagonal elements of the coefficient matrix to obtain better preconditioners for the original linear systems, which is different from [14–17]. This modification does not require any significant computational cost as compared with the original SSOR preconditioner and also requires no additional storage cost.The remainder of this paper is organized as follows. In Section2, the MSSOR preconditioner for solving the resulting linear system is presented. In Section 3, numerical experiments are given to illustrate the efficiency of the presented preconditioner. Finally, in Section 4 some conclusions are drawn.


## 2. Modified SSOR Preconditioner

To improve the convergence rate of iterative methods, an appropriate preconditioner should be incorporated. That is, it is often preferable to solve the preconditioned linear system as follows:(2.1)P-1Ax=P-1b,
where P, called the preconditioner, is a nonsingular matrix. The choice of the preconditioner P plays an important role in actual implements. In general, the preconditioner P is chosen such that the condition number of the preconditioned matrix P-1𝒜 is less than that of the original matrix 𝒜. Based on the excellent survey of [13] by Benzi, a good preconditioner should meet the following requirements:(1)
the preconditioned system should be easy to solve;(2)
the preconditioner should be cheap to construct and apply.Certainly, the best choice forP-1 is the inverse of 𝒜. However, it is unpractical in actual implements because the cost of the computation of 𝒜-1 may be high. If 𝒜 is a symmetric positive definite matrix, the approximation of 𝒜-1 can be replaced by SSOR or multigrid. However, in fact, the Helmholtz equation often results in an indefinite linear system, for which SSOR or multi-grid may be not guaranteed to converge.To introduce the modified SSOR preconditioner, a brief review of the classical and well-known SSOR preconditioner is needed. The SSOR preconditioner is established by the SSOR iterative method, which is a symmetric version of the well-known SOR iterative method. Based on matrix splitting, the coefficient matrix𝒜 is split as follows:
(2.2)A=D+L+LT,
where 𝒟 and ℒ are the diagonal parts and strictly lower triangular of 𝒜. According to the foregoing matrix splitting (2.2), the standard SSOR preconditioner [6, 20] is defined by
(2.3)PSSOR=(D+L)D-1(D+LT).
It is difficult to show theoretically the behavior of a preconditioner when the coefficient matrix 𝒜 is a large, sparse, and symmetric indefinite. The SSOR iterative method is not convergent, but PSSOR may be still used as a preconditioner. By a simple modification on the original indefinite linear systems (1.3), we establish the following coefficient matrix:
(2.4)A̅=A+h2D+h2Ei=D̃+L+LT,
where
(2.5)D̃=diag⁡(A)+h2D+h2Ei.
Obviously, 𝒜̅ is a symmetric and positive stable H-matrix. To increase the values of diagonal elements of the coefficient matrix to obtain better preconditioners for the original linear systems and reduce computation times and amount of memory, based on (2.4), the MSSOR preconditioner is defined by
(2.6)PMSSOR=(D̅+L)D̅-1(D̅+LT),
with(2.7)D̅=diag⁡(|D̃|),|D̃|=|d̃ii|(i=1,2,…,n).
This idea is based on an absolute diagonal scaling technique, which is cheap and easy to implement.Since the coefficient matrix of the linear systems (1.3) is neither positive definite nor Hermitian with p being a sufficiently large positive number, the Conjugate Gradient (CG) method [21] may breakdown. To solve the complex symmetric linear systems, van der Vorst and Melissen [22] proposed the conjugate orthogonal conjugate gradient (COCG) method, which is regarded as an extension of CG method.To solve the linear systems (1.3) efficiently, (1.3) is transformed into the following form with the preconditioner PMSSOR, that is,
(2.8)PMSSOR(A-h2D+ih2E)x=PMSSORb.
Then the MSSOR preconditioned COCG (PCOCG) method can be employed to solve the preconditioned linear systems (2.8).In the following, we give the MSSOR preconditioned COCG method for solving the linear systems (2.8). The MSSOR preconditioned COCG (PCOCG) algorithm is described as followsAlgorithm PCOCG [22]: given an initial guess x0(1)
v0=b-𝒜x0;(2)
p-1=0;β-1=0;(3)
w0=PMSSOR-1v0;ρ0=v0Tw0;(4)
fori=0,1,2,…(5)
pi=wi+βi-1pi-1;(6)
ui=𝒜pi;(7)
μi=ujTpj; if μi=0 then quit (failure);(8)
α=ρi/μi;(9)
xi+1=xi+αpi;vi+1=vi-αui(10)
ifxi+1 is accurate enough, then quit (convergence);(11)
wi+1=PMSSOR-1vi+1;(12)
ρi+1=vi+1Twi+1; if |ρi+1| too small, then quit (failure);(13)
βi=ρi+1/ρi;(14)
End fori.It is not difficult to find that the main computation of algorithm PCOCG involves one matrix-vector multiplication and two triangular linear systems. These computations are very easy to implement. The main advantage is no extra computational cost in construction of MSSOR preconditioner.Note that thetranspose in all dot products in this algorithm is essential [23]. Meanwhile, note that two different breakdowns of this algorithm may occur: one is that if vi+1TPMSSOR-1vi+1 is too small, but vi+1 exists (line 12), algorithm PCOCG breaks down and the other is that when the search direction pi≠0, but piT𝒜pi=0 (line 7), algorithm PCOCG breaks down. The breakdown can be fixed to some extent by restarting the process [22], such as the restarted process in GMRES [24]. However, breakdown scarcely happens in the actual computation of the Helmholtz equation.


## 3. Numerical Experiments

In this section, some numerical experiments are given to demonstrate the performance of both preconditionerPSSOR and preconditioner PMSSOR for solving the Helmholtz equation.Example 3.1 (see [25]).
Consider the following complex Helmholtz equation:(3.1)-Δu-pu+iqu=finΩ,u=gon⁡∂Ω,
where Ω=[0,1]×[0,1] and p≥0 and q≥0 are real constants. Discretizing (3.1) with the approach above in introduction, we obtain the complex symmetric indefinite linear systems 𝒜x=(A-h2D+ih2E)x=b, and f and g are adjusted such that b=𝒜e(e=(1,1,…,1)T).All tests are started from the zero vector, preformed in MATLAB with machine precision10-16. The COCG iteration terminates if the relative residual error satisfies ∥r(k)∥2/∥r(0)∥2<10-6 or the iteration number is more than 500.In Tables1, 2, 3, and 4, we present some iteration results to illustrate the convergence behaviors of the COCG method preconditioned by PSSOR and PMSSOR to solve the complex symmetric indefinite linear systems with the different values of p and q. In Tables 1–4, (·,·) denotes the values of p and q. “CPU(s)’’ denotes the time (in seconds) required to solve a problem. “IT’’ denotes the number of iteration.Table 1
Iteration numbers for the two-dimensional Helmholtz equations forh=1/19 meshes, using COCG with the preconditioner PSSOR and PMSSOR for solving the complex symmetric indefinite linear systems.
(p, q)hPreconditioner(800, 10)(800, 20)(800, 30)(800, 40)1/19PSSORIT246242239214CPU(s)0.84380.81250.78130.7344PMSSORIT138133126117CPU(s)0.250.21880.20310.1875Table 2
Iteration numbers for the two-dimensional Helmholtz equations forh=1/34 meshes, using COCG with the preconditioner PSSOR and PMSSOR for solving the complex symmetric indefinite linear systems.
(p, q)hPreconditioner(1400, 40)(1500, 40)(1600, 40)(1700, 40)1/34PSSORIT230274310336CPU(s)3.28133.93754.84384.7188PMSSORIT214228237249CPU(s)1.43751.51561.65631.6250Table 3
Iteration numbers for the two-dimensional Helmholtz equations forh=1/64 and p=4100, using COCG with the preconditioner PSSOR and PMSSOR.
qhPreconditioner1001201501601/64PSSORIT471456374349CPU(s)32.934430.968827.984426.0156PMSSORIT373366347326CPU(s)12.468812.250011.609410.9063Table 4
Iteration numbers for the two-dimensional Helmholtz equations forh=1/119 and q=2000, using COCG with the preconditioner PSSOR and PMSSOR.
phPreconditioner150001550016000165001/119PSSORIT148153147165CPU(s)55.359457.453157.859452.4844PMSSORIT134137139144CPU(s)22.62522.906323.281323.7344From Tables1–4, it is not difficult to find that when the COCG method preconditioned by PSSOR and PMSSOR is used to solve the complex symmetric indefinite linear systems, the convergence rate of the preconditioner PMSSOR is more efficient than that of the preconditioner PSSOR by the iteration numbers and CPU time. That is, the preconditioner PMSSOR outperforms the preconditioner PSSOR under certain conditions. Compared with the preconditioner PSSOR, the preconditioner PMSSOR may be the “preferential” choice under certain conditions.


## 4. Conclusions

In this paper, MSSOR preconditioned COCG algorithm has been applied for solving the complex symmetric indefinite systems arising from Helmholtz equations. Due to the reduction of the iteration numbers and CPU time, the MSSOR preconditioner presented is feasible and effective. Without extra costs, MSSOR preconditioner is more efficient than SSOR preconditioner.



---

*Source: 365124-2011-12-28.xml*