北海道大学

Performance Evaluation of the Mixed Precision GMRES(m)Method using FP64 and FP32
The international Conference on High Performance Computing in Asia-Pacific Region(HPC Asia 2022)
January 12, 2022 | online

Presenter: Yingqi Zhao (Hokkaido University)
Collaborator：Takeshi Fukaya (Hokkaido University, JST PRESTO)
Takeshi Iwashita (Hokkaido University)

Overview◆ Background➢ In scientific computations, exploiting the high ability of low precision computing is now getting crucial.
➢ In order to accelerate linear solvers, it is necessary to developing a new method that can exploit low
precision computing and provide computed results with the same accuracy as by traditional methods
using only FP64.

◆ Main study➢ Target problem : 𝐴𝑥 = 𝑏, 𝐴 ∈ ℝ𝑛×𝑛 , 𝑥, 𝑏 ∈ ℝ𝑛 , where 𝐴 is large, sparse and non symmetric.
➢ Mixed Precision GMRES(𝑚) (called MP-GMRES(𝑚)) using FP64 and FP32 was investigated through thenumerical experiments and compared with GMRES(𝑚) using only FP64(called FP64-GMRES(𝑚)).
➢ Numerical behaviors ( e.g. Number of the total iterations, achievable accuracy) were examined andevaluated.
2Outline◼ Introduction➢Background➢

Main study◼ Method➢

GMRES(m) method➢

Mixed precision GMRES(m) method◼ Numerical experiments➢

Experimental settings➢

Number of iterations for attaining the convergence criterion➢

Achievable accuracy within the limit of iterations➢

Execution time for attaining convergence criterion◼ Conclusion◼ References

3

Introduction

Background◆ Traditional scientific computing➢ FLOPS(Floating-point Operations Per Second) for Double-precision floating-point number (FP64)
operation is metric for performance of CPU or computer system.
➢ FP64 has been used in scientific computations as standard.

◆ Current situation➢ Due to several reasons(e.g. power budget), it is getting difficult to improve FP64 FLOPS.
➢ Recent hardware trend: Some new applications accept low precision, for example, AI applications.

◆ Advantages of using low precision computing➢ Reduce the cost of data transfer.

➢ Increase the arithmetic performance, e.g. wider SIMD operations.

Exploit the advantages of low precision computing➔ Key idea is mixed precision computing, in which both FP64 and low precision computing
4are combined.

Introduction

Main study◆ Target Problem➢ 𝐴𝑥 = 𝑏, 𝐴 ∈ ℝ𝑛×𝑛 , 𝑥, 𝑏 ∈ ℝ𝑛

𝐴 : coefficient matrix
large, sparse, and non symmetric◆ Method➢ Original algorithm: GMRES(𝑚) method (Saad & Schultz, 1986; Saad, 2003).
➢ Target algorithm: mixed precision GMRES(𝑚) method, a mixed precision variant of the GMRES(𝑚)method using FP64 and FP32 (single-precision floating-point number).

◆ Objective➢ Investigate the numerical behavior of mixed precision GMRES(𝑚) method using FP64 and FP32.

➢ Make a comparison between traditional GMRES(𝑚) method using only FP64 and mixed precision
GMRES(𝑚) method.
5

Method : Iterative refinement◆ Iterative refinement
The solution 𝑥 can be iteratively improved and reach the required accuracy by iterative refinement.
Let 𝑥෤ be the current approximate solution for 𝐴𝑥 = 𝑏,
Step 1 : Compute the residual 𝑟 = 𝑏 − 𝐴𝑥.
෤Step 2 : Solve the error equation 𝐴𝑒 = 𝑟 and obtain an approximate solution 𝑒ǁStep 3 : Update the solution 𝑥෤ = 𝑥෤ + 𝑒ǁ

◆ A mixed precision variant of the iterative refinement scheme➢ Solution 𝑥෤ is stored in high precision.
➢ Step 1 and 3 use high precision.
➢ Step 2 uses low precision.
6

Method : GMRES(𝒎) method
Algorithm 1 GMRES(𝑚)Input : 𝑥0 : initial guess, 𝜖: convergence criterion, 𝐴: coefficient matrix, 𝑏: right-hand side vector, maximum
number of iterations
1:repeat2:𝑟0 = 𝑏 − 𝐴𝑥03:

𝛽 = 𝑟0 2

Step 1: Solve 𝐴𝑥 = 𝑏 by 𝑚-iteration

4:

If 𝛽/ 𝑏 2 ≤ 𝜖 then return 𝑥0

GMRES with initial guess 𝑥0 and

5:𝑣0 = 𝑟0 /𝛽

obtain approximate solution 𝑥𝑚 .

6:

Step 2: Update the initial guess :𝑥0 =7:

෩𝑚 by the 𝑚-step Arnoldi process with 𝐴 and 𝑣0 .
Compute 𝑉𝑚 and 𝐻෩𝑚 .
Compute 𝑦𝑚 from 𝛽 and 𝐻

8:𝑥𝑚 = 𝑥0 + 𝑉𝑚 𝑦𝑚

9:𝑥0 = 𝑥𝑚

10:𝑥𝑚 .

until attain the maximum number of iterations

Output: 𝑥07

Method : Mixed precision GMRES(𝒎) method
In GMRES(𝑚) method
Step 1 : Solve 𝐴𝑥 = 𝑏 by 𝑚-iteration GMRES withinitial guess 𝑥0 and obtain approximate solution 𝑥𝑚 .
Step 2 : Update the initial guess :𝑥0 = 𝑥𝑚 .
Mathematically equivalent(Imakura, A. et al., 2012 )Step 1 : Solve 𝐴𝑒 = 𝑟 by 𝑚-iteration GMRES
with initial guess 𝑒0 = 0 and obtain approximate
solution 𝑒𝑚 , where 𝑟 = 𝑏 − 𝐴𝑥0 .
Step 2 : Update the initial guess : 𝑥0 = 𝑥0 + 𝑒𝑚 .

GMRES(𝑚) has the iterative refinement structure.
➔ Obtain a mixed precision variant of GMRES(𝑚) method based on iterative refinement(mixed precision GMRES(𝑚) method).

8

Method : Mixed precision GMRES(𝒎) method
Algorithm 2 Mixed Precision GMRES(𝑚)Input : 𝑥0 : initial guess, 𝜖: convergence criterion, 𝐴: coefficient matrix, 𝑏: right-hand side vector, maximum number ofiterations𝐿𝑜𝑤() : precision converts from standard to low
1:
𝐴(𝐿) = 𝐿𝑜𝑤 𝐴
2:repeat3:𝑟0 = 𝑏 − 𝐴𝑥0 ,𝛽 = 𝑟0 2

4:

If 𝛽/ 𝑏 2 ≤ 𝜖 then return 𝑥0

5:𝑣0 = 𝑟0 /𝛽6:𝛽 (𝐿) = 𝐿𝑜𝑤 (𝛽) , 𝑣0

7:

(𝐿)
(𝐿)
(𝐿)
෩𝑚
Compute 𝑉𝑚 and 𝐻
by the 𝑚-step Arnoldi process with 𝐴(𝐿) and 𝑣0 .

8:

෩𝑚 .
Compute 𝑦𝑚 from 𝛽 (𝐿) and 𝐻

9:𝑧𝑚 = 𝑉𝑚 𝑦𝑚

10:𝑧𝑚 = 𝑆𝑡𝑑(𝑧𝑚 )11:𝑥𝑚 = 𝑥0 + 𝑧𝑚 , 𝑥0 = 𝑥𝑚

12:

(𝐿)

= 𝐿𝑜𝑤 (𝑣0 )

(𝐿)

(𝐿)

(𝐿)

Using low precision arithmetic(𝐿) (𝐿)(𝐿)

𝑆𝑡𝑑() : precision converts from low to standard

until attain the maximum number of iterations

Output: 𝑥09

Numerical experiments: Experimental settings◆ Environment➢ CPU : Xeon Gold 6148 (20 cores, 2.4GHz)×2➢ Program : C language, OpenMP, CRS format for sparse matrix 𝐴➢ Compiler : “Icc” (ver. 19.1.3.304), “-xCORE-AVX512 –qopenmp”(40 threads)◆ Test matrices➢ Database : SuiteSparse Matrix Collection➢ Matrix kind : Real matrices appear in scientific simulation➢ Condition number 𝜅2 𝐴 (= 𝐴 2 / 𝐴−1 2 ) : Range from 𝑂(101 ) to 𝑂(108 )◆ Settings➢ 𝑏 = 1,1, … , 1 𝑇 , and initial guess 𝑥0 = 0,0, … , 0 𝑇➢ Convergence criterion : 𝑏 − 𝐴𝑥 2 / 𝑏 2 ≤ 10−10➢ Iteration limit : (total # of inner iterations) ≤ n➢ Candidates of 𝑚 : 50, 100, 200, 300, 400, 500

10Numerical experiments—— Number of iterations for attaining the convergence criterion
FP64-GMRES(m) attained within one restart𝜅2 (𝐴)

Matrix𝒎 = 𝟓𝟎𝒎 = 𝟏𝟎𝟎

𝒎 = 𝟐𝟎𝟎𝒎 = 𝟑𝟎𝟎𝒎 = 𝟒𝟎𝟎𝒎 = 𝟓𝟎𝟎

FP64MPFP64MPFP64MPFP64MPFP64MPFP64MP𝑂(101 ) cage10265926

110

26

211

26

3112641126

511

𝑂(102 ) Zhao1

437343

123

43223

43323

43423

43523

𝑂(103 ) epb1

1994 1881 1781 1786 1516 1517 1290 1302 1153 1153 1095 1272

𝑂(104 ) af23560

—

—

—

—

—

—

4799 4777 4305 4283 4442 4351

𝑂(105 ) memplus 5900 4521 3030 2943 1942 2044 1532 2404 1459 3265 1352 3108
𝑂(106 ) inlet————————————

𝑂(107 ) garon2—

—

—

—

—

—

—

—

—

—

—

—

𝑂(108 ) rajat15

—

—

—

—

—

—

—

—

—

—

—

—

Results marked in red are the minimum # of iterations, “—” means the method didn’t attain the convergence within limit of iterations

Analysis:➢ When 𝑚 is small, the number of iterations of the two methods are almost equal.
➢ When 𝑚 becomes larger, the number of iterations of FP64-GMRES(𝑚) basically decreases, but that ofMP-GMRES(𝑚) tends to increase.

11

Numerical experiments—— Case study : matrix 𝒎𝒆𝒎𝒑𝒍𝒖𝒔◆ The number of iterations:➢ FP64-GMRES(𝑚): monotonically decrease➢ MP-GMRES(𝑚): first decrease and thenincrease◆ MP-GMRES(𝑚) attained the convergence criterionwith smaller iterations than FP64-GMRES(𝑚) when𝑚 = 50, 100.
Fig.1 The change of the number of the total inner iterations in FP64GMRES(𝑚) and MP-GMRES(𝑚) for 𝑚𝑒𝑚𝑝𝑙𝑢𝑠.

12Numerical experiments—— Achievable accuracy within the limit of iterations
𝒎 = 𝟓𝟎𝒎 = 𝟏𝟎𝟎

𝒎 = 𝟐𝟎𝟎𝒎 = 𝟑𝟎𝟎𝒎 = 𝟒𝟎𝟎𝒎 = 𝟓𝟎𝟎

FP64MPFP64MPFP64MPFP64MPFP64MPFP64MP𝑂(101 ) cage10-11-11-11-11-11-11-11-11-11-11-11-11𝑂(102 ) Zhao1-11-11-11-11-11-11-11-11-11-11-11-11𝑂(103 ) epb1-11-11-11-11-11-11-11-11-11-11-11-11𝑂(104 ) af23560-1-1-1-1-1-1-11-11-11-11-11-11𝑂(105 ) memplus-11-11-11-11-11-11-11-11-11-11-11-11𝑂(106 ) inlet-1-1-1-1-1-1-1-1-1-1-1-1𝑂(107 ) garon2-1-1-1-1-2-2-3-3-3-3-3-3𝑂(108 ) rajat15-1-1-1-1-1-1-1-1-1-1-10

𝜅2 (𝐴)Matrix

log10 ( 𝑏 − 𝐴𝑥 2 / 𝑏 2 ) when attained the iteration limit or the convergence criterion(“-11” corresponds to this case)

Analysis:➢ In most cases, the final accuracy of both FP64-GMRES(𝑚) and MP-GMRES(𝑚) is “-11”, which meansthat both methods attained the convergence criterion𝑏−𝐴𝑥 2≤ 10−10 .
𝑏 2➢ Generally, if both two methods did not converge, there is almost no difference in the final accuracy.

13

Numerical experiments—— Execution time for attaining convergence criterion
𝜅2 (𝐴)
1𝑂(10 )

𝑂(102 )𝑂(103 )

𝑂(104 )

Matrixcage10Zhao1
epb1

af23560

𝑂(105 ) memplus

Method𝒎 = 𝟓𝟎𝒎 = 𝟏𝟎𝟎𝒎 = 𝟐𝟎𝟎𝒎 = 𝟑𝟎𝟎𝒎 = 𝟒𝟎𝟎𝒎 = 𝟓𝟎𝟎

FP641.03 × 10−2

1.22 × 10−2

1.10 × 10−2

1.25 × 10−2

9.49 × 10−3

1.14 × 10−2MP

1.43 × 10−2

3.40 × 10−2

9.85 × 10−2

2.01 × 10−1

4.09 × 10−1

5.43 × 10−1FP64

3.38 × 10−2

3.77 × 10−2

3.72 × 10−2

3.79 × 10−2

3.27 × 10−2

3.49 × 10−2MP

2.90 × 10−2

7.75 × 10−2

3.03 × 10−1

7.62 × 10−1

1.39 × 100

2.18 × 100

FP64

6.09 × 10−1

9.65 × 10−1

1.59 × 100

2.19 × 100

2.81 × 100

3.50 × 100

MP3.24 × 10−1

5.19 × 10−1

8.09 × 10−1

9.96 × 10−1

1.24 × 100

1.50 × 100

FP64

—

—

—

1.47 × 101

1.78 × 101

2.46 × 101

MP

—

—

—

5.89 × 100

7.35 × 100

9.91 × 100

FP64

2.32 × 100

2.12 × 100

2.69 × 100

3.49 × 100

4.36 × 100

5.16 × 100

MP1.30 × 100

1.23 × 100

1.30 × 100

2.18 × 100

3.37 × 100

4.14 × 100

Results marked in red are the shortest execution time, “—” means that the method did not attain the convergence criterion.

Analysis:➢ When 𝑚 is small, the execution time of both methods tends to be short.
➢ For many test matrices, MP-GMRES(𝑚) has shorter execution time than MP-GMRES(𝑚).

14

Numerical experiments—— Case study : matrix 𝒘𝒂𝒏𝒈𝟒
Even in the case that larger 𝑚 reduce the number of iterations in FP64-GMRES(𝑚),MP-GMRES(𝑚) is still faster.
𝜅2 (𝐴)

Matrix𝑂(105 ) wang4

ItemMethod 𝒎 = 𝟓𝟎 𝒎 = 𝟏𝟎𝟎 𝒎 = 𝟐𝟎𝟎𝒎 = 𝟑𝟎𝟎𝒎 = 𝟒𝟎𝟎𝒎 = 𝟓𝟎𝟎

ExecutiontimeFP64——

3.53 × 100

3.62 × 100

3.34 × 100

3.65 × 100

MP——1.69 × 100

2.00 × 100

2.85 × 100

4.55 × 100# ofiterations

FP64——

1741

1101

738

688

MP——

1896

1423

1477

1792

For matrix 𝑤𝑎𝑛𝑔4, about three times the iterations are required,MP-GMRES(𝑚) is still faster.
15

Numerical experiments—— Performance modeling◆ Our approach to performance modeling.
➢In order to confirm the acceleration effect of MP-GMRES(𝑚), speedup of MP-GMRES(𝑚)over FP64-GMRES(𝑚) is considered. Firstly, a model is established to analyze thespeedup theoretically.
➢GMRES(𝑚) is memory bounded computation and Arnoldi process is used. The 𝑚-step
Arnoldi process mainly consists of 𝑚 SpMVs and MGS for 𝑚 vectors.

➢The estimation of the memory access cost in the 𝑚-step Arnoldi process in FP64 andFP32 is summarized.
16

Numerical experiments—— Expected speedup◆ Model construction and expected speedup

SpMVsMGSTotal

FP64

MP12𝑁𝑛𝑧 + 20𝑛 𝑚
1 28∙ 𝑚 𝑛
24𝑚2 𝑛 + 12𝑚𝑁𝑛𝑧 + 20𝑚𝑛

8𝑁𝑛𝑧 + 12𝑛 𝑚
1 24∙ 𝑚 𝑛
22𝑚2 𝑛 + 8𝑚𝑁𝑛𝑧 + 12𝑚𝑛

Let 𝛾 be the ratio of the number of the iterations in MP-GMRES(𝑚) over that in FP64-GMRES(𝑚)The speedup of MP-GMRES(𝑚) over FP64-GMRES(𝑚) is estimated as
4𝑚2 𝑛 + 12𝑚𝑁𝑛𝑧 + 20𝑚𝑛
2𝑚 + 6𝑁𝑛𝑧 /𝑛 + 102𝑐 + 6𝑆𝑝𝑒𝑒𝑑𝑢𝑝 ==≃𝛾(2𝑚2 𝑛 + 8𝑚𝑁𝑛𝑧 + 12𝑚𝑛 ) 𝛾(𝑚 + 4𝑁𝑛𝑧 /𝑛 + 6)
𝛾 𝑐+4

where 𝑐 = 𝑚Τ(𝑁𝑛𝑧 Τ𝑛)17Numerical experiments—— Comparison between actual speedup and estimated speedup
For the cases in which the setting of 𝑚 that provides the shortest execution time is not different between

FP64-GMRES(𝑚) and MP-GMRES(𝑚), the actual speedup and the estimated speedup is compared.
Actual Estimatedspeedup speedup𝜅2 (𝐴)

Matrix𝒎𝑵𝒏𝒛 Τ𝒏𝒄𝜸

𝑂(102 )

ns3Da50

82

0.61

1.00

1.26

1.56

𝑂(103 )

epb150

6

8.33

0.94

1.88

1.94

𝑂(104 )

af23560

300

21

14.29

1.00

2.50

1.90

𝑂(105 )

memplus

100

7

14.29

1.02

1.73

1.86

➢ Since the estimated speedup is calculated based on the rough modeling, the estimated speedup

seems to be close to the actual speedup.
➢ However, as SpMV is more complicated than MGS, it might bring difficulties to estimation which
may cause a gap between actual and estimated speedup.

18

Conclusion◆ What we did➢ MP-GMRES(𝑚) using FP64 and FP32 was investigated through the numerical experiments and comparedwith FP64-GMRES(𝑚) using only FP64.
➢ The numerical behaviors(e.g. number of the total iterations and achievable accuracy) were examined.
➢ The execution time of the methods using 40 threads was also evaluated on a standard CPU platform.

◆ Main observations➢ If a problem is solved by FP64-GMRES(𝑚), basically, MP-GMRES(𝑚) can solve it.
➢ When 𝑚 is small, the number of iterations of the two methods are almost equal.

➢ As 𝑚 grows, the number of the total iterations often decreases in FP64-GMRES(𝑚) but tends to increasein MP-GMRES(𝑚).
➢ On the standard CPU platform, it can be expected that MP-GMRES(𝑚) outperforms FP64-GMRES(𝑚) interms of the execution time.

19

References[1] Abdelfattah, A., Anzt, H., Boman, E. G., Carson, E., Cojean, T., Dongarra, J., Fox, A., Gates, M., Higham, N. J., Li, X. S., Loe, J., Luszczek, P.,Pranesh, S., Rajamanickam, S., Ribizel, T., Smith, B. F., Swirydowicz, K., Thomas, S., Tomov, S., Tsai, Y. M. and Yang, U. M.: A survey ofnumerical linear algebra methods utilizing mixed-precision arithmetic, The International Journal of High Performance Computing
Applications, Vol. 35, No. 4, pp. 344–369 (2021).
[2] Saad, Y. and Schultz, M. H.: GMRES: A Generalized Minimal Residual Algorithm for Solving Nonsymmetric Linear Systems, SIAM Journal
on Scientic and Statistical Computing, Vol. 7, No. 3, pp. 856-869 (1986).
[3] Saad, Y.: Iterative Methods for Sparse Linear Systems, Society for Industrial and Applied Mathematics (2003).
[4] Imakura, A., Sogabe, T. and Zhang, S.-L.: An Effecient Variant of the GMRES(𝑚) Method Based on the Error Equations, East Asian
Journal on Applied Mathematics, Vol. 2, No. 1, pp. 19-32 (2012).

[5] Lindquist, N., Luszczek, P. and Dongarra, J.: Improving the Performance of the GMRES Method Using Mixed Precision Techniques,
Driving Scientific and Engineering Discoveries Through the Convergence of HPC, Big Data and AI, pp. 51–66 (2020).
[6] Lindquist, N., Luszczek, P. and Dongarra, J.: Accelerating Restarted GMRES with Mixed Precision Arithmetic, IEEE Transactions onParallel and Distributed Systems, Vol. 33, pp. 1027–1037 (2022).
[7] Demmel, J. W.: Applied Numerical Linear Algebra, Society for Industrial and Applied Mathematics (1997).
[8] Higham, N. J.: Accuracy and Stability of Numerical Algorithms, Society for Industrial and Applied Mathematics, second edition (2002).

20

THANKS