Hi,
I have implemented a linear operator as a method of LinearAlgebra.mul!(x, A, b) to use with IterativeSolvers.jl, but when profiling my code I noticed a huge performance difference due to unaliasing arrays when testing my operator on its own versus in a GMRES call. I realized that this is because my mul! method reshapes x and b and then does broadcasting with them, which matters because in GMRES x and b are views of columns of the same array whereas in my testing I had allocated x and b in different arrays. So my problem is that may reshaped views of non-overlapping memory may alias, shown here:
N=3
a = rand(N*N,2)
a1 = view(a,:,1)
a2 = view(a,:,2)
Base.mightalias(a1, a2) # false
aa1 = reshape(a1, N, N)
aa2 = reshape(a2, N, N)
Base.mightalias(aa1, aa2) # true
So I would like to ask if there is any way to tell the broadcasting machinery that these reshaped views wouldn’t alias.
Thanks,
Lorenzo
P.S. Since I want the most performance possible for my linear operator, I will replace my problematic broadcast expressions with hand-written loops that improve memory locality, but I am still curious if there is a possibility of a more precise Base.mightalias method that can improve the performance of my naive broadcast expressions.
I might be wrong, but my intuition is that AbstractArray wrappers of non-aliasing arrays wouldn’t alias either, and mightalias just doesn’t know how to check.
julia> @which Base.mightalias(a1, a2)
mightalias(A::SubArray, B::SubArray)
@ Base multidimensional.jl:1049
julia> @which Base.mightalias(aa1, aa2)
mightalias(A::AbstractArray, B::AbstractArray)
@ Base abstractarray.jl:1537
julia> methods(Base.mightalias)
# 3 methods for generic function "mightalias" from Base:
[1] mightalias(A::SubArray, B::SubArray)
@ multidimensional.jl:1049
[2] mightalias(A::AbstractArray, B::AbstractArray)
@ abstractarray.jl:1537
[3] mightalias(x, y)
@ abstractarray.jl:1538
The ::Any fallback assumes no aliasing, but we can ignore that danger because we’re dealing with AbstractArrays. The AbstractArray fallback guesses no aliasing if:
either array isbits, so they have no references to share data
either array isempty, so they have no data to share at the moment. Odd to me because that can change real fast: a = []; println(Base.mightalias(a, a)); push!(a, 1); println(Base.mightalias(a, a))
the arrays don’t reference the same memory regions, as reported by Base.dataids
One possible workaround is to replace broadcast expressions with map!, which doesn’t check for aliasing. Schematically, replace
y .= f.(x)
with
map!(f, y, x)
This is less flexible in that it doesn’t expand singleton dimensions like broadcasting does, but if it fits your use case it can help you save on the explicit loops.
I agree that the logic in the generic mightalias method doesn’t know how to unwrap the ReshapedArray wrapper to apply the method for SubArray. Perhaps the logic can be improved like this?
mightalias(A::AbstractArray, B::AbstractArray) = (parent(A) === A || A isa SubArray) && (parent(B) === B || B isa SubArray) ? !isbits(A) && !isbits(B) && !isempty(A) && !isempty(B) && !_isdisjoint(dataids(A), dataids(B)) : mightalias(A isa SubArray ? A : parent(A), B isa SubArray ? B : parent(B))
My only second thoughts about this are that it seems difficult to extend to other methods. I haven’t tried doing this with multiple dispatch.
I like the suggestion for map! since that gives the right semantics, so I’ve marked it as a solution. In my case I ended up just writing loops as it gave clearer code and since it was better for memory locality since I was originally breaking up many 2x2 matrix multiplications into a broadcast expression.