A C implementation of the Smith massager algorithm

TL;DR

This paper presents a highly efficient C implementation of a Smith normal form algorithm for integer matrices, achieving near-matrix multiplication speed through advanced optimization techniques.

Contribution

The authors introduce practical implementation strategies that significantly improve the performance of the Smith massager algorithm, bridging theory and real-world application.

Findings

Implementation scales proportionally to matrix multiplication time.

Optimization techniques enable collapsing multiple iterations into one.

Experimental results up to 10,007 dimensions demonstrate practical efficiency.

Abstract

We describe a C implementation of the Las Vegas algorithm of Birmpilis, Labahn and Storjohann from 2020 for computing the Smith normal form of a nonsingular integer matrix. The algorithm computes a Smith massager for the input matrix using $O(n^{\omega}\, \B(\log n + \log \|A\|)\, (\log n)^2)$ bit operations, which is softly equivalent to the cost of multiplying two matrices of the same dimension and entry size. We describe the key implementation techniques that bridge the gap between the theoretical algorithm and practical performance, including BLAS-accelerated modular arithmetic via the Residue Number System and an adaptive batching scheme that collapses the theoretical $O(\log n)$ iterations to $O(1)$ in practice. Experiments on matrices of dimension up to $n = 10007$ show that the implementation's running time scales proportionally to that of a single BLAS matrix multiplication, with both exhibiting the same effective growth rate on a log-log plot.