Performance enhancement of the Ozaki Scheme on integer matrix multiplication unit

Abstract

This study was aimed at simultaneously achieving sufficient accuracy and high performance for general matrix multiplications. Recent architectures, such as NVIDIA GPUs, feature high-performance units designed for low-precision matrix multiplications in machine learning models, and next-generation architectures are expected to follow the same design principle. The key to achieving superior performance is to fully leverage such architectures. The Ozaki scheme, a highly accurate matrix multiplication algorithm using error-free transformations, enables higher-precision matrix multiplication to be performed through multiple lower-precision matrix multiplications and higher-precision matrix additions. Ootomo et al. implemented the Ozaki scheme on high-performance matrix multiplication units with the aim of achieving both sufficient accuracy and high performance. This paper proposes alternative approaches to improving performance by reducing the numbers of lower-precision matrix multiplications and higher-precision matrix additions. Numerical experiments demonstrate the accuracy of the results and conduct performance benchmarks of the proposed approaches. These approaches are expected to yield more efficient results in next-generation architectures.

1. Introduction

The field of machine learning, including AI, is evolving daily, and the scale and complexity of machine learning models are continually increasing. Recent architectures are designed to process these models rapidly and with high energy efficiency. For machine learning models, matrix multiplications using low-precision floating-point systems and integers are essential. Therefore, recent architectures are equipped with high-performance low-precision floating-point and integer matrix multiplication units, and high-performance mixed-precision matrix multiplication units that leverage these capabilities. A prime example of this is the NVIDIA Tensor Cores (see NVIDIA Corporation (2024)). Table 1 shows the specifications of the NVIDIA GPUs equipped with Tensor Core technology. Note that the specification of FP16 TC on RTX 4090 (165 TFLOPS) is for FP16 input and FP32 output, the specifications for the H100 and GH200 are the same as those for the H200, and the specifications for the B200 in the table is based on NVIDIA’s Technical Brief. In the future, numerical computation algorithms that leverage the performance of cutting-edge architectures will be essential. This study focused on researching high-performance matrix multiplication algorithms that maximize the potential of the latest architectures.

Table 1.

Specifications in TFLOPS/TOPS of NVIDIA GPUs for dense data NVIDIA Corporation (2024).

	Ampere	ADA	Hopper	Blackwell
	A100 SXM4	RTX 4090	H200 SXM5	B200
	GA100	AD102-300
FP64	9.7	1.29	34	40
FP64 TC	19.5	–	67	40
FP32	19.5	82.6	67	80
TF32 TC	156	82.6	494	1100
BF16 TC	312	165	989	2250
FP16 TC	312	165	989	2250
INT8 TC	624	661	1979	4500
FP8 TC	–	661	1979	4500
FP6 TC	–	–	–	4500
INT4 TC	–	1321	–	–
FP4 TC	–	–	–	9000

A highly accurate matrix multiplication scheme via the error-free transformation of matrix products was proposed by Ozaki et al. (2012). The scheme is called the Ozaki scheme and it converts a matrix product into a sum of multiple matrix products. It is also possible to convert a matrix product into a sum of lower-precision matrix products to take full advantage of the immense computational power of recent architectures. In order to compute D ← AB for higher-precision matrices $A\in {\mathbb{R}}^{m\times n}$ and $B\in {\mathbb{R}}^{n\times p}$ , the Ozaki scheme using lower/mixed-precision, provided by Mukunoki et al. (2020), is constructed as the following four steps:

(i) Extract the lower-precision matrices A₁, A₂, …, A _k from the higher-precision matrix A, where k is specified by the user, shifting A _i to prevent overflow and underflow at lower-precision arithmetic.

(ii) Extract k lower-precision matrices B _i from B in the similar way as A _i .

(iii) Compute A _i B _j for i + j ≤ k + 1 using lower/mixed-precision arithmetic.

(iv) Reverse (shift) A _i B _j and accumulate the results into D using higher-precision arithmetic.

When emulating the GEMM routine, the following fifth step is also performed:

• (v) Compute C ← αD + βC for scalars α, β and a matrix C.

Ootomo et al. (2024) implemented the Ozaki scheme using the INT8 Tensor Cores and evaluated the performance for emulating DGEMM, an FP64 matrix multiplication routine. The Ozaki scheme implemented by Ootomo is named the “Ozaki Scheme on Integer Matrix Multiplication Unit” (ozIMMU), and the code is available from the GitHub repository by Ootomo (2024). Figure 1 shows the accuracy of ozIMMU. Herein, ϕ specifies the tendency of the difficulty in terms of the accuracy of matrix multiplications. At least 7 or 8 slices are required to obtain sufficiently accurate results, even for well-conditioned matrix multiplications. As ϕ increases, more slices are required to obtain sufficient accuracy. Ootomo’s implementation for splitting matrices offers bit masking. Thus, the extracted matrices A _i , B _i may not be optimal and the splitting method can be improved (see Section 3.1 for details). Improvement of the splitting method should contribute to improving the accuracy of results.

Figure 1.

Accuracy of ozIMMU. Matrix A has entries a _ij ≔(U _ij − 0.5) ⋅ exp(ϕ ⋅ N _ij ), where U _ij ∈ (0, 1) are uniformly distributed and N _ij are drawn from standard normal distribution for 1 ≤ i, j ≤ m and m = n = p. Matrix B is composed similarly.

Figures 2 and 3 show the time breakdown of ozIMMU for double-precision matrices $A\in {\mathbb{R}}^{m\times n}$ and $B\in {\mathbb{R}}^{n\times p}$ , defined by IEEE Computer Society (2019), on GeForce RTX 4090 and GH200 Grace Hopper Superchip, respectively. Note that “split A”, “split B”, “Gemm in INT8-TC”, “accumulation in FP64”, and “copy” correspond to the steps (i), (ii), (iii), (iv), and (v) in the Ozaki scheme for emulating the GEMM routine, respectively. It can be seen that the computing time for the accumulation of matrix products in FP64 is not negligible even though the computation cost is $\mathcal{O}(mp)$ operations. This is because the performance of the INT8 Tensor Core is nearly 512 times and 60 times higher than that of FP64 on the RTX 4090 and GH200, respectively. Because this ratio increases significantly on the B200 and future architectures, accelerating the accumulation process is critical for the Ozaki scheme.

Figure 2.

Time breakdown of ozIMMU on NVIDIA GeForce RTX 4090.

Figure 3.

Time breakdown of ozIMMU on NVIDIA GH200 Grace Hopper Superchip.

In this paper, we propose new implementation methods for accelerating the accumulation of the Ozaki scheme using the INT8 Tensor Core. In addition, an alternative splitting method is also applied to improve the accuracy of results. For the number of slices, the splitting method contributes to obtaining more accurate results than those of Ootomo’s ozIMMU. The remainder of this paper is organized as follows: Section 2 overviews previous studies about the Ozaki scheme using an integer matrix multiplication unit; Section 3 presents the proposed method for accelerating the accumulation and optimizing the splitting method; Section 4 shows numerical results to illustrate the efficiency of the proposed methods; Section 5 provides a rounding error analysis of the Ozaki scheme with the proposed methods; and Section 6 presents final remarks.

2. Previous study

In this section, we briefly summarize previous studies. Let u be the relative error unit, e.g., u = 2⁻⁵³ for double-precision floating-point numbers defined in IEEE 754. Define $\mathbb{F}$ as a set of binary floating-point numbers with u. Let O_m,n be an m × n zero matrix. Suppose that $A\in {\mathbb{F}}^{m\times n}$ and $B\in {\mathbb{F}}^{n\times p}$ . Ozaki et al. (2012) proposed an error-free transformation of a matrix product AB. For a user-specified constant $k\in \mathbb{N}$ , the Ozaki scheme transforms each of A and B into k + 1 matrices

$\begin{array}{l}\hfill A≕A_1+A_2+\cdots +A_k+V_k,\\ \hfill B≕B_1+B_2+\cdots +B_k+W_k,\end{array}$

where $A_i\in {\mathbb{F}}^{m\times n}$ , $B_i\in {\mathbb{F}}^{n\times p}$ ,

$V_i≔A-\sum _{j=1}^i{A}_j,\text{and}{W}_i≔B-\sum _{j=1}^i{B}_j,$

(1)

such that $V_i\in {\mathbb{F}}^{m\times n}$ and $W_i\in {\mathbb{F}}^{n\times p}$ for i ≤ k, and $|{(A_s)}_{ij}|\ge |{(A_{s+1})}_{ij}|$ holds if ${(A_s)}_{ij}\ne 0$ for s = 1, …, k − 2, and B _s satisfies the corresponding relation. If the technique provided by Minamihata et al. (2016) is used, then $|{(A_s)}_{ij}|>|{(A_{s+1})}_{ij}|$ holds when ${(A_s)}_{ij}\ne 0$ . Note that ${(A_s)}_{ij}$ and ${(B_s)}_{ij}$ are obtained depending on ${\max }_{1\le j\le n}|{(V_{s-1})}_{ij}|$ and ${\max }_{1\le i\le m}|{(W_{s-1})}_{ij}|$ , respectively.

Let ${\mathbb{F}}_N$ be a set of N-bit binary floating-point numbers. For the Ozaki scheme using the FP16 Tensor Core with accumulation in FP32 and outputs in FP32, the mixed-precision splitting method with shifting to prevent overflow and underflow was utilized by Mukunoki et al. (2020). For $A\in {\mathbb{F}}_{64}^{m\times n}$ and $B\in {\mathbb{F}}_{64}^{n\times p}$ , the splitting method produces

$\begin{array}{l}\hfill A≕\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}({\mu }^{\prime (1)})A_1^{\prime }+\cdots +\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}({\mu }^{\prime (k)})A_k^{\prime }+V_k,\end{array}$

(2)

$\begin{array}{l}B≕B_1^{\prime }\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}({\nu }^{\prime (1)})+\cdots +B_k^{\prime }\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}({\nu }^{\prime (k)})+W_k\hfill \end{array}$

(3)

with double-precision shift values ${\mu }^{\prime (i)}\in {\mathbb{F}}_{64}^m$ , ${\nu }^{\prime (i)}\in {\mathbb{F}}_{64}^p$ and half-precision matrices $A_i^{\prime }\in {\mathbb{F}}_{16}^{m\times n}$ , $B_i^{\prime }\in {\mathbb{F}}_{16}^{n\times p}$ for i = 1, …, k. Algorithm 1 represents the splitting method for obtaining (2) using Minamihata’s technique. Equation (3) can be obtained by transposing the results of Algorithm 1 executed for B ^T . Algorithm 1 can be described as a loop that is executed until s = k for simplicity; however, the loop terminates when A = O_m,n in practice. Note that the binary logarithm is used at the 3rd line in Algorithm 1; however, the calculation using the binary logarithm occasionally returns erroneous results. Therefore, it is better to use a calculation method without the binary logarithm, such as a bitwise operation or a technique leveraging rounding error in floating-point arithmetic as follows:

${\mu }_i^{\prime (s)}\leftarrow u^{-1}{\alpha }_i+(1-u^{-1}){\alpha }_i,$

where α _i ≔max _j |a _ij |. This technique was developed by Rump and provided by Ozaki et al. (2013).

Algorithm 1 Mixed-precision splitting method from Mukunoki et al. (2020) for Ozaki scheme between FP64 and FP16 using floating-point arithmetic in roundto-nearest-even mode with Minamihata’s technique from Minamihata et al. (2016)
Input: $A\in {\mathbb{F}}_{64}^{m\times n}$ , $k\in \mathbb{N}$ Output: $A_s^{\prime }\in {\mathbb{F}}_{16}^{m\times n}$ , ${\mu }^{\prime (s)}\in {\mathbb{F}}_{64}^m$ , s = 1, …, k  1: β ←⌈(29 − log2n)/2⌉
2: for s = 1, …, k do
3:   ${\mu }_i^{\prime (s)}\leftarrow 2^{\lceil {\log }_2{\max }_j|a_{ij}|\rceil }$ for i = 1, …, m
4:   ${\sigma }_i\leftarrow 0.75\cdot {\mu }_i^{\prime (s)}\cdot 2^{\beta }$ for i = 1, …, m
5:   ${(A_s)}_{ij}\leftarrow (a_{ij}+{\sigma }_i)-{\sigma }_i\forall i,j$ {Extract in round-to-nearest-even mode}
6:   ${(A_s^{\prime })}_{ij}\leftarrow \mathrm{F}\mathrm{P}\mathrm{16}({({\mu }_i^{\prime (s)})}^{-1}\cdot {(A_s)}_{ij})\forall i,j$ {Convert to FP16}
7:   $a_{ij}\leftarrow a_{ij}-{(A_s)}_{ij}\forall i,j$
8: end for

After splitting matrices as in (2) and (3),

$\begin{array}{c}AB=\sum _{s+t\le k+1}\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}({\mu }^{\prime (s)})A_s^{\prime }{B}_t^{\prime }\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}({\nu }^{\prime (t)})\hfill \\ +\sum _{s+t>k+1}\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}({\mu }^{\prime (s)})A_s^{\prime }{B}_t^{\prime }\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}({\nu }^{\prime (t)})\hfill \\ +\sum _{s=1}^k\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}({\mu }^{\prime (s)})A_s^{\prime }{W}_k\hfill \\ +\sum _{s=1}^k{V}_k{B}_s^{\prime }\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}({\nu }^{\prime (s)})\hfill \\ +V_k{W}_k\hfill \end{array}$

holds and there is no rounding error in $A_s^{\prime }{B}_t^{\prime }$ for 1 ≤ s, t ≤ k on the FP16 Tensor Core. In the method implemented by Mukunoki et al. (2020), the approximation of AB was computed only using matrix multiplications with the FP16 Tensor Core, as

$AB\approx \sum _{s+t\le k+1}\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}({\mu }^{\prime (s)})A_s^{\prime }{B}_t^{\prime }\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}({\nu }^{\prime (t)})$

shown in Algorithm 2. This method was referred to as the “Fast Mode” defined by Mukunoki et al. (2020). For larger k, $|{(W_k)}_{ij}|$ and $|{(V_k)}_{ij}|$ are smaller and the approximation is more accurate.

Algorithm 2 Mixed-precision matrix multiplication method for Ozaki scheme using FP16 Tensor Core in fast mode from Mukunoki et al. (2020)
Input: $A_s^{\prime }\in {\mathbb{F}}_{16}^{m\times n}$ , ${\mu }^{\prime (s)}\in {\mathbb{F}}_{64}^m$ , $B_s^{\prime }\in {\mathbb{F}}_{16}^{n\times p}$ , ${\nu }^{\prime (s)}\in {\mathbb{F}}_{64}^n$ , s = 1, …, k Output: $C\in {\mathbb{F}}_{64}^{m\times p}$
1: C ← Om,p { $C^{\prime }\in {\mathbb{F}}_{32}^{m\times p}$ }
2: for g = 2, …, k + 1 do
3:  for s = 1, …, g − 1 do
4:    $C^{\prime }\leftarrow A_s^{\prime }{B}_{g-s}^{\prime }$ {Compute using GEMM on FP16 Tensor Core}
5:   C ← C + diag(μ′(s))FP64(C′)diag(ν′(t)) {Compute in FP64}
6:  end for
7: end for

Let ${\mathbb{I}}_N$ be a set of N-bit signed integers. It holds that −2^N−1 ≤ i ≤ 2^N−1 − 1 for all $i\in {\mathbb{I}}_N$ . Remember that no error occurs in integer arithmetic, barring overflow. For the Ozaki scheme using the INT8 Tensor Core with accumulation in INT32, the mixed-precision splitting method via bit masking with shifting is offered in Ootomo et al. (2024) and provided by Ootomo (2024). Algorithm 3 shows the splitting method. Let

$\beta ≔\min \left(7,\lfloor \frac{31-{\log }_{2}n}{2}\rfloor \right).$

(4)

Assume that β ≥ 1, i.e., n ≤ 2²⁹. That splitting method produces

$\begin{array}{l}A≕\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}({\mu }^{\prime \prime })(2^{-\beta +1}{A}_1^{\prime \prime }+\cdots +2^{-k\beta +1}{A}_k^{\prime \prime })+V_k,\hfill \end{array}$

(5)

$\begin{array}{l}B≕(2^{-\beta +1}{B}_1^{\prime \prime }+\cdots +2^{-k\beta +1}{B}_k^{\prime \prime })\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}({\nu }^{\prime \prime })+W_k\hfill \end{array}$

(6)

with double-precision shift values ${\mu }^{\prime \prime }\in {\mathbb{F}}_{64}^m$ , ${\nu }^{\prime \prime }\in {\mathbb{F}}^p$ and 8-bit integer matrices $A_i^{\prime \prime }\in {\mathbb{I}}_8^{m\times n}$ , $B_i^{\prime \prime }\in {\mathbb{I}}_8^{n\times p}$ for i = 1, …, k.

Algorithm 3 Mixed-precision splitting method for Ozaki scheme via bit masking from Ootomo et al. (2024). This algorithm is specialized for the Ozaki scheme using the INT8 Tensor Core.
Input: $A\in {\mathbb{F}}_{64}^{m\times n}$ , $k\in \mathbb{N}$ Output: $A_s^{\prime \prime }\in {\mathbb{I}}_8^{m\times n}$ , ${\mu }^{\prime \prime }\in {\mathbb{F}}_{64}^m$ , s = 1, …, k
1: β ← min(7, ⌊(31 − log2n)/2⌋)
2: ${\mu }_i^{\prime \prime }\leftarrow 2^{\lfloor {\log }_2{\max }_j|a_{ij}|\rfloor }$ for i = 1, …, m
3: for s = 1, …, k do
4:  Extract sth β bits of mantissa of a ij ∀i, j via bit masking and hold those as INT8 variable ${(A_s^{\prime \prime })}_{ij}$
5: end for

After splitting matrices as in (5) and (6),

$\begin{array}{ll}\hfill AB& =\sum _{s+t\le k+1}\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}({\mu }^{\prime \prime })2^{-\beta s+1}{2}^{-\beta t+1}{A}_s^{\prime \prime }{B}_t^{\prime \prime }\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}({\nu }^{\prime \prime })\hfill \\ \hfill & +\sum _{s+t>k+1}\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}({\mu }^{\prime \prime })2^{-\beta s+1}{2}^{-\beta t+1}{A}_s^{\prime \prime }{B}_t^{\prime \prime }\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}({\nu }^{\prime \prime })\hfill \\ \hfill & +\sum _{s=1}^k\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}({\mu }^{\prime \prime })2^{-\beta s+1}{A}_s^{\prime \prime }{W}_k\hfill \\ \hfill & +\sum _{s=1}^k{V}_k{2}^{-\beta s+1}{B}_s^{\prime \prime }\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}({\nu }^{\prime \prime })\hfill \\ \hfill & +V_k{W}_k\hfill \end{array}$

holds and there is no error in $A_s^{\prime \prime }{B}_t^{\prime \prime }$ for 1 ≤ s, t ≤ k on the INT8 Tensor Core barring overflow. In the method implemented by Ootomo et al. (2024), the approximation of AB is computed using only matrix multiplications with the INT8 Tensor Core, as

$AB\approx \sum _{s+t\le k+1}\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}({\mu }^{\prime \prime })2^{-\beta s+1}{2}^{-\beta t+1}{A}_s^{\prime \prime }{B}_t^{\prime \prime }\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}({\nu }^{\prime \prime }),$

(7)

shown in Algorithm 4. For larger k, $|{(W_k)}_{ij}|$ and $|{(V_k)}_{ij}|$ are smaller and the approximation is more accurate.

Algorithm 4 Mixed-precision matrix multiplication method for Ozaki scheme using INT8 Tensor Core from Ootomo et al. (2024)
Input: $A_s^{\prime \prime }\in {\mathbb{I}}_8^{m\times n}$ , ${\mu }^{\prime \prime (s)}\in {\mathbb{F}}_{64}^m$ , $B_s^{\prime \prime }\in {\mathbb{I}}_8^{n\times p}$ , ${\nu }^{\prime \prime (s)}\in {\mathbb{F}}_{64}^p$ , s = 1, …, k Output: $C\in {\mathbb{F}}_{64}^{m\times p}$
1: β ← min(7, ⌊(31 − log2n)/2⌋)
2: C ← Om,p { $C^{\prime \prime }\in {\mathbb{I}}_{32}^{m\times p}$ }
3: for g = 2, …, k + 1 do
4:  for s = 1, …, g − 1 do
5:    $C^{\prime \prime }\leftarrow A_s^{\prime \prime }{B}_{g-s}^{\prime \prime }$ {Compute using GEMM on INT8 Tensor Core}
6:   C ← C + diag(μ″)2−βs+12−β(g−s)+1FP64(C″)diag(ν″) {Compute in FP64}
7:  end for
8: end for

Figure 4 represents images of Algorithms 2 and 4. In Figure 4, $A_i≔\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}({\mu }^{\prime (i)})A_i^{\prime }$ and $B_i≔B_i^{\prime }\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}({\nu }^{\prime (i)})$ for Algorithms 1 and 2, and $A_i≔\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}({\mu }^{\prime \prime })2^{-i\beta +1}{A}_i^{\prime \prime }$ and $B_i≔2^{-i\beta +1}{B}_i^{\prime \prime }\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}({\nu }^{\prime \prime })$ for Algorithms 3 and 4. The exponent of A _i B _j is larger for smaller indices of i, j, so those terms have stronger effects on the accuracy of the final result.

Figure 4.

Images of matrix multiplications in Ozaki scheme for k = four.

3. Proposed methods

3.1. Splitting with rounding to nearest

We begin by explaining our motivation with a simple example. Algorithm 3 described in Section 2 used bitmask for the splitting. We assume that each number keeps three bits. Let a floating-point number be x≔377, represented by (101011111)₂ in binary format. The bitmask strategy divides x into three numbers x = x₁ + x₂ + x₃:

$\left\{\begin{array}{c}{x}_1={(\underline{101}000000)}_2,\hfill \\ x_2={(000\underline{011}000)}_2,\hfill \\ x_3={(000000\underline{111})}_2.\hfill \end{array}\right.$

(8)

In decimal format,

$x_1=256+0+64,x_2=0+16+8,x_3=4+2+1.$

Here x₁ is not the nearest number to x.

On the other hand, let x₁ and x₂ be the nearest number to x and x − x₁, respectively. Then,

$\left\{\begin{array}{c}{x}_1={(\underline{101}000000)}_2,\hfill \\ x_2={(000\underline{100}000)}_2,\hfill \\ x_3=-{(00000000\underline{1})}_2.\hfill \end{array}\right.$

(9)

In decimal format,

$x_1=256+0+64,x_2=32,x_3=-1.$

There are two drawbacks in the bitmask strategy. One drawback is that even if the leading bit is zero, we keep it, so that the number of significant figures decreases as in x₂ in (8). The other drawback is that the truncation error can be reduced. For example, if we set two slices, there is a difference in the truncation error |x₃| in (8) and (9). It is expected that the nearest strategy makes the truncation error small. Summarizing, even if we obtain a computed result with k slices using the rounding to nearest strategy, the accuracy may be comparable to that with k + 1 slices using the bitmask strategy described in Section 2 as the previous study. In this case, we must carefully monitor the following via numerical examples:

• • the increase in the splitting cost because the cost for the bitmask is low;

• • the reduction in the number of matrix multiplications; and

• • the reduction in the term for the accumulation.

Note that in the worst case, the results of the splitting with bitmask and the nearest strategy are the same. It can be explained using the following example:

$x≔{(100010001)}_2.$

For the bit mask strategy, we have

$\left\{\begin{array}{c}{x}_1={(\underline{100}000000)}_2,\hfill \\ x_2={(000\underline{010}000)}_2,\hfill \\ x_3={(000000\underline{001})}_2.\hfill \end{array}\right.$

For the rounding to the nearest strategy, we obtain

$\left\{\begin{array}{c}{x}_1={(\underline{100}000000)}_2,\hfill \\ x_2={(0000\underline{100}00)}_2,\hfill \\ x_3={(00000000\underline{1})}_2.\hfill \end{array}\right.$

The following is an algorithm for the splitting with the round to nearest strategy.

Algorithm 5 Proposed method for mixed-precision splitting for Ozaki scheme between FP64 and INT8 using floatingpoint arithmetic in round-to-nearest-even mode with Minamihata’s technique from Minamihata et al. (2016)
Input: $A\in {\mathbb{F}}_{64}^{m\times n}$ , $k\in \mathbb{N}$ Output: $A_s^{\prime }\in {\mathbb{I}}_8^{m\times n}$ , ${\mu }^{\prime (s)}\in {\mathbb{F}}_{64}^m$ , s = 1, …, k
1: β ← min(7, ⌊(31 − log2n)/2⌋)
2: for s = 1, …, k do
3:  ${\mu }_i^{\prime (s)}\leftarrow 2^{\lceil {\log }_2{\max }_j|a_{ij}|\rceil }\cdot 2^{1-\beta }$ for i = 1, …, m
4:  ${\sigma }_i=0.75\cdot 2^{53}\cdot {\mu }_i^{\prime (s)}$ for i = 1, …, m
5:  ${(A_s)}_{ij}=(a_{ij}+{\sigma }_i)-{\sigma }_i$ for i = 1, …, m {Extract in round-to-nearest-even mode}
6:  ${(A_s^{\prime })}_{ij}\leftarrow \mathrm{I}\mathrm{N}\mathrm{T}\mathrm{8}({({\mu }_i^{\prime (s)})}^{-1}\cdot {(A_s)}_{ij})\forall i,j$ {Convert to INT8}
7:  $a_{ij}\leftarrow a_{ij}-{(A_s)}_{ij}\forall i,j$
8: end for

3.2. Group-wise error-free accumulation

Next, we propose a method for accelerating the accumulation in FP64 for ozIMMU. Algorithm 4 requires the computation of a sum of k(k + 1)/2 FP64 matrices at the 6th line. The accumulation accounts for a large ratio of the total computation time of ozIMMU as shown in Figures 2 and 3. For this challenge, we propose a method for accelerating ozIMMU by reducing the number of additions in FP64.

Recall that Algorithm 4 performs

$\begin{array}{cc}\hfill & C\leftarrow \hfill \\ \hfill & \sum _{g=2}^{k+1}\sum _{s=1}^{g-1}\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}({\mu }^{\prime \prime })2^{-\beta s+1}{2}^{-\beta (g-s)+1}{A}_s^{\prime \prime }{B}_{g-s}^{\prime \prime }\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}({\nu }^{\prime \prime }).\hfill \end{array}$

(10)

The inner sum of (10) can be expressed as follows:

$\begin{array}{cc}\hfill & \sum _{s=1}^{g-1}\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}({\mu }^{\prime \prime })2^{-\beta s+1}{2}^{-\beta (g-s)+1}{A}_s^{\prime \prime }{B}_{g-s}^{\prime \prime }\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}({\nu }^{\prime \prime })\hfill \\ \hfill & =2^{-\beta g+2}\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}({\mu }^{\prime \prime })\left(\sum _{s=1}^{g-1}{A}_s^{\prime \prime }{B}_{g-s}^{\prime \prime }\right)\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}({\nu }^{\prime \prime })\hfill \end{array}$

(11)

If no overflow occurs in $A_s^{\prime \prime }{B}_{g-s}^{\prime \prime }+A_{s+1}^{\prime \prime }{B}_{g-s-1}^{\prime \prime }$ (1 ≤ s ≤ g − 2) in INT32, the summation can be performed on the accumulator of GEMM on the INT8 Tensor Core. Therefore, applying the above transformation reduces the numbers of slow conversions from INT32 to FP64, scalings, and slow summations in FP64. Let $r\in \mathbb{N}$ be

$r≔\max (1,2^{31-2\beta -\lceil {\log }_{2}n\rceil }).$

(12)

Then, the summation

$\sum _{s=1}^{\min (r,g-1)}{A}_s^{\prime \prime }{B}_{g-s}^{\prime \prime }$

can be computed without error using GEMM on the INT8 Tensor Core. The summation of r instances of $A_s^{\prime \prime }{B}_{g-s}^{\prime \prime }$ for s = 1, 2, …, r can be computed without error; however, the summation of all instances of $A_s^{\prime \prime }{B}_{g-s}^{\prime \prime }$ for s = 1, 2, …, g − 1 cannot always be computed without error. The proof validating error-free summation is provided in Section 5. Finally, we present Algorithm 6, which has a reduced number of summands in FP64 accumulation.

Algorithm 6 Proposed method for groupwise error-free accumulation for Ozaki scheme using INT8 Tensor Core
Input: $A_s^{\prime \prime }\in {\mathbb{I}}_8^{m\times n}$ , ${\mu }^{\prime \prime }\in {\mathbb{F}}_{64}^m$ , $B_s^{\prime \prime }\in {\mathbb{I}}_8^{n\times p}$ , ${\nu }^{\prime \prime (s)}\in {\mathbb{F}}_{64}^p$ , s = 1, …, k Output: $C\in {\mathbb{F}}_{64}^{m\times p}$
1: β ← min(7, ⌊(31 − log2n)/2⌋)
2: $r\leftarrow \max (1,2^{31-2\beta -\lceil {\log }_{2}n\rceil })$ {#addends for error-free accum.}
3: C ← Om,p
4: for g = 2, …, k + 1 do
5:  q ← 0
6:  C″ ← Om,p { $C^{\prime \prime }\in {\mathbb{I}}_{32}^{m\times p}$ }
7:  for s = 1, …, g − 1 do
8:   q ← q + 1
9:    $C^{\prime \prime }\leftarrow C^{\prime \prime }+A_s^{\prime \prime }{B}_{g-s}^{\prime \prime }$ {Compute using GEMM on INT8 Tensor Core}
10:   if q = r or s = g − 1 then
11:    C ← C + 2−βg+2diag(μ″)FP64(C″)diag(ν″) {Compute in FP64}
12:    q ← 0
13:    C″ ← O
14:   end if
15:  end for
16: end for

Assuming r ≥ k for simplification, all $A_s^{\prime \prime }{B}_{g-s}^{\prime \prime }$ for s = 1, 2, …, g − 1 can be accumulated without overflow for g = 3, …, k + 1. Algorithm 7 shows the Ozaki scheme with groupwise error-free accumulation for r ≥ k.

Algorithm 7 Simple version of proposed method for groupwise error-free accumulation for Ozaki scheme using INT8 Tensor Core for r ≥ k
Input: $A_s^{\prime \prime }\in {\mathbb{I}}_8^{m\times n}$ , ${\mu }^{\prime \prime }\in {\mathbb{F}}_{64}^m$ , $B_s^{\prime \prime }\in {\mathbb{I}}_8^{n\times p}$ , ${\nu }^{\prime \prime (s)}\in {\mathbb{F}}_{64}^p$ , s = 1, …, k Output: $C\in {\mathbb{F}}_{64}^{m\times p}$
1: β ← min(7, ⌊(31 − log2n)/2⌋) = 7
2: C ← Om,p
3: for g = 2, …, k + 1 do
4:  C″ ← Om,p { $C^{\prime \prime }\in {\mathbb{I}}_{32}^{m\times p}$ }
5:  for s = 1, …, g − 1 do
6:    $C^{\prime \prime }\leftarrow C^{\prime \prime }+A_s^{\prime \prime }{B}_{g-s}^{\prime \prime }$ {Compute using GEMM on INT8 Tensor Core}
7:  end for
8:  C ← C + diag(μ″)2−βg+2FP64(C″)diag(ν″) {Compute in FP64}
9: end for

3.3. Combination of proposals

Next, we accelerate ozIMMU by reducing the number of summands in FP64 accumulation and improving the splitting method. For this purpose, we combine the methods proposed in Sections 3.1 and 3.2. In order to use group-wise error-free accumulation as in Algorithm 6, A and B need to be expressed as in (5) and (6). Thus, we determine

$\begin{array}{l}\hfill {\mu }_i^{\prime \prime }\leftarrow 2^{\lceil {\log }_2{\max }_j|a_{ij}|\rceil }\cdot 2^{1-\beta },\\ \hfill {\sigma }_i\leftarrow 0.75\cdot 2^{53-\beta (s-1)}\cdot {\mu }_i^{\prime \prime }\end{array}$

as the 3rd and 4th lines in Algorithm 5. Then, A _s are extracted using rounding to nearest floating-point arithmetic as

${(A_s)}_{ij}\leftarrow (a_{ij}+2^{\beta (s-1)}{\sigma }_i)-2^{\beta (s-1)}{\sigma }_i$

for the constant 2^β(s−1)σ _i with a common ratio of 2 ^β . Finally, we obtain Algorithm 8. The shift values for A _s are ${\mu }_i^{\prime \prime (s)}\cdot 2^{\beta (s-1)}$ ; thus, ${\mu }_i^{\prime \prime (s)}$ are determined before the “for” statement and the shift values for A _s for s ≥ 1 are calculated by shifting ${\mu }_i^{\prime \prime (s)}$ by a constant ratio 2 ^β . Hence, the algorithm finds the maximum absolute values max _j |a _ij | for i = 1, …, m only once before the “for” statement. On the other hand, the shift values ${\mu }_i^{\prime (s)}$ are determined inside the “for” statement in Algorithm 5. Therefore, the algorithm finds the maximum absolute values max _j |a _ij | for i = 1, …, m at most k times.

Algorithm 8 Another proposed method for mixed-precision splitting for Ozaki scheme between FP64 and INT8 using floating-point arithmetic in round-to-nearest-even mode with Minamihata’s technique from Minamihata et al. (2016). The results can be expressed as in (5) for error-free group-wise accumulation.
Input: $A\in {\mathbb{F}}_{64}^{m\times n}$ , $k\in \mathbb{N}$ Output: $A_s^{\prime \prime }\in {\mathbb{I}}_8^{m\times n}$ , ${\mu }^{\prime \prime }\in {\mathbb{F}}_{64}^m$ , s = 1, …, k
1: β ← min(7, ⌊(31 − log2n)/2⌋)
2: ${\mu }_i^{\prime \prime }\leftarrow 2^{\lceil {\log }_2{\max }_j|a_{ij}|\rceil }\cdot 2^{1-\beta }$ for i = 1, …, m
3: for s = 1, …, k do
4:   ${\sigma }_i\leftarrow 0.75\cdot 2^{53-\beta (s-1)}\cdot {\mu }_i^{\prime \prime }$ for j = 1, …, m
5:   ${(A_s)}_{ij}\leftarrow (a_{ij}+{\sigma }_i)-{\sigma }_i\forall i,j$ {Extract in round-to-nearest-even mode}
6:   ${(A_s^{\prime \prime })}_{ij}\leftarrow \mathrm{I}\mathrm{N}\mathrm{T}\mathrm{8}({{\mu }_i^{\prime \prime }}^{-1}\cdot 2^{\beta (s-1)}\cdot {(A_s)}_{ij})\forall i,j$ {Convert to INT8}
7:   $a_{ij}\leftarrow a_{ij}-{(A_s)}_{ij}\forall i,j$
8: end for

4. Numerical examples

We have implemented the proposed methods in the GitHub repository by Uchino (2024), replacing specific code in Ootomo’s ozIMMU library with our code. All numerical experiments were conducted on NVIDIA GH200 Grace Hopper Superchip and NVIDIA GeForce RTX 4090 with the GNU C++ Compiler 11.4.1 and NVIDIA CUDA Toolkit 12.5.82. The tested methods will be denoted as follows:

• ozIMMU-k: Ootomo’s implementation with k slices

• ozIMMU_RN-k: Proposed method in Section 3.1 with k slices

• ozIMMU_EF-k: Proposed method in Section 3.2 with k slices

• ozIMMU_H-k: Proposed method in Section 3.3 with k slices.

4.1. Accuracy

Figure 5 shows the accuracies of ozIMMU-k, ozIMMU_RN-k, ozIMMU_EF-k, and ozIMMU_H-k for k = 7, …, 11. We set $A\in {\mathbb{F}}_{64}^{n\times n}$ as a _ij ≔(U _ij − 0.5) ⋅ exp(ϕ ⋅ N _ij ), where U _ij ∈ (0, 1) are uniformly distributed random numbers and N _ij are drawn from the standard normal distribution, for 1 ≤ i, j ≤ n. The constant ϕ specifies the tendency of the difficulty in terms of matrix multiplication accuracy. A matrix $B\in {\mathbb{F}}_{64}^{n\times n}$ is set similarly to A.

Figure 5.

Comparison of accuracy between ozIMMU and proposed methods. Matrix A has entries a _ij ≔(U _ij − 0.5) ⋅ exp(ϕ ⋅ N _ij ), where U _ij ∈ (0, 1) are uniformly distributed random numbers and N _ij are drawn from the standard normal distribution, for 1 ≤ i, j ≤ m and m = n = p. Matrix B has a similar composition.

The accuracy deteriorated for n > 2¹⁷ because the maximum number of significant bits β = min(7, ⌊(31 − log₂n)/2⌋) of the elements of the split matrices is less than 7. The accuracies of ozIMMU and ozIMMU_EF are comparable, as are those of, ozIMMU_RN and ozIMMU_H. ozIMMU_RN and ozIMMU_H, which use a splitting method via rounding to nearest floating-point arithmetic, are more accurate than ozIMMU and ozIMMU_EF, which use a splitting method via bit masking. Occasionally, to obtain results comparable with that of FP64, ozIMMU_RN-k, ozIMMU_H-k, ozIMMU-(k + 1), and ozIMMU_EF-(k + 1) are required; i.e., ozIMMU_RN and ozIMMU_H require fewer splits than ozIMMU and ozIMMU_EF. In particular, for ϕ = 2, ozIMMU_RN-9 and ozIMMU_H-9 produce comparable results to FP64; however, the accuracies of the results of ozIMMU-9 and ozIMMU_EF-9 are worse than that of FP64. In such cases, k = 10 is required for ozIMMU and ozIMMU_EF to attain more accurate results than FP64.

4.2. Time breakdown

Figures 6–8 and Figures 9–11 show time breakdowns of the proposed methods on RTX 4090 and GH200, respectively. Note that “split A”, “split B”, “Gemm in INT8-TC”, “accumulation in FP64”, and “copy” correspond to the steps (i), (ii), (iii), (iv), and (v) in the Ozaki scheme for emulating the GEMM routine as shown in Section 1.

Figure 6.

Time breakdown of ozIMMU_RN on NVIDIA GeForce RTX 4090.

Figure 7.

Time breakdown of ozIMMU_EF on NVIDIA GeForce RTX 4090.

Figure 8.

Time breakdown of ozIMMU_H on NVIDIA GeForce RTX 4090.

Figure 9.

Time breakdown of ozIMMU_RN on NVIDIA GH200 Grace Hopper Superchip.

Figure 10.

Time breakdown of ozIMMU_EF on NVIDIA GH200 Grace Hopper Superchip.

Figure 11.

Time breakdown of ozIMMU_H on NVIDIA GH200 Grace Hopper Superchip.

The execution time of FP64 accumulation in ozIMMU_RN is not negligible, because the number of additions in FP64 is not reduced. The execution times of FP64 accumulation in ozIMMU_EF and ozIMMU_H are less than those in ozIMMU and ozIMMU_RN. From Figures 6–8, the computation time of splitting via rounding to nearest floating-point arithmetic is not so fast because FP64 is much slower than the lower-precision arithmetic on RTX 4090. From Figures 10 and 11, the ratios of the computation time of splitting in ozIMMU_EF and ozIMMU_H are comparable on GH200.

4.3. Performance

Figures 12 and 13 show throughput in TFLOPS and ratio to ozIMMU. A smaller number of slices corresponds to better performance because the total computation cost is less.

Figure 12.

Throughput comparison (left) and speedups to ozIMMU (right) on NVIDIA GeForce RTX 4090.

Figure 13.

Throughput comparison (left) and speedups to ozIMMU (right) on NVIDIA GH200 Grace Hopper Superchip.

On RTX 4090, all methods are faster than or comparable to matrix multiplication in FP64 for n ≥ 1024 because FP64 is much slower than the lower-precision arithmetic. ozIMMU_EF and ozIMMU_H are faster than ozIMMU almost everywhere. In particular, ozIMMU_EF-12 and ozIMMU_H-12 are respectively 1.6 and 1.4 times faster than ozIMMU for n = 4096. ozIMMU_RN is slower than ozIMMU. For k = 3, ozIMMU, ozIMMU_RN, ozIMMU_EF, and ozIMMU_H are 38.8, 33.8, 44.4, and 40.0 times faster than matrix multiplication in FP64 for n = 16384, respectively. For k = 12, ozIMMU, ozIMMU_RN, ozIMMU_EF, and ozIMMU_H are 1.4, 1.0, 1.9, and 1.6 times faster than matrix multiplication in FP64 for n = 1024, respectively. In addition, ozIMMU, ozIMMU_RN, ozIMMU_EF, and ozIMMU_H are 9.4, 8.5, 12.0, and 10.9 times faster than FP64 for k = 7 and n = 16384, and 7.4, 6.8, 9.5, and 8.6 times faster than FP64 for k = 8 and n = 16384.

On GH200, the methods with small k are faster than matrix multiplication in FP64 for n ≥ 4096; however, the methods are much slower than FP64 otherwise. ozIMMU_EF and ozIMMU_H are faster than ozIMMU. In particular, ozIMMU_EF-11 and ozIMMU_H-11 are respectively 1.5 and 1.6 times faster than ozIMMU for n = 2048. ozIMMU_RN is not as slow on GH200 as it is on RTX 4090. For k = 3, ozIMMU, ozIMMU_RN, ozIMMU_EF, and ozIMMU_H are 2.7, 2.2, 3.0, and 3.0 times faster than matrix multiplication in FP64 for n = 16384, respectively. ozIMMU-5, ozIMMU_RN-5, ozIMMU_EF-6, and ozIMMU_H-6 have comparable computation times to matrix multiplication in FP64 for n = 16384; however, ozIMMU-k, ozIMMU_RN-k, ozIMMU_EF-(k + 1), and ozIMMU_H-(k + 1) are slower than FP64 for n < 16384 or k > 5.

4.4. Performance vs. Accuracy

Figures 14 and 15 show throughput in a TFLOPS vs. Accuracy plot, illustrating relationships between performance and accuracy for the Ozaki scheme and for matrix multiplication in FP64 for n = 4096 and ϕ = 0. The numbers inside the symbols are the numbers of splices for the Ozaki scheme.

Figure 14.

Relationship between performance and accuracy for n = 4096 and ϕ = 0 on NVIDIA GeForce RTX 4090.

Figure 15.

Relationship between performance and accuracy for n = 4096 and ϕ = 0 on NVIDIA GH200 Grace Hopper Superchip.

On RTX 4090, all methods are better than matrix multiplication in FP64 with respect to both performance and accuracy for k ≥ 8. ozIMMU_H and ozIMMU_EF produce results with comparable performances that are more accurate than those of ozIMMU. In particular, ozIMMU_H-8 and ozIMMU_EF-9 produce results with comparable performances that are more accurate than those of ozIMMU-7. On GH200 as well, ozIMMU_H and ozIMMU_EF produce results with comparable performances that are more accurate than those of ozIMMU. In particular, ozIMMU_H-7 and ozIMMU_EF-7 produce results with comparable performances that are more accurate than those of ozIMMU-6. These results also indicate that ozIMMU_H and ozIMMU_EF can be computed with less memory consumption than ozIMMU.

5. Rounding error analysis

In this section, we give a rounding analysis for a fixed number of slices. We have described Algorithms 3 and 4 (ozIMMU), and proposed Algorithms 5 and 6, Algorithms 3 and 6, and Algorithms 8 and 6. Even if we used the rounding to nearest strategy as in Algorithms 5 and 8, their error bounds would be the same as that of the bitmask strategy used in Algorithm 5. Therefore, we focus on Algorithms 3 and 4 (ozIMMU) and Algorithms 3 and 6. Below, absolute value notation applied to a matrix means the matrix from applying absolute value element-wise. We assume that neither overflow nor underflow occurs. For $A\in {\mathbb{F}}^{m\times n}$ and $B\in {\mathbb{F}}^{n\times p}$ , the following deterministic error bound is given by Jeannerod and Rump (2013):

$|AB-\mathtt{f}\mathtt{l}(AB)|\le nu|A||B|.$

The following alternative probabilistic error bound comes from Higham and Mary (2019) (which has the details on the assumptions and probabilities):

$|AB-\mathtt{f}\mathtt{l}(AB)|\lesssim \sqrt{n}u|A||B|.$

In this section, our aim is to derive an error bound on the computed result. Let $T\in {\mathbb{F}}^{m\times p}$ be a computed result using k slices, such as

$A\approx A_1+A_2+\cdots +A_k,B\approx B_1+B_2+\cdots +B_k.$

Note that for i = 1, 2, …, k,

$A_i≔\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}({\mu }^{\prime (i)})A_i^{\prime },B_i≔B_i^{\prime }\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}({\nu }^{\prime (i)})$

for Algorithms 1, 2, and 5, and

$A_i≔\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}({\mu }^{\prime \prime })2^{-i\beta +1}{A}_i^{\prime \prime },B_i≔2^{-i\beta +1}{B}_i^{\prime \prime }\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}({\nu }^{\prime \prime })$

for Algorithms 3, 4, 6, 7, and 8. Let T _k be an approximation of AB with k slices. Our aim is to obtain the upper bound on |AB − T| as follows:

$\begin{array}{cc}\hfill & |AB-T_k|\hfill \\ \hfill & \le \left|AB-\sum _{i=1}^k\sum _{j=1}^{k-i+1}{A}_i{B}_j\right|+\left|\sum _{i=1}^k\sum _{j=1}^{k-i+1}{A}_i{B}_j-T_k\right|.\hfill \end{array}$

(13)

Note that we can obtain the exact product A _i B _j in the above by using GEMM in the INT8 Tensor Core and scaling by powers of two. The first term of the bound in (13) indicates the truncation error, and the second term shows an error arising in the accumulation process. As matrix scaling does not affect rounding errors, a discussion on it is omitted.

In the following subsections, we use an upper bound provided by Jeannerod and Rump (2013) for a sum:

$\left|\sum _{i=1}^{n}p-\mathtt{f}\mathtt{l}\left(\sum _{i=1}^{n}p\right)\right|\le (n-1)u\sum _{i=1}^n|p_i|,p\in {\mathbb{F}}^n.$

(14)

Let $\mathtt{u}\mathtt{f}\mathtt{p}(c)$ for $c\in \mathbb{R}$ be defined as

$\mathtt{u}\mathtt{f}\mathtt{p}(c)≔\left\{\begin{array}{cc}0\hfill & \text{if}c=0,\hfill \\ 2^{\lfloor {\log }_2|c|\rfloor }\hfill & \text{otherwise}.\hfill \end{array}\right.$

The following condition is used to check whether a real number is a floating-point number. For $c\in \mathbb{R}$ , if |c| is smaller than the maximum floating-point number and c is an integral multiple of the minimum positive floating-point number, then it holds that

$c\in 2u\cdot \mathtt{u}\mathtt{f}\mathtt{p}(c)\mathbb{Z}\Rightarrow c\in \mathbb{F}.$

(15)

5.1. Error bound for Algorithms 3 and 4 (ozIMMU)

Let matrices $V_k\in {\mathbb{F}}^{m\times n}$ and $W_k\in {\mathbb{F}}^{n\times p}$ be defined as in (1) where these show the truncation errors of A and B for the case of k slices, respectively. Notation e indicates the vector $e={(\mathrm{1,1},\dots ,1)}^T\in {\mathbb{F}}^n$ . Define two vectors $g\in {\mathbb{F}}^m$ and $h\in {\mathbb{F}}^p$ related to row-wise maximums in A and column-wise maximums in B in the sense of the unit in the first place as

$g_i≔\mathtt{u}\mathtt{f}\mathtt{p}\left(\underset{j}{\max }|a_{ij}|\right),f_j≔\mathtt{u}\mathtt{f}\mathtt{p}\left(\underset{i}{\max }|b_{ij}|\right).$

Then, from Figure 16 and (1), we have

$|V_i|\le 2^{-\beta i+1}g{e}^T,|W_i|\le 2^{-\beta i+1}e{f}^T,$

(16)

and matrices |A _i | and |B| are bounded as

$|A_i|\le 2^{-\beta (i-1)+1}g{e}^T,|B|\le 2e{f}^T.$

(17)

Because

$AB=(A_1+\cdots +A_k+V_k)(B_1+\cdots +B_k+W_k),$

we have

$AB-\sum _{i=1}^k\sum _{j=1}^{k-i+1}{A}_i{B}_j=\sum _{i=1}^k{A}_i{W}_{k-i+1}+V_{k}B.$

From the last equation, (16), and (17), the truncation error is bounded as

$\begin{array}{ll}\hfill & \left|AB-\sum _{i=1}^k\sum _{j=1}^{k-i+1}{A}_i{B}_j\right|\hfill \\ \hfill & \le \sum _{i=1}^k\left|A_i\right|\left|W_{k-i+1}\right|+\left|V_k\right|\left|B\right|\hfill \\ \hfill & \le \sum _{i=1}^k{2}^{-\beta (i-1)+1}g{e}^T\cdot 2^{-\beta (k-i+1)+1}e{f}^T\hfill \\ \hfill & +2\cdot 2^{-\beta k+1}g{e}^T\cdot e{f}^T\hfill \\ \hfill & =4n\sum _{i=1}^k{2}^{-\beta k}g{f}^T+4n\cdot 2^{-\beta k}g{f}^T\hfill \\ \hfill & =4(k+1)n{2}^{-\beta k}g{f}^T.\hfill \end{array}$

(18)

Figure 16.

An image of (24).

If |g _i − |a _ij ‖ for 1 ≤ j ≤ n and |f _j − |b _ij ‖ for 1 ≤ i ≤ n are very small, we can assume that

$4ng{f}^T\approx |A||B|.$

(19)

In this case, we have

$\left|AB-\sum _{i=1}^k\sum _{j=1}^{k-i+1}{A}_i{B}_j\right|\lesssim (k+1)2^{-\beta k}|A\Vert B|.$

(20)

We next focus on an error in the accumulation, the second term in (13). If A _i and B _i are obtained by Algorithm 3,

$|A|=\sum _{i=1}^k|A_i|,|B|=\sum _{i=1}^k|B_i|$

(21)

are satisfied. Because we compute the sum of k(k + 1)/2 floating-point matrices, the following bound is immediately obtained from (14) and (21):

$\begin{array}{ll}\hfill & \left|\sum _{i=1}^k\sum _{j=1}^{k-i+1}{A}_i{B}_j-T_k\right|\hfill \\ \hfill & \le \left(\frac{1}{2}k(k+1)-1\right)u\sum _{i=1}^k\sum _{j=1}^{k-i+1}|A_i{B}_j|\hfill \\ \hfill & \le \left(\frac{1}{2}k(k+1)-1\right)u|A||B|.\hfill \end{array}$

(22)

We will show that k′ (≤k) exists such that

$\sum _{s=1}^{k^{\prime }}\sum _{t=1}^{k^{\prime }-s+1}{A}_s{B}_t\in {\mathbb{F}}^{m\times n}.$

(23)

Then, the bound (22) is improved. From Figure 16, we have

${\left(A_s\right)}_{ij}\in 2^{-s\beta +1}{g}_i\mathbb{Z},{\left(B_t\right)}_{ij}\in 2^{-t\beta +1}{f}_j\mathbb{Z},$

(24)

and these derive

${\left(A_s\right)}_{i\ell }{\left(B_t\right)}_{\ell j}\in 2^{-\beta (s+t)+2}{g}_i{f}_j\mathbb{Z}.$

(25)

Figure 16 shows an image of (24). It follows from this that

${\left(\sum _{s=1}^{k^{\prime }}\sum _{t=1}^{k^{\prime }-s+1}{A}_s{B}_t\right)}_{ij}\in 2^{-\beta (k^{\prime }+1)+2}{g}_i{f}_j\mathbb{Z}.$

(26)

From Figure 16, the elements of the matrices are bounded as

$|A_s{|}_{ij}\le 2^{-\beta (s-1)+1}{g}_i,|B_t{|}_{ij}\le 2^{-\beta (t-1)+1}{f}_j,$

so that we have

$|A_s{B}_t{|}_{ij}\le n{2}^{-\beta (s+t-2)+2}{g}_i{f}_j.$

(27)

Therefore, from (27), we have the following bound on a dot product:

$\begin{array}{ll}\hfill & {\left|\sum _{s=1}^{k^{\prime }}\sum _{t=1}^{k^{\prime }-s+1}{A}_s{B}_t\right|}_{ij}\hfill \\ \hfill & ={\left|\sum _{s=2}^{k^{\prime }+1}\sum _{i_2+j_2=s}{A}_{i_2}{B}_{j_2}\right|}_{ij}\hfill \\ \hfill & \le {\left(\sum _{s=2}^{k^{\prime }+1}\sum _{i_2+j_2=s}|A_{i_2}||B_{j_2}|\right)}_{ij}\hfill \\ \hfill & \le \sum _{s=2}^{k^{\prime }+1}\sum _{i_2+j_2=s}n{2}^{-(i_2+j_2-2)\beta +2}{g}_i{f}_j\hfill \\ \hfill & =\sum _{s=2}^{k^{\prime }+1}\sum _{i_2+j_2=s}n{2}^{-(s-2)\beta +2}{g}_i{f}_j\hfill \\ \hfill & =\sum _{s=2}^{k^{\prime }+1}(s-1)n{2}^{-(s-2)\beta +2}{g}_i{f}_j.\hfill \end{array}$

If β ≥ 3,

$\sum _{s=3}^{k^{\prime }+1}(s-1)n{2}^{-(s-2)\beta +2}{g}_i{f}_j\le n{2}^2{g}_i{f}_j$

Finally, for β ≥ 3, we have

${\left|\sum _{s=1}^{k^{\prime }}\sum _{t=1}^{k^{\prime }-s+1}{A}_s{B}_t\right|}_{ij}\le 8n\cdot g_i{f}_j$

(28)

From (26), (28), and (15), if

$2u\cdot \mathtt{u}\mathtt{f}\mathtt{p}(8n)\le 2^{-\beta (k^{\prime }+1)+2}$

(29)

is satisfied, then (23) holds.

Let $k_{\max }^{\prime }$ be the maximum k′ satisfying (29). No rounding error occurs for $i=\mathrm{2,3},\dots ,k_{\max }^{\prime }$ . Therefore, we have

$\begin{array}{cc}\hfill & \left|\sum _{k=1}^k\sum _{j=1}^{k-i+1}{A}_i{B}_j-T_k\right|\hfill \\ \hfill & \le \left(\frac{1}{2}k(k+1)-\frac{1}{2}{k}_{\max }^{\prime }(k_{\max }^{\prime }+1)-1\right)u|A||B|.\hfill \end{array}$

(30)

Summarizing, from (18) and (30), we have

$\begin{array}{cc}\hfill & |AB-T_k|\hfill \\ \hfill & \le 4(k+1)n{2}^{-\beta k}g{f}^T\hfill \\ \hfill & +\left(\frac{1}{2}k(k+1)-\frac{1}{2}{k}_{\max }^{\prime }(k_{\max }^{\prime }+1)-1\right)u|A||B|.\hfill \end{array}$

If (19) is satisfied, then we have

$\begin{array}{cc}\hfill & |AB-T_k|\hfill \\ \hfill & \lesssim (k+1)2^{-\beta k}|A\Vert B|\hfill \\ \hfill & +\left(\frac{1}{2}k(k+1)-\frac{1}{2}{k}_{\max }^{\prime }(k_{\max }^{\prime }+1)-1\right)u|A||B|.\hfill \end{array}$

5.2. Error bound for Algorithms 3 and 6

We focus on

$\sum _{i_2+j_2=s}{A}_{i_2}^{\prime \prime }{B}_{j_2}^{\prime \prime },A_{i_2}^{\prime \prime }\in {\mathbb{I}}_{\beta +1}^{m\times n},B_{j_2}^{\prime \prime }\in {\mathbb{I}}_{\beta +1}^{n\times p}.$

(31)

We consider r such that the following holds

$\sum _{\begin{array}{c}{i}_2+j_2=s\\ s\le \min (r,k+1)\end{array}}{A}_{i_2}^{\prime \prime }{B}_{j_2}^{\prime \prime }\in {\mathbb{I}}_{32}^{m\times p}.$

(32)

Let E_m,n be the matrix ${\mathbb{I}}_2^{m\times n}$ whose elements are all ones. From the definitions in (31),

$|A_{i_2}^{\prime \prime }|\le (2^{\beta }-1)E_{m,n},|B_{i_2}^{\prime \prime }|\le (2^{\beta }-1)E_{n,p},$

and we have

$\begin{array}{ll}\hfill \left|\sum _{i_2+j_2=s}{A}_{i_2}^{\prime \prime }{B}_{j_2}^{\prime \prime }\right|& \le \sum _{i_2+j_2=s}|A_{i_2}^{\prime \prime }||B_{j_2}^{\prime \prime }|\hfill \\ \hfill & \le (s-1)n{(2^{\beta }-1)}^2{E}_{m,p}.\hfill \end{array}$

Since the largest number in ${\mathbb{I}}_{32}$ is 2³¹ − 1, if s ≤ r for $r≔\max (1,2^{31-2\beta -\lceil {\log }_{2}n\rceil })$ in (12), we have

$\begin{array}{ll}\hfill (s-1)n{(2^{\beta }-1)}^2& \le (r-1)n{(2^{\beta }-1)}^2\hfill \\ \hfill & \le (r-1)\cdot 2^{\lceil {\log }_{2}n\rceil }\cdot 2^{2\beta }\hfill \\ \hfill & =r\cdot 2^{\lceil {\log }_{2}n\rceil }\cdot 2^{2\beta }-2^{\lceil {\log }_{2}n\rceil }\cdot 2^{2\beta }\hfill \\ \hfill & \le 2^{31}-1,\hfill \end{array}$

which implies that (32) holds. Note that

$(s-1)n{(2^{\beta }-1)}^2<2^{31}-1$

for β ≥ 1. This is why constant r was set as in (12).

If we define the number of the terms, w, for the accumulation as

$w≔\lceil \frac{k}{r}\rceil \left(k-\frac{r}{2}⌊\frac{k-1}{r}⌋\right),$

we have

$|AB-T|\le 4(k+1)n{2}^{-\beta k}g{f}^T+(w-1)u|A||B|.$

If (19) is satisfied, then we have

$|AB-T|\lesssim (k+1)2^{-\beta k}|A\Vert B|+(w-1)u|A||B|.$

6. Conclusion

We proposed three implementation methods for accelerating the Ozaki scheme using the INT8 Tensor Core. In the original implementation, ozIMMU, provided by Ootomo et al., the ratio of the accumulation in FP64 is not negligible, as it accounts for approximately 40–50 % of the total computation time. The proposed methods ozIMMU_EF and ozIMMU_H reduce the computation time ratio of the accumulation to approximately 10–20 % and achieve a 1.2- to 1.6-fold speedup. With future architectures expected to achieve high speed in lower-precision matrix multiplication, the computation time ratio of the accumulation in FP64 is expected to increase. Thus, these proposed methods contribute to leveraging the performance of future architectures more effectively than ozIMMU. The proposed methods ozIMMU_RN and ozIMMU_H offer alternative splitting methods using floating-point arithmetic in round-to-nearest-even mode and produce more accurate results than ozIMMU for the same number of slices.