- Email: [email protected]
Gauss type preconditioning techniques for linear systems
Applied Mathematics and Computation 188 (2007) 612–633 www.elsevier.com/locate/amc
Gauss type preconditioning techniques for linear systems Yong Zhang a, Ting-Zhu Huang
a,*
, Xing-Ping Liu
q
b
a
b
School of Applied Mathematics, University of Electronic Science and Technology of China, Chengdu, Sichuan 610054, PR China National Key Lab of Computational Physics, Institute of Applied Physics and Computational Mathematics, Beijing 100088, PR China
Abstract Many researchers have considered preconditioners chosen to eliminate the off-diagonal elements of the coefficient matrix of a linear system. In this work, we generalize the left Gauss type preconditioners [Y. Zhang, T.Z. Huang, X.P. Liu, Modified iterative methods for nonnegative matrices and M-matrices linear systems, Comput. Math. Appl. 50 (2005) 1587–1602] which eliminate the strictly lower triangular elements. Right Gauss type preconditioners that eliminate strictly upper triangular elements are proposed in this paper. These Gauss type preconditioners are partly derived from the LU factorization method. Theoretic analysis on spectral radii of the two kinds of Gauss type preconditioners is given. Numerical experiments are used to show the performance of the improved inbuilt left and right Gauss type preconditioning algorithms associated with Jacobi type and Gauss–Seidel type iterative methods. 2006 Elsevier Inc. All rights reserved. Keywords: Preconditioning; Jacobi and Gauss–Seidel type iterative methods; LU factorization; Gauss elimination; M-Matrix
1. Introduction Consider the numerical solution of the linear system Ax ¼ b; where A = (aij) 2 R
ð1Þ n·n
is a known nonsingular M-matrix.
Definition 1 [6]. A matrix A is a Z-matrix if aij 6 0, i, j = 1, 2, . . . , n; i 5 j. Definition 2 [5]. A matrix A is an L-matrix if aii > 0, i = 1, 2, . . . , n and aij 6 0 for all i, j = 1, 2, . . . , n; i 5 j. Definition 3 [6]. A matrix A is an M-matrix if A has the form A = sI B, s > 0, B P 0 and s > q(B) (the spectral radius of B). q
Supported by NCET of China (NCET-04-0893), foundation of the National Key Lab of Computational Physics, Beijing Appl. Phys. Comput. Math. Institute and applied basic research foundations of Sichuan province (05JY029-068-2). * Corresponding author. E-mail addresses: [email protected] (Y. Zhang), [email protected] (T.-Z. Huang). 0096-3003/$ - see front matter 2006 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2006.10.053
Y. Zhang et al. / Applied Mathematics and Computation 188 (2007) 612–633
613
Since A is an M-matrix, aii > 0, i = 1, 2, . . . , n, which shows that A is an L-matrix [1]. Without loss of generality, we can write A ¼ I L U;
ð2Þ
where I is the identity matrix, L and U are strictly lower triangular and strictly upper triangular parts of A respectively. M = I, N = L + U and M = I L, N = U leads to the classical Jacobi and classical Gauss–Seidel iterative methods respectively. Then the corresponding iteration matrices of the classical Jacobi and Gauss– Seidel methods are given by B = L + U and T = (I L)1U respectively. One can transform the original system (1) into the preconditioned form PAx ¼ Pb
ð3Þ
or AQy ¼ b;
x ¼ Qy:
ð4Þ
Here P in (3) and Q in (4) are called the left preconditioner and the right preconditioner, respectively. Then one can define two kinds of basic iterative schemes respectively ðkÞ xðkþ1Þ ¼ M 1 þ M 1 p N px p Pb;
k ¼ 0; 1; . . . ;
ð5Þ
where PA = Mp Np and Mp is nonsingular, and ðkÞ y ðkþ1Þ ¼ M 1 þ M 1 p N py p b;
k ¼ 0; 1; . . . ;
ð6Þ
where AQ = Mp Np and Mp is nonsingular. Many researchers have considered the left preconditioners applied to system (1) that make the associated Jacobi and Gauss–Seidel type methods converge faster than the original ones. In a simpler form, Milaszewicz [2] considered the preconditioner 1 0 1 C B C B a21 1 C B C B a31 1 C B C B . . C B . . ð7Þ . P1 ¼ B . C; C B C B ai1 1 C B C B . .. C B . . A @ . an1 1 which eliminates the elements of the first column of A below the diagonal. Hadjidimos et al. [1], for M-matrices systems, generalized Milaszewicz’s preconditioning technique and considered the preconditioner 1 0 1 C B C B a2 a21 1 C B C B a3 a31 1 C B C B . . C B . . P 1 ðaÞ ¼ B . ð8Þ . C: C B C B ai ai1 1 C B C B . . C B . . . A @ . an an1
1
Kohno et al. [7] parameterized the preconditioner presented by Gunawardena et al. [8] and considered a preconditioner
614
Y. Zhang et al. / Applied Mathematics and Computation 188 (2007) 612–633
0 B B B B B B B P 2 ðaÞ ¼ B B B B B B @
1
a1 a12 1
1 a2 a23 .. .
..
.
1
ai ai;iþ1 .. .
C C C C C C C C: C C .. C . C C 1 an1 an1;n A 1
ð9Þ
Recently, Zhang et al. [3] proposed the left Gauss type preconditioning techniques based on Hadjidimos et al. [1], Milaszewicz [2] and LU factorization method [4]. The Gauss type preconditioning technique utilizes the Gauss transformation matrices [4] as the base of the Gauss type preconditioner. The construction of Gauss transformation matrices [4] is as follows: 1 0 1 0 0 0 B .. .. C .. .. B. .C . . C B C B 1 0 0C B0 C; k ¼ 1; 2; . . . ; n 1; Mk ¼ B ð10Þ B0 skþ1 1 0C C B B. .. C .. .. C B. @. .A . . 0
sn
0
1
where sT ¼ ð0; . . . ; 0; skþ1 ; . . . ; sn Þ;
si ¼
aik ; i ¼ k þ 1; . . . ; n: akk
The fore k elements of the vector s 2 Rn are zero, and s is called a Gauss vector. In some cases the two preconditioners (7) and (8) have the same expression with one of Gauss transformation matrices, denoted by M1, of LU factorization methods [4], that is to say M 1 ¼ P 1 ¼ P 1 ðaÞ;
T
T
a ¼ ða1 ; a2 ; . . . ; an1 Þ ¼ ð1; 1; . . . ; 1Þ :
Due to the similarity between the preconditioner (7) and (8) and Gauss transformation matrix (10) (i.e. M1), Zhang et al. [3] consider the following left preconditioners: P 1 ¼ M 1; P 2 ¼ M 2M 1; .. . P m ¼ M m M m1 M 2 M 1 :
ð11Þ
Then, • preconditioners Pm (m = 1, 2, . . . , n 1) are all lower triangular matrices with unity diagonal entries; • the entries of the fore m columns of PmA (m = 1, 2, . . . , n 1) are all zero below the diagonal; • Pn1A is an upper triangular matrix in particular. This work is organized as follows. In Section 2, we extend the left Gauss preconditioners by giving a family of preconditioners based on the elimination of elements at the strictly upper triangular part (i.e. up the diagonal), and, theoretic analysis is presented which theoretically confirm the results in [3]. In Section 3, improved Gauss preconditioning algorithms associated with inbuilt technique are given. In Section 4, numerical experiments are presented to show the performance and efficiency of improved Gauss preconditioning algorithms.
Y. Zhang et al. / Applied Mathematics and Computation 188 (2007) 612–633
615
2. Generalized Gauss type preconditioning techniques 2.1. Right Gauss type preconditioners From the fore section, we know that the left Gauss type preconditioner (11) could increase the number of zero elements in the strictly lower triangular part of the preconditioned matrix (PmA) column by column. Similarly, we can construct the preconditioners which eliminate the elements up the diagonal row by row if possible. First, let Qm be such kind of preconditioners with order m. The following expression shows the way to construct Qm: ðAQm ÞT ¼ QTm AT ;
ð12Þ QTm ,
then, in order to find the kind of wanted preconditioners, we only need to construct the matrix just by the fore mentioned left preconditioning techniques which could make the entries of the fore m columns of the strictly lower triangular part of the preconditioned matrix QTm AT all zero. We note that these preconditioning matrices Qm could make the entries of the fore m rows of the strictly upper triangular part of the preconditioned matrix AQm all zero. For brevity, these preconditioners Qm is named as right Gauss preconditioners which could be applied to (1) as right preconditioners. Similarly, we define the corresponding right Gauss transformation matrices 0
1 B. B .. B B B0 Mk ¼ B B0 B B. B. @.
0 .. . 1 0 .. .
1 .. .
0
0
0
0 .. . ukþ1
1 0 .. C . C C C un C C; 0 C C .. C C . A 1
k ¼ 1; 2; . . . ; n 1;
ð13Þ
where uT ¼ ð0; . . . ; 0; ukþ1 ; . . . ; un Þ;
uj ¼
akj ; j ¼ k þ 1; . . . ; n: akk
The fore k elements of u 2 Rn are all 0, and u is also called a Gauss vector. We note that Qm ¼ M 1 M 2 M m1 M m ;
1 6 m 6 n 1:
Now we begin to analyze the elements of the preconditioned matrices ALm and ARm which satisfy ALm ¼ P m A
and ARm ¼ AQm :
First, let us get the expression of the elements of MkA and AM k , which help us to build the corresponding preconditioning algorithms presented in the next section, from the previous notations and considerations. Defining B ¼ ðbij Þ ¼ M k A;
j ¼ 1; 2; . . . ; n;
we have bij ¼ i.e.
( bij ¼
aij ;
i ¼ 1; 2; . . . ; k;
aij si akj ;
i ¼ k þ 1; . . . ; n;
aij ; i ¼ 1; 2; . . . ; k; aik akj aij akk ; i ¼ k þ 1; . . . ; n:
ð14Þ
616
Y. Zhang et al. / Applied Mathematics and Computation 188 (2007) 612–633
Defining B ¼ ð bij Þ ¼ AM k ; we have bij ¼ i.e.
i ¼ 1; 2; . . . ; n; j ¼ 1; 2; . . . ; k;
aij ;
aij uj aik ; j ¼ k þ 1; . . . ; n; (
bij ¼
aij ; j ¼ 1; 2; . . . ; k; aik akj aij akk ; j ¼ k þ 1; . . . ; n:
ð15Þ
Define ðlmÞ
ðrmÞ
ALm ¼ ðaij Þ;
ARm ¼ ðaij Þ;
m ¼ 0; 1; 2; . . . ; n 1;
where ðl0Þ
ðr0Þ
A ¼ AL0 ¼ AR0 ¼ ðaij Þ ¼ ðaij Þ: By the fore known properties of Gauss preconditioners ALm and ARm can be written as 0 ðl;mÞ ðl;mÞ ðl;mÞ ðl;mÞ a12 a1m a1;mþ1 a11 B ðl;mÞ ðl;mÞ ðl;mÞ B 0 a22 a2m a2;mþ1 B B .. .. .. .. B .. B . . . . . B B ðl;mÞ ðl;mÞ 0 0 a a B m;mþ1 mm ALm ¼ B B ðl;mÞ 0 amþ1;mþ1 B 0 B B . .. .. .. .. B .. . . . . B B . .. .. .. .. B . . . . . @ . ¼
0 Lm1
Lm2
Lm3
Lm4
ðl;mÞ
0
an;mþ1
and considerations, the preconditioned matrices ðl;mÞ
a1n
1
C C C C C C C C aðl;mÞ C mn C ðl;mÞ C amþ1;n C C .. C . C C .. C C . A aðl;mÞ nn ðl;mÞ
a2n .. .
ð16Þ
ð17Þ
and 0
ARm
ðr;mÞ
a11
B B ðr;mÞ B a21 B B . B .. B B ðr;mÞ B am1 ¼B B ðr;mÞ B amþ1;1 B B . B . B . B . B . @ . ðr;mÞ
¼
an1 Rm1 Rm3
0 ðr;mÞ
a22 .. .
ðr;mÞ am2 ðr;mÞ amþ1;2
.. . .. .
.. . .. . ðr;mÞ
an2 Rm2 ; Rm4
0 .. .
.. . .. .
.. . .. .
0
0 .. . .. .
aðr;mÞ mm ðr;mÞ amþ1;m
0
ðr;mÞ amþ1;mþ1
.. . .. .
.. . .. .
aðr;mÞ nm
an;mþ1
ðr;mÞ
0 .. . .. .
1
C C C C C C C C 0 C C ðr;mÞ C amþ1;n C C .. C C . C .. C C . A
ð18Þ
aðr;mÞ nn ð19Þ
Y. Zhang et al. / Applied Mathematics and Computation 188 (2007) 612–633
617
where the block form in (17) and (19) is the same, and square matrices Lm1 and Rm1 are all of dimension m · m and square matrices Lm4 and Rm4 are all of dimension (n m) · (n m). Now, we can get the analysis expressions of two kinds of preconditioned matrices in terms of (14) and (15): • For m = 1, we have ðARm Þ11 ¼ ðALm Þ11 ¼ a11 ; ai1 a1j ; a11
ðALm Þij ¼ ðARm Þij ¼ aij
i ¼ m þ 1; . . . ; n; j ¼ m þ 1; . . . ; n:
ð20Þ
• For m = 2, we have ðl;mÞ
ðALm Þ22 ¼ a22
ðl;m1Þ
¼ a22
;
ðr;mÞ
ðr;m1Þ
ðARm Þ22 ¼ a22
¼ a22
ðr;m1Þ
¼ a22 ;
:
By (20) we know ðl;m1Þ
a22
ðl;21Þ
¼ a22
ðl;1Þ
¼ a22 ;
a22
ðr;21Þ
¼ a22
ðr;1Þ
ðl;1Þ
ðr;1Þ
a22 ¼ a22 ;
then (ALm)22 = (ARm)22. • In a general way, by (14), (15), (20), for any 1 6 m 6 n 1, we have ðALm Þij ¼ ðARm Þij ¼
ðl;m1Þ aij
ðl;m1Þ ðl;m1Þ amj ðl;m1Þ amm
aim
¼
ðr;m1Þ aij
ðr;m1Þ ðr;m1Þ amj ðr;m1Þ amm
aim
;
ð21Þ
i = m + 1, . . . , n; j = m + 1, . . . , n. Naturally, from the analysis expression (21), we conclude that (1) Preconditioned matrices ALm and ARm have the same diagonal elements, i.e. ðl;mÞ
ðALm Þii ¼ ðARm Þii ¼ aii
ðr;mÞ
¼ aii
;
1 6 i 6 n:
ð22Þ
(2) The sub matrices of the (m + 1)th row to the nth row, the (m + 1)th column to the nth column of ALm and ARm are equal, i.e. Lm4 ¼ Rm4 :
ð23Þ
That is to say ðl;mÞ
ðALm Þij ¼ ðARm Þij ¼ aij
ðr;mÞ
¼ aij
;
i; j ¼ m þ 1; . . . ; n; 1 6 m 6 n 1:
ð24Þ
(3) In particular, ðALm Þij ¼ 0
ðj ¼ 1; 2; . . . ; m; i ¼ j þ 1; j þ 2; . . . ; nÞ;
ðARm Þij ¼ 0 ði ¼ 1; 2; . . . ; m; j ¼ i þ 1; i þ 2; . . . ; nÞ: We consider the usual triangular splitting of ALm and ARm ALm ¼ DLm ELm F Lm ; where DLm, ELm, FLm are diagonal, strictly lower triangular and strictly upper triangular parts of ALm respectively, and ARm ¼ DRm ERm F Rm ; where DRm, ERm, FRm are diagonal, strictly lower triangular and strictly upper triangular parts of ARm respectively.
618
Y. Zhang et al. / Applied Mathematics and Computation 188 (2007) 612–633
Below we define the Jacobi type iteration associated with the above splittings B ¼ L þ U; BLm ¼ D1 Lm ðE Lm þ F Lm Þ; BRm ¼
D1 Rm ðE Rm
ð25Þ
þ F Rm Þ;
as well as the splittings that define the Gauss–Seidel type iteration matrices T ¼ ðI LÞ1 U ; T Lm ¼ ðDLm ELm Þ1 F Lm ; 1
T Rm ¼ ðDRm ERm Þ F Rm : In terms of the expressions of ALm (16) and ARm (18), we can express these splitting matrices in the following block form: ! ! ! 0 DL1 0 0 F L1 F L2 DLm ¼ ; ELm ¼ ; F Lm ¼ ; 0 DL4 0 EL4 0 F L4 ! ! ! 0 0 DR1 ER1 0 0 DRm ¼ ; ERm ¼ ; F Rm ¼ : 0 DR4 ER3 ER4 0 F R4 Thus, we have BLm ¼ D1 Lm ðELm þ F Lm Þ !1 DL1 0 F L1 ¼ 0 DL4 0 ¼
EL4 þ F L4 !
D1 L1 F L1
D1 L1 F L2
0
D1 L4 ðE L4 þ F L4 Þ
BRm ¼ D1 Rm ðE Rm þ F Rm Þ !1 DR1 0 ER1 ¼ 0 DR4 ER3 ¼
!
F L2
ð26Þ
;
!
0 ER4 þ F R4 !
D1 R1 ER1
0
D1 R4 ER3
D1 R4 ðE R4 þ F R4 Þ
ð27Þ
;
T Lm ¼ ðDLm ELm Þ1 F Lm ¼
¼
0
0
DL4 EL4
D1 L1
F L1
F L2
0
F L4
!
0 ðDL4 EL4 Þ
0 ¼
!1
DL1
1
F L1
F L2
0
F L4
D1 L1 F L1
D1 L1 F L2
0
ðDL4 EL4 Þ F L4
1
! !
! ;
ð28Þ
Y. Zhang et al. / Applied Mathematics and Computation 188 (2007) 612–633
619
and 1
T Rm ¼ ðDRm ERm Þ F Rm 1 DR1 ER1 0 0 ¼ DR4 ER4 ER3 0 ðDR1 ER1 Þ
¼ ¼
F R4
1
1
0
ðDR4 ER4 Þ ER3 ðDR1 ER1 Þ 0 0 : 1 0 ðDR4 ER4 Þ F R4
!
0 1
ðDR4 ER4 Þ
1
0
0
0 F R4 ð29Þ
Then we have Lemma 1. If A is an M-matrix, then for any 1 6 m 6 n 1 (1) ALm exists; (2) Mm P 0 and Pm P 0; (3) ALm and Lm4 both are M-matrices. Proof. For the case m = 1: Since A is an M-matrix, a11 > 0. Then there exist M1 and P1 such that M 1 P 0;
P 1 P 0:
Hence, by AL1 = P1A, AL1 exists. And for a Z-matrix A, ‘‘A is an M-matrix’’ is equivalent to ‘‘there exists a positive vector y (>0) 2 Rn such that Ay > 0’’ (see [6, p. 136]. P1 P 0, implies AL1 y ¼ P 1 Ay > 0: Consequently, AL1, which is a Z-matrix, is an M-matrix. Since for a Z-matrix the statement ‘‘A is an M-matrix’’ is equivalent to ‘‘all the principal minors of A are positive’’ [6, p. 134], the submatrix L14 is also an M-matrix. Similarly, we could use the same way for the case m = 1 step by step to prove the results for 2 6 m 6 n 1. The proof is completed. h Lemma 2. If A is a nonsingular M-matrix, then, for any 1 6 m 6 n 1, (1) ARm exists; (2) M m P 0 and Qm P 0; (3) ARm and Rm4 both are nonsingular M-matrices. Proof. Since A is an M-matrix, AT is also an M-matrix. For the case m = 1: For A being an M-matrix, so a11 > 0. Then there exist M 1 and Q1 such that M 1 P 0;
Q1 P 0:
Hence, for AR1 = A Q1, AR1 exists. And for a Z-matrix A, ‘‘A is an M-matrix’’ is equivalent to ‘‘there exists a positive vector y(>0) 2 Rn such that Ay > 0’’ (see [6, p. 136]). So there exists a positive vector y(>0) 2 Rn such that ATy > 0. Q1 P 0, implies T
ATR1 y ¼ ðAQ1 Þ y ¼ QT1 AT y > 0: Consequently, ATR1 , which is a Z-matrix, is an M-matrix. Then AR1 is an M-matrix. Since for a Z-matrix the statement ‘‘A is an M-matrix’’ is equivalent to the statement ‘‘all the principal minors of A are positive’’ [6, p. 134], then the submatrix R14 is also an M-matrix. Similarly, we could use the same way for the case m = 1 step by step to prove the results for 2 6 m 6 n 1. The proof is completed. h
620
Y. Zhang et al. / Applied Mathematics and Computation 188 (2007) 612–633
Lemma 3. If A is an M-matrix, then Lm3 ¼ 0; Rm2 ¼ 0; diagðLm1 Þ ¼ diagðRm1 Þ; Lm4 ¼ Rm4 :
ð30Þ ð31Þ ð32Þ ð33Þ
Proof. Since A is an M-matrix, from Lemmas 1 and 2, there exist ALm and ARm for any 1 6 m 6 n 1. From (17) and (19), we immediately get Lm3 = 0 and Rm2 = 0. And from (21), (22) and (24), we have diagðLm1 Þ ¼ diagðRm1 Þ;
ð34Þ
Lm4 ¼ Rm4 :
ð35Þ
Lemma 4. If A is an M-matrix, then DL1 ¼ DR1 ; DL4 ¼ DR4 ; EL4 ¼ ER4 ; F L4 ¼ F R4 :
ð36Þ ð37Þ ð38Þ ð39Þ
Proof. Since A is an M-matrix, from Lemmas 1 and 2, there exist ALm and ARm for any 1 6 m 6 n 1. From the Eq. (21), the result follows. h Lemma 5. The existence of Gauss preconditioners Pm and Qm is guaranteed by the same condition (1) a11 5 0 for m = 1; ðl;m2Þ ðl;m2Þ ðl;m2Þ ðr;m2Þ ðr;m2Þ ðr;m2Þ ðr;m2Þ (2) aðl;m2Þ am1;m1 6¼ am;m1 am1;m or am;m am1;m1 6¼ am;m1 am1;m for 2 6 m 6 n 1. m;m Proof. Notice that, for any 1 6 m 6 n 1, two kinds of Gauss preconditioners Pm and Qm with order m exists if and only if the fore mentioned preconditioned matrices AL,m1 and AR,m1 have the nonzero diagonal elements i.e., ðl;m1Þ
ðr;m1Þ
am1;m1 ¼ am1;m1 6¼ 0: For the case m = 1, it only need a11 5 0 which keep the existence of P1 and Q1. For the case m = 2, it only needs M1A and AM 1 having nonzero diagonal elements, that is to say a22
a21 a12 6¼ 0; a11
i.e. a22a11 5 a21a12. Totally, for any 2 6 m 6 n 1, left Gauss preconditioners Pm with order m exist if and only if ðl;m2Þ ðl;m2Þ
¼ aðl;m2Þ aðl;m1Þ m;m m;m
am;m1 am1;m ðl;m2Þ
am1;m1
ðl;m2Þ
6¼ 0;
ð40Þ
ðl;m2Þ ðl;m2Þ
which implies aðl;m2Þ am1;m1 6¼ am;m1 am1;m . m;m For the same m as Pm, right Gauss preconditioners Qm with order m exists if and only if ðr;m2Þ ðr;m2Þ
ðr;m1Þ ðr;m2Þ ¼ am;m am;m
am;m1 am1;m ðr;m2Þ
am1;m1
ðr;m2Þ
6¼ 0;
ðr;m2Þ ðr;m2Þ
ðr;m2Þ which implies am;m am1;m1 6¼ am;m1 am1;m .
ð41Þ
Y. Zhang et al. / Applied Mathematics and Computation 188 (2007) 612–633
621
In terms of (21), we obtain ðr;m2Þ aðl;m2Þ ¼ am;m ; m;m ðl;m2Þ
ðr;m2Þ
am1;m1 ¼ am1;m1 ; ðl;m2Þ
ðr;m2Þ
ðl;m2Þ
ðr;m2Þ
ð42Þ
am;m1 ¼ am;m1 ; am1;m ¼ am1;m : Thus, by expression (40), expression (41) and (42), the proof is completed. h For convenience, the following notations are needed: • B(M): The Jacobi iteration matrix of the matrix M. • T(M): The Gauss–Seidel iteration matrix of the matrix M. • BLm4, BRm4: the Jacobi iteration matrices of Lm4, Rm4, respectively, where, BLm4 ¼ D1 L4 ðE L4 þ F L4 Þ;
ð43Þ
D1 R4 ðER4
ð44Þ
BRm4 ¼
þ F R4 Þ:
• TLm4, TRm4: the Gauss–Seidel iterative matices of Lm4, Rm4, respectively, where, T Lm4 ¼ ðDL4 EL4 Þ1 F L4 ; 1
T Rm4 ¼ ðDR4 ER4 Þ F R4 :
ð45Þ ð46Þ
Lemma 6. If A is an M-matrix, then (1) BLm4 = BRm4; (2) TLm4 = TRm4. Proof. Since A is an M-matrix, from Lemmas 1 and 2, there exist ALm and ARm for any 1 6 m 6 n 1. ALm 1 and ARm are M-matrices with nonzero diagonal elements, then diagonal matrices D1 L4 and DR4 exist. From Lemma 4, the Eqs. (43)–(46), the results are followed. h Now, we will discuss the properties of Gauss type preconditioners associated with Jacobi, Gauss–Seidel iterative methods, mainly including: (1) The preconditioned matrices of same order left and right Gauss preconditioning methods having the same spectral radii, i.e. q(ARm) = q(ALm). (2) The Jacobi iteration matrices of preconditioned matrices of same order left and right Gauss preconditioning methods having the same spectral radii, i.e. q(BLm) = q(BRm). (3) The Gauss–Seidel iteration matrices of preconditioned matrices of same order left and right Gauss preconditioning methods having the same spectral radii, i.e. q(TLm) = q(TRm). (4) The Jacobi iteration matrices of left (right) Gauss type preconditioned system and its submatrix Lm4 (Rm4) having the same spectral radii, i.e. q(BLm) = q(BLm4) (q(BRm) = q(BRm4)). (5) The Gauss–Seidel iteration matrices of left (right) Gauss type preconditioned system and its submatrix Lm4 (Rm4) having the same spectral radii, i.e. q(TLm) = q(TLm4) (q(TRm) = q(TRm4)). (6) Spectral radii properties on the Jacobi, Gauss–Seidel iteration matrices applied with different order left (right) Gauss preconditioners. Theorem 1. If A is an M-matrix, then qðALm Þ ¼ qðARm Þ;
i:e:; qðP m AÞ ¼ qðAQm Þ:
ð47Þ
622
Y. Zhang et al. / Applied Mathematics and Computation 188 (2007) 612–633
Proof. For A being an M-matrix, from Lemmas 1 and 2, there exist ALm and ARm for any 1 6 m 6 n 1. !! F L2 DL1 F L1 detðkI ALm Þ ¼ det kI 0 DL4 EL4 F L4 !! DL1 0 ¼ det kI 0 DL4 EL4 F L4 !! DL1 0 ¼ det kI ð48Þ 0 Lm4 and DR1 EL1 0 detðkI ARm Þ ¼ det kI ER3 DR4 ER4 F R4 DR1 0 ¼ det kI 0 DR4 ER4 F R4 DR1 0 ¼ det kI : 0 Rm4
ð49Þ
From Lemmas 3 and 4, we have that DL1 = DR1 and Lm4 = Rm4. Then the two equations det(kI ALm) = 0 and det(kI ARm) = 0 have the same characteristic roots which implies qðARm Þ ¼ qðALm Þ; i.e. q(PmA) = q(AQm).
h
Theorem 2. If A is an M-matrix, then qðBLm Þ ¼ qðBRm Þ;
ð50Þ
qðT Rm Þ ¼ qðT Lm Þ:
ð51Þ
Proof. For A being an M-matrix, from Lemmas 1 and 2, there exist ALm and ARm for any 1 6 m 6 n 1. ALm 1 1 and ARm are M-matrices with nonzero diagonal elements, then the four diagonal matrices D1 L1 , DL4 , DR1 and 1 DR4 exist. (1) detðkI BLm Þ ¼ det kI
D1 L1 F L1 0
D1 L1 F L2 1 DL4 ðEL4 þ F L4 Þ
D1 R1 ER1 D1 R4 ER3
0
Similarly, we can get detðkI BRm Þ ¼ det kI
D1 R4 ðE R4 þ F R4 Þ
!!
!!
0 0 : ¼ det kI 0 D1 L4 ðE L4 þ F L4 Þ
ð52Þ
0 ¼ det kI 0
ð53Þ
0 D1 R4 ðE R4 þ F R4 Þ
:
From Lemma 4, it gives DL4 = DR4, EL4 = ER4 and FL4 = FR4. Then detðkI BLm Þ ¼ detðkI BRm Þ: Hence, the two equations det(kI BLm) = 0 and det(kI BRm) = 0 have the same characteristic roots. So qðBRm Þ ¼ qðBLm Þ: (2) Since the triangular matrices DL4 EL4 and DR4 ER4 have nonzero diagonal elements, (DL4 EL4)1 and (DR4 ER4)1 exist. From Lemma 4, we have
Y. Zhang et al. / Applied Mathematics and Computation 188 (2007) 612–633 1
1
ðDR4 ER4 Þ F R4 ¼ ðDL4 EL4 Þ F L4 :
623
ð54Þ
Defining 1
1
T 4 ¼ ðDR4 ER4 Þ F R4 ¼ ðDL4 EL4 Þ F L4 ; then we can write TLm and TRm as ! D1 D1 L1 F L1 L1 F L2 T Lm ¼ ; 0 T4 0 0 T Rm ¼ 1 0 ðDR4 ER4 Þ F R4 0 0 ¼ : 0 T4
ð55Þ
ð56Þ
1 Since diagðD1 L1 F L1 Þ ¼ 0 and DL1 F L1 is a strictly upper triangular matrix, we have 0 0 detðkI T Lm Þ ¼ det kI ¼ detðkI T Rm Þ: 0 T4
Then the two equations det(kI TLm) = 0 and det(kI TRm) = 0 have the same roots. So qðT Lm Þ ¼ qðT Rm Þ:
Remark 1. By Theorem 2, we know the same order left and right Gauss preconditioners have the same convergence rate in Jacobi and Gauss–Seidel iterative methods, respectively. Theorem 3. If A is an M-matrix, then qðBLm Þ ¼ qðBLm4 Þ;
ð57Þ
qðT Lm Þ ¼ qðT Lm4 Þ:
ð58Þ
Proof. (1) We know the Jacobi iteration matrices of ALm and Lm4 are BLm and BLm4. For A being an M-matrix, from Lemmas 1 and 2, there exist ALm and ARm for any 1 6 m 6 n 1. ALm and ARm are M-matri1 ces with nonzero diagonal elements, then D1 Lm , DL4 exist. Hence, BLm and BLm4 exists. By (52), we have 0 0 detðkI BLm Þ ¼ det kI ¼ kk detðkI nk D1 L4 ðEL4 þ F L4 ÞÞ ðE þ F Þ 0 D1 L4 L4 L4 k ¼ kk detðkI nk D1 L4 ðEL4 þ F L4 ÞÞ ¼ k detðkI nk BLm4 Þ:
Then det(kI BLm) = 0 and det(kInk BLm4) = 0 have the same nonzero characteristic roots. So q(BLm) = q(BLm4). (2) For the matrix DL4 EL4 being a lower triangular matrix and the diagonal elements are all nonzero, then (DL4 EL4)1 exists. By (55), we have !! D1 D1 L1 F L1 L1 F L2 detðkI T Lm Þ ¼ det kI : 1 0 ðDL4 EL4 Þ F L4 Since the matrix D1 L1 F L1 is a strictly upper triangular matrix, so we have 0 0 ¼ kk detðkI nk ðDL4 EL4 Þ1 F L4 Þ ¼ kk detðkI nk T Lm4 Þ; detðkI T Lm Þ ¼ det kI 1 0 ðDL4 EL4 Þ F L4 ð59Þ where Ink is the identity matrix with order n k.
624
Y. Zhang et al. / Applied Mathematics and Computation 188 (2007) 612–633
Then we conclude that det(kI TLm) = 0 and det(kInk TLm4) = 0 have the same nonzero eigenvalues. Then q(TLm) = q(TLm4). The proof is completed. h Theorem 4. If A is an M-matrix, then qðBRm Þ ¼ qðBRm4 Þ;
ð60Þ
qðT Rm Þ ¼ qðT Rm4 Þ:
ð61Þ
Proof. The result (60) is an immediate consequence of Eqs. (50), (57) and Lemma 6. The result (61) is an immediate consequence of Eqs. (51), (58) and Lemma 6.
h
Theorem 5. If A is an M-matrix, then for any 1 6 j 6 i 6 n 1, 0 ¼ qðBL;n1 Þ 6 qðBLi Þ 6 qðBLj Þ 6 qðBL1 Þ 6 qðBÞ 6 1;
ð62Þ
0 ¼ qðBR;n1 Þ 6 qðBRi Þ 6 qðBRj Þ 6 qðBR1 Þ 6 qðBÞ 6 1;
ð63Þ
0 ¼ qðT L;n1 Þ 6 qðT Li Þ 6 qðT Lj Þ 6 qðT L1 Þ 6 qðT Þ 6 1;
ð64Þ
0 ¼ qðT R;n1 Þ 6 qðT Ri Þ 6 qðT Rj Þ 6 qðT R1 Þ 6 qðT Þ 6 1:
ð65Þ
Proof. For A being an M-matrix, from Lemmas 1 and 2, there exist ALm and ARm for any 1 6 m 6 n 1 and ALm, ARm, Lm4 and Rm4 are M-matrices. For A being an M-matrix, from [1, Theorem 2.1], we have qðBÞ 6 1;
qðT Þ 6 1:
For (62): For m = 1, by [1, Theorem 2.1], we have qðBLm Þ ¼ qðBL1 Þ 6 qðBÞ: For m = 2, we still use the Theorem 2.1 in [1] for the submatrix L14 with order n 1, then qðBðL24 ÞÞ 6 qðBðL14 ÞÞ; i.e. qðBLm4 Þ 6 qðBLðm1Þ4 Þ:
ð66Þ
By Theorem 3, we have qðBL1 Þ ¼ qðBLðm1Þ Þ ¼ qðBLðm1Þ4 Þ
ð67Þ
qðBL2 Þ ¼ qðBLm Þ ¼ qðBLm4 Þ:
ð68Þ
and Then, by (66)–(68), we have q(BL2) 6 q(BL1), i.e., qðBLm Þ 6 qðBLðm1Þ Þ;
m ¼ 2:
For 3 6 m 6 n 1: we use the same method proving the cases m = 2 and m = 1 step by step, then qðBL;n1 Þ 6 qðBLi Þ 6 qðBLj Þ 6 qðBL1 Þ 6 qðBÞ 6 1;
m 6 j 6 i 6 n 1:
In addition, since BL,n1 is an upper triangular matrix, then q(BL,n1) = 0. So the result (62) follows. For (63): By Theorem 2, we know that qðBLm Þ ¼ qðBRm Þ;
m ¼ 1; 2; . . . ; n 1:
Then, by (62), (63) follows. For (64): Since the result of Theorem 2.1 in [1] hold for Gauss–Seidel iterative method, we can use the same proving method used to prove the result (62), to prove the result (64).
Y. Zhang et al. / Applied Mathematics and Computation 188 (2007) 612–633
625
For (65): By Theorem 2, we know that qðT Lm Þ ¼ qðT Rm Þ;
m ¼ 1; 2; . . . ; n 1:
Then by the result (64), the result (65) follows. h It should be mentioned that Theorem 5 theoretically confirm the statistical results in [3]. 3. Preconditioning algorithms In this section we study the construction of the left and right Gauss preconditioners in detail. Matlab-like notation is applied in algorithm descriptions. Firstly, we discuss the construction of two kinds of Gauss transformation matrices which are essential to the corresponding Gauss preconditioners. For brevity, we denote one function by GetML with two input parameters A and k and a return parameter Mk which implement the computation to get the matrix Mk by expression (10). Similarly, the function GetMR with the same input parameters as GetML and a return parameter M k is used to implement the computation to get the matrix M k by (13). So these two procedures used to compute the Gauss transformation matrices are given in Algorithms 1 and 2, respectively. Algorithm 1 (Construction of left Gauss transformation matrices) Input: A, m Output: Mm 1. for k = 1 to m, 2. M = GetML(A, k) 3. A=M*A 4. end for 5. Mm = M Algorithm 2 (Construction of right Gauss transformation matrices) Input: A, m Output: M m 1. for k = 1 to m, 2. M = GetMR(A, k) 3. A=A*M 4. end for 5. M m ¼ M Secondly, the normal algorithms computing Pm and Qm are given in the following. We note that Pm = MmMm1 . . . M2M1 and Qm ¼ M 1 M 2 M m1 M m . Algorithm 3 [3], (Constructing left Gauss type preconditioners) Input: A, m Output: P 1. Initialize: 2. P = I 3. T = A 4. Computation for Pm: 5. for k = 1 to m, 6. Set v be zero row vector. 7. v(k + 1: n) = T(k + 1: n, k)/T(k, k)
626
8. 9. 10. 11. 12. 13. 14. 15. 16.
Y. Zhang et al. / Applied Mathematics and Computation 188 (2007) 612–633
if k = 1 then P(k + 1: n, 1) = v(k + 1: n) else for r = k + 1 to n, P(r, 1: k) = P(r, 1: k) + v(r) * P(k, 1: k) end for end if T=P*A end for.
The algorithm computing the right Gauss type preconditioners is based on expression (12) as follows. Algorithm 4 (Constructing right Gauss type preconditioner) Input: A, m Output: Qm 1. Initialize: 2. P = I 3. T = AT 4. Computation for Pm: 5. for k = 1 to m, 6. Set v be n dimension zero column vector. 7. v(k + 1: n) = T(k + 1: n, k)/T(k, k) 8. if k = 1 then 9. P(k + 1: n, 1) = v(k + 1: n) 10. else 11. for r = k + 1 to n, 12. P(r, 1: k) = P(r,1: k) + v(r) * P(k, 1: k) 13. end for 14. end if 15. T = P * AT 16. end for. 17. Qm = PT We note that in general these two kinds of Gauss preconditioners Pm and Qm need not be constructed if we could get the left preconditioned matrices on ALm = PmA, bLm = Pmb for left Gauss preconditioning techniques and the right preconditioned matrices ARm = AQm and Qm. After the left Gauss preconditioning it remains to solve the preconditioned system ALmx = bLm. In addition, after the right Gauss preconditioning, it remains to solve the preconditioned system ARmy = b and x = Qmy. Thirdly, therefore, we could give the improved Gauss type preconditioning algorithms by expression (14) and expression (15) to the preconditioned linear systems. Algorithm 5 (Inbuilt left Gauss type preconditioning techniques) Input: A, b, m Output: ALm, bL, Pm 1. Set n be the order of A 2. P = I 3. for k = 1 to m 4. Set v be n dimension zero column vector. 5. v(k + 1: n) = A(k + 1: n, k)/A(k, k) 6. for r = k + 1 to n 7. A(r, k + 1: n) = A(r, k + 1: n) + v(r) * A(k, k + 1: n)
Y. Zhang et al. / Applied Mathematics and Computation 188 (2007) 612–633
8. 9. 10. 11. 12. 13. 14. 15. 16. 17.
627
b(r) = b(r) + v(r) * b(k) end for A(k + 1: n, k) = 0 for i = 1 to k P(k + 1: n, i) = P(k + 1: n, i) + P(k, i) * v(k + 1: n) end for end for. ALm = A bLm = b Pm = P
The lines 2, 11 to 13, 17 in Algorithm 5 are only used to show the computation procedure of the left preconditioners Pm. We should delete these lines in the actual application. Note that, for m = n 1, Algorithm 5 constructing the left preconditioner of order n 1 become a procedure of Gauss elimination. Algorithm 6 (Inbuilt right Gauss type preconditioning techniques) Input: A, m Output: ARm, Qm 1. Set n be the order of A. 2. Q = I 3. for k = 1 to m 4. Set v be n dimension zero row vector. 5. v(k + 1: n) = A(k, k + 1: n)/A(k, k) 6. for r = k + 1 to n 7. A(k + 1: n, r) = A(k + 1: n, r) + v(r) * A(k + 1: n, k) 8. end for 9. A(k, k + 1: n) = 0 10. for i = 1 to k 11. Q(i, k + 1: n) = Q(i, k + 1: n) + Q(i, k) * v(k + 1: n) 12. end for 13. end for. 14. ARm = A 15. Qm = Q Finally, we study floating point operations of Algorithms 5 and 6. Proposition 1. The number of floating point operations of Algorithm 5 (m order left Gauss preconditioning technique) is 1 mð4m2 mð3 þ 12nÞ þ 12n2 þ 6n 7Þ: 6
ð69Þ
Proof. For computing Gauss vector v, the number of floating point operations in line 5, denoted by Nv, is N v ¼ n k: For computing preconditioned matrix ALm, the number of floating point operations in lines 6–9 except line 8, denoted by NA, is n X 2 2ðn kÞ ¼ 2ðn kÞ : NA ¼ r¼kþ1
For the computing preconditioned right hand vector b, the number of floating point operations in lines 6–9 except line 7, denoted by Nb, is
628
Y. Zhang et al. / Applied Mathematics and Computation 188 (2007) 612–633
N b ¼ 2ðn kÞ: Now the total cost of the kth step of the main loop is N k ¼ N v þ N A þ N b:
ð70Þ
Summing Eq. (70) over k = 1 ! m yields the total number of floating point operation of Algorithm 5 as m m m X X X Nk ¼ Nv þ NA þ Nb ¼ ½ðn kÞ þ 2ðn kÞ2 þ 2ðn kÞ k¼1
k¼1
k¼1
¼ mð2n2 þ 3nÞ ð4n þ 3Þ
mðm þ 1Þ 2mðm þ 1Þð2m þ 1Þ þ 2 6
1 ¼ mð4m2 mð3 þ 12nÞ þ 12n2 þ 6n 7Þ: 6
ð71Þ
Proposition 2. The number of floating point operations of Algorithm 5 (left Gauss preconditioning technique) with m = n 1 is 2 3 1 2 7 n þ n n: 3 2 6
ð72Þ
Proof. Substituting m = n 1 in (69) gives n1 X k¼1
2 1 7 N k ¼ n3 þ n2 n: 3 2 6
ð73Þ
Proposition 3. The number of floating point operations of Algorithm 6 (m order right Gauss preconditioning technique) is 2n2 m nm2
m m2 : 2 2
ð74Þ
Proof. For computing Gauss vector v, the number of floating point operations in line 5, denoted by Nv, is N v ¼ n k: For computing preconditioned matrix ARm, the number of floating point operations in lines 6–8, denoted by NA, is n X 2 2ðn kÞ ¼ 2ðn kÞ : NA ¼ r¼kþ1
For computing preconditioner Q, the number of floating point operations in lines 10–12, denoted by NQ, is N Q ¼ 2kðn kÞ: Now the total cost of the kth step of the main loop is N k ¼ N v þ N A þ N Q:
ð75Þ
Summing Eq. (75) over k = 1 ! m yields the total number of floating point operation of Algorithm 5 as m m m X X X 2 Nk ¼ Nv þ NA þ NQ ¼ ½ðn kÞ þ 2ðn kÞ þ 2kðn kÞ k¼1
k¼1
k¼1
m X mðm þ 1Þ ½ðn þ 2n2 Þ ð1 þ 2nÞk ¼ nð1 þ 2nÞm ð1 þ 2nÞ ¼ 2 k¼1 2 mþ1 m m ¼ mð1 þ 2nÞ n ¼ 2n2 m nm2 : 2 2 2
ð76Þ
Y. Zhang et al. / Applied Mathematics and Computation 188 (2007) 612–633
629
Proposition 4. The number of floating point operations of Algorithm 6 (right Gauss preconditioning technique) with m = n 1 is n3 12 n2 12 n. Proof. Substituting m = n 1 in (74) gives n1 X k¼1
Nk ¼
nð1 þ 2nÞðn 1Þ n2 n ¼ n3 : 2 2 2
Proposition 5. The number of floating point operations of Algorithm 5 (left Gauss preconditioning technique) is always less than that of Algorithm 6 (right Gauss preconditioning technique) for any m 6 n 1 when n P 2. Proof. Subtracting (69) by (74), for any m 6 n 1, then the result, denoted by Nd, yields 1 1 1 N d ¼ mð5 3m þ 2m2 þ 6n 6mnÞ ¼ mð2mðm nÞ þ ð6n 3m 5 4mnÞÞ < mð6n 4mnÞ: 6 6 6 ð77Þ For the case m = 1, we have Nd = 1 < 0. For the cases m P 2, we have Nd < 0 immediately from (77). So, the result is proved. h Remark 2. Comparing Propositions 2 and 4 shows that Algorithm 5 costs n3 ðn2 3n þ 2Þ less floating point operations. Therefore, we conclude that the left Gauss preconditioning techniques are superior to the right Gauss preconditioning techniques concerning the proposed algorithms. Numerical results of comparisons on the flops of Algorithms 5 and 6 are given in Figs. 1 and 2. 4. Numerical experiments Numerical experiments are presented to see the performance of these inbuilt Gauss type preconditioning algorithms (Algorithms 5 and 6) associated with the Jacobi, Gauss–Seidel iterative methods. In all the numerical methods we tested, the initial approximation x(0) is taken as the zero-vector, and the right hand side vector b is chosen so that x ¼ ð1; 1; . . . ; 1ÞT is the solution vector for the considered linear system. We have used kb Ax(k)k/k bk 6 106 as the stopping criteria. The max iteration in iterative methods is 8
10
x 10
Left Gauss Right Gauss 9
The number of floating point operations
8
7
6
5
4
3
2
1
0
0
100
200
300
400
500 m (n=1000)
600
700
800
900
1000
Fig. 1. Comparison on flops of Algorithms 5 and 6 on the case n = 1000.
630
Y. Zhang et al. / Applied Mathematics and Computation 188 (2007) 612–633 11
10
x 10
Left Gauss Right Gauss 9
The number of floating point operations
8
7
6
5
4
3
2
1
0
0
1000
2000
3000
4000
5000 6000 m (n=10000)
7000
8000
9000
10000
Fig. 2. Comparison on flops of Algorithms 5 and 6 on the case n = 10,000.
limited to 500. All experiments were executed on a PC using the Matlab programming package. Matlab carries out all calculations in double precision by default. In all tables, we report the spectral radius of the corresponding iteration matrix (Radius), the time in seconds for computing the preconditioner (P-time), the time for the iterative part only (It-time), and the number of iterations (Its). Left (1), Right (1), Left (2) and Right (2) are used to denote Algorithms 3–6 respectively. Example 1. Consider the two-dimensional elliptic partial differential equation (Helmholtz) [9]
o2 u o2 u þ uðx; yÞ ¼ 0; ox2 oy 2
on the rectangular region R defined by the inequalities 0 6 x 6 2 and 0 6 y 6 1 with the Dirichlet boundary conditions. We have used two uniform meshes of Dx = 2/21, Dy = 1/21 which leads to a matrix of order n = 20 · 20 where Dx and Dy refer to the mesh sizes in the x-direction and y-direction, respectively. Table 1 Comparisons of left Gauss preconditioning methods with Jacobi iterative method for Example 1 Preconditioner
Classical GS(1) GS(10) GS(50) GS(90) GS(130) GS(170) GS(210) GS(250) GS(290) GS(330) GS(370) GS(399)
Radius
0.987935 0.987935 0.987921 0.987719 0.987074 0.985779 0.983534 0.979719 0.972898 0.959152 0.924344 0.779000 0.302955
Its
815 815 814 800 762 698 611 505 388 266 150 64 43
‘‘Classical’’: the unpreconditioned iterative method.
Left (1)
Left (2)
P-time
It-time
P-time
It-time
0.13 0.81 3.84 6.95 10.03 13.14 16.22 19.31 22.34 25.34 28.27 30.36
7.02 7.34 7.31 7.13 6.83 6.30 5.47 4.64 3.61 2.42 1.41 0.64 0.45
0.02 0.20 0.92 1.45 1.89 2.30 2.53 2.73 2.89 2.92 2.94 3.00
7.02 6.87 6.78 6.73 6.34 5.89 5.08 4.25 3.27 2.22 1.33 0.55 0.36
Y. Zhang et al. / Applied Mathematics and Computation 188 (2007) 612–633
631
Note that in Tables 1–4, for the highest order left and right Gauss preconditioned linear system, here we use the classical Jacobi and Gauss–Seidel to solve the triangular linear systems. Table 2 Comparisons of right Gauss preconditioning methods with Jacobi iterative method for Example 1 Preconditioner
Radius
Its
Classical GS(1) GS(10) GS(50) GS(90) GS(130) GS(170) GS(210) GS(250) GS(290) GS(330) GS(370) GS(399)
0.987935 0.987935 0.987921 0.987719 0.987074 0.985779 0.983534 0.979719 0.972898 0.959152 0.924344 0.779000 0.302955
815 815 814 802 764 700 614 509 391 270 154 64 43
Left (1)
Left (2)
P-time
It-time
P-time
It-time
0.14 0.83 3.88 6.95 10.05 13.17 16.25 19.33 22.38 25.38 28.30 30.39
7.22 7.00 7.02 6.88 6.50 6.05 5.25 4.39 3.41 2.34 1.38 0.61 0.44
0.02 0.11 0.55 0.98 1.42 1.83 2.17 2.61 2.88 3.19 3.41 3.53
7.22 6.94 6.88 6.81 6.44 5.94 5.20 4.33 3.45 2.28 1.31 0.55 0.38
Table 3 Comparisons of left Gauss preconditioning methods with Gauss–Seidel iterative method for Example 1 Preconditioner
Radius
Its
Classical GS(1) GS(10) GS(50) GS(90) GS(130) GS(170) GS(210) GS(250) GS(290) GS(330) GS(370) GS(399)
0.976015 0.976015 0.975988 0.975590 0.974315 0.971762 0.967342 0.959858 0.946551 0.920040 0.854735 0.611080 0.297231
437 437 437 430 410 376 329 273 210 144 87 49 43
Left (1)
Left (2)
P-time
It-time
P-time
It-time
0.14 0.81 3.88 6.94 10.95 13.19 16.27 19.34 22.39 25.33 28.34 30.41
7.06 6.81 6.86 6.55 6.27 6.09 5.02 4.19 3.28 2.25 1.38 0.83 0.73
0.03 0.20 0.89 1.52 2.02 2.30 2.52 2.73 2.89 2.94 2.95 3.00
7.06 6.66 6.69 6.53 7.06 7.45 5.03 4.17 3.22 2.16 1.30 0.75 0.66
Table 4 Comparisons of right Gauss preconditioning methods with Gauss–Seidel iterative method for Example 1 Preconditioner
Classical GS(1) GS(10) GS(50) GS(90) GS(130) GS(170) GS(210) GS(250) GS(290) GS(330) GS(370) GS(399)
Radius
0.976015 0.976015 0.975988 0.975590 0.974315 0.971762 0.967342 0.959858 0.946551 0.920040 0.854735 0.611080 0.000000
Its
437 437 437 430 410 376 329 273 210 144 82 30 2
Left (1)
Left (2)
P-time
It-time
P-time
It-time
0.14 0.81 3.88 6.95 10.06 13.20 16.30 19.38 22.44 25.47 28.73 30.67
6.92 6.72 6.80 6.69 6.39 5.88 5.13 4.27 3.30 2.27 1.33 0.53 0.09
0.02 0.11 0.56 1.00 1.45 1.84 2.19 2.61 2.89 3.19 3.38 3.52
6.92 6.72 6.73 6.61 6.30 5.83 5.06 4.19 3.22 2.20 1.27 0.47 0.03
632
Y. Zhang et al. / Applied Mathematics and Computation 188 (2007) 612–633
From Tables 1–4, we find that the improved inbuilt Gauss preconditioning algorithms (Algorithms 5 and 6) are greatly superior to the non-inbuilt ones (Algorithms 2 and 4). Example 2. Problems set called LANPRO (Lanczos with partial reorthogonalization) from the Harwell-Boeing Collection. Problems matrices Nos1, Nos2 and Nos5 all with nonzero diagonal elements are involved in the experiments. Note that the highest order preconditioned systems are solved by Gauss–Seidel iterative method and direct triangular solving method (Marked with ‘‘Tri’’), respectively. Numerical results are given in Tables 5–7.
Table 5 Comparisons of left Gauss Preconditioning methods with Gauss–Seidel iterative method (Example 2. Nos1 with order n = 237) Preconditioner
Radius
Its
Left (1)
Left (2)
P-time
It-time
Min
Max
Classical GS(1) GS(10) GS(57) GS(104) GS(151) GS(198) GS(236) GS(236)Tri
0.999999 0.999999 0.999999 0.999999 0.999998 0.999996 0.999945 0.000000 0.000000
500 500 500 500 500 500 500 91 –
– 0.02 0.13 0.41 0.74 0.69 0.84 0.97 0.97
3.92 3.62 3.66 3.70 4.02 4.00 3.86 0.70 0.02
0.240140 0.240140 0.124329 0.010632 0.010627 0.010178 0.002500 1.000000 1.000000
2.242217 2.242217 2.242217 2.242217 2.242220 2.242260 2.244213 1.000000 1.000000
‘‘Min’’ : the min entry of the final numerical solution vector. ‘‘Max’’ : the max entry of the final numerical solution vector.
Table 6 Comparisons of Left Gauss preconditioning methods with Gauss–Seidel iterative method (Example 2. Nos5 with order n = 468) Preconditioner
Radius
Its
Left (1)
Left (2)
P-time
It-time
Min
Max
Classical GS(1) GS(10) GS(103) GS(196) GS(289) GS(382) GS(467) GS(467)Tri
0.997904 0.997904 0.997903 0.997893 0.997851 0.997744 0.996816 0.000000 0.000000
500 500 500 500 500 500 500 33 –
– 0.09 0.34 2.39 4.11 4.45 4.78 4.84 4.84
10.55 10.19 10.20 10.23 10.20 10.16 10.25 0.67 0.02
0.417626 0.417129 0.416974 0.413352 0.386359 0.440179 0.650278 1.000000 1.000000
1.008461 1.008470 1.008475 1.008573 1.017740 1.016211 1.053018 1.000000 1.000000
Table 7 Comparisons of left Gauss Preconditioning methods with Gauss–Seidel iterative method (Example 2. Nos2 with order n = 957) Preconditioner
Radius
Its
Classical GS(1) GS(10) GS(201) GS(392) GS(583) GS(774) GS(956) GS(956)Tri
1.000000 1.000000 1.000000 1.000000 1.000000 1.000000 1.000000 0.000000 0.000000
500 500 500 500 500 500 500 334 –
Left (1)
Left (2)
P-time
It-time
Min
Max
– 0.19 1.22 17.02 29.45 34.28 35.81 36.33 36.33
32.88 31.64 31.80 34.39 32.25 31.88 31.44 22.69 0.03
4.208663 4.208663 3.746223 0.061133 0.061133 0.061133 0.061133 1.000000 1.000000
6.216589 6.216589 6.216589 6.216589 6.216589 6.216589 6.216589 1.000000 1.000000
Y. Zhang et al. / Applied Mathematics and Computation 188 (2007) 612–633
633
Finally, this paper presents some methods and numerical improvements for preconditioning. Parameterized m-order Gauss preconditioning technique is a natural problem. Further research is required. References [1] A. Hadjidimos, D. Noutsos, M. Tzoumas, More on modifications and improvements of classical iterative schemes for M-matrices, Linear Algebra Appl. 364 (2003) 253–279. [2] J.P. Milaszewicz, Improving Jacobi and Gauss–Seidel iterations, Linear Algebra Appl. 93 (1987) 161–170. [3] Y. Zhang, T.Z. Huang, X.P. Liu, Modified iterative methods for nonnegative matrices and M-matrices linear systems, Comput. Math. Appl. 50 (2005) 1587–1602. [4] G.H. Golub, C.F. Van Loan, Matrix Computations, third ed., The Johns Hopkins University Press, 1996. [5] D.J. Evans, M.M. Martins, M.E. Trigo, The AOR iterative method for new preconditioned linear systems, J. Comput. Appl. Math. 132 (2001) 461–466. [6] A. Berman, P.J. Plemmons, Nonnegative Matrices in the Mathematical Sciences, Academic Press, 1979. [7] T. Kohno, H. Kotakemori, H. Niki, M. Usui, Improving the Gauss–Seidel method for Z-matrices, Linear Algebra Appl. 267 (1997) 113–123. [8] A.D. Gunawardena, S.K. Jain, L. Snyder, Modified iterative methods for consistent linear systems, Linear Algebra Appl. 154–156 (1991) 123–143. ¨ zlu¨, On the numerical solution of Helmholtz equation by alternating group explicit (AGE) methods, Appl. Math. [9] Necdet Bildik, Sevil O Comput. 163 (2005) 505–518.
Gauss type preconditioning techniques for linear systems Yong Zhang a, Ting-Zhu Huang
a,*
, Xing-Ping Liu
q
b
a
b
School of Applied Mathematics, University of Electronic Science and Technology of China, Chengdu, Sichuan 610054, PR China National Key Lab of Computational Physics, Institute of Applied Physics and Computational Mathematics, Beijing 100088, PR China
Abstract Many researchers have considered preconditioners chosen to eliminate the off-diagonal elements of the coefficient matrix of a linear system. In this work, we generalize the left Gauss type preconditioners [Y. Zhang, T.Z. Huang, X.P. Liu, Modified iterative methods for nonnegative matrices and M-matrices linear systems, Comput. Math. Appl. 50 (2005) 1587–1602] which eliminate the strictly lower triangular elements. Right Gauss type preconditioners that eliminate strictly upper triangular elements are proposed in this paper. These Gauss type preconditioners are partly derived from the LU factorization method. Theoretic analysis on spectral radii of the two kinds of Gauss type preconditioners is given. Numerical experiments are used to show the performance of the improved inbuilt left and right Gauss type preconditioning algorithms associated with Jacobi type and Gauss–Seidel type iterative methods. 2006 Elsevier Inc. All rights reserved. Keywords: Preconditioning; Jacobi and Gauss–Seidel type iterative methods; LU factorization; Gauss elimination; M-Matrix
1. Introduction Consider the numerical solution of the linear system Ax ¼ b; where A = (aij) 2 R
ð1Þ n·n
is a known nonsingular M-matrix.
Definition 1 [6]. A matrix A is a Z-matrix if aij 6 0, i, j = 1, 2, . . . , n; i 5 j. Definition 2 [5]. A matrix A is an L-matrix if aii > 0, i = 1, 2, . . . , n and aij 6 0 for all i, j = 1, 2, . . . , n; i 5 j. Definition 3 [6]. A matrix A is an M-matrix if A has the form A = sI B, s > 0, B P 0 and s > q(B) (the spectral radius of B). q
Supported by NCET of China (NCET-04-0893), foundation of the National Key Lab of Computational Physics, Beijing Appl. Phys. Comput. Math. Institute and applied basic research foundations of Sichuan province (05JY029-068-2). * Corresponding author. E-mail addresses: [email protected] (Y. Zhang), [email protected] (T.-Z. Huang). 0096-3003/$ - see front matter 2006 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2006.10.053
Y. Zhang et al. / Applied Mathematics and Computation 188 (2007) 612–633
613
Since A is an M-matrix, aii > 0, i = 1, 2, . . . , n, which shows that A is an L-matrix [1]. Without loss of generality, we can write A ¼ I L U;
ð2Þ
where I is the identity matrix, L and U are strictly lower triangular and strictly upper triangular parts of A respectively. M = I, N = L + U and M = I L, N = U leads to the classical Jacobi and classical Gauss–Seidel iterative methods respectively. Then the corresponding iteration matrices of the classical Jacobi and Gauss– Seidel methods are given by B = L + U and T = (I L)1U respectively. One can transform the original system (1) into the preconditioned form PAx ¼ Pb
ð3Þ
or AQy ¼ b;
x ¼ Qy:
ð4Þ
Here P in (3) and Q in (4) are called the left preconditioner and the right preconditioner, respectively. Then one can define two kinds of basic iterative schemes respectively ðkÞ xðkþ1Þ ¼ M 1 þ M 1 p N px p Pb;
k ¼ 0; 1; . . . ;
ð5Þ
where PA = Mp Np and Mp is nonsingular, and ðkÞ y ðkþ1Þ ¼ M 1 þ M 1 p N py p b;
k ¼ 0; 1; . . . ;
ð6Þ
where AQ = Mp Np and Mp is nonsingular. Many researchers have considered the left preconditioners applied to system (1) that make the associated Jacobi and Gauss–Seidel type methods converge faster than the original ones. In a simpler form, Milaszewicz [2] considered the preconditioner 1 0 1 C B C B a21 1 C B C B a31 1 C B C B . . C B . . ð7Þ . P1 ¼ B . C; C B C B ai1 1 C B C B . .. C B . . A @ . an1 1 which eliminates the elements of the first column of A below the diagonal. Hadjidimos et al. [1], for M-matrices systems, generalized Milaszewicz’s preconditioning technique and considered the preconditioner 1 0 1 C B C B a2 a21 1 C B C B a3 a31 1 C B C B . . C B . . P 1 ðaÞ ¼ B . ð8Þ . C: C B C B ai ai1 1 C B C B . . C B . . . A @ . an an1
1
Kohno et al. [7] parameterized the preconditioner presented by Gunawardena et al. [8] and considered a preconditioner
614
Y. Zhang et al. / Applied Mathematics and Computation 188 (2007) 612–633
0 B B B B B B B P 2 ðaÞ ¼ B B B B B B @
1
a1 a12 1
1 a2 a23 .. .
..
.
1
ai ai;iþ1 .. .
C C C C C C C C: C C .. C . C C 1 an1 an1;n A 1
ð9Þ
Recently, Zhang et al. [3] proposed the left Gauss type preconditioning techniques based on Hadjidimos et al. [1], Milaszewicz [2] and LU factorization method [4]. The Gauss type preconditioning technique utilizes the Gauss transformation matrices [4] as the base of the Gauss type preconditioner. The construction of Gauss transformation matrices [4] is as follows: 1 0 1 0 0 0 B .. .. C .. .. B. .C . . C B C B 1 0 0C B0 C; k ¼ 1; 2; . . . ; n 1; Mk ¼ B ð10Þ B0 skþ1 1 0C C B B. .. C .. .. C B. @. .A . . 0
sn
0
1
where sT ¼ ð0; . . . ; 0; skþ1 ; . . . ; sn Þ;
si ¼
aik ; i ¼ k þ 1; . . . ; n: akk
The fore k elements of the vector s 2 Rn are zero, and s is called a Gauss vector. In some cases the two preconditioners (7) and (8) have the same expression with one of Gauss transformation matrices, denoted by M1, of LU factorization methods [4], that is to say M 1 ¼ P 1 ¼ P 1 ðaÞ;
T
T
a ¼ ða1 ; a2 ; . . . ; an1 Þ ¼ ð1; 1; . . . ; 1Þ :
Due to the similarity between the preconditioner (7) and (8) and Gauss transformation matrix (10) (i.e. M1), Zhang et al. [3] consider the following left preconditioners: P 1 ¼ M 1; P 2 ¼ M 2M 1; .. . P m ¼ M m M m1 M 2 M 1 :
ð11Þ
Then, • preconditioners Pm (m = 1, 2, . . . , n 1) are all lower triangular matrices with unity diagonal entries; • the entries of the fore m columns of PmA (m = 1, 2, . . . , n 1) are all zero below the diagonal; • Pn1A is an upper triangular matrix in particular. This work is organized as follows. In Section 2, we extend the left Gauss preconditioners by giving a family of preconditioners based on the elimination of elements at the strictly upper triangular part (i.e. up the diagonal), and, theoretic analysis is presented which theoretically confirm the results in [3]. In Section 3, improved Gauss preconditioning algorithms associated with inbuilt technique are given. In Section 4, numerical experiments are presented to show the performance and efficiency of improved Gauss preconditioning algorithms.
Y. Zhang et al. / Applied Mathematics and Computation 188 (2007) 612–633
615
2. Generalized Gauss type preconditioning techniques 2.1. Right Gauss type preconditioners From the fore section, we know that the left Gauss type preconditioner (11) could increase the number of zero elements in the strictly lower triangular part of the preconditioned matrix (PmA) column by column. Similarly, we can construct the preconditioners which eliminate the elements up the diagonal row by row if possible. First, let Qm be such kind of preconditioners with order m. The following expression shows the way to construct Qm: ðAQm ÞT ¼ QTm AT ;
ð12Þ QTm ,
then, in order to find the kind of wanted preconditioners, we only need to construct the matrix just by the fore mentioned left preconditioning techniques which could make the entries of the fore m columns of the strictly lower triangular part of the preconditioned matrix QTm AT all zero. We note that these preconditioning matrices Qm could make the entries of the fore m rows of the strictly upper triangular part of the preconditioned matrix AQm all zero. For brevity, these preconditioners Qm is named as right Gauss preconditioners which could be applied to (1) as right preconditioners. Similarly, we define the corresponding right Gauss transformation matrices 0
1 B. B .. B B B0 Mk ¼ B B0 B B. B. @.
0 .. . 1 0 .. .
1 .. .
0
0
0
0 .. . ukþ1
1 0 .. C . C C C un C C; 0 C C .. C C . A 1
k ¼ 1; 2; . . . ; n 1;
ð13Þ
where uT ¼ ð0; . . . ; 0; ukþ1 ; . . . ; un Þ;
uj ¼
akj ; j ¼ k þ 1; . . . ; n: akk
The fore k elements of u 2 Rn are all 0, and u is also called a Gauss vector. We note that Qm ¼ M 1 M 2 M m1 M m ;
1 6 m 6 n 1:
Now we begin to analyze the elements of the preconditioned matrices ALm and ARm which satisfy ALm ¼ P m A
and ARm ¼ AQm :
First, let us get the expression of the elements of MkA and AM k , which help us to build the corresponding preconditioning algorithms presented in the next section, from the previous notations and considerations. Defining B ¼ ðbij Þ ¼ M k A;
j ¼ 1; 2; . . . ; n;
we have bij ¼ i.e.
( bij ¼
aij ;
i ¼ 1; 2; . . . ; k;
aij si akj ;
i ¼ k þ 1; . . . ; n;
aij ; i ¼ 1; 2; . . . ; k; aik akj aij akk ; i ¼ k þ 1; . . . ; n:
ð14Þ
616
Y. Zhang et al. / Applied Mathematics and Computation 188 (2007) 612–633
Defining B ¼ ð bij Þ ¼ AM k ; we have bij ¼ i.e.
i ¼ 1; 2; . . . ; n; j ¼ 1; 2; . . . ; k;
aij ;
aij uj aik ; j ¼ k þ 1; . . . ; n; (
bij ¼
aij ; j ¼ 1; 2; . . . ; k; aik akj aij akk ; j ¼ k þ 1; . . . ; n:
ð15Þ
Define ðlmÞ
ðrmÞ
ALm ¼ ðaij Þ;
ARm ¼ ðaij Þ;
m ¼ 0; 1; 2; . . . ; n 1;
where ðl0Þ
ðr0Þ
A ¼ AL0 ¼ AR0 ¼ ðaij Þ ¼ ðaij Þ: By the fore known properties of Gauss preconditioners ALm and ARm can be written as 0 ðl;mÞ ðl;mÞ ðl;mÞ ðl;mÞ a12 a1m a1;mþ1 a11 B ðl;mÞ ðl;mÞ ðl;mÞ B 0 a22 a2m a2;mþ1 B B .. .. .. .. B .. B . . . . . B B ðl;mÞ ðl;mÞ 0 0 a a B m;mþ1 mm ALm ¼ B B ðl;mÞ 0 amþ1;mþ1 B 0 B B . .. .. .. .. B .. . . . . B B . .. .. .. .. B . . . . . @ . ¼
0 Lm1
Lm2
Lm3
Lm4
ðl;mÞ
0
an;mþ1
and considerations, the preconditioned matrices ðl;mÞ
a1n
1
C C C C C C C C aðl;mÞ C mn C ðl;mÞ C amþ1;n C C .. C . C C .. C C . A aðl;mÞ nn ðl;mÞ
a2n .. .
ð16Þ
ð17Þ
and 0
ARm
ðr;mÞ
a11
B B ðr;mÞ B a21 B B . B .. B B ðr;mÞ B am1 ¼B B ðr;mÞ B amþ1;1 B B . B . B . B . B . @ . ðr;mÞ
¼
an1 Rm1 Rm3
0 ðr;mÞ
a22 .. .
ðr;mÞ am2 ðr;mÞ amþ1;2
.. . .. .
.. . .. . ðr;mÞ
an2 Rm2 ; Rm4
0 .. .
.. . .. .
.. . .. .
0
0 .. . .. .
aðr;mÞ mm ðr;mÞ amþ1;m
0
ðr;mÞ amþ1;mþ1
.. . .. .
.. . .. .
aðr;mÞ nm
an;mþ1
ðr;mÞ
0 .. . .. .
1
C C C C C C C C 0 C C ðr;mÞ C amþ1;n C C .. C C . C .. C C . A
ð18Þ
aðr;mÞ nn ð19Þ
Y. Zhang et al. / Applied Mathematics and Computation 188 (2007) 612–633
617
where the block form in (17) and (19) is the same, and square matrices Lm1 and Rm1 are all of dimension m · m and square matrices Lm4 and Rm4 are all of dimension (n m) · (n m). Now, we can get the analysis expressions of two kinds of preconditioned matrices in terms of (14) and (15): • For m = 1, we have ðARm Þ11 ¼ ðALm Þ11 ¼ a11 ; ai1 a1j ; a11
ðALm Þij ¼ ðARm Þij ¼ aij
i ¼ m þ 1; . . . ; n; j ¼ m þ 1; . . . ; n:
ð20Þ
• For m = 2, we have ðl;mÞ
ðALm Þ22 ¼ a22
ðl;m1Þ
¼ a22
;
ðr;mÞ
ðr;m1Þ
ðARm Þ22 ¼ a22
¼ a22
ðr;m1Þ
¼ a22 ;
:
By (20) we know ðl;m1Þ
a22
ðl;21Þ
¼ a22
ðl;1Þ
¼ a22 ;
a22
ðr;21Þ
¼ a22
ðr;1Þ
ðl;1Þ
ðr;1Þ
a22 ¼ a22 ;
then (ALm)22 = (ARm)22. • In a general way, by (14), (15), (20), for any 1 6 m 6 n 1, we have ðALm Þij ¼ ðARm Þij ¼
ðl;m1Þ aij
ðl;m1Þ ðl;m1Þ amj ðl;m1Þ amm
aim
¼
ðr;m1Þ aij
ðr;m1Þ ðr;m1Þ amj ðr;m1Þ amm
aim
;
ð21Þ
i = m + 1, . . . , n; j = m + 1, . . . , n. Naturally, from the analysis expression (21), we conclude that (1) Preconditioned matrices ALm and ARm have the same diagonal elements, i.e. ðl;mÞ
ðALm Þii ¼ ðARm Þii ¼ aii
ðr;mÞ
¼ aii
;
1 6 i 6 n:
ð22Þ
(2) The sub matrices of the (m + 1)th row to the nth row, the (m + 1)th column to the nth column of ALm and ARm are equal, i.e. Lm4 ¼ Rm4 :
ð23Þ
That is to say ðl;mÞ
ðALm Þij ¼ ðARm Þij ¼ aij
ðr;mÞ
¼ aij
;
i; j ¼ m þ 1; . . . ; n; 1 6 m 6 n 1:
ð24Þ
(3) In particular, ðALm Þij ¼ 0
ðj ¼ 1; 2; . . . ; m; i ¼ j þ 1; j þ 2; . . . ; nÞ;
ðARm Þij ¼ 0 ði ¼ 1; 2; . . . ; m; j ¼ i þ 1; i þ 2; . . . ; nÞ: We consider the usual triangular splitting of ALm and ARm ALm ¼ DLm ELm F Lm ; where DLm, ELm, FLm are diagonal, strictly lower triangular and strictly upper triangular parts of ALm respectively, and ARm ¼ DRm ERm F Rm ; where DRm, ERm, FRm are diagonal, strictly lower triangular and strictly upper triangular parts of ARm respectively.
618
Y. Zhang et al. / Applied Mathematics and Computation 188 (2007) 612–633
Below we define the Jacobi type iteration associated with the above splittings B ¼ L þ U; BLm ¼ D1 Lm ðE Lm þ F Lm Þ; BRm ¼
D1 Rm ðE Rm
ð25Þ
þ F Rm Þ;
as well as the splittings that define the Gauss–Seidel type iteration matrices T ¼ ðI LÞ1 U ; T Lm ¼ ðDLm ELm Þ1 F Lm ; 1
T Rm ¼ ðDRm ERm Þ F Rm : In terms of the expressions of ALm (16) and ARm (18), we can express these splitting matrices in the following block form: ! ! ! 0 DL1 0 0 F L1 F L2 DLm ¼ ; ELm ¼ ; F Lm ¼ ; 0 DL4 0 EL4 0 F L4 ! ! ! 0 0 DR1 ER1 0 0 DRm ¼ ; ERm ¼ ; F Rm ¼ : 0 DR4 ER3 ER4 0 F R4 Thus, we have BLm ¼ D1 Lm ðELm þ F Lm Þ !1 DL1 0 F L1 ¼ 0 DL4 0 ¼
EL4 þ F L4 !
D1 L1 F L1
D1 L1 F L2
0
D1 L4 ðE L4 þ F L4 Þ
BRm ¼ D1 Rm ðE Rm þ F Rm Þ !1 DR1 0 ER1 ¼ 0 DR4 ER3 ¼
!
F L2
ð26Þ
;
!
0 ER4 þ F R4 !
D1 R1 ER1
0
D1 R4 ER3
D1 R4 ðE R4 þ F R4 Þ
ð27Þ
;
T Lm ¼ ðDLm ELm Þ1 F Lm ¼
¼
0
0
DL4 EL4
D1 L1
F L1
F L2
0
F L4
!
0 ðDL4 EL4 Þ
0 ¼
!1
DL1
1
F L1
F L2
0
F L4
D1 L1 F L1
D1 L1 F L2
0
ðDL4 EL4 Þ F L4
1
! !
! ;
ð28Þ
Y. Zhang et al. / Applied Mathematics and Computation 188 (2007) 612–633
619
and 1
T Rm ¼ ðDRm ERm Þ F Rm 1 DR1 ER1 0 0 ¼ DR4 ER4 ER3 0 ðDR1 ER1 Þ
¼ ¼
F R4
1
1
0
ðDR4 ER4 Þ ER3 ðDR1 ER1 Þ 0 0 : 1 0 ðDR4 ER4 Þ F R4
!
0 1
ðDR4 ER4 Þ
1
0
0
0 F R4 ð29Þ
Then we have Lemma 1. If A is an M-matrix, then for any 1 6 m 6 n 1 (1) ALm exists; (2) Mm P 0 and Pm P 0; (3) ALm and Lm4 both are M-matrices. Proof. For the case m = 1: Since A is an M-matrix, a11 > 0. Then there exist M1 and P1 such that M 1 P 0;
P 1 P 0:
Hence, by AL1 = P1A, AL1 exists. And for a Z-matrix A, ‘‘A is an M-matrix’’ is equivalent to ‘‘there exists a positive vector y (>0) 2 Rn such that Ay > 0’’ (see [6, p. 136]. P1 P 0, implies AL1 y ¼ P 1 Ay > 0: Consequently, AL1, which is a Z-matrix, is an M-matrix. Since for a Z-matrix the statement ‘‘A is an M-matrix’’ is equivalent to ‘‘all the principal minors of A are positive’’ [6, p. 134], the submatrix L14 is also an M-matrix. Similarly, we could use the same way for the case m = 1 step by step to prove the results for 2 6 m 6 n 1. The proof is completed. h Lemma 2. If A is a nonsingular M-matrix, then, for any 1 6 m 6 n 1, (1) ARm exists; (2) M m P 0 and Qm P 0; (3) ARm and Rm4 both are nonsingular M-matrices. Proof. Since A is an M-matrix, AT is also an M-matrix. For the case m = 1: For A being an M-matrix, so a11 > 0. Then there exist M 1 and Q1 such that M 1 P 0;
Q1 P 0:
Hence, for AR1 = A Q1, AR1 exists. And for a Z-matrix A, ‘‘A is an M-matrix’’ is equivalent to ‘‘there exists a positive vector y(>0) 2 Rn such that Ay > 0’’ (see [6, p. 136]). So there exists a positive vector y(>0) 2 Rn such that ATy > 0. Q1 P 0, implies T
ATR1 y ¼ ðAQ1 Þ y ¼ QT1 AT y > 0: Consequently, ATR1 , which is a Z-matrix, is an M-matrix. Then AR1 is an M-matrix. Since for a Z-matrix the statement ‘‘A is an M-matrix’’ is equivalent to the statement ‘‘all the principal minors of A are positive’’ [6, p. 134], then the submatrix R14 is also an M-matrix. Similarly, we could use the same way for the case m = 1 step by step to prove the results for 2 6 m 6 n 1. The proof is completed. h
620
Y. Zhang et al. / Applied Mathematics and Computation 188 (2007) 612–633
Lemma 3. If A is an M-matrix, then Lm3 ¼ 0; Rm2 ¼ 0; diagðLm1 Þ ¼ diagðRm1 Þ; Lm4 ¼ Rm4 :
ð30Þ ð31Þ ð32Þ ð33Þ
Proof. Since A is an M-matrix, from Lemmas 1 and 2, there exist ALm and ARm for any 1 6 m 6 n 1. From (17) and (19), we immediately get Lm3 = 0 and Rm2 = 0. And from (21), (22) and (24), we have diagðLm1 Þ ¼ diagðRm1 Þ;
ð34Þ
Lm4 ¼ Rm4 :
ð35Þ
Lemma 4. If A is an M-matrix, then DL1 ¼ DR1 ; DL4 ¼ DR4 ; EL4 ¼ ER4 ; F L4 ¼ F R4 :
ð36Þ ð37Þ ð38Þ ð39Þ
Proof. Since A is an M-matrix, from Lemmas 1 and 2, there exist ALm and ARm for any 1 6 m 6 n 1. From the Eq. (21), the result follows. h Lemma 5. The existence of Gauss preconditioners Pm and Qm is guaranteed by the same condition (1) a11 5 0 for m = 1; ðl;m2Þ ðl;m2Þ ðl;m2Þ ðr;m2Þ ðr;m2Þ ðr;m2Þ ðr;m2Þ (2) aðl;m2Þ am1;m1 6¼ am;m1 am1;m or am;m am1;m1 6¼ am;m1 am1;m for 2 6 m 6 n 1. m;m Proof. Notice that, for any 1 6 m 6 n 1, two kinds of Gauss preconditioners Pm and Qm with order m exists if and only if the fore mentioned preconditioned matrices AL,m1 and AR,m1 have the nonzero diagonal elements i.e., ðl;m1Þ
ðr;m1Þ
am1;m1 ¼ am1;m1 6¼ 0: For the case m = 1, it only need a11 5 0 which keep the existence of P1 and Q1. For the case m = 2, it only needs M1A and AM 1 having nonzero diagonal elements, that is to say a22
a21 a12 6¼ 0; a11
i.e. a22a11 5 a21a12. Totally, for any 2 6 m 6 n 1, left Gauss preconditioners Pm with order m exist if and only if ðl;m2Þ ðl;m2Þ
¼ aðl;m2Þ aðl;m1Þ m;m m;m
am;m1 am1;m ðl;m2Þ
am1;m1
ðl;m2Þ
6¼ 0;
ð40Þ
ðl;m2Þ ðl;m2Þ
which implies aðl;m2Þ am1;m1 6¼ am;m1 am1;m . m;m For the same m as Pm, right Gauss preconditioners Qm with order m exists if and only if ðr;m2Þ ðr;m2Þ
ðr;m1Þ ðr;m2Þ ¼ am;m am;m
am;m1 am1;m ðr;m2Þ
am1;m1
ðr;m2Þ
6¼ 0;
ðr;m2Þ ðr;m2Þ
ðr;m2Þ which implies am;m am1;m1 6¼ am;m1 am1;m .
ð41Þ
Y. Zhang et al. / Applied Mathematics and Computation 188 (2007) 612–633
621
In terms of (21), we obtain ðr;m2Þ aðl;m2Þ ¼ am;m ; m;m ðl;m2Þ
ðr;m2Þ
am1;m1 ¼ am1;m1 ; ðl;m2Þ
ðr;m2Þ
ðl;m2Þ
ðr;m2Þ
ð42Þ
am;m1 ¼ am;m1 ; am1;m ¼ am1;m : Thus, by expression (40), expression (41) and (42), the proof is completed. h For convenience, the following notations are needed: • B(M): The Jacobi iteration matrix of the matrix M. • T(M): The Gauss–Seidel iteration matrix of the matrix M. • BLm4, BRm4: the Jacobi iteration matrices of Lm4, Rm4, respectively, where, BLm4 ¼ D1 L4 ðE L4 þ F L4 Þ;
ð43Þ
D1 R4 ðER4
ð44Þ
BRm4 ¼
þ F R4 Þ:
• TLm4, TRm4: the Gauss–Seidel iterative matices of Lm4, Rm4, respectively, where, T Lm4 ¼ ðDL4 EL4 Þ1 F L4 ; 1
T Rm4 ¼ ðDR4 ER4 Þ F R4 :
ð45Þ ð46Þ
Lemma 6. If A is an M-matrix, then (1) BLm4 = BRm4; (2) TLm4 = TRm4. Proof. Since A is an M-matrix, from Lemmas 1 and 2, there exist ALm and ARm for any 1 6 m 6 n 1. ALm 1 and ARm are M-matrices with nonzero diagonal elements, then diagonal matrices D1 L4 and DR4 exist. From Lemma 4, the Eqs. (43)–(46), the results are followed. h Now, we will discuss the properties of Gauss type preconditioners associated with Jacobi, Gauss–Seidel iterative methods, mainly including: (1) The preconditioned matrices of same order left and right Gauss preconditioning methods having the same spectral radii, i.e. q(ARm) = q(ALm). (2) The Jacobi iteration matrices of preconditioned matrices of same order left and right Gauss preconditioning methods having the same spectral radii, i.e. q(BLm) = q(BRm). (3) The Gauss–Seidel iteration matrices of preconditioned matrices of same order left and right Gauss preconditioning methods having the same spectral radii, i.e. q(TLm) = q(TRm). (4) The Jacobi iteration matrices of left (right) Gauss type preconditioned system and its submatrix Lm4 (Rm4) having the same spectral radii, i.e. q(BLm) = q(BLm4) (q(BRm) = q(BRm4)). (5) The Gauss–Seidel iteration matrices of left (right) Gauss type preconditioned system and its submatrix Lm4 (Rm4) having the same spectral radii, i.e. q(TLm) = q(TLm4) (q(TRm) = q(TRm4)). (6) Spectral radii properties on the Jacobi, Gauss–Seidel iteration matrices applied with different order left (right) Gauss preconditioners. Theorem 1. If A is an M-matrix, then qðALm Þ ¼ qðARm Þ;
i:e:; qðP m AÞ ¼ qðAQm Þ:
ð47Þ
622
Y. Zhang et al. / Applied Mathematics and Computation 188 (2007) 612–633
Proof. For A being an M-matrix, from Lemmas 1 and 2, there exist ALm and ARm for any 1 6 m 6 n 1. !! F L2 DL1 F L1 detðkI ALm Þ ¼ det kI 0 DL4 EL4 F L4 !! DL1 0 ¼ det kI 0 DL4 EL4 F L4 !! DL1 0 ¼ det kI ð48Þ 0 Lm4 and DR1 EL1 0 detðkI ARm Þ ¼ det kI ER3 DR4 ER4 F R4 DR1 0 ¼ det kI 0 DR4 ER4 F R4 DR1 0 ¼ det kI : 0 Rm4
ð49Þ
From Lemmas 3 and 4, we have that DL1 = DR1 and Lm4 = Rm4. Then the two equations det(kI ALm) = 0 and det(kI ARm) = 0 have the same characteristic roots which implies qðARm Þ ¼ qðALm Þ; i.e. q(PmA) = q(AQm).
h
Theorem 2. If A is an M-matrix, then qðBLm Þ ¼ qðBRm Þ;
ð50Þ
qðT Rm Þ ¼ qðT Lm Þ:
ð51Þ
Proof. For A being an M-matrix, from Lemmas 1 and 2, there exist ALm and ARm for any 1 6 m 6 n 1. ALm 1 1 and ARm are M-matrices with nonzero diagonal elements, then the four diagonal matrices D1 L1 , DL4 , DR1 and 1 DR4 exist. (1) detðkI BLm Þ ¼ det kI
D1 L1 F L1 0
D1 L1 F L2 1 DL4 ðEL4 þ F L4 Þ
D1 R1 ER1 D1 R4 ER3
0
Similarly, we can get detðkI BRm Þ ¼ det kI
D1 R4 ðE R4 þ F R4 Þ
!!
!!
0 0 : ¼ det kI 0 D1 L4 ðE L4 þ F L4 Þ
ð52Þ
0 ¼ det kI 0
ð53Þ
0 D1 R4 ðE R4 þ F R4 Þ
:
From Lemma 4, it gives DL4 = DR4, EL4 = ER4 and FL4 = FR4. Then detðkI BLm Þ ¼ detðkI BRm Þ: Hence, the two equations det(kI BLm) = 0 and det(kI BRm) = 0 have the same characteristic roots. So qðBRm Þ ¼ qðBLm Þ: (2) Since the triangular matrices DL4 EL4 and DR4 ER4 have nonzero diagonal elements, (DL4 EL4)1 and (DR4 ER4)1 exist. From Lemma 4, we have
Y. Zhang et al. / Applied Mathematics and Computation 188 (2007) 612–633 1
1
ðDR4 ER4 Þ F R4 ¼ ðDL4 EL4 Þ F L4 :
623
ð54Þ
Defining 1
1
T 4 ¼ ðDR4 ER4 Þ F R4 ¼ ðDL4 EL4 Þ F L4 ; then we can write TLm and TRm as ! D1 D1 L1 F L1 L1 F L2 T Lm ¼ ; 0 T4 0 0 T Rm ¼ 1 0 ðDR4 ER4 Þ F R4 0 0 ¼ : 0 T4
ð55Þ
ð56Þ
1 Since diagðD1 L1 F L1 Þ ¼ 0 and DL1 F L1 is a strictly upper triangular matrix, we have 0 0 detðkI T Lm Þ ¼ det kI ¼ detðkI T Rm Þ: 0 T4
Then the two equations det(kI TLm) = 0 and det(kI TRm) = 0 have the same roots. So qðT Lm Þ ¼ qðT Rm Þ:
Remark 1. By Theorem 2, we know the same order left and right Gauss preconditioners have the same convergence rate in Jacobi and Gauss–Seidel iterative methods, respectively. Theorem 3. If A is an M-matrix, then qðBLm Þ ¼ qðBLm4 Þ;
ð57Þ
qðT Lm Þ ¼ qðT Lm4 Þ:
ð58Þ
Proof. (1) We know the Jacobi iteration matrices of ALm and Lm4 are BLm and BLm4. For A being an M-matrix, from Lemmas 1 and 2, there exist ALm and ARm for any 1 6 m 6 n 1. ALm and ARm are M-matri1 ces with nonzero diagonal elements, then D1 Lm , DL4 exist. Hence, BLm and BLm4 exists. By (52), we have 0 0 detðkI BLm Þ ¼ det kI ¼ kk detðkI nk D1 L4 ðEL4 þ F L4 ÞÞ ðE þ F Þ 0 D1 L4 L4 L4 k ¼ kk detðkI nk D1 L4 ðEL4 þ F L4 ÞÞ ¼ k detðkI nk BLm4 Þ:
Then det(kI BLm) = 0 and det(kInk BLm4) = 0 have the same nonzero characteristic roots. So q(BLm) = q(BLm4). (2) For the matrix DL4 EL4 being a lower triangular matrix and the diagonal elements are all nonzero, then (DL4 EL4)1 exists. By (55), we have !! D1 D1 L1 F L1 L1 F L2 detðkI T Lm Þ ¼ det kI : 1 0 ðDL4 EL4 Þ F L4 Since the matrix D1 L1 F L1 is a strictly upper triangular matrix, so we have 0 0 ¼ kk detðkI nk ðDL4 EL4 Þ1 F L4 Þ ¼ kk detðkI nk T Lm4 Þ; detðkI T Lm Þ ¼ det kI 1 0 ðDL4 EL4 Þ F L4 ð59Þ where Ink is the identity matrix with order n k.
624
Y. Zhang et al. / Applied Mathematics and Computation 188 (2007) 612–633
Then we conclude that det(kI TLm) = 0 and det(kInk TLm4) = 0 have the same nonzero eigenvalues. Then q(TLm) = q(TLm4). The proof is completed. h Theorem 4. If A is an M-matrix, then qðBRm Þ ¼ qðBRm4 Þ;
ð60Þ
qðT Rm Þ ¼ qðT Rm4 Þ:
ð61Þ
Proof. The result (60) is an immediate consequence of Eqs. (50), (57) and Lemma 6. The result (61) is an immediate consequence of Eqs. (51), (58) and Lemma 6.
h
Theorem 5. If A is an M-matrix, then for any 1 6 j 6 i 6 n 1, 0 ¼ qðBL;n1 Þ 6 qðBLi Þ 6 qðBLj Þ 6 qðBL1 Þ 6 qðBÞ 6 1;
ð62Þ
0 ¼ qðBR;n1 Þ 6 qðBRi Þ 6 qðBRj Þ 6 qðBR1 Þ 6 qðBÞ 6 1;
ð63Þ
0 ¼ qðT L;n1 Þ 6 qðT Li Þ 6 qðT Lj Þ 6 qðT L1 Þ 6 qðT Þ 6 1;
ð64Þ
0 ¼ qðT R;n1 Þ 6 qðT Ri Þ 6 qðT Rj Þ 6 qðT R1 Þ 6 qðT Þ 6 1:
ð65Þ
Proof. For A being an M-matrix, from Lemmas 1 and 2, there exist ALm and ARm for any 1 6 m 6 n 1 and ALm, ARm, Lm4 and Rm4 are M-matrices. For A being an M-matrix, from [1, Theorem 2.1], we have qðBÞ 6 1;
qðT Þ 6 1:
For (62): For m = 1, by [1, Theorem 2.1], we have qðBLm Þ ¼ qðBL1 Þ 6 qðBÞ: For m = 2, we still use the Theorem 2.1 in [1] for the submatrix L14 with order n 1, then qðBðL24 ÞÞ 6 qðBðL14 ÞÞ; i.e. qðBLm4 Þ 6 qðBLðm1Þ4 Þ:
ð66Þ
By Theorem 3, we have qðBL1 Þ ¼ qðBLðm1Þ Þ ¼ qðBLðm1Þ4 Þ
ð67Þ
qðBL2 Þ ¼ qðBLm Þ ¼ qðBLm4 Þ:
ð68Þ
and Then, by (66)–(68), we have q(BL2) 6 q(BL1), i.e., qðBLm Þ 6 qðBLðm1Þ Þ;
m ¼ 2:
For 3 6 m 6 n 1: we use the same method proving the cases m = 2 and m = 1 step by step, then qðBL;n1 Þ 6 qðBLi Þ 6 qðBLj Þ 6 qðBL1 Þ 6 qðBÞ 6 1;
m 6 j 6 i 6 n 1:
In addition, since BL,n1 is an upper triangular matrix, then q(BL,n1) = 0. So the result (62) follows. For (63): By Theorem 2, we know that qðBLm Þ ¼ qðBRm Þ;
m ¼ 1; 2; . . . ; n 1:
Then, by (62), (63) follows. For (64): Since the result of Theorem 2.1 in [1] hold for Gauss–Seidel iterative method, we can use the same proving method used to prove the result (62), to prove the result (64).
Y. Zhang et al. / Applied Mathematics and Computation 188 (2007) 612–633
625
For (65): By Theorem 2, we know that qðT Lm Þ ¼ qðT Rm Þ;
m ¼ 1; 2; . . . ; n 1:
Then by the result (64), the result (65) follows. h It should be mentioned that Theorem 5 theoretically confirm the statistical results in [3]. 3. Preconditioning algorithms In this section we study the construction of the left and right Gauss preconditioners in detail. Matlab-like notation is applied in algorithm descriptions. Firstly, we discuss the construction of two kinds of Gauss transformation matrices which are essential to the corresponding Gauss preconditioners. For brevity, we denote one function by GetML with two input parameters A and k and a return parameter Mk which implement the computation to get the matrix Mk by expression (10). Similarly, the function GetMR with the same input parameters as GetML and a return parameter M k is used to implement the computation to get the matrix M k by (13). So these two procedures used to compute the Gauss transformation matrices are given in Algorithms 1 and 2, respectively. Algorithm 1 (Construction of left Gauss transformation matrices) Input: A, m Output: Mm 1. for k = 1 to m, 2. M = GetML(A, k) 3. A=M*A 4. end for 5. Mm = M Algorithm 2 (Construction of right Gauss transformation matrices) Input: A, m Output: M m 1. for k = 1 to m, 2. M = GetMR(A, k) 3. A=A*M 4. end for 5. M m ¼ M Secondly, the normal algorithms computing Pm and Qm are given in the following. We note that Pm = MmMm1 . . . M2M1 and Qm ¼ M 1 M 2 M m1 M m . Algorithm 3 [3], (Constructing left Gauss type preconditioners) Input: A, m Output: P 1. Initialize: 2. P = I 3. T = A 4. Computation for Pm: 5. for k = 1 to m, 6. Set v be zero row vector. 7. v(k + 1: n) = T(k + 1: n, k)/T(k, k)
626
8. 9. 10. 11. 12. 13. 14. 15. 16.
Y. Zhang et al. / Applied Mathematics and Computation 188 (2007) 612–633
if k = 1 then P(k + 1: n, 1) = v(k + 1: n) else for r = k + 1 to n, P(r, 1: k) = P(r, 1: k) + v(r) * P(k, 1: k) end for end if T=P*A end for.
The algorithm computing the right Gauss type preconditioners is based on expression (12) as follows. Algorithm 4 (Constructing right Gauss type preconditioner) Input: A, m Output: Qm 1. Initialize: 2. P = I 3. T = AT 4. Computation for Pm: 5. for k = 1 to m, 6. Set v be n dimension zero column vector. 7. v(k + 1: n) = T(k + 1: n, k)/T(k, k) 8. if k = 1 then 9. P(k + 1: n, 1) = v(k + 1: n) 10. else 11. for r = k + 1 to n, 12. P(r, 1: k) = P(r,1: k) + v(r) * P(k, 1: k) 13. end for 14. end if 15. T = P * AT 16. end for. 17. Qm = PT We note that in general these two kinds of Gauss preconditioners Pm and Qm need not be constructed if we could get the left preconditioned matrices on ALm = PmA, bLm = Pmb for left Gauss preconditioning techniques and the right preconditioned matrices ARm = AQm and Qm. After the left Gauss preconditioning it remains to solve the preconditioned system ALmx = bLm. In addition, after the right Gauss preconditioning, it remains to solve the preconditioned system ARmy = b and x = Qmy. Thirdly, therefore, we could give the improved Gauss type preconditioning algorithms by expression (14) and expression (15) to the preconditioned linear systems. Algorithm 5 (Inbuilt left Gauss type preconditioning techniques) Input: A, b, m Output: ALm, bL, Pm 1. Set n be the order of A 2. P = I 3. for k = 1 to m 4. Set v be n dimension zero column vector. 5. v(k + 1: n) = A(k + 1: n, k)/A(k, k) 6. for r = k + 1 to n 7. A(r, k + 1: n) = A(r, k + 1: n) + v(r) * A(k, k + 1: n)
Y. Zhang et al. / Applied Mathematics and Computation 188 (2007) 612–633
8. 9. 10. 11. 12. 13. 14. 15. 16. 17.
627
b(r) = b(r) + v(r) * b(k) end for A(k + 1: n, k) = 0 for i = 1 to k P(k + 1: n, i) = P(k + 1: n, i) + P(k, i) * v(k + 1: n) end for end for. ALm = A bLm = b Pm = P
The lines 2, 11 to 13, 17 in Algorithm 5 are only used to show the computation procedure of the left preconditioners Pm. We should delete these lines in the actual application. Note that, for m = n 1, Algorithm 5 constructing the left preconditioner of order n 1 become a procedure of Gauss elimination. Algorithm 6 (Inbuilt right Gauss type preconditioning techniques) Input: A, m Output: ARm, Qm 1. Set n be the order of A. 2. Q = I 3. for k = 1 to m 4. Set v be n dimension zero row vector. 5. v(k + 1: n) = A(k, k + 1: n)/A(k, k) 6. for r = k + 1 to n 7. A(k + 1: n, r) = A(k + 1: n, r) + v(r) * A(k + 1: n, k) 8. end for 9. A(k, k + 1: n) = 0 10. for i = 1 to k 11. Q(i, k + 1: n) = Q(i, k + 1: n) + Q(i, k) * v(k + 1: n) 12. end for 13. end for. 14. ARm = A 15. Qm = Q Finally, we study floating point operations of Algorithms 5 and 6. Proposition 1. The number of floating point operations of Algorithm 5 (m order left Gauss preconditioning technique) is 1 mð4m2 mð3 þ 12nÞ þ 12n2 þ 6n 7Þ: 6
ð69Þ
Proof. For computing Gauss vector v, the number of floating point operations in line 5, denoted by Nv, is N v ¼ n k: For computing preconditioned matrix ALm, the number of floating point operations in lines 6–9 except line 8, denoted by NA, is n X 2 2ðn kÞ ¼ 2ðn kÞ : NA ¼ r¼kþ1
For the computing preconditioned right hand vector b, the number of floating point operations in lines 6–9 except line 7, denoted by Nb, is
628
Y. Zhang et al. / Applied Mathematics and Computation 188 (2007) 612–633
N b ¼ 2ðn kÞ: Now the total cost of the kth step of the main loop is N k ¼ N v þ N A þ N b:
ð70Þ
Summing Eq. (70) over k = 1 ! m yields the total number of floating point operation of Algorithm 5 as m m m X X X Nk ¼ Nv þ NA þ Nb ¼ ½ðn kÞ þ 2ðn kÞ2 þ 2ðn kÞ k¼1
k¼1
k¼1
¼ mð2n2 þ 3nÞ ð4n þ 3Þ
mðm þ 1Þ 2mðm þ 1Þð2m þ 1Þ þ 2 6
1 ¼ mð4m2 mð3 þ 12nÞ þ 12n2 þ 6n 7Þ: 6
ð71Þ
Proposition 2. The number of floating point operations of Algorithm 5 (left Gauss preconditioning technique) with m = n 1 is 2 3 1 2 7 n þ n n: 3 2 6
ð72Þ
Proof. Substituting m = n 1 in (69) gives n1 X k¼1
2 1 7 N k ¼ n3 þ n2 n: 3 2 6
ð73Þ
Proposition 3. The number of floating point operations of Algorithm 6 (m order right Gauss preconditioning technique) is 2n2 m nm2
m m2 : 2 2
ð74Þ
Proof. For computing Gauss vector v, the number of floating point operations in line 5, denoted by Nv, is N v ¼ n k: For computing preconditioned matrix ARm, the number of floating point operations in lines 6–8, denoted by NA, is n X 2 2ðn kÞ ¼ 2ðn kÞ : NA ¼ r¼kþ1
For computing preconditioner Q, the number of floating point operations in lines 10–12, denoted by NQ, is N Q ¼ 2kðn kÞ: Now the total cost of the kth step of the main loop is N k ¼ N v þ N A þ N Q:
ð75Þ
Summing Eq. (75) over k = 1 ! m yields the total number of floating point operation of Algorithm 5 as m m m X X X 2 Nk ¼ Nv þ NA þ NQ ¼ ½ðn kÞ þ 2ðn kÞ þ 2kðn kÞ k¼1
k¼1
k¼1
m X mðm þ 1Þ ½ðn þ 2n2 Þ ð1 þ 2nÞk ¼ nð1 þ 2nÞm ð1 þ 2nÞ ¼ 2 k¼1 2 mþ1 m m ¼ mð1 þ 2nÞ n ¼ 2n2 m nm2 : 2 2 2
ð76Þ
Y. Zhang et al. / Applied Mathematics and Computation 188 (2007) 612–633
629
Proposition 4. The number of floating point operations of Algorithm 6 (right Gauss preconditioning technique) with m = n 1 is n3 12 n2 12 n. Proof. Substituting m = n 1 in (74) gives n1 X k¼1
Nk ¼
nð1 þ 2nÞðn 1Þ n2 n ¼ n3 : 2 2 2
Proposition 5. The number of floating point operations of Algorithm 5 (left Gauss preconditioning technique) is always less than that of Algorithm 6 (right Gauss preconditioning technique) for any m 6 n 1 when n P 2. Proof. Subtracting (69) by (74), for any m 6 n 1, then the result, denoted by Nd, yields 1 1 1 N d ¼ mð5 3m þ 2m2 þ 6n 6mnÞ ¼ mð2mðm nÞ þ ð6n 3m 5 4mnÞÞ < mð6n 4mnÞ: 6 6 6 ð77Þ For the case m = 1, we have Nd = 1 < 0. For the cases m P 2, we have Nd < 0 immediately from (77). So, the result is proved. h Remark 2. Comparing Propositions 2 and 4 shows that Algorithm 5 costs n3 ðn2 3n þ 2Þ less floating point operations. Therefore, we conclude that the left Gauss preconditioning techniques are superior to the right Gauss preconditioning techniques concerning the proposed algorithms. Numerical results of comparisons on the flops of Algorithms 5 and 6 are given in Figs. 1 and 2. 4. Numerical experiments Numerical experiments are presented to see the performance of these inbuilt Gauss type preconditioning algorithms (Algorithms 5 and 6) associated with the Jacobi, Gauss–Seidel iterative methods. In all the numerical methods we tested, the initial approximation x(0) is taken as the zero-vector, and the right hand side vector b is chosen so that x ¼ ð1; 1; . . . ; 1ÞT is the solution vector for the considered linear system. We have used kb Ax(k)k/k bk 6 106 as the stopping criteria. The max iteration in iterative methods is 8
10
x 10
Left Gauss Right Gauss 9
The number of floating point operations
8
7
6
5
4
3
2
1
0
0
100
200
300
400
500 m (n=1000)
600
700
800
900
1000
Fig. 1. Comparison on flops of Algorithms 5 and 6 on the case n = 1000.
630
Y. Zhang et al. / Applied Mathematics and Computation 188 (2007) 612–633 11
10
x 10
Left Gauss Right Gauss 9
The number of floating point operations
8
7
6
5
4
3
2
1
0
0
1000
2000
3000
4000
5000 6000 m (n=10000)
7000
8000
9000
10000
Fig. 2. Comparison on flops of Algorithms 5 and 6 on the case n = 10,000.
limited to 500. All experiments were executed on a PC using the Matlab programming package. Matlab carries out all calculations in double precision by default. In all tables, we report the spectral radius of the corresponding iteration matrix (Radius), the time in seconds for computing the preconditioner (P-time), the time for the iterative part only (It-time), and the number of iterations (Its). Left (1), Right (1), Left (2) and Right (2) are used to denote Algorithms 3–6 respectively. Example 1. Consider the two-dimensional elliptic partial differential equation (Helmholtz) [9]
o2 u o2 u þ uðx; yÞ ¼ 0; ox2 oy 2
on the rectangular region R defined by the inequalities 0 6 x 6 2 and 0 6 y 6 1 with the Dirichlet boundary conditions. We have used two uniform meshes of Dx = 2/21, Dy = 1/21 which leads to a matrix of order n = 20 · 20 where Dx and Dy refer to the mesh sizes in the x-direction and y-direction, respectively. Table 1 Comparisons of left Gauss preconditioning methods with Jacobi iterative method for Example 1 Preconditioner
Classical GS(1) GS(10) GS(50) GS(90) GS(130) GS(170) GS(210) GS(250) GS(290) GS(330) GS(370) GS(399)
Radius
0.987935 0.987935 0.987921 0.987719 0.987074 0.985779 0.983534 0.979719 0.972898 0.959152 0.924344 0.779000 0.302955
Its
815 815 814 800 762 698 611 505 388 266 150 64 43
‘‘Classical’’: the unpreconditioned iterative method.
Left (1)
Left (2)
P-time
It-time
P-time
It-time
0.13 0.81 3.84 6.95 10.03 13.14 16.22 19.31 22.34 25.34 28.27 30.36
7.02 7.34 7.31 7.13 6.83 6.30 5.47 4.64 3.61 2.42 1.41 0.64 0.45
0.02 0.20 0.92 1.45 1.89 2.30 2.53 2.73 2.89 2.92 2.94 3.00
7.02 6.87 6.78 6.73 6.34 5.89 5.08 4.25 3.27 2.22 1.33 0.55 0.36
Y. Zhang et al. / Applied Mathematics and Computation 188 (2007) 612–633
631
Note that in Tables 1–4, for the highest order left and right Gauss preconditioned linear system, here we use the classical Jacobi and Gauss–Seidel to solve the triangular linear systems. Table 2 Comparisons of right Gauss preconditioning methods with Jacobi iterative method for Example 1 Preconditioner
Radius
Its
Classical GS(1) GS(10) GS(50) GS(90) GS(130) GS(170) GS(210) GS(250) GS(290) GS(330) GS(370) GS(399)
0.987935 0.987935 0.987921 0.987719 0.987074 0.985779 0.983534 0.979719 0.972898 0.959152 0.924344 0.779000 0.302955
815 815 814 802 764 700 614 509 391 270 154 64 43
Left (1)
Left (2)
P-time
It-time
P-time
It-time
0.14 0.83 3.88 6.95 10.05 13.17 16.25 19.33 22.38 25.38 28.30 30.39
7.22 7.00 7.02 6.88 6.50 6.05 5.25 4.39 3.41 2.34 1.38 0.61 0.44
0.02 0.11 0.55 0.98 1.42 1.83 2.17 2.61 2.88 3.19 3.41 3.53
7.22 6.94 6.88 6.81 6.44 5.94 5.20 4.33 3.45 2.28 1.31 0.55 0.38
Table 3 Comparisons of left Gauss preconditioning methods with Gauss–Seidel iterative method for Example 1 Preconditioner
Radius
Its
Classical GS(1) GS(10) GS(50) GS(90) GS(130) GS(170) GS(210) GS(250) GS(290) GS(330) GS(370) GS(399)
0.976015 0.976015 0.975988 0.975590 0.974315 0.971762 0.967342 0.959858 0.946551 0.920040 0.854735 0.611080 0.297231
437 437 437 430 410 376 329 273 210 144 87 49 43
Left (1)
Left (2)
P-time
It-time
P-time
It-time
0.14 0.81 3.88 6.94 10.95 13.19 16.27 19.34 22.39 25.33 28.34 30.41
7.06 6.81 6.86 6.55 6.27 6.09 5.02 4.19 3.28 2.25 1.38 0.83 0.73
0.03 0.20 0.89 1.52 2.02 2.30 2.52 2.73 2.89 2.94 2.95 3.00
7.06 6.66 6.69 6.53 7.06 7.45 5.03 4.17 3.22 2.16 1.30 0.75 0.66
Table 4 Comparisons of right Gauss preconditioning methods with Gauss–Seidel iterative method for Example 1 Preconditioner
Classical GS(1) GS(10) GS(50) GS(90) GS(130) GS(170) GS(210) GS(250) GS(290) GS(330) GS(370) GS(399)
Radius
0.976015 0.976015 0.975988 0.975590 0.974315 0.971762 0.967342 0.959858 0.946551 0.920040 0.854735 0.611080 0.000000
Its
437 437 437 430 410 376 329 273 210 144 82 30 2
Left (1)
Left (2)
P-time
It-time
P-time
It-time
0.14 0.81 3.88 6.95 10.06 13.20 16.30 19.38 22.44 25.47 28.73 30.67
6.92 6.72 6.80 6.69 6.39 5.88 5.13 4.27 3.30 2.27 1.33 0.53 0.09
0.02 0.11 0.56 1.00 1.45 1.84 2.19 2.61 2.89 3.19 3.38 3.52
6.92 6.72 6.73 6.61 6.30 5.83 5.06 4.19 3.22 2.20 1.27 0.47 0.03
632
Y. Zhang et al. / Applied Mathematics and Computation 188 (2007) 612–633
From Tables 1–4, we find that the improved inbuilt Gauss preconditioning algorithms (Algorithms 5 and 6) are greatly superior to the non-inbuilt ones (Algorithms 2 and 4). Example 2. Problems set called LANPRO (Lanczos with partial reorthogonalization) from the Harwell-Boeing Collection. Problems matrices Nos1, Nos2 and Nos5 all with nonzero diagonal elements are involved in the experiments. Note that the highest order preconditioned systems are solved by Gauss–Seidel iterative method and direct triangular solving method (Marked with ‘‘Tri’’), respectively. Numerical results are given in Tables 5–7.
Table 5 Comparisons of left Gauss Preconditioning methods with Gauss–Seidel iterative method (Example 2. Nos1 with order n = 237) Preconditioner
Radius
Its
Left (1)
Left (2)
P-time
It-time
Min
Max
Classical GS(1) GS(10) GS(57) GS(104) GS(151) GS(198) GS(236) GS(236)Tri
0.999999 0.999999 0.999999 0.999999 0.999998 0.999996 0.999945 0.000000 0.000000
500 500 500 500 500 500 500 91 –
– 0.02 0.13 0.41 0.74 0.69 0.84 0.97 0.97
3.92 3.62 3.66 3.70 4.02 4.00 3.86 0.70 0.02
0.240140 0.240140 0.124329 0.010632 0.010627 0.010178 0.002500 1.000000 1.000000
2.242217 2.242217 2.242217 2.242217 2.242220 2.242260 2.244213 1.000000 1.000000
‘‘Min’’ : the min entry of the final numerical solution vector. ‘‘Max’’ : the max entry of the final numerical solution vector.
Table 6 Comparisons of Left Gauss preconditioning methods with Gauss–Seidel iterative method (Example 2. Nos5 with order n = 468) Preconditioner
Radius
Its
Left (1)
Left (2)
P-time
It-time
Min
Max
Classical GS(1) GS(10) GS(103) GS(196) GS(289) GS(382) GS(467) GS(467)Tri
0.997904 0.997904 0.997903 0.997893 0.997851 0.997744 0.996816 0.000000 0.000000
500 500 500 500 500 500 500 33 –
– 0.09 0.34 2.39 4.11 4.45 4.78 4.84 4.84
10.55 10.19 10.20 10.23 10.20 10.16 10.25 0.67 0.02
0.417626 0.417129 0.416974 0.413352 0.386359 0.440179 0.650278 1.000000 1.000000
1.008461 1.008470 1.008475 1.008573 1.017740 1.016211 1.053018 1.000000 1.000000
Table 7 Comparisons of left Gauss Preconditioning methods with Gauss–Seidel iterative method (Example 2. Nos2 with order n = 957) Preconditioner
Radius
Its
Classical GS(1) GS(10) GS(201) GS(392) GS(583) GS(774) GS(956) GS(956)Tri
1.000000 1.000000 1.000000 1.000000 1.000000 1.000000 1.000000 0.000000 0.000000
500 500 500 500 500 500 500 334 –
Left (1)
Left (2)
P-time
It-time
Min
Max
– 0.19 1.22 17.02 29.45 34.28 35.81 36.33 36.33
32.88 31.64 31.80 34.39 32.25 31.88 31.44 22.69 0.03
4.208663 4.208663 3.746223 0.061133 0.061133 0.061133 0.061133 1.000000 1.000000
6.216589 6.216589 6.216589 6.216589 6.216589 6.216589 6.216589 1.000000 1.000000
Y. Zhang et al. / Applied Mathematics and Computation 188 (2007) 612–633
633
Finally, this paper presents some methods and numerical improvements for preconditioning. Parameterized m-order Gauss preconditioning technique is a natural problem. Further research is required. References [1] A. Hadjidimos, D. Noutsos, M. Tzoumas, More on modifications and improvements of classical iterative schemes for M-matrices, Linear Algebra Appl. 364 (2003) 253–279. [2] J.P. Milaszewicz, Improving Jacobi and Gauss–Seidel iterations, Linear Algebra Appl. 93 (1987) 161–170. [3] Y. Zhang, T.Z. Huang, X.P. Liu, Modified iterative methods for nonnegative matrices and M-matrices linear systems, Comput. Math. Appl. 50 (2005) 1587–1602. [4] G.H. Golub, C.F. Van Loan, Matrix Computations, third ed., The Johns Hopkins University Press, 1996. [5] D.J. Evans, M.M. Martins, M.E. Trigo, The AOR iterative method for new preconditioned linear systems, J. Comput. Appl. Math. 132 (2001) 461–466. [6] A. Berman, P.J. Plemmons, Nonnegative Matrices in the Mathematical Sciences, Academic Press, 1979. [7] T. Kohno, H. Kotakemori, H. Niki, M. Usui, Improving the Gauss–Seidel method for Z-matrices, Linear Algebra Appl. 267 (1997) 113–123. [8] A.D. Gunawardena, S.K. Jain, L. Snyder, Modified iterative methods for consistent linear systems, Linear Algebra Appl. 154–156 (1991) 123–143. ¨ zlu¨, On the numerical solution of Helmholtz equation by alternating group explicit (AGE) methods, Appl. Math. [9] Necdet Bildik, Sevil O Comput. 163 (2005) 505–518.