Multi-limb integer matrix multiplication using AMX and AVX512-IFMA instructions Joris van der HoevenA, Marc MezzarobbaB, Laboratoire d'informatique de l'École polytechnique Preliminary version of October 9, 2026

. This work has been supported by an ERC-2023-ADG grant for the ODELIX project (number 101142171).

Funded by the European Union. Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or the European Research Council Executive Agency. Neither the European Union nor the granting authority can be held responsible for them.

We study the hardware-accelerated multiplication of matrices with integer entries of 8 to about 800 bits. Our focus is on two relatively new hardware features in Intel CPUs: the IFMA “integer FMA” instruction and the AMX matrix extensions. We survey a number of approaches to the multiplication of matrices of multiple-precision integers using these instruction sets and analyze to what extent our implementations on top of the JIL library approach theoretical peak performance.

1Introduction

Motivation

Arguably the most significant recent trend in computer hardware is the eruption of tensor cores in both GPUs and CPUs. While this trend is driven by applications in artificial intelligence, is there a potential for tensor cores to be used for other purposes? After all, GPUs have found applications well beyond their initial ones in the gaming industry.

The main motivation of this paper is to study this question for applications in scientific computation and more particularly in computer algebra and computer arithmetic. A basic question is that of matrix multiplication with -bit integer entries. Tensor cores are typically designed to accelerate such computations when the matrix dimensions are large and . Can this be used to efficiently deal with other bit precisions, ranging from small ones like to larger ones like ? We focus here on the case of square matrices. This is indeed most favorable for tensor cores, but also most significant from a complexity perspective, since the complexity of many problems in linear algebra and beyond can be expressed in terms of the complexity of square matrix products [5]. This restriction will also simplify comparisons, because varying both and will already give rise to numerous different cases, as we shall see.

To be fair, we should compare implementations that exploit tensor cores with alternatives based on more conventional vector arithmetic. For this reason, we focus on recent Intel Xeon processors which support both AMX matrix products and AVX512 vector arithmetic with IFMA. Using AMX, each core can do a matrix FMA (MFMA) in about cycles, where and (of dimensions and ) have -bit entries and (of dimension ) has 32-bit entries. Roughly speaking, Intel's AVX512 integer FMA (IFMA) can compute in one or half a cycle, where and are -wide vectors of -bit integers and is an -wide vector of -bit integers. Note that we preferred the terminology MFMA over GeMM (general matrix multiplication, ) because AMX only supports positive accumulation.

Recent AMD processors also support AVX-512 and the IFMA instruction, and part of our work also applies to these processors. However, AMD processors do not yet support the AMX instruction set. Current ARM processors only support 128-bit-wide Neon vector arithmetic (although approximately twice as efficient as on X86 architectures) and there is no equivalent of the IFMA instruction. The matrix coprocessors integrated in Apple's M1–M5 CPUs are based on fast pipelined outer products [10], not dot products like AMX-style matrix multiplication, so that adapting our work to such processors would require significant changes.

In most of the article, we focus on the performance on a single CPU core. This simplification allows us to compare implementation strategies that exploit different hardware features without being distracted by largely orthogonal considerations on how to scale up these base implementations on multi-core architectures. On request of one of the referees, though, we briefly present very preliminary work on a multi-threaded implementation at the end of the paper.

The timings presented in the main body of the article are relative to theoretical peak performance. For instance, for a naive implementation of matrix multiplication using IFMAs, we give the measured cost in cycles divided by the theoretical cost of IFMAs based on published instruction throughputs. These ratios are fairly independent of the particular machine and indicate the quality of our implementation. Quotients close to one mean that we could not have done much better for the algorithm at hand. Higher ratios indicate potential problems such as suboptimal memory access patterns, asymptotically subdominant costs that become high for certain sizes, or insurmountable hardware obstacles. In an ideal world in which all ratios are one, the theoretical complexity analysis of our methods would correspond perfectly to their practical performance. For completeness, we also provide an appendix with absolute timings, both for our new implementation and some existing software.

All implementations in this paper are done inside or on top of the C++ library JIL [2, 41] for computations with straight-line programs (SLPs). This work is the first step of a project to equip the JIL library with fast linear algebra over the integers and fixed point numbers. We would like to cover all ranges of matrix dimensions and bit-sizes of the integers, while taking advantage of recent hardware features. Conversely, matrix multiplication is an excellent problem for extending JIL beyond the mere compilation of SLPs: another important motivation is to use it as a case study for how to integrate SLPs inside larger pieces of code such as multiple loops. Integrating AMX instructions into the SLP framework is also a challenge.

Previous work

Multiple precision integer multiplication has a long history and there exist many algorithms for various bit precisions, such as the traditional schoolbook method, Karatsuba's algorithm [44], Toom's method [65], and modular methods [25, Chapter 5]. For high precisions, there are also methods based on the fast Fourier transform [11, 59, 35], but we will not use them in this paper. GMP [29] and MPFR [32] are the reference libraries for multiple precision integer and floating point arithmetic. We also refer to [39, 30, 19, 8, 40] for some recent work on this topic in the HPC context.

In the context of integer matrix multiplication, fast multiple precision methods tend to become useful already at very low precisions. This is due to the fact that all known methods can be regarded as modular or evaluation-interpolation methods. Roughly speaking, multiplying two matrices with -bit entries is reduced to multiplications of matrices with -bit entries, where the quantity depends on the method (this will be detailed in Section 3 and Table 1 below). The cost of this reduction scales with , whereas the cost of the low precision matrix multiplications scales with . Consequently, the cost of the reduction becomes negligible for large and, for a given , it is best to pick the method for which is minimal. Good practical implementations should reflect this fact and this is indeed what our work aims at.

Theoretically speaking, matrices can be multiplied in time with [12] for large . Among the algorithms with , only Strassen's algorithm [63] with has allowed for some modest practical gains so far. However, due to the high arithmetic intensity of matrix multiplication and now to the availability of dedicated hardware accelerators, the practical complexity of matrix multiplication behaves as if were well below , making matrix multiplication a crucial building block for other fast algorithms, both in theory and in practice.

Classical reference implementations for matrix multiplication and other routines from linear algebra include Blas [6], Lapack [48], Atlas [3], Blis [7] and OpenBlas [64]. We refer to [49, 16, 15, 31, 69, 28, 67, 52] for some of the underlying scientific background. For matrices with multiple precision integer entries, the Chinese remainder theorem (CRT) is typically used to reduce matrix multiplication to modular matrix multiplication, which is then performed using a standard BLAS [17, 14]. When working with very high precisions, FFT-based techniques are also used [36, 34]. The Flint library [33] performs well on a large spectrum of matrix dimensions and integer sizes. We refer to [42, 22] for a discussion of SIMD-accelerated modular arithmetic and to [4, 66] for CRT-based matrix multiplication on modern hardware. Doubling the precision of a BLAS has been considered in [56].

Several works [13, 20] and research implementations [51, 50] specifically address the challenges raised by Intel's AMX units. OpenBlas [64] also contains a bf16 kernel (sbgemm_kernel_16x16_spr_tmpl.c) using AMX instructions, and GGML [27] an int8 kernel interleaved with encoding and decoding of quantized tensors (ggml-cpu/amx/mmq.cpp). Separate techniques have been developed for other tensor accelerator architectures, including the matrix coprocessors of Apple M1–M3 CPUs [21]. We refer to [53] for an overview of the types of matrix multiplication accelerators available in current CPUs and to [43, Chapter 20] for details about Intel's AMX units. The use of tensor cores to multiply matrices of floating-point numbers of precisions beyond the word length of the accelerator is studied in [55, 54, 66, 1]. Georganas et al. [26] propose a batch-reduce kernel that may also be useful for this purpose.

Finally, our work on JIL draws inspiration from the compilation of codelets for critical tasks pioneered in Fftw3 [23, 24] and Spiral [57]. Matrix multiplication has also been studied from a similar compilation-oriented perspective [68, 46, 47].

Contributions

We present an experimental high performance implementation of integer matrix multiplication for a wide range of matrix dimensions and bit sizes of integers. We have spent a comparable amount of effort on exploiting two types of hardware acceleration, namely the IFMA and AMX extensions in recent Intel processors. Our implementation is available at https://sourcesup.renater.fr/projects/jil/ on the amx branch and distributed under the GNU GPLv3. The version described in the present paper corresponds to revision a04ea451, based on JIL revision 999ee20d.

In Table 1, we present an idealized complexity model for large matrix dimensions as a function of the bit size of the coefficients and the integer multiplication method being used. Although the underlying methods are not new, we expect this synthesis to be useful as a guideline for implementors. Our IFMA-based kernels, and to a lesser extent the AMX-based ones, often come close to the theoretical peak performance for which the idealized model is most relevant. In absolute terms, we observe that AMX-based algorithms outperform the IFMA-based ones for matrices of dimensions about 128 and higher. For many matrix dimensions and bit precisions, our sequential implementations outperform the current reference libraries Flint and NTL.

Our preliminary multi-threaded implementation scales well up to approximately 32 threads for the IFMA-based algorithms. For 64 threads, it might be an issue that our computer is equipped with two physical CPUs. The multi-threaded versions of our AMX-based algorithms often lose a factor two with respect to the theoretically maximal parallel speed-up. It seems that the performance is penalized by memory operations, and further work is needed to better approach the theoretical peak performance.

Note on AI usage

For the present revised version of the paper, we improved and completed the AMX-accelerated CRT-based matrix multiplication kernel. Claude Opus 5 was used as an assistant for this work, as detailed in Section 6.4.

2Preliminaries

2.1Notations

Throughout this paper, we write for the set of nonnegative integers, and for the set of unsigned -bit integers.

Given a ring and , we write for the set of vectors of length with entries in . We will also write for the set of polynomials of degree in over . We will assume the natural dense representation for such polynomials by their vectors of coefficients.

We formally extend the polynomial notation to the case when is replaced by an integer radix , typically of the form for some . This allows us to rewrite unsigned integers in as “polynomials” in the radix. Following terminology from GMP [29], it is customary to call the coefficients limbs. If is a machine integer type, then we may represent an integer by its vector of limbs.

For dimensions , we also write for the set of matrices with entries in . Unless stated otherwise, we will always assume that matrices are represented in row-major order, as vectors . By “ matrix multiplication”, we understand the multiplication of an matrix with an matrix into an matrix.

2.2SIMD arithmetic

Let be the machine precision (or one of the available machine precisions). Standard arithmetic operations on naturally extend to SIMD vectors of width in . Whenever a function that only involves basic arithmetic operations needs to be evaluated exactly times, this can be done replacing all operations by their SIMD versions.

For instance, given a matrix and a vector , consider the function that computes the matrix-vector product . If we need to evaluate this function times, then the simplest solution is to present the input as a matrix and a vector , in which case the exact same formula yields all matrix-vector products in SIMD format.

We call this way of vectorizing the pure SIMD mode. This works best whenever a given function needs to be evaluated a multiple of times. This also assumes that all data types are replaced functorially by their standard vectorized versions. For instance, the standard vectorization of the matrix type is . The standard vectorization of a multiple precision integer type consists of working with “limbs” that are vectors in . Hence SIMD multiple precision integers are represented as with and . Note that this may not seem to be the most natural choice since one might prefer to represent as with .

SIMD arithmetic can often be used as well for evaluating a function only once. For instance, the multiplication of two matrices and can benefit from SIMD acceleration whenever is a multiple of the SIMD width . Indeed, it suffices to reinterpret as a matrix and as a matrix , while casting to a matrix by repeating each of its entries times. After that, we can compute in the pure SIMD mode. However, these reinterpretations require some gymnastics at least, and exploiting SIMD arithmetic becomes even more delicate when none of the dimensions are a multiple of . We will come back to this tinkered SIMD mode in Section 4.3 below.

2.3Benchmarks

All experiments were conducted on a system equipped with two Intel Xeon Gold 6538Y+ (Emerald Rapids) processors (each with 32 cores running at a base frequency of 2.2 GHz and maximum frequency of 4GHz, 2 MB of L2 cache per core, and a 60 MB shared L3 cache). The system was running Ubuntu 24.04.4 with Linux 6.8.0 and we compiled our code with gcc 13.3.0.

We use the rdtsc instruction for measuring time in terms of TSC cycles. For the relative timings with respect to theoretical single-core peak performance, we multiply the number of elapsed TSC cycles by the ratio between the maximum core frequency and the TSC frequency to obtain an upper bound on the number of core cycles. (Compared to disabling dynamic frequency scaling, this avoids artificially slowing down the processor relative to memory accesses.) In the case of absolute timings, we directly convert the number of TSC cycles to a duration in seconds.

The reported values are means or median figures over a number of runs depending on the problem scale. For full details, we refer to the source code of the programs used for producing our timings tables, which are distributed with our implementations. We stress however that the timings provided in this article are meant to provide a first validation of each implementation strategy as well as some elements of comparison between them but not intended as definitive performance figures.

3Multiple precision arithmetic

3.1The standard representation

Let be our favorite (and supported) machine bit precision for integer arithmetic, like or . The standard machine representation of an integer is the vector of “its” limbs.

The hardware is assumed to provide an instruction for multiplying two limbs , with a result in (sometimes, the low and high parts need to be computed separately). For small precisions, multiple precision libraries typically implement integer multiplications by cross multiplying all limbs and summing the results while handling carries appropriately.

3.2Redundant representations

Unfortunately, carry propagation is unpleasant to vectorize and SIMD units typically do not provide “vector carry registers” (except for certain GPUs that handle this via mask registers). In a SIMD world, it is therefore more convenient to work with redundant representations: we still represent integers as polynomials in some radix , but the limbs are stored with a higher precision (which may be different for multiplicands and products). If , then and the representation is said to be normalized. We will call the bit precision and the bit capacity. The difference is the number of nail bits in the terminology of GMP.

On recent Intel Xeon processors there are two instructions and for integer FMAs (IFMA). Formally, given and , we have

where and denote the quotient and the remainder of the euclidean division of by . The IFMA representation is the redundant representation with and . Given a two-limb integer , we can accumulate the product of two limbs using one and one . As long as , this computation does not overflow, but typically does not remain normalized.

On even more recent Intel Xeon processors with AMX support, there is an instruction for matricial FMAs , where , , and . For individual integer coefficients of the product, this corresponds to working with a redundant AMX representation with and . Note that, with our convention, products of limbs are two limbs long for the IFMA representation, but only one limb long for the AMX representation. Indeed, when using AMX, products and multiplicands are represented using limbs of different -bit and -bit capacities, respectively.

3.3Naive multiple precision multiplication

Redundant representations have the advantage that the nail bits can be used to accumulate carries, while postponing normalization to the very end.

Consider for instance the multiplication of and using integer FMAs, where and . Assume also that . Then we first compute the non-normalized product with , and then normalize it, as follows:

This algorithm has been represented graphically in Figure 1; the lower and upper triangles correspond to and operations, respectively.

If , then the mere highest limbs of the result can be computed similarly. With the naive algorithm, the right picture in Figure 1 shows that this can be done at roughly half the cost of a full multiplication. This is very useful for fixed point arithmetic, by representing unsigned -limb fixed point numbers as -limb integers divided by . It is interesting to note that the nail bits can be exploited to represent small integer parts in this case.

However, our current implementation is limited to unsigned fixed point numbers. Signed variants can be obtained by representing signed integers by and in , so that . Alternatively, one may use a two's-complement representation for the highest limb and adapt multiple precision multiplication accordingly. But it becomes harder to exploit the extra nail bits in both cases.

Figure 1. Illustration of a full multiple precision integer product (on the left) and a truncated high product (on the right) using integer FMAs (the lower yellow triangles and the upper orange triangles correspond to and , respectively).

3.4Vector and matrix limbs

We briefly discussed pure SIMD vectorization of multiple precision integers in Section 2.2. The idea to represent such integers as polynomials in with vectorial limbs also works for redundant representations. This allows us to apply the algorithms from the previous subsection directly to pure SIMD multiple precision integers in .

In a similar way, matrices of multiple precision integers in can be rewritten into polynomials in with matrix limbs in . Computing a non-normalized product with , , and can then be done naively using AMX-accelerated matrix multiplications . Note that we do not compute low and high parts separately in this case. When doing high -limb multiplication in this way, this means that we need to do multiplications, which is slightly more than one half of the number for full products.

3.5Karatsuba multiplication

For large matrix multiplications with -limb coefficients, using a naive algorithm for polynomial multiplication becomes suboptimal. If is small, then it becomes better to use Karatsuba's algorithm [44].

Given a ring , the traditional Karatsuba trick reduces one multiplication with into three multiplications in and several additions and subtractions using the formulas , , and . This trick generalizes to the case where is replaced by and by , in which case the multiplication of two polynomials essentially reduces to three multiplications of polynomials in . These latter polynomials can be computed recursively using the same trick.

When working with multiple precision integers one technical difficulty is that and are in and not necessarily in . This makes it better to work with limbs of bit-size instead of . If our input integers used -bit limbs, then this requires them to be rewritten into the -bit limb representation and vice versa for the output. When using Karatsuba's algorithm recursively for -limb integers, we may either directly use limbs of bit-size (which yields a multiplication algorithm for -bit integers), or do the limb rewritings recursively as well (which yields a multiplication algorithm for -bit integers).

Karatsuba's trick can also be applied to matrices . In that case, we need to do multiplications in , but only additions and subtractions, which becomes negligible for large . Asymptotically, we thus gain a factor with respect to naive integer matrix multiplication.

Toom generalized Karatsuba's trick [65] and gave a way to reduce the multiplication of two polynomials in to multiplications in at the cost of an increased loss of bit-precision for the integer variant. Interestingly, Karatsuba's trick can also be used directly for , while losing a single bit of precision.

The idea (see also [37, Section 4.4.1], [45], [9, Exercise 1.4], and [60]) is to observe that the contribution of to can be computed using , for all . In this way, we need to first compute the products and then further products of the form , which yields a total number of products. If we just want the high coefficients , then we only need to do multiplications. This is just as good as Toom's algorithm for and almost as good for and .

3.6Chinese remaindering

A well known method for multiplying two matrices is based on the Chinese remainder theorem. We first pick pairwise coprime moduli . Best is to choose them as large as possible while fitting into bits. Let and assume that for some small with .

In order to compute , we first reduce and modulo for . There is an elegant way to evaluate the map using a vector-matrix product [14]:

(1)

Taking the limbs of the entries of as our rows , the computation of for thus essentially reduces to an matrix product over and similarly for . We next compute

for . We finally recover from by applying the Chinese remainder theorem. Consider the map with and for . We have

where is the inverse of modulo for . After precomputing the -adic expansions , it is again possible to compute using a vector-matrix product [14]:

			(2)

In a similar way as above, the recovery of from thus essentially boils down to an matrix product over . In order to ease the computation of the remainder with respect to , it is sometimes preferred to first compute the remainders of modulo and then to work with the -adic expansions of instead of . In this way, the quotient is always bounded by .

3.7Summary and asymptotic complexity

Each of the schemes that we described so far can be regarded as a way to reduce the multiplication of two matrices in to multiplications of matrices in , through the use of an appropriate encoding:

(3)

The cost of the reduction may depend on , but scales linearly with the size of the matrices. Since the cost of multiplications of matrices in is proportional to , the cost of the reduction is negligible for large . For large matrices, it is therefore best to select the scheme for which is minimal. The encoding/decoding process may also induce a loss of precision, which should be kept minimal as well.

Table 1 summarizes the factors and the loss of precision for the methods from the previous sections. For our implementations, we only considered the naive method, Karatsuba's method, and Chinese remaindering. In certain cases, Toom's algorithm or FFT-based algorithms may also be of interest.

Method 1 2 3 4 5 6 8 10 12 16 Bit loss Naive 1 4 9 16 25 36 64 100 144 256 0 Karatsuba 1 3 6 10 15 21 36 55 78 136 Karatsuba2 9 18 30 45 63 108 or CRT 2 4 6 8 10 12 16 20 24 32 Naive, high 1 3 6 10 15 21 36 55 78 136 Karatsuba, high 1 3 5 8 11 15 24 35 48 80 CRT, high 2 4 6 8 10 12 16 20 24 32 Naive, IFMA, high 1 4 9 16 25 36 64 100 144 256 0 Karatsuba, IFMA, high 1 4 8 13 19 26 43 64 89 151 CRT, IFMA, high 4 8 12 16 20 24 32 40 48 64

Table 1. Overview of the factor as a function of for various multiplication schemes, as well as the corresponding losses in bit precision. For CRTs, the bit loss is for or and , but for and . The first four rows correspond to full multiplication, whereas the other rows correspond to high multiplication. The last three rows count the number of IFMA operations that are required if low and high machine products need to be computed separately. For Karatsuba multiplication, we do not recurse, except for Karatsuba2, which recurses once.

4Multiplying small matrices using IFMA

In this section, we consider the problem of multiplying “small” integer matrices. Since AMX-acceleration is mainly interesting for quite large multiplications (of size at least ), we focus on SIMD-accelerated IFMA-based algorithms.

4.1Implementation strategy

Since we assume both and to be “small”, we will use naive algorithms for both matrix and polynomial multiplication. Here the qualifier “naive” is meant to capture the family of all algorithms that perform pairs of / integer FMAs in total. However, there are many such “naive” algorithms and the precise order in which we do perform the integer FMAs matters, since this crucially influences how well we can keep coefficients in registers and how much pressure we put on registers.

Instead of implementing various multiplication strategies directly in a language like C++, our strategy is to rely on the JIL library [2, 41] for computations and JIT compilation of straight-line programs (SLPs): we implement a large number of strategies and ways to combine strategies as symbolic programs. The JIL library provides highly efficient mechanisms to record the execution of these programs as SLPs and to JIT compile the resulting SLPs into highly efficient machine code. The recording plus compilation time is typically about 1000 times larger than a single execution of the generated code, which is one to three orders of magnitude faster than traditional compilers.

Let us illustrate this with a short example. We recall that an SLP is just a sequence of arithmetic operations. The operations belong to a fixed signature like and we operate on a finite number of data fields that are all of the same type. This type could either be a scalar type like or FP32, or a SIMD vector type like . A finite number of the data fields are designated to be “input” and “output” fields.

Now given a generic function operating on data fields of the same type, but which may involve loops and function calls, JIL provides a way to record the trace of its execution as an SLP. In Figure 2, we illustrated the result of recording a C++ template for matrix multiplication as an SLP.

When implementing multiplication algorithms in this way, the focus is not on the design of efficient C++ code templates, but rather on the design of clean templates, whose recording produces good SLPs. This makes it easy and efficient to build new strategies on top of other ones, like a multiplication algorithm for polynomials in on top of a multiplication algorithm for matrices in . JIL also provides standard mechanisms for composition and tensor products of maps, reindexation, etc. We refer to [ 41 ] for more information.

template<typename C> void matmul_222 (C* c, C* a, C* b) { for (int k=0; k<2; k++) for (int i=0; i<2; i++) for (int j=0; j<2; j++) if (j == 0) c[2*i+j] = a[2*i+k] * b[2*k+j]; else c[2*i+j]+= a[2*i+k] * b[2*k+j]; }

Figure 2. Illustration of the recording process of a program as an SLP. For the purpose of readability, we have chosen nice names for the data fields of the SLP. Inside JIL, the data fields are referred to by indices .

4.2Pure SIMD mode

Let us first consider the problem of multiplying integer matrices in the pure SIMD mode, using IFMA-based arithmetic. This means that our multiplicands are matrices in with coefficients that are SIMD multiple precision integers in , whose limbs are SIMD vectors in .

Two things matter for the design of efficient “naive” algorithms. We first should decide whether we view our multiplicands as matrices in with SIMD integer entries or as polynomials in with SIMD matrix “limbs”. The first point of view tends to be faster when is large with respect to , whereas the second one tends to be faster when is large with respect to .

Secondly, divide and conquer algorithms tend to reduce register spilling. For instance, matrix products can recursively be reduced to smaller matrix products by halving , , or . However, recursing too far down may increase register pressure. We implemented an ad hoc strategy that recurses until a certain threshold is reached. In the future, we plan to automatically generate many strategies of this kind and select the best one.

After implementing the desired multiplication strategy as a map in JIL, we record one evaluation of this map as an SLP, and JIT compile the SLP into machine code.

Remark 1. The SLP compilation framework allows us to physically work with a specific memory layout, say matrices with coefficients in , while applying an algorithm that works logically with another representation, like polynomials in with coefficients in . Indeed, the alternative representation just corresponds to a permutation of coefficients in memory. The implementation described in this section uses an abstraction provided by JIL (the aliased type in domain.hpp) to automate this kind of re-interpretations.

The timings in Tables 2, 3, and 4 were obtained by running the program mat_mul_bench in the JIL library (Git hash tag 999ee20) with the options –delay=5 –peak –scale=0.55.

Table 2 shows the ratios of our timings for an limb long matrix multiplication divided by the theoretical time to perform integer FMAs over (in theory, our processor has a throughput of two integer FMAs per cycle). The left-hand table shows the ratios for the standard representation; the right-hand one concerns the alternative representation of our multiplicands as polynomials in the radix with matrix coefficients.

We also recall that the generated codelet unrolls the entire computation. For , this means that the generated program performs IFMAs. The codelet therefore takes more than 12 MB in memory and does not even fit in the L2 cache. This explains the dramatic increase in the ratios for large and/or .

Multiplicands in \ 1 2 3 4 6 8 12 16 1 6.55 2.41 2.09 1.86 1.61 1.47 1.41 1.88 2 1.78 1.70 1.60 1.64 1.67 1.53 1.53 2.56 3 1.54 2.02 1.91 1.81 1.65 1.95 2.23 2.78 4 1.79 2.09 1.94 1.83 2.18 1.97 2.43 2.93 6 1.51 2.18 2.74 2.58 2.37 2.08 2.29 6.13 8 1.71 3.21 3.09 2.78 2.30 2.03 5.69 9.99 12 1.92 3.16 2.93 2.62 5.17 5.85 8.08 10.1 16 2.48 3.41 3.92 7.07 7.60 6.96 8.10 10.6 Multiplicands in \ 1 2 3 4 6 8 12 16 1 6.54 2.39 2.07 1.86 1.61 1.47 1.41 1.88 2 1.78 1.70 1.56 1.60 1.65 1.69 2.25 2.60 3 1.54 1.76 2.05 2.11 2.03 2.66 2.73 2.79 4 2.15 2.16 2.28 2.22 3.09 3.08 3.30 3.22 6 1.83 1.70 2.42 2.64 2.68 2.72 2.74 3.35 8 2.12 3.18 3.27 3.30 3.28 3.85 6.58 9.80 12 2.33 2.57 2.51 2.48 4.53 6.03 7.38 7.42 16 2.99 3.16 3.87 6.45 9.35 9.18 9.20 9.41

Table 2. Pure SIMD integer matrix multiplication using a naive algorithm and IFMA. We show the ratios with respect to the theoretical peak performance of two IFMAs per cycle.

4.3Tinkered SIMD mode

Let us now turn to the computation of a single matrix product over . In Section 2.2, we already described a reduction to the computation of a matrix product over , whenever is divisible by . We may then use an algorithm in pure SIMD mode to compute this product.

In fact the situation is even a bit better than that: instead of casting into a matrix , it is better to regard directly as an SIMD matrix over , and to replicate individual coefficients of across lanes whenever needed. This more compact representation of has the advantage that we may typically keep wholly in registers, after which obtaining coefficients from does not involve memory accesses.

If is not divisible by , then more work is required to benefit from SIMD acceleration: for an product such that is divisible by , we may write , , and for some with . We next redundantly encode , , and by matrices , , by taking

(4)

where , , and , , . Every SIMD multiplication in then computes a matrix product with coefficients in . If , then some of the lanes of need to be summed in order to retrieve .

The standard SLP signature in JIL includes instructions for permutations of lanes for SIMD base types (mathematically speaking, really maps on lanes, because replication of lanes is also allowed). Consequently, the tinkered SIMD multiplication algorithms can still be recorded into SLPs and compiled within the JIL framework. As a consequence of Remark 1, we note in addition that any of the pure SIMD multiplication strategies can naturally be used in conjunction with any of the tinkered encoding strategies.

Table 3 shows timings that we obtained, with similar conventions as for Table 2 . This time, we divided the total number of cycles by instead of .

Multiplicands in \ 1 2 3 4 6 8 12 16 2 9.31 4.95 4.24 3.54 3.89 3.17 2.56 2.51 4 3.27 2.63 2.55 2.33 2.25 2.03 2.03 3.20 6 3.44 3.41 3.17 2.77 2.99 2.92 3.62 3.69 8 2.42 2.21 2.32 2.19 2.74 2.39 2.87 3.42 12 1.97 2.53 3.13 2.75 2.42 2.16 2.48 6.91 16 1.77 3.36 3.21 2.93 2.48 2.14 6.06 10.3 24 1.84 3.56 3.15 2.78 4.57 6.25 8.49 10.4 32 2.71 3.63 4.81 6.65 7.99 7.36 8.77 11.6 Multiplicands in \ 1 2 3 4 6 8 12 16 2 7.84 4.27 3.65 3.03 3.39 2.63 2.17 2.12 4 2.76 2.22 2.17 1.95 1.71 1.59 1.43 2.23 6 2.92 2.43 2.65 2.40 3.50 3.22 3.41 3.14 8 2.05 1.89 1.86 1.67 2.70 2.59 2.60 2.50 12 1.74 1.70 2.65 2.52 2.24 2.14 2.11 2.68 16 1.58 2.52 2.43 2.22 2.17 2.06 4.36 7.78 24 1.85 2.31 2.28 2.18 3.89 7.24 6.62 6.71 32 2.79 3.07 3.34 8.22 9.87 9.57 9.88 10.4

Table 3. Tinkered SIMD integer matrix multiplication using a naive algorithm and IFMA. We show the ratios with respect to the theoretical peak performance of two IFMAs per cycle.

4.4Fixed point entries

One important application of integer matrix multiplication is linear algebra over fixed point numbers. In that case, we only need to compute the highest limbs of the product, which saves about half of the computations. Table 4 reports on the obtained timings. We recall that a shortcoming of our implementation is that it is currently limited to unsigned fixed point numbers.

\ 1 2 3 4 6 8 12 16 2 4.73 3.30 3.37 3.28 2.42 1.98 4 3.47 2.76 2.45 2.13 2.15 1.98 1.74 1.61 6 2.57 1.99 2.19 1.85 2.60 2.23 8 2.60 2.23 2.10 1.91 1.82 1.89 1.87 1.72 12 2.10 1.76 1.96 2.38 2.41 2.18 2.07 1.94 16 1.93 1.77 2.37 2.26 2.20 2.06 1.98 2.31 24 1.72 2.12 2.50 2.40 2.26 2.52 6.02 5.96 32 2.10 2.89 2.58 2.51 3.06 6.01 6.26 5.98

Table 4. Tinkered SIMD matrix multiplication for limb fixed point entries and using a naive algorithm and IFMA. We show the ratios with respect to the theoretical peak performance of two IFMAs per cycle.

5Multiplying large matrices using IFMA

5.1Implementation strategy

The simplest way to accelerate a large matrix multiplication via JIL is to consider , , and as block matrices. We use JIL for the compilation of efficient SLPs to multiply or multiply-accumulate blocks. The block matrices themselves are multiplied using conventional C++ code. The last possibly incomplete row and column are treated using zero-padding or recursive decomposition. The block sizes should be carefully chosen and may differ for , , and . By taking sufficiently large blocks, we intend to make the overhead of the C++ part of the code as small as possible, so that the total cost essentially boils down to multiplications of blocks.

One goal of the above strategy is that JIL can focus on compiling SLPs of a fixed size and that we can use a higher level language like C++ for the rest, with minimal overhead. However, for practical block sizes, the overhead of the C++ code is small, but not necessarily negligible, which makes it desirable to further improve the interface between JIT compiled SLPs and their use within high level C++ programs. For the purposes of this paper, we therefore extended JIL with a mechanism (called “promises”, but unrelated to the classical concurrency primitive of the same name) that allows it to compile slightly more complex programs than mere SLPs.

For instance, a common scenario is to execute the same SLP on a vector of successive input and output chunks. This happens for example when encoding all entries of an integer matrix in into CRT format. Our implementation of this step uses a “loop promise” whose compilation transparently rewrites the SLP into a new one that can handle entries at a time in SIMD fashion (this requires transpositions in memory to work in the SIMD representation), or even a multiple of entries at a time via unrolling. In the future, we also plan to automatically determine and factor out loop invariants.

Back to matrix multiplication, another common scenario is that the input of the SLP has first to be gathered from memory and/or that the output needs to be scattered to memory, as we explained in Section 3.7 and (3). For instance, Chinese remaindering essentially reduces the multiplication of two matrices in to multiplications of matrices in . Now it is more natural to CRT encode an input matrix in as a vector of matrices in rather than a matrix with vector entries in . The latter case corresponds to the loop scenario that we described above, whereas the first one requires us to submit the output to an additional matrix transposition over . In order to reduce unnecessary memory accesses, we preferred to introduce a new scenario that scatters entries in directly into the desired matrices while processing the CRTs.

For the present revised version of our paper, we further optimized the critical loop of the main IFMA-based kernel using AVX512 intrinsics. This deviation from our original design strategy allowed us to test some ideas more quickly than with a pure JIL-based implementation. Through these experiments, we also identified a few features that are currently missing in JIL. We plan to extend JIL accordingly in the future and hope that this will further improve our timings.

5.2Experiments

We have tested most of the strategies of Table 1 for the multiplication of large matrices using integer FMAs. For the “Multiply” step of (3), we use the blocking technique that we described above, and also a tinkered SIMD as in Section 4.3 for the multiplication of blocks. We report on relative timings with respect to the theoretical peak performance. Tables with absolute timings are given in Appendix A.3. The timings were obtained using the command matmul_ifma_bench –threads=1 –peak –scale=0.55.

In Table 6, we consider long limb multiplications of matrices for the naive and the Chinese remaindering (CRT) strategies. More precisely, we show the number of cycles that it takes to compute the matrix product, divided by , the theoretical number of cycles to perform all IFMA operations in the multiplication step of (3). In Table 5, we show timings for the non-recursive variant of Karatsuba's algorithm, both in the case of long limb multiplication and high limb multiplication.

The tables show that we are not far from the theoretical peak performance. The factors from Table 1 can therefore guide the design of efficient algorithms when gets large. In particular, it is interesting to note that Karatsuba's algorithm is most efficient for high multiplication until five limbs.

1 2 3 4 6 8 12 16 32 3.12 2.42 2.41 2.34 2.34 2.54 2.83 3.29 64 1.91 1.63 1.67 1.77 1.92 2.00 2.09 5.89 96 1.58 1.62 1.70 1.79 1.80 1.81 4.24 4.72 128 1.45 1.60 1.65 1.68 1.63 3.44 3.88 4.66 192 1.52 1.52 1.49 1.47 2.66 3.02 3.57 3.81 256 1.46 1.43 1.42 2.32 2.62 2.93 3.06 3.19 384 1.50 1.47 1.47 1.86 2.90 3.16 3.53 3.76 512 1.46 1.39 1.37 2.41 2.56 2.86 3.04 3.14 768 1.44 1.38 1.38 2.40 2.55 2.83 3.02 3.13 1024 1.44 1.37 1.34 2.40 2.54 2.83 3.00 3.15 1536 1.43 1.37 1.34 2.40 2.53 2.81 2.99 3.13 2048 1.42 1.37 1.34 2.39 2.51 2.80 2.98 3.12 1 2 3 4 6 8 12 16 32 4.58 2.71 2.27 2.35 2.32 2.37 2.66 2.87 64 2.47 1.76 1.52 1.69 1.86 1.98 2.06 2.09 96 1.91 1.70 1.50 1.69 1.81 1.82 1.80 4.31 128 1.68 1.70 1.52 1.62 1.64 1.63 3.51 3.95 192 1.66 1.61 1.38 1.43 1.42 2.78 3.27 3.60 256 1.64 1.50 1.28 1.32 2.38 2.70 3.03 3.09 384 1.80 1.61 1.38 1.44 2.52 3.03 3.37 3.59 512 1.70 1.44 1.25 1.31 2.49 2.67 3.01 3.07 768 1.67 1.43 1.25 1.30 2.48 2.63 2.99 3.04 1024 1.67 1.43 1.24 1.29 2.47 2.65 2.97 3.05 1536 1.66 1.46 1.24 1.29 2.46 2.62 2.96 3.03 2048 1.65 1.44 1.25 1.29 2.46 2.60 2.95 3.02

Table 5. Relative timings for one long or one high limb matrix product (on the left-hand and right-hand sides, respectively) using IFMA and the non-recursive variant of Karatsuba's algorithm.

There is still quite some room for further improvements. Our implementation of the naive algorithm was done first and based on a suboptimal implementation of matrix transposition for the encoding and decoding steps; this explains the poor performance for small . In general, with more work, we expect that it should be possible to reach factors that are very close to (at least for, say, and ).

1 2 3 4 6 8 12 16 32 3.78 2.95 2.14 1.79 1.65 1.51 1.54 1.83 64 2.37 1.89 1.57 1.44 1.42 1.36 1.34 1.84 96 1.88 1.71 1.51 1.40 1.32 1.26 1.26 1.63 128 1.69 1.60 1.42 1.33 1.28 1.25 1.60 1.62 192 1.62 1.47 1.35 1.28 1.25 1.30 1.56 1.51 256 1.54 1.41 1.31 1.26 1.29 1.27 1.45 1.44 384 1.63 1.45 1.34 1.28 1.32 1.29 1.41 1.43 512 1.54 1.38 1.29 1.32 1.28 1.27 1.45 1.43 768 1.54 1.39 1.28 1.32 1.37 1.29 1.45 1.43 1024 1.53 1.37 1.28 1.31 1.37 1.28 1.45 1.42 1536 1.52 1.37 1.27 1.31 1.28 1.26 1.44 1.42 2048 1.52 1.36 1.26 1.31 1.28 1.26 1.44 1.42 1 2 3 4 6 8 12 16 32 2.64 2.69 2.90 3.15 3.69 4.24 5.94 8.21 64 1.74 1.80 1.95 2.18 2.58 2.95 3.92 4.77 96 1.49 1.69 1.84 1.98 2.26 2.54 3.07 3.62 128 1.47 1.62 1.74 1.84 2.01 2.21 2.57 2.98 192 1.41 1.49 1.54 1.59 1.68 1.80 2.13 4.05 256 1.35 1.39 1.43 1.45 1.56 2.87 3.37 3.84 384 1.41 1.45 1.52 1.56 1.67 2.24 2.94 4.16 512 1.35 1.37 1.40 1.43 1.64 2.79 3.21 3.63 768 1.36 1.36 1.40 1.42 1.78 2.77 3.21 3.60 1024 1.34 1.35 1.39 1.42 1.85 2.77 3.17 3.59 1536 1.34 1.35 1.41 1.44 1.61 2.76 3.17 3.56 2048 1.34 1.35 1.40 1.44 1.61 2.75 3.14 3.52

Table 6. Relative timings for one long limb matrix product using IFMA. The left-hand and right-hand sides respectively indicate the relative timings (with respect to theoretical peak performance) for the naive and the CRT strategies.

6Accelerating matrix multiplication using AMX

6.1Implementation strategy

The AMX unit can be configured for computing the matrix FMA where , , and for any , , and such that is a multiple of four. The processor can issue one such MFMA every 16 cycles when , and somewhat faster if . The total number of scalar FMAs per cycle is highest when using the maximal allowed dimensions , , and , which is the main configuration used in our implementation.

For , the matrices , , and need to be loaded from memory and stored in dedicated tile registers. There are eight available tile registers of 1 KB each and there are dedicated strided load and store operations for tiles. In addition, the second multiplicand is stored in a special VNNI format, which coincides neither with the row-major nor column-major formats.

Since the 1 KB size of tiles is very different from the 64 byte size of AVX512 registers, AMX and AVX512 instructions do not fit well together into the SLP framework of JIL. For this reason, we decided to start with a base C++ implementation of the matrix FMA (MFMA) for arbitrary sizes , using Intel intrinsics. In the future, we plan to extend JIL with special “promises” (supplementing those described in Section 5.1) for building code that gracefully mixes AMX and AVX512 instructions.

The core of our work on AMX is an efficient implementation of MFMA for bit entries. On top of it, we implemented multiple precision matrix multiplication using various schemes as in (3). The encoding and decoding stages are mostly done using AVX512 instructions and one of our main concerns will be to keep the cost of these stages reasonable. This is challenging due to the fact that the throughput of AVX512 instructions is 8 to 16 times lower than that of AMX instructions (for eight bit entries, we can do FMAs per 16 cycles with AMX versus 64 operations per cycle and per execution unit with AVX512). We again heavily rely on SLPs from JIL and a similar extra layer as described in Section 5.1 for loops.

6.2Int8 matrix FMAs

Our 8-bit MFMA kernel follows a relatively standard design informed by the classic paper by Goto and van de Geijn [28]; see also [67, Sec. 5.1] and the recommendations of the Intel Optimization Reference Manual [43, Chapter 20].

Like many implementations of general matrix multiplication on AMX units, it is based on a microkernel that multiplies a tile block from with a tile block from and accumulates the result in a tile block of , using all eight tile registers. Assuming that this microkernel is run repeatedly with the four accumulation tiles fixed, it requires about four loads and four tile FMAs per iteration. While, in ideal circumstances, the architecture we target can sustain up to two tile loads from L1 cache (and one store) per tile FMA, the pattern makes it easier to hide memory movements behind computations even when not all loads hit L1 cache.

The typical scenario we aim for is for to reside in L1 cache while resides in L2 cache, as it is still possible in that case to fully overlap the loads with the FMAs in the steady state. We access the matrix using non-temporal load instructions so as not to needlessly evict data from from L1 cache. In addition, we use the strided load instructions mentioned above to read the 16 rows making up a tile of from their original, typically non-contiguous memory locations without first packing them together, and similarly for writing tiles back to . Our experiments suggest that, for data aligned on a cache line and in the absence of cache conflicts, strided loads and stores are not significantly slower than contiguous ones (see however [51] for more on the limits of this strategy). We deal with odd numbers of tiles in either dimension using a variant of the same microkernel, and with incomplete - or -tiles using a temporary buffer and AVX512 masked loads and stores.

The microkernel is wrapped in loops over (from outermost to innermost) that multiply a series of tile blocks of each fitting in L1 cache by a single block of that fits in L2 cache. On a single core, we can typically process blocks of about tiles in this fashion. The resulting macrokernel is itself wrapped in a cache-blocking kernel, using the same loop order, that repacks roughly L2-sized blocks of into contiguous VNNI-encoded full tiles (padded with zeros if necessary) before processing them using the macrokernel. The encoder operates on blocks of up to entries of , corresponding to two adjacent tiles, using AVX512 instructions. Slightly surprisingly, we found this block format to perform better in our context than the alternative of entries, even though the latter makes it possible to load full cache lines (rows of the block, to be split among four tiles) at once. We speculate that this may be due to slightly different pressures on execution ports or differences between the speed of 256 and 512-bit loads from L1 cache.

We also implemented some microkernels dedicated to the special cases where one of the dimensions of the product is less than or equal to one tile.

Remark 2. Endo, Ohshima and Nanri [20] recently proposed an alternative microkernel that reserves six tile registers for accumulating a block of and uses the remaining two to hold the current tiles of and . The tiles of are loaded once per block and those of are loaded twice. While this results in an fma/load ratio of only , the idea is that two of the loads are virtually guaranteed to hit L1 cache, whereas, depending on the blocking strategy, all four of the tiles accessed by the kernel may have been evicted by the time they are reused. While Endo et al. report reaching better performance with their kernel than with the standard one, in our experiments, it performed worse except for . Endo et al.'s code does not appear to be publicly available, but a plausible explanation would be that the use of non-temporal tile loads makes their technique unnecessary.

6.3Naive multiple precision multiplications

The simplest approach to build an algorithm for -limb matrix multiplication on top of the work from the previous subsection is to rewrite input matrices as elements in . Then a long multiplication boils down to matrix FMAs and a short (low or high) multiplication to matrix FMAs.

The N and L columns in Table 7 show the factors between the performance of our implementation and the theoretical peak performance. These timings can be reproduced with the command matmul_alt_bench –peak –scale=0.55 –threads=no using the provided software distribution. As the variance of the timings obtained with AMX is quite large, with frequent outliers that can significantly alter mean timings, this table shows median timings over a large number of runs. If one is satisfied with a factor of 2 to 3, then our current implementation will do, but approaching the theoretical peak performance as well as we did for IFMA-based methods seems to require significantly more effort. In Table 21 in the appendix, we can see that the naive AMX-based implementation is competitive with the IFMA-based one for large matrices.

The cost of encoding and decoding using AVX512 is far from negligible. It seems hard to entirely get rid of this cost, since the algorithms from the previous subsection benefit from the fact that most of the matrices have a natural memory layout. Nonetheless, by interleaving the AMX instructions appropriately with the encoding and decoding, we expect that the factors can be reduced significantly. As increases, the cost of decoding and encoding should a priori become smaller with respect to the cost of the inner matrix multiplications. This is indeed what we observe for small values (say ). However, for larger values of , we start to suffer from cache misses. In the future, we plan to reduce these misses through a more careful interleaving of memory loads and stores with the actual computations. But again, part of the slow-down is probably unavoidable.

AMX N16 N24 N32 N48 N64 L16 L24 L32 L48 L64 64 2.71 12.9 9.78 7.97 6.39 5.56 12.1 8.75 7.22 6.10 5.31 96 4.43 9.69 8.72 7.85 7.06 6.73 8.87 7.94 7.40 6.89 6.66 128 3.79 6.99 6.25 5.82 5.37 5.15 6.53 5.80 5.47 5.31 5.20 192 2.74 4.61 4.18 4.08 4.11 4.35 4.28 4.17 4.16 4.56 4.41 256 2.34 3.84 3.77 4.03 4.01 3.89 3.78 4.29 4.19 4.00 3.81 384 1.75 4.02 3.71 3.55 3.30 3.15 4.09 3.74 3.51 3.24 3.08 512 1.80 4.26 3.40 3.19 2.93 2.82 3.81 3.35 3.12 2.90 2.77 768 1.82 2.87 2.62 2.49 2.36 2.31 2.83 2.57 2.44 2.36 3.98 1024 1.75 2.52 2.35 2.32 3.91 3.53 2.48 2.40 4.46 3.94 3.58 1536 2.05 2.60 2.59 4.28 3.86 3.67 2.72 4.84 4.44 3.86 3.51 2048 2.00 4.76 4.40 4.11 3.62 3.50 5.11 4.49 4.12 3.77 3.38 3072 1.96 4.24 3.87 3.70 3.32 3.06 4.84 4.09 3.72 3.23 2.97 4096 2.75 4.38 4.10 3.97 3.65 3.46 4.71 4.23 3.92 3.57 3.39

Table 7. Median factors with respect to theoretical peak performance for multiplying two matrices of unsigned integers, measured using the matmul_alt_bench executable provided with our software. The AMX column corresponds to the raw bit product from Section 6.2. The N columns correspond to a bit low product, and the L columns to a bit product, in both cases computed with the long multiplication strategy of Section 6.3.

6.4Chinese remaindering

We also implemented an AMX-accelerated version of the Chinese remaindering approach from section 3.6 in which both the modular products and the CRT transforms are based on AMX tile multiplication. As mentioned in the introduction, we used large language models, including in particular Claude Opus 5, to develop this kernel. More precisely, after deciding on a general implementation strategy (modular products using the AMX kernel presented above, CRT transforms by dedicated AMX-accelerated kernels in the style of [14], cache blocking), we used Claude Opus to prototype several variants of the algorithm, interactively refine certain details, and develop a first draft of the implementation we settled on. We then reviewed and significantly modified this first draft to obtain the version presented here.

Our implementation currently supports up to 32 moduli of bits. The transforms (1) and (2) are done in large batches by replacing the left-hand rows in (1), (2) by matrices with one row per integer being transformed. This amounts to the multiplication of a tall and skinny matrix by a fixed encoding or decoding matrix of more balanced dimensions. Thus, while the multiplication itself benefits from AMX acceleration, memory traffic and the remaining reduction steps dominate the cost. In addition, typical integer multiplication formats like bits lead to matrix products with an inner dimension significantly smaller than , which cannot use the full power of the AMX unit and are best performed using adapted tile configurations such as or .

In some more detail, our CRT-AMX kernel works as follows. The output matrix is divided into square blocks to be computed sequentially as the product of a full-width horizontal panel of and a full-height vertical panel of . We determined experimentally that a block side of 256 to 512 works well in practice but have no complete theoretical explanation for this choice. The modular reductions of each horizontal panel of are computed just before starting the work on the blocks of that depend on it. The modular reductions of , in contrast, are computed upfront for the whole matrix and immediately repacked into VNNI-encoded panels. In both cases, the multiplication by the encoding matrix yields a matrix of integers of at most 21 bits that need to be reduced modulo the 8-bit moduli and scattered in memory so that each reduced matrix is stored contiguously. This is done as the tiles of the product by the encoding matrix become available, using AVX512 single precision floating-point arithmetic for the 21-bit to 8-bit modular reduction. The VNNI-encoding step is a separate pass over blocks of .

For each block, we loop over the moduli and compute each modular product by a variant of the AMX Int8 matrix multiplication kernel of section 6.2. Compared to that kernel, the variant accepts a right operand pre-packed as a matrix of VNNI-encoded tiles instead of reorganizing it on the fly. In addition, it reduces each of the 32-bit entries of its output matrix modulo the current 8-bit modulus and narrows it to 8 bits before storing it. For the reduction of a -bit integer , we precompute , and , after which can be computed in a single AVX512-VNNI instruction.

As soon as all the modular reductions of a given block of have been computed, they are lifted to recover the corresponding integer entries. The inverse CRT map is implemented as a left multiplication by a fixed matrix, namely the transpose of the encoding matrix from (2), again followed by some postprocessing. An advantage of multiplying to the left is that the result comes out in a layout that is convenient for further processing using AVX512 instructions. The entries of the modular reduction of can also be read by full cache lines in their natural layout, but must be VNNI-encoded before they are multiplied by the decoding matrix. As in the direct CRT phase, the tiles holding the result of that multiplication are stored to a temporary buffer and postprocessed on the fly.

Writing for the product of all the moduli in use as in section 3.6, the postprocessing needed to recover one entry of the final output matrix consists in reducing modulo the sum where the are the entries of a given column of the matrix . This is done for 16 columns of in parallel, using vector instructions. Each is the dot product of two vectors of length at most 32 of 8-bit integers, so and we need to subtract from a multiple of with . We estimate from the most significant terms of the sum evaluated in floating-point arithmetic, then compute in -limb integer arithmetic with full carry propagation. If our estimated was correct, this yields the normalized representation of . In rare cases, the estimate will be off by one unit, but provided that is strictly larger than the bit length of (which can always be arranged), this can be detected from the result of the reduction and fixed. Finally, the matrix (represented as AVX512 vectors) of limbs obtained this way is transposed and each entry of can be written to the destination buffer.

Table 8 shows the performance of our implementation relative to theoretical peak performance. For matrices and moduli, we compute this peak performance as the cost in cycles of AMX matrix multiplications in dimensions ; in other words, the normalization factor does not take into account the cost of Chinese remaindering. We see that the measured overhead is large for small matrices but comes within a factor of two of that of the naive method from the previous section when and . In view of the relative costs of the underlying algorithms, this means that the CRT-AMX kernel outperforms the schoolbook multiplication kernel in these cases. This is confirmed by the absolute timings in Table 22 of the appendix, which also shows that the CRT-AMX kernel is faster than our IFMA-based implementation for comparable bit sizes when .

6 8 10 16 19 28 32 16 32 AMX C16 C24 C32 C55 C64 C96 C112 M124 M237 64 2.71 121. 94.8 59.4 57.2 51.4 39.9 45.3 40.2 45.4 96 4.43 79.3 62.3 38.8 37.4 33.8 26.7 31.1 26.5 32.4 128 3.79 59.3 46.7 28.8 28.2 25.5 21.4 24.6 19.9 25.8 192 2.74 39.8 31.7 20.1 20.4 18.5 14.8 16.9 14.3 17.3 256 2.34 31.1 25.6 16.3 16.1 14.6 11.9 13.3 11.7 13.6 384 1.75 21.5 17.9 11.6 11.3 10.4 8.56 9.60 8.48 9.82 512 1.80 16.7 14.0 9.07 8.84 8.15 6.85 7.52 6.76 7.77 768 1.82 11.9 9.72 6.53 6.53 6.06 5.21 5.84 5.21 8.95 1024 1.75 9.21 7.63 5.38 5.32 5.04 5.97 6.63 4.46 7.39 1536 2.05 7.33 6.13 4.63 5.70 6.49 6.07 6.58 4.94 7.50 2048 2.00 6.56 6.89 5.43 5.65 5.48 5.13 5.54 5.14 6.14 3072 1.96 5.35 5.18 4.36 4.44 4.33 4.17 4.45 4.09 4.69 4096 2.75 5.45 5.00 4.40 4.47 4.37 4.25 4.43 4.13 4.63

Table 8. This table complements Table 7 with timings for the CRT-based algorithms. We first repeated the column with the raw AMX times. The C columns correspond to a bit multiplication (that is, multimodular multiplication of matrices with -bit entries using enough moduli for the integer result to be exact for all considered matrix dimensions). The M columns correspond to multiplication modulo a fixed -bit integer. Above these columns, we indicate the number of -bit moduli that we used.

Again, we have not investigated in detail how to best schedule the various stages of the algorithm to limit memory transfers as much as possible and overlap work using different CPU resources, but we expect that this could lead to some speedups for small matrices. Finally, as mentioned, our implementation is limited to 32 moduli, allowing the computation of results with entries of up to 237 bits. The same design could be extended to at most 54 coprime moduli and 362 bits, after which we run out of 8-bit primes. For larger bit-sizes, we plan to resort to 16-bit primes, for which we pay a factor two, or to 14-bit primes with Karatsuba multiplication, which still leads to an overhead of .

7Multi-threaded matrix multiplication

7.1Implementation strategy and limitations

In order to measure the interest of IFMA and AMX instructions for matrix multiplication in the cleanest possible way, it is natural to first focus on their behavior on a single core, as we did in most of this paper. However, for real-world applications, we will also need multi-threaded versions, and we should examine how well the proposed strategies scale in that context.

In this section, we describe our preliminary multi-threaded implementation for the multiplication of very large integer matrices. The strategy is simple: for an matrix multiplication or matrix FMA , we subdivide the matrix in blocks and (roughly) use threads to execute matrix multiplications in parallel. Each multiplication is itself subdivided into multiplications, where and are chosen such that the destination matrix and the and multiplicands fit into the local L2 cache. In fact, our implementation copies these basic blocks into hot thread-local memory, as this turned out to be more efficient than working directly on blocks of the original matrices.

For one large matrix multiplication, reading the and blocks from global memory requires about coefficients to be transferred. These blocks typically do not fit into the L3 cache, so they really come from global memory. On the other hand, the design requires the output coefficients to be written only once to the global memory. In our implementation, all blocks are transferred using massive block matrix copy operations, and we usually pay the full price for cache misses. If the bandwidth of the global memory is small, then our implementation will become memory bound and the read operations will dominate the total cost.

Our strategy is a compromise between performance and ease of implementation. It has the advantage that the base multiplications can be done directly using any abstract single-threaded implementation. However, if the number of available threads becomes large, then it would be better to interleave the block transfers with the actual computations. But this improvement requires the single-threaded implementation to be adapted accordingly to allow for this kind of interleaving. For our first implementation, we have decided to ignore this intricate optimization, which doubles the efficiency at best.

It should be noted that the multiplication of small or medium size integer matrices may also be of interest in certain HPC applications in computer algebra, like sparse interpolation [58] or numerical homotopy continuation [62]. In that case, parallelization is trivial (although memory issues may still play a role), and it is relevant to optimize single-threaded algorithms, as we did in this work.

7.2Experiments

In Table 9, we give the speed-ups for naive matrix multiplication using IFMA instructions for various numbers of threads, for limb long multiplication with . Table 10 is similar, but for AMX instructions. For threads, our CPU runs at its maximal boosted frequency of 4 GHz. For threads, the frequency drops to 3.92, 3.42, and 2.65 GHz, allowing for theoretically maximal parallel speed-ups of 15.7, 27.4, and 42.4, respectively.

At present, we only implemented multi-threaded versions of the naive multi-limb algorithms. For variants based on Karatsuba multiplication and CRTs, it is natural to distribute the transformed products among the available cores, but to obtain a significant speedup this way it is necessary to parallelize the encoding and decoding phases as well.

2 4 8 16 32 64 256 1.9 3.0 2.9 3.3 4.4 4.4 512 1.9 3.6 6.2 10.8 9.1 15.2 1024 2.0 3.7 7.0 13.2 21.0 23.5 2048 2.1 3.7 7.2 13.9 23.6 27.7 4096 2.0 3.7 7.3 14.0 23.0 29.3 8192 2.0 3.7 7.6 14.5 25.2 25.5 16384 2.0 3.8 7.4 14.9 25.0 26.8 32768 2.0 3.8 7.5 14.5 25.2 30.5 bits 2 4 8 16 32 64 256 1.9 3.2 5.4 8.9 8.4 8.4 512 2.0 3.1 6.7 12.2 17.1 11.4 1024 1.9 3.1 6.9 11.1 19.7 12.8 2048 1.9 3.2 7.2 12.0 18.9 23.9 4096 2.0 3.2 7.3 13.1 19.8 25.6 8192 2.0 3.3 7.8 14.1 21.3 27.1 16384 2.0 3.9 7.7 14.2 21.1 26.9 bits 2 4 8 16 32 64 256 1.3 2.0 3.5 5.1 4.9 4.6 512 1.9 3.2 7.4 14.6 22.1 27.8 1024 2.0 2.6 6.9 12.1 20.1 26.0 2048 2.1 2.6 7.4 12.6 21.3 19.1 4096 2.1 2.6 7.5 13.6 21.1 20.4 8192 2.1 2.6 7.3 13.7 20.6 19.7 16384 1.3 2.4 6.2 8.6 20.1 19.4 bits

Table 9. Speed-ups for naive multi-threaded matrix multiplication using IFMA instructions and threads. The three tables show the parallel speed-ups for limb long multiplication with and , respectively.

2 4 8 16 32 64 1024 0.7 0.8 1.3 2.1 3.1 3.5 2048 1.0 1.4 2.5 4.6 6.8 10.0 4096 1.3 1.8 3.6 7.1 12.4 19.7 8192 1.4 1.9 3.7 7.4 13.3 21.4 16384 1.4 1.8 3.6 7.0 13.6 20.9 32768 1.5 2.1 3.7 7.3 13.7 20.8 65536 1.6 2.5 4.5 8.5 15.2 23.6 bits 2 4 8 16 32 64 1024 1.8 2.7 4.7 8.0 11.4 15.2 2048 1.7 5.0 7.7 10.7 15.6 15.8 4096 1.7 2.5 4.6 12.7 15.1 18.4 8192 1.7 2.9 4.6 5.9 6.7 6.7 16384 1.9 3.3 5.5 7.1 7.6 6.1 32768 1.9 3.7 6.7 10.5 12.5 11.5 bits 2 4 8 16 32 64 1024 3.0 5.2 9.9 17.3 23.7 29.6 2048 1.8 3.0 7.0 10.9 17.4 18.8 4096 1.7 3.0 5.5 7.9 7.3 21.2 8192 1.8 3.4 5.5 7.6 8.8 7.0 16384 2.0 3.9 7.2 10.8 13.6 11.3 32768 1.9 3.2 7.4 13.5 21.0 23.4 bits

Table 10. Speed-ups for naive multi-threaded matrix multiplication using AMX instructions and threads. The three tables show the parallel speed-ups for limb long multiplication with and , respectively.

8Conclusion and perspectives

We have described a first implementation of accelerated integer matrix multiplication using Intel's IFMA and AMX instructions. Most of the theory behind our work is well known, but we carefully describe how to adapt known algorithms in order to exploit the IFMA and AMX instructions. For matrix multiplication with -bit coefficients, the best strategy and its efficiency greatly depend on , , , and . For many regimes, our single-threaded implementation is within a factor 1.5 (IFMA) or 3 (AMX) of the theoretical peak performance of the best known multiplication strategy.

The IFMA-based variant is faster than state-of-the-art implementations based on GMP or double precision BLAS in almost all regimes, often by a substantial margin. In particular, it provides a good base strategy whenever IFMA instructions are available. Our implementation is very efficient for small matrices as well, thanks to the use of dedicated codelets for various sizes that are generated via JIL. The performance of the preliminary multi-threaded version scales reasonably well up to about or threads, but more work is needed to reach the theoretical peak performance for high numbers of threads. Nonetheless, we are confident that it should eventually be possible to make it as efficient as traditional double precision BLAS libraries, since IFMA is really a variant of the usual FMA instruction.

On a single core, our AMX-based implementation further improves upon the IFMA-based one for large matrices (e.g. with best performance if , , and are multiples of ), even when the behavior of the former is close to optimal. It also provides much finer control over the bit-precision . However, obtaining these results required a good deal of trial and error. The extremely high efficiency of the AMX instruction also requires the asymptotically subdominant AVX512 parts of the algorithms to be implemented with care. AMX-based implementations finally tend to become memory-bound for a lower number of threads, which makes multi-threading challenging.

Our code is still work in progress and there are multiple directions for future improvements. From a software engineering point of view, our ultimate aim is to provide a mixed-precision BLAS library that can easily be integrated into other software. For this, a better high-level interface is needed which automatically selects the best strategy for every , after an auto-tuning phase for the specific machine. We believe that there is room for further lowering the overhead of our AMX-based multiplication kernels by experimenting with other strategies and local optimizations. Last, but not least, we plan to extend the CRT-AMX kernel to 14 and 16 bit moduli, which will allow us to better take advantage of AMX accelerators for bit precisions .

Bibliography

AAbsolute timings

As motivated in the introduction, the timings that we reported so far are relative to the theoretical peak performance. This appendix provides some of the corresponding absolute timings. We also added absolute reference timings for existing software.

Of course, the running times are typically best for bit lengths that are multiples of a limb's bit-length, and our present implementation is limited to such bit lengths. The limb size is typically 64 bits for traditional implementations, but 52 or 51 bits for our implementation based on the IFMA instruction. The tables for existing implementations therefore show the timings for multiples of both 52 and 64 bits and they need to be interpreted with the necessary care.

A.1Pre-existing software

Tables 11 to 16 present some timings obtained with Flint [ 33 ], NTL [ 61 ], and Fflas [ 18 ], three well-known algebraic computation libraries. We use the last development versions available at the time we ran our experiment, namely revision 1e74b110 of the Flint source tree, NTL release 11.6.0, and Fflas revision 06aec3f8 installed as part of LinBox revision bb590e22 . Flint and Fflas are configured to link with OpenBlas (revision 9e857d00 ). The code we use to build the external libraries and measure their running times is available at

Since modular matrix multiplication has often benefited from more optimization effort than plain integer matrix multiplication, we include implementations working over or in addition to implementations working over .

\ 8 16 24 32 52 64 104 128 208 256 416 512 832 8 0.41μs 0.41μs 0.41μs 0.61μs 0.60μs 2.05μs 2.07μs 6.21μs 8.51μs 8.34μs 18.4μs 28.6μs 53.7μs 16 1.66μs 1.65μs 1.63μs 3.52μs 3.91μs 12.6μs 13.0μs 65.3μs 86.7μs 96.6μs 0.18ms 0.23ms 0.36ms 32 9.21μs 9.36μs 9.24μs 28.9μs 29.7μs 87.6μs 87.0μs 0.29ms 0.46ms 0.53ms 0.92ms 2.11ms 3.00ms 64 21.3μs 19.7μs 61.4μs 0.20ms 0.20ms 0.64ms 0.63ms 1.43ms 2.46ms 2.61ms 4.50ms 5.52ms 10.3ms 96 92.1μs 90.6μs 0.33ms 0.68ms 0.66ms 2.06ms 2.01ms 4.03ms 6.39ms 7.21ms 12.1ms 14.9ms 27.1ms 128 0.20ms 0.19ms 0.79ms 1.52ms 1.50ms 4.77ms 4.75ms 8.70ms 13.4ms 15.1ms 26.2ms 31.8ms 56.9ms 192 0.51ms 0.51ms 2.15ms 3.52ms 3.47ms 14.2ms 15.5ms 25.3ms 39.9ms 44.1ms 76.7ms 91.8ms 0.16 s 256 1.07ms 1.06ms 4.79ms 7.81ms 7.97ms 12.2ms 19.2ms 53.6ms 84.6ms 96.4ms 0.17 s 0.20 s 0.36 s 384 2.84ms 2.69ms 15.3ms 25.9ms 25.6ms 33.5ms 61.5ms 73.3ms 0.13 s 0.13 s 0.23 s 0.28 s 0.50 s 512 5.95ms 5.53ms 36.5ms 61.0ms 62.7ms 83.0ms 0.15 s 0.15 s 0.25 s 0.25 s 0.42 s 0.53 s 0.96 s 768 16.3ms 15.5ms 0.12 s 0.21 s 0.23 s 0.25 s 0.38 s 0.38 s 0.64 s 0.69 s 1.11 s 1.40 s 2.37 s 1024 39.4ms 33.6ms 0.17 s 0.48 s 0.49 s 0.49 s 0.75 s 0.78 s 1.29 s 1.37 s 2.25 s 2.72 s 5.00 s 2048 0.28 s 0.28 s 1.00 s 1.49 s 2.00 s 2.36 s 3.49 s 4.36 s 6.18 s 7.52 s 11.8 s 15.2 s 27.6 s

Table 11. Matrix multiplication in using the fmpz_mat_mul function from Flint.

\ 8 16 24 32 52 64 104 128 208 256 416 512 832 8 0.25μs 0.38μs 0.38μs 0.38μs 0.38μs 0.38μs 2.13μs 2.17μs 6.39μs 6.01μs 13.7μs 19.2μs 36.8μs 16 0.55μs 2.13μs 2.13μs 2.11μs 2.12μs 2.14μs 11.5μs 11.5μs 39.7μs 41.9μs 90.1μs 0.12ms 0.25ms 32 1.65μs 7.69μs 7.77μs 7.70μs 10.6μs 10.2μs 71.0μs 70.2μs 0.27ms 0.37ms 0.65ms 0.91ms 1.90ms 64 6.31μs 38.9μs 39.0μs 39.0μs 62.3μs 60.4μs 0.50ms 0.47ms 2.10ms 2.51ms 5.02ms 7.06ms 8.49ms 96 15.6μs 0.12ms 0.12ms 0.12ms 0.25ms 0.25ms 1.56ms 1.50ms 5.75ms 6.85ms 11.8ms 13.1ms 22.8ms 128 38.1μs 97.4μs 97.4μs 0.32ms 0.61ms 0.61ms 3.61ms 3.45ms 13.2ms 14.7ms 24.9ms 28.6ms 49.4ms 192 0.10ms 0.28ms 0.27ms 0.96ms 1.97ms 1.96ms 11.9ms 20.6ms 35.8ms 41.2ms 71.1ms 82.5ms 0.14 s 256 0.23ms 0.66ms 0.63ms 2.37ms 4.38ms 4.37ms 36.2ms 46.4ms 79.9ms 92.0ms 0.15 s 0.18 s 0.32 s 384 0.74ms 1.90ms 1.78ms 7.42ms 14.0ms 14.0ms 0.11 s 0.14 s 0.24 s 0.27 s 0.46 s 0.55 s 0.93 s 512 1.72ms 3.93ms 17.6ms 23.5ms 33.6ms 41.5ms 0.17 s 0.21 s 0.37 s 0.42 s 0.70 s 0.82 s 1.36 s 768 5.45ms 12.3ms 48.9ms 65.2ms 90.7ms 0.11 s 0.48 s 0.60 s 0.98 s 1.11 s 1.84 s 2.15 s 3.56 s 1024 12.8ms 29.6ms 0.10 s 0.16 s 0.20 s 0.30 s 1.03 s 1.30 s 2.10 s 2.40 s 4.02 s 4.76 s 7.79 s 2048 91.4ms 0.24 s 0.75 s 0.99 s 1.58 s 1.88 s 6.41 s 8.03 s 12.9 s 14.8 s 24.9 s 25.4 s 43.4 s

Table 12. Matrix multiplication in using the nmod_mat_mul_ui (for ), nmod_mat_mul (for ), or mpn_mod_mat_mul (otherwise) function from Flint.

\ 8 16 24 32 52 64 104 128 208 256 416 512 832 8 5.50μs 5.45μs 5.49μs 6.48μs 5.71μs 7.01μs 8.01μs 9.23μs 13.4μs 14.6μs 23.7μs 26.9μs 54.9μs 16 43.8μs 43.4μs 43.1μs 59.6μs 45.3μs 60.7μs 65.3μs 80.8μs 0.10ms 0.12ms 0.20ms 0.23ms 0.45ms 32 0.33ms 0.33ms 0.33ms 0.49ms 0.35ms 0.50ms 0.51ms 0.64ms 0.84ms 0.98ms 1.63ms 1.83ms 3.55ms 64 2.58ms 2.57ms 2.56ms 3.91ms 2.68ms 4.02ms 4.21ms 5.44ms 8.19ms 9.83ms 15.6ms 17.2ms 30.4ms 96 9.21ms 9.10ms 9.10ms 14.1ms 9.48ms 13.5ms 15.4ms 18.0ms 36.5ms 42.2ms 68.3ms 81.2ms 0.13 s 128 32.1ms 25.2ms 31.2ms 46.3ms 35.1ms 44.9ms 53.0ms 59.5ms 0.11 s 0.12 s 0.20 s 0.20 s 0.34 s 192 0.13 s 0.13 s 0.13 s 0.19 s 0.13 s 0.19 s 0.20 s 0.23 s 0.32 s 0.33 s 0.55 s 0.66 s 1.10 s 256 0.32 s 0.32 s 0.33 s 0.46 s 0.34 s 0.46 s 0.48 s 0.59 s 0.74 s 0.80 s 1.31 s 1.41 s 2.65 s 384 1.16 s 1.15 s 1.18 s 1.60 s 1.23 s 1.62 s 1.77 s 2.02 s 2.41 s 2.58 s 4.48 s 5.31 s 8.50 s

Table 13. Multiplication of matrices in represented as NTL objects of type mat_ZZ.

\ 8 16 24 32 52 64 104 128 208 256 416 512 832 8 7.47μs 7.04μs 7.44μs 8.86μs 7.25μs 9.78μs 9.95μs 12.7μs 17.5μs 19.6μs 30.4μs 37.2μs 72.8μs 16 56.5μs 49.0μs 48.8μs 60.4μs 51.6μs 65.7μs 74.0μs 92.6μs 0.12ms 0.14ms 0.21ms 0.26ms 0.51ms 32 93.6μs 0.12ms 0.14ms 0.14ms 0.18ms 0.21ms 0.29ms 0.34ms 0.51ms 0.60ms 1.02ms 1.34ms 2.40ms 64 0.36ms 0.45ms 0.49ms 0.48ms 0.64ms 0.75ms 1.06ms 1.28ms 1.94ms 2.28ms 3.96ms 5.11ms 9.31ms 96 0.80ms 0.98ms 1.08ms 1.09ms 1.42ms 1.66ms 2.40ms 2.87ms 4.43ms 5.23ms 8.92ms 11.2ms 20.6ms 128 1.39ms 1.75ms 1.92ms 1.93ms 2.57ms 2.97ms 4.30ms 5.11ms 7.91ms 9.13ms 15.5ms 19.7ms 36.6ms 192 3.20ms 4.12ms 4.68ms 4.64ms 6.06ms 7.01ms 10.2ms 12.1ms 18.2ms 21.3ms 37.5ms 47.5ms 87.6ms 256 5.91ms 7.65ms 8.59ms 8.58ms 11.4ms 13.3ms 19.3ms 22.9ms 34.8ms 40.8ms 71.4ms 93.1ms 0.18 s 384 14.1ms 18.5ms 21.4ms 21.4ms 28.8ms 33.6ms 50.0ms 59.8ms 92.6ms 0.11 s 0.20 s 0.28 s 0.50 s 512 27.7ms 36.3ms 41.7ms 42.0ms 57.3ms 67.8ms 0.11 s 0.13 s 0.21 s 0.25 s 0.41 s 0.53 s 0.96 s 768 75.6ms 0.10 s 0.12 s 0.13 s 0.18 s 0.21 s 0.32 s 0.37 s 0.55 s 0.65 s 1.09 s 1.37 s 2.45 s 1024 0.17 s 0.23 s 0.28 s 0.28 s 0.38 s 0.45 s 0.67 s 0.79 s 1.20 s 1.42 s 2.34 s 2.93 s 5.12 s 2048 0.95 s 1.33 s 1.67 s 1.66 s 2.40 s 2.79 s 4.27 s 5.05 s 7.58 s 9.12 s 14.9 s 18.6 s 30.7 s

Table 14. Multiplication of matrices in represented as NTL objects of type mat_ZZ_p.

\ 8 16 24 32 52 8 0.58μs 0.57μs 0.58μs 0.74μs 0.73μs 16 0.84μs 0.85μs 0.85μs 3.97μs 4.03μs 32 3.09μs 3.04μs 2.84μs 12.4μs 16.3μs 64 14.5μs 14.4μs 14.3μs 92.0μs 0.12ms 96 38.8μs 38.4μs 38.1μs 0.30ms 0.41ms 128 78.1μs 78.2μs 77.3μs 0.71ms 0.95ms 192 0.22ms 0.22ms 0.26ms 2.38ms 3.19ms 256 0.49ms 0.49ms 0.55ms 5.62ms 7.61ms 384 1.53ms 1.53ms 1.84ms 19.1ms 25.6ms 512 3.54ms 3.54ms 4.35ms 41.1ms 54.8ms 768 11.6ms 11.4ms 14.3ms 0.14 s 0.18 s 1024 32.0ms 30.2ms 35.6ms 0.29 s 0.39 s 2048 0.27 s 0.26 s 0.30 s 2.10 s 2.81 s

Table 15. Multiplication of matrices in represented as NTL objects of type mat_zz_p.

\ 8 16 24 32 52 64 104 128 208 256 416 512 832 8 1.59ms 1.60ms 1.61ms 1.62ms 1.63ms 1.65ms 1.68ms 1.69ms 1.75ms 1.80ms 1.98ms 2.04ms 2.44ms 16 1.63ms 1.65ms 1.67ms 1.68ms 1.70ms 1.72ms 1.78ms 1.81ms 1.96ms 2.04ms 2.37ms 2.51ms 3.14ms 32 1.82ms 1.82ms 1.86ms 1.92ms 2.03ms 2.07ms 2.25ms 2.38ms 2.77ms 3.00ms 3.82ms 4.24ms 6.16ms 64 2.56ms 2.56ms 2.80ms 2.93ms 3.35ms 3.51ms 4.13ms 4.71ms 6.20ms 7.06ms 10.3ms 12.1ms 18.1ms 96 3.78ms 3.72ms 4.31ms 4.59ms 5.62ms 5.95ms 7.77ms 8.64ms 11.8ms 14.7ms 21.6ms 25.9ms 39.2ms 128 5.52ms 5.44ms 6.03ms 6.91ms 8.37ms 9.27ms 12.9ms 14.5ms 21.0ms 25.5ms 39.9ms 44.6ms 72.4ms 192 10.4ms 10.4ms 11.9ms 13.7ms 17.1ms 19.0ms 26.0ms 28.9ms 42.7ms 51.5ms 82.8ms 0.11 s 0.18 s 256 16.5ms 16.5ms 19.3ms 22.9ms 28.8ms 33.0ms 47.0ms 54.1ms 85.2ms 0.10 s 0.18 s 0.23 s 0.38 s 384 36.9ms 44.1ms 45.0ms 55.4ms 71.0ms 78.6ms 0.13 s 0.14 s 0.22 s 0.29 s 0.44 s 0.55 s 0.89 s 512 66.5ms 79.3ms 79.9ms 0.10 s 0.14 s 0.17 s 0.25 s 0.27 s 0.42 s 0.53 s 0.86 s 1.04 s 1.69 s 768 0.16 s 0.20 s 0.24 s 0.24 s 0.34 s 0.41 s 0.60 s 0.71 s 1.06 s 1.28 s 2.15 s 2.60 s 4.21 s 1024 0.32 s 0.40 s 0.49 s 0.50 s 0.70 s 0.87 s 1.25 s 1.43 s 2.17 s 2.64 s 4.22 s 5.08 s 8.36 s 2048 1.59 s 2.04 s 2.56 s 3.06 s 4.15 s 4.64 s 6.71 s 8.22 s 12.8 s 15.3 s 24.2 s 29.6 s 48.4 s

Table 16. Multiplication of matrices in with entries of type Givaro::ZRing<Givaro::Integer> using the fgemm function from Fflas.

A.2Multiplying small matrices using IFMA

Tables 17 and 18 show absolute counterparts of Tables 2, 3, and 4 from Section 4. Note that the left-hand of Table 17 corresponds to the left-hand of Table 2 (for multiplicands in ), whereas the right-hand corresponds to the right-hand of Table 3 (for multiplicands that use the alternative representation in ). The timings were obtained using mat_mul_bench –delay=5.

52 104 156 208 312 416 624 832 52 104 156 208 312 416 624 832 1 2 3 4 6 8 12 16 1 2 3 4 6 8 12 16 1 1.86ns 2.76ns 5.37ns 8.48ns 16.7ns 26.8ns 57.9ns 0.14μs 2 4.07ns 15.6ns 32.7ns 60.2ns 0.14μs 0.22μs 0.51μs 1.50μs 2.15ns 4.61ns 8.91ns 13.2ns 34.3ns 46.4ns 84.2ns 0.15μs 3 11.9ns 61.5ns 0.13μs 0.22μs 0.46μs 1.03μs 2.49μs 5.53μs 4 32.8ns 0.15μs 0.32μs 0.53μs 1.43μs 2.30μs 6.32μs 13.5μs 5.95ns 19.2ns 42.0ns 67.4ns 0.13μs 0.22μs 0.43μs 1.20μs 6 91.1ns 0.53μs 1.50μs 2.51μs 5.19μs 8.19μs 20.6μs 88.9μs 20.7ns 68.7ns 0.17μs 0.27μs 0.90μs 1.49μs 3.60μs 5.86μs 8 0.25μs 1.81μs 3.94μs 6.38μs 11.8μs 18.7μs 0.10ms 0.36ms 35.6ns 0.13μs 0.29μs 0.47μs 1.71μs 2.90μs 6.51μs 11.1μs 12 0.93μs 5.99μs 12.6μs 20.1μs 75.7μs 0.18ms 0.56ms 1.24ms 99.9ns 0.39μs 1.31μs 2.22μs 4.81μs 8.20μs 17.7μs 35.2μs 16 2.77μs 15.1μs 40.7μs 0.11ms 0.31ms 0.51ms 1.36ms 3.27ms 0.22μs 1.39μs 3.07μs 5.09μs 10.8μs 18.2μs 84.6μs 0.23ms 24 0.86μs 4.17μs 9.09μs 16.7μs 53.3μs 0.17ms 0.45ms 0.78ms 32 2.99μs 13.1μs 37.7μs 0.12ms 0.37ms 0.65ms 1.45ms 2.54ms multiplicands in multiplicands in

Table 17. Absolute timings for pure and tinkered SIMD multiplication of two matrices in or .

52 104 156 208 312 416 624 832 1 2 3 4 6 8 12 16 2 2.34ns 6.57ns 15.1ns 25.6ns 42.8ns 63.1ns 4 3.47ns 10.8ns 21.5ns 33.3ns 75.4ns 0.12μs 0.24μs 0.40μs 6 33.4ns 0.10μs 0.25μs 0.39μs 1.22μs 1.87μs 8 20.4ns 69.0ns 0.15μs 0.24μs 0.50μs 0.93μs 2.08μs 3.40μs 12 54.3ns 0.18μs 0.45μs 1.03μs 2.19μs 3.50μs 7.52μs 12.4μs 16 0.11μs 0.42μs 1.27μs 2.18μs 4.73μs 7.89μs 17.1μs 39.0μs 24 0.35μs 1.71μs 4.53μs 7.71μs 16.4μs 54.0μs 0.17ms 0.35ms 32 0.99μs 5.39μs 13.8μs 18.6μs 72.7μs 0.17ms 0.47ms 0.80ms

Table 18. Absolute timings for high multiplication of matrices in in tinkered SIMD mode.

A.3Multiplying large matrices using IFMA

Tables 19 presents the absolute counterparts of the timings in Table 6 from Section 5. Similarly, Table 20 provides absolute counterparts for the timings in Table 5.

52 104 156 208 312 416 624 832 51 103 155 207 311 415 623 831 1 2 3 4 6 8 12 16 1 2 3 4 6 8 12 16 32 3.87μs 11.8μs 20.1μs 29.4μs 61.2μs 98.7μs 0.22ms 0.48ms 5.41μs 11.2μs 18.8μs 25.5μs 44.7μs 69.6μs 0.15ms 0.27ms 64 19.4μs 61.5μs 0.11ms 0.19ms 0.42ms 0.73ms 1.58ms 3.89ms 28.5μs 59.0μs 97.0μs 0.15ms 0.26ms 0.39ms 0.77ms 1.26ms 96 51.9μs 0.19ms 0.38ms 0.63ms 1.32ms 2.25ms 5.03ms 11.8ms 82.1μs 0.19ms 0.31ms 0.44ms 0.77ms 1.14ms 2.05ms 3.22ms 128 0.11ms 0.42ms 0.84ms 1.41ms 3.04ms 5.25ms 15.3ms 27.6ms 0.19ms 0.43ms 0.70ms 0.97ms 1.60ms 2.32ms 4.08ms 6.27ms 192 0.36ms 1.30ms 2.70ms 4.57ms 10.1ms 18.7ms 51.0ms 86.5ms 0.63ms 1.33ms 2.10ms 2.83ms 4.51ms 6.42ms 11.5ms 28.7ms 256 0.82ms 2.98ms 6.21ms 10.7ms 24.8ms 43.4ms 0.11 s 0.19 s 1.44ms 2.94ms 4.53ms 6.14ms 10.1ms 24.3ms 42.5ms 63.8ms 384 2.89ms 10.3ms 21.5ms 36.4ms 85.7ms 0.15 s 0.37 s 0.66 s 5.01ms 10.4ms 16.4ms 22.3ms 35.8ms 64.2ms 0.13 s 0.24 s 512 6.49ms 23.3ms 49.3ms 91.2ms 0.20 s 0.35 s 0.89 s 1.55 s 11.4ms 23.2ms 35.7ms 48.3ms 84.5ms 0.19 s 0.32 s 0.49 s 768 21.8ms 78.1ms 0.17 s 0.31 s 0.68 s 1.17 s 3.01 s 5.22 s 39.4ms 77.7ms 0.12 s 0.16 s 0.29 s 0.64 s 1.09 s 1.64 s 1024 51.8ms 0.19 s 0.39 s 0.72 s 1.94 s 3.07 s 7.12 s 12.3 s 90.3ms 0.18 s 0.28 s 0.38 s 0.70 s 1.51 s 2.58 s 3.91 s 1536 0.17 s 0.62 s 1.30 s 2.42 s 5.28 s 9.25 s 23.9 s 41.5 s 0.30 s 0.61 s 0.93 s 1.29 s 2.14 s 4.93 s 8.48 s 12.7 s 2048 0.41 s 1.47 s 3.08 s 5.74 s 12.5 s 21.9 s 56.6 s 98.2 s 0.71 s 1.43 s 2.22 s 3.04 s 5.15 s 11.6 s 19.9 s 29.9 s Naive CRT

Table 19. Absolute timings for a long limb matrix product using IFMA, both for the naive algorithm and the CRT method. The top row indicates the corresponding bit-precisions .

51 102 153 204 306 408 612 816 51 102 153 204 306 408 612 816 1 2 3 4 6 8 12 16 1 2 3 4 6 8 12 16 32 3.19μs 7.58μs 15.3μs 23.6μs 50.8μs 91.5μs 0.22ms 0.45ms 2.32μs 5.60μs 9.48μs 15.6μs 31.0μs 52.0μs 0.12ms 0.22ms 64 15.7μs 40.1μs 85.1μs 0.15ms 0.33ms 0.58ms 1.30ms 6.49ms 10.1μs 28.8μs 49.9μs 92.0μs 0.19ms 0.34ms 0.74ms 1.25ms 96 43.9μs 0.13ms 0.28ms 0.49ms 1.03ms 1.75ms 9.01ms 17.7ms 26.4μs 87.8μs 0.17ms 0.31ms 0.64ms 1.06ms 2.18ms 8.89ms 128 96.2μs 0.31ms 0.64ms 1.08ms 2.20ms 8.01ms 19.5ms 40.7ms 55.9μs 0.22ms 0.39ms 0.68ms 1.37ms 2.27ms 10.2ms 19.2ms 192 0.33ms 1.00ms 1.98ms 3.21ms 12.3ms 23.7ms 60.7ms 0.11 s 0.18ms 0.70ms 1.22ms 2.05ms 4.02ms 13.0ms 31.7ms 59.3ms 256 0.76ms 2.24ms 4.38ms 12.0ms 28.3ms 54.2ms 0.12 s 0.22 s 0.43ms 1.54ms 2.65ms 4.44ms 16.1ms 29.8ms 69.7ms 0.12 s 384 2.64ms 7.73ms 15.5ms 32.8ms 0.11 s 0.20 s 0.48 s 0.89 s 1.58ms 5.63ms 9.61ms 16.3ms 59.9ms 0.11 s 0.26 s 0.47 s 512 6.01ms 17.3ms 33.9ms 99.7ms 0.22 s 0.42 s 0.98 s 1.77 s 3.49ms 12.0ms 20.8ms 35.3ms 0.13 s 0.24 s 0.55 s 0.96 s 768 20.3ms 57.9ms 0.12 s 0.34 s 0.75 s 1.42 s 3.29 s 5.96 s 11.8ms 40.2ms 70.0ms 0.12 s 0.45 s 0.79 s 1.86 s 3.22 s 1024 47.6ms 0.14 s 0.27 s 0.79 s 1.76 s 3.36 s 7.74 s 14.1 s 27.7ms 94.9ms 0.17 s 0.28 s 1.06 s 1.87 s 4.37 s 7.64 s 1536 0.16 s 0.46 s 0.90 s 2.68 s 5.93 s 11.2 s 26.1 s 47.6 s 93.2ms 0.33 s 0.55 s 0.94 s 3.58 s 6.28 s 14.7 s 25.6 s 2048 0.38 s 1.09 s 2.13 s 6.33 s 14.0 s 26.6 s 61.6 s 112. s 0.22 s 0.78 s 1.33 s 2.22 s 8.43 s 14.8 s 34.7 s 60.3 s long short

Table 20. Absolute timings for long and short limb matrix products using IFMA and the non-recursive variant of Karatsuba's algorithm. The top row also indicates the corresponding bit-precisions .

A.4Multiplying matrices using AMX

Table 21 shows absolute timings for naive AMX-accelerated matrix multiplication at various precisions; absolute timings for the CRT strategy are presented in Table 22. The corresponding relative timings with respect to the theoretical peak performance were given in Tables 7 and 8, respectively.

AMX N16 N24 N32 N48 N64 L16 L24 L32 L48 L64 I52 64 0.17μs 2.48μs 3.76μs 4.78μs 8.32μs 12.3μs 3.20μs 4.78μs 7.14μs 13.5μs 21.5μs 19.4μs 96 0.93μs 5.90μs 11.2μs 16.9μs 31.5μs 52.3μs 7.61μs 15.3μs 25.4μs 54.1μs 91.9μs 51.9μs 128 1.96μs 10.4μs 19.0μs 29.7μs 58.0μs 95.5μs 13.1μs 26.9μs 44.8μs 99.1μs 0.18ms 0.11ms 192 4.40μs 23.6μs 43.5μs 70.5μs 0.15ms 0.27ms 29.6μs 64.5μs 0.12ms 0.29ms 0.49ms 0.36ms 256 9.57μs 47.6μs 96.3μs 0.16ms 0.35ms 0.58ms 65.6μs 0.15ms 0.28ms 0.61ms 1.01ms 0.82ms 384 24.1μs 0.16ms 0.31ms 0.49ms 0.95ms 1.57ms 0.23ms 0.47ms 0.78ms 1.63ms 2.75ms 2.89ms 512 62.2μs 0.37ms 0.67ms 1.05ms 2.02ms 3.32ms 0.50ms 0.99ms 1.65ms 3.41ms 5.79ms 6.49ms 768 0.21ms 0.96ms 1.75ms 2.79ms 5.56ms 9.37ms 1.26ms 2.59ms 4.36ms 9.67ms 28.8ms 21.8ms 1024 0.45ms 1.96ms 3.61ms 5.94ms 21.9ms 34.6ms 2.73ms 5.45ms 19.8ms 39.3ms 62.4ms 51.8ms 1536 1.81ms 7.00ms 25.0ms 38.8ms 73.2ms 0.12 s 9.99ms 39.5ms 63.0ms 0.13 s 0.21 s 0.17 s 2048 4.24ms 30.5ms 55.8ms 84.9ms 0.16 s 0.28 s 43.9ms 85.7ms 0.14 s 0.28 s 0.46 s 0.41 s 3072 14.4ms 91.9ms 0.17 s 0.26 s 0.50 s 0.79 s 0.13 s 0.26 s 0.42 s 0.84 s 1.37 s 1.41 s 4096 48.0ms 0.22 s 0.41 s 0.67 s 1.29 s 2.11 s 0.31 s 0.64 s 1.06 s 2.18 s 3.67 s 3.36 s

Table 21. Time to multiply two matrices with bit positive integer coefficients using naive AMX-accelerated arithmetic. The columns Np corresponds to a bit low product, whereas the columns Lp corresponds to a bit long product. The column AMX corresponds to the raw bit product from Section 6.2. For reference, the column I52 indicates the timings of limb products using IFMA.