Verified Solutions of Large Sparse Linear
Systems Arising from 3D Poisson
Equation
OGITA, Takeshi
Waseda University
joint work withNAKAJIMA, Kengo, The University of Tokyo

45th ASE Seminar
November 29, 2023

Outline• Develop verified solutions of linear systems
discretized from 3D Poisson equation
Frequently used for benchmarking in HPC

𝜆!
~cond𝜆"(Dirichlet boundary condition)
𝜆 : thermal conductivity

𝜆! is considered asthermal insulation.

2What we focus on
PDE (3D Poisson equation)
Discretization by FDM, FEM, etc.
Discretization error
Rounding error

Linear system Ax = bIterative solver, such as
CG, GMRES, etc.
Truncation error
Rounding error

In this talk,we consider this part.

Approximate solution3Condition• 7-point FDM mesh
• Discretized by finite volume method• Matrices have M-property. (𝐴"# ≥ 𝑂)
• We can apply fast verification methods [1,2] forlinear systems with M-matrices.
[1] T. Ogita, S. Oishi, Y. Ushiro: Fast verification of solutions for sparsemonotone matrix equations, Computing, Supplement 15 (2001), 175-187.
[2] A. Minamihata, K. Sekine, T. Ogita, S. M. Rump, S. Oishi: Improved error
bounds for linear systems with H-matrices, Nonlinear Theory and Its
Applications, IEICE, 6:3 (2015), 377-382.

4

Verification algorithm (1)1. Solve a discretized linear system 𝐴𝑥 = 𝑏.
Ø 𝑥:
' a computed solution

2. Solve a linear system 𝐴𝑦 = 𝑒 where all
elements of 𝑒 are 1's.
Ø 𝑦:
' a computed solution

3. Verify M-property of 𝐴 using 𝑦.
) (𝑦' > 0 ⇒ 𝐴𝑦' > 0)
4. Compute an error bound using𝑦) $ 𝑏 − 𝐴𝑥) $𝑥 − 𝑥) $ ≤1 − 𝑒 − 𝐴𝑦) $if 𝑒 − 𝐴𝑦) $ < 1.
𝑦' $𝐴#! $ ≤1 − 𝑒 − 𝐴𝑦' $5

Numerical results (1)• Computer: Reedbush-U (1 node)– Intel Xeon E5-2695v4 (2.1GHz, 18 cores) x 2 sockets– 1.21 TFLOP/s per socket– 256 GiB (153.6GB/s)• Solver: ICCG with CM-RCM, MC(20)• Stopping criteria:ØFor 𝐴𝑥 = 𝑏,%"&(' !
% !

< 10"#)

ØFor 𝐴𝑦 = 𝑒, 𝑒 − 𝐴𝑦) $ < 10")
• FP64 (double precision), OpenMP (36 threads)
6

Result (1-1): 𝜆! = 𝜆" = 1.0

NX=NY=NZ=128 (n = 2,097,152)• Upper bounds of maximum relative error andrelative residual norm:– max–(" "('"≤ 3.38×10"-#*+*,("%"&(' !

< 3.66×10"##% !

• Computing time
Approximation Verification(Solve 𝐴𝑥 = 𝑏) (Solve 𝐴𝑦 = 𝑒)Elapsed time (sec.)

3.75

1.60

#iterations for ICCG

415

176
7

Result (1-2): NX=NY=NZ=128
Vary 𝜆!/𝜆"~cond between 1 and 106
Error bound106

better103

𝑥& − 𝑥'&max!%&%'𝑥&Residual norm𝑏 − 𝐴𝑥' "𝑏 "It is difficult to estimate the error of a computed solution
only from residual norm!

8Improved error bound• To reduce the overestimation, we replace𝑦) $ 𝑏 − 𝐴𝑥) $𝑥 − 𝑥) $ ≤1 − 𝑒 − 𝐴𝑦) $by𝑦) $ 𝑏 − 𝐴 𝑥) + 𝑧̂ $𝑥 − 𝑥) $ ≤ 𝑧̂ $ +.
1 − 𝑒 − 𝐴𝑦) $𝑥 − 𝑥) = 𝐴"# 𝑏 − 𝐴𝑥) = 𝑧̂ + 𝐴"# 𝑏 − 𝐴(𝑥) + 𝑧̂ )• The correction term 𝑧̂ is obtained by solving alinear system 𝐴𝑧 = 𝑟 with 𝑟 = 𝑏 − 𝐴𝑥.
)[3] T. Ogita, S. Oishi, Y. Ushiro: Fast inclusion and residual iteration forsolutions of matrix equations, Computing, Supplement 16 (2002), 171-184.

9

Verification algorithm (2)1. Solve a discretized linear system 𝐴𝑥 = 𝑏.
2. Solve a linear system 𝐴𝑦 = 𝑒.
3. Verify M-property of 𝐴 using 𝑦.
) (𝑦' > 0 ⇒ 𝐴𝑦' > 0)
4. Compute 𝑟 = 𝑏 − 𝐴𝑥) with an error bound.
Ø 𝑟:̂ a computed residual, 𝑒( : an error bound of 𝑟̂

5. Solve a linear system 𝐴𝑧 = 𝑟.̂
6. Compute an error bound using𝑦) $ ( 𝑟̂ − 𝐴𝑧̂ $ + 𝑒4 $ )𝑥 − 𝑥) $ ≤ 𝑧̂ $ +.
1 − 𝑒 − 𝐴𝑦) $[4] T. Ogita, S. M. Rump, S. Oishi: Accurate sum and dot product,
SIAM Journal on Scientific Computing, 26:6 (2005), 1955-1988.

10

Numerical results (2)• Computer: Reedbush-U (1 node)– Intel Xeon E5-2695v4 (2.1GHz, 18 cores) x 2 sockets– 1.21 TFLOP/s per socket, 256 GiB (153.6GB/s)• Solver: ICCG with CM-RCM, MC(20)• Stopping criteria:ØFor 𝐴𝑥 = 𝑏,%"&(' !
% !

< 10"#)

ØFor 𝐴𝑦 = 𝑒, 𝑒 − 𝐴𝑦) $ < 10")ØFor 𝐴𝑧 = 𝑟,̂

̂
4"&
6̂ !
4̂ !

< 10"7• FP64 (double precision), OpenMP (36 threads)
11

Result (2-1): NX=NY=NZ=128
Vary 𝜆!/𝜆"~cond between 1 and 106
1.E+00Previous boundImproved error bound
1.E-02New error boundMaximum relative error
Relative residual norm

1.E-0410-5

1.E-06𝑥& − 𝑥'&max!%&%'𝑥&

1.E-08better

1.E-1010-4

1.E-12
1.E+0 1.E+1 1.E+2 1.E+3 1.E+4 1.E+5 1.E+6
lambda1 / lambda2

Residual norm𝑏 − 𝐴𝑥' "𝑏 "Computed error bounds are significantly improved!
12

Result (2-2): Computing time• Elapsed time (sec.) and #iterations for ICCGNX=NY=NZ=128 Approximation
Verification(n = 2,097,152)(Solve 𝐴𝑥 = 𝑏) (Solve 𝐴𝑦 = 𝑒 and 𝐴𝒛 = 𝑟)̂
𝜆" /𝜆! = 1
𝜆" /𝜆! = 10#

3.75(415)
6.83(686)

4.47(493 = 176 + 317)
7.85(777 = 328 + 449)

NX=NY=NZ=256 Approximation
Verification(n = 16,777,216) (Solve 𝐴𝑥 = 𝑏) (Solve 𝐴𝑦 = 𝑒 and 𝐴𝒛 = 𝑟)̂
𝜆" /𝜆! = 1
𝜆" /𝜆! = 10#

76.64
(903)
110.15(1377)

80.89(964 = 387 + 577)
120.07(1504 = 639 + 865)
13

Future work1) Accelerate verification process:––

Apply mixed-precision approach
Optimization with double right-hand side2) Consider the discretization error3) Develop efficient verification algorithms forother matrix classes:–For example, positive definite matrices

14