# Preparing an antisymmetric matrix of normal numbers as efficiently as possible

**URL:** https://discourse.julialang.org/t/preparing-an-antisymmetric-matrix-of-normal-numbers-as-efficiently-as-possible/138954
**Category:** Performance
**Created:** [August 20, 2026, 5:40pm UTC](https://discourse.julialang.org/t/preparing-an-antisymmetric-matrix-of-normal-numbers-as-efficiently-as-possible/138954 "2026-08-20T17:40:08Z")
**Posts on this page:** 1
**Showing post:** 5

<div class="post-metadata">

### Author: ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)
#### Post date: [August 20, 2026, 7:44pm UTC](https://discourse.julialang.org/t/preparing-an-antisymmetric-matrix-of-normal-numbers-as-efficiently-as-possible/138954/5 "2026-08-20T19:44:36Z")

</div>

> [@stevengj](#):
>
> You really want to do neither, in principle — it seems similar to optimizing a transpose operation, in which the best approaches involve some form of blocking/tiling in order to maximize cache-line utilization. See e.g. [WIP: cache oblivious linear algebra algorithms - Pull Request #6690 - JuliaLang/julia - GitHub](https://github.com/JuliaLang/julia/pull/6690)

Here is a little [recursive cache-oblivious implementation](https://en.wikipedia.org/wiki/Cache-oblivious_algorithm) of this idea. ([Recursion can be fast!](https://discourse.julialang.org/t/recursion-in-julia-bad-idea/99304/5))

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-auto
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 😉. (Caveat: it’s only lightly tested.)

> **Code**
>
> ```julia-auto
> 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-auto
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.

---

_[View the full topic](https://discourse.julialang.org/t/preparing-an-antisymmetric-matrix-of-normal-numbers-as-efficiently-as-possible/138954)._