Journal of Mathematical Research with Applications
Jan., 2021, Vol. 41, No. 1, pp. 87–98
DOI:10.3770/j.issn:2095-2651.2021.01.009
Http://jmre.dlut.edu.cn

An Improvement of Multigrid Methods Using Multiple
Grids on Each Layer for Parallel Computing
Akira IMAKURA
University of Tsukuba, 1-1-1 Tennodai, Tsukuba, Ibaraki 305-8573, Japan
Abstract Multigrid methods are widely used and well studied for linear solvers and preconditioners of Krylov subspace methods. The multigrid method is one of the most powerful approaches for solving large scale linear systems; however, it may show low parallel efficiency on coarse
grids. There are several kinds of research on this issue. In this paper, we intend to overcome
this difficulty by proposing a novel multigrid algorithm that has multiple grids on each layer.
Numerical results indicate that the proposed method shows a better convergence rate comparedwith the existing multigrid method.
Keywords

linear systems; multigrid; preconditioner

MR(2020) Subject Classification 65F08; 65F10

1. Introduction
This paper focuses on algebraic multigrid (AMG) methods [1,2] for solving linear systems of
the form
Ax = b, A ∈ Rn×n , b, x ∈ Rn .

(1.1)

Algebraic multigrid methods and multigrid preconditioned Krylov subspace methods are a powerful approach for solving large scale linear systems. However, it has a difficulty in parallel
computing. That is, using the methods, we have that parallel efficiency decreases in computations of problems on coarse layers. To improve the parallel efficiency in coarse grids, several typesof research have been proposed such as the use of efficient parallel smoothers and the restriction
techniques of parallelism in coarse grids [3–8].
In this paper, to improve the parallel efficiency and convergence rate of algebraic multigrid
methods, we propose a novel multigrid algorithm that has multiple grids on each layer. The
proposed method increases computational costs due to the use of multiple grids in coarse layers
and thus, improves the parallel efficiency. The proposed method is expected to be more efficientthan the existing multigrid methods regarding computation time in parallel computing because
the computation for coarse grids can be parallelized. We also propose a practical algorithm basedon the smoothed aggregation AMG (SA-AMG) method [9]. Notice that, since the main purpose
Received December 18, 2019; Accepted April 23, 2020
Supported by the Japan Science and Technology Agency (JST), ACT-I (Grant No. JPMJPR16U6) and theJapan Society for the Promotion of Science (JSPS), Grants-in-Aid for Scientific Research (Grant Nos. 17K12690;
18H03250; 19KK0255).
E-mail address: imakura@cs.tsukuba.ac.jp

88

Akira IMAKURA

Algorithm 1 Algebraic two-grid method
Function: x = A2G(A, x(0) , b, A1 , P1 , R1 )1: Pre-smoother:x := smoothν1 (A, x(0) , b)
2: Get residual:r = b − Ax
3: Coarsen:
b1 = R1 r
4: Solve:
Solve A1 y1 = b1 for y1
5: Correct:
x := x + P1 y1
6: Post-smoother:x := smoothν2 (A, x, b)of this paper is to propose a basic algorithm of the proposed method and demonstrate that the
proposed method can improve the convergence rate. The parallel efficiency evaluation of the
proposed method is open to research.
The remainder of this paper is organized as follows. Section 2 briefly reviews the algebraic
multigrid methods. In Section 3, we propose a basic concept of the proposed method anda practical algorithm based on the SA-AMG method [9]. Also, we analyse the relationshipbetween the proposed method and the existing multigrid methods. In Section 4, we describe thenumerical experiments used to evaluate the convergence performance of the proposed method.
Finally, conclusions are given in Section 5.

2. Algebraic multigrid methods
Let x(0) be an approximate solution of the linear system (1). Then, a two-grid method isthe algorithm that considers updating the approximate solution x(0) of the linear systems (1) bysolving one coarse grid linear system with small matrix A1 ∈ Rn1 ×n1 (n1 < n), that is,
A1 y 1 = b 1
as shown in Algorithm 1, where R1 ∈ Rn1 ×n and P1 ∈ Rn×n1 are called as a restriction matrix andan interpolation matrix, respectively. Here, the function smoothν (A, x(0) , b) applies ν iterationsof smoother with an initial guess x(0) for solving Ax = b. For a smoother, stationary iterative
methods, the CG method [10] and the GMRES method [11] are usually used.
Multigrid methods are the algorithms that apply the two-grid methods recursively. Let Lmax
be the number of layers and n` , ` = 1, 2, . . . , Lmax be the size of each grid and n0 = n. LetR` ∈ Rn` ×n`−1 be restriction matrices and P` ∈ Rn`−1 ×n` be interpolation matrices. Then, thealgorithm of the algebraic multigrid method is written as Algorithm 2. Algorithm 2 is called theV-cycle version. Instead of the above V-cycle version, the W-cycle version is also well used forachieving a better convergence rate [1,2].
In algebraic multigrid methods, a coefficient matrix in each grid is set as
A` = R` A`−1 P`
and A0 = A. Therefore, updating of the solution in the two-grid method (Algorithm 1 Step 2–5)
is written as−1
x := x + P1 A−1
R1 (b − Ax).
1 R1 r = x + P1 (R1 AP1 )An improvement of multigrid methods using multiple grids on each layer for parallel computing

89

Algorithm 2 Algebraic multigrid method(0)

Function: x` = AMG(`, A` , x` , b` , A`+1 , P`+1 , R`+1 )(0)
1: Pre-smoother:x` := smoothν1 (A` , x` , b` )2: Get residual:r` = b` − A` x`
3: Coarsen:b`+1 = R`+1 r`
4: if ` + 1 = Lmax then5:
Solve:
Solve A`+1 y`+1 = b`+1 for y`+1
6: else7:
Recursion:y` = AMG(` + 1, A`+1 , 0, b`+1 , A`+2 , P`+2 , R`+2 )8: end if
9: Correct:x` := x` + P`+1 y`+1
10: Post-smoother:x` := smoothν2 (A` , x` , b` )This is regarded as Petrov-Galerkin-type projection method using subspaces R(P1 ) and R(R1T ),x ∈ x + R(P1 ), r := b − Ax ⊥ R(R1T ).

3. Proposed method
Here, we propose a basic concept of the proposed method and a practical algorithm basedon the SA-AMG method
3.1. Basic concept
In multigrid methods, the computation on each layer is performed sequentially. Since thesize of problems on coarse grids is relatively small, then the parallel efficiency decreases incomputations on coarse grids. To improve the parallel efficiency in coarse grids, several typesof research have been proposed such as usage of efficient parallel smoothers and restriction
techniques of parallelism in coarse grids [3–8].
In this paper, our goal is to improve the parallel efficiency and convergence rate of algebraic
multigrid methods. We propose a novel multigrid algorithm that has multiple grids on each layer.
(1)

(2)

Here, we consider a two-grid method with two couples of restriction matrices R1 , R1 ∈ Rn1 ×n
(1)

(2)

and interpolation matrices P1 , P1
(1)

∈ Rn×n1 and small matrices
(1)

(1)

(2)

(2)

(2)

A1 = R1 AP1 , A1 = R1 AP1 .
Then, ignoring smoother, the updating is written as
(1)

(1)

(1)

(1)

(2)

(2)

(2)

(2)

x(1) = x(0) + P1 (R1 AP1 )−1 R1 (b − Ax(0) ),x(2) = x(1) + P1 (R1 AP1 )−1 R1 (b − Ax(1) ).
This can be regarded as applying two different Petrov-Galerkin-type methods sequentially.
Here, x(2) can be rewritten as
(2)

(2)

(2)

(2)

(2)

(2)

(2)

(2)

x(2) = x(1) + P1 (R1 AP1 )−1 R1 (b − Ax(1) )
(1)

(1)

(1)

(1)

= x(1) + P1 (R1 AP1 )−1 R1 {b − A[x(0) + P1 (R1 AP1 )−1 R1 (b − Ax(0) )]}.

90

Akira IMAKURA

Algorithm 3 Proposed two-grid method
(1)

(2)

(1)

(2)

(1)

(2)

Function: x = PA2G(A, x(0) , b, A1 , A1 , P1 , P1 , R1 , R1 )1: Pre-smoother:x := smoothν1 (A, x(0) , b)
2: Get residual:r = b − Ax
(1)
(1)
3: Coarsen:
b1 = R1 r(2)(2)
4: Coarsen:
b1 = R1 r(1) (1)(1)(1)
5: Solve:
Solve A1 y1 = b1 for y1(2) (2)(2)(2)
6: Solve:
Solve A1 y1 = b1 for y1(1) (1)(2) (2)
7: Correct:
x := x + P1 y1 + P1 y1
8: Post-smoother:x := smoothν2 (A, x, b)

Algorithm 4 Proposed multigrid method
(0)

(1)

(2)

(1)

(2)

(1)

(2)

Function: x` = PAMG(`, A` , x` , b` , A`+1 , A`+1 , P`+1 , P`+1 , R`+1 , R`+1 )1: Pre-smoother:
2: Get residual:
3: Coarsen:(0)x` := smoothν1 (A` , x` , b` )r` = b` − A` x`(1)(1)
b`+1 = R`+1 r`
(2)

(2)

4: Coarsen:b`+1 = R`+1 r`
5: if ` + 1 = Lmax then
(1)
(1)
(1)
(1)
6:
Solve:
Solve A`+1 y`+1 = b`+1 for y`+1
(2)

10:Recursion:
11: end if
12: Correct:
13: Post-smoother:

(2)

(2)

(2)

Solve A`+1 y`+1 = b`+1 for y`+1

7:
Solve:
8: else9:
Recursion:

(1)

(1)
(1)
(1,1)
(1,2)
(1,1)
(1,2)
(1,1)
(1,2)
= PAMG(` + 1, A`+1 , 0, b`+1 , A`+2 , A`+2 , P`+2 , P`+2 , R`+2 , R`+2 )
(2)
(2)
(2)
(2,1)
(2,2)
(2,1)
(2,2)
(2,1)
(2,2)
y` = PAMG(` + 1, A`+1 , 0, b`+1 , A`+2 , A`+2 , P`+2 , P`+2 , R`+2 , R`+2 )

y`

(1)

(1)

(2)

(2)

x` := x` + P`+1 y`+1 + P`+1 y`+1x` := smoothν2 (A` , x` , b` )(2)Now, we assume that the matrices R1(1)

and P1(2)

satisfy(1)

R1 AP1

= O.

(3.1)

In this case, the computation of x(2) becomes
(2)

(2)

(2)

(2)

(1)

(1)

(1)

(1)

x(2) = x(1) + P1 (R1 AP1 )−1 R1 (b − Ax(0) )
(2)

(2)

(2)

(2)

= x(0) + P1 (R1 AP1 )−1 R1 (b − Ax(0) ) + P1 (R1 AP1 )−1 R1 (b − Ax(0) ).
Here, we can observe that the 2nd and 3rd terms,
(1)

(1)

(1)

(1)

(2)

(2)

(2)

(2)

P1 (R1 AP1 )−1 R1 (b − Ax(0) ),P1 (R1 AP1 )−1 R1 (b − Ax(0) )can be computed in parallel. This idea follows the proposed two-grid and multigrid methods asshown in Algorithms 3 and 4.
Comparing Algorithms 2 and 4, we have that the convergence rate of the proposed method
is expected to be improved, although the proposed method increases computational costs forcomputing multiple coarse grids. On the other hand, since the computation of problems oncoarse grids of the proposed method can be parallelized, then the proposed method may beefficient regarding computation time in parallel computing.
3.2. Relationship with existing methods

An improvement of multigrid methods using multiple grids on each layer for parallel computing

91

Here, we analyse the relationship between the proposed method (Algorithm 4) and theexisting multigrid method (Algorithm 2) with two times larger subspace than the proposed
method. We observe that from the condition (3.1) that updating the solution using the subspace
(1)

(2)

(1)

(2)

R([P1 , P2 ]) and R([(R1 )T , (R2 )T ]) is written as" (1) #! "#hihi −1 R(1) R1(1) (2)(1) (2)1(1)(0)
x = x + P1 P1
A P1 P 1
b − Ax(0)
(2)
(2)
R1
R1
" (1)
#
"
#
h

i R AP (1) R(1) AP (2) −1 R(1) (1) (2)
1
1
1
1
1
(0)
= x + P1 P1
b − Ax(0)
(2)
(1)
(2)
(2)
(2)
R1 AP1R1 AP1
R1
" (1)
#
"
#

h
i R AP (1) R(1) AP (2) −1 R(1) (1) (2)
1
1
1
1
1
(0)
= x + P1 P1
b − Ax(0) .
(2)
(2)
(2)
O
R1 AP1
R1
In general, for the case of symmetric positive definite A, the restriction matrix R` is set as
R` = cP`T , c ∈ R.
In this case, if the condition (3.1) is satisfied, then we have
(1)

(2)

R1 AP1

=O

and
x

(1)

=x

(0)

=x

(0)

h
i
(1) (2)+ P1 P1

"

(1)

(1)

R1 AP1

O

O

(2)
(2)
R1 AP1

#−1 "

(1)

R1

(2)
R1

#



b − Ax(0)



(1)
(1)
(1)
(1)
(2)
(2)
(2)
(2)
+ P1 (R1 AP1 )−1 R1 (b − Ax(0) ) + P1 (R1 AP1 )−1 R1 (b − Ax(0) ).

Thus, for the two-grid methods, the proposed method is mathematically equivalent to the existing
two-grid method with two times larger subspace. However, we note that the relationship of themultigrid methods is more complicated due to smoother.
From the comparison with the existing methods, we obtained the following properties of the
proposed multigrid method:• The computational cost of the proposed multigrid method is two times higher than that of
the existing method using the same size subspace. But the proposed multigrid method is almost
the same as that of the existing method using the two times larger subspace.
• The convergence rate of the proposed multigrid method is better than that of the existing
method using the same size subspace. But the proposed multigrid method is almost the same asthat of the existing method using the two times larger subspace.
• The proposed method has a better parallel efficiency than the existing methods.
3.3. Practical algorithm based on the SA-AMG method
Here, we then propose a practical algorithm based on the typical algebraic multigrid method,
the smoothed aggregation AMG (SA-AMG) method [9]. In the SA-AMG method, we set the
restriction matrix and the coefficient matrix for each grid asR` = P`T , A` = P`T A`−1 P` .

92

Akira IMAKURA

Then, we recursively develop the interpolation matrix P` using the following five steps:[SA-AMG]
Step 1. Let L be the number of initial vectors, and set an initial matrix B1 ∈ Rn×L .
Step 2. Separate the index set into N1a subsets, called “aggregates”, such that indexes ineach aggregate are strongly connected.
Step 3. Using the QR factorization of each sub-matrix corresponding to each aggregate ofT
T
T
T
B1 = [B1,1, B1,2, . . . , B1,Na] ,1

B1,i = Q1,i R1,i , i = 1, 2, . . . , N1a ,the matrix B1 is decomposed into
 

Q1,1
B1,1
 

 B1,2   
B1 = 
=
..
 

.
 

B1,N1a  R1,2 




..


.


R1,N1aQ1,N1a

Q1,2
..

R1,1

.

=: Pe1 B2 .
aaThen, we obtain matrices Pe1 ∈ Rn×(LN1 ) and B2 ∈ R(LN1 )×L .
Step 4. Applying a smoothing technique, we obtain the interpolation matrix P1 . A simple
and appropriate choice for the prolongation smoother is Richardson’s method, that is,a4A)Pe1 ∈ Rn×(LN1 ) ,3λwhere λ is an upper bound on the spectral radius of the matrix, λ ≥ max λ(A) (see [12]).

P1 = (I −Step 5. Applying Steps 2–4 with B2 as the initial matrix, we get the interpolation matrixaaaP2 ∈ R(LN1 )×(LN2 ) and the matrix B3 ∈ R(LN1 )×L . Repeating these steps inductively, we obtain
the interpolation matrices P1 , P2 , . . . , PLmax .
Following the same approach, the proposed multigrid method generates interpolation matrices
(1)

(2)

by splitting the matrix B` into B` = [B` , B` ] in each layer. The algorithm of the proposed
method based on the SA-AMG method can be shown below.
[The proposed method]
(1)

(2)

Step 1. Let L be the number of initial vectors, and set an initial matrix B1 = [B1 , B1 ] ∈
(1)

(2)

Rn×L (B1 , B1

(1)

∈ Rn×(L/2) ). Note that B1

is the matrix whose columns are the first L/2columns of the matrix B1 .
Step 2. Separate the index set into N1a subsets, called “aggregates”, such that indexes ineach aggregate are strongly connected.
(1)

(2)

Step 3. Apply Steps 3 and 4 of the SA-AMG method to the matrix B1 , B1
(2)
(1)
interpolation matrices P1 , P1

∈ R

n×(LN1a /2)(1)(2)
and the matrices B2 , B2

and obtain thea∈ R(LN1 /2)×(L/2) ,respectively.
(1)

(2)

Step 4. For satisfying (P1 )T AP1
(2)

P1

(2)

= P1

(2)

as follows:

(1)

(2)

= O, update P1(1)

T− P1 A−1
1 (P1 ) AP1 .

(3.2)

An improvement of multigrid methods using multiple grids on each layer for parallel computing
(1)

(2)

Step 5. Applying Steps 2–4 with B2 , B2

93

as the initial matrices, we get the interpolation

a
a
(1)
(2)
(3)
(4)
matrices P2 , P2 , P2 , P2 ∈ R(LN1 /2)×(LN2 /4) . Repeating these steps inductively, we obtain
(1)
(2)
(1)
(2Lmax )the interpolation matrices P1 , P1 , . . . , PLmax , . . . , PLmax .

Since the proposed method solves the linear systems with A1 in (3.2), the computational costsof the proposed method increase for constructing interpolation matrices. This computational cost
is expected to be reduced using the rough solution of the linear system with A1 . We will evaluate
the potential in Section 4.2.

4. Numerical experiments
In this section, we evaluate the performance of the proposed method introduced in Section 3.3.
We also compare it with the two types of SA-AMG method which are SA-AMG(1) and SAAMG(2). SA-AMG(1) has the same grid size as the proposed method for each layer, while the
SA-AMG(2) has the grid size as the same as some of the grid sizes in each layer of the proposed
method. The details are shown in Table 1, where L is the number of initial vectors. We test
both the V-cycle and W-cycle.

Number of aggregates
Size of each grid
Number of each grid

SA-AMG(1)SA-AMG(2)proposed16, 4, 1

32, 16, 8

32, 16, 8

16L, 4L, L

32L, 16L, 8L

16L, 4L, L

1, 1, 1

1, 1, 12, 4, 8

Table 1 Number of aggregate, size of each grid and number of each grid of SA-AMG preconditioners

We remark that the main purpose of this paper is to propose a basic algorithm of the proposed
method and demonstrate that the proposed method can improve the convergence rate. In thispaper, we evaluate the convergence rate of the proposed method and compare it with the SAAMG method. Parallel efficiency of the proposed method will be evaluated in the future.
All numerical experiments were performed by MATLAB2018b.
4.1. Experiment I
Here, we evaluate the convergence rate and the eigenvalue distribution that affects the efficiency of the preconditioned Krylov subspace methods. For the model problem, we used thefinite differential two dimensional Laplacian obtained from the MATLAB command
A = delsq(numgrid(’N’,52));.
The problem size is n = 2500. We set b as random vector and x0 = [0, 0, . . . , 0]T .
For parameters, let the number of initial vectors L = 8 and the number of grids Lmax = 3.
We used the weighted Jacobi method for pre/post-smoother and the number of iterations as 1.
We tuned the weight parameter of the weighted Jacobi method using the auto-tuning techniquein [13,14]. We also set B1 as the rough solution of the homogeneous system
Ax = 094

Akira IMAKURA

102
SA-AMG(1)SA-AMG(2)Proposed010

Relative residual 2-norm

Relative residual 2-norm

102

10-2
-4

10

-6

10

10-8
10-10

0

100

200

300

Number of iterations

400

500

10

SA-AMG(1)SA-AMG(2)Proposed

0

10-2
-4

10

-6

10

10-8
10-10

0

100

(a) V-cycle

200

300

Number of iterations

400

500

(b) W-cycle

Figure 1 Convergence history of multigrid methods(a) V-cycle(b) W-cycle

Figure 2 Eigenvalue distributions of the preconditioned matrices with the multigrid preconditioners (1:non-preconditioner, 2: SA-AMG(1), 3: SA-AMG(2), 4: proposed)

obtained by 10 iterations of the weighted Jacobi method. The stopping criterion was set as
krk k2 /kbk2 ≤ 10−10 .
Figures 1 and 2 show the convergence history of the multigrid methods and the eigenvalue
distribution of the preconditioned matrix with the multigrid methods for both V-cycle and Wcycle. Figure 1 demonstrates that the proposed method achieves a better convergence rate thanthe SA-AMG(1) method as well as the SA-AMG(2) method. This numerical result supports thediscussion in Section 3.2. We can also see from Figure 2 that the proposed method improves the
eigenvalue distribution as well as the SA-AMG(2) method.
4.2. Experiment II
In this experiment, we evaluate the convergence rate of the proposed method using the roughsolutions of (3.2). For this purpose, we add the perturbation in the small matrix A1 likeA1 → A1 + εkA1 kE
and evaluate the convergence rate. Let ε denote the amount of the perturbation and E be arandom matrix whose entries are in [−1, 1].

An improvement of multigrid methods using multiple grids on each layer for parallel computing(a) V-cycle

95(b) W-cycle

Figure 3 Number of iterations of the proposed method with rough linear solutions

We used the same problem and the same parameters as in Experiment I. Figure 3 shows
the number of iterations for each ε. From the numerical result, we observed that the number of
iterations increases for large perturbation ε ≥ 10−3 , but does not change for small perturbationε ≤ 10−4 . This result indicates that we can reduce the computational costs using the roughsolution for computing (3.2). In practice, we can use some iterative solvers for efficiently obtaining
an approximation solution.
4.3. Experiment III
We evaluate the convergence of the preconditioned Krylov subspace method with multigrid
methods for the following eight model problems.
Problem 1 [SPD] The finite differential two dimensional Laplacian (n = 10000) obtained
from the MATLAB command
A = delsq(numgrid(’N’,102));.
Problem 2 [SPD] Pres Poission (Computational fluid dynamics problem, n = 14822).
Problem 3 [SPD] S1RMQ4M1 (Finite element analysis of cylindrical shells, n = 5489).
Problem 4 [SPD] S2RMQ4M1 (Finite element analysis of cylindrical shells, n = 5489).
Problem 5 [Non Symm] CAVITY05 (Computational fluid dynamics problem, n = 1182).
Problem 6 [Non Symm] PDE2961 (2D/3D problem, n = 2961).
Problem 7 [Non Symm] RAEFSKY1 (Computational fluid dynamics problem, n = 3242).
Problem 8 [Non Symm] RAEFSKY2 (Computational fluid dynamics problem, n = 3242).
Here, Problems 2–8 are obtained from Matrix Market [15]. Let SPD and Non Symm denotesymmetric positive definite and non-symmetric problems. We set b as random vector and x0 =[0, 0, . . . , 0]T .
To evaluate the performance of the preconditioned Krylov subspace method with multigrid
method, we test the V-cycle of the proposed method, the SA-AMG(1) and SA-AMG(2) methods. We also test the weighted Jacobi preconditioner. For the parameters of the multigrid
preconditioners, we set L = 32 for Problems 1–4 and L = 16 for Problems 5–8. For other pa-

96

Akira IMAKURA

104
non-preconditioner
Weighted Jacobi
SA-AMG(1)SA-AMG(2)Proposed

2

10

100
-2

10

-4

10

-6

10

10-8
10-10

0

100

200

300

Number of iterations

400

Relative residual 2-norm

Relative residual 2-norm

104

10

100
10-2
-4

10

-6

10

10-8
10-10

500

non-preconditioner
Weighted Jacobi
SA-AMG(1)SA-AMG(2)Proposed

20(a) Problem 1: 2D-Laplace

300

400

500

104
non-preconditioner
Weighted Jacobi
SA-AMG(1)SA-AMG(2)Proposed

2

10

0

10

-2

10

-4

10

10-6
-8

10

0

100

200

300

Number of iterations

400

500

(c) Problem 3: S1RMQ4M1

Relative residual 2-norm

Relative residual 2-norm

200

Number of iterations(b) Problem 2: Pres Poisson

104

10-10

100

non-preconditioner
Weighted Jacobi
SA-AMG(1)SA-AMG(2)Proposed

102
10

0

-2

10

10-4
10-6
10-8
10-10

0

100

200

300

Number of iterations

400

500

(d) Problem 4: S2RMQ4M1

Figure 4 Convergence histories for symmetric positive definite problems

rameters of the multigrid preconditioners, we used the same as Experiment I. We set the numberof iterations of the weighted Jacobi preconditioner as 5. For the Krylov subspace method, weused the CG method [10] for SPD problems (Problems 1–4) and the BiCGSTAB method [16] for
non-symmetric problems (Problems 5–8). The stopping criterion was set as krk k2 /kbk2 ≤ 10−10 .

MethodPreconditioner
nonJacobi

SA-AMG(1)SA-AMG(2)Proposed

2D-Laplacian

0.10

3.22

9.90

9.21

11.38

Pre Poisson

——

12.64

42.77

57.72

37.79

S1RMQ4M1

——

4.80

22.14

41.67

32.43

S2RMQ4M1

——

——

——

103.0052.22

CAVITY05

——

0.72

4.31

5.78

4.86

PDE2981

0.03

0.66

2.27

3.07

3.26

RAEFSKY1

——

1.02

2.09

1.83

2.01

RAEFSKY2

——

1.39

1.85

2.50

2.90

Table 2 Computation time [sec.] of the preconditioned CG and BiCGSTAB methods

An improvement of multigrid methods using multiple grids on each layer for parallel computing

104
non-preconditioner
Weighted Jacobi
SA-AMG(1)SA-AMG(2)Proposed

2

10

100
-2

10

-4

10

-6

10

10-8
10-10

0

100

200

300

Number of iterations

400

Relative residual 2-norm

Relative residual 2-norm

104

10

100
10-2
-4

10

-6

10

10-8
10-10

500

non-preconditioner
Weighted Jacobi
SA-AMG(1)SA-AMG(2)Proposed

20(a) Problem 5: CAVITY05

200

300

Number of iterations

400

500

104
non-preconditioner
Weighted Jacobi
SA-AMG(1)SA-AMG(2)Proposed

2

10

0

10

-2

10

-4

10

10-6
-8

10

0

100

200

300

Number of iterations

400

500

(c) Problem 7: RAEFSKY1

Relative residual 2-norm

Relative residual 2-norm

100(b) Problem 6: PDE2961

104

10-10

97

non-preconditioner
Weighted Jacobi
SA-AMG(1)SA-AMG(2)Proposed

102
10

0

-2

10

10-4
10-6
10-8
10-10

0

100

200

300

Number of iterations

400

500

(d) Problem 8: RAEFSKY2

Figure 5 Convergence histories for non-symmetric problems

Figures 4 and 5 show the convergence history. These results demonstrate that the proposedpreconditioned method achieves a better convergence rate than the SA-AMG(1) preconditioner
for both SPD and non-symmetric problems. The proposed preconditioned method shows almost
the same convergence rate as the SA-AMG(2) preconditioner for Pres Poisson and RAEFSKY1,and a better convergence rate for S2RMQ4M1.
Although the main purpose of this paper is to propose a basic algorithm of the proposed
method and demonstrate that the proposed method can improve the convergence rate, here weevaluate the computation time in serial computing. From the numerical results shown in Table 2,
we observed that the proposed preconditioner shows competitive results compared with the SAAMG preconditioners, even in serial computing. Considering also higher parallel efficiency of the
proposed method, the proposed method has a competitive advantage over the existing SA-AMG
methods in parallel computing.

5. Conclusion
In this paper, our primary goal is to improve the parallel efficiency and the convergence
rate of algebraic multigrid methods. To achieve this, we proposed a novel multigrid algorithm

98

Akira IMAKURA

that has multiple grids on each layer. The proposed method increases computational costs
due to the usage of the multiple grids in coarse layers; however, it improves parallel efficiency.
We also proposed a practical algorithm based on the SA-AMG method. From the numericalexperiments, the proposed multigrid method and preconditioner show a better convergence rate
than the existing multigrid method and preconditioner. We note that the proposed method cangenerate more than three grids in each layer.
It is now required to evaluate the parallel efficiency for solving large scale problems. The
practical approximation technique for computing (3.2) is also open for investigation.
Acknowledgements The author would like to thank anonymous referees for useful comments.

References[1] W. L. BRIGGS, V. E. HENSON, S. F. MCCORMICK. A Multigrid Tutorial, 2nd ed. SIAM Publications,
2000.
[2] Y. SAAD. Iterative Methods for Sparse Linear Systems. 2nd edition. SIAM, Philadelphia, PA, 2003.
[3] A. BAKER, T. GAMBLIN, M. SCHULTZ, et al. Challenge of scaling algebraic multigrid across modern
multicore architectures. in: Proceedings of the 2011 IEEE International Parallel & Distributed Processing
Symposium (IPDPS’11), 2011, 275–286.
[4] H. SUNDAR, G. BIROS, C. BURSTEDDE, et al. Parallel Geometric-Algebraic Multigrid on Unstructured
Forests of Octrees. in: ACM/IEEE Proceedings of the 2012 International Conference for High Performance
Computing, Networking, Storage and Analysis (SC12), 2012.
[5] P. T. LIN. Improving multigrid performance for unstructured mesh drift-diffusion simulations on 147,000
cores. Internat. J. Numer. Methods Engrg., 2012, 91(9): 971–989.
[6] V. E. HENSON, U. M. YANG. BoomerAMG: a parallel algebraic multigrid solver and preconditioner. Appl.
Numer. Math., 2002, 41(1): 155–177.
[7] K. NAKAJIMA. Automatic tuning of parallel multigrid solvers using OpenMP/MPI hybrid parallel programming models. Lecture Notes in Computer Science, 2013, 7851: 435–450.
[8] K. NAKAJIMA. Optimization of serial and parallel communications for parallel geometric multigrid method.
in: Proceedings of the 20th IEEE International Conference for Parallel and Distributed Systems (ICPADS2014), 2014, 25–32.
[9] P. VANĚK, J. MANDEL, M. BREAINA. Algebraic multigrid by smoothed aggregation for second and fourthorder elliptic problems. Computing, 1996, 56(3): 179–196.
[10] M. R. HESTENES, E. STIEFEL. Methods of conjugate gradients for solving linear systems. J. ResearchNat. Bur. Standards, 1952, 49: 409–436.
[11] Y. SAAD, M. H. SCHULTZ. GMRES: a generalized minimal residual algorithm for solving nonsymmetriclinear systems. SIAM J. Sci. Stat. Comput., 1986, 7(3): 856–869.
[12] M. BREZINA, R. FALGOUT, S. MACLACHLAN, et al. Adaptive Smoothed Aggregation (αSA). SIAM J.
Sci. Comput., 2004, 25(6): 1896–1920.
[13] A. IMAKURA, T. SAKURAI, K. SUMIYOSHI, et al. An auto-tuning technique of the weighed Jacobi-type
Iteration used for preconditioners of Krylov subspace methods. in: IEEE 6th International Symposium on
Embedded Multicore SoCs (MCSoC-12), 2012, 183–190.
[14] A. IMAKURA, T. SAKURAI, K. SUMIYOSHI, et al. A parameter optimization technique for a weighted
Jacobi-type preconditioner. JSIAM Lett., 2012, 4: 41–44.
[15] Matrix Market. http://math.nist.gov/MatrixMarket/.
[16] H. A. VAN DER VORST. Bi-CGSTAB: a fast and smoothly converging variant of Bi-CG for the solution of
nonsymmetric linear systems. SIAM J. Sci. Statist. Comput., 1992, 13(2): 631–644.