AICS TECHNICAL Report
N0. 2014-001

Block BiCGStab for lattice QCD on the K
computer

By

H. Suno*, Y. Nakamura, K.-I. Ishikawa, andY. Kuramashi
RIKEN Advanced Institute for Computational Science* Corresponding Author E-mail: suno@riken.jp

Submitted on 11/08/2014
Accepted on 03/10/2014
Published and copyrighted by
RIKEN Advanced Institute for Computational Science (AICS)
7-1-26 Minatojima-minami-machi, Chuo-ku, Kobe, 650-0047, Japan

Block BiCGStab for lattice QCD on the K computer
H. Suno∗RIKEN Advanced Institute for Computational Science, Kobe 650-0047, Japanand RIKEN Nishina Center, Wako 351-0198, Japan
Y. Nakamura
RIKEN Advanced Institute for Computational Science, Kobe 650-0047, Japan
K.-I. Ishikawa
Graduate School of Science, Hiroshima University, Hiroshima 739-8526, Japan
Y. Kuramashi
Graduate School of Pure and Applied Science,
University of Tsukuba, Ibaraki 305-8571, Japan,
Center for Computational Science, University of Tsukuba, Ibaraki 305-8577, Japan,
RIKEN Advanced Institute for Computational Science, Kobe 650-0047, Japan(Dated: August 11, 2014)

Abstract
We have developed a block version of the Bi-Conjugate Gradient Stabilized (BiCGStab) solver
for treating the Wilson-Dirac equation with multiple right-hand sides in lattice QCD employing
O(a)-improved Wilson fermions. Code development and optimization are carried out on the K computer at RIKEN Advanced Institute for Computational Science (AICS). We implement the QR
decomposition, a preconditioner based on the domain decomposed Schwarz alternating procedure(SAP) and the symmetric successive over-relaxation (SSOR) iteration with locally lexicographic ordering for the inversion within the domain. The block BiCGStab solver, written in single precision,is included as a preconditioning step for an outer BiCGStab solver in double precision. Numericaltest is performed on the lattice size of 964 using 2048 nodes of the K computer. The performanceof a computational kernel exceeds 50% eﬃciency, and the single precision block BiCGStab solver
shows more than ∼ 30% sustained eﬃciency.

∗

Electronic address: suno@riken.jp

1

I.

INTRODUCTION

The physical point simulation has been one of the primary issues in the ﬁrst principle
calculation of lattice QCD in the recent years. A great eﬀort has been paid to reduce thecomputational cost for the conﬁguration generation with light quark masses. The algorithmic
improvements together with rapid increase of the computational power actually allow us tocarry out a direct full QCD simulation on the physical up and down quark masses [1, 2].
While previous improvements mostly concerns the conﬁguration generation, growing eﬀort
has been also paid to algorithmic improvements for the measurement of physical observables,
so that we are now able to obtain the results for physical quantities at the physical up, down
and strange quark masses [3]. Yet it still remains necessary to reﬁne the results reducing
the systematic errors and to challenge computationally demanding problems, for example,
calculation of the disconnected diagrams.
The main diﬃculty in lattice QCD simulations resides in solving the Wilson-Dirac equation, which is computationally expensive near the physical up and down quark masses andmust be repeated many times both in the conﬁguration generation and measurement ofphysical observables on given conﬁgurations. This can be done in an eﬃcient way only byusing a Krylov subspace method, such as the Bi-Conjugate Gradient Stabilized (BiCGStab)method. The performance of the BiCGStab is usually improved with a suitable preconditioning technique such as the Schwarz alternating procedure (SAP) preconditioner proposedby Lüscher, which is applied to the domain-decomposed lattice [4]. The performance isfurther improved by the nested BiCGStab algorithm with an inner-outer strategy, in whichanother Krylov subspace method is implemented as a preconditioner[5, 6].
In addition, the Wilson-Dirac equation must be solved with multiple sources in the measurement of physical observables on given conﬁgurations: twelve in the simplest case andO(10 − 100) for the stochastic technique. These are typical examples for diﬀerential equations with multiple right-hand sides. For this type of equations, it is well known that theblock Krylov subspace solvers works successfully in reducing the computational cost [7, 8].
Since the Wilson-Dirac matrix in lattice QCD is non-Hermitian, we might expect the block
BiCGStab algorithm [9] to be applicable in a straightforward manner. Recently, new algorithms, the block BiCGGR and block BiCGStab were found to show signiﬁcant improvements for this problem [10, 11].

2

In this work we ﬁrst develop and optimize an inner block BiCGStab solver in single
precision for the K computer at the RIKEN Advanced Institute for Computational Science.
Systematic performance tests are carried out for diﬀerent computational kernels. The solver
is used as a preconditioning step for an outer BiCGStab solver in double precision, where
diﬀerent right-hand sides are treated independently. Numerical tests are performed with a
964 lattice using 2048 nodes of the K computer.
This report is organized as follows. We explain the block BiCGStab algorithm in detail in
Sec. II. Code development and tuning are explained in Sec. III. The results of the numericaltest are presented in Sec. IV. Conclusions and discussions are summarized in Sec. V.

II.

ALGORITHM

The lattice QCD is deﬁned on a hypercubic four-dimensional lattice with the threedimensional spatial extent Lx × Ly × Lz and the temporal extent Lt . The ﬁelds are deﬁned
on the four-dimensional lattice sites n with periodic boundary conditions. We deﬁne two
types of ﬁelds on the lattice. One is the gauge ﬁeld represented by (Uµ (n))a,b with fourdimensional direction indices µ = 1, 2, 3, 4 and color indices a, b = 1, 2, 3, thus which is a
3 × 3 matrix assigned on each link. The other is the quark ﬁeld (q(n))aα that locates on eachsite carrying the Dirac indices α = 1, 2, 3, 4. The O(a)-improved Wilson-Dirac operator
preconditioned with the clover term is represented by the matrix elements
a,b
a,c
Dα,β
(n, m) = δ a,b δα,β δ(n, m) − κFα,γ(n)4∑[(1 − γµ )γ,β (Uµ (n))c,b δ(n + µ̂, m)
µ=1

b,c ∗

+ (1 + γµ )γ,β ((Uµ (m)) ) δ(n − µ̂, m)],(1)

where µ̂ denotes the unit vector in the µ direction in the four-dimensional space-time, F (n)the inverse clover term (1 − (cSW κ/2)σµν Fµν (n))−1 with cSW being an adjusting parameterfor the O(a)-improvement. Our gamma matrices γµ are deﬁned as

0 0 0 −i

0 0 0 −1







 0 0 1 0 
 0 0 −i 0 
γ1 = , , γ2 =  0 1 0 0 
0 i 0 0 
−1 0 0 0
i 0 0 03

0 0 −i 0





1 0 0

0







0 1 0 0 
0 0 0 i 
γ3 = 
 , γ4 = .
 i 0 0 0
 0 0 −1 0 
0 −i 0 0
0 0 0 −1

(2)

The O(a)-improved Wilson-Dirac operator deﬁned above is a complex non-Hermitian squarematrix, where only 51 out of Lx ×Ly ×Lz ×Lt ×3×4 entries in each row have nonzero values.
The matrix is fairly sparse in current numerical simulations with Lx,y,z,t ∼ O(10 − 100).
Calculations of the physical quantities requires the solution of linear equations with multiple right-hand sides:

∑

a,b
Dα,β
(n, m)x(i) (m)bβ = b(i) (n)aα ,

(3)

m,β,b

where b(i) (i = 1, · · · L) represents the i-th source vector. This is expressed in a simple form:
DX = B,(4)

where D is an N × N complex sparse non-Hermitian matrix, X and B are N × L complex
rectangular matrices given by
X = [x(1) , . . . , x(i) , . . . , x(L) ],(5)

B = [b(1) , . . . , b(i) , . . . , b(L) ].

(6)

In the Wilson-Dirac equation, the matrix dimension is given by N = Lx ×Ly ×Lz ×Lt ×3×4,
whereas we designate as L the number of the right-hand sides corresponding to the diﬀerent
source vectors in lattice QCD.
Pseudocode for the block BiCGStab algorithm is described in Algorithm 1 below, which is
improved with two ingredients. One is the QR decomposition with a modiﬁed Gram-Schmidt
orthogonalization method, which improves numerical accuracy since each span works eﬀectively to search approximative solutions. The other is Lüscher’s SAP preconditioner [4],which is applied to the domain-decomposed lattice: The Wilson-Dirac matrix D is decomposed in the 2 × 2 blocked matrix form



D=

DEE DEODOE DOO

4

,

(7)

where the subscripts E and O denote the even and odd domains, respectively. The SAP
preconditioner MSAP is computed as
MSAP = KNSAP∑−1(1 − DK)j

(8)

j=0



with

K=


AEE

0

−AOO DOE AEE AOO

,

(9)

where AEE (AOO ) can be any approximation for (DEE )−1 ((DOO )−1 ). When AEE and AOO
are exact, DK becomes block triangular and is expected to be well conditioned. In thecase |DK| < 1, MSAP converges to D−1 when NSAP → ∞. The approximate inverse of theoperator AEE (AOO ) is the most important part in the SAP. The exactness is not required
for the SAP, however, the better approximation with less computational cost is preferredfor AEE (AOO ). We adopt the symmetric successive over-relaxation (SSOR) method [12],which is derived by decomposing the original operator into the sum of an upper and a lower
triangle matrices. We have then:
AEE = (1 − ωUEE )−1

[N
SSOR
∑

]
j
(1 − ASSOR(1 − ωLEE )−1
EE )

(10)

j=0

with
ASSOR
=
EE

1
[(1 − ωLEE )−1 + (1 − ωUEE )−1 + (ω − 2)(1 − ωLEE )−1 (1 − ωUEE )−1 ],ω(11)where LEE and UEE are the lower and upper triangle matrices corresponding to the forward
and backward hopping terms, respectively. To extract 8 core parallelism for the K computer(for the computer architecture, see Sec. III), we further divide the block into 16 sub-blocks via
the locally-lexicographical ordering (ll -ordering) as described in detail in Refs. [12, 13]. Since(1 − ωLEE ) and (1 − ωUEE ) are triangular matrices, their inverse can be directly calculated
via forward and backward substitutions, respectively. Since more than 80% of the CPU timeis spent for the forward/backward substitutions, optimization must concentrate upon theseparts.

III.

CODE DEVELOPMENT AND OPTIMIZATION FOR THE K COMPUTER

Our code development and optimization are carried out for the K computer at the RIKEN
Advanced Institute for Computational Science. The machine consists of 82944 computa5

Algorithm 1 Block BiCGStab algorithm (D, MSAP , B, ϵi ).
1: initial guess X ∈ CN ×L ,2: compute R = B − DX,
3: set P = R,
4: choose R̃ ∈ CN ×L ,5: while (|r (i) |/|b(i) |) > ϵi do
6:

QR decomposition P = Qγ, P ← Q,

7:

U = MSAP P ,8:

V = DU ,9:solve (R̃H V )α = R̃H R for α,

10:

R ← R − V α,11:

X ← X + U α,

12:

S = MSAP R,

13:

Z = DS,14:

ζ = Tr(ZiH Ri )/Tr(ZiH Zi ),15:

X ← X + ζS,

16:

R ← R − ζZ,

17:solve (R̃H V )β = −R̃H Z for β,18:

P ← R + (P − ζV )β,19: end while

tional nodes and 5184 I/O nodes connected by the so-called “Tofu” network, providing
11.28 Pﬂops of computing capability. The Tofu network topology is six-dimensional one
with 3D-mesh times 3D-torus shape. Each node has a single 2.0GHz SPARC64 VIIIfx processor equipping 8 cores with SIMD enabled 256 registers, 6MB shared L2 cache and 16GB ofmemory. The L1 cache sizes per each core are 32KB/2WAY (instruction) and 32KB/2WAY(data).
We have developed a block version of the forward and backward solvers constructing AEE(or AOO ) in the SSOR for the multiple right-hand sides, based on the non-block version[13].
The forward and backward kernels constitute the most important part in the SSOR. Bothconsists of several four-fold nested loops: three spatial direction x, y, z-loops and one temporal t-loop. They are typically written as
6for ix = 0 to Nx − 1 dofor iy = 0 to Ny − 1 dofor iz = 0 to Nz − 1 do
OpenMP Barrier
for it = 0 to Nt − 1 do
if it > 0 theny ′′ (n)aα = (1 + γ4 )α,β ((U (n − t̂))b,a )∗ y(n − t̂)bβend if...
′′by ′′ (n)aα ⇐ F (n)a,bα,β y (n)βy ′ (n)aα = y(n)aα + κy ′′ (n)aα
end forend forend forend for,or as the same way with all the loops in the reversed order, where Nx,y,z,t denote the numberof lattice sites in x, y, z, t directions inside the SAP block. Note that an explicit OpenMP
barrier is inserted because of the recurrence dependency among the OpenMP threads. There
are several conditional branches for site location like it > 0, ix > 0 and ix < Nx − 1, etc...
inside the t-loop. The question arises as to where to put the loop running over diﬀerent
right-hand side vectors in Eq. (4). Our choice is just inside the t-loop like (we designate ascase (a)):
for it = 0 to Nt − 1 dofor i = 0 to L − 1 do
if it > 0 theny ′′ (n)aα = (1 + γ4 )α,β ((U (n − t̂))b,a )∗ y (i) (n − t̂)bβend if...
′′by ′′ (n)aα ⇐ F (n)a,bα,β y (n)βy ′(i) (n)aα = y(n)aα + κy ′′ (n)aα
end for7

end for,or even inside the if statements (we designate as case (b)):
for it = 0 to Nt − 1 do
if it > 0 thenfor i = 0 to L − 1 doy ′′(i) (n)aα = (1 + γ4 )α,β ((U (n − t̂))b,a )∗ y (i) (n − t̂)bβ
end forend if...
for i = 0 to L − 1 do′′(i)y ′′(i) (n)aα ⇐ F (n)a,b
(n)bβ
α,β y

y ′(i) (n)aα = y (i) (n)aα + κy ′′(i) (n)aα
end forend for.
We have carried out performance tests for these two cases with lattice sizes 123 × 24 and243 × 48, where the SAP block sizes correspond to 64 and 124 , respectively, on 16 nodes of
the K computer choosing L = 1, 4 and 12 for the number of the right-hand sides. For all
the performance tests presented in this Section, we choose the hopping and improvement
parameters (κ, cSW ) = (0.123, 1.0). Table I shows eﬃciency, SIMD rate and several time
consuming parts. Although it is diﬃcult to observe the general trend because of the complexity of their behavior with respect to L, the eﬃciency and SIMD rate increase well with
L. We notice that (OpenMP) barrier synchronization wait and integer cache wait decrease
with L. This is because the L-loop does not include any OpenMP barrier statement inside
it. For all the values of L the better performance is found for case (a) than for case (b);We obtain a better SIMD rate and less time for ﬂoating-point cache wait. We should also
note that case (a) has less loop overhead than case (b) so that SIMDization and software
pipelining should work better. We have hence chosen to put the L-loop inside the t-loop andoutside the if statements, as in case (a). Since the DEE and DOO kernels involve a similar
four-fold nested loop with several if statements inside, though without OpenMP barrier
synchronization, we have inserted the L-loop inside the t-loop and outside if statements.
We have measured the performance of the AEE and DEE kernels as well as the DMSAP89

4 instr.

2/3 instr.

1 instr.

34

36

19

25

26

17

21

24

4

12

(b) 123 × 24 1

4

12

243 × 48 1

4

12

40

12

29

39

4

243 × 48 1

31

75.0

73.6

60.9

74.7

73.0

59.4

87.8

84.9

77.0

88.0

86.8

78.1

8.8

7.4

8.1

7.8

7.7

8.5

13.0

11.7

10.6

13.8

14.1

11.4

19.9

19.8

17.2

23.4

23.9

20.3

23.9

23.5

21.7

27.6

25.1

23.2

20.5

18.0

16.1

22.6

21.8

18.0

19.0

17.7

15.6

20.4

22.2

15.2

[%] rate[%] commit[%] commit[%] commit[%](a) 123 × 24 1

size

Case Lattice L Eﬃc. SIMD

5.3

7.7

7.6

8.8

10.4

10.8

5.7

9.5

9.3

9.5

10.3

10.5

wait[%]

8.3

8.3

8.8

8.9

9.3

10.0

13.3

11.5

10.1

12.1

15.5

11.6

24.2

22.5

14.8

16.8

14.7

14.6

18.1

16.1

13.3

11.4

6.8

14.3

0.7

1.8

3.6

1.1

2.7

5.8

0.9

2.0

4.6

1.1

2.3

5.4

LA wait[%] LA wait[%] LA wait[%]

4.7

6.4

10.5

0.9

0.2

0.0

1.6

2.6

9.5

1.1

0.2

0.0

LA wait[%]Barrier sync. Fl. pt. op. Fl. pt. cache Int. cache Fl.pt. memory

L-loop is located just inside the t-loop in case (a), while inside the if statements in case (b). LA stands for load access.

TABLE I: Results for the performance tests of the block AEE (AOO ) kernel on 16 nodes of the K computer. We have tested two cases: Themultiplication and the block BiCGStab algorithm with two lattice sizes of 123 × 24 and243 × 48 on 16 nodes of the K computer. Figure 1 presents the performance eﬃciencies withrespect to the ﬂoating point number operations for these kernels as a function of L. All theperformance eﬃciencies are increased with L. We should note that the diﬀerence between
the block BiCGStab and DMSAP slightly becomes larger as L increases. This is mostly due
to the QR decomposition, which requires more MPI Allreduce calls for increased numberof L. We have also carried out weak-scaling performance tests by increasing the numberof nodes with the SAP block size ﬁxed to 64 and 124 at L = 12: We have chosen lattices
sizes of 243 × 48 on 256 nodes, 484 on 2048 nodes, 483 × 96 on 256 nodes, and 964 on 2048
nodes. The results are shown in Fig. 2. We observe that the performance drops little withincreasing number of nodes for most of the computational kernels. Even in the case of theblock BiCGStab with the SAP block size of 64 , which shows the worst scaling behavior, theeﬃciency is still better than 30% for the largest lattice size. The numerical results of these
performance tests are summarized in Table II, together with the measured SIMD rate andseveral time consuming parts.
The nested BiCGStab algorithm is designed to be ﬂexible against the preconditioner
changing iteration by iteration [5, 6], where the outer BiCGStab contains an inner BiCGStab
as the ﬂexible preconditioner. The ﬂexibility ensures the double precision accuracy of thesolution vector even if we used the single precision for the inner BiCGStab solver. In
our work, the block BiCGStab solver described above is used as a preconditioner or innersolver for the outer solver, the former in single precision and the latter in double precision.
Pseudocode for the outer solver is described in Algorithm 2. Note that the preconditioning
steps at lines 8 and 14 in Algorithm 2 consist of solving approximatively the Dirac equation
with multiple right-hand sides DV = P and DV = R by using the single precision block
BiCGStab, which is described in Algorithm 1. We can quickly obtain an approximative
solution to the Wilson-Dirac equation by using the block BiCGStab algorithm.

10Performance Efficiency [%]

80
70
60
50

3

Bl-BiCGStab (12 x24)
DEE
AEE
DMSAP
3
Bl-BiCGStab (24 x48)
DEE
AEE
DMSAP

40
30
20
0

2

4

8

6

10

12

14

L
FIG. 1: Performance eﬃciencies of DEE , AEE , DMSAP kernels and the block BiCGStab algorithm
as a function of L with 123 ×24 (closed) and 243 ×48 (open) lattices on 16 nodes of the K computer.

Performance Efficiency [%]

90
80
70
60

4

Bl-BiCGStab (SAP bl. size 6 )
DEE
AEE
DMSAP
4
Bl-BiCGStab (SAP bl. size 12 )
DEE
AEE
DMSAP

50
40
30
20 1
10

10

2

10

3

4

10

# of nodes
FIG. 2: Performance eﬃciencies of DEE , AEE , DMSAP kernels and the block BiCGStab algorithm
as a function of the number of nodes for two diﬀerent SAP block size 64 (closed) and 124 (open).
The number of right-hand sides is ﬁxed to L = 12.

11

12

4 instr.

2/3 instr.

1 instr.

Calculated for L = 4

2048

964b

Calculated for L = 1

2048

964a

b

2048

964

a

256

483 × 96

2048

484

16

256

243 × 48(b) 243 × 48

16

31

27

32

32

33

30

33

35

85.0

76.0

84.9

85.1

85.1

83.6

80.7

84.6

15.6

12.4

12.0

13.0

13.0

17.1

9.7

13.9

21.0

19.5

21.1

21.3

21.5

22.1

24.5

24.4

16.6

15.5

17.9

17.7

17.7

18.1

20.1

19.2

of nodes [%] rate[%] commit[%] commit[%] commit[%](a) 123 × 24

size

Case Lattice number Eﬃc. SIMD

6.2

8.3

4.7

4.6

4.7

6.2

7.2

7.6

wait[%]

10.7

10.3

11.0

11.1

11.2

8.9

9.7

10.5

17.6

13.2

17.8

18.1

17.9

11.2

11.9

13.0

3.0

5.3

1.8

1.3

1.0

7.2

4.3

2.7

LA wait[%] LA wait[%] LA wait[%]

4.2

8.6

7.0

7.2

7.2

1.7

1.8

1.9

LA wait[%]Barrier sync. Fl. pt. op. Fl. pt. cache Int. cache Fl.pt. memory

the number L = 12 of right-hand sides, but for 964 lattice we also include the results for L = 1 and L = 4 for the sake of comparison.

K computer. The SAP block sizes are 64 in case (a) and 124 in case (b). LA stands for load access. These tests have been carried out withTABLE II: Results for the performance tests of the block BiCGStab kernel with diﬀerent lattice sizes on diﬀerent numbers of nodes of theAlgorithm 2 Outer BiCGStab algorithm
1: for i = 0 to L − 1 do
2:initial guess x(i) ∈ CN ,3:compute r(i) = b(i) − Dx(i) ,4:

set p(i) = r(i) ,5:choose r̃(i) such that (r̃(i) , r(i) ) ̸= 0,

6: end for
7: while maxi (|r (i) |/|b(i) |) > ϵouter do
8:solve D[v (1) , . . . , v (L) ] = [p(1) , . . . , p(L) ] by the block BiCGStab

9:for i = 0 to L − 1 do

10:α(i) = (r̃(i) , r(i) )/(r̃(i) , Dv (i) ),11:x(i) ← x(i) + α(i) v (i) ,12:r(i) ← r(i) − α(i) Dv (i) ,13:

end for14:

solve D[v (1) , . . . , v (L) ] = [r(1) , . . . , r(L) ] by the block BiCGStab

15:

for i = 0 to L − 1 do

16:

ζ (i) = (Dv (i) , r(i) )/(Dv (i) , Dv (i) ),17:x(i) ← x(i) + ζ (i) v (i) ,18:r(i) ← r(i) − ζ (i) Dv (i) ,19:

β (i) = (α(i) /ζ (i) ) · (r̃(i) , r(i) )/(r̃(i) , r′(i) ),20:p(i) ← r(i) + βk (p(i) − ζ (i) Dv (i) ),21:r′(i) = r(i)(i)

end for22:

23: end whileIV.

NUMERICAL TESTOur numerical tests are performed using 2048 nodes of the K computer.

We usea thermalized conﬁguration generated at almost the physical point with (κud , κs ) =(0.126117, 0.124790) on 964 lattice in 2 + 1 ﬂavor lattice QCD employing the nonpertur-

13

batively O(a)-improved Wilson quark action with cSW = 1.11 and the Iwasaki gauge action at β = 1.82. Stout smearing procedure is performed with the smearing parameters
Nstout = 6 and α = 0.1. We choose the hopping parameter κ = 0.126125 for the WilsonDirac equation, and NSAP = 5 with 124 SAP block size following Ref. [1]. Parameters
for the SSOR method are ﬁxed with NSSOR = 1 and ω = 1.24. The stopping criterion
is set to be maxi (|r(i) |/|b(i) |) ≤ 10−14 for the outer BiCGStab, while the stopping criterion (|r(i) |/|b(i) |) ≤ ϵinner(i = 1, 2, ..., L) for the inner block BiCGStab is set by the outeri≈ 10−7 to 10−3 . We have used
BiCGStab from iteration to iteration ranging from ϵinneriso-called “noise sources” generated with random numbers for multiple right-hand sides.
Figure 3 presents a representative case for the residual norm as a function of the total
number of inner matrix-vector multiplications (IMVMs) for the block BiCGstab with L =1, 2, 3, 4, 6 and 12 in addition to the residual norm of the outer BiCGStab solver. The
residual norm of the inner block BiCGStab decreases well along with the iteration, and oncethe stopping condition is satisﬁed, the algorithm is switched from the inner BiCGStab tothe outer one and returns again to the inner one in order to resume the next preconditioning
step with the residual norm initialized. Although for some cases the inner residual norm
decreases more quickly as L increases, the total number of iterations for convergence does notdiminish so drastically as observed in Refs. [6, 11], where the block BiCGStab is employed as
the outer solver. A possible reason is a rather loose stopping criterion compared to the oneused in Refs. [6, 11], where the required number of iterations seems to decrease much more
conspicuously as a function of L with a more stringent stopping criterion. Table III shows
the L dependence of the computational cost to solve the Wilson-Dirac equation as well as
the maximum and minimum values of the true residuals. We present the averaged values
over the (12/L) samples of L noise source vectors for the number of inner matrix-vector
multiplications and the execution time. The computational cost divided by L is found to
decrease as L increases. We do not have a problem of the convergence failure. In Table III
we also show the performance and eﬃciency of the block BiCGStab solver. It turns outthat we obtain more than 30% eﬃciency with L ≥ 2 and the performance goes beyond
85TFLOPS for several values of L. This implies that we have a clear gain in using the block
BiCGStab instead of its non-blocked version, where the maximum eﬃciency is about 26%.

14

10

2

10
10

-2

Residual

10
10

-4
-6

10
10
10

0

L=1
L=2
L=3
L=4
L=6
L=12

-8

-10

-12

10

-14

0

κ=0.126125

200

400

600 800#IMVM

1000 1200

FIG. 3: Respresentative behavior of the residual norm as a function of the number of inner matrixvector multiplications (IMVMs), together with the residual norm of the outer solver. The numberof right-hand side vectors are chosen to be L = 1, 2, 3, 4, 6, and 12.
TABLE III: L dependence of the number of inner matrix-vector multiplications (IMVMs) and theexecution time to solve the Wilson-Dirac equation. True residual is evaluated after the relative
residual in the outer solver reaches the tolerance of 10−14 .
L

#IMVM

Time [s]Time/L [s]True residual

PerformanceEﬃciency

Max

Min[TFLOPS]

[%]

1

915.3

65.2

65.2

4.99 × 10−15

7.24 × 10−16

71

27

2

913.0

128.4

64.2

2.69 × 10−15

4.59 × 10−16

81

31

3

876.8

193.8

64.6

4.49 × 10−15

4.20 × 10−16

86

33

4

869.7

235.6

58.9

2.21 × 10−15

4.09 × 10−16

82

31

6

856.5

349.3

58.2

4.62 × 10−15

4.48 × 10−16

85

32

12

770.0

669.0

55.8

4.25 × 10−15

6.52 × 10−16

83

32

V.

SUMMARY

In this work we have developed the single precision block BiCGStab solver for treating
the O(a)-improved Wilson-Dirac equation with multiple right-hand sides in lattice QCD.
The code development is carried out for the K computer and the optimization for the15

computation kernels shows the improvement in performance depending on the SAP block
size, the number of compute nodes, and the number of right-hand sides L. The performanceof the block BiCGStab solver is found to be more than 30% eﬃciency. With the block
BiCGStab solver used as the preconditioning step for the outer BiCGStab solver in double
precision, we have succeeded in test calculations with 964 lattice size on 2048 nodes of the K
computer. The computational cost reduces with increasing number of L. It may be possible
to make a further improvement for the QR decomposition in future.

Acknowledgments

All the results are obtained by using the K computer at the RIKEN Advanced Institutefor Computational Science.

[1] S. Aoki, K.-I. Ishikawa, N. Ishizuka, T. Izubuchi, D. Kadoh, K. Kanaya, Y. Kuramashi,
Y. Namekawa, M. Okawa, Y. Taniguchi, et al. (PACS-CS Collaboration), Phys. Rev. D 79,
034503 (2009).
[2] Y. Kuramashi, PoS Lattice 2008, 018 (2008).
[3] S. Aoki, K.-I. Ishikawa, N. Ishizuka, T. Izubuchi, D. Kadoh, K. Kanaya, Y. Kuramashi,
Y. Namekawa, M. Okawa, Y. Taniguchi, et al. (PACS-CS Collaboration), Phys. Rev. D 81,
074503 (2010).
[4] M. Lüscher, JHEP 0305, 052 (2003).
[5] J. A. Vogel, App. Math. Comput. 188, 226 (2007).
[6] H. Tadano and T. Sakurai, LSSC’07, Lec. Notes Comput. Sci. 4818, 721 (2008).
[7] D. O’Leary, Lenear Algebra Appl. 29, 293 (1980).
[8] A. Nikishin and A. Yu. Yeremin, SIAM J. Matrix Anal. Appl. 16, 1135 (1995).
[9] A. E. Guennouni, K. Jbilou, and H. Sadok, Electron. Trans. Numer. Anal. 16, 129 (2003).
[10] T. Sakurai, H. Tadano, and Y. Kuramashi, Comput. Phys. Commun. 181, 113 (2010).
[11] Y. Nakamura, K.-I. Ishikawa, Y. Kuramashi, T. Sakurai, and H. Tadano, Comput. Phys.
Commun. 183, 34 (2012).
[12] S. Fischer, A. Frommer, U. Glässner, T. Lippert, G. Ritzenhöfer, and K. Schilling, Comput.

16

Phys. Commun. 98, 20 (1996).
[13] T. Boku, K.-I. Ishikawa, Y. Kuramashi, K. Minami, Y. Nakamura, F. Shoji, D. Takahashi,
M. Terai, A. Ukawa, and T. Yoshié, PoS Lattice 2012, 188 (2012).

17