A preconditioned SSOR iteration method for solving complex symmetric system of linear equations

We present a preconditioned version of the symmetric successive overrelaxation (SSOR) iteration method for a class of complex symmetric linear systems. The convergence results of the proposed method are established and conditions under which the spectral radius of the iteration matrix of the method is smaller than that of the SSOR method are analyzed. Numerical experiments illustrate the theoretical results and depict the efficiency of the new iteration method.

1. Introduction. Many scientific and engineering applications often require the solution of large and sparse complex symmetric system of linear equations where W, T ∈ R n×n . We assume that W and T are symmetric positive semidefinite with at least one of them, e.g., W , being positive definite. Several effective iteration methods have been proposed in the literature for solving the system (1). For solving such symmetric systems, Bunse-Gerstner et al. presented a new conjugate gradient-type iterative method in [11]. Splitting iteration methods such as generalized successive overrelaxation (GSOR) [18], preconditioned GSOR (PGSOR) [16], scale splitting (SCSP) [15], two-parameter two-step SCSP (TTSCSP) in [19] and parameterized splitting (PS) in [20] proposed for the linear system (1). Edalatpour et al. in [13] solved inexactly the linear system (1) by two efficient methods. Using the Hermitian and skew-Hermitian splitting [5] of the coefficient matrix A, Bai et al. proposed the modified HSS (MHSS) iteration method in [2] and preconditioned MHSS (PMHSS) iteration method in [6].
In [1] it has been shown that by assuming A = W + iT , u = x + iy and b = p + iq and separating the real and imaginary parts of the system (1), it can be transformed into the real equivalent form This linear system can be formally regarded as a special case of the generalized saddle point problem [9,10]. Bai et al. in [7], applied the GSOR method for the block 2-by-2 linear system of different sub-blocks rather than the original complex linear system (1). Recently Salkuyeh et al. in [18] solved the real equivalent system (2) by the GSOR iterative method. Under our hypotheses, it can be easily proved that the matrix A is nonsingular. In [18] it has been shown that if W and T are symmetric positive definite and symmetric, respectively, then the GSOR method is convergent.
In [17], Liang et al. have applied the symmetric SOR (SSOR) method to the following preconditioned linear system and introduced the accelerated variant of the SSOR (ASSOR) method, in order to increase the convergence rate of the method. They proved that, under certain conditions on ω, the ASSOR method is convergent and obtained the optimal value of the iteration parameter ω.
For improving the proposed method, we use the preconditioner where 0 < α ∈ R and instead of solving the linear system (2), we solve the preconditioned linear system by the SSOR method and define the PSSOR method (for preconditioned version of SSOR). Evidently, when α = 1 the PSSOR method reduces to the ASSOR method. We study conditions under which the spectral radius of the PSSOR iteration matrix becomes smaller than that of the SSOR method. It is noteworthy that the motivation of choosing this preconditioner stems from [8] in which Benzi et al. have presented some examples of preconditioners for the real formulation (2) of the system (1). Throughout this paper, for a square matrix B, ρ(B) and σ(B) stand for the spectral radius and spectrum of B, respectively.
The rest of the paper is organized as follows. In Section 2, we present the PSSOR method and discuss its convergence properties. Section 3 is devoted to numerical experiments to verify the theoretical results given in Section 2. Finally, in Section 4 we present a brief concluding remarks to end the paper.
2. The PSSOR method and its convergence. In this section, we present the PSSOR iterative method and study its convergence properties. It is easy to see that the system (5) is equivalent to where W α = αW + T , T α = αT − W and p α = αp + q, q α = αq − p. Now, we use the SSOR method to solve the preconditioned linear system (6) and split the matrix According to this splitting, the SSOR method for solving the system (6) is given by Obviously, the above system can be rewritten as in which is the iteration matrix of the PSSOR method, with and Therefore the matrix M ω,α can be used as a preconditioner for the system (5). From (7), the PSSOR method for the solving the system (2) can be stated as follows.
Theorem 2.1. Let W and T ∈ R n×n be symmetric positive definite and symmetric positive semidefinite, respectively. Let also W α = αW + T and T α = αT − W , where α is a positive constant. Then, the PSSOR method is convergent if one of the following conditions holds: The optimal value of the iteration parameters ω and α which minimize the spectral radius of the PSSOR method are given by where µ min and µ max are the smallest and largest eigenvalues of G, respectively. The corresponding optimal convergence factor of the PSSOR method is . Therefore from (9), we can easily write it follows from (9) that Since the eigenvalues of the matrices H ω,α and H * ω,α are the same, we study the eigenvalues of matrix H * ω,α , instead of H ω,α . Using Remark 1, W α and T α are symmetric positive definite and symmetric, respectively. Then the matrix S α is symmetric and has the spectral decomposition S α = V ΓV −1 , where V ∈ R n×n is an invertible matrix and Γ = diag(µ 1 , µ 2 , . . . , µ n ) in which µ i 's are the eigenvalues of S α . By defining Hence the latter matrix is similar to H ω,α . Therefore, if λ is an eigenvalue of the above matrix, then we have for all µ ∈ σ(S α ). Letting β = ω(2 − ω), Eq. (13) can be written as Then, from Theorem 1 in [18], we get the following convergence condition After some simplifications, we conclude that the latter equation is equivalent to (i) if ρ(S α ) 1, then 0 < ω < 2 should be satisfied; ii) if ρ(S α ) > 1, then 0 < ω < 1 − ρ(Sα)−1 ρ(Sα)+1 or 1 + ρ(Sα)−1 ρ(Sα)+1 < ω < 2 should be satisfied. Also from Theorem 2 in [18], we see that for a given α the optimal value of β is obtained from β opt = 2 and therefore from On the other hand, using Theorem 2 in [18], we have This equation shows that to make ρ(H ωopt,α ) minimum we should minimize the value of ρ(S α ). According to Theorem 2 in [16], the optimal value of the parameter α which minimizes the spectral radius ρ(S α ) is given by where µ min and µ max are the smallest and largest eigenvalues of G, respectively. Now, using ρ(H ωopt,αopt ) = min α ρ(H ωopt,α ), Eq. (12) is proved.
The next theorem simplifies the convergence conditions of the PSSOR method.
Theorem 2.2. Let W, T ∈ R n×n be symmetric positive definite and symmetric positive semidefinite, respectively and α be a positive constant. Suppose that µ min and µ max be the smallest and largest eigenvalues of G = W − 1 2 T W − 1 2 , respectively, and S α = (αI + G) −1 (αG − I). The PSSOR method is convergent if one of the following conditions holds: (i) if µ max > 1, then a < α < d and 0 < ω < 2 should be satisfied; (ii) if µ max > 1, then α < a or α > d and 0 < ω < 1 − c or 1 + c < ω < 2 should be satisfied; (iii) if µ max 1, then α > a and 0 < ω < 2 should be satisfied; (iv) if µ max 1, then α < a and 0 < ω < 1 − c or 1 + c < ω < 2 should be satisfied; Proof. We have seen that the PSSOR method is convergent, if one of the following conditions holds true: (a) if ρ(S α ) 1, then 0 < ω < 2 should be satisfied; It is easy to see that the eigenvalues of S α = (αI + G) −1 (αG − I) are of the form where µ is an eigenvalue of the G. Since, the function h is increasing, the spectral radius of S α is equal to Hence ρ(S α ) 1 if and only if It is easy to see that the latter relation is equivalent to The first inequality of (18) holds if and only if Obviously, if µ max 1, then the second inequality of (18) always holds. When µ max > 1, the second inequality of (18) is equivalent to Therefore, from (a) we see that the PSSOR method is convergent in each of the cases (i) and (iii).
In the same manner, considering the case (b), it can be seen that the PSSOR method is convergent in the cases (ii) and (iv).

Numerical experiments.
In this section, we use three test problems from [2] to illustrate the feasibility and effectiveness of the PSSOR method to solve the system (2). We also compare the performance of this method with the PMHSS, GSOR, SSOR and ASSOR methods in terms of both the number of iterations (denoted by IT) and the total computing times denoted by CPU in the tables. The PMHSS iteration method are employed to solve the complex system (1) and other methods to solve the equivalent real system (2). In each iteration of the PMHSS, GSOR, SSOR, ASOR and PSSOR iteration methods, we use the Cholesky factorization of the coefficient matrices to solve the sub-systems.
All runs are implemented in Matlab R2014b with a Laptop with 2.40 GHz central processing unit (Intel(R) Core(TM) i7-5500), 8 GB memory and Windows 10 operating system.
We use a zero vector as an initial guess and the stopping criterion where u (k) = x (k) + iy (k) and f = (p; q).

Example 1. (See [2]) Consider the system of linear equations
where τ is the time step-size and K is the five-point centered difference matrix approximating the negative Laplacian operator L = −∆ with homogeneous Dirichlet boundary conditions, on a uniform mesh in the unit square [0, 1] × [0, 1] with the mesh-size h = 1 m+1 . The matrix K ∈ R n×n possesses the tensor-product form K = I ⊗ V m + V m ⊗ I, with V m = h −2 tridiag(−1, 2, −1) ∈ R m×m . Hence, K is an n × n block-tridiagonal matrix, with n = m 2 . We take I, and T = K + 3 + √ 3 τ I and the right-hand side vector b with its jth entry b j being given by b j = (1 − i)j τ (j + 1) 2 , j = 1, 2, . . . , n. In our tests, we take τ = h. Furthermore, we normalize coefficient matrix and right-hand side by multiplying both by h 2 .

Example 2. (See [2]) Consider the system of linear equations (1) as following
where M and K are the inertia and the stiffness matrices, C V and C H are the viscous and the hysteretic damping matrices, respectively, and ω is the driving circular frequency. We take C H = µK with µ a damping coefficient, M = I, C V = 10I, and K the five-point centered difference matrix approximating the negative Laplacian operator with homogeneous Dirichlet boundary conditions, on a uniform mesh in the unit square [0, 1] × [0, 1] with the mesh-size h = 1/(m + 1). The matrix K ∈ R n×n possesses the tensor-product form K = I ⊗ V m + V m ⊗ I, with V m = h −2 tridiag(−1, 2, −1) ∈ R m×m . Hence, K is an n × n block-tridiagonal matrix, with n = m 2 . In addition, we set µ = π, ω = π, and the right-hand side vector b to be b = (1 + i)A1, with 1 being the vector of all entries equal to 1. As before, we normalize the system by multiplying both sides through by h 2 .
Example 3. (See [2]) Consider the linear system of equations (1) as following m − e m e T 1 ∈ R m×m and e 1 and e m are the first and last unit vectors in R m , respectively. We take the right-hand side vector b to be b = (1 + i)A1, with 1 being the vector of all entries equal to 1.
Here T and W correspond to the five-point centered difference matrices approximating the negative Laplacian operator with homogeneous Dirichlet boundary conditions and periodic boundary conditions, respectively, on a uniform mesh in the unit square [0, 1] × [0, 1] with the mesh-size h = 1/(m + 1). In Table 1, the optimal parameters of the PMHSS, GSOR, SSOR, ASSOR and PSSOR methods are presented. The optimal parameters of the GSOR method are chosen based on Theorem 2 of [18] and the ones of the other methods are experimentally found.
In Tables 2, 3 and 4, we report the numerical results for Examples 1, 2 and 3. We present the iteration numbers (IT) and CPU time (CPU) for the PMHSS, GSOR, SSOR, ASSOR and PSSOR iteration methods for different grids. As we observe, for all the test examples the iteration counts of the PSSOR method is always less that that of the other methods. From the CPU time of view, this is the case for large problems. In fact, we see the the PSSOR method outperforms the other methods for large systems. 4. Conclusion. In this paper we have presented a preconditioned version of the symmetric successive overrelaxation (PSSOR) iterative method to solve the equivalent real system of linear system (2), where W is symmetric positive definite and T   is symmetric positive semidefinite. Convergence properties of this method have also been investigated. We have compared the numerical results of the PSSOR method with PMHSS, GSOR, SSOR and ASSOR. Numerical results show that PSSOR method is superior to the other methods in terms of both iteration number and CPU time for large problems.