Takeshi Iwashita(Kyoto Univ. / Hokkaido Univ.)

1Iwashita lab.

• Mixed precision computing

1. Development of mixed precision GMRES(m) solver• The solver is based on the iterative refinement method. An m-step
FP32 GMRES solver is used for the inner iteration. In the outer
iteration, the residual vector is calculated using FP64 and theconvergence of the relative residual norm is checked. Thus theaccuracy of the solution is the same as that of a conventional FP64
GMRES solver.

• Yingqi Zhao, Takeshi Fukaya, Linjie Zhang, Takeshi Iwashita:
Numerical Investigation into the Mixed Precision GMRES(m) Method Using
FP64 and FP32. J. Inf. Process. 30: 525-537 (2022)

2. Investigation of mixed precision IR-based Bi-CGSTAB solver• Methods based on the Lanczos principle look unsuitable for the innersolver of the iterative refinement method. However, we implemented
it and got some good results.
• To be published

2Iwashita lab.

3. Introduction of mixed precision computing to H-matrices• A part of low-rank submatrices are expressed using lower precisiondata. The proposed method was evaluated in a BEM analysis forelectric field.
• Rise Ooi, Takeshi Iwashita, Takeshi Fukaya, Akihiro Ida, Rio Yokota:
Effect of Mixed Precision Computing on H-Matrix Vector Multiplication in
BEM Analysis. HPC Asia 2020: 92-101.

4. Integer arithmetic based GMRES solver• An integer arithmetic based preconditioned GMRES solver wasdeveloped. The solver is based on the iterative refinement method
and only integer arithmetic are used in the inner iteration.
• Takeshi Iwashita, Kengo Suzuki, Takeshi Fukaya: An Integer
Arithmetic-Based Sparse Linear Solver Using a GMRES Method andIterative Refinement. ScalA@SC 2020: 1-8.

3Iwashita lab.

• ILU-type preconditioning techniques effectively using
SIMD vectorization
1. Proposal of hierarchical block multi-color ordering• The proposed parallel ordering method attains the same
convergence rate as the block multi-color ordering method, while
the forward and backward substitutions are efficiently parallelizedand vectorized. (The preconditioning matrix involves diagonalmatrices.)• Takeshi Iwashita, Senxi Li, Takeshi Fukaya: Hierarchical block multicolor ordering: a new parallel ordering method for vectorization andparallelization of the sparse triangular solver in the ICCG method. CCF
Trans. High Perform. Comput. 2(2): 84-97 (2020)

4Iwashita lab.

2. Proposal of ILUB preconditioning• A new fill-in strategy was proposed for ILU preconditioning. Nonzero
blocks are defined and all fill-ins in a nonzero block are permitted inILU factorization. Then, the forward and backward substitutions
consist of small dense matrix operations which is efficiently
vectorized.
• Kengo Suzuki, Takeshi Fukaya, Takeshi Iwashita: A novel ILU
preconditioning method with a block structure suitable for SIMD
vectorization. J. Comput. Appl. Math. 419: 114687 (2023)

5Iwashita lab.

•GPU oriented preconditioners

1. Proposal of a new AINV preconditioner• A modification of the AIV algorithm was presented. A newapproximation method for the preconditioner matrix greatly reduces
the preconditioner setup time with preserving the preconditioning
effect.
• Kengo Suzuki, Takeshi Fukaya, Takeshi Iwashita: A New AINVPreconditioner for the CG Method in Hybrid CPU-GPU Computing
Environment. J. Inf. Process. 30: 755-765 (2022)6Iwashita lab.

Aim of the research: Speedup of the convergence of asparse linear iterative solver• Subspace correction method•• J. Xu, “Iterative Methods by Space Decomposition and Subspace Correction”, SIAM Rev., Vol.
34, 1992.

•Deflation method• R. A. Nicolaides, “Deflation of conjugate gradients with application to boundary valueproblems”, SIAM J. Numer. Anal., Vol. 24, 1987.
• Y. Saad et al., “A deflated version of the conjugate gradient algorithm”, SIAM J. Sci. Comput.,Vol. 21, 2000.

•These techniques can collaborate with standardpreconditioning techniques and use a specific subspace.

• It is important to set an appropriate subspace.
• Here we consider a linear system of equations having a S.P.D. matrix. In many
practical simulations, a few isolated eigenvalues with very small magnitudes
cause a convergence problem. If we identify the eigenspaces corresponding tothese eigenvalues, we can accelerate the convergence of the iterative (CG)solver.
7

Iwashita lab.

•The way to obtain (approximate) eigenvectors withsmall eigenvalues• Based on knowledge of problems• Multigrid method, finite elements with high aspect ratio•

In this research, we focus on algebraic auxiliary matrixconstruction.
• Reduction of programming cost
• (Hidden) unapparent characteristics of a problem may be exploited.

•

We focus on a situation in which a series of linear
systems with an identical or similar coefficient matrix
are solved.
• In practical simulations, these linear systems often arise.

• Time dependent calculations, non-linear analysis (Newton method)

8Iwashita lab.

•Problem:
The right hand side vector bk
are functions of x1, x2,…,xk-1.

Proposed method, the following 4 steps• In the first solution step forare performed. For simple explanation, letanddenoted by and .
1. Sampling of approximation vectors
2. Calculation of error vectors
3. Rayleigh-Ritz method
4. Construction of auxiliary matrix•

9Iwashita lab.

❶Sampling of approximation vectors◦ Number of samples:◦ Sampling is done with a certain interval◦ ObtainSampling𝟏𝟐𝟑❷ Calculation of error vectors◦ After the solution process completes,the error vectors corresponding to thesampled approximation vectors arecalculated:.
◦ Apply the Gram-Schmidt process to thevectors and get orthonormal basis
vectors:
10𝟒(𝑠 = 1, … , 𝑚)

Iwashita lab.

❸ Rayleigh-Ritz method◦ Solve an eigenvalue problem:◦ Select Ritz vectors with small Ritz values less thanshould be much less than 1 when the coefficient matrix
is diagonally (or properly) scaled.
 : the number of selected Ritz vectors
 Selected Ritz vectors:❹ Auxiliary matrix◦

11

Iwashita lab.



We used a node of Fujitsu CX2550 (2 Xeon Skylake processors)
We tested three solvers◦ ICCG◦ ES-SC-ICCG: additive Schwarz SC and IC preconditioned CG solver◦ ES-D-ICCG: deflated ICCG solverWe also tested parallelized versions of three solvers

Using theproposedtechnique◦ Block Jacobi IC method is used for the parallelization of IC preconditioningThe convergence criterion: relative residual norm less than
Test problems: 30 matrices from SuiteSparse Matrix collection◦ We selected 30 largest datasets of symmetric positive-definite matrices.



We solved a linear system 6 times. In / after the first solution step,
vector sampling and the auxiliary matrix construction are performed.
Two types of right-hand side vectors are used: (1) a vector of onesand (2) a random vector.
The number of sampling vectors, : 20
The threshold for Ritz vectors, :,,The number of threads in tests for parallel solvers: 40
12Iwashita lab.

Results of tests using random vectors (dataset: Hook_1498)

Convergence acceleration

ES-SC-ICCG
Convergence acceleration

ES-D-ICCG

Significant convergence acceleration was achieved.
When we use a larger , which typically increases the number of selected Ritz
vector, , reduction in the number of iterations is enlarged.
Although the mechanism of convergence acceleration is different between
subspace correction and deflation, the convergence behaviors are mostly identical.
(This phenomenon was observed in all datasets.) It implies that the difference ofboth methods are marginally when the used subspace is identical.
13

Iwashita lab.

Results of tests using random vectors

The convergence rate was doubled for 20 out of 30 datasets.
◦When a vector of ones was used, the improvement in convergence rate was more significant.

For 12 datasets, the best result was obtained when(# sampled vectors).
For these datasets, the increase inpossibly improve the performance.

14

Iwashita lab.

Results of tests using random vectors

The time to solution was reduced for 28 out of 30 datasets using the proposedtechnique.
For 12 datasets, the best result was obtained when(# sampled vectors).
For these datasets, the increase inpossibly improve the performance.

15

Iwashita lab.

••••

We newly introduce an algebraic auxiliary matrix (subspace)construction method for subspace correction and deflation
based on error vector sampling in a series of linear systems.
In the method, limited number of approximation vectors arepreserved. The error vectors are calculated after the solution
step finishes. The auxiliary matrix is constructed using theseerror vectors.
Numerical tests using matrices from SuiteSparse confirm
that the subspace correction preconditioning and deflation
using the proposed technique accelerate the convergenceand the solution of a preconditioned CG solver.
Numerical test also indicates that the proposed method canbe used for the condition number estimation.

16Iwashita lab.

Basic concept◦ T. Iwashita, S. Kawaguchi, T. Mifune, and T. Matsuo, Automatic mapping
operator construction for subspace correction method to solve a series of linear
systems, JSIAM Letters, vol. 9, pp. 25-28, 2017.

Preliminary analysis on computational electromagnetics◦ T. Iwashita, S. Kawaguchi, T. Mifune, and T. Matsuo, Acceleration of Transient
Non-Linear Electromagnetic Field Analyses Using an Automated Subspace
Correction Method, IEEE Trans. Magn., vol. 55, pp. 1-4, 2019.

Regular paper◦ T. Iwashita, K. Ikehara, T. Fukaya, and T. Mifune, Convergence acceleration ofpreconditioned conjugate gradient solver based on error vector sampling for asequence of linear systems, Numer. Linear Algebra Appl. (2023), e2512.

17

Iwashita lab.

We will welcome your visit to Hokkaido & Kyoto Univ.

18

Iwashita lab.