Hi,
I’m trying to optimize a loop and I get a warning about `LoopVectorization.check_args` on your inputs failed.
Here is my MWE:
psi = rand(ComplexF64,4,4)
dpsi = zeros(ComplexF64,4,4)
function test_turbo(dpsi::Array{ComplexF64,2}, psi::Array{ComplexF64,2}, t::Float64)
@turbo for pos in 1:length(psi)
dpsi[pos] = -2psi[pos]
end
end
julia> test_turbo(dpsi, psi , 0.0)
┌ Warning: #= /path/to/file/test.jl:70 =#:
│ `LoopVectorization.check_args` on your inputs failed; running fallback `@inbounds @fastmath` loop instead.
│ Use `warn_check_args=false`, e.g. `@turbo warn_check_args=false ...`, to disable this warning.
└ @ Main ~/.julia/packages/LoopVectorization/ndGJi/src/condense_loopset.jl:825
It looks like I’m missing something very basic.
Any help will be welcome!
This is a little bit inelegant. But we don’t have better tools at present.
In general, SOA also helps LLVM to generate better simd code. If your data could be stored in StructArray from the very beginning. You might get better performance even without LV.jl. Just a example:
julia> a = randn(ComplexF64, 4096); b = randn(ComplexF64, 4096); c = similar(a, Float64);
julia> f(c, a, b) = c .= abs2.(a .* b) .+ abs2.(a .+ b).^3
f (generic function with 2 methods)
julia> @btime f($c, $a, $b);
3.650 μs (0 allocations: 0 bytes)
julia> sa = StructArray(a); sb = StructArray(b);
julia> @btime f($c, $sa, $sb);
2.267 μs (0 allocations: 0 bytes)
Here is my example,
I worked with ghost cells so I also used OffsetArrays.jl for helping me with the indices.
There is a least factor 10 between using @turbo or not.
Unfortunately, it all requires reinterpret and manually writing the complex operations.
The LoopVectorization rewrite will eliminate the need for this, but that’s probably at least a year away.
In the mean time, PRs for better documentation are always welcome.
Depends on the rest of your code.
All else equal, StructArrays is better because that is the optimal data layout for SIMD. LV will need extra shuffles to permute the data otherwise.
But is the rest of your code prefers regular arrays, shuffling a bit within the loop could be faster.
Hi, I’ve been having the same issue and thought it was also because of the complex arrays I’m dealing with. However, I found the same kind of error (check_args fail) in loops that deal with only real input. For example, this loop doesn’t work.
Can you give me code to copy/paste so I can run the above example?
@tturbo for I in CartesianIndices(v) #This loop raises a warning and falls back to @inbounds @fastmath
ix0, iy0, iz0 = Tuple(I)
I’m not sure if that’d work. I don’t remember adding support for it Tuple(I), anyway.
So you could try replacing it with the individual loops manually, i.e.
@tturbo for iz0 = axes(z,3), iy0 = axes(z,2), ix0 = axes(z,1)
...
r[ix0, iy0, iz0] = ...
end
It is a bit long to put here completely I think but I’ll add initialization with some random values
grid_size = (128,128,128)
function x_vec(field::AbstractArray, box_size::SVector{3,T}, box_min::SVector{3,T}) where T<:Real
dims = [size(field)...]
cell_size = box_size ./ dims
Tuple(map(T, collect(box_min[i] + 0.5 * cell_size[i]:cell_size[i]:box_min[i] + box_size[i])) for i in 1:3)
end #func
r = ones(Float32, grid_size...)
v = similar(r)
f = similar(r)
box_size = @SVector [ 1000f0, 1000f0, 1000f0]
box_min = @SVector [ 0f0, 0f0, 0f0]
x_vec = x_vec(f, box_size, box_min) # I realize the variable naming is not good, this is normally not an issue since x_vec is called differently outside the function.
β = 1f0
damping_factor = 0.5f0
niterations = 4
los = nothing # or [0,0,1] should work with both
residual!(r, v, f, x_vec, box_size, box_min, β, damping_factor niterations; los=los)
I will also test looping as you suggest, but I wanted to use CartesianIndices so I make sure I don’t mix things up and loop the dimensions in the wrong (slow) order.
Would
LoopVectorization.jl picks the order for you. If you place them in a slow order, it’ll swap them to a fast one. The original order they were in doesn’t actually matter/impact what LoopVectorization.jl does.
Note that CartesianIndices will often go in the wrong order, e.g.
but now when using @tturbo I get an out of memory error that I do not get if I uncomment the @threads for and unpack the ix0,, iy0, iz0, have you encountered this before? Should the vectorized code use more memory? Or has it something to do with the compilation?
Can you try again with VectorizationBase 0.12.54 (which was just released)?
I get
julia> @benchmark prolong!($v1h, $v2h)
BenchmarkTools.Trial: 303 samples with 1 evaluation.
Range (min … max): 16.214 ms … 20.038 ms ┊ GC (min … max): 0.00% … 0.00%
Time (median): 16.333 ms ┊ GC (median): 0.00%
Time (mean ± σ): 16.511 ms ± 556.707 μs ┊ GC (mean ± σ): 0.00% ± 0.00%
▅█▆▂
▄█████▇▆▄▄▄▅▁▁▁▄▄▁▅▅▄▁▅▁▆▁▁▁▁▅▄▁▁▄▁▁▁▄▁▁▁▁▁▁▁▁▁▄▁▁▁▅▁▆▄▄▄▄▄▄ ▆
16.2 ms Histogram: log(frequency) by time 18.8 ms <
Memory estimate: 0 bytes, allocs estimate: 0.
julia> @benchmark prolong_threads!($v1h, $v2h)
BenchmarkTools.Trial: 199 samples with 1 evaluation.
Range (min … max): 24.969 ms … 27.071 ms ┊ GC (min … max): 0.00% … 0.00%
Time (median): 25.063 ms ┊ GC (median): 0.00%
Time (mean ± σ): 25.136 ms ± 266.319 μs ┊ GC (mean ± σ): 0.00% ± 0.00%
▄▇██▆▄
▇██████▇▅▆▅▅▅▁▁▇▆▁▁▅▁▅▁▁▅▅▁▁▁▆▅▁▅▁▁▆▁▁▅▁▁▁▁▁▅▁▁▁▁▁▁▁▁▁▁▁▁▁▅▅ ▅
25 ms Histogram: log(frequency) by time 26.3 ms <
Memory estimate: 10.59 KiB, allocs estimate: 109.
After moving the @inbounds inside of the @threads for the @threads version (otherwise, it’s ignored).
You could probably improve performance by splitting the loop up to avoid those branches on indexes.
In theory, LoopVectorization should do that automatically, but it’s not something I’d gotten around to implementing.
I’ll very likely implement it in the rewrite, since it is a pattern that shows up so often.
Vectorizing those branches is very expensive, because it makes the loads/stores that use them go from using fairly cheap vector moves, to requiring expensive gather/scatter instructions.
julia> using LinuxPerf
julia> foreachf(f::F, N::Int, arg1::A, args::Vararg{Any,K}) where {F,A,K} = foreach(_ -> f(arg1, args...), 1:N)
foreachf (generic function with 1 method)
julia> @time @pstats "(cpu-cycles,task-clock),(instructions,branch-instructions,branch-misses), (L1-dcache-load-misses, L1-dcache-loads, cache-misses, cache-references)" begin
foreachf(prolong!, 100, v1h, v2h)
end
1.689482 seconds (32.32 k allocations: 1.924 MiB, 2.14% compilation time)
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━
┌ cpu-cycles 1.01e+11 33.3% # 3.4 cycles per ns
└ task-clock 2.93e+10 33.3% # 29.3 s
┌ instructions 7.86e+09 66.7% # 0.1 insns per cycle
│ branch-instructions 5.10e+07 66.7% # 0.6% of insns
└ branch-misses 1.68e+05 66.7% # 0.3% of branch insns
┌ L1-dcache-load-misses 1.29e+09 33.3% # 53.2% of dcache loads
│ L1-dcache-loads 2.42e+09 33.3%
│ cache-misses 9.56e+08 33.3% # 96.3% of cache refs
└ cache-references 9.93e+08 33.3%
aggregated from 18 threads
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━
julia> @time @pstats "(cpu-cycles,task-clock),(instructions,branch-instructions,branch-misses), (L1-dcache-load-misses, L1-dcache-loads, cache-misses, cache-references)" begin
foreachf(prolong_threads!, 100, v1h, v2h)
end
2.604119 seconds (43.48 k allocations: 2.966 MiB, 1.39% compilation time)
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━
┌ cpu-cycles 1.71e+11 33.3% # 3.7 cycles per ns
└ task-clock 4.58e+10 33.3% # 45.8 s
┌ instructions 3.03e+11 66.7% # 1.8 insns per cycle
│ branch-instructions 6.96e+09 66.7% # 2.3% of insns
└ branch-misses 5.65e+06 66.7% # 0.1% of branch insns
┌ L1-dcache-load-misses 1.07e+09 33.3% # 1.1% of dcache loads
│ L1-dcache-loads 1.01e+11 33.3%
│ cache-misses 1.02e+09 33.3% # 95.6% of cache refs
└ cache-references 1.06e+09 33.3%
aggregated from 18 threads
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━
Look at the difference in number of instructions needed (and the number of instructions per clock cycle)!
That is, the @threads version required 38.5 times as many instructions as the @tturbo version!
This is running it 100 times.
Howe do I make sure I use the correct VectorizationBase? I just did using VectorizationBase and I got
│ Package VectorizationBase not found, but a package named VectorizationBase is available from a registry.
│ Install package?
│ (v1.7) pkg> add VectorizationBase
└ (y/n) [y]:
Resolving package versions...
Updating `/global/u1/d/dforero/.julia/nersc/perlmutter/environments/v1.7/Project.toml`
[3d5dd08c] + VectorizationBase v0.21.49
No Changes to `/global/u1/d/dforero/.julia/nersc/perlmutter/environments/v1.7/Manifest.toml`
ERROR: Unsatisfiable requirements detected for package VectorizationBase [3d5dd08c]:
VectorizationBase [3d5dd08c] log:
├─possible versions are: 0.1.0-0.21.49 or uninstalled
└─restricted to versions 0.21.54 by an explicit requirement — no versions left