Generalized Preconditioned MHSS Method for a Class of Complex Symmetric Linear Systems

Abstract

Based on the modified Hermitian and skew-Hermitian splitting (MHSS) and preconditioned MHSS (PMHSS) methods, a generalized preconditioned MHSS (GPMHSS) method for a class of complex symmetric linear systems is presented. Theoretical analysis gives an upper bound for the spectral radius of the iteration matrix. From a practical point of view, we have analyzed and implemented inexact GPMHSS (IGPMHSS) iteration, which employs Krylov subspace methods as its inner processes. Numerical experiments are reported to confirm the efficiency of the proposed methods.

1. Introduction

Consider an iterative solution of the system of linear equations as follows: 𝐴𝑥=𝑏,𝐴∈ℂ^𝑛×𝑛,𝑥,𝑏∈ℂ^𝑛, (1) where 𝐴∈ℂ^𝑛×𝑛 is a complex symmetric matrix of the following form: 𝐴=𝑊+𝑖𝑇, (2) with 𝑊∈ℝ^𝑛×𝑛 being symmetric positive definite and 𝑇∈ℝ^𝑛×𝑛 being symmetric positive semidefinite. Here and in the sequel, we use 𝑖=√−1 as the imaginary unit. One can readily verify that 𝐴 is non-Hermitian, that is to say, the linear system (1) is a non-Hermitian linear system. System such as (1) is important and arises in a variety of scientific and engineering applications, including structural dynamics [1–4], diffuse optical tomography [5, 6], FFT-based solution of certain time-dependent PDEs [7], lattice quantum chromodynamics [8], molecular dynamics and fluid dynamics [9], quantum chemistry, and eddy current problem [10, 11]. One can see [12, 13] for more examples and additional references. In order to solve (1) more effectively, many efficient numerical algorithms have been proposed in the literature, see [14–20].

Based on the specific structure of the coefficient matrix 𝐴, one can verify that the Hermitian and skew-Hermitian parts of the coefficient matrix 𝐴, respectively, are 1𝐻=2𝐴+𝐴^∗1=𝑊,𝑆=2𝐴−𝐴^∗=𝑖𝑇. (3) Obviously, the above Hermitian and skew-Hermitian splitting of the coefficient matrix 𝐴 is in line with the real and imaginary parts splitting of the coefficient matrix 𝐴. Based on the HSS method [21], Bai et al. [2] skillfully designed a modified HSS (MHSS) method to solve the complex symmetric linear system (1) which is described below.

The MHSS Method. Let 𝑥⁽⁰⁾∈ℂ^𝑛 be an arbitrary initial guess. For 𝑘=0,1,2,…, until the sequence of iterates {𝑥^(𝑘)}^∞_𝑘=0 converges, compute the next iterate 𝑥^(𝑘+1) according to the following procedure: (𝛼𝐼+𝑊)𝑥^(𝑘+1/2)=(𝛼𝐼−𝑖𝑇)𝑥^(𝑘)()𝑥+𝑏,𝛼𝐼+𝑇^(𝑘+1)=(𝛼𝐼+𝑖𝑊)𝑥^(𝑘+1/2)−𝑖𝑏, (4) where 𝛼 is a given positive constant and 𝐼 is an identity matrix.

The potential advantage of the MHSS method over the HSS method [21] for solving the complex symmetric linear system (1) is that only two linear subsystems with coefficient matrices 𝛼𝐼+𝑊 and 𝛼𝐼+𝑇, both being real and symmetric positive definite, need to be solved at each step. Therefore, in this case, these two linear subsystems can be solved efficiently using mostly real arithmetic either exactly by a sparse Cholesky factorization or inexactly by conjugated gradient scheme. That is to say, the MHSS method successfully avoids solving a shifted skew-Hermitian linear subsystem with coefficient matrix 𝛼𝐼+𝑖𝑇.

Theoretical analysis in [2] shows that the MHSS method converges unconditionally to the unique solution of the complex symmetric linear system (1) when 𝑊∈ℝ^𝑛×𝑛 is real symmetric positive definite and 𝑇∈ℝ^𝑛×𝑛 is real symmetric positive semidefinite. The corresponding optimum parameter 𝛼=𝜆_min(𝑊)𝜆_max(𝑊) is obtained to minimize an upper bound on the spectral radius of the iteration matrix associated with (4).

The MHSS method immediately attracted considerable attention and resulted in many papers devoted to various aspects of the new algorithms. For instance, preconditioned modified Hermitian and skew-Hermitian splitting (PMHSS) iteration in [3], lopsided preconditioned modified Hermitian and skew-Hermitian splitting (LPMHSS) iteration in [22], amongst others. On the other hand, the MHSS method was successfully extended to the solution of control problems in [23].

In this paper, based on the splitting (3), we generalize the MHSS iterative scheme into a new approach, called generalized preconditioned MHSS (GPMHSS) iteration. By introducing two symmetric positive definite matrices 𝑃₁∈ℝ^𝑛×𝑛 and 𝑃₂∈ℝ^𝑛×𝑛, the GPMHSS iterative scheme works as follows.

The GPMHSS Method. Let 𝑥⁽⁰⁾∈ℂ^𝑛 be an arbitrary initial guess. For 𝑘=0,1,2,…, until the sequence of iterates {𝑥^(𝑘)}^∞_𝑘=0 converges, compute the next iterate 𝑥^(𝑘+1) according to the following procedure: 𝛼𝑃₁+𝑊𝑥^(𝑘+1/2)=𝛼𝑃₁−𝑖𝑇𝑥^(𝑘)+𝑏,𝛽𝑃₂+𝑇𝑥^(𝑘+1)=𝛽𝑃₂+𝑖𝑊𝑥^(𝑘+1/2)−𝑖𝑏, (5) where 𝛼 is a nonnegative constant and 𝛽 is a positive constant.

Note that the GPMHSS iteration method can cover many existing variants of the standard MHSS iteration. For instance, when 𝛼=𝛽 and 𝑃₁=𝑃₂=𝐼, the GPMHSS iteration method is equivalent to the standard MHSS iteration method in [2]; when 𝛼=𝛽 and 𝑃₁=𝑃₂, the GPMHSS iteration method is equivalent to the standard PMHSS iteration method in [3]; when 𝛼=0 and 𝑃₂=𝐼, it leads to the LPMHSS iteration method in [22].

Theoretical analysis gives an upper bound about the contraction factor of GPMHSS iteration method, which shows the relations among GPMHSS, MHSS, and other existing variants. From a practical point of view, we also discuss the inexact variants of the GPMHSS iteration method and their implementation. A number of numerical experiments are presented to illustrate the advantages of the GPMHSS and IGPMHSS methods.

This paper is organized as follows. In Section 2, we study the convergence properties of the GPMHSS iteration method. In Section 3, we discuss the implementation of GPMHSS iteration method and the corresponding inexact GPMHSS (IPGMHSS) iteration method. Numerical experiments are reported to confirm the efficiency of the proposed methods in Section 4. Finally, we end the paper with concluding remarks in Section 5.

2. Convergence Analysis for the GPMHSS Method

In this section, the convergence of the GPMHSS method is studied and an upper bound for the contraction factor of the GPMHSS method is derived. The GPMHSS iteration method can be generalized into a two-step splitting iteration framework. The following lemma is required to study the convergence properties of the GPMHSS method.

The spectral radius of the matrix 𝐴 is the nonnegative real number 𝜌(𝐴)=max{|𝜆|∶𝜆∈𝜎(𝐴)}, where 𝜎(𝐴) denotes the spectrum of matrix 𝐴. In fact, we have the basic property on spectral radius of the product of two matrices that 𝜌(𝐴𝐵)=𝜌(𝐵𝐴), which is used in the proof of the following theorem.

Lemma 1 (see [21]). Let 𝐴∈ℂ^𝑛×𝑛, 𝐴=𝑀_𝑖−𝑁_𝑖(𝑖=1,2) be two splittings of 𝐴, and 𝑥⁽⁰⁾∈ℂ^𝑛 be a given initial vector. If {𝑥^(𝑘)} is a two-step iteration sequence defined by 𝑀₁𝑥^(𝑘+1/2)=𝑁₁𝑥^(𝑘)𝑀+𝑏,₂𝑥^(𝑘+1)=𝑁₂𝑥^(𝑘+1/2)+𝑏, (6)𝑘=0,1,…, then 𝑥^(𝑘+1)=𝑀₂⁻¹𝑁₂𝑀₁⁻¹𝑁₁𝑥^(𝑘)+𝑀₂⁻¹𝐼+𝑁₂𝑀₁⁻¹𝑏,𝑘=0,1,…. (7) Moreover, if the spectral radius 𝜌(𝑀₂⁻¹𝑁₂𝑀₁⁻¹𝑁₁)<1, then the iterative sequence {𝑥^(𝑘)} converges to the unique solution 𝑥_∗∈ℂ^𝑛 of the system (1) for all initial vectors 𝑥⁽⁰⁾∈ℂ^𝑛.

Applying this lemma to the GPMHSS method, we get the convergence property in the following theorem.

Theorem 2. Let 𝑃₁∈ℝ^𝑛×𝑛 and 𝑃₂∈ℝ^𝑛×𝑛 be two symmetric positive definite matrices. Let 𝐴=𝑊+𝑖𝑇∈ℂ^𝑛×n, with 𝑊∈ℝ^𝑛×𝑛 and 𝑇∈ℝ^𝑛×𝑛 symmetric positive definite and symmetric positive semidefinite, respectively, and let 𝛼 be a nonnegative constant and let 𝛽 be a positive constant. Then the iteration matrix 𝑀_𝛼,𝛽 of GPMHSS method is 𝑀_𝛼,𝛽=𝛽𝑃₂+𝑇⁻¹𝛽𝑃₂+𝑖𝑊𝛼𝑃₁+𝑊⁻¹𝛼𝑃₁−𝑖𝑇. (8) Denote 𝛾(𝑃𝑃^𝑇),𝜆(𝑊), and 𝜇(𝑇) to be the spectral sets of the matrices 𝑃𝑃^𝑇, 𝑊, and 𝑇, respectively, where 𝑃=𝑃₁^−1/2𝑃₂^1/2, 𝑊=𝑃₁^−1/2𝑊𝑃₁^−1/2, and 𝑇=𝑃₂^−1/2𝑇𝑃₂^−1/2. Then, 𝜌𝑀_𝛼,𝛽≤𝛿=𝜅𝑃𝛿₁𝛿₂, (9) where 𝜅(𝑃) denotes the spectral condition number of the matrix 𝑃, and 𝛿₁=𝛽max_{𝛾𝑖∈𝛾}𝑃𝑃^𝑇𝛾_𝑖−1+𝛽−𝛼max_{𝜆𝑖∈𝜆}𝑊1𝛼+𝜆_𝑖+max_{𝜆𝑖∈𝜆(}𝑊₎𝛼²+𝜆²_𝑖𝛼+𝜆_𝑖,𝛿₂=𝛼max_{𝛾𝑖∈𝜆}𝑃𝑃^𝑇1𝛾_𝑖−1+𝛽−𝛼max_{𝜇𝑖∈𝜇(}𝑇₎1𝛽+𝜇_𝑖+1. (10)

Proof. Let 𝑀₁=𝛼𝑃₁+𝑊, 𝑁₁=𝛼𝑃₁−𝑖𝑇, 𝑀₂=𝛽𝑃₂+𝑇, and 𝑁₂=𝛽𝑃₂+𝑖𝑊. Obviously, 𝛼𝑃₁+𝑊 and 𝛽𝑃₂+𝑇 are nonsingular for any nonnegative constants 𝛼 and positive constants 𝛽. So formula (8) is valid.
Let 𝑊=𝑃₁^−1/2𝑊𝑃₁^−1/2, 𝑇=𝑃₂^−1/2𝑇𝑃₂^−1/2, and 𝑃=𝑃₁^−1/2𝑃₂^1/2. Then, 𝛼𝑃₁+𝑊=𝑃₁^1/2𝛼𝐼+𝑊𝑃₁^1/2,𝛽𝑃₂+𝑖𝑊=𝑃₁^1/2𝛽𝑃𝑃^𝑇+𝑖𝑊𝑃₁^1/2,𝛽𝑃₂+𝑇=𝑃₂^1/2𝛽𝐼+𝑇𝑃₂^1/2,𝛼𝑃₁−𝑖𝑇=𝑃₂^1/2𝛼𝑃⁻¹𝑃^−𝑇−𝑖𝑇𝑃₂^1/2. (11) Hence, 𝛽𝑃₂+𝑖𝑊𝛼𝑃₁+𝑊⁻¹=𝑃₁^1/2𝛽𝑃𝑃^𝑇+𝑖𝑊𝛼𝐼+𝑊⁻¹𝑃₁^−1/2=𝑃₁^1/2𝛽𝑃𝑃^𝑇−𝛼𝐼𝛼𝐼+𝑊⁻¹+𝛼𝐼+𝑖𝑊𝛼𝐼+𝑊⁻¹𝑃₁^−1/2=𝑃₁^1/2𝛽𝑃𝑃^𝑇−𝐼𝛼𝐼+𝑊⁻¹+𝛽−𝛼𝛼𝐼+𝑊⁻¹+𝛼𝐼+𝑖𝑊𝛼𝐼+𝑊⁻¹𝑃₁^−1/2,𝛼𝑃₁−𝑖𝑇𝛽𝑃₂+𝑇⁻¹=𝑃₂^1/2𝛼𝑃⁻¹𝑃^−𝑇−𝑖𝑇𝛽𝐼+𝑇⁻¹𝑃₂^−1/2=𝑃₂^1/2𝛼𝑃⁻¹𝑃^−𝑇−𝛽𝐼𝛽𝐼+𝑇⁻¹+𝛽𝐼−𝑖𝑇𝛽𝐼+𝑇⁻¹𝑃₂^−1/2=𝑃₂^1/2𝛼𝑃⁻¹𝑃^−𝑇−𝐼𝛽𝐼+𝑇⁻¹+×𝑇𝛼−𝛽𝛽𝐼+⁻¹+𝛽𝐼−𝑖𝑇𝛽𝐼+𝑇⁻¹𝑃₂^−1/2. (12) Further, we have 𝜌𝑀_𝛼,𝛽=𝜌𝛽𝑃₂+𝑇⁻¹𝛽𝑃₂+𝑖𝑊𝛼𝑃₁+𝑊⁻¹𝛼𝑃₁−𝑖𝑇=𝜌𝛽𝑃₂+𝑖𝑊𝛼𝑃₁+𝑊⁻¹𝛼𝑃₁−𝑖𝑇𝛽𝑃₂+𝑇⁻¹=𝜌𝑃₁^1/2𝛽𝑃𝑃^𝑇+𝑖𝑊𝛼𝐼+𝑊⁻¹×𝑃𝛼𝑃⁻¹𝑃^−𝑇−𝑖𝑇𝛽𝐼+𝑇⁻¹𝑃₂^−1/2=𝜌𝛽𝑃𝑃^𝑇+𝑖𝑊𝛼𝐼+𝑊⁻¹𝑃𝛼𝑃⁻¹𝑃^−𝑇𝑇×𝑇−𝑖𝛽𝐼+⁻¹𝑃⁻¹≤𝛽𝑃𝑃^𝑇+𝑖𝑊𝛼𝐼+𝑊⁻¹𝑃𝛼𝑃⁻¹𝑃^−𝑇𝑇×𝑇−𝑖𝛽𝐼+⁻¹𝑃⁻¹₂≤𝜅𝑃𝛽𝑃𝑃^𝑇−𝐼₂𝛼𝐼+𝑊⁻¹₂+×𝑊𝛽−𝛼𝛼𝐼+⁻¹₂+𝛼𝐼+𝑖𝑊𝛼𝐼+𝑊⁻¹₂⋅𝛼𝑃⁻¹𝑃^−𝑇−𝐼₂𝛽𝐼+𝑇⁻¹₂+×𝑇𝛼−𝛽𝛽𝐼+⁻¹₂+𝛽𝐼−𝑖𝑇𝛽𝐼+𝑇⁻¹₂. (13) Through simple calculations, we can get that 𝛼𝐼+𝑖𝑊𝛼𝐼+𝑊⁻¹₂≤max_{𝜆𝑖∈𝜆(}𝑊₎𝛼²+𝜆²_𝑖𝛼+𝜆_𝑖,𝛽𝐼−𝑖𝑇𝛽𝐼+𝑇⁻¹₂≤max_{𝜇𝑖∈𝜇(}𝑇₎𝛽²+𝜇²_𝑖𝛽+𝜇_𝑖𝑃𝑃≤1,_𝑇−𝐼₂≤max_{𝛾𝑖∈𝛾(}𝑃𝑃^𝑇₎𝛾_𝑖,𝑃−1₋₁𝑃^−𝑇−𝐼₂≤max_{𝛾𝑖∈𝛾(}𝑃𝑃^𝑇₎1𝛾_𝑖−1. (14) The above inequalities give the upper bound for 𝜌(𝑀_𝛼,𝛽) in (9).

Some remarks on Theorem 2 are given below.

(i) When 𝑃₁=𝑃₂, GPMHSS reduces to a two-parameter PMHSS method. In this case, 𝑃=𝑃₁^−1/2𝑃₂^1/2=𝐼 and 𝜅(𝑃)=1. In the meantime, the upper bound in (9) results in 𝛿=max_{𝜆𝑖∈𝜆}𝑊𝛽−𝛼𝛼+𝜆_𝑖+max_{𝜆𝑖∈𝜆}𝑊𝛼²+𝜆²_𝑖𝛼+𝜆_𝑖×max_{𝜇𝑖∈𝜇}𝑇𝛽−𝛼𝛽+𝜇_𝑖+1, (15) which not only gives the upper bound for the PMHSS method and complements the theoretical results of paper [24], but also includes the special case in [2] as follows: 𝛿=max_{𝜆𝑖∈𝜆(}𝑊₎𝛼²+𝜆²_𝑖𝛼+𝜆_𝑖, (16) when 𝛼=𝛽.

(ii) Based on (9), it is obvious that the convergence rate not only depends on the choices of two parameters 𝛼 and 𝛽, but also depends on the choices of two auxiliary matrices 𝑃₁ and 𝑃₂. Obviously, the efficiency of GPMHSS is best when one can obtain two iteration parameters and auxiliary matrices to minimize an upper bound on the spectral radius of the iteration matrix. That is to say, it is necessary to discuss the choices of 𝛼, 𝛽, 𝑃₁, and 𝑃₂. In fact, under certain conditions, in theory, some optimal iteration parameters 𝛼 and 𝛽 can be derived. For example, if 𝑃₁=0, then the GPMHSS method reduces to the LPMHSS method and the corresponding optimal parameter is 𝜆²_min/𝜇_max, where 𝜆_min and 𝜇_max, respectively, denote the smallest eigenvalue of 𝑃₂⁻¹𝑊 and the largest eigenvalue of 𝑃₂⁻¹𝑇 [22]. If 𝑃₁=𝑃₂ and 𝛼=𝛽, then the GPMHSS method reduces to the PMHSS method and the corresponding optimal parameter is 𝜆_min𝜆_max, where 𝜆_min and 𝜆_max, respectively, denote the smallest and largest eigenvalues of 𝑃₂⁻¹𝑊 [3]. Whereas, if 𝑃₁=𝑃₂ and 𝛼≠𝛽, in theory, one cannot find the optimal values of the iteration parameters 𝛼 and 𝛽 in general. With respect to this point, one can see [24] for more details. Based on this fact, in theory, a conclusion may be obtained, that is, one cannot find the optimal values of the iteration parameters 𝛼 and 𝛽 for GPMHSS in general. But even so, the optimal values of the iteration parameters 𝛼 and 𝛽 for GPMHSS may be obtained experimentally (to see the fourth section). In practice, to further improve the efficiency of the GPMHSS method, it is desirable to determine or find a good estimate of the optimal parameter that minimizes the convergence factor. In fact, to find the actual optimal estimates 𝛼 and 𝛽 for the GPMHSS method is a hard task because its solution strongly depends on the particular structures and properties of the coefficient matrix 𝐴, as well as the splittings matrices 𝑊 and 𝑇, and needs further in-depth study, both from theory and computation point of view. When the optimal parameter cannot be derived in theory, the value of parameter 𝛼=𝛽 with 𝑃₁=𝑃₂ is selected by the statement on the choice of the iteration parameter in [25], that is to say, experience suggests that in most applications and for an appropriate scaling of the problem, a “small” value of parameter (usually between 0.01 and 0.5) may give good results. Whereas choosing a parameter so as to minimize the spectral radius of the iteration matrix is not necessarily the best choice when the algorithm is used as a preconditioner for a Krylov subspace method. Remarkably, it can be shown that, for certain problems, the alternating iteration results in an ℎ-independent preconditioner for GMRES when parameter is chosen sufficiently small, correspond to a spectral radius very close to 1. In fact, if we define the optimal value of parameter as the one that minimizes the total amount of work needed to compute an approximate solution, this will not necessarily be the same as the parameter that minimizes the number of (outer) iterations. Overall, the analytic determination of such an optimal value for parameter appears to be daunting. With respect to the choices of two auxiliary matrices 𝑃₁ and 𝑃₂, 𝑃₁ and 𝑃₂ may have a better degree of freedom. For instance, if 𝑊+𝑇 is positive definite, then 𝑃₁=𝑇 and 𝑃₂=𝑊 with 𝛼>0 would be a natural choice for the preconditioners.

(iii) When 𝑃₁≠𝑃₂ and 𝛼=𝛽, the upper bound in (9) results in 𝛿=max_{𝛾𝑖∈𝛾}𝑃𝑃^𝑇𝛾_𝑖−1max_{𝜆𝑖∈𝜆}𝑊𝛼𝛼+𝜆_𝑖+max_{𝜆𝑖∈𝜆}𝑊𝛼²+𝜆²_𝑖𝛼+𝜆_𝑖⋅max_{𝛾𝑖∈𝜆}𝑃𝑃^𝑇1𝛾_𝑖−1max_{𝜇𝑖∈𝜇}𝑇𝛼𝛼+𝜇_𝑖+1. (17)

The approach to minimize the upper bound is very important in theoretical viewpoint. However, it is not practical since the corresponding spectral radius of the iteration matrix 𝑀_𝛼,𝛽 is not optimal. How to choose the suitable preconditioners and parameters for practical problem is still a great challenge.

3. The IGPMHSS Iteration

In the GPMHSS method, it is required to solve two systems of linear equations whose coefficient matrices are 𝛼𝑃₁+𝑊 and 𝛽𝑃₂+𝑇, respectively. However, this may be very costly and impractical in actual implementations. To overcome this disadvantage and improve the computational efficiency of the GPMHSS iteration method, we propose to solve the two subproblems iteratively [21, 26], which leads to the inexact GPMHSS (IGPMHSS) iterative scheme. Its convergence can be shown in a similar way to that of the IHSS iteration method, using Theorem 3.1 of [21]. Since 𝛼𝑃₁+𝑊 and 𝛽𝑃₂+𝑇 are symmetric positive definite, there can employ some Krylov subspace methods (such as CG) to gain its solution easily. Of course, if good preconditioners for matrices 𝛼𝑃₁+𝑊 and 𝛽𝑃₂+𝑇 are available, we can use the preconditioned conjugate gradient (PCG) method instead of CG for the two inner systems, this yields a better performance of IGPMHSS method. If either 𝛼𝑃₁+𝑊 or 𝛽𝑃₂+𝑇 (or both) is Toeplitz, we can use fast algorithms for solution of the corresponding subsystems [27]. Here, just like the IHSS iteration, the IGPMHSS iterative scheme is presented in Algorithm 1.

𝑘=0; while (not convergent)  𝑟(𝑘)=𝑏−𝐴𝑥(𝑘);  approximately solve (𝛼𝑃1+𝑊)𝑧(𝑘)=𝑟(𝑘) by employing CG method,  such that the residual 𝑝(𝑘)=𝑟(𝑘)−(𝛼𝑃1+𝑊)𝑧(𝑘) of the iteration satisfies  ‖𝑝(𝑘)‖2≤𝜂𝑘‖𝑟(𝑘)‖2;  𝑥(𝑘+1/2)=𝑥(𝑘)+𝑧(𝑘);  𝑟(𝑘+1/2)=𝑏+𝑖𝐴𝑥(𝑘+1/2);  approximately solve (𝛽𝑃2+𝑇)𝑧(𝑘+1/2)=𝑟(𝑘+1/2) by employing CG method,  such that the residual 𝑞(𝑘+1/2)=𝑟(𝑘+1/2)−(𝛽𝑃2+𝑇)𝑧(𝑘+1/2) of the  iteration satisfies ‖𝑞(𝑘+1/2)‖2≤𝜏𝑘‖𝑟(𝑘+1/2)‖2;  𝑥(𝑘+1)=𝑥(𝑘+1/2)+𝑧(𝑘+1/2);  𝑘=𝑘+1; end

Algorithm 1 

Remark 3. It is not difficult to find the fact that if the inner systems can be solved exactly with 𝜂_𝑘=0 and 𝜏_𝑘=0, then the IGPMHSS iteration reduces to the GPMHSS iteration. In fact, to guarantee the convergence of the IGPMHSS iteration, it is not necessary for 𝜂_𝑘 and 𝜏_𝑘 to go to zero when 𝑘 increases.
To derive the convergence properties for IGPMHSS iteration, the following lemma is required, which is presented by Bai et al. [21]. Here ‖|𝑥|‖_𝑀=‖𝑀𝑥‖ for any 𝑥∈ℂ^𝑛, which immediately induces the matrix norm ‖|𝑋|‖_𝑀=‖𝑀𝑋𝑀⁻¹‖ for any 𝑋∈ℂ^𝑛×𝑛.

Lemma 4. Let 𝐴∈ℂ^𝑛×𝑛, 𝐴=𝑀_𝑖−𝑁_𝑖(𝑖=1,2) be two splittings of 𝐴. If {𝑥^(𝑘)} is an iteration sequence defined as follows: 𝑥^(𝑘+1/2)=𝑥^(𝑘)+𝑧^(𝑘)𝑤𝑖𝑡ℎ𝑀₁𝑧^(𝑘)=𝑟^(𝑘)+𝑝^(𝑘), (18) satisfying ‖𝑝^(𝑘)‖₂≤𝜂_𝑘‖𝑟^(𝑘)‖₂, where 𝑟^(𝑘)=𝑏−𝐴𝑥^(𝑘), and 𝑥^(𝑘+1)=𝑥^(𝑘+1/2)+𝑧^(𝑘+1/2)𝑤𝑖𝑡ℎ𝑀₂𝑧^(𝑘+1/2)=𝑟^(𝑘+1/2)+𝑞^(𝑘+1/2), (19) satisfying ‖𝑞^(𝑘+1/2)‖₂≤𝜏_𝑘‖𝑟^(𝑘+1/2)‖₂, where 𝑟^(𝑘+1/2)=𝑏−𝐴𝑥^(𝑘+1/2), then {𝑥^(𝑘)} is of the form 𝑥^(𝑘+1)=𝑀₂⁻¹𝑁₂𝑀₁⁻¹𝑁₁𝑥^(𝑘)+𝑀₂⁻¹𝐼+𝑁₂𝑀₁⁻¹𝑏+𝑀₂⁻¹𝑁₂𝑀₁⁻¹𝑝^(𝑘)+𝑞^(𝑘+1/2). (20) Moreover, if 𝑥_∗∈ℂ^𝑛 is the exact solution of the system of linear equations (1), then we have 𝑥^(𝑘+1)−𝑥_∗𝑀2≤𝜁+𝜇𝜃𝜂_𝑘+𝜃𝜌+𝜃𝜈𝜂_𝑘𝜏_𝑘𝑥^(𝑘)−𝑥_∗𝑀2,𝑘=1,2,…, (21) where 𝜁=𝑁₂𝑀₁⁻¹𝑁₁𝑀₂⁻¹₂,𝜌=𝑀₂𝑀₁⁻¹𝑁₁𝑀₂⁻¹₂,𝜇=𝑁₂𝑀₁⁻¹₂,𝜃=𝐴𝑀₂⁻¹₂,𝜈=𝑀₂𝑀₁⁻¹₂. (22) In particular, if 𝜁+𝜇𝜃𝜂_max+𝜃𝜌+𝜃𝜈𝜂_max𝜏_max<1, (23) then the iteration sequence {𝑥^(𝑘)} converges to 𝑥_∗∈ℂ^𝑛, where 𝜂_max=max_𝑘{𝜂_𝑘} and 𝜏_max=max_𝑘{𝜏_𝑘}.

Applying Lemma 4, the following theorem gives a convergence analysis of the IGPMHSS iteration.

Let 𝛿_𝛼,𝛽=𝛼𝑃₁−𝑖𝑇𝛼𝑃₁+𝑊⁻¹𝛽𝑃₂+𝑖𝑊𝛽𝑃₂+𝑇⁻¹₂. (24) Then we have the following result.

Theorem 5. Let 𝑃₁∈ℝ^𝑛×𝑛 and 𝑃₂∈ℝ^𝑛×𝑛 be two symmetric positive definite matrices. Let 𝐴=𝑊+𝑖𝑇∈ℂ^𝑛×𝑛, with 𝑊∈ℝ^𝑛×𝑛 and 𝑇∈ℝ^𝑛×𝑛 symmetric positive definite and symmetric positive semidefinite, respectively, and let 𝛼 be a nonnegative constant and let 𝛽 be a positive constant. If {𝑥^(𝑘)} is an iteration sequence generated by the IGPMHSS iteration method (Algorithm 1) and if 𝑥_∗∈ℂ^𝑛 is the exact solution of the system of linear equations (1), then it holds that 𝑥^(𝑘+1)−𝑥_∗≤𝛿_𝛼,𝛽+𝜇𝜃𝜂_𝑘+𝜃𝜌+𝜃𝜐𝜂_𝑘𝜏_𝑘𝑥^(𝑘)−𝑥_∗,𝑘=1,2,…, (25) where 𝜌=𝛽𝑃₂+𝑇𝛼𝑃₁+𝑊⁻¹𝛼𝑃₁−𝑖𝑇𝛽𝑃₂+𝑇⁻¹₂,𝜃=𝐴𝛽𝑃₂+𝑇⁻¹₂,𝜇=𝛽𝑃₂+𝑖𝑊𝛼𝑃₁+𝑊⁻¹₂,𝜐=𝛽𝑃₂+𝑇𝛼𝑃₁+𝑊⁻¹₂. (26) In particular, if 𝛿_𝛼,𝛽+𝜇𝜃𝜂_max+𝜃𝜌+𝜃𝜈𝜂_max𝜏_max<1, (27) then the iteration sequence {𝑥^(𝑘)} converges to 𝑥_∗, where 𝜂_max=max_𝑘{𝜂_𝑘} and 𝜏_max=max_𝑘{𝜏_𝑘}.

In fact, replacing 𝑀_𝑖 and 𝑁_𝑖(𝑖=1,2) in Lemma 4 with 𝑀₁=𝛼𝑃₁+𝑊, 𝑁₁=𝛼𝑃₁−𝑖𝑇, 𝑀₂=𝛽𝑃₂+𝑇, and 𝑁₂=𝛽𝑃₂+𝑖𝑊, we straightforwardly obtain the proof of Theorem 5.

Theorem 5 shows that the choices of the tolerances 𝜂_𝑘 and 𝜏_𝑘 are to make the IGPMHSS convergent. In fact, as mentioned previously, to get the convergence of the IGPMHSS iteration, the tolerances 𝜂_𝑘 and 𝜏_𝑘 are not required to approach zero as 𝑘 increases. However, the optimal tolerances 𝜂_𝑘 and 𝜏_𝑘 are not easy to analyze.

4. Numerical Experiments

In this section, we give some numerical experiments to demonstrate the performance of the GPMHSS and IGPMHSS methods for solving the linear system (1)-(2). Numerical comparisons with the MHSS and PMHSS methods are also presented to show the advantage of the GPMHSS method.

In our implementations, the initial guess is chosen to be 𝑥⁽⁰⁾=0 and the stopping criteria for outer iterations (when MHSS, PMHSS, and GPMHSS methods are used as solvers) are 𝑏−𝐴𝑥^(𝑘)₂‖𝑏‖₂≤10⁻⁶. (28)

The preconditioner 𝑃 used in PMHSS method is chosen to be 𝑃=𝑊. For the sake of comparing, the corresponding preconditioners 𝑃₁ and 𝑃₂ used in GPMHSS method are chosen to be 𝑃₁=𝑃₂=𝑊. Similarly, if the preconditioner 𝑃 used in PMHSS method is chosen to be 𝑃=𝑇, then corresponding preconditioners 𝑃₁ and 𝑃₂ used in GPMHSS method are chosen to be 𝑃₁=𝑃₂=𝑇. Since the numerical results in [2, 3] show that the PMHSS method outperforms the MHSS and HSS methods when they are employed as preconditioners for the GMRES method or its restarted variants [28], we just examine the efficiency of the GPMHSS method as a solver for solving complex symmetric linear system (1)-(2) by comparing the iteration numbers (denoted as IT) and CPU times (in seconds, denoted as CPU) of the GPMHSS method with those of the MHSS and PMHSS methods.

Example 6 (see [2, 3]). The complex symmetric linear system (1)-(2) is of the following form: 𝐴𝑥≡(𝑊+𝑖𝑇)𝑥=𝑏, (29) with 𝑇=𝐼⊗𝑉+𝑉⊗𝐼,𝑊=10𝐼⊗𝑉_𝑐+𝑉_𝑐⊗𝐼+9𝑒₁𝑒^𝑇_𝑚+𝑒_𝑚𝑒^𝑇₁⊗𝐼, (30)(31) where 𝑉= tridag (−1,2,−1)∈ℝ^𝑚×𝑚, 𝑉_𝑐=𝑉−𝑒₁𝑒^𝑇_𝑚−𝑒_𝑚𝑒^𝑇₁∈ℝ^𝑚×𝑚, and 𝑒₁ and 𝑒_𝑚 are the first and the last unit vectors in ℝ^𝑚, respectively. We take the right-hand side vector 𝑏 to be 𝑏=(1+𝑖)𝐴𝟏, with 1 being the vector of all entries equal to 1. Here 𝑇 and 𝑊 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 ℎ=1/(𝑚+1).

As is known, the spectral radius of the iteration matrix may be decisive for the convergence of the iteration method. The spectral radius corresponding to the iteration method is necessary to consider. The comparisons of the spectral radius of the three different iteration matrices derived by MHSS, PMHSS, and GPMHSS methods with different mesh-sizes are performed in Tables 1 and 2. In Tables 1 and 2, we use the optimal values of the parameters 𝛼 and 𝛽, denoted by 𝛼^∗_𝑀 for MHSS method, 𝛼^∗_𝑃 for PMHSS method, and 𝛼^∗_𝐺, 𝛽^∗_𝐺 for GPMHSS method. These parameters are obtained experimentally with the least spectral radius for the iteration matrices of the three methods.

Table 1 

The experimentally optimal parameters and the spectral radii for the iteration matrices of MHSS, PMHSS, and GPMHSS with 𝑃=𝑃₁=𝑃₂=𝑊.

Table 2 

The experimentally optimal parameters and the spectral radii for the iteration matrices of MHSS, PMHSS, and GPMHSS with 𝑃=𝑃₁=𝑃₂=𝑇.

From Tables 1 and 2, one can see that with the mesh-size 𝑚 increasing, the trend of the experimentally optimal parameter of the MHSS method decreases. Whereas, with the mesh-size 𝑚 increasing, the experimentally optimal parameter of the PMHSS method may be changeless. Fixing the parameter 𝛼 in GPMHSS method, the trend of the experimentally optimal parameter 𝛽 of the MHSS method decreases. No matter what 𝑃=𝑃₁=𝑃₂=𝑊 or 𝑃=𝑃₁=𝑃₂=𝑇, we observe that the optimal spectral radius for the iteration matrices of the three methods grows with problem size, and the optimal spectral radius of GPMHSS method is still smaller than those of MHSS and PMHSS methods. In this case, the GPMHSS method may outperform the MHSS and PMHSS methods. To this end, we need to examine the efficiencies of the MHSS, PHSS, and GPMHSS methods for solving the systems of linear equations 𝐴𝑥=𝑏, where 𝐴 is described above.

In Tables 3 and 4, we list the iteration numbers and computational times for MHSS, PMHSS, and GPMHSS iteration methods by using the optimal parameters in Tables 1 and 2. In Tables 3 and 4, “RES” denotes the relative residual error.

Table 3 

RES, CPUs, and IT for MHSS, PMHSS, and GPMHSS with 𝑃=𝑃₁=𝑃₂=𝑊.

Table 4 

RES, CPUs, and IT for MHSS, PMHSS, and GPMHSS with 𝑃=𝑃₁=𝑃₂=𝑇.

From Tables 3 and 4, we see that GPMHSS iteration method is the best among three methods in terms of the iteration numbers and computational time, and the PMHSS scheme requires fewer iteration numbers than the MHSS scheme. For the MHSS and GPMHSS methods, the number of iterations grows with the problem size. For the PMHSS method, its iteration numbers are relatively stable. That is to say, the PMHSS method does not have any growth in iterations numbers by increasing grid dimension. The presented results in Tables 3 and 4 show that in all cases GPMHSS is superior to another two methods in terms of the CPU time. That is to say, compared with the MHSS and PMHSS methods, the GPMHSS method may be given priority under certain conditions. In fact, the speed-up of GDMHSS/PMHSS with respect to MHSS is quite noticeable, where we define it by speed-up=CPUofHSSmethod.CPUofGPMHSS(orPMHSS)method (32) Obviously, the efficiency of GPMHSS is superior to that of MHSS and PMHSS.

In particular, here we test the efficiency of the GPMHSS method when 𝑃₁=𝑇 and 𝑃₂=𝑊. In this case, some numerical results are obtained in Table 5. In Table 5, it is not difficult to find that the spectral radius grows with problem size. Simultaneously, the number of iterations grows with the problem size. Numerical results in Table 5 show that, under certain conditions, the GPMHSS method is feasible and efficient, compared with the MHSS and PMHSS methods.

Table 5 

Numerical results for GPMHSS with 𝑃₁=𝑇 and 𝑃₂=𝑊.

As already noted, in the two half steps of the GPMHSS iteration, it is necessary to solve two systems of linear equations, whose coefficient matrices are 𝛼𝑃₁+𝑊 and 𝛽𝑃₂+𝑇, respectively. This can be very costly and impractical in actual implementations. We use the IGPMHSS method to solve the systems of linear equations (1)-(2) in the actual implementations. That is, it is necessary to solve two systems of linear equations with 𝛼𝑃₁+𝑊 and 𝛽𝑃₂+𝑇 by using the IGPMHSS iteration. It is easy to know that 𝛼𝑃₁+𝑊 and 𝛽𝑃₂+𝑇 are symmetric positive definite. So, solving the above two subsystems, there can employ the CG method.

In our computations, the inner CG iteration is terminated if the current residual of the inner iterations satisfies 𝑝^(𝑗)₂𝑟^(𝑘)₂≤0.1𝜏^(𝑘),𝑞^(𝑗)₂𝑟^(𝑘)₂≤0.1𝜏^(𝑘), (33) (cf. Algorithm 1) where 𝑝^(𝑗) and 𝑞^(𝑗) are, respectively, the residuals of the 𝑗th inner CG for 𝛼𝑃₁+𝑊 and 𝛽𝑃₂+𝑇. 𝑟^(𝑘) is the 𝑘th outer IGPMHSS iteration, and 𝜏 is a tolerance.

Some results are listed in Tables 6 and 7, which are the numbers of outer IGPMHSS iteration (it.s), the average numbers (avg1) of inner CG iteration for 𝛼𝑃₁+𝑊, and the average numbers (avg2) of CG iteration for 𝛽𝑃₂+𝑇.

Table 6 

Convergence results for IGPMHSS with 𝑃₁=𝑃₂=𝑊.

Table 7 

Convergence results for IGPMHSS with 𝑃₁=𝑃₂=𝑇.

In our numerical computations, it is easy to find the fact that the choice of 𝜏 is important to the convergence rate of the IGPMHSS method. According to Tables 6 and 7, the iteration numbers of the IGPMHSS method generally increases when 𝜏 decreases. This increase can probably be eliminated using a suitable preconditioner. It is noted that no preconditioning is used for these inner iterations in our numerical computations.

5. Conclusion

In this paper, we have generalized the MHSS method into the GPMHSS method for a class of complex symmetric linear systems. Theoretical analysis shows that for any initial guess the GPMHSS method converges to the unique solution of the linear system for a wide range of the parameters. Then, an inexact version of GPMHSS (IGPMHSS) has been presented and implemented for saving the computational cost. Numerical experiments show that the GPMHSS and IGPMHSS methods are efficient and competitive.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgments

The authors would like to thank the anonymous referees for their helpful suggestions, which greatly improved the paper. This research was supported by the NSFC (no. 11301009), the Science and Technology Development Plan of Henan Province (no. 122300410316), and the Natural Science Foundations of Henan Province (no. 13A110022).