diff --git a/docs/src/getting-started.md b/docs/src/getting-started.md
index 263e964..547d640 100644
--- a/docs/src/getting-started.md
+++ b/docs/src/getting-started.md
@@ -6,6 +6,9 @@ CurrentModule = Octavian
 
 ## Multi-threaded matrix multiplication: `matmul!` and `matmul`
 
+Octavian exports the functions `matmul!` and `matmul`, which provide
+multithreaded matrix multiplication in pure Julia.
+
 Remember to start Julia with multiple threads with e.g. one of the following:
 - `julia -t auto`
 - `julia -t 4`
@@ -20,7 +23,7 @@ B = [7 8 9 10; 11 12 13 14; 15 16 17 18]
 
 C = Matrix{Int}(undef, 2, 4)
 
-Octavian.matmul!(C, A, B) # (multi-threaded) multiply A×B and store the result in C (overwriting the contents of C)
+matmul!(C, A, B) # (multi-threaded) multiply A×B and store the result in C (overwriting the contents of C)
 
 C
 
@@ -34,7 +37,7 @@ A = [1 2 3; 4 5 6]
 
 B = [7 8 9 10; 11 12 13 14; 15 16 17 18]
 
-C = Octavian.matmul(A, B) # (multi-threaded) multiply A×B and return the result
+C = matmul(A, B) # (multi-threaded) multiply A×B and return the result
 
 C
 
diff --git a/test/_matmul.jl b/test/_matmul.jl
new file mode 100644
index 0000000..027254a
--- /dev/null
+++ b/test/_matmul.jl
@@ -0,0 +1,101 @@
+# The following variables need to be defined before `include`-ing this file:
+# `testset_name_suffix`
+# `n_values`
+# `k_values`
+# `m_values`
+
+@time @testset "Matrix Multiply Float32 $(testset_name_suffix)" begin
+    T = Float32
+    for n ∈ n_values
+        for k ∈ k_values
+            for m ∈ m_values
+                A = rand(T, m, k)
+                B = rand(T, k, n)
+                A′ = permutedims(A)'
+                B′ = permutedims(B)'
+                AB = A * B; A′B = A′ * B; AB′ = A * B′; A′B′ = A′ * B′;
+                @info "" T n k m
+                @test @time(Octavian.matmul(A, B)) ≈ AB
+                @test @time(Octavian.matmul(A′, B)) ≈ A′B
+                @test @time(Octavian.matmul(A, B′)) ≈ AB′
+                @test @time(Octavian.matmul(A′, B′)) ≈ A′B′
+                @test @time(Octavian.matmul_serial(A, B)) ≈ AB
+                @test @time(Octavian.matmul_serial(A′, B)) ≈ A′B
+                @test @time(Octavian.matmul_serial(A, B′)) ≈ AB′
+                @test @time(Octavian.matmul_serial(A′, B′)) ≈ A′B′
+            end
+        end
+    end
+end
+
+@time @testset "Matrix Multiply Float64 $(testset_name_suffix)" begin
+    T = Float64
+    for n ∈ n_values
+        for k ∈ k_values
+            for m ∈ m_values
+                A = rand(T, m, k)
+                B = rand(T, k, n)
+                A′ = permutedims(A)'
+                B′ = permutedims(B)'
+                AB = A * B; A′B = A′ * B; AB′ = A * B′; A′B′ = A′ * B′;
+                @info "" T n k m
+                @test @time(Octavian.matmul(A, B)) ≈ AB
+                @test @time(Octavian.matmul(A′, B)) ≈ A′B
+                @test @time(Octavian.matmul(A, B′)) ≈ AB′
+                @test @time(Octavian.matmul(A′, B′)) ≈ A′B′
+                @test @time(Octavian.matmul_serial(A, B)) ≈ AB
+                @test @time(Octavian.matmul_serial(A′, B)) ≈ A′B
+                @test @time(Octavian.matmul_serial(A, B′)) ≈ AB′
+                @test @time(Octavian.matmul_serial(A′, B′)) ≈ A′B′
+            end
+        end
+    end
+end
+
+@time @testset "Matrix Multiply Int32 $(testset_name_suffix)" begin
+    T = Int32
+    for n ∈ n_values
+        for k ∈ k_values
+            for m ∈ m_values
+                A = rand(T, m, k)
+                B = rand(T, k, n)
+                A′ = permutedims(A)'
+                B′ = permutedims(B)'
+                AB = A * B; A′B = A′ * B; AB′ = A * B′; A′B′ = A′ * B′;
+                @info "" T n k m
+                @test @time(Octavian.matmul(A, B)) == AB
+                @test @time(Octavian.matmul(A′, B)) == A′B
+                @test @time(Octavian.matmul(A, B′)) == AB′
+                @test @time(Octavian.matmul(A′, B′)) == A′B′
+                @test @time(Octavian.matmul_serial(A, B)) == AB
+                @test @time(Octavian.matmul_serial(A′, B)) == A′B
+                @test @time(Octavian.matmul_serial(A, B′)) == AB′
+                @test @time(Octavian.matmul_serial(A′, B′)) == A′B′
+            end
+        end
+    end
+end
+
+@time @testset "Matrix Multiply Int64 $(testset_name_suffix)" begin
+    T = Int64
+    for n ∈ n_values
+        for k ∈ k_values
+            for m ∈ m_values
+                A = rand(T, m, k)
+                B = rand(T, k, n)
+                A′ = permutedims(A)'
+                B′ = permutedims(B)'
+                AB = A * B; A′B = A′ * B; AB′ = A * B′; A′B′ = A′ * B′;
+                @info "" T n k m
+                @test @time(Octavian.matmul(A, B)) == AB
+                @test @time(Octavian.matmul(A′, B)) == A′B
+                @test @time(Octavian.matmul(A, B′)) == AB′
+                @test @time(Octavian.matmul(A′, B′)) == A′B′
+                @test @time(Octavian.matmul_serial(A, B)) == AB
+                @test @time(Octavian.matmul_serial(A′, B)) == A′B
+                @test @time(Octavian.matmul_serial(A, B′)) == AB′
+                @test @time(Octavian.matmul_serial(A′, B′)) == A′B′
+            end
+        end
+    end
+end
diff --git a/test/matmul.jl b/test/matmul.jl
index 90ed726..f6d2005 100644
--- a/test/matmul.jl
+++ b/test/matmul.jl
@@ -1,108 +1,7 @@
-@time @testset "Matrix Multiply Float32" begin
-    Mc, Kc, Nc = map(Int, Octavian.block_sizes(Float32))
-    for logn ∈ range(log(1), log(2Nc+1), length = 7)
-        n = round(Int, exp(logn))
-        for logk ∈ range(log(1), log(2Kc+1), length = 7)
-            k = round(Int, exp(logk))
-            B = rand(Float32, k, n)
-            B′ = permutedims(B)'
-            for logm ∈ range(log(1), log(2Mc+1), length = 7)
-                m = round(Int, exp(logm))
-                A = rand(Float32, m, k)
-                A′ = permutedims(A)'
-                @show m, k, n
-                AB = A * B; A′B = A′ * B; AB′ = A * B′; A′B′ = A′ * B′;
-                @test @time(Octavian.matmul(A, B)) ≈ AB
-                @test @time(Octavian.matmul(A′, B)) ≈ A′B
-                @test @time(Octavian.matmul(A, B′)) ≈ AB′
-                @test @time(Octavian.matmul(A′, B′)) ≈ A′B′
-                @test @time(Octavian.matmul_serial(A, B)) ≈ AB
-                @test @time(Octavian.matmul_serial(A′, B)) ≈ A′B
-                @test @time(Octavian.matmul_serial(A, B′)) ≈ AB′
-                @test @time(Octavian.matmul_serial(A′, B′)) ≈ A′B′
-            end
-        end
-    end
-end
+n_values  = [200, 300, 400]
+k_values  = [200, 300, 400]
+m_values  = [200, 300, 400]
 
-@time @testset "Matrix Multiply Float64" begin
-    Mc, Kc, Nc = map(Int, Octavian.block_sizes(Float64))
-    for logn ∈ range(log(1), log(2Nc+1), length = 7)
-        n = round(Int, exp(logn))
-        for logk ∈ range(log(1), log(2Kc+1), length = 7)
-            k = round(Int, exp(logk))
-            B = rand(Float64, k, n)
-            B′ = permutedims(B)'
-            for logm ∈ range(log(1), log(2Mc+1), length = 7)
-                m = round(Int, exp(logm))
-                A = rand(Float64, m, k)
-                A′ = permutedims(A)'
-                @show m, k, n
-                AB = A * B; A′B = A′ * B; AB′ = A * B′; A′B′ = A′ * B′;
-                @test @time(Octavian.matmul(A, B)) ≈ AB
-                @test @time(Octavian.matmul(A′, B)) ≈ A′B
-                @test @time(Octavian.matmul(A, B′)) ≈ AB′
-                @test @time(Octavian.matmul(A′, B′)) ≈ A′B′
-                @test @time(Octavian.matmul_serial(A, B)) ≈ AB
-                @test @time(Octavian.matmul_serial(A′, B)) ≈ A′B
-                @test @time(Octavian.matmul_serial(A, B′)) ≈ AB′
-                @test @time(Octavian.matmul_serial(A′, B′)) ≈ A′B′
-            end
-        end
-    end
-end
-
-@time @testset "Matrix Multiply Int32" begin
-    Mc, Kc, Nc = map(Int, Octavian.block_sizes(Int32))
-    for logn ∈ range(log(1), log(1.5Nc+1), length = 7)
-        n = round(Int, exp(logn))
-        for logk ∈ range(log(1), log(1.5Kc+1), length = 7)
-            k = round(Int, exp(logk))
-            B = rand(Int32, k, n)
-            B′ = permutedims(B)'
-            for logm ∈ range(log(1), log(1.5Mc+1), length = 7)
-                m = round(Int, exp(logm))
-                A = rand(Int32, m, k)
-                A′ = permutedims(A)'
-                @show m, k, n
-                AB = A * B; A′B = A′ * B; AB′ = A * B′; A′B′ = A′ * B′;
-                @test @time(Octavian.matmul(A, B)) ≈ AB
-                @test @time(Octavian.matmul(A′, B)) ≈ A′B
-                @test @time(Octavian.matmul(A, B′)) ≈ AB′
-                @test @time(Octavian.matmul(A′, B′)) ≈ A′B′
-                @test @time(Octavian.matmul_serial(A, B)) ≈ AB
-                @test @time(Octavian.matmul_serial(A′, B)) ≈ A′B
-                @test @time(Octavian.matmul_serial(A, B′)) ≈ AB′
-                @test @time(Octavian.matmul_serial(A′, B′)) ≈ A′B′
-            end
-        end
-    end
-end
-
-@time @testset "Matrix Multiply Int64" begin
-    Mc, Kc, Nc = map(Int, Octavian.block_sizes(Int64))
-    for logn ∈ range(log(1), log(1.5Nc+1), length = 7)
-        n = round(Int, exp(logn))
-        for logk ∈ range(log(1), log(1.5Kc+1), length = 7)
-            k = round(Int, exp(logk))
-            B = rand(Int64, k, n)
-            B′ = permutedims(B)'
-            for logm ∈ range(log(1), log(1.5Mc+1), length = 7)
-                m = round(Int, exp(logm))
-                A = rand(Int64, m, k)
-                A′ = permutedims(A)'
-                @show m, k, n
-                AB = A * B; A′B = A′ * B; AB′ = A * B′; A′B′ = A′ * B′;
-                @test @time(Octavian.matmul(A, B)) ≈ AB
-                @test @time(Octavian.matmul(A′, B)) ≈ A′B
-                @test @time(Octavian.matmul(A, B′)) ≈ AB′
-                @test @time(Octavian.matmul(A′, B′)) ≈ A′B′
-                @test @time(Octavian.matmul_serial(A, B)) ≈ AB
-                @test @time(Octavian.matmul_serial(A′, B)) ≈ A′B
-                @test @time(Octavian.matmul_serial(A, B′)) ≈ AB′
-                @test @time(Octavian.matmul_serial(A′, B′)) ≈ A′B′
-            end
-        end
-    end
-end
+testset_name_suffix = "(main)"
 
+include("_matmul.jl")
diff --git a/test/matmul_coverage.jl b/test/matmul_coverage.jl
index 2596fe9..4dd9482 100644
--- a/test/matmul_coverage.jl
+++ b/test/matmul_coverage.jl
@@ -1,75 +1,7 @@
-@time @testset "Matrix Multiply Float32 (coverage)" begin
-    m = 20
-    n = 20
-    k = 20
-    A = rand(Float32, m, k)
-    B = rand(Float32, k, n)
-    A′ = permutedims(A)'
-    B′ = permutedims(B)'
-    @show m, k, n
-    @test @time(Octavian.matmul(A, B)) ≈ A * B
-    @test @time(Octavian.matmul(A′, B)) ≈ A′ * B
-    @test @time(Octavian.matmul(A, B′)) ≈ A * B′
-    @test @time(Octavian.matmul(A′, B′)) ≈ A′ * B′
-end
+n_values  = [10, 20, 50, 100, 150, 200]
+k_values  = [10, 20, 50, 100, 150, 200]
+m_values  = [10, 20, 50, 100, 150, 200]
 
-@time @testset "Matrix Multiply Float64 (coverage)" begin
-    m = 20
-    n = 20
-    k = 20
-    A = rand(Float64, m, k)
-    B = rand(Float64, k, n)
-    A′ = permutedims(A)'
-    B′ = permutedims(B)'
-    @show m, k, n
-    @test @time(Octavian.matmul(A, B)) ≈ A * B
-    @test @time(Octavian.matmul(A′, B)) ≈ A′ * B
-    @test @time(Octavian.matmul(A, B′)) ≈ A * B′
-    @test @time(Octavian.matmul(A′, B′)) ≈ A′ * B′
-end
-
-@time @testset "Matrix Multiply Int32 (coverage)" begin
-    m = 20
-    n = 20
-    k = 20
-    A = rand(Int32, m, k)
-    B = rand(Int32, k, n)
-    A′ = permutedims(A)'
-    B′ = permutedims(B)'
-    @show m, k, n
-    @test @time(Octavian.matmul(A, B)) == A * B
-    @test @time(Octavian.matmul(A′, B)) == A′ * B
-    @test @time(Octavian.matmul(A, B′)) == A * B′
-    @test @time(Octavian.matmul(A′, B′)) == A′ * B′
-end
-
-@time @testset "Matrix Multiply Int64 (coverage)" begin
-    m = 20
-    n = 20
-    k = 20
-    A = rand(Int32, m, k)
-    B = rand(Int32, k, n)
-    A′ = permutedims(A)'
-    B′ = permutedims(B)'
-    @show m, k, n
-    @test @time(Octavian.matmul(A, B)) == A * B
-    @test @time(Octavian.matmul(A′, B)) == A′ * B
-    @test @time(Octavian.matmul(A, B′)) == A * B′
-    @test @time(Octavian.matmul(A′, B′)) == A′ * B′
-end
-
-@time @testset "A not-a-StrideArray" begin
-    m = 20
-    n = 20
-    k = 20
-    A = view(rand(Float64, 2m, 2k), 1:2:2m, 2:2:2k)
-    B = rand(Float64, k, n)
-    A′ = view(permutedims(parent(A))', 1:2:2m, 2:2:2k)
-    B′ = permutedims(B)'
-    @show m, k, n
-    @test @time(Octavian.matmul(A, B)) ≈ A * B
-    @test @time(Octavian.matmul(A′, B)) ≈ A′ * B
-    @test @time(Octavian.matmul(A, B′)) ≈ A * B′
-    @test @time(Octavian.matmul(A′, B′)) ≈ A′ * B′
-end
+testset_name_suffix = "(coverage)"
 
+include("_matmul.jl")
diff --git a/test/runtests.jl b/test/runtests.jl
index 384e49b..e69738c 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,5 +1,6 @@
 import Octavian
 
+import Aqua
 import BenchmarkTools
 import InteractiveUtils
 import LinearAlgebra
@@ -19,9 +20,6 @@ include("test_suite_preamble.jl")
 
 Random.seed!(123)
 
-using Aqua: test_all
-test_all(Octavian)
-
 include("utils.jl")
 include("block_sizes.jl")
 include("init.jl")
@@ -31,3 +29,7 @@ include("matmul_coverage.jl")
 if !coverage
     include("matmul.jl")
 end
+
+@testset "Aqua.jl" begin
+    Aqua.test_all(Octavian)
+end