Support Symmetric Matrices

glennmoy commented 4 years ago

The problem is not resolved by defining Base.convert(::Type{Symmetric}, x::Array{Float64, 2}).

I believe the problem is the non-commutativity of Symmetric and BlockArray methods. Specifically, I think the sub-block types are being defined as Symmetric in this line and then setindex! complains when it tries to assign an Array{T, 2} to the blocks.

At least that's my educated guess. Folks more familiar with what those functions do will know more.

julia> using BlockArrays, LinearAlgebra

julia> Symmetric(BlockArray(rand(5, 5), [3, 2], [2, 3]))
5×5 Symmetric{Float64,BlockArray{Float64,2,Array{Array{Float64,2},2},Tuple{BlockedUnitRange{Array{Int64,1}},BlockedUnitRange{Array{Int64,1}}}}} with indices 1:1:5×1:1:5:
 0.80611   0.980725  │  0.396233  0.218845  0.126326
 0.980725  0.519793  │  0.370195  0.88799   0.996251
 ────────────────────┼──────────────────────────────
 0.396233  0.370195  │  0.344366  0.171352  0.478706
 0.218845  0.88799   │  0.171352  0.789426  0.58824 
 0.126326  0.996251  │  0.478706  0.58824   0.565382

julia> BlockArray(Symmetric(rand(5, 5)), [3, 2], [2, 3])
ERROR: MethodError: Cannot `convert` an object of type Array{Float64,2} to an object of type Symmetric{Float64,Array{Float64,2}}
Closest candidates are:
  convert(::Type{#s664} where #s664<:Symmetric, ::Union{Hermitian, Symmetric}) at /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.4/LinearAlgebra/src/symmetric.jl:193
  convert(::Type{T}, ::T) where T<:AbstractArray at abstractarray.jl:14
  convert(::Type{T}, ::Factorization) where T<:AbstractArray at /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.4/LinearAlgebra/src/factorization.jl:55
  ...
Stacktrace:
 [1] setindex! at ./array.jl:828 [inlined]
 [2] setblock! at /Users/glenn/.julia/packages/BlockArrays/Iqf2f/src/blockarray.jl:400 [inlined]
 [3] setindex! at /Users/glenn/.julia/packages/BlockArrays/Iqf2f/src/abstractblockarray.jl:203 [inlined]
 [4] macro expansion at /Users/glenn/.julia/packages/BlockArrays/Iqf2f/src/blockarray.jl:184 [inlined]
 [5] macro expansion at ./cartesian.jl:64 [inlined]
 [6] macro expansion at /Users/glenn/.julia/packages/BlockArrays/Iqf2f/src/blockarray.jl:181 [inlined]
 [7] BlockArray{Float64,N,R,BS} where BS<:Tuple{Vararg{AbstractUnitRange{Int64},N}} where R<:(AbstractArray{#s13,N} where #s13<:AbstractArray{Float64,N}) where N(::Symmetric{Float64,Array{Float64,2}}, ::Tuple{BlockedUnitRange{Array{Int64,1}},BlockedUnitRange{Array{Int64,1}}}) at /Users/glenn/.julia/packages/BlockArrays/Iqf2f/src/blockarray.jl:179
 [8] BlockArray{Float64,N,R,BS} where BS<:Tuple{Vararg{AbstractUnitRange{Int64},N}} where R<:(AbstractArray{#s13,N} where #s13<:AbstractArray{Float64,N}) where N(::Symmetric{Float64,Array{Float64,2}}, ::Array{Int64,1}, ::Array{Int64,1}) at /Users/glenn/.julia/packages/BlockArrays/Iqf2f/src/blockarray.jl:172
 [9] BlockArray(::Symmetric{Float64,Array{Float64,2}}, ::Array{Int64,1}, ::Array{Int64,1}) at /Users/glenn/.julia/packages/BlockArrays/Iqf2f/src/blockarray.jl:175
 [10] top-level scope at REPL[199]:1
 [11] eval(::Module, ::Any) at ./boot.jl:331
 [12] eval_user_input(::Any, ::REPL.REPLBackend) at /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.4/REPL/src/REPL.jl:86
 [13] run_backend(::REPL.REPLBackend) at /Users/glenn/.julia/packages/Revise/tV8FE/src/Revise.jl:1165
 [14] top-level scope at none:0

dlfivefifty commented 4 years ago

Thanks for the issue. Don't have time to fix this right now but if you make a PR I can review it

glennmoy commented 4 years ago

No problem - I'm deep diving into the package as we speak to figure out how to fix it. I'll post a PR if I manage to but just wanted to open the issue first :)

nickrobinson251 commented 4 years ago

I think the sub-block types are being defined as Symmetric

Yeah, i think that right: the blocks are being converted to the same type T as the input array. But that doesn't make sense for "structured" array types (like Symmetric).

The errors we hit are

calling convert(T, block): for most structured array type convert(T, ::Array) is not defined (MethodError)
if it were defined (e.g. if we were instead calling T(block), it may nevertheless error: e.g. trying to convert a non-square matrix to Symmetric throws a DimensionMismatch

I guess we need a different approach for handling Structured matrix types? Store the blocks as Matrixes? But I might be missing something, as i'd have expected the working with Structured matrix types as input would have come up before (but maybe not). And i'm not too familiar with the BlockArrays 😊

dlfivefifty commented 4 years ago

It looks like the issue is that it is assuming the storage type is the same as the array type:

https://github.com/JuliaArrays/BlockArrays.jl/blob/abf38ef8a495b87f6ca7b260e636c0346dc64e35/src/blockarray.jl#L180

Probably best to default to Array so change this to something like:

...
block_arr = _BlockArray(Array{_block_typeof(arr),N}, baxes)
...
_block_typeof(::AbstractArray{T,N}) where {T,N} = Array{T,N}
_block_typeof(A::AbstractSparseArray) = typeof(A)

JuliaArrays / BlockArrays.jl

Support Symmetric Matrices #128