Preparing an antisymmetric matrix of normal numbers as efficiently as possible

General Usage Performance
5

Here is a little recursive cache-oblivious implementation of this idea. (Recursion can be fast!)

It first generates all of the N (N-1)/2 random numbers with a single SIMD-vectorized randn! call in a contiguous temporary buffer v, and then copies this data anti-symmetrically into the matrix with a cache-oblivious algorithm that divides the array recursively into smaller and smaller blocks until ending up at tiny blocks (that should fit into L1 cache but big enough to amortize the recursion overhead) where it uses a simple loop. (I experimented a little and 16 seems like a good block size on my machine.)

julia> N = 1024; A = zeros(N, N); v = zeros((N * (N-1)) >> 1);

julia> @btime prepare!($A); # your original version, fixed to take A as an argument
  1.633 ms (0 allocations: 0 bytes)

julia> @btime prepare3!($A, $v);
  817.416 μs (0 allocations: 0 bytes)

julia> @btime randn!(v); # lower bound: just the raw random-number generation cost
  553.708 μs (0 allocations: 0 bytes)

I didn’t put many comments in my code, so I’ll leave deciphering it as an exercise :wink:. (Caveat: it’s only lightly tested.)

Code 
using LinearAlgebra, Random, BenchmarkTools

function _prepare3_diag!(A, i0, n, v, k)
    if n <= 16
        for i = i0:i0+(n-1)
            @inbounds A[i,i] = 0
            for j = i0:i-1
                @inbounds A[i,j] = -(A[j,i] = v[k += 1])
            end
        end
    else
        n2 = n >> 1
        k += _prepare3_diag!(A, i0, n2, v, k)
        k += _prepare3_diag!(A, i0+n2, n-n2, v, k)
        _prepare3_offdiag!(A, i0, i0+n2, n2, n-n2, v, k)
    end
    return (n * (n-1)) >> 1 # number of values consumed
end

function _prepare3_offdiag!(A, i0, j0, nr, nc, v, k)
    if nr * nc <= 256
        for i = i0:i0+(nr-1), j = j0:j0+(nc-1)
            @inbounds A[i,j] = -(A[j,i] = v[k += 1])
        end
    else
        nr2 = nr >> 1
        nc2 = nc >> 1
        k += _prepare3_offdiag!(A, i0, j0, nr2, nc2, v, k)
        k += _prepare3_offdiag!(A, i0+nr2, j0, nr-nr2, nc2, v, k)
        k += _prepare3_offdiag!(A, i0, j0+nc2, nr2, nc-nc2, v, k)
        _prepare3_offdiag!(A, i0+nr2, j0+nc2, nr-nr2, nc-nc2, v, k)
    end
    return nr * nc
end

function prepare3!(A::AbstractMatrix, v::AbstractVector)
    N = LinearAlgebra.checksquare(A)
    length(v) == ((N * (N-1)) >> 1) || throw(DimensionMismatchError())
    randn!(v) # SIMD-vectorized random-number generation
    _prepare3_diag!(A, 1, N, v, 0)
    return A
end

Note, by the way, that a power-of-two size like N=1024 above is usually the worst case for cache-associativity conflicts. N=1023 is substantially faster for both the original code and the cache-oblivious version:

julia> N = 1023; A = zeros(N, N); v = zeros((N * (N-1)) >> 1);

julia> @btime prepare!($A);  # original version, fixed to take A as an argument
  1.387 ms (0 allocations: 0 bytes)

julia> @btime prepare3!($A, $v);
  758.750 μs (0 allocations: 0 bytes)

(This is all on an Apple M4 from 2024.)

Of course, this doesn’t get into threading, but in principle you could spawn threads from the top-level recursive subdivisions to do them in parallel while still maintaining locality with each thread. Not sure how big N has to be for that to be beneficial, though.