Batches and SIMD

Batches

Batch routines work on many points or matrices at once, stored one after another in a flat float-array, a specialised array of single floats: points as x, y and z, matrices as 16 floats each, column by column.

Function	Does
`(mat4-transform-points matrix source destination &key start end)`	Transforms the points of `source` by an affine matrix, and stores them at the same places in `destination`.
`(mat4-transform-directions matrix source destination &key start end)`	The same for directions, without the translation.
`(mat4*-batch matrix matrices destination &key start end)`	Multiplies `matrix` by each matrix of `matrices`, as `mat4*` does, such as a parent's transform applied to its children's.

destination may be the source array, to work in place.
start and end count items, points or matrices, not floats: by default, all of them.
Each returns destination, and allocates nothing.

The batches of doubles, dmat4-transform-points, dmat4-transform-directions and dmat4*-batch, work the same way on a double-array, a specialised array of doubles.

(let ((points (make-array 9 :element-type 'single-float
                            :initial-contents '(0.0 0.0 0.0  1.0 0.0 0.0  0.0 1.0 0.0)))
      (world (make-array 9 :element-type 'single-float :initial-element 0.0)))
  (gm:mat4-transform-points (gm:mat4-translation (gm:make-vec3 :z -5)) points world))
;; => #(0.0 0.0 -5.0 1.0 0.0 -5.0 0.0 1.0 -5.0)

With a job system

A range of items is exactly what a chunk of cl-threads' parallel-for gives, so the workers of a job pool can share a batch:

(asdf:load-system :cl-threads)

(defpackage #:example
  (:use #:cl)
  (:local-nicknames (#:gm #:cl-gmath)
                    (#:th #:cl-threads)))
(in-package #:example)

(th:with-job-pool (pool)
  (let* ((count 100000)
         (points (make-array (* 3 count) :element-type 'single-float :initial-element 1.0))
         (world (make-array (* 3 count) :element-type 'single-float :initial-element 0.0))
         (model (gm:mat4-translation (gm:make-vec3 :z -5))))
    (th:parallel-for pool 0 count
                     (lambda (start end)
                       (gm:mat4-transform-points model points world :start start :end end)))
    (subseq world 0 6)))
;; => #(1.0 1.0 -4.0 1.0 1.0 -4.0)

The lambda is a new closure each time it is made; to run a batch every frame without allocating, make the chunk function once, as cl-threads' Parallel For shows.

SIMD

Where SBCL has sb-simd, on x86-64 and on ARM64 from SBCL 2.6.7, the 4x4 operations and the batches use SIMD instructions. Each has a kernel for each instruction set that measured faster than the one below it, and runs the best one the processor has. sb-simd's instruction-set-case makes the choice when cl-gmath is loaded and again whenever an image starts, so an image saved on a machine with AVX still runs on one without it.

Every kernel computes exactly what the scalar code computes: the same operations in the same order, sums added in pairs, and no fused multiply-add, which rounds once where a multiplication and an addition round twice. So every processor gets the same bits. That is what a game needs to compute the same on every player's machine, as lockstep multiplayer and replays do.

The kernels, and the time each takes on a Ryzen 9 9950X3D with SBCL 2.6.9, per call or per item of a batch:

Operation	Scalar	SSE, SSE2	AVX
`mat4*`	24 ns	3.7 ns	1.5 ns
`mat4*vec4`	3.6 ns	1.3 ns
`mat4-transpose`	2.3 ns	1.3 ns
`mat4*-batch`, a matrix	6.8 ns	4.6 ns	2.4 ns
`mat4-transform-points`, a point	1.1 ns	0.93 ns	0.44 ns
`mat4-transform-directions`, a direction	1.0 ns	0.88 ns	0.41 ns
`dmat4*`	24 ns	7.3 ns	2.7 ns
`dmat4*dvec4`	3.6 ns	2.2 ns	1.1 ns
`dmat4-transpose`	2.3 ns	1.2 ns	0.92 ns
`dmat4*-batch`, a matrix	6.9 ns		4.7 ns

A call through a function, such as mat4*, adds about 2 ns of checking its arguments and choosing its kernel; a batch adds it once. Batches of points work on four points at a time with SSE and eight with AVX, their x, y and z gathered into registers of their own. Everything else is scalar. SIMD values never appear in the API: every function takes and returns the same arrays either way.

Instruction sets left out

SSE4.1 adds blends, rounding, and inserting and extracting lanes, of which only the extraction could serve, for transforming points one at a time. That took 1.78 ns a point, slower than SSE's four at a time, and than the scalar code.
AVX2 adds, of what cl-gmath could use, a permute of doubles across a whole register, which made dmat4* 7% faster and dmat4*dvec4 a third slower. Its broadcast from a register is AVX2's too, which is why the AVX kernels spread values with AVX's own permutes.
AVX-512, which sb-simd has from SBCL 2.6.9, has no permutes there yet, which a 4x4 product needs to spread over 16 lanes, and AVX already holds a whole column of either precision.
FMA would change the results, and in the AVX kernels made dmat4* and mat4*-batch about 9% faster, dmat4*dvec4 and dmat4*-batch slower, and the rest no faster.

ARM64

On ARM64, NEON kernels do mat4*, mat4*vec4, mat4-transpose, mat4*-batch, dmat4*, dmat4*dvec4 and dmat4-transpose, the counterparts of the x86-64 kernels that were faster. They are tested under QEMU's emulation, which computes exactly but cannot be timed; see Getting Started.

Which path a build uses

Where sb-simd exists and the build was not told to stay scalar, loading cl-gmath pushes the feature :cl-gmath-simd:

(find :cl-gmath-simd *features*)   ; => :CL-GMATH-SIMD on SBCL on x86-64 and ARM64

The kernels of every instruction set are compiled in then, and the choice between them is made when an image starts. To build without SIMD, push :cl-gmath-scalar onto *features* before cl-gmath is compiled, and compile it again, since compiled files built the other way are not rebuilt on their own:

(pushnew :cl-gmath-scalar *features*)
(asdf:load-system :cl-gmath :force '("cl-gmath"))

The scalar code stays the reference: the tests compare every kernel the processor has with it, exactly, also when the result is one of the arguments.

cl-gmath

Reference

API Reference