Fast parameterized inexact Uzawa method for complex symmetric linear systems

Fast parameterized inexact Uzawa method for complex symmetric linear systemsqQing-Qing Zheng,Chang-Feng Ma ⇑School of Mathematics and Computer Science,Fujian Normal University,Fuzhou 350007,PR Chinaa r t i c l e i n f o Keywords:Complex symmetric linear system Iterative methods Correction technique The PIU methodConvergence analysis Numerical experimentsa b s t r a c tIn previous years,Bai and Wang presented a class of parameterized inexact Uzawa (PIU)methods for solving the generalized saddle point problems.In this paper,we consider the same method for iteratively solving the complex symmetric linear systems.Our main contribution is accelerating the convergence of the parameterized inexact Uzawa method by correction technique.First,the corrected model for the PIU method is established and the corrected PIU method is presented.Then we study the convergence property of the cor-rected PIU method.In fact,the corrected PIU method can converge faster than some Uza-wa-type and HSS-like methods.Finally,numerical experiments on a few model problems are presented to illustrate the theoretical results and examine the numerical effectiveness of the new method.Ó2015Elsevier Inc.All rights reserved.1.IntroductionLet n be a positive integer.We consider the iterative solution of systems of linear equations of the formAx ¼b ;A 2C n Ân and x ;b 2C n ;ð1:1Þwhere A 2C n Ân is a complex symmetric matrix of the formA ¼W þi T ð1:2Þand W ;T 2R n Ân are real,symmetric,and positive definite matrices.Here and in the sequel,we use i ¼ffiffiffiffiffiffiffiÀ1p to denote theimaginary unit.Complex symmetric linear systems of this kind arise in many important problems in scientific computing and engineering applications.For example,FFT-based solution of certain time-dependent PDEs [1],diffuse optical tomography [2],algebraic eigenvalue problems [3,4],molecular scattering [5],numerical solutions of the complex Helmholtz equation and numerical computations in eddy current problems.For more examples,we refer to [6–12]and the references therein.Hence,there is a strong need for the fast solution of complex symmetric linear systems.For solving the complex symmetric linear system (1.1)efficiently,van derVorst and Mellissen [9]proposed the conjugate orthogonal conjugate gradient (COCG)method,which is regarded as an extension of the Conjugate Gradient (CG)method [13].Relatively complicated but robust methods such as QMR [14],CSYM [15],and Bi-CGCR [16]are also useful.QMR is/10.1016/j.amc.2015.01.0230096-3003/Ó2015Elsevier Inc.All rights reserved.qThe Project supported by National Natural Science Foundation of China (Grant Nos.11071041,11201074),Fujian Natural Science Foundation (Grant No.2013J01006),The University Special Fund Project of Fujian (Grant No.JK2013060)and R&D of Key Instruments and Technologies for Deep Resources Prospecting (the National R&D Projects for Key Scientific Instruments)under Grant No.ZDYZ2012-1-02-04.⇑Corresponding author.E-mail address:macf88@ (C.-F.Ma).derived from the complex symmetric Lanczos method,CSYM is obtained from the idea of QMR and tridiagonalization of A by Householder reflections,and Bi-CGCR is derived from a particular case in Bi-CG [17]for solving non-Hermitian linear sys-tems.In [18]Sogabe and Zhang extended the conjugate residual (CR)method described in [19,20]to complex symmetric linear systems based on an observation of deriving CG,CR,and COCG.Moreover,Based on the Hermitian and skew-Hermi-tian splitting (HSS)A ¼H þS ;of the matrix A 2C n Ân ,withH ¼12ðA þA ÃÞand S ¼12ðA ÀA ÃÞbeing the Hermitian and skew-Hermitian parts and A Ãbeing the conjugate transpose of the matrix A 2C n Ân ,we can apply theHSS iteration method [21]or its preconditioned variant PHSS (i.e.,the preconditioned HSS,see [22])which were proposed by Bai and his co-authors to compute an approximate solution of the linear system (1.1).In addition,the convergence properties of the PHSS method can be found in [23].In [24],Bai,Golub and Ng further generalized the technique for constructing HSS iteration method for solving large sparse non-Hermitian positive definite system of linear equations to the normal/skew-Hermitian (NS)splitting obtaining a class of normal/skew-Hermitian splitting (NSS)iteration methods.Theoretical analyses shown that the NSS iteration method converges unconditionally to the exact solution of the system of linear Eq.(1.1).A potential difficulty with the HSS iteration approach is the need to solve the shifted skew-Hermitian sub-system of linear equations at each iteration step.Hence,Bai,Benzi and Chen presented a modification of the HSS iteration scheme in [25]and some of its basic properties are studied.In [26],the authors proposed a preconditioned variant of the modified HSS (PMHSS)iteration method for solving the complex symmetric systems of linear equations.Moreover,the authors in [27]pre-sented analytical and extensive numerical comparisons of some available numerical solution methods for the complex val-ued linear algebraic systems (1.1).In this paper,we study the corrected parameterized inexact Uzawa method for solving the complex symmetric linear sys-tem (1.1).We call this new method CPIU method for the sake of simplicity.We also introduce the overall reduction coeffi-cient a ði Þto measure the effectiveness of the correction process.If the overall reduction coefficient a ði Þ<1is satisfied in each iteration step,then the corrected PIU method will converge faster than the PIU method.Moreover,sufficient conditions for the convergence of the CPIU method for the linear system (1.1)are also provided in the paper.The paper is organized as follows.In Section 2,we introduce the parameterized inexact Uzawa method and establish the corrected PIU model for the Eq.(1.1),then the CPIU method is presented.In Section 3,The corrected PIU model is solved.Moreover,we analyze the convergence property of the CPIU method.In Section 4,some numerical experiments are given to show the efficiency of the CPIU method.Finally,some conclusion remarks are proposed in Section 5.The following notations will be used throughout this paper.We denote the identity matrix and the 0-matrix by I and O ,respectively.½u ;v denote the inner product of vectors u and v .For a matrix A ,we denote the spectral norm of A by k A k 2.And denote the range and the null spaces of A by R ðA Þand N ðA Þ,respectively.The conjugate transpose of A is denoted by A H .More-over,the Moore–Penrose inverse of A is defined as the unique matrix A þwhich satisfies the following four matrix equations:AA þA ¼A ;A þAA þ¼A þ;ðAA þÞH¼AA þ;ðA þA ÞH¼A þA :2.The corrected PIU methodLet x ¼y þiz and b ¼p þiq ,then from (1.1)and (1.2),we can get ðW þiT Þðy þiz Þ¼p þiq ,which implies that we can obtain the following block two-by-two systems of linear equationDX ¼W ÀT TWyz¼p q¼g :ð2:1ÞConversely,from the linear Eq.(2.1),we can get the complex symmetric linear system (1.1).So the complex symmetric linear system (1.1)is formally identical to the above block two-by-two systems of linear Eq.(2.1).Based on the PMHSS precondi-tioning matrix,Bai in [28]constructed a class of rotated block triangular preconditioners for linear system (2.1),and analyze the eigen-properties of the corresponding preconditioned matrices.Moreover,the block two-by-two systems of linear Eq.(2.1)can be formally regarded as a special case of the generalized saddle point problem [29,30].It frequently arises from finite element discretizations of elliptic partial differential equation (PDE)-constrained optimization problems such as dis-tributed control problems [31–35]and so on.Based on the parameterized inexact Uzawa (PIU)iteration method for solving the following generalized saddle point problemA B B >ÀCe y e z¼e p e q;ð2:2Þ12Q.-Q.Zheng,C.-F.Ma /Applied Mathematics and Computation 256(2015)11–19in this section we derive a corrected PIU iteration method for solving the block two-by-two systems of linear Eq.(2.1).To this end,we first introduce the PIU iteration method proposed in Bai and Wang [36]for the Eq.(2.2).This PIU iteration method is methodically described in the following.Method 2.1.Given initial vectors e y ð0Þ2R n and ez ð0Þ2R n and two relaxation factors x ;s with x ;s –0.For k ¼0;1;2;...,until the iteration sequence fðe y ðk Þ>;ez ðk Þ>Þ>g converges to the exact solution of the saddle point problem (2.2),compute e y ðk þ1Þ¼e y ðk Þþx P À1ðe p ÀA e y ðk ÞÀB e z ðk ÞÞ;e z ðk þ1Þ¼e z ðk Þþs Q À1ðB >e y ðk þ1ÞÀC ez ðk ÞÀe q Þ;(where P 2R n Ân and Q 2R n Ân are prescribed symmetric positive definite matrices.If P –A and x ¼s ¼1then the PIU iteration Method 2.1yields inexact Uzawa method [37–39]for solving the saddle point problems.Here,in this paper we consider the case that P ¼Q ¼A and x ¼1.Based on the above Method 2.1,we can get the following iteration method.Method 2.2.Given initial vectors y ð0Þ2R n ;z ð0Þ2R n and the relaxation factor s with s –0.For k ¼0;1;2;...,until the iter-ation sequence fðy ðk Þ>;z ðk Þ>Þ>g converges to the exact solution of the saddle point problem (2.1),computey ðk þ1Þ¼W À1ðp þTz ðk ÞÞ;z ðk þ1Þ¼z ðk ÞÀs W À1ðT >y ðk þ1ÞþWz ðk Þþq Þ:(Denote X ðk Þ¼y ðk Þz ðk Þ,then the PIU Method 2.2can be rewritten as X ðk þ1Þ¼HX ðk ÞþM À1g ;ð2:3ÞwhereH ¼OW À1TO I Às W À1ðC þTW À1T Þ!and M ¼W O T1s W!:Let N ¼M ÀD ¼O T O1s À1ÀÁW,thenD ¼M ÀN ð2:4Þdefines a splitting of the coefficient matrix D of the block two-by-two systems of linear Eq.(2.1),and the PIU Method 2.2can also be induced by the matrix splitting (2.4).Easily,we see that H ¼M À1N is the iteration matrix of the PIU Method 2.2.Therefore,the PIU Method 2.2is convergent if and only if the spectral radius of the iteration matrix H is less than one,i.e.,q ðH Þ<1.See [40–42].Assume X 1;X 2;...;X m are approximate solutions of Eq.(2.1),where m >1is a positive integer and DX i –g ði ¼1;2;...;m Þ.The basic idea of the correction model for the parameterized inexact Uzawa Method 2.2is:construct a vector b Xwhich is a better approximate solution of Eq.(2.1)than X i ði ¼1;2;...;m Þ.That is,k g ÀD b X k 2<min 16i 6mk g ÀDX i k 2:ð2:5ÞThe vector b Xthat satisfies (2.5)is called the corrected solution of X 1;X 2;...;X m ,and the progress to determine the vector X is called CPIU process.The corrected model for the PIU is presented as follows.The CPIU model.b X¼a 1X 1þa 2X 2þÁÁÁþa c X m ;a 1þa 2þÁÁÁþa m ¼1;k g ÀD b Xk 2¼min :ð2:6ÞWith the above CPIU model we can obtain the following detailed methodic description of the CPIU method for the blocktwo-by-two system of linear Eq.(2.1).The CPIU method.Give an initial vector X ð0Þ2R 2n ,let k ¼0. WhileX ðk Þ1¼X ðk À1ÞFor ðj ¼2;3;...;m ÞX ðk Þj¼HX ðk Þj À1þM À1gend ForX ðk Þ¼a ðk Þ1X ðk Þ1þa ðk Þ2X ðk Þ2þÁÁÁþa ðk Þm X ðk Þm k ¼k þ1End WhileQ.-Q.Zheng,C.-F.Ma /Applied Mathematics and Computation 256(2015)11–1913Here aðkÞiði¼1;2;...;mÞis determined by the k th CPIU progress.How to solve the above CPIU model will be illustrated detailedly in Section3.From the Eq.(2.5),we can see that if the iteration(2.3)is convergent,i.e.,qðHÞ<1,then the sequence f XðkÞg will converge.3.Analysis for the CPIU methodIn this section,we will solve the CPIU model(2.6).Moreover,the sufficient conditions for the convergence of the CPIU method are also provided.When the linear system EX¼r is consistent,then the vector X02C2n that satisfies k X0k2¼min EX¼r k X k2is called the min-imum norm solution of linear system EX¼r.When the linear system EX¼r is not consistent,then the vector X02C2n thatsatisfies k X0k2¼min min k EXÀr k2k X k2is called the minimal norm least squares solution of linear system EX¼r.Lemma3.1.If the linear system EX¼r is consistent,then the minimum norm solution X0of EX¼r is unique and X0¼Eþr.Proof.Firstly,we prove X02RðE HÞholds true.If X0R RðE HÞ,then we can obtain the decomposition X0¼X1þX2;X12RðE HÞ;X22NðEÞ¼½RðE HÞ ?;X2–0:Note X1?X2,we obtaink X0k22¼k X1k22þk X2k22>k X1k22:Because EX1¼EX1þEX2¼EX0¼r,then we get X0is not the minimum norm solution of EX¼r.This is a contradiction,so X02RðE HÞ.Then we illustrate the uniqueness of minimum norm solution X0.If Z0is another minimum norm solution of EX¼r,then Z02RðE HÞ.So we haveX0ÀZ02RðE HÞ¼½NðDÞ ?:Moreover,because EðX0ÀZ0Þ¼0,then we have X0ÀZ02NðEÞ.So X0ÀZ0¼0,which implies X0¼Z0.Now we illustrate X0¼Eþr.Obviously,the general solution of the equation EX¼r isX¼EþrþðIÀEþEÞy;ðy2C nÞ:Take y¼0,then we obtain that X0¼Eþr is a solution of EX¼r.Because½ðIÀEþEÞy;X0 ¼y HðIÀEþEÞH X0¼y HðIÀEþEÞEþr¼y HðEþÀEþEEþÞr¼0;then we havek X k22¼k X0k22þkðIÀEþEÞy k22P k X0k22;which implies that X0¼Eþr is the minimum norm solution of linear system EX¼r.The proof is completed.hLemma3.2.If the linear system EX¼r is not consistent,then the minimal norm least squares solution X0of EX¼r is unique and X0¼Eþr.Proof.From Lemma3.1,we can demonstrate that the result holds true.hTogether with Lemmas3.1and3.2we prove the following result which show how to solve the CPIU model(2.6).Theorem3.1.If X1;X2;...;X m are approximate solutions of Eq.(2.1),for k¼1;2;3;...,let r k¼gÀDX k.Denoteb i¼r iÀr1ði¼2;3;...;mÞ;Q¼ðb2;b3;...;b mÞ;ð3:1Þy 1¼ða2;a3;...;a mÞ>¼ÀQþr1;a1¼1ÀX mi¼2a i:ð3:2ÞThen b X¼a1X1þa2X2þÁÁÁþa m X m is the corrected solution of X1;X2;...;X m.Proof.Assume b X¼b1X1þb2X2þÁÁÁþb m X m is the corrected solution of X1;X2;...;X m of CPIU model(2.6),where b1þb2þÁÁÁþb m¼1.Then we haveb X¼X1þX mi¼2b iðX iÀX1Þ;14Q.-Q.Zheng,C.-F.Ma/Applied Mathematics and Computation256(2015)11–19which impliesrðb XÞ¼gÀD b X¼r1þX mi¼2b iðr iÀr1Þ:From(3.1)and(3.2),we obtain that the above equation can be rewritten as rðb XÞ¼r1þQðb2;b3;...;b mÞ>.This implies thatfinding a vector b X which can minimize k rðb XÞk2is equivalent tofinding y2C mÀ1which can minimize k r1þQ y k2.By making use of Lemmas 3.1and 3.2,the minimum norm solution satisfies k y1k2¼min k r1þQ y k2k y k2isy1¼ða2;a3;...;a mÞ>¼Qþr1.Let a1¼1ÀP mi¼2a i,then we can obtain thatb X¼a1X1þa2X2þÁÁÁþa m X m is the correctedsolution of X1;X2;...;X m of the CPIU model(2.6).The proof is completed.hRemark3.1.Take j2f1;2;...;m g,and denotebi¼r iÀr jði–jÞ;Q j¼ðb1;...;b jÀ1;b jþ1;...;b mÞ;y j ¼ða1;...;a jÀ1;a jþ1;...;a mÞ>¼ÀQþjr j;a1¼1ÀX mi¼2a i:Then from the proof of Theorem 3.1,we can obtain that b X¼a1X1þa2X2þÁÁÁþa m X m is the corrected solution of X1;X2;...;X m of the CPIU model(2.6).Lemma3.3.Let A2R nÂn;P2R nÂn and Q2R nÂn be symmetric positive definite,and B2R nÂn is nonsingular.When C¼d Q with d>0a real constant,then the PIU Method2.1for the generalized saddle point problem(2.2)is convergent,provided that x sat-isfies0<x<2gmaxand s satisfies0<s<2dand s½2dþðl maxÀd g maxÞx <2ð2Àg max xÞ;ð3:3Þwhere g max is the largest eigenvalue of the matrix PÀ1A,and l max is the largest eigenvalue of the matrix J¼QÀ1B>PÀ1B.Proof.See Theorem2.2of[36].hTheorem3.2.If W2R nÂn and T2R nÂn are symmetric positive definite,then the PIU Method2.2for solving the block two-by-two systems of linear Eq.(2.1)is convergent,provided that s satisfies0<s<2and sð1þl maxÞ<2;ð3:4Þwhere l max is the largest eigenvalue of the matrix WÀ1T>WÀ1T.Proof.Let A¼C¼W;B¼ÀT;P¼Q¼A and x¼1,then we can obtain d¼1;PÀ1A¼I and J¼QÀ1B>PÀ1B¼WÀ1T>WÀ1T. Also from PÀ1A¼I we have g max¼1,then from Lemma3.3,we can get the result directly.hIn order to measure the effect of the CPIU correction process,we introduce the overall reduction coefficient a¼k rðXÞk2min16i6m k r i k2:ð3:5ÞFrom(2.5)we have06a<1.If a¼1,then the CPIU process(2.6)is failed.If06a<1,then the CPIU process is successful. And the smaller a is,the better of the effect of the CPIU progress is.In particular,if a¼0,then the CPIU process isfinished, i.e.,the corrected solution is the exact solution of Eq.(2.1).Denote the overall reduction coefficient of the k th correction progress as aðkÞ,thenaðkÞ¼k rðX ðkÞÞk2min16j6m k rðkÞj k 2;where rðkÞj ¼gÀDXðkÞjðj¼1;2;...;mÞand rðXðkÞÞ¼gÀDXðkÞ,with XðkÞiði¼1;2;...;mÞis the approximate solutions of k thCPIU progress and XðkÞis the corrected solution of XðkÞ1;XðkÞ2;...;XðkÞm.With the overall reduction coefficient of the k th correction progress aðkÞ,we can get the following result.Theorem3.3.If W2R nÂn and T2R nÂn are symmetric positive definite,then the CPIU method for solving the block two-by-two systems of linear Eq.(2.1)is convergent,provided that s satisfiesQ.-Q.Zheng,C.-F.Ma/Applied Mathematics and Computation256(2015)11–19150<s<2and sð1þl maxÞ<2;ð3:6Þwhere l max is the largest eigenvalue of the matrix WÀ1T>WÀ1T.Moreover,if aðkÞ<1is satisfied,then the CPIU method will con-verge faster than the PIU method.Proof.If0<s<2and sð1þl maxÞ<2,then from Theorem3.2and the CPIU method which we have presented in Section2 for solving the block two-by-two system of linear Eq.(2.1)is convergent.Moreover,if aðkÞ<1is satisfied,then the CPIU method will converge faster than the PIU method.This complete the proof.h4.Numerical experimentsIn this section,we perform some numerical examples to illustrate the theoretical results and show the effectiveness of the CPIU iteration method for solving the complex symmetric linear system(1.1)in terms of both iteration steps(denoted as IT) and computing time(in s,denoted as CPU),and the norm of the residual(denoted as‘‘RES’’)defined byRES¼ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffik pÀWyðkÞþTzðkÞk22þk qþTyðkÞþWzðkÞk22 q:In actual computations,the iteration schemes are started from the zero vector and terminated if the current iterations satisfy ERR610À6or the number of the prescribed iteration steps k¼500are exceeded,whereERR¼ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffik pÀWyðkÞþTzðkÞk22þk qþTyðkÞþWzðkÞk22 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffik pÀWyð0ÞþTzð0Þk22þk qþTyð0ÞþWzð0Þk22 q:All experiments are performed in MATLAB(version7.4.0.336(R2010b))with machine precision10À16,and all experiments are implemented on a personal computer with2.20GHz central processing unit(Intel(R)Core(TM)i3–2310M),2.00G mem-ory and Win7operating system.The CPIU iteration method is compared with the preconditioned Uzawa(PU)method[39]and PMHSS[26]methods.For the tests reported in this section we used the optimal values of the parameter s(denoted by s opt)for the CPIU,PU and PMHSS iteration methods.The experimentally found optimal parameters s opt are the ones resulting in the least numbers of iterations for the three methods for each of the numerical examples and for each choice of the spatial mesh-sizes.Example4.1(See[10,25,26]).The system of linear Eq.(1.1)is of the formKþ3Àffiffiffi3ps I!þi Kþ3þffiffiffi3ps I!"#x¼b;ð4:1Þwhere s is the time step-size and K is thefive-point centered difference matrix approximating the negative Laplacian oper-ator L¼ÀD with homogeneous Dirichlet boundary conditions,on a uniform mesh in the unit square½0;1 ½0;1 with the mesh-size h¼1.The matrix K2R nÂn possesses the tensor-product form K¼I V lþV l I,with V l¼hÀ2tridiag ðÀ1;2;À1Þ2R lÂl.Hence,K is an nÂn block-tridiagonal matrix,with n¼l2.We takeW¼Kþ3Àffiffiffi3ps I;and T¼Kþ3þffiffiffi3ps Iand the right-hand side vector b with its j th entry½bjbeing given by½bj ¼ð1ÀiÞjsðjþ1Þ2;j¼1;2;...;n:In our tests,we take s¼h.Furthermore,we normalize coefficient matrix and right-hand side by multiplying both by h2.For more details,we refer to[10].With respect to different sizes of the coefficient matrix,we list IT,CPU and RES about the PU,PMHSS and CPIU(m¼8) methods in Table1for Example4.1.By comparing the results in Table1,we observe that the CPIU iteration method outper-forms the PU and PMHSS methods,as it requires much less time and iteration steps to achieve the stopping criterion.Moreover,we also take m¼2;3;4;5;6;7in this numerical experiment to illustrate the special property of the CPIU method.The results for different parameters m of the CPIU iteration method for Example4.1are presented in Table2.From the results of Table2we observe that the iteration steps of the CPIU iteration method become smaller when the parameter m becomes larger which means that the CPIU iteration method converges faster if the parameter m becomes larger.Moreover, the solution that our CPIU iteration method found is more accurate than the other two methods,because the norm of the residual(RES)is more smaller than the PU and PMHSS methods.So we can get a good approximate solutions by taking a proper parameter m in our actual calculation.That illustrates that CPIU iteration method is a very good method for solving the complex symmetric linear systems.16Q.-Q.Zheng,C.-F.Ma/Applied Mathematics and Computation256(2015)11–19Example 4.2(See [11,25,26]).The system of linear Eqs.(1.1)is of the form½ðÀx 2M þK Þþi ðx C V þC H Þ x ¼b ;ð4:2Þ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 x is the driving circular frequency.We take C H ¼l K with l a damping coefficient,M ¼I ;C V ¼10I ,and K the five-point centered difference matrix approximating the negative Laplacian operator with homogeneous Dirichletboundary conditions,on a uniform mesh in the unit square [0,1]Â[0,1]with the mesh-size h ¼1.The matrix K 2R n Ân pos-sesses the tensor-product form K ¼I B l þB l I ,with B l ¼h À2tridiag ðÀ1;2;À1Þ2R l Âl .Hence,K is an n Ân block-tridiagonalmatrix,with n ¼l 2.In addition,we set x ¼p ;l ¼0:02,and the right-hand side vector b to be b ¼ð1þi ÞA 1,with 1being thevector of all entries equal to 1.As before,we normalize the system by multiplying both sides through by h 2.In fact,this com-plex symmetric system of linear equation arises in direct frequency domain analysis of an n -degree-of-freedom (n -DOF)lin-ear system.For more details,we refer to [10,40].With respect to different sizes of the coefficient matrix,we list IT,CPU and RES about the PU,PMHSS and CPIU (m ¼8)methods in Table 3for Example 4.2.By comparing the results in Tables 2,we observe that the CPIU iteration method out-performs the PU and PMHSS methods,as it requires much less time and iteration steps to achieve the stopping criterion.Moreover,we also take m ¼2;3;4;5;6;7in this numerical experiment.The results for different parameters m of the CPIU iteration method for Example 4.2are presented in Table 4.As we can see from the results of Table 4,we can always get a good approximate solutions by taking a proper parameter m in our actual calculation.Example 4.3(See [25,26]).The system of linear Eq.(1.1)is of the form ðW þi T Þx ¼b ,withT ¼I V þV Iand W ¼10ðI V c þV c I Þþ9ðe 1e >l þe l e >1Þ I ;where V ¼tridiag ðÀ1;2;À1Þ2R l Âl ;V c ¼V Àe 1e >l Àe l e >12R l Âl ;e 1and e l are the first and the last unit vectors in R l ,respec-tively.We take the right-hand side vector b to be b ¼ð1þi ÞA 1,with 1being the vector of all entries equal to 1.Here T and W correspond to the five-point centered difference matrices approximat-ing the negative Laplacian operator with homogeneous Dirichlet boundary conditions and periodic boundary conditions,respectively,on a uniform mesh in theunit square ½0;1 ½0;1 with the mesh-size h ¼1.This problem is an artificially constructed one,but it is quite challenging for iterative solvers.Table 1IT,CPU and RES for PU,PMHSS and CPIU methods for Example 4.1.Method l8162432PUs opt0.1960.2160.2260.210IT 252384164CPU 0.0080.085 1.1677.485RES5.09e À8 3.94e-8 3.15e À8 2.30e À8PMHSSs opt1.091 1.356 1.350 1.120IT 21212121CPU 0.0090.104 1.245 3.099RES4.33e À8 2.02e À8 3.04e À8 3.20e À7CPIUs opt0.1960.2160.2260.210IT 2222CPU 0.0070.0790.056 2.771RES1.53e À133.29e À113.57e À101.29e À9Table 2Numerical results for different parameters m of CPIU method for Example 4.1.l IT CPU RES IT CPU RES IT CPU RES m ¼2m ¼3m ¼48190.0947.80e À870.0598.18e À840.042 1.50e À816230.126 3.84e À890.095 1.71e À850.087 3.85e À924240.875 2.94e À890.633 1.83e À860.606 1.92e À93231 4.007 1.82e À810 3.103 1.12e À86 2.9867.79e À9m ¼5m ¼6m ¼7830.030 1.63e À920.007 3.44e À920.021 1.97e À111630.079 2.68e À930.094 1.46e À1020.089 1.93e À92440.655 1.71e À930.6519.59e À1020.561 1.74e À93242.9526.34e À932.8693.36e À933.0343.28e À11Q.-Q.Zheng,C.-F.Ma /Applied Mathematics and Computation 256(2015)11–1917。