A Polynomial Preconditioned Global CMRH Method for Linear Systems with Multiple Right-Hand Sides

Abstract

The restarted global CMRH method (Gl-CMRH(m)) (Heyouni, 2001) is an attractive method for linear systems with multiple right-hand sides. However, Gl-CMRH(m) may converge slowly or even stagnate due to a limited Krylov subspace. To ameliorate this drawback, a polynomial preconditioned variant of Gl-CMRH(m) is presented. We give a theoretical result for the square case that assures that the number of restarts can be reduced with increasing values of the polynomial degree. Numerical experiments from real applications are used to validate the effectiveness of the proposed method.

1. Introduction

We consider the linear systems of the form 𝐴𝑋=𝐵, (1) where 𝐴 is a nonsingular and sparse matrix of order 𝑛, 𝐵∈ℝ^𝑛×𝑠 with usually 𝑠≤𝑛. Such a situation arises from, for instance, wave-scattering problems, image restorations, recursive least squares computations, numerical methods for integral equations, and structural mechanics problems; see [1, 2] and references therein.

Numerical solvers for (1) can be roughly divided into two categories. The first is the direct method; for instance, the LU factorization is competent since 𝐴 is only factorized once to recast (1) into a triangular system which is easy to solve. However, if the coefficient matrix 𝐴 is large and sparse or sometimes not readily available, then iterative solvers may become the only choice and possibly fall into the following three classes.

The first class is the block method. For symmetric problems, the first block methods are due to O’Leary, including the block conjugate gradient method and block biconjugate gradient method [3]. For nonsymmetric cases, the block generalized minimal residual method [4–6], the block BiCGSTAB method [7], the block quasiminimal residual method [8], the block least squares method [9], the block Lanczos method [10], and the block IDR(𝑠) method [11] have been proposed recently. In general, the block solvers are more suitable for dense systems with precondition.

The second class is the seed method. The main idea of this kind of methods is briefed below. We first select a single system (seed system) and develop the corresponding Krylov subspace. Then we project all the residuals of the other systems onto the same Krylov subspace to find new approximations as initial guess; see [2, 12, 13] for details.

The last class is the global method. To our knowledge, the term global is at least due to Saad [14, Chapter 10] and has been further populated by Jbilou et al. [15] with the global FOM and global GMRES methods for matrix equations. Following the work [15], many other global methods have been developed, including, to name just a few of them, the global BiCG and global BiCGSTAB methods [16, 17], the global Hessenberg and global CMRH (changing minimal residual method based on the Hessenberg process) methods [18] and their weighted variants [19], the skew-symmetric methods [20], and the global SCD method [21]. Generally, the global methods are more appropriate for large and sparse systems.

It is well known that the performance of the above Krylov subspace methods can be reinforced with a suitable preconditioner [14] or through effective matrix splitting techniques [22, 23]. In this paper, we are interested in preconditioning the global methods. Specifically, we aim at improving the convergence behavior of the restarted global CMRH method (Gl-CMRH(𝑚)) [18], which is originally proposed to reduce the increasing storage requirement in its full version. However, because of the use of a small subspace (say 𝑚≪𝑛), Gl-CMRH(𝑚) is likely to slow down or even stalls out. Heyouni and Essai give a weighted version of Gl-CMRH(𝑚) (WGl-CMRH(𝑚)) [19] to alleviate such disadvantage. Instead, we propose a different approach, that is, by polynomial preconditioner to improve Gl-CMRH(𝑚) in this paper.

The remainder of this work is organized as follows. In Section 2, we first recall some notations and properties of the global method, and then we sketch the Gl-CMRH(𝑚) method. In Section 3, we construct the polynomial preconditioner tailored to Gl-CMRH(𝑚) by exploiting the relation between the Krylov matrix and the global basis. For square right-hand side matrices, we also give a theoretical result that justifies the use of the proposed polynomial preconditioner. In Section 4, several numerical examples are employed to substantiate the effectiveness of the proposed method. Some concluding remarks and potential future work are briefed in the last section.

2. Notations and the Global CMRH Method

In this section, we first give some notations and properties used in the global methods, which will henceforth be adopted extensively in deriving the main results. Then we present a brief introduction of the Gl-CMRH(𝑚) method [18]. More details about the global methods can be found, for instance, in [15, 18, 19].

2.1. Notations and Properties

Throughout this paper, the following notations will be used. The norms ‖⋅‖₂ and ‖⋅‖_𝐹 represent the vector 2-norm and matrix 𝐹-norm, respectively. Let 𝕄 be the set of 𝑛×𝑠 rectangular matrices. If 𝑋∈𝕄, then 𝑋^𝑇 stands for its transpose. For a square matrix 𝐴, 𝐴⁻¹ indicates the inverse of 𝐴 if existed. Unless otherwise stated, subscripts denote the corresponding iteration step; for example, 𝑋_𝑘 denotes the 𝑘th iterate of the matrix (vector) 𝑋. Moreover, the (𝑖,𝑗) entries of matrices 𝑌 and 𝑋_𝑘 are denoted by (𝑌)_𝑖,𝑗 and (𝑋_𝑘)_𝑖,𝑗, respectively. If a column or a row of a matrix is invoked, then we denote it in a dot format; that is, (𝑋_𝑘)_⋅,𝑗 and (𝑋_𝑘)_𝑖,⋅ mean correspondingly the 𝑗th column and the 𝑖th row of 𝑋_𝑘. Besides, (𝑋_𝑘)_{𝑖∶𝑗,𝑠∶𝑡} extracts the submatrix from 𝑖 to 𝑗 rows and from 𝑠 to 𝑡 columns of 𝑋_𝑘.

Next we present some notations and basic properties used in the global methods [15]. Given the 𝑛×𝑚𝑠 block matrix 𝒱_𝑚=[𝑉₁,𝑉₂,…,𝑉_𝑚], where 𝑉_𝑖∈𝕄, 𝑖=1,2,…,𝑚, then we define the matrix product ∗ as 𝒱_𝑚∗𝑓=^𝑚_𝑖=1𝑓_𝑖𝑉_𝑖, (2) where 𝑓∈ℝ^𝑚. For any matrix 𝐻∈ℝ^𝑚×𝑚, we define analogously the ∗ product by 𝒱_𝑚∗𝐻=𝒱_𝑚∗(𝐻)_⋅,1,𝒱_𝑚∗(𝐻)_⋅,2,…,𝒱_𝑚∗(𝐻)_⋅,𝑚. (3)

It can be verified that such matrix product satisfies the following properties: 𝒱_𝑚∗𝑓+𝑔=𝒱_𝑚∗𝑓+𝒱_𝑚,𝒱∗𝑔_𝑚∗𝐻∗𝑓=𝒱_𝑚∗𝐻𝑓, (4) where 𝑓,𝑔∈ℝ^𝑚.

2.2. The Global CMRH Method

The Gl-CMRH method [18] is an efficient extension of the CMRH method [24] for solving (1). It is based on the global Hessenberg process [18]. As for the numerical performance, Gl-CMRH is in general competitive with the classic global GMRES method (Gl-GMRES) [15].

Now we give a brief sketch of Gl-CMRH. Let 𝑋₀∈𝕄 be the initial guess of (1) with the associated residual matrix 𝑅₀=𝐵−𝐴𝑋₀. The 𝑚th iteration 𝑋_𝑚 is searched in the affine subspace 𝑋₀+𝒦_𝑚(𝐴,𝑅₀); that is, 𝑋_𝑚−𝑋₀=𝑊_𝑚∈𝒦_𝑚(𝐴,𝑅₀), where 𝑊_𝑚 is the 𝑚th correction matrix. The matrix Krylov subspace is defined as 𝒦_𝑚(𝐴,𝑅₀)=span{𝑅₀,𝐴𝑅₀,…,𝐴^𝑚−1𝑅₀}. Using the basis 𝒱_𝑚 of 𝒦_𝑚(𝐴,𝑅₀) given by the global Hessenberg process [18], we get 𝑋_𝑚=𝑋₀+𝒱_𝑚∗𝑦_𝑚, (5) where 𝑦_𝑚∈ℝ^𝑚. The global Hessenberg process in [18] also yields 𝐴𝒱_𝑚=𝒱_𝑚+1∗𝐻_𝑚=𝒱_𝑚∗𝐻_𝑚+𝐻_{𝑚𝑚+1,𝑚}𝟎,𝟎,…,𝟎,𝑉_𝑚+1, (6) where 𝐻_𝑚 is an (𝑚+1)×𝑚 upper Hessenberg matrix, 𝐻_𝑚 is obtained by deleting the last row of 𝐻_𝑚, and 0 is the zero matrix in 𝕄. Thus it follows immediately that 𝑅_𝑚=𝑅₀−𝐴𝒱_𝑚∗𝑦_𝑚=𝒱_𝑚+1∗𝛽𝑒₁^(𝑚+1)−𝐻_𝑚𝑦_𝑚, (7) where 𝛽=max_{1≤𝑖≤𝑛,1≤𝑗≤𝑠}{|(𝑅₀)_𝑖,𝑗|} and 𝑒₁^(𝑚+1)=[1,0,…,0]^𝑇∈ℝ^𝑚+1. To obtain the vector 𝑦_𝑚, a restriction is imposed on the Gl-CMRH method; that is, 𝑅_𝑚𝐹=min_{𝑊∈𝒦𝑚𝐴,𝑅0}𝑅₀−𝐴𝑊_𝐹. (8)

Relations (7) and (8) yield 𝑦_𝑚=argmin_{𝑦∈ℝ𝑚}𝒱_𝑚+1∗𝛽𝑒₁^(𝑚+1)−𝐻_𝑚𝑦_𝐹. (9)

Instead of solving (9), which requires 𝒪(𝑛𝑚²) operations and 𝒪(𝑛𝑚) storage, we solve a smaller problem min_{𝑦∈ℝ𝑚}𝛽𝑒₁^(𝑚+1)−𝐻_𝑚𝑦₂, (10) which leads to 𝑦_𝑚=𝛽(𝐻^𝑇_𝑚𝐻_𝑚)⁻¹𝐻^𝑇_𝑚𝑒₁^(𝑚+1) by assuming that 𝐻_𝑚 is of full rank. From (5) and (10), the 𝑚th iterate 𝑋_𝑚 can be updated by 𝑋_𝑚=𝑋₀+𝛽𝒱_𝑚∗𝐻^𝑇_𝑚𝐻_𝑚⁻¹𝐻^𝑇_𝑚𝑒₁^(𝑚+1). (11)

As in the Gl-GMRES method [15], a restarting strategy is used to address the problem that the computational and storage requirements increase with iterations. Algorithm 1 gives a framework of the restarted version of Gl-CMRH (Gl-CMRH(𝑚)). We refer to [18, 19] for elaborate explanation for the Gl-CMRH method.

(1) Choose an initial guess 𝑋0, the restarting frequency 𝑚 and the tolerance tol. Compute 𝑅0=𝐵−𝐴𝑋0. Determine 𝑖0,𝑗0 such that (𝑅0)𝑖0,𝑗0=max1≤𝑖≤𝑛,1≤𝑗≤𝑠{(𝑅0)𝑖0,𝑗0}. Set 𝛽=(𝑅0)𝑖0,𝑗0, 𝑉1=𝑅0/𝛽, 𝑝1,1=𝑖0 and 𝑝1,2=𝑗0. (2) Construct the matrix basis 𝒱𝑚=[𝑉1,𝑉2,…,𝑉𝑚] and 𝐻𝑚 by the global Hessenberg process [18]. (3) Solve 𝑦𝑚 by (10) and update 𝑋𝑚 by (5). (4) Compute 𝑅𝑚=𝐵−𝐴𝑋𝑚. If 𝑅𝑚𝐹/𝑅0𝐹≤tol, then stop; otherwise set 𝑋0=𝑋𝑚, 𝑅0=𝑅𝑚.  Choose 𝑖0,𝑗0 such that (𝑅0)𝑖0,𝑗0=max1≤𝑖≤𝑛,1≤𝑗≤𝑠{(𝑅0)𝑖0,𝑗0}.  Set 𝛽=(𝑅0)𝑖0,𝑗0, 𝑉1=𝑅0/𝛽, 𝑝1,1=𝑖0 and 𝑝1,2=𝑗0. Go to Step >2.

Algorithm 1 

The Gl-CMRH(𝑚) method for 𝐴𝑋=𝐵 [18].

3. A New Polynomial Preconditioned Gl-CMRH(𝑚) Method

In many cases, the accuracy of Gl-CMRH(𝑚) is sufficient. Due to the limited dimension of the matrix Krylov subspace 𝒦_𝑚, however, Gl-CMRH(𝑚) may suffer from slow convergence or even stagnation in practice, just like in GMRES(𝑚) [25] and Gl-GMRES(𝑚) [15]. To remedy this drawback, some accelerating techniques are demanded, for instance, a weighting strategy exploited in [19]. Besides, a polynomial preconditioner can also be adapted to improve the convergence [26]. In this section we focus on constructing an efficient polynomial preconditioner pertinent to the Gl-CMRH(𝑚) method.

The essence of the polynomial preconditioned method is to devise a polynomial 𝑄(𝐴)≈𝐴⁻¹ such that an easier system 𝑄(𝐴)𝐴𝑋=𝑄(𝐴)𝐵 (12) is solved instead of solving the original system (1). In what follows, we obtain a polynomial preconditioner 𝑄(𝑥) by extracting some useful information from Gl-CMRH(𝑚).

Now suppose that the block Krylov matrix 𝐾_𝑘 is of the form 𝐾_𝑘=𝑅₀,𝐴𝑅₀,…,𝐴^𝑘−1𝑅₀, (13) where 𝑅₀ is the initial residual matrix. By comparing the last 𝑠 columns of the second equation in (6), we have 𝐴𝑉_𝑘=𝒱_𝑘∗𝐻_𝑘⋅,𝑘+𝐻_{𝑘𝑘+1,𝑘}𝑉_𝑘+1. (14)

The equality (14) can be rearranged as 𝑉_𝑘+1=𝐻_{𝑘𝑘+1,𝑘}⁻¹𝐴𝑉_𝑘−𝒱_𝑘∗𝐻_𝑘⋅,𝑘. (15)

Let us consider the relationship between 𝐾_𝑘 and the basis 𝒱_𝑘. Since 𝒱_𝑘 and 𝐾_𝑘 span the same space, it follows that 𝒱_𝑘=𝐾_𝑘∗𝑈_𝑘, (16) where 𝑈_𝑘 is an upper triangular matrix. The relation (16), however, does not shed too much light because how to compute 𝑈_𝑘 still remains unclear. Fortunately, an explicit recurrence for 𝑈_𝑘 can be derived in terms of 𝑈_𝑘−1 and 𝐻_𝑘−1. By combining (16) and (4), we get 𝒱_𝑘∗𝐻_𝑘⋅,𝑘=𝐾_𝑘∗𝑈_𝑘𝐻_𝑘⋅,𝑘=𝐾_𝑘,𝐴^𝑘𝑅₀∗𝑈_𝑘𝐻_𝑘⋅,𝑘0. (17)

Since 𝑉_𝑘=(𝒱_𝑘)_{⋅,(𝑘−1)𝑠+1∶𝑘𝑠}=𝐾_𝑘∗(𝑈_𝑘)_⋅,𝑘, then 𝐴𝑉_𝑘=𝐴𝑅₀,𝐴²𝑅₀,…,𝐴^𝑘𝑅₀∗𝑈_𝑘⋅,𝑘=𝐾_𝑘+1∗0𝑈_𝑘⋅,𝑘. (18)

Substituting (17) and (18) into (15) gives rise to 𝑉_𝑘+1=𝐻_{𝑘𝑘+1,𝑘}⁻¹𝐾_𝑘+1∗0𝑈_𝑘⋅,𝑘−𝑈_𝑘𝐻_𝑘⋅,𝑘0. (19)

Besides, the relation (16) gives 𝑉_𝑘+1=𝐾_𝑘+1∗(𝑈_𝑘+1)_⋅,𝑘+1. By combining it with (19), we obtain a recurrence for the (𝑘+1)st column of 𝑈_𝑘+1; that is, 𝑈_{𝑘+1⋅,𝑘+1}=𝐻_{𝑘𝑘+1,𝑘}⁻¹0𝑈_𝑘⋅,𝑘−𝑈_𝑘𝐻_𝑘⋅,𝑘0. (20)

Therefore, 𝑈_𝑘 in (16) can be updated recursively by (20). Recall that in (5) 𝑋_𝑘 is updated on the basis 𝒱_𝑘. Here we will show another way to update 𝑋_𝑘 which is based on the block Krylov matrix 𝐾_𝑘. It follows from (5), (16), and (4) that 𝑋_𝑘=𝑋₀+𝐾_𝑘∗𝑈_𝑘∗𝑦_𝑘=𝑋₀+𝐾_𝑘∗𝑈_𝑘𝑦_𝑘=𝑋₀+^𝑘−1_𝑖=0𝛼_𝑖𝐴^𝑖𝑅₀, (21) where 𝑈_𝑘𝑦_𝑘=[𝛼₀,𝛼₁,…,𝛼_𝑘−1]^𝑇 and 𝑦_𝑘 is solved from (10). Denote by 𝑄_𝑘−1 a polynomial in 𝐴 of degree 𝑘−1; that is, 𝑄_𝑘−1(𝐴)=∑^𝑘−1_𝑖=0𝛼_𝑖𝐴^𝑖. Hence (21) can be recast as 𝑋_𝑘=𝑋₀+𝑄_𝑘−1(𝐴)𝑅₀. (22)

The matrix polynomial 𝑄_𝑘−1(𝐴) in (22) can be regarded as the approximation to 𝐴⁻¹ in some sense. This is justified for the case 𝑛=𝑠 in (1) by the following result.

Theorem 1. Let 𝑄_𝑘−1(𝐴), 𝑋₀ and 𝑅₀=𝐵−𝐴𝑋₀ be the square matrices defined in (22), and lot 𝑋_∗ be the true solution of (1). Suppose that 𝑋_∗−𝑋₀ is nonsingular. Then one has 𝐼−𝑄_𝑘−1(𝐴)𝐴_𝐹≤𝐸_𝑘𝐹𝐸₀⁻¹_𝐹, (23) where 𝐸_𝑘∶=𝑋_∗−𝑋_𝑘 and 𝐸₀∶=𝑋_∗−𝑋₀.

Proof. The inequality (23) follows immediately from an arrangement of (22).

Remark 2. In (23), the term ‖𝐸_𝑘‖_𝐹 becomes smaller with growing 𝑘 and hence the upper bound diminishes correspondingly, which in turn implies that 𝑄_𝑘−1(𝐴) approximates 𝐴⁻¹ asymptotically. This justifies the use of the polynomial preconditioner. In general, (23) assures that the number of restarts will be reduced correspondingly with increasing 𝑘. Yet this does not necessarily mean that the CPU time will be reduced simultaneously since the time saved from the reduction of restarts may be offset by the extra time spent in constructing the polynomial. In practice, we are often more concerned with the CPU time than the restarting number. Therefore, we restrict ourselves to small values of 𝑘. For 𝑠<𝑛, an inequality similar to (23) is generally unavailable. Nevertheless, numerical examples seem to demonstrate that the asymptotical property of (23) is also shared by the case 𝑠<𝑛; see Example 3 in Section 4 for more discussions.

By putting all together, we propose the new polynomial preconditioned global CMRH method (PGl-CMRH(𝑚,deg)) that is shown in Algorithm 2.

(1) Input: 𝑚, deg and 𝑋0. (2) Phase I: Compute 𝑅0=𝐵−𝐴𝑋0. Determine 𝑖0,𝑗0 such that (𝑅0)𝑖0,𝑗0=max1≤𝑖≤𝑛,1≤𝑗≤𝑠|(𝑅0)𝑖,𝑗|. Set 𝛽=(𝑅0)𝑖0,𝑗0, 𝑉1=𝑅0/𝛽, 𝑝1,1=𝑖0 and 𝑝1,2=𝑗0. Let 𝑈1,1=(𝑅0)𝑖0,𝑗0 (the upper triangular matrix defined in (16)). (3) for  𝑘=1∶deg (4)  𝑀=𝐴𝑉𝑘 (5)  for 𝑗=1∶𝑘 (6)   (𝐻)𝑗,𝑘=(𝑀)𝑝𝑗,1,𝑝𝑗,2 (7)   𝑀=𝑀−(𝐻)𝑗,𝑘𝑉𝑗 (8)  end for (9)  Determine 𝑖0,𝑗0 such that (𝑀)𝑖0,𝑗0=max1≤𝑖≤𝑛,1≤𝑗≤𝑠(𝑀)𝑖,𝑗. (10)  Set (𝐻)𝑘+1,𝑘=(𝑀)𝑖0,𝑗0, 𝑉𝑘+1=𝑀/(𝐻)𝑘+1,𝑘, 𝑝𝑘+1,1=𝑖0 and 𝑝𝑘+1,2=𝑗0. (11)  Compute (𝑈𝑘+1).,𝑘+1 from (20). (12) end for (13) 𝑋deg=𝑋0+𝒱deg∗𝑦deg, where 𝑦deg=argmin𝑦∈ℝdeg(𝑅0)𝑖0,𝑗0𝑒1(deg+1)−𝐻deg𝑦2. (14) Construct the polynomial preconditioner 𝑄deg−1(𝐴) with its coefficients decided by entries of 𝑈deg𝑦deg. (15)  Phase II: Solve 𝑄deg−1(𝐴)𝐴𝑋=𝑄deg−1(𝐴)𝐵 using Algorithm 1.

Algorithm 2 

PGl-CMRH(𝑚,deg) method for 𝐴𝑋=𝐵.

4. Numerical Examples

In this section, we present some numerical experiments which are coded with MATLAB 7.8.0. For fair comparisons, some other global solvers mentioned earlier like Gl-CMRH(𝑚) [18], WGl-CMRH(𝑚) [19], and Gl-GMRES(𝑚) [15] have also been implemented. From now on we drop the parameters 𝑚 and deg in brackets without ambiguity. In all examples, we assume that 𝑋₀=𝟎. The terminating criterion for the 𝑘th iteration is tol=‖𝑅_𝑘‖_𝐹/‖𝑅₀‖_𝐹≤10⁻¹⁰. Though other alternatives are possible, we use 𝐷=√𝑛𝑠|(𝑅₀)_𝑖,𝑗|/‖𝑅₀‖_𝐹 as the the weighting matrix for WGl-CMRH, which is also preferred in [19]. The coefficient matrices 𝐴 in the first two examples are derived from the discretizations of the Poisson’s equation and convection-diffusion problems which occur frequently in applied science and engineering. The coefficient matrices 𝐴 in the third example are quoted from the Matrix Market [27].

Example 1. We consider the linear systems of (1) in which its coefficient matrix 𝐴 is obtained from the discretization of ℒ𝑢=𝑢_𝑥𝑥+𝑢_𝑦𝑦 (24) on the unit square (0,1)×(0,1) with 𝑢=0 on the boundary. It can be discretized through the centered difference scheme at the gird points (𝑥_𝑖,𝑦_𝑖) with 𝑥_𝑖=𝑖ℎ,𝑦_𝑗=𝑗ℎ, where the mesh size ℎ=1/(𝑁+1) for 𝑖,𝑗=0,…,𝑁+1. This yields a block tridiagonal matrix of size 𝑛=𝑁². The right-hand side matrix 𝐵 is chosen with entries uniformly distributed on [0,1]; see [14, Chapter 2] for more details about (24). Related parameters are given by 𝑠=2, 𝑚=20 and deg=5. The number of restarts and CPU time for matrices 𝐴 of different sizes are given in Table 1. As observed from Table 1, PGl-CMRH improves the original CMRH method by time ratios from 17.8% to 55.4%. Compared with WGl-CMRH, PGl-CMRH requires less number of restarts and CPU time to achieve the required accuracy. Note that WGl-CMRH does not speed up the convergence of Gl-CMRH. This indicates that a different weighting matrix should be used. To find the optimal weighting matrix, however, remains an open problem [19].

Table 1 

Number of restarts and CPU time (in brackets) for Example 1.

Example 2. Consider the linear systems of (1) where its coefficient matrix 𝐴 is obtained from the discretization of the three-dimensional convection-diffusion problem 𝒯𝑢=−𝑢_𝑥𝑥+𝑢_𝑦𝑦+𝑢_𝑧𝑧+𝑞𝑢_𝑥+𝑢_𝑦+𝑢_𝑧 (25) on the unit cube Ω=[0,1]×[0,1]×[0,1]. Here 𝑞 is a constant coefficient and (25) subjects to Dirichlet-type boundary conditions. This equation can be discretized by applying seven-point finite difference discretizations. For instance, we use the centered difference to the diffusive terms and the first-order upwind approximations to the convective terms. This approach yields a coefficient matrix 𝐴 of size 𝑛=𝑁³, where the equidistant mesh size ℎ=1/(𝑁+1) is used, and the natural lexicographic ordering is adopted to the unknowns; we refer to [22, Section 4] for more details. The right-hand side matrix 𝐵 is chosen with entries uniformly distributed on [0,1]. Here, 𝑠=2, 𝑚=15, and deg=5. The number of restarts and CPU time for 𝑞=0.1 and 𝑞=1 is given in Table 2. For this large problem, as expected, PGl-CMRH performs better than CMRH and other variants concerning CPU time.

Table 2 

Number of restarts and CPU time for Example 2; 𝑞=0.1 (top) and 𝑞=1 (bottom).

Example 3. In practice, the degree of the polynomial preconditioner 𝑄_deg−1 has a great impact on the numerical performance of PGl-CMRH. Thus it deserves our attention to investigate how to choose the “optimal” degree (if existed) for generic matrices. Nevertheless, theoretical analysis to this end can be very hard. Instead, we show empirically how to choose a range of degrees for the polynomial 𝑄_deg−1 such that PGl-CMRH at least yields a modest performance. To this end, we use ten unsymmetrical testing matrices from [27] and illustrate how PGl-CMRH performs for each matrix with deg varying from 2 to 15; see Figure 1. Some properties of these testing matrices are listed in Table 3. The right-hand side matrix 𝐵 is chosen with entries uniformly distributed on [0,1]. Since we are only concerned with the value of deg that makes PGl-CMRH performs stably with the shortest CPU time, we have normalized values of CPU time by dividing the maximum value of CPU time for a certain curve. Take the matrix pde2961 for example. The longest time is 4.2 seconds (with deg=3); then we divide all values of CPU time by 4.2 for pde2961 and plot the result in Figure 1. This approach facilitates our comparison since different curves become more clustered now. Some remarks can be made from Figure 1. First, the curves seem rather problem-dependent and are not necessarily nonincreasing with increasing values of deg; for instance, the curve of rdb2048l is rather irregular and hence unpredictable. However, this does not contradict Theorem 1 where it is stated that 𝑄_deg−1(𝐴) can approximate 𝐴⁻¹ better with growing values of deg. In other words, Theorem 1 explains theoretically that the total number of restarts will be reduced with increasing values of deg. However, this does not apply to the change of CPU time. In fact, it is likely that PGl-CMRH with high degree preconditioner takes more CPU time in generating the polynomial preconditioner (even with less number of restarts) and hence uses more time to converge than that of its low degree counterpart. Second, most curves locate the corresponding shortest CPU time point with deg between 2 and 10. This can be the first reason for favoring small values of deg. Finally, more rounding errors can be introduced in developing high-degree polynomial preconditioners from the numerical point of view. This is the second reason for the approval of low-degree polynomial preconditioners. Therefore it is useful to test with deg from 2 to 10. Under extreme situations, however, higher degree may be demanded if a low-degree preconditioner fails to bring the required accuracy.

Table 3 

Properties of testing matrices in Example 3.

Figure 1 

Normalized CPU time against deg (from 2 to 15) for ten testing matrices.

5. Conclusion

To remedy the slow convergence of the original Gl-CMRH(𝑚) method, a new variant of Gl-CMRH(𝑚) for linear systems with multiple right-hand sides is developed. The proposed method often yields better performance than its predecessor Gl-CMRH(𝑚) and other global variants in terms of CPU time. We show experimentally that polynomial preconditioners with degree lower than 10 should be considered if no prior knowledge is known.

Acknowledgments

The authors would like to thank Professor Jinyun Yuan and the referees for their valuable remarks that improved this paper. The work is supported by the National Natural Science Foundation (11371243), the Key Disciplines of Shanghai Municipality (S30104), the Innovation Program of Shanghai Municipal Education Commission (13ZZ068), and the Anhui Provincial Natural Science Foundation (1308085QF117).