You signed in with another tab or window. Reload to refresh your session.You signed out in another tab or window. Reload to refresh your session.You switched accounts on another tab or window. Reload to refresh your session.Dismiss alert
The plan is for `LoopVectorization.jl` to be maintained through the [SciML Small Grants](https://sciml.ai/small_grants/#update_loopvectorization_to_support_changes_in_julia_v112_200) program.
18
18
If you would like to see an issue fixed, or support extended to another Julia patch release, please consider:
19
19
1. Donating to the SciML Small Grants program, with a note on the purpose of your donation.
@@ -58,7 +58,7 @@ For simple loops like a dot product, LoopVectorization.jl's most important optim
58
58
<details>
59
59
<summaryClick me! ></summary>
60
60
<p>
61
-
61
+
62
62
```julia
63
63
julia>using LoopVectorization, BenchmarkTools
64
64
@@ -314,7 +314,7 @@ f = KwargCall(round, (digits = 3,));
314
314
<p>
315
315
316
316
The key to the `@turbo` macro's performance gains is leveraging knowledge of exactly how data like `Float64`s and `Int`s are handled by a CPU. As such, it is not strightforward to generalize the `@turbo` macro to work on arrays containing structs such as `Matrix{Complex{Float64}}`. Instead, it is currently recommended that users wishing to apply `@turbo` to arrays of structs use packages such as [StructArrays.jl](https://github.com/JuliaArrays/StructArrays.jl) which transform an array where each element is a struct into a struct where each element is an array. Using StructArrays.jl, we can write a matrix multiply (gemm) kernel that works on matrices of `Complex{Float64}`s and `Complex{Int}`s:
317
-
```julia
317
+
```julia
318
318
using LoopVectorization, LinearAlgebra, StructArrays, BenchmarkTools, Test
319
319
320
320
BLAS.set_num_threads(1); @show BLAS.vendor()
@@ -373,7 +373,7 @@ julia> @test C1 β C2
373
373
Test Passed
374
374
```
375
375
376
-
Similar approaches can be taken to make kernels working with a variety of numeric struct types such as [dual numbers](https://github.com/JuliaDiff/DualNumbers.jl), [DoubleFloats](https://github.com/JuliaMath/DoubleFloats.jl), etc.
376
+
Similar approaches can be taken to make kernels working with a variety of numeric struct types such as [dual numbers](https://github.com/JuliaDiff/DualNumbers.jl), [DoubleFloats](https://github.com/JuliaMath/DoubleFloats.jl), etc.
Copy file name to clipboardExpand all lines: docs/src/examples/dot_product.md
+4-6Lines changed: 4 additions & 6 deletions
Display the source diff
Display the rich diff
Original file line number
Diff line number
Diff line change
@@ -11,7 +11,7 @@ function jdotavx(a, b)
11
11
s
12
12
end
13
13
```
14
-
To execute the loop using SIMD (Single Instruction Multiple Data) instructions, you have to unroll the loop. Rather than evaluating the loop as written -- adding element-wise products to a single accumulator one after the other -- you can multiply short vectors loaded from `a` and `b` and add their results to a vector of accumulators.
14
+
To execute the loop using SIMD (Single Instruction Multiple Data) instructions, you have to unroll the loop. Rather than evaluating the loop as written -- adding element-wise products to a single accumulator one after the other -- you can multiply short vectors loaded from `a` and `b` and add their results to a vector of accumulators.
15
15
16
16
Most modern CPUs found in laptops or desktops have the AVX instruction set, which allows them to operate on 256 bit vectors -- meaning the vectors can hold 4 double precision (64 bit) floats. Some have the AVX512 instruction set, which increases the vector size to 512 bits, and also adds many new instructions that make vectorizing easier. To be general across CPUs and data types, I'll refer to the number of elements in the vectors with `W`. I'll also refer to unrolling a loop by a factor of `W` and loading vectors from it as "vectorizing" that loop.
17
17
@@ -21,20 +21,20 @@ This means that if we used a single vector to accumulate a product, we'd only ge
21
21
If we had 8 accumulators, then theoretically we could perform two per clock cycle, and after the 4th cycle, our first operations are done so that we can reuse them.
22
22
23
23
However, there is another bottle neck: we can only perform 2 aligned loads per clock cycle (or 1 unaligned load). [Alignment here means with respect to a memory address boundary, if your vectors are 256 bits, then a load/store is aligned if it is with respect to a memory address that is an integer multiple of 32 bytes (256 bits = 32 bytes).]
24
-
Thus, in 4 clock cycles, we can do up to 8 loads. But each `fma` requires 2 loads, meaning we are limited to 4 of them per 4 clock cyles, and any unrolling beyond 4 gives us no benefit.
24
+
Thus, in 4 clock cycles, we can do up to 8 loads. But each `fma` requires 2 loads, meaning we are limited to 4 of them per 4 clock cycles, and any unrolling beyond 4 gives us no benefit.
25
25
26
26
Double precision benchmarks pitting Julia's builtin dot product, and code compiled with a variety of compilers:
What we just described is the core of the approach used by all these compilers. The variation in results is explained mostly by how they handle vectors with lengths that are not an integer multiple of `W`. I ran these on a computer with AVX512 so that `W = 8`. LLVM, the backend compiler of both Julia and Clang, shows rapid performance degradation as `N % 4W` increases, where `N` is the length of the vectors.
29
-
This is because, to handle the remainder, it uses a scalar loop that runs as written: multiply and add single elements, one after the other.
29
+
This is because, to handle the remainder, it uses a scalar loop that runs as written: multiply and add single elements, one after the other.
30
30
31
31
Initially, GCC (gfortran) stumbled in throughput, because it does not use separate accumulation vectors by default except on Power, even with `-funroll-loops`.
32
32
I compiled with the flags `-fvariable-expansion-in-unroller --param max-variable-expansions-in-unroller=4` to allow for 4 accumulation vectors, yielding good performance.
33
33
34
34
The Intel compilers have a secondary vectorized loop without any additional unrolling that masks off excess lanes beyond `N` (for when `N` isn't an integer multiple of `W`).
35
35
LoopVectorization uses `if/ifelse` checks to determine how many extra vectors are needed, the last of which is masked.
36
36
37
-
Neither GCC nor LLVM use masks (without LoopVectorization's assitance).
37
+
Neither GCC nor LLVM use masks (without LoopVectorization's assistance).
38
38
39
39
I am not certain, but I believe Intel and GCC check for the vector's alignment, and align them if necessary. Julia guarantees that the start of arrays beyond a certain size are aligned, so this is not an optimization I have implemented. But it may be worthwhile for handling large matrices with a number of rows that isn't an integer multiple of `W`. For such matrices, the first column may be aligned, but the next will not be.
40
40
@@ -56,5 +56,3 @@ For this reason, LoopVectorization now unrolls by 8 -- it decides how much to un
56
56
This algorithm may need refinement, because Julia (without LoopVectorization) only unrolls by 4, yet achieves roughly the same performance as LoopVectorization at multiples of `4W = 32`, although performance declines rapidly from there due to the slow scalar loop. Performance for most is much higher -- more GFLOPS -- than the normal dot product, but still under half of the CPU's potential 131.2 GFLOPS, suggesting that some other bottlenecks are preventing the core from attaining 2 fmas per clock cycle.
57
57
Note also that `8W = 64`, so we don't really have enough iterations of the loop to amortize the overhead of performing the reductions of all these vectors into a single scalar.
58
58
By the time the vectors are long enough to do this, we'll start running into memory bandwidth bottlenecks.
Copy file name to clipboardExpand all lines: docs/src/examples/matrix_multiplication.md
+3-7Lines changed: 3 additions & 7 deletions
Display the source diff
Display the rich diff
Original file line number
Diff line number
Diff line change
@@ -15,7 +15,7 @@ function A_mul_B!(C, A, B)
15
15
end
16
16
end
17
17
```
18
-
and this can handle all transposed/not-tranposed permutations. LoopVectorization will change loop orders and strategy as appropriate based on the types of the input matrices. For each of the others, I wrote separate functions to handle each case.
18
+
and this can handle all transposed/not-transposed permutations. LoopVectorization will change loop orders and strategy as appropriate based on the types of the input matrices. For each of the others, I wrote separate functions to handle each case.
19
19
Letting all three matrices be square and `Size` x `Size`, we attain the following benchmark results:
The optimal pattern for `π = π * πα΅` is almost identical to that for `π = π * π`. Yet, gfortran's `matmul` intrinsic stumbles, surprisingly doing much worse than gfortran + loops, and almost certainly worse than allocating memory for `πα΅` and creating the explicit copy.
26
26
27
-
ifort did equally well whethor or not `π` was transposed, while LoopVectorization's performance degraded slightly faster as a function of size in the transposed case, because strides between memory accesses are larger when `π` is transposed. But it still performed best of all the compiled loops over this size range, losing out to MKL and eventually OpenBLAS.
27
+
ifort did equally well whether or not `π` was transposed, while LoopVectorization's performance degraded slightly faster as a function of size in the transposed case, because strides between memory accesses are larger when `π` is transposed. But it still performed best of all the compiled loops over this size range, losing out to MKL and eventually OpenBLAS.
28
28
icc interestingly does better when it is transposed.
29
29
30
-
GEMM is easiest when the matrix `π` is not tranposed (assuming column-major memory layouts), because then you can sum up columns of `π` to store into `π`. If `π` were transposed, then we cannot efficiently load contiguous elements from `π` that can best stored directly in `π`. So for `π = πα΅ * π`, contiguous vectors along the `k`-loop have to be reduced, adding some overhead.
30
+
GEMM is easiest when the matrix `π` is not transposed (assuming column-major memory layouts), because then you can sum up columns of `π` to store into `π`. If `π` were transposed, then we cannot efficiently load contiguous elements from `π` that can best stored directly in `π`. So for `π = πα΅ * π`, contiguous vectors along the `k`-loop have to be reduced, adding some overhead.
Packing is critical for performance here. LoopVectorization does not pack, therefore it is well behind MKL and OpenBLAS, which do. Eigen packs, but is poorly optimized for this CPU architecture.
33
33
34
34
When both `π` and ` π` are transposed, we now have `π = πα΅ * πα΅ = (π * π)α΅`.
Julia, Clang, and gfortran all struggled to vectorize this, because none of the matrices share a contiguous access: `M` for `π`, `K` for `πα΅`, and `N` for `πα΅`. However, LoopVectorization and all the specialized matrix multiplication functions managed to do about as well as normal; transposing while storing the results takes negligible amounts of time relative to the matrix multiplication itself.
0 commit comments