# A Parallel Preconditioned Modified Conjugate Gradient Method for Large Sylvester Matrix Equation

Authors: Junxia Hou; Quanyi Lv; Manyu Xiao
Journal: Mathematical Problems in Engineering
 (2014)

---

## Abstract

Computational effort of solving large-scale Sylvester equationsA
X
+
X
B
+
F
=
O is frequently hindered in dealing with many complex control problems. In this work, a parallel preconditioned algorithm for solving it is proposed based on combination of a parameter iterative preconditioned method and modified form of conjugate gradient (MCG) method. Furthermore, Schur’s inequality and modified conjugate gradient method are employed to overcome the involved difficulties such as determination of parameter and calculation of inverse matrix. Several numerical results finally show that high performance of proposed parallel algorithm is obtained both in convergent rate and in parallel efficiency.


---

## Body

## 1. Introduction

Sylvester equations play vital roles in a number of applications such as control theory [1], matrix eigendecompositions [2], model reduction [3–5], and image processing [6]. Particularly, when solving parabolic and elliptic PDEs by the discretization method, the issue of solving a large scale Sylvester matrix equation is often encountered. Nowadays, there are two types of methods for solving these matrix equations: the first common one is the standard direct method proposed by Bartels and Stewart [7]. However, this is not applicable in large-scale settings. In order to overcome this limitation, the other concerning iterative form is discussed such as a few simple iterative schemes proposed by Ding and Chen [8]. ADI method was first introduced by Peaceman and Rachford [9] for solving parabolic and elliptic PDEs and was later adapted solving the Sylvester equation by Wachspress [10]. Although using ADI method needs to choose optimal parameters, it is too difficult to find them for a large scale problem. Therefore, the parameter iterative method was first proposed in [11], and the available choice of parameter was also given for a large scale problem. Recently, a new effective method based on Krylov subspace [12–14] is discussed. In particular, conjugate gradient (CG) method [15] is its essential one, which has the advantages of small store, small computing, and the fact that it is suitable for parallel. But it is only applied for the case that the coefficient matrix is positive definite symmetric matrices. So we can choose the modified conjugate gradient method (MCG) [16, 17] in other cases. The convergence speed of MCG has been closely related to the singular values of the coefficient matrix, and the more concentrated the singular values are, the better the convergence speed of MCG is. In order to enhance the speed of convergence, the idea of preconditioning [18] is introduced. By this way, it will concentrate the singular values of the coefficient matrix. All of these brought about a need for the development of efficient parallel algorithm as a computational time is prohibitive, which motivates us to propose a parallel preconditioned algorithm for solving large Sylvester matrix equation.The remainder of the paper is organized as follows. In Section2, there is the introduction of MCG algorithm for solving large matrix equation A
X
+
X
B
+
F
=
O, and then our proposed preconditioned MCG algorithm is presented in Section 3, followed by the implementation of parallel algorithm (in Section 4). Convergence Analysis is given in Section 5. Finally, some numerical examples are discussed with the results analysis and comparison, followed by concluding comments.Throughout this paper we use the following notations. LetA and B stand for two n
×
n matrices and A
T, ρ
(
A
), and I represent the transpose of matrix A, spectral radius of matrix, A and the n
×
n identity matrix, respectively. At the same time, we define the inner product as [
A
,
B
]
=
tr
⁡
(
A
T
B
); then the matrix norm of A induced by the inner product is Frobenius norm and denoted by ∥
A
∥; that is, ∥
A
∥
=
[
A
,
A
]
=
tr
⁡
(
A
T
A
).


## 2. Introduction of MCG Algorithm for SolvingA

X
+
X

B
+
F
=
O

Ann
×
n Sylvester matrix equation takes the form
(1)
A
X
+
X
B
+
F
=
O
,
where A
∈
ℝ
n
×
n, B
∈
ℝ
n
×
n, and F
∈
ℝ
n
×
n are given matrices, and X
∈
ℝ
n
×
n is unknown matrix. MCG algorithm [15, 19, 20] can be expressed as follows.Algorithm 1.
(
1
) Choose initial vector X
(
0
)
∈
ℝ
n
×
n; let k
:
=
0, and calculate
(2)
R
(
0
)
=
-
F
-
A
X
(
0
)
-
X
(
0
)
B
,
Q
(
0
)
=
A
T
R
(
0
)
+
R
(
0
)
+
R
(
0
)
B
T
,
Z
(
0
)
=
Q
(
0
)
.
(
2
) If R
(
k
)
=
O, or R
(
k
)
≠
O but Z
(
k
)
=
O stop; otherwise, calculate
(3)
γ
k
=
[
R
(
k
)
,
R
(
k
)
]
[
Z
(
k
)
,
Z
(
k
)
]
,
X
(
k
+
1
)
=
X
(
k
)
+
γ
k
Z
(
k
)
.
(
3
) Then perform
(4)
R
(
k
+
1
)
=
-
F
-
A
X
(
k
+
1
)
-
X
(
k
+
1
)
B
,
Q
(
k
+
1
)
=
A
T
R
(
k
+
1
)
+
R
(
k
+
1
)
+
R
(
k
+
1
)
B
T
.Furthermore, the convergent condition of MCG algorithm referred in [15, 19, 20] is illustrated as Theorem 2.Theorem 2.
Assume thatA
X
+
X
B
+
F
=
O is consistent; then, for any arbitrary initial matrix X
(
0
), a solution of A
X
+
X
B
+
F
=
O can be obtained with finite iteration steps in the absence of round-off error.It is important to accelerate the speed of the convergence. Form the matrix equationA
X
+
X
B
+
F
=
O, we know that the speed of convergence depends on the singular values centralization of the coefficient matrix (
A
⊗
I
+
I
⊗
B
T
). With the point of theorem, it is difficult to calculate its singular values. In order to easily enhance the speed of convergence, we propose a preconditioned MCG algorithm mentioned in the following section.


## 3. Proposed Preconditioned MCG Algorithm

In order to choose preconditioner for MCG, we introduce the parameter iterative method in [11] for solving (1), which is formed as follows:
(5)
X
(
k
+
1
)
=
A
~
X
(
k
)
B
~
+
F
~
,
where
(6)
A
~
=
(
α
A
-
I
)
-
1
(
α
A
+
I
)
=
I
+
2
(
α
A
-
I
)
-
1
,
B
~
=
(
α
B
-
I
)
-
1
(
α
B
+
I
)
=
I
+
2
(
α
B
-
I
)
-
1
,
F
~
=
2
α
(
α
A
-
I
)
-
1
F
(
α
B
-
I
)
-
1
,
and α is a parameter. The speed of convergence depends on the value of ρ
(
A
~
⊗
B
~
)
T, since A
~
=
(
α
A
-
I
)
-
1
(
α
A
+
I
), B
~
=
(
α
B
-
I
)
-
1
(
α
B
+
I
). It is much simple to obtain the best convergence effect after decision of an optimal value of α. However, two inverse matrices of α
A
+
I and α
B
+
I must be solved to get the matrices A
~, B
~. Let Y
1
=
(
α
A
-
I
)
-
1 and Y
2
=
(
α
B
-
I
)
-
1; we can find Y
1 and Y
2 by solving the matrix equations (
α
A
-
I
)
Y
1
=
I and (
α
B
-
I
)
Y
2
=
I.Next subsection will give the methods determining the parameterα and solving the matrices A
~, B
~, and F
~, respectively.


### 3.1. Technique for Getting the Parameterα and the Inverse of Matrix (
α
A
-
I
)

#### 3.1.1. Determination of the Parameterα

If the negative real numbersλ
1
≥
λ
2
≥
⋯
≥
λ
n are the eigenvalues of A and the negative real numbers μ
1
≥
μ
2
≥
⋯
≥
μ
n are the eigenvalues of B, based on the optimal parameter measurement in [11], we have either the optimal parameter α
=
1
/
λ
1
λ
n or α
=
1
/
μ
1
μ
n.On basis ofSchur’s inequality in [21], we determine α approximate value.Theorem 3 (Schur’s inequality).
LetA
=
(
a
r
s
)
∈
ℂ
n
×
n, and let the eigenvalues of A be λ
1
,
λ
2
,
…
,
λ
n; then it holds that:
(7)
∑
r
=
1
n
|
λ
r
|
2
≤
∑
r
,
s
=
1
n
|
a
r
s
|
2
=
∥
A
∥
F
2
.
Ifλ
1
,
λ
2
,
…
,
λ
n satisfy |
λ
1
|
≥
|
λ
2
|
≥
⋯
≥
|
λ
n
|, fromSchur’s inequality, one has
(8)
|
λ
1
|
2
≥
∥
Α
∥
F
2
n
,
|
λ
n
|
2
≤
∥
A
∥
F
2
n
;
then we choose either α
=
n
/
∥
Α
∥
F or α
=
n
/
∥
B
∥
F; it is very easy for us to find α.


#### 3.1.2. Algorithm for Solving Matrix Equationα
(
A
-
I
)
Y
1
=
I

We use the MCG algorithm to solve the matrix equations(
α
A
-
I
)
Y
1
=
I and (
α
B
-
I
)
Y
2
=
I. The specific algorithm for it will be given as follows.Algorithm 4.
(
1
) Give an initial matrix Y
1
(
0
)
∈
ℝ
n
×
n, let k
:
=
0 calculate
(9)
R
(
0
)
=
I
-
(
α
A
-
I
)
Y
1
(
0
)
,
Q
(
0
)
=
(
α
A
-
I
)
T
R
(
0
)
,
Z
(
0
)
=
Q
(
0
)
.
(
2
) If R
(
k
)
=
O, or R
(
k
)
≠
O but Z
(
k
)
=
O stop; otherwise, do the form
(10)
γ
k
=
[
R
(
k
)
,
R
(
k
)
]
[
Z
(
k
)
,
Z
(
k
)
]
,
Y
1
(
k
+
1
)
=
Y
1
(
k
)
+
γ
k
Z
(
k
)
.
(
3
) Then perform
(11)
R
(
k
+
1
)
=
I
-
(
α
A
-
I
)
Y
1
(
k
+
1
)
,
Q
(
k
+
1
)
=
(
α
A
-
I
)
T
R
(
k
+
1
)
,
η
k
=
[
R
(
k
+
1
)
,
R
(
k
+
1
)
]
[
R
(
k
)
,
R
(
k
)
]
,
Z
(
k
+
1
)
=
Q
(
k
+
1
)
+
η
k
Z
(
k
)
.
(
4
) Let k
:
=
k
+
1 and turn back to (
2
).Through Algorithm4, we can get the matrix (
α
A
-
I
)
-
1. In similar manner, the inverse of α
B
-
I can be obtained. Therefore, the matrices A
~, B
~, and F
~ can be obtained.


### 3.2. Algorithm for Solving the New Matrix EquationA
~ 
X

B
~
- X
+ F
~
= O

We still use MCG algorithm for solving the matrix equationA
~
X
B
~
-
X
+
F
~
=
O, which is obtained from the parameter iterative method of the above introduced, for solving (5). Thereby the precondition MCG is developed; it can be expressed as follows.Algorithm 5.
(
1
) Given initial matrix X
(
0
)
∈
ℝ
n
×
n, calculate
(12)
R
(
0
)
=
-
F
~
+
X
(
0
)
-
A
~
X
(
0
)
B
~
,
Q
(
0
)
=
A
~
T
R
(
0
)
B
~
T
-
R
(
0
)
;
Z
(
0
)
=
Q
(
0
)
.
(
2
) For k
=
0,1
,
…, calculateStep  1: If∥
R
(
k
)
∥
<
ε, or ∥
R
(
k
)
∥
>
ε but ∥
Z
(
k
)
∥
<
ε stop; otherwise, turn to Step  2;
Step  2:γ
k
=
[
R
(
k
)
,
R
(
k
)
]
/
[
Z
(
k
)
,
Z
(
k
)
];
Step  3:X
(
k
+
1
)
=
X
(
k
)
+
γ
k
Z
(
k
);
Step  4:R
(
k
+
1
)
=
-
F
~
+
X
(
k
+
1
)
-
A
~
X
(
k
+
1
)
B
~, Q
(
k
+
1
)
=
A
~
T
R
(
k
+
1
)
B
~
T
-
R
(
k
+
1
);
Step  5:η
k
=
[
R
(
k
+
1
)
,
R
(
k
+
1
)
]
/
[
R
(
k
)
,
R
(
k
)
];
Step  6:Z
(
k
+
1
)
=
Q
(
k
+
1
)
+
η
k
Z
(
k
).From every step of the previously mentioned algorithms, we can find that there are only matrix multiplication, addition, and inner product, which shows that the algorithm has very good parallelism. Furthermore, at each iteration, amount of calculation with preconditioned modified conjugate gradient method is the same as that with modified conjugate gradient method. There is only computational cost on obtaining the inverse of two matrices. The numerical results will also confirm its advantage both in convergent rate and parallel efficiency.


## 3.1. Technique for Getting the Parameterα and the Inverse of Matrix (
α
A
-
I
)

### 3.1.1. Determination of the Parameterα

If the negative real numbersλ
1
≥
λ
2
≥
⋯
≥
λ
n are the eigenvalues of A and the negative real numbers μ
1
≥
μ
2
≥
⋯
≥
μ
n are the eigenvalues of B, based on the optimal parameter measurement in [11], we have either the optimal parameter α
=
1
/
λ
1
λ
n or α
=
1
/
μ
1
μ
n.On basis ofSchur’s inequality in [21], we determine α approximate value.Theorem 3 (Schur’s inequality).
LetA
=
(
a
r
s
)
∈
ℂ
n
×
n, and let the eigenvalues of A be λ
1
,
λ
2
,
…
,
λ
n; then it holds that:
(7)
∑
r
=
1
n
|
λ
r
|
2
≤
∑
r
,
s
=
1
n
|
a
r
s
|
2
=
∥
A
∥
F
2
.
Ifλ
1
,
λ
2
,
…
,
λ
n satisfy |
λ
1
|
≥
|
λ
2
|
≥
⋯
≥
|
λ
n
|, fromSchur’s inequality, one has
(8)
|
λ
1
|
2
≥
∥
Α
∥
F
2
n
,
|
λ
n
|
2
≤
∥
A
∥
F
2
n
;
then we choose either α
=
n
/
∥
Α
∥
F or α
=
n
/
∥
B
∥
F; it is very easy for us to find α.


### 3.1.2. Algorithm for Solving Matrix Equationα
(
A
-
I
)
Y
1
=
I

We use the MCG algorithm to solve the matrix equations(
α
A
-
I
)
Y
1
=
I and (
α
B
-
I
)
Y
2
=
I. The specific algorithm for it will be given as follows.Algorithm 4.
(
1
) Give an initial matrix Y
1
(
0
)
∈
ℝ
n
×
n, let k
:
=
0 calculate
(9)
R
(
0
)
=
I
-
(
α
A
-
I
)
Y
1
(
0
)
,
Q
(
0
)
=
(
α
A
-
I
)
T
R
(
0
)
,
Z
(
0
)
=
Q
(
0
)
.
(
2
) If R
(
k
)
=
O, or R
(
k
)
≠
O but Z
(
k
)
=
O stop; otherwise, do the form
(10)
γ
k
=
[
R
(
k
)
,
R
(
k
)
]
[
Z
(
k
)
,
Z
(
k
)
]
,
Y
1
(
k
+
1
)
=
Y
1
(
k
)
+
γ
k
Z
(
k
)
.
(
3
) Then perform
(11)
R
(
k
+
1
)
=
I
-
(
α
A
-
I
)
Y
1
(
k
+
1
)
,
Q
(
k
+
1
)
=
(
α
A
-
I
)
T
R
(
k
+
1
)
,
η
k
=
[
R
(
k
+
1
)
,
R
(
k
+
1
)
]
[
R
(
k
)
,
R
(
k
)
]
,
Z
(
k
+
1
)
=
Q
(
k
+
1
)
+
η
k
Z
(
k
)
.
(
4
) Let k
:
=
k
+
1 and turn back to (
2
).Through Algorithm4, we can get the matrix (
α
A
-
I
)
-
1. In similar manner, the inverse of α
B
-
I can be obtained. Therefore, the matrices A
~, B
~, and F
~ can be obtained.


## 3.1.1. Determination of the Parameterα

If the negative real numbersλ
1
≥
λ
2
≥
⋯
≥
λ
n are the eigenvalues of A and the negative real numbers μ
1
≥
μ
2
≥
⋯
≥
μ
n are the eigenvalues of B, based on the optimal parameter measurement in [11], we have either the optimal parameter α
=
1
/
λ
1
λ
n or α
=
1
/
μ
1
μ
n.On basis ofSchur’s inequality in [21], we determine α approximate value.Theorem 3 (Schur’s inequality).
LetA
=
(
a
r
s
)
∈
ℂ
n
×
n, and let the eigenvalues of A be λ
1
,
λ
2
,
…
,
λ
n; then it holds that:
(7)
∑
r
=
1
n
|
λ
r
|
2
≤
∑
r
,
s
=
1
n
|
a
r
s
|
2
=
∥
A
∥
F
2
.
Ifλ
1
,
λ
2
,
…
,
λ
n satisfy |
λ
1
|
≥
|
λ
2
|
≥
⋯
≥
|
λ
n
|, fromSchur’s inequality, one has
(8)
|
λ
1
|
2
≥
∥
Α
∥
F
2
n
,
|
λ
n
|
2
≤
∥
A
∥
F
2
n
;
then we choose either α
=
n
/
∥
Α
∥
F or α
=
n
/
∥
B
∥
F; it is very easy for us to find α.


## 3.1.2. Algorithm for Solving Matrix Equationα
(
A
-
I
)
Y
1
=
I

We use the MCG algorithm to solve the matrix equations(
α
A
-
I
)
Y
1
=
I and (
α
B
-
I
)
Y
2
=
I. The specific algorithm for it will be given as follows.Algorithm 4.
(
1
) Give an initial matrix Y
1
(
0
)
∈
ℝ
n
×
n, let k
:
=
0 calculate
(9)
R
(
0
)
=
I
-
(
α
A
-
I
)
Y
1
(
0
)
,
Q
(
0
)
=
(
α
A
-
I
)
T
R
(
0
)
,
Z
(
0
)
=
Q
(
0
)
.
(
2
) If R
(
k
)
=
O, or R
(
k
)
≠
O but Z
(
k
)
=
O stop; otherwise, do the form
(10)
γ
k
=
[
R
(
k
)
,
R
(
k
)
]
[
Z
(
k
)
,
Z
(
k
)
]
,
Y
1
(
k
+
1
)
=
Y
1
(
k
)
+
γ
k
Z
(
k
)
.
(
3
) Then perform
(11)
R
(
k
+
1
)
=
I
-
(
α
A
-
I
)
Y
1
(
k
+
1
)
,
Q
(
k
+
1
)
=
(
α
A
-
I
)
T
R
(
k
+
1
)
,
η
k
=
[
R
(
k
+
1
)
,
R
(
k
+
1
)
]
[
R
(
k
)
,
R
(
k
)
]
,
Z
(
k
+
1
)
=
Q
(
k
+
1
)
+
η
k
Z
(
k
)
.
(
4
) Let k
:
=
k
+
1 and turn back to (
2
).Through Algorithm4, we can get the matrix (
α
A
-
I
)
-
1. In similar manner, the inverse of α
B
-
I can be obtained. Therefore, the matrices A
~, B
~, and F
~ can be obtained.


## 3.2. Algorithm for Solving the New Matrix EquationA
~ 
X

B
~
- X
+ F
~
= O

We still use MCG algorithm for solving the matrix equationA
~
X
B
~
-
X
+
F
~
=
O, which is obtained from the parameter iterative method of the above introduced, for solving (5). Thereby the precondition MCG is developed; it can be expressed as follows.Algorithm 5.
(
1
) Given initial matrix X
(
0
)
∈
ℝ
n
×
n, calculate
(12)
R
(
0
)
=
-
F
~
+
X
(
0
)
-
A
~
X
(
0
)
B
~
,
Q
(
0
)
=
A
~
T
R
(
0
)
B
~
T
-
R
(
0
)
;
Z
(
0
)
=
Q
(
0
)
.
(
2
) For k
=
0,1
,
…, calculateStep  1: If∥
R
(
k
)
∥
<
ε, or ∥
R
(
k
)
∥
>
ε but ∥
Z
(
k
)
∥
<
ε stop; otherwise, turn to Step  2;
Step  2:γ
k
=
[
R
(
k
)
,
R
(
k
)
]
/
[
Z
(
k
)
,
Z
(
k
)
];
Step  3:X
(
k
+
1
)
=
X
(
k
)
+
γ
k
Z
(
k
);
Step  4:R
(
k
+
1
)
=
-
F
~
+
X
(
k
+
1
)
-
A
~
X
(
k
+
1
)
B
~, Q
(
k
+
1
)
=
A
~
T
R
(
k
+
1
)
B
~
T
-
R
(
k
+
1
);
Step  5:η
k
=
[
R
(
k
+
1
)
,
R
(
k
+
1
)
]
/
[
R
(
k
)
,
R
(
k
)
];
Step  6:Z
(
k
+
1
)
=
Q
(
k
+
1
)
+
η
k
Z
(
k
).From every step of the previously mentioned algorithms, we can find that there are only matrix multiplication, addition, and inner product, which shows that the algorithm has very good parallelism. Furthermore, at each iteration, amount of calculation with preconditioned modified conjugate gradient method is the same as that with modified conjugate gradient method. There is only computational cost on obtaining the inverse of two matrices. The numerical results will also confirm its advantage both in convergent rate and parallel efficiency.


## 4. Parallel Implementation of Proposed Algorithm

For describing convenience, let the processor number bep, n
=
p
l, where l is integer; p
i
(
i
=
1,2
,
…
,
p
) represents ith processor. And mark
(13)
A
=
(
A
1
T
,
A
2
T
,
…
,
A
p
T
)
T
,
A
~
=
(
A
~
1
T
,
A
~
2
T
,
…
,
A
~
p
T
)
T
,
C
~
=
(
C
~
1
,
C
~
2
,
…
,
C
~
p
)
,
B
=
(
B
1
T
,
B
2
T
,
…
,
B
p
T
)
T
,
B
~
=
(
B
~
1
T
,
B
~
2
T
,
…
,
B
~
p
T
)
T
,
I
=
(
I
1
T
,
I
2
T
,
…
,
I
p
T
)
T
,
(
α
A
-
I
)
-
1
=
(
D
1
T
,
D
2
T
,
…
,
D
p
T
)
T
,
(
α
B
-
I
)
-
1
=
(
E
1
T
,
E
2
T
,
…
,
E
p
T
)
T
,
F
=
(
F
1
T
,
F
2
T
,
…
,
F
p
T
)
T
,
X
(
k
)
=
(
(
X
1
(
k
)
)
T
,
(
X
2
(
k
)
)
T
,
…
,
(
X
p
(
k
)
)
T
)
T
,
R
(
k
)
=
(
(
R
1
(
k
)
)
T
,
(
R
2
(
k
)
)
T
,
…
,
(
R
p
(
k
)
)
T
)
T
,
Z
(
k
)
=
(
(
Z
1
(
k
)
)
T
,
(
Z
2
(
k
)
)
T
,
…
,
(
Z
p
(
k
)
)
T
)
T
,
where A
i, A
~
i, B
i, B
~
i, I
i, D
i, E
i, F
i, X
i
(
k
), R
i
(
k
), and Z
i
(
k
) are l
×
n subblock matrices, and C
~
i is n
×
l subblock matrix (
i
=
1,2
,
…
,
p
). Because of subblock matrices A
~
i and C
~
i of A
~ need be stored on p
i
(
i
=
1,2
,
…
,
p
); then, the data of A
~ will be redistributed during the precondition process. And the details can be seen in precondition process. In addition, since the storage manner directly influences the ways of matrix multiplication in parallel, two parallel computing algorithms of matrix multiplication proposed in [22] are used in the paper:Multiplication 1: matrix multiplication with block row-row matrices stored;
Multiplication 2: matrix multiplication with block row-column matrices stored.Its detailed implementation can be expressed as follows.(1)  Storage Method. Store A
i, B
i, I
i, F
i, X
i
(
0
), α on processor p
i
(
i
=
1,2
,
…
,
p
).(2)  Precondition Process. Solve matrix equations (
α
A
-
I
)
Y
1
=
I, (
α
B
-
I
)
Y
2
=
I by Algorithm 4 in parallel, then we can get the matrices (
α
A
-
I
)
-
1 and (
α
B
-
I
)
-
1. Calculate I
i
+
2
D
i on p
i
(
i
=
1,2
,
…
,
p
), the subblock matrices A
~
i of A
~ can be obtained by p
i
(
i
=
1,2
,
…
,
p
); all the p processors pass p times in turn, and the sub-block matrices C
~
i of A
~ can be obtained by p
i
(
i
=
1,2
,
…
,
p
). Thus, the two storage ways of A
~ are realized. Calculate I
i
+
2
E
i on p
i
(
i
=
1,2
,
…
,
p
); the subblock matrices B
~
i of B
~ can be obtained by p
i
(
i
=
1,2
,
…
,
p
).(3)  Calculate Process. Calculate matrix multiplication A
~
X
(
0
) by Multiplication 1; then calculate matrix multiplication (
A
~
X
(
0
)
)
B
~ by Multiplication 1. Finally, we can obtain R
(
0
); then, get matrix multiplication A
~
T
R
(
0
) by Multiplication 1, and calculate (
A
~
T
R
(
0
)
)
B
~
T by Multiplication 2, so Q
(
0
) can be obtained.(4)  Cycle Process:
(1)
Onp
i
(
i
=
1,2
,
…
,
p
), calculate [
R
i
(
k
)
,
R
i
(
k
)
] and [
Z
i
(
k
)
,
Z
i
(
k
)
]; we can get [
R
(
k
)
,
R
(
k
)
] and [
Z
(
k
)
,
Z
(
k
)
] after all-reduce. If ∥
R
(
k
)
∥
<
ε, or ∥
R
(
k
)
∥
>
ε but ∥
Z
(
k
)
∥
<
ε stop; otherwise, turn back to (
2
).
(2)
Calculateγ
k
=
[
R
(
k
)
,
R
(
k
)
]
/
[
Z
(
k
)
,
Z
(
k
)
] in each processor.
(3)
Onp
i
(
i
=
1,2
,
…
,
p
) calculate X
i
(
k
+
1
)
=
X
i
(
k
)
+
γ
k
Z
i
(
k
).
(4)
Then getA
~
X
(
k
+
1
) and (
A
~
X
(
k
+
1
)
)
B
~; Processor p
i
(
i
=
1,2
,
…
,
p
) calculates R
i
(
k
+
1
) and [
R
i
(
k
+
1
)
,
R
i
(
k
+
1
)
], and finally obtain [
R
(
k
+
1
)
,
R
(
k
+
1
)
] after all-reduce.
(5)
After thatA
~
T
R
(
k
+
1
) and (
A
~
T
R
(
k
+
1
)
)
B
~
T are done; then calculate Q
i
(
k
+
1
) on p
i
(
i
=
1,2
,
…
,
p
).
(6)
Calculate the factorη
k
=
[
R
(
k
+
1
)
,
R
(
k
+
1
)
]
/
[
R
(
k
)
,
R
(
k
)
] in each processor; then on processor p
i
(
i
=
1,2
,
…
,
p
), renew the value Z
i
(
k
+
1
)
=
Q
i
(
k
+
1
)
+
η
k
Z
i
(
k
+
1
).
(7)
Letk
:
=
k
+
1; turn back to (
1
).Note that if n is not divisible by p, some processors store [
n
/
p
]
+
1 rows-block of A, B and I, sequentially, and others store [
n
/
p
] rows-block to make the load balance of each processor and short wait-time.


## 5. Convergent Analysis

Let the matricesA and B be both positive definite symmetric matrices. λ
i, μ
j
(
i
,
j
=
1,2
,
…
,
n
) are the eigenvalues of A and B, respectively. λ
i,  μ
j
(
i
,
j
=
1,2
,
…
,
n
) satisfy λ
n
≥
⋯
≥
λ
2
≥
λ
1
>
0, μ
n
≥
⋯
≥
μ
2
≥
μ
1
>
0.For the matrix equationsA
X
+
X
B
+
F
=
O, the speed of the convergence depends on
(14)
max
⁡
|
λ
i
+
μ
j
|
min
⁡
|
λ
i
+
μ
j
|
=
λ
n
+
μ
n
λ
1
+
μ
1
.
While our proposed matrix forms X
=
A
~
X
B
~
+
F
~, its convergent rate is decided by
(15)
max
⁡
i
,
j
⁡
|
1
-
(
(
α
λ
i
+
1
)
/
(
α
λ
i
-
1
)
)
(
(
α
μ
j
+
1
)
/
(
α
μ
j
-
1
)
)
|
min
⁡
i
,
j
⁡
|
1
-
(
(
α
λ
i
+
1
)
/
(
α
λ
i
-
1
)
)
(
(
α
μ
j
+
1
)
/
(
α
μ
j
-
1
)
)
|
.
The following comparison could indicate that our algorithm is more effective.Theorem 6.
If the eigenvalues of matricesA and B satisfy
(16)
λ
1
λ
n
≤
(
5
-
1
2
)
2
,
μ
1
μ
n
≤
5
-
1
2
,
then there exists α such that the inequality
(17)
max
⁡
i
,
j
|
1
-
(
(
α
λ
i
+
1
)
/
(
α
λ
i
-
1
)
)
(
(
α
μ
j
+
1
)
/
(
α
μ
j
-
1
)
)
|
min
⁡
i
,
j
⁡
|
1
-
(
(
α
λ
i
+
1
)
/
(
α
λ
i
-
1
)
)
(
(
α
μ
j
+
1
)
/
(
α
μ
j
-
1
)
)
|
<
λ
n
+
μ
n
λ
1
+
μ
1
holds.Proof.
Since|
(
(
α
λ
i
+
1
)
/
(
α
λ
i
-
1
)
)
(
(
α
μ
j
+
1
)
/
(
α
μ
j
-
1
)
)
|
<
1, for α
<
0, we have
(18)
max
⁡
i
,
j
|
1
-
(
(
α
λ
i
+
1
)
/
(
α
λ
i
-
1
)
)
(
(
α
μ
j
+
1
)
/
(
α
μ
j
-
1
)
)
|
min
⁡
i
,
j
⁡
|
1
-
(
(
α
λ
i
+
1
)
/
(
α
λ
i
-
1
)
)
(
(
α
μ
j
+
1
)
/
(
α
μ
j
-
1
)
)
|
≤
2
min
⁡
i
,
j
⁡
|
1
-
(
(
α
λ
i
+
1
)
/
(
α
λ
i
-
1
)
)
(
(
α
μ
j
+
1
)
/
(
α
μ
j
-
1
)
)
|
≤
2
1
-
ρ
(
A
~
⊗
B
~
T
)
.
From [17], we have known that ρ
(
A
~
⊗
B
~
T
)
=
(
(
λ
n
-
λ
1
)
/
(
λ
n
+
λ
1
)
)
(
(
λ
1
λ
n
-
μ
1
)
/
(
λ
1
λ
n
+
μ
1
)
), for α
=
-
1
/
λ
1
λ
n. Since λ
1
/
λ
n
≤
(
(
5
-
1
)
/
2
)
2 andμ
1
/
μ
n
≤
(
5
-
1
)
/
2, we have that
(19)
λ
n
(
λ
n
-
λ
1
λ
n
-
λ
1
)
+
λ
n
(
μ
n
-
μ
1
-
μ
1
λ
1
λ
n
)
>
0
.
So
(20)
λ
n
(
λ
n
+
λ
1
)
(
λ
1
λ
n
+
μ
1
)
≥
1
λ
n
+
μ
n
,
that is,
(21)
2
1
-
ρ
(
A
~
⊗
B
~
T
)
≤
λ
n
+
μ
n
λ
1
+
μ
1
,
for
α
=
-
1
λ
1
λ
n
.
Therefore
(22)
max
⁡
i
,
j
⁡
|
1
-
(
(
α
λ
i
+
1
)
/
(
α
λ
i
-
1
)
)
(
(
α
μ
j
+
1
)
/
(
α
μ
j
-
1
)
)
|
min
⁡
i
,
j
⁡
|
1
-
(
(
α
λ
i
+
1
)
/
(
α
λ
i
-
1
)
)
(
(
α
μ
j
+
1
)
/
(
α
μ
j
-
1
)
)
|
<
λ
n
+
μ
n
λ
1
+
μ
1
for
α
=
-
1
λ
1
λ
n
.From previously mentioned theorem, we can observe that proposed preconditioned MCG is better than traditional MCG ifλ
1
/
λ
n and μ
1
/
μ
n are less, that is, if the singular values of A and B are centralized slightly.


## 6. Numerical Examples and Results Analysis

In this section, the comparison with three methods the parameter iterative method, the MCG algorithm, and the proposed algorithm (PMCG) is discussed for solving two examples. The parallel platform is the parallel machine Lenovo Shenteng 1800 cluster. In the following test cases,X
(
0
)
=
O and termination condition ε
=
1
0
-
10.Example 7.
Consider the elliptical partial differential equation(23)
∂
2
u
∂
x
2
+
∂
2
u
∂
y
2
+
sin
2
π
x
∂
u
∂
x
+
sin
2
π
y
∂
u
∂
y
+
u
=
0
,
∂
2
u
∂
x
2
+
∂
2
u
∂
y
2
+
sin
2
π
x
∂
u
∂
x
00000
+
0
≤
x
,
y
≤
1
,
the boundary condition are
(24)
u
|
x
=
0
=
u
|
x
=
1
=
10
+
cos
⁡
π
y
,
u
|
y
=
0
=
u
|
y
=
1
=
10
+
cos
⁡
π
x
.
Let the step size beh
=
1
/
1001, and five-point difference form is used to discrete above equation, which is finally transformed to solving a Sylvester matrix equation seen in (1). Since the eigenvalues of matrix A and matrix B are all greater than zero, then we have α
=
-
n
/
∥
Α
∥
F
=
-
0.3015 and α
=
-
n
/
∥
B
∥
F
=
-
0.4083, respectively. The numerical results are shown in Table 1.Table 1
Numerical results of Example7.
p
8
12
16
20
Parameter iterative method (α
=
-
0.3015)
T (s)
571.2402
411.2114
310.7856
251.4239
K
33
33
33
33
E
1
0.53
0.49
0.48
0.48
Δ
7.21
e
-
11
7.21
e
-
11
7.21
e
-
11
7.21
e
-
11
Parameter iterative method (α
=
-
0.4083)
T (s)
564.0080
394.3398
301.6541
246.2671
K
31
31
31
31
E
1
0.53
0.51
0.50
0.49
Δ
7.21
e
-
11
7.21
e
-
11
7.21
e
-
11
7.21
e
-
11
MCG algorithm
T (s)
641.1601
427.8546
332.7490
271.4428
K
79
79
79
79
E
1
0.47
0.47
0.46
0.45
Δ
9.80
e
-
11
9.80
e
-
11
9.80
e
-
11
9.80
e
-
11
The algorithm in the paper (α
=
-
0.3015)
T (s)
303.8564
220.1811
176.6607
151.9282
K
20
20
20
20
E
1
1
0.92
0.86
0.80
Δ
9.82
e
-
11
9.82
e
-
11
9.82
e
-
11
9.82
e
-
11
The algorithm in the paper (α
=
-
0.4083)
T (s)
301.6773
216.2561
167.5985
141.9657
K
18
18
18
18
E
1
1
0.93
0.90
0.85
Δ
6.39
e
-
11
6.39
e
-
11
6.39
e
-
11
6.39
e
-
11Example 8.
If matricesA
=
(
a
i
j
)
n
×
n and B
=
(
b
i
j
)
n
×
n, where
(25)
a
i
j
=
{
3
+
sin
(
i
+
j
)
i
=
j
cos
⁡
(
i
+
j
)
others
,
b
i
j
=
{
3
+
2
sin
(
i
+
j
)
i
=
j
2
cos
⁡
(
i
+
j
)
others
,
F is any given n
×
n matrix, n
=
500. The comparison of three different methods for solving Sylvester matrix equation A
X
+
X
B
+
F
=
O is discussed here.
Since the eigenvalues of matrixA and matrix B have not only negative real numbers but also positive ones, so there are four choices for α. We have α
=
-
n
/
∥
Α
∥
F
=
-
0.0621, α
=
n
/
∥
Α
∥
F
=
0.0621, α
=
-
n
/
∥
B
∥
F
=
-
0.0314, and α
=
n
/
∥
B
∥
F
=
0.0314, respectively. The numerical results are shown in Table 2.Table 2
Numerical results of Example8.
p
4
8
12
16
Parameter iterative method (α
=
-
0.0621, 0.0621, -
0.0314, 0.0314)
T (s)
—
—
—
—
K
—
—
—
—
E
1
—
—
—
—
Δ
—
—
—
—
MCG algorithm
T (s)
612.0169
336.6452
237.2380
179.1261
K
331
331
331
331
E
1
0.23
0.21
0.20
0.20
Δ
8.35
e
-
11
8.35
e
-
11
8.35
e
-
11
8.35
e
-
11
The algorithm in the paper (α
=
-
0.0621)
T (s)
144.6455
79.1208
54.7878
41.6944
K
80
79
81
81
E
1
1
0.91
0.88
0.87
Δ
3.64
e
-
11
5.31
e
-
11
9.44
e
-
11
9.44
e
-
11
The algorithm in the paper (α
=
0.0621)
T (s)
221.1802
120.1056
84.8724
64.2965
K
134
133
133
133
E
1
1
0.92
0.87
0.86
Δ
7.70
e
-
11
9.87
e
-
11
9.87
e
-
11
9.87
e
-
11
The algorithm in the paper (α
=
-
0.0314)
T (s)
169.8777
91.0232
65.0475
48.8154
K
104
103
103
103
E
1
1
0.93
0.87
0.87
Δ
8.07
e
-
11
6.30
e
-
11
6.30
e
-
11
6.30
e
-
11
The algorithm in the paper (α
=
0.0314)
T (s)
245.1565
132.7123
93.1277
70.7914
K
157
154
154
154
E
1
1
0.92
0.88
0.87
Δ
7.03
e
-
11
5.91
e
-
11
5.91
e
-
11
5.91
e
-
11
Note:
(1) Herep, T (s),  K, E1, and Δ represent the process number, time, iteration number, absolute parallel efficiency, and error, respectively.
(2) Because the computational time of one processor is too long, we take the minimal time with one of multiprocessors as a reference to calculate the relative speedup and relative parallel efficiency.


### 6.1. Results Analysis

From the results shown in the above tables, we observed the following.(i)
Compared with the modified conjugate gradient method and parameter iterative method, the iteration number and the computational time by our algorithm are significantly reduced. It is shown that the preconditioned method enhances its validity and convergent rate. Obviously, better parallel efficiency can be obtained with our algorithm.
(ii)
When the eigenvalues of matrixA and matrix B are all positive real numbers (seen in Example 7), the parallel efficiency of our algorithm is much higher than it with traditional MCG. Furthermore, the numerical results are consistent with the theory; in the case that the eigenvalues have not only negative real numbers but also positive ones in Example 8, the parallel efficiency of our algorithm is still higher than that with MCG. That means the preconditioned factor α is efficient.
(iii)
From the results in the above Tables, we can observe that the parameter method is ineffective in Example8. This illustrates that our algorithm has a wider range of applications than parameter method.


## 6.1. Results Analysis

From the results shown in the above tables, we observed the following.(i)
Compared with the modified conjugate gradient method and parameter iterative method, the iteration number and the computational time by our algorithm are significantly reduced. It is shown that the preconditioned method enhances its validity and convergent rate. Obviously, better parallel efficiency can be obtained with our algorithm.
(ii)
When the eigenvalues of matrixA and matrix B are all positive real numbers (seen in Example 7), the parallel efficiency of our algorithm is much higher than it with traditional MCG. Furthermore, the numerical results are consistent with the theory; in the case that the eigenvalues have not only negative real numbers but also positive ones in Example 8, the parallel efficiency of our algorithm is still higher than that with MCG. That means the preconditioned factor α is efficient.
(iii)
From the results in the above Tables, we can observe that the parameter method is ineffective in Example8. This illustrates that our algorithm has a wider range of applications than parameter method.



---

*Source: 598716-2014-03-26.xml*