LoopVectorization.jl

https://juliasimd.github.io/LoopVectorization.jl/latest/
Admin
LoopVectorization.jl provides macros and functions for Julia that optimize performance through SIMD...
Tokens:30,441
Snippets:222
Trust Score:9.7
Update:1 month ago
Show doc for...
Tokens:
Raw
### Install LoopVectorization.jl

Source: https://juliasimd.github.io/LoopVectorization.jl/latest/getting_started

Use the Julia package manager to add the library.

```julia
using Pkg
Pkg.add("LoopVectorization")
```

--------------------------------

### Benchmarking vmapnt! Performance

Source: https://juliasimd.github.io/LoopVectorization.jl/latest/api

Example comparing standard map!, vmap!, and vmapnt! performance using BenchmarkTools.

```julia
julia> f(x, y) = exp(-0.5abs2(x - y))
using LoopVectorization, BenchmarkTools

julia> @benchmark map!(f, $z, $x, $y)
x = rand(10^8); y = rand(10^8); z = similar(x);

julia> @benchmark vmap!(f, $z, $x, $y)
f (generic function with 1 method)

julia> @benchmark vmapnt!(f, $z, $x, $y)
BenchmarkTools.Trial:
  memory estimate:  0 bytes
  allocs estimate:  0
  --------------
  minimum time:     439.613 ms (0.00% GC)
  median time:      440.729 ms (0.00% GC)
  mean time:        440.695 ms (0.00% GC)
  maximum time:     441.665 ms (0.00% GC)
  --------------
  samples:          12
  evals/sample:     1
```

--------------------------------

### Optimize Loops with @turbo

Source: https://juliasimd.github.io/LoopVectorization.jl/latest/getting_started

Example of using @turbo for a matrix-vector product with rectangular iteration spaces.

```julia
using LoopVectorization 

function mvp(P, basis, coeffs::Vector{T}) where {T}
    C = length(coeffs)
    A = size(P, 1)
    p = zero(T)
    @turbo for c ∈ 1:C
        pc = coeffs[c]
        for a = 1:A
            pc *= P[a, basis[a, c]]
        end
        p += pc
    end
    p
end

maxdeg = 20; nbasis = 1_000; dim = 15;
r = 1:maxdeg+1
basis = rand(r, (dim, nbasis));
coeffs = rand(T, nbasis);
P = rand(T, dim, maxdeg+1);

mvp(P, basis, coeffs)
```

--------------------------------

### Example of a load operation

Source: https://juliasimd.github.io/LoopVectorization.jl/latest/devdocs/lowering

A basic representation of a load operation within the lowering process.

```julia
vload(vptr_A, (i,j,k))
```

--------------------------------

### Optimizing matrix multiplication with @turbo

Source: https://juliasimd.github.io/LoopVectorization.jl/latest/api

Example of using @turbo to optimize a matrix multiplication function. It assumes loop iterations are independent and does not perform bounds checks.

```julia
function AmulB!(C, A, B)
    @turbo for m ∈ indices((A,C), 1), n ∈ indices((B,C), 2) # indices((A,C),1) == axes(A,1) == axes(C,1)
        Cₘₙ = zero(eltype(C))
        for k ∈ indices((A,B), (2,1)) # indices((A,B), (2,1)) == axes(A,2) == axes(B,1)
            Cₘₙ += A[m,k] * B[k,n]
        end
        C[m,n] = Cₘₙ
    end
end
```

--------------------------------

### LoopSet loop bounds

Source: https://juliasimd.github.io/LoopVectorization.jl/latest/devdocs/constructing_loopsets

Example of loop bounds representation in a LoopSet summary.

```julia
(LoopVectorization.StaticLowerUnitRange{0}(M), LoopVectorization.StaticLowerUnitRange{0}(N), LoopVectorization.StaticLowerUnitRange{0}(K))
```

--------------------------------

### Matrix multiplication with HybridArrays and LoopVectorization

Source: https://juliasimd.github.io/LoopVectorization.jl/latest/examples/array_interface

Setup of hybrid arrays containing 3x3 static matrices and a custom bmul! function using @turbo for optimized multiplication.

```julia
using HybridArrays, StaticArrays, LoopVectorization, BenchmarkTools

A_static = [@SMatrix(rand(3,3)) for i in 1:32, j in 1:32];
B_static = [@SMatrix(rand(3,3)) for i in 1:32, j in 1:32];
C_static = similar(A_static);

A_hybrid = HybridArray{Tuple{StaticArrays.Dynamic(),StaticArrays.Dynamic(),3,3}}(permutedims(reshape(reinterpret(Float64, A_static), (3,3,size(A_static)...)), (3,4,1,2)));
B_hybrid = HybridArray{Tuple{StaticArrays.Dynamic(),StaticArrays.Dynamic(),3,3}}(permutedims(reshape(reinterpret(Float64, B_static), (3,3,size(B_static)...)), (3,4,1,2)));
C_hybrid = HybridArray{Tuple{StaticArrays.Dynamic(),StaticArrays.Dynamic(),3,3}}(permutedims(reshape(reinterpret(Float64, C_static), (3,3,size(C_static)...)), (3,4,1,2)));

# C is M x N x I x J
# A is M x K x I x L
# B is K x N x L x J
function bmul!(C, A, B)
    @turbo for n in axes(C,2), m in axes(C,1), j in axes(C,4), i in axes(C,3)
        Cmnji = zero(eltype(C))
        for k in axes(B,1), l in axes(B,3)
            Cmnji += A[m,k,i,l] * B[k,n,l,j]
        end
        C[m,n,i,j] = Cmnji
    end
end
```

--------------------------------

### LoopSet operation definition

Source: https://juliasimd.github.io/LoopVectorization.jl/latest/devdocs/constructing_loopsets

Example of an operation definition within a LoopSet summary.

```julia
:LoopVectorization, :getindex, LoopVectorization.OperationStruct(0x0000000000000013, 0x0000000000000000, 0x0000000000000000, 0x0000000000000000, LoopVectorization.memload, 0x02, 0x03)
```

--------------------------------

### Relative prime count with @tturbo

Source: https://juliasimd.github.io/LoopVectorization.jl/latest/examples/multithreading

Example from ThreadsX.jl README using @tturbo to count relative primes. Demonstrates multithreading with a reduction over a loop.

```julia
function relative_prime_count(x, N)
    c = 0
    @tturbo for i ∈ 1:N
        c += gcd(x, i) == 1
    end
    c
end
```

--------------------------------

### LoopSet array reference structure

Source: https://juliasimd.github.io/LoopVectorization.jl/latest/devdocs/constructing_loopsets

Example of an array reference structure used in LoopSet summaries.

```julia
LoopVectorization.ArrayRefStruct(0x0000000000000101, 0x0000000000000102, 0xffffffffffffe03b)
```

--------------------------------

### Julia: Matrix multiplication kernel with @tturbo

Source: https://juliasimd.github.io/LoopVectorization.jl/latest/examples/multithreading

Computes C = A * B using @tturbo with explicit loops over indices. Uses a local accumulator Cmn and assigns to C[m,n]. This is a compute-bound example compared against MKL, OpenBLAS, Tullio.jl, and Octavian.jl.

```julia
function A_mul_B!(C, A, B)
    @tturbo for n ∈ indices((C,B), 2), m ∈ indices((C,A), 1)
        Cmn = zero(eltype(C))
        for k ∈ indices((A,B), (2,1))
            Cmn += C[m,k] * B[k,n]
        end
        C[m,n] = Cmn
    end
end
```

--------------------------------

### Benchmark vmap!, vmapnt!, and vmapntt! against map!

Source: https://juliasimd.github.io/LoopVectorization.jl/latest/vectorized_convenience_functions

Benchmarks map!, vmap!, vmapnt!, and vmapntt! on a large array. vmapnt! uses non-temporal stores to avoid cache pollution, and vmapntt! is threaded. The output shows performance improvements and the number of threads.

```julia
julia> using LoopVectorization, BenchmarkTools

julia> f(x,y) = exp(-0.5abs2(x - y))
f (generic function with 1 method)

julia> x = rand(10^8); y = rand(10^8); z = similar(x);

julia> @benchmark map!(f, $z, $x, $y)
BenchmarkTools.Trial:
  memory estimate:  0 bytes
  allocs estimate:  0
  --------------
  minimum time:     442.614 ms (0.00% GC)
  median time:      443.750 ms (0.00% GC)
  mean time:        443.664 ms (0.00% GC)
  maximum time:     444.730 ms (0.00% GC)
  --------------
  samples:          12
  evals/sample:     1

julia> @benchmark vmap!(f, $z, $x, $y)
BenchmarkTools.Trial:
  memory estimate:  0 bytes
  allocs estimate:  0
  --------------
  minimum time:     177.257 ms (0.00% GC)
  median time:      177.380 ms (0.00% GC)
  mean time:        177.423 ms (0.00% GC)
  maximum time:     177.956 ms (0.00% GC)
  --------------
  samples:          29
  evals/sample:     1

julia> @benchmark vmapnt!(f, $z, $x, $y)
BenchmarkTools.Trial:
  memory estimate:  0 bytes
  allocs estimate:  0
  --------------
  minimum time:     143.521 ms (0.00% GC)
  median time:      143.639 ms (0.00% GC)
  mean time:        143.645 ms (0.00% GC)
  maximum time:     143.821 ms (0.00% GC)
  --------------
  samples:          35
  evals/sample:     1

julia> Threads.nthreads()
36

julia> @benchmark vmapntt!(f, $z, $x, $y)
BenchmarkTools.Trial:
  memory estimate:  25.69 KiB
  allocs estimate:  183
  --------------
  minimum time:     30.065 ms (0.00% GC)
  median time:      30.130 ms (0.00% GC)
  mean time:        30.146 ms (0.00% GC)
  maximum time:     31.277 ms (0.00% GC)
  --------------
  samples:          166
  evals/sample:     1
```

--------------------------------

### Inspect LoopSet operations

Source: https://juliasimd.github.io/LoopVectorization.jl/latest/devdocs/reference

Demonstrates how to create a LoopSet from a Julia expression and inspect the resulting operations.

```julia
julia> using LoopVectorization

julia> AmulBq = :(for m ∈ 1:M, n ∈ 1:N
           C[m,n] = zero(eltype(B))
           for k ∈ 1:K
               C[m,n] += A[m,k] * B[k,n]
           end
       end);

julia> lsAmulB = LoopVectorization.LoopSet(AmulBq);

julia> LoopVectorization.operations(lsAmulB)
6-element Vector{LoopVectorization.Operation}:
 var"##RHS#245" = var"##zero#246"
 C[m, n] = var"##RHS#245"
 var"##tempload#248" = A[m, k]
 var"##tempload#249" = B[k, n]
 var"##RHS#245" = LoopVectorization.vfmadd(var"##tempload#248", var"##tempload#249", var"##RHS#245")
 var"##RHS#245" = LoopVectorization.identity(var"##RHS#245")
```

--------------------------------

### Benchmark vmapreduce against mapreduce

Source: https://juliasimd.github.io/LoopVectorization.jl/latest/vectorized_convenience_functions

Compares vmapreduce and mapreduce for reducing with hypot and +. vmapreduce applies f(a[i], b[i], c[i]) and reduces with op.

```julia
julia> using LoopVectorization, BenchmarkTools

julia> x = rand(127); y = rand(127);

julia> @btime vmapreduce(hypot, +, $x, $y)
  191.420 ns (0 allocations: 0 bytes)
96.75538300513509

julia> @btime mapreduce(hypot, +, $x, $y)
  1.777 μs (5 allocations: 1.25 KiB)
96.75538300513509
```

--------------------------------

### Sample Input and Expected Output for DateTime Scaling

Source: https://juliasimd.github.io/LoopVectorization.jl/latest/examples/datetime_arrays

Defines the sample input vector of DateTime objects and the corresponding expected output vector of scaled Float64 values.

```julia
using Dates

sample_input = [
    Dates.DateTime(2021, 5, 5, 10, 0, 0),
    Dates.DateTime(2021, 5, 5, 10, 5, 15),
    Dates.DateTime(2021, 5, 6, 10, 0, 0),
    Dates.DateTime(2021, 5, 6, 10, 5, 15),
    Dates.DateTime(2021, 5, 7, 10, 0, 20),
]

expected_output = [
    0.0,
    0.0018227057053581761,
    0.499942136326814,
    0.5017648420321722,
    1.0,
]
```

--------------------------------

### Benchmark Broadcasting Performance

Source: https://juliasimd.github.io/LoopVectorization.jl/latest/getting_started

Compare standard broadcasting performance against @turbo using BenchmarkTools.

```julia
julia> using LoopVectorization, BenchmarkTools

julia> M, K, N = 47, 73, 7;

julia> A = rand(M, K);

julia> b = rand(K);

julia> c = rand(M);

julia> d = rand(1,K,N);

julia> #You can use a LowDimArray when you have a leading dimension of size 1.
       ldad = LowDimArray{(false,true,true)}(d);

julia> E1 = Array{Float64}(undef, M, K, N);

julia> E2 = similar(E1);

julia> @benchmark @. $E1 = exp($A - $b' +    $d) * $c
BenchmarkTools.Trial: 
  memory estimate:  112 bytes
  allocs estimate:  5
  --------------
  minimum time:     224.142 μs (0.00% GC)
  median time:      225.773 μs (0.00% GC)
  mean time:        229.146 μs (0.00% GC)
  maximum time:     289.601 μs (0.00% GC)
  --------------
  samples:          10000
  evals/sample:     1

julia> @benchmark @turbo @. $E2 = exp($A - $b' + $ldad) * $c
BenchmarkTools.Trial: 
  memory estimate:  0 bytes
  allocs estimate:  0
  --------------
  minimum time:     19.666 μs (0.00% GC)
  median time:      19.737 μs (0.00% GC)
  mean time:        19.759 μs (0.00% GC)
  maximum time:     29.906 μs (0.00% GC)
  --------------
  samples:          10000
  evals/sample:     1

julia> E1 ≈ E2
true
```

--------------------------------

### Dot product implementation with @turbo

Source: https://juliasimd.github.io/LoopVectorization.jl/latest/examples/dot_product

A basic implementation of a dot product using the @turbo macro to enable SIMD vectorization.

```julia
function jdotavx(a, b)
    s = zero(eltype(a))
    @turbo for i ∈ eachindex(a, b)
        s += a[i] * b[i]
    end
    s
end
```

--------------------------------

### Benchmark sequential vs turbo performance

Source: https://juliasimd.github.io/LoopVectorization.jl/latest/examples/datetime_arrays

Compares the execution time of sequential and turbo-optimized scaling functions using BenchmarkTools.

```julia
julia> using BenchmarkTools

julia> data_100000 = generate_timestamps(100000);

julia> data_200000 = generate_timestamps(200000);

julia> @benchmark scale_timeseries_sequential(data_100000)
BenchmarkTools.Trial: 10000 samples with 1 evaluation.
 Range (min … max):  318.864 μs … 967.760 μs  ┊ GC (min … max): 0.00% … 40.41%
 Time  (median):     321.291 μs               ┊ GC (median):    0.00%
 Time  (mean ± σ):   332.503 μs ±  52.040 μs  ┊ GC (mean ± σ):  1.97% ±  6.98%

  █▆▅▂▂▂▁                                                       ▁
  █████████▆▆▆▅▅▅▅▅▄▄▄▄▁▁▄▁▁▃▁▁▄▄▃▃▁▃▄▁▃▁▁▁▃▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▅█▇ █
  319 μs        Histogram: log(frequency) by time        701 μs <

 Memory estimate: 781.33 KiB, allocs estimate: 2.

julia> @benchmark scale_timeseries_turbo(data_100000)
BenchmarkTools.Trial: 10000 samples with 1 evaluation.
 Range (min … max):   71.942 μs … 933.400 μs  ┊ GC (min … max):  0.00% … 71.93%
 Time  (median):      87.926 μs               ┊ GC (median):     0.00%
 Time  (mean ± σ):   100.082 μs ±  89.095 μs  ┊ GC (mean ± σ):  11.63% ± 11.43%

  ▄█▃▁                                                        ▁ ▁
  ████▇▄▁▁▁▁▁▁▁▁▃▄▆▄▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁█ █
  71.9 μs       Histogram: log(frequency) by time        764 μs <

 Memory estimate: 781.33 KiB, allocs estimate: 2.

julia> @benchmark scale_timeseries_sequential(data_200000)
BenchmarkTools.Trial: 7153 samples with 1 evaluation.
 Range (min … max):  637.692 μs …   2.277 ms  ┊ GC (min … max): 0.00% … 65.01%
 Time  (median):     640.729 μs               ┊ GC (median):    0.00%
 Time  (mean ± σ):   694.282 μs ± 184.965 μs  ┊ GC (mean ± σ):  3.69% ±  8.68%

  █▆▅▃▂▁                                                        ▁
  ███████▇▅▆▆▄▄▅▁▁▁▁▁▁▆██▇▅▄▄▁▃▁▃▃▄▁▁▁▁▁▁▁▁▁▁▁▇█▇▇▅▅▄▃▄▄▁▁▁▄▄█▇ █
  638 μs        Histogram: log(frequency) by time       1.71 ms <

 Memory estimate: 1.53 MiB, allocs estimate: 2.

julia> @benchmark scale_timeseries_turbo(data_200000)
BenchmarkTools.Trial: 10000 samples with 1 evaluation.
 Range (min … max):  159.023 μs …   2.092 ms  ┊ GC (min … max):  0.00% … 50.30%
 Time  (median):     176.559 μs               ┊ GC (median):     0.00%
 Time  (mean ± σ):   230.513 μs ± 189.542 μs  ┊ GC (mean ± σ):  11.86% ± 12.80%

  █▇▅▄▄▃▂▂▁          ▁▁                                         ▂
  ██████████▇▅▅▃▁▄▁▃▁██▇▆▄▅▄▄▅▄▅▃▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁███▇▅▅▅▅▄▅▄▅█▇▇ █
  159 μs        Histogram: log(frequency) by time       1.22 ms <

 Memory estimate: 1.53 MiB, allocs estimate: 2.
```

--------------------------------

### Benchmark: dot product variants

Source: https://juliasimd.github.io/LoopVectorization.jl/latest/examples/multithreading

Benchmark comparing LinearAlgebra.dot, dot_turbo, dot_tturbo, dot_baseline, and cdot for N=10,000. Shows relative performance.

```julia
julia> N = 10_000; x = rand(N); y = rand(N);

julia> @btime dot($x, $y) # LinearAlgebra
  1.114 μs (0 allocations: 0 bytes)
2480.296446711209

julia> @btime dot_turbo($x, $y)
  761.621 ns (0 allocations: 0 bytes)
2480.296446711209

julia> @btime dot_tturbo($x, $y)
  622.723 ns (0 allocations: 0 bytes)
2480.296446711209

julia> @btime dot_baseline($x, $y)
  1.294 μs (0 allocations: 0 bytes)
2480.2964467112097

julia> @btime cdot($x, $y)
  6.109 μs (0 allocations: 0 bytes)
2480.2964467112092
```

--------------------------------

### Julia: Wrapper for C dot product using @ccall

Source: https://juliasimd.github.io/LoopVectorization.jl/latest/examples/multithreading

Calls the C function cdot3 from the OPENMPTEST library using @ccall. Passes references to complex arrays and dimensions as Clong. Returns the result from a Ref.

```julia
function cdot(x::AbstractVector{Complex{Float64}}, A::AbstractMatrix{Complex{Float64}}, y::AbstractVector{Complex{Float64}})
    c = Ref{Complex{Float64}}()
    M, N = size(A)
    @ccall OPENMPTEST.cdot3(c::Ref{Complex{Float64}}, x::Ref{Complex{Float64}}, A::Ref{Complex{Float64}}, y::Ref{Complex{Float64}}, M::Clong, N::Clong)::Cvoid
    c[]
end
```

--------------------------------

### vmapt!

Source: https://juliasimd.github.io/LoopVectorization.jl/latest/api

A threaded variant of vmap!.

```APIDOC
## vmapt!

### Description
A threaded variant of `vmap!`.

### Function Signature
```julia
vmapt!(::Function, dest, args...)
```

### Parameters
- `f` (Function) - Required - The function to apply to batches of elements.
- `dest` (AbstractArray) - Required - The destination array.
- `args...` (AbstractArray) - Required - Input arrays.
```

--------------------------------

### Julia REPL: Cache size calculations

Source: https://juliasimd.github.io/LoopVectorization.jl/latest/examples/multithreading

Calculates the number of doubles that fit in L2 cache per core and total, then divides among three matrices to estimate a matrix size for benchmarking. Uses Sys.CPU_THREADS and assumes 18 cores.

```julia
julia> doubles_per_l2 = (2 ^ 20) ÷ 8
131072

julia> total_doubles_in_l2 = doubles_per_l2 * (Sys.CPU_THREADS ÷ 2) # doubles_per_l2 * 18
2359296

julia> doubles_per_mat = total_doubles_in_l2 ÷ 3 # divide up among 3 matrices
786432

julia> sqrt(ans)
886.8100134752651
```

--------------------------------

### Matrix-Vector Multiplication with jgemvavx!

Source: https://juliasimd.github.io/LoopVectorization.jl/latest/examples/matrix_vector_ops

Implements matrix-vector multiplication using the @turbo macro for loop optimization.

```julia
function jgemvavx!(𝐲, 𝐀, 𝐱)
    @turbo for i ∈ eachindex(𝐲)
        𝐲i = zero(eltype(𝐲))
        for j ∈ eachindex(𝐱)
            𝐲i += 𝐀[i,j] * 𝐱[j]
        end
        𝐲[i] = 𝐲i
    end
end
```

--------------------------------

### vmapreduce

Source: https://juliasimd.github.io/LoopVectorization.jl/latest/api

Vectorized version of mapreduce.

```APIDOC
## vmapreduce

### Description
Vectorized version of `mapreduce`.
```

--------------------------------

### vmap and vmap!

Source: https://juliasimd.github.io/LoopVectorization.jl/latest/api

Vectorized versions of map and map!.

```APIDOC
## vmap and vmap!

### Description
Vectorized versions of `map` and `map!`.
```

--------------------------------

### Implement matrix multiplication with @turbo for StaticArrays

Source: https://juliasimd.github.io/LoopVectorization.jl/latest/examples/array_interface

Defines AmulB! and AmulB methods that use @turbo to multiply matrices, supporting MMatrix and SMatrix inputs. SMatrix inputs are converted to MMatrix to enable pointer access, and the result is converted back. This avoids writing hundreds of lines of unrolled code.

```julia
using StaticArrays, LoopVectorization

@inline function AmulB!(C, A, B)
    @turbo for n ∈ axes(C,2), m ∈ axes(C,1)
        Cmn = zero(eltype(C))
        for k ∈ axes(B,1)
            Cmn += A[m,k] * B[k,n]
        end
        C[m,n] = Cmn
    end
    C
end
@inline AmulB(A::MMatrix{M,K,T}, B::MMatrix{K,N,T}) where {M,K,N,T} = AmulB!(MMatrix{M,N,T}(undef), A, B)
@inline AmulB(A::SMatrix, B::SMatrix) = SMatrix(AmulB(MMatrix(A), MMatrix(B)))
```

--------------------------------

### OpenMP dot product in C

Source: https://juliasimd.github.io/LoopVectorization.jl/latest/examples/multithreading

C implementation of dot product using OpenMP parallel for reduction. Compilation command is provided in a comment.

```c
#include<omp.h>
//  gcc -Ofast -march=native -mprefer-vector-width=512 -fopenmp -shared -fPIC openmp.c -o libomptest.so
double dot(double* a, double* b, long N){
  double s = 0.0;
  #pragma omp parallel for reduction(+: s)
  for(long n = 0; n < N; n++){
    s += a[n]*b[n];
  }
  return s;
}
```

--------------------------------

### Benchmark: ThreadsX.sum vs relative_prime_count

Source: https://juliasimd.github.io/LoopVectorization.jl/latest/examples/multithreading

Benchmark comparing ThreadsX.sum with the @tturbo version. Shows significant performance difference due to SIMD and multithreading.

```julia
julia> @btime ThreadsX.sum(gcd(42, i) == 1 for i in 1:10_000)
  130.928 μs (3097 allocations: 240.39 KiB)
2857

julia> @btime relative_prime_count(42, 10_000)
  3.376 μs (0 allocations: 0 bytes)
2857
```

--------------------------------

### Dot product with @tturbo and baseline

Source: https://juliasimd.github.io/LoopVectorization.jl/latest/examples/multithreading

Two Julia functions for dot product: one using @tturbo and a baseline using @fastmath @inbounds @simd. Used to compare threading overhead.

```julia
function dot_tturbo(a::AbstractArray{T}, b::AbstractArray{T}) where {T <: Real}
    s = zero(T)
    @tturbo for i ∈ eachindex(a,b)
        s += a[i] * b[i]
    end
    s
end
function dotbaseline(a::AbstractArray{T}, b::AbstractArray{T}) where {T}
    s = zero(T)
    @fastmath @inbounds @simd for i ∈ eachindex(a,b)
        s += a[i]' * b[i]
    end
    s
end
```

--------------------------------

### Customizing @turbo with keyword arguments

Source: https://juliasimd.github.io/LoopVectorization.jl/latest/api

Advanced users can customize the implementation of the @turbo-annotated block using specific keyword arguments.

```julia
@turbo inline = false unroll = 2 thread = 4 body
```

--------------------------------

### Show parent operations of a reduction

Source: https://juliasimd.github.io/LoopVectorization.jl/latest/devdocs/loopset_structure

Selects the sixth operation (the vfmadd_fast reduction) and then shows its parent operations using LoopVectorization.parents. The parents are the two loads and the reduction zero initialization.

```julia
LoopVectorization.operations(lsAmulB)[6]
var"##reduction#260" = LoopVectorization.vfmadd_fast(var"##tempload#258", var"##tempload#259", var"##reduction#260")

LoopVectorization.parents(ans)
3-element Array{LoopVectorization.Operation,1}:
 var"##tempload#258" = A[m, k]
 var"##tempload#259" = B[k, n]
 var"##reduction#260" = var"##reductzero#261"
```

--------------------------------

### vmapntt and vmapntt!

Source: https://juliasimd.github.io/LoopVectorization.jl/latest/api

Threaded variants of vmapnt and vmapnt!.

```APIDOC
## vmapntt and vmapntt!

### Description
Threaded variants of `vmapnt` and `vmapnt!`.
```

--------------------------------

### Benchmark LoopVectorization vs StaticArrays matrix multiplication

Source: https://juliasimd.github.io/LoopVectorization.jl/latest/examples/array_interface

Benchmarks the custom AmulB! against StaticArrays' mul! for both SMatrix and MMatrix types, measuring time and GFLOPS. Uses BenchmarkTools, LinearAlgebra, DataFrames, and VegaLite to generate a plot saved as sarraymatmul.svg.

```julia
using BenchmarkTools, LinearAlgebra, DataFrames, VegaLite
BLAS.set_num_threads(1);

matdims(x::Integer) = (x, x, x)
matdims(x::NTuple{3}) = x
matflop(x::Integer) = 2x^3
matflop(x::NTuple{3}) = 2prod(x)

function runbenches(sr, ::Type{T}, fa = identity, fb = identity) where {T}
    bench_results = Matrix{Float64}(undef, length(sr), 4);
    for (i,s) ∈ enumerate(sr)
        M, K, N = matdims(s)
        Am = @MMatrix rand(T, M, K)
        Bm = @MMatrix rand(T, K, N)
        As = Ref(SMatrix(Am));
        Bs = Ref(SMatrix(Bm));
        Css = fa(As[]) * fb(Bs[]);
        Csl = AmulB(fa(As[]), fb(Bs[]))
        Cms = similar(Css); mul!(Cms, fa(Am), fb(Bm));
        Cml = similar(Css); AmulB!(Cml, fa(Am), fb(Bm));
        @assert Array(Css) ≈ Array(Csl) ≈ Array(Cms) ≈ Array(Cml) # Once upon a time Julia crashed on ≈ for large static arrays
        bench_results[i,1] = @belapsed $fa($As[]) * $fb($Bs[])
        bench_results[i,2] = @belapsed AmulB($fa($As[]), $fb($Bs[]))
        bench_results[i,3] = @belapsed mul!($Cms, $fa($Am), $fb($Bm))
        bench_results[i,4] = @belapsed AmulB!($Cml, $fa($Am), $fb($Bm))
        @show s, bench_results[i,:]
    end
    gflops = @. 1e-9 * matflop(sr) / bench_results
    array_type = append!(fill("Static", 2length(sr)), fill("Mutable", 2length(sr)))
    sa = fill("StaticArrays", length(sr)); lv = fill("LoopVectorization", length(sr));
    matmul_lib = vcat(sa, lv, sa, lv);
    sizes = reduce(vcat, (sr for _ ∈ 1:4))
    DataFrame(
        Size = sizes, Time = vec(bench_results), GFLOPS = vec(gflops),
        ArrayType = array_type, MatmulLib = matmul_lib, MulType = array_type .* ' ' .* matmul_lib
    )
end

df = runbenches(1:24, Float64);
df |> @vlplot(:line, x = :Size, y = :GFLOPS, color = :MulType, height=640,width=960) |> save("sarraymatmul.svg")
```

--------------------------------

### vmap!

Source: https://juliasimd.github.io/LoopVectorization.jl/latest/api

Vectorized map!, applying f to batches of elements of a (or paired batches of a, b, ...) and storing the result in destination. The function f must accept VectorizationBase.AbstractSIMD inputs.

```APIDOC
## vmap!

### Description
Vectorized-`map!`, applying `f` to batches of elements of `a` (or paired batches of `a`, `b`, ...) and storing the result in `destination`. The function `f` must accept `VectorizationBase.AbstractSIMD` inputs. Ideally, all this requires is making sure that `f` is defined to be agnostic with respect to input types, but if the function `f` contains branches or loops, more work will probably be needed.

### Function Signature
```julia
vmap!(f, destination, a::AbstractArray)
vmap!(f, destination, a::AbstractArray, b::AbstractArray, ...)
```

### Parameters
- `f` (Function) - Required - The function to apply to batches of elements. Must accept `VectorizationBase.AbstractSIMD` inputs.
- `destination` (AbstractArray) - Required - The array where results are stored.
- `a` (AbstractArray) - Required - The first input array.
- `b...` (AbstractArray) - Optional - Additional input arrays for element-wise operations.

### Notes
If `f` contains branches or loops, it may need to be rewritten to be SIMD-friendly. For example, a function `f(x) = x > 0 ? log(x) : inv(x)` can be rewritten using `IfElse.ifelse`:
```julia
using IfElse
f(x) = IfElse.ifelse(x > 0, log(x), inv(x))
```
```

--------------------------------

### vmapt

Source: https://juliasimd.github.io/LoopVectorization.jl/latest/api

A threaded variant of vmap.

```APIDOC
## vmapt

### Description
A threaded variant of `vmap`.

### Function Signature
```julia
vmapt(f, a::AbstractArray)
vmapt(f, a::AbstractArray, b::AbstractArray, ...)
```

### Parameters
- `f` (Function) - Required - The function to apply to each element or paired elements.
- `a` (AbstractArray) - Required - The first input array.
- `b...` (AbstractArray) - Optional - Additional input arrays for element-wise operations.

### Returns
A new array containing the results of applying `f` to each element or paired elements, computed in parallel.
```

--------------------------------

### Julia wrapper for OpenMP dot

Source: https://juliasimd.github.io/LoopVectorization.jl/latest/examples/multithreading

Julia wrapper using @ccall to call the compiled OpenMP dot function. Requires the shared library to be compiled first.

```julia
using Libdl; const OPENMPTEST = joinpath(pkgdir(LoopVectorization), "benchmark", "libomptest.$(Libdl.dlext)");
cdot(a::AbstractVector{Float64},b::AbstractVector{Float64}) = @ccall OPENMPTEST.dot(a::Ref{Float64}, b::Ref{Float64}, length(a)::Clong)::Float64
```

--------------------------------

### Benchmark convolution performance

Source: https://juliasimd.github.io/LoopVectorization.jl/latest/examples/multithreading

Compares the performance of the custom convlayer! function against NNlib's default conv! implementation.

```julia
julia> NNlib.conv!(out2, img, kern, dcd);

julia> convlayer!(out1, img, kern, dcd);

julia> out1 ≈ out2
true

julia> @benchmark convlayer!($out1, $img, $kern, $dcd)
BenchmarkTools.Trial:
  memory estimate:  0 bytes
  allocs estimate:  0
  --------------
  minimum time:     5.377 ms (0.00% GC)
  median time:      5.432 ms (0.00% GC)
  mean time:        5.433 ms (0.00% GC)
  maximum time:     5.682 ms (0.00% GC)
  --------------
  samples:          920
  evals/sample:     1

julia> @benchmark NNlib.conv!($out2, $img, $kern, $dcd)
BenchmarkTools.Trial:
  memory estimate:  675.02 MiB
  allocs estimate:  195
  --------------
  minimum time:     182.749 ms (0.00% GC)
  median time:      190.472 ms (0.60% GC)
  mean time:        197.527 ms (4.98% GC)
  maximum time:     300.536 ms (35.82% GC)
  --------------
  samples:          26
  evals/sample:     1
```

--------------------------------

### Benchmarking HybridArrays vs StaticArrays multiplication

Source: https://juliasimd.github.io/LoopVectorization.jl/latest/examples/array_interface

Benchmark results comparing the performance of the custom bmul! function with HybridArrays against standard mul! with StaticArrays.

```julia
julia> @benchmark bmul!($C_hybrid, $A_hybrid, $B_hybrid)
BenchmarkTools.Trial:
  memory estimate:  0 bytes
  allocs estimate:  0
  --------------
  minimum time:     15.550 μs (0.00% GC)
  median time:      15.663 μs (0.00% GC)
  mean time:        15.685 μs (0.00% GC)
  maximum time:     50.286 μs (0.00% GC)
  --------------
  samples:          10000
  evals/sample:     1
  
julia> @benchmark mul!($C_static, $A_static, $B_static)
BenchmarkTools.Trial:
  memory estimate:  336 bytes
  allocs estimate:  6
  --------------
  minimum time:     277.736 μs (0.00% GC)
  median time:      278.035 μs (0.00% GC)
  mean time:        278.310 μs (0.00% GC)
  maximum time:     299.259 μs (0.00% GC)
  --------------
  samples:          10000
  evals/sample:     1

julia> all(I -> C_hybrid[Tuple(I)[1],Tuple(I)[2],:,:] ≈ C_static[I], CartesianIndices(C_static))
true

julia> length(C_hybrid) * size(B_hybrid,1) * size(B_hybrid,3) * 2e-9 / 15.55e-6 # GFLOPS loops + hybrid arrays
113.79241157556271

julia> length(C_hybrid) * size(B_hybrid,1) * size(B_hybrid,3) * 2e-9 / 277.736e-6 # GFLOPS LinearAlgebra.mul! + StaticArrays
6.371057407034018
```

--------------------------------

### vreduce

Source: https://juliasimd.github.io/LoopVectorization.jl/latest/api

Vectorized version of reduce.

```APIDOC
## vreduce

### Description
Vectorized version of `reduce`.
```

--------------------------------

### Benchmark vfilter against filter

Source: https://juliasimd.github.io/LoopVectorization.jl/latest/vectorized_convenience_functions

Compares Base.filter and LoopVectorization.vfilter on a random vector. vfilter requires LLVM 7 or greater and is likely faster only on AVX512 CPUs due to the compressed store intrinsic.

```julia
julia> using LoopVectorization, BenchmarkTools

julia> x = rand(997);

julia> y1 = filter(a -> a > 0.7, x);

julia> y2 = vfilter(a -> a > 0.7, x);

julia> y1 == y2
true

julia> @benchmark filter(a -> a > 0.7, $x)
BenchmarkTools.Trial:
  memory estimate:  7.94 KiB
  allocs estimate:  1
  --------------
  minimum time:     955.389 ns (0.00% GC)
  median time:      1.050 μs (0.00% GC)
  mean time:        1.191 μs (9.72% GC)
  maximum time:     82.799 μs (94.92% GC)
  --------------
  samples:          10000
  evals/sample:     18

julia> @benchmark vfilter(a -> a > 0.7, $x)
BenchmarkTools.Trial:
  memory estimate:  7.94 KiB
  allocs estimate:  1
  --------------
  minimum time:     477.487 ns (0.00% GC)
  median time:      575.166 ns (0.00% GC)
  mean time:        711.526 ns (17.87% GC)
  maximum time:     9.257 μs (79.17% GC)
  --------------
  samples:          10000
  evals/sample:     193
```

--------------------------------

### vmapntt!

Source: https://juliasimd.github.io/LoopVectorization.jl/latest/api

A threaded variant of vmapnt!.

```APIDOC
## vmapntt!

### Description
A threaded variant of `vmapnt!`.

### Function Signature
```julia
vmapntt!(::Function, dest, args...)
```

### Parameters
- `f` (Function) - Required - The function to apply to batches of elements.
- `dest` (AbstractArray) - Required - The destination array.
- `args...` (AbstractArray) - Required - Input arrays.
```

--------------------------------

### Writing SIMD-compatible functions with IfElse

Source: https://juliasimd.github.io/LoopVectorization.jl/latest/api

Functions containing branches must be rewritten to be agnostic to input types. Use IfElse.ifelse to handle conditional logic.

```julia
f(x) = x > 0 ? log(x) : inv(x)
```

```julia
using IfElse
f(x) = IfElse.ifelse(x > 0, log(x), inv(x))
```

--------------------------------

### Sequential DateTime Scaling Implementation

Source: https://juliasimd.github.io/LoopVectorization.jl/latest/examples/datetime_arrays

A sequential implementation of the timeseries scaling function that iterates over a Vector of Dates.DateTime.

```julia
using Dates

function scale_timeseries_sequential(data::Vector{Dates.DateTime})
  out = similar(data, Float64)
  ϕ = (data[lastindex(data)] - data[1]).value

  @inbounds for i ∈ eachindex(data)
      out[i] = (data[i] - data[1]).value / ϕ
  end

  return out
end
```

--------------------------------

### vmap

Source: https://juliasimd.github.io/LoopVectorization.jl/latest/api

SIMD-vectorized map, applying f to each element of a (or paired elements of a, b, ...) and returning a new array.

```APIDOC
## vmap

### Description
SIMD-vectorized `map`, applying `f` to each element of `a` (or paired elements of `a`, `b`, ...) and returning a new array.

### Function Signature
```julia
vmap(f, a::AbstractArray)
vmap(f, a::AbstractArray, b::AbstractArray, ...)
```

### Parameters
- `f` (Function) - Required - The function to apply to each element or paired elements.
- `a` (AbstractArray) - Required - The first input array.
- `b...` (AbstractArray) - Optional - Additional input arrays for element-wise operations.

### Returns
A new array containing the results of applying `f` to each element or paired elements.
```
Tokens:
Raw