`lu!` does not support `SubArray` indexed by `Vector` rather `UnitRange`

[wsshin]

lu! works for A[ii,jj] when ii and jj are UnitRange{Int64}:

julia> VERSION
v"0.7.1-pre.0"

julia> A = rand(3,3)
3×3 Array{Float64,2}:
 0.433089  0.964333  0.921431
 0.736251  0.108516  0.147767
 0.435144  0.707014  0.36105

julia> Asub1 = @view(A[2:3,2:3])
2×2 view(::Array{Float64,2}, 2:3, 2:3) with eltype Float64:
 0.108516  0.147767
 0.707014  0.36105

julia> lu!(Asub1)
LU{Float64,SubArray{Float64,2,Array{Float64,2},Tuple{UnitRange{Int64},UnitRange{Int64}},false}}
L factor:
2×2 Array{Float64,2}:
 1.0       0.0
 0.153484  1.0
U factor:
2×2 Array{Float64,2}:
 0.707014  0.36105
 0.0       0.092351

However, lu! does not work for A[ii,jj] if ii and jj are Vector{Int64}:

julia> Asub2 = @view(A[[1,3],[1,3]])
2×2 view(::Array{Float64,2}, [1, 3], [1, 3]) with eltype Float64:
 0.433089  0.921431
 0.435144  0.092351

julia> lu!(Asub2)
ERROR: MethodError: no method matching lu!(::SubArray{Float64,2,Array{Float64,2},Tuple{Array{Int64,1},Array{Int64,1}},false})
Closest candidates are:
  lu!(::Union{DenseArray{T<:Union{Complex{Float32}, Complex{Float64}, Float32, Float64},2}, ReinterpretArray{T<:Union{Complex{Float32}, Complex{Float64}, Float32, Float64},2,S,A} where S where A<:Union{SubArray{T,N,A,I,true} where I<:Tuple{AbstractUnitRange,Vararg{Any,N} where N} where A<:DenseArray where N where T, DenseArray}, ReshapedArray{T<:Union{Complex{Float32}, Complex{Float64}, Float32, Float64},2,A,MI} where MI<:Tuple{Vararg{SignedMultiplicativeInverse{Int64},N} where N} where A<:Union{ReinterpretArray{T,N,S,A} where S where A<:Union{SubArray{T,N,A,I,true} where I<:Tuple{AbstractUnitRange,Vararg{Any,N} where N} where A<:DenseArray where N where T, DenseArray} where N where T, SubArray{T,N,A,I,true} where I<:Tuple{AbstractUnitRange,Vararg{Any,N} where N} where A<:DenseArray where N where T, DenseArray}, SubArray{T<:Union{Complex{Float32}, Complex{Float64}, Float32, Float64},2,A,I,L} where L where I<:Tuple{Vararg{Union{Int64, AbstractRange{Int64}, AbstractCartesianIndex},N} where N} where A<:Union{ReinterpretArray{T,N,S,A} where S where A<:Union{SubArray{T,N,A,I,true} where I<:Tuple{AbstractUnitRange,Vararg{Any,N} where N} where A<:DenseArray where N where T, DenseArray} where N where T, ReshapedArray{T,N,A,MI} where MI<:Tuple{Vararg{SignedMultiplicativeInverse{Int64},N} where N} where A<:Union{ReinterpretArray{T,N,S,A} where S where A<:Union{SubArray{T,N,A,I,true} where I<:Tuple{AbstractUnitRange,Vararg{Any,N} where N} where A<:DenseArray where N where T, DenseArray} where N where T, SubArray{T,N,A,I,true} where I<:Tuple{AbstractUnitRange,Vararg{Any,N} where N} where A<:DenseArray where N where T, DenseArray} where N where T, DenseArray}}) where T<:Union{Complex{Float32}, Complex{Float64}, Float32, Float64} at /Users/wsshin/packages/julia/julia-0.7/usr/share/julia/stdlib/v0.7/LinearAlgebra/src/lu.jl:37
  lu!(::Union{DenseArray{T<:Union{Complex{Float32}, Complex{Float64}, Float32, Float64},2}, ReinterpretArray{T<:Union{Complex{Float32}, Complex{Float64}, Float32, Float64},2,S,A} where S where A<:Union{SubArray{T,N,A,I,true} where I<:Tuple{AbstractUnitRange,Vararg{Any,N} where N} where A<:DenseArray where N where T, DenseArray}, ReshapedArray{T<:Union{Complex{Float32}, Complex{Float64}, Float32, Float64},2,A,MI} where MI<:Tuple{Vararg{SignedMultiplicativeInverse{Int64},N} where N} where A<:Union{ReinterpretArray{T,N,S,A} where S where A<:Union{SubArray{T,N,A,I,true} where I<:Tuple{AbstractUnitRange,Vararg{Any,N} where N} where A<:DenseArray where N where T, DenseArray} where N where T, SubArray{T,N,A,I,true} where I<:Tuple{AbstractUnitRange,Vararg{Any,N} where N} where A<:DenseArray where N where T, DenseArray}, SubArray{T<:Union{Complex{Float32}, Complex{Float64}, Float32, Float64},2,A,I,L} where L where I<:Tuple{Vararg{Union{Int64, AbstractRange{Int64}, AbstractCartesianIndex},N} where N} where A<:Union{ReinterpretArray{T,N,S,A} where S where A<:Union{SubArray{T,N,A,I,true} where I<:Tuple{AbstractUnitRange,Vararg{Any,N} where N} where A<:DenseArray where N where T, DenseArray} where N where T, ReshapedArray{T,N,A,MI} where MI<:Tuple{Vararg{SignedMultiplicativeInverse{Int64},N} where N} where A<:Union{ReinterpretArray{T,N,S,A} where S where A<:Union{SubArray{T,N,A,I,true} where I<:Tuple{AbstractUnitRange,Vararg{Any,N} where N} where A<:DenseArray where N where T, DenseArray} where N where T, SubArray{T,N,A,I,true} where I<:Tuple{AbstractUnitRange,Vararg{Any,N} where N} where A<:DenseArray where N where T, DenseArray} where N where T, DenseArray}}, ::Union{Val{false}, Val{true}}; check) where T<:Union{Complex{Float32}, Complex{Float64}, Float32, Float64} at /Users/wsshin/packages/julia/julia-0.7/usr/share/julia/stdlib/v0.7/LinearAlgebra/src/lu.jl:37
  lu!(::Union{Hermitian{T,S}, Symmetric{T,S}} where S where T) at /Users/wsshin/packages/julia/julia-0.7/usr/share/julia/stdlib/v0.7/LinearAlgebra/src/lu.jl:45
  ...
Stacktrace:
 [1] top-level scope at none:0

Still, lu (instead of lu!) works even in that case:

julia> lu(Asub2)
LU{Float64,Array{Float64,2}}
L factor:
2×2 Array{Float64,2}:
 1.0       0.0
 0.995276  1.0
U factor:
2×2 Array{Float64,2}:
 0.435144  0.092351
 0.0       0.829516

[andreasnoack]

We could probably make this work for the general fallback. It's currently restricted to StridedArray. However, I suspect it will typically be faster to copy. What is the application?

[wsshin]

I need to set certain rows and columns of A indicated by ind and then solve A[ind,ind] \ b[ind]. I need to do this many times (> a thousand). Every time I solve the equation, ind changes. I wanted to reduce the number of allocations by

ind = # some code
for i = ind, j = ind
    A[i,j] = # some code
end

LU = lu!(@view(A[ind,ind]))
ldiv(x[ind], LU, @view(b[ind]))

instead of

ind = # some code
for i = ind, j = ind
    A[i,j] = # some code
end

x = @view(A[ind,ind]) \ @view(b[ind])

[wsshin]

Further performance test to justify the needs for lu!(@view(A[ind,ind])). I combined (lu, lu!) with (A[ind,ind], @view(A[ind,ind])), and here are the results:

julia> A = rand(40,40);

julia> @btime lu($A[1:10,1:10]);
  1.122 μs (4 allocations: 1.94 KiB)

julia> @btime lu($@view(A[1:10,1:10]));
  1.189 μs (3 allocations: 1.06 KiB)

julia> @btime lu!($A[1:10,1:10]);
  1.119 μs (3 allocations: 1.06 KiB)

julia> @btime lu!($@view(A[1:10,1:10]));
  839.973 ns (2 allocations: 192 bytes)

[timholy]

The error should be fixed, but @andreasnoack 's point is that once you are using non-ranges the fastest implementation will almost certainly be lu!(A[ind,ind]). LU decomposition is O(N^3) and copying is O(N^2), so the copy should be insignificant in terms of the overall cost. In contrast, slowing down the element access---which is inevitable with the view since you require an indirect lookup in ind---would be a performance hit you'd feel, no matter how smart your implementation.

If you're using very small matrices, you should consider StaticArrays which may allow you to avoid stack allocation altogether.