Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

mul!(::Matrix, ::ArrayPartition, ::Adjoint(ArrayPartition)) #186

Open
vpuri3 opened this issue Jan 31, 2022 · 4 comments
Open

mul!(::Matrix, ::ArrayPartition, ::Adjoint(ArrayPartition)) #186

vpuri3 opened this issue Jan 31, 2022 · 4 comments

Comments

@vpuri3
Copy link
Member

vpuri3 commented Jan 31, 2022

JuliaNLSolvers/Optim.jl#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
@vpuri3
Copy link
Member Author

vpuri3 commented Jan 31, 2022

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

@ChrisRackauckas
Copy link
Member

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.

@vpuri3
Copy link
Member Author

vpuri3 commented Feb 1, 2022

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

@ChrisRackauckas
Copy link
Member

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.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

2 participants