SciML
`mul!(::Matrix, ::ArrayPartition, ::Adjoint(ArrayPartition))`
https://github.com/JuliaNLSolvers/Optim.jl/pull/972
MWE
using RecursiveArrayTools, LinearAlgebra
u = ArrayPartition(rand(2), rand(1))
v = copy(u)
H = u * v' # works
mul!(H, u, v') # error
julia> mul!(H,u,v')
ERROR: BoundsError: attempt to access Tuple{Vector{Float64}, Vector{Float64}} at index [3]
Stacktrace:
[1] getindex
@ ~/.julia/dev/RecursiveArrayTools.jl/src/array_partition.jl:151 [inlined]
[2] copy_transpose!(B::Matrix{Float64}, ir_dest::UnitRange{Int64}, jr_dest::UnitRange{Int64}, A::ArrayPartition{Float64, Tuple{Vector{Float64}, Vector{Float64}}}, ir_src::UnitRange{Int64}, jr_src::UnitRange{Int64})
@ LinearAlgebra /Applications/Julia-1.7.app/Contents/Resources/julia/share/julia/stdlib/v1.7/LinearAlgebra/src/transpose.jl:197
[3] copy_transpose!
@ /Applications/Julia-1.7.app/Contents/Resources/julia/share/julia/stdlib/v1.7/LinearAlgebra/src/matmul.jl:693 [inlined]
[4] _generic_matmatmul!(C::Matrix{Float64}, tA::Char, tB::Char, A::ArrayPartition{Float64, Tuple{Vector{Float64}, Vector{Float64}}}, B::ArrayPartition{Float64, Tuple{Vector{Float64}, Vector{Float64}}}, _add::LinearAlgebra.MulAddMul{true,
true, Bool, Bool})
@ LinearAlgebra /Applications/Julia-1.7.app/Contents/Resources/julia/share/julia/stdlib/v1.7/LinearAlgebra/src/matmul.jl:834
[5] generic_matmatmul!
@ /Applications/Julia-1.7.app/Contents/Resources/julia/share/julia/stdlib/v1.7/LinearAlgebra/src/matmul.jl:801 [inlined]
[6] mul!
@ /Applications/Julia-1.7.app/Contents/Resources/julia/share/julia/stdlib/v1.7/LinearAlgebra/src/matmul.jl:478 [inlined]
[7] mul!(C::Matrix{Float64}, A::ArrayPartition{Float64, Tuple{Vector{Float64}, Vector{Float64}}}, B::Adjoint{Float64, ArrayPartition{Float64, Tuple{Vector{Float64}, Vector{Float64}}}})
@ LinearAlgebra /Applications/Julia-1.7.app/Contents/Resources/julia/share/julia/stdlib/v1.7/LinearAlgebra/src/matmul.jl:275
[8] top-level scope
@ REPL[134]:1
in LinearAlgebra/transpose.jl, v is being treated as an AbstractVecOrMat and is being accessed as follows v[3,1]
julia> checkbounds(Bool,v,3,1)
true
julia> v[3,1]
ERROR: BoundsError: attempt to access Tuple{Vector{Float64}, Vector{Float64}} at index [3]
Stacktrace:
[1] getindex(A::ArrayPartition{Float64, Tuple{Vector{Float64}, Vector{Float64}}}, i::Int64, j::Int64)
@ RecursiveArrayTools ~/.julia/dev/RecursiveArrayTools.jl/src/array_partition.jl:151
it would be good to define a checkbounds(::ArrayPartition)
function Base.checkbounds(A::ArrayPartition, i, j...)
@boundscheck checkbounds(A.x, i)
b = A.x[i]
@boundscheck checkbounds(b, j)
end
though that still doesn't solve my problem
No, it's a vector not a matrix so the indexing A[i,j] does not make sense. There's not a clear matrix definition because there's no constraint that the sizes of each underlying array align.
i understand that ArrayPartitions are vectors and A[i,j] isn't the same as standard cartesian indexing.
however in https://github.com/JuliaLang/julia/blob/master/stdlib/LinearAlgebra/src/transpose.jl#L181-L203,
A::ArrayPartition{Float64, Tuple{Vector{Float64}, Vector{Float64}}} is treated as an AbstractVecOrMat and is accessed like A[isrc,jsrc] == A[3,1] in line 197 with boundscheck @boundscheck checkbounds(A, ir_src, jr_src) passing in line 192. (see stack trace above)
IDK how to avoid this behaviour
I think mul! and such need an overload to handle this well. But honestly, if you're using mul!, you should probably use ComponentArrays.jl instead.