  We develop an algorithm for i) computing generalized regular $k$-point grids,
ii) reducing the grids to their symmetrically distinct points, and iii) mapping
the reduced grid points into the Brillouin zone. The algorithm exploits the
connection between integer matrices and finite groups to achieve a
computational complexity that is linear with the number of $k$-points. The
favorable scaling means that, at a given $k$-point density, all possible
commensurate grids can be generated (as suggested by Moreno and Soler) and
quickly reduced to identify the grid with the fewest symmetrically unique
$k$-points. These optimal grids provide significant speed-up compared to
Monkhorst-Pack $k$-point grids; they have better symmetry reduction resulting
in fewer irreducible $k$-points at a given grid density. The integer nature of
this new reduction algorithm also simplifies issues with finite precision in
current implementations. The algorithm is available as open source software.
