Strange behavior of BlockDiagonal constructor with Vector of BandedMatrices

As you can see from the MWE I have to call D = BlockDiagonal([Dbv[1], Dbv[2], Dbv[3]]) instead of D = BlockDiagonal(Dbv) . It seems that some type information is lost when the Vector is banded matrices is created, but then it is recovered by calling the constructor with [Dbv[1], Dbv[2], Dbv[3]] …

using BandedMatrices
using BlockDiagonals

ix = [9, 16, 22]
Dbv = []
iold = 0
for i in ix
    Db = BandedMatrix{Float64}(undef, (i-iold-1, i-iold), (0,1))
    Db[band(0)]  .= -1
    Db[band(1)] .=  1
    iold = i
    push!(Dbv, copy(Db))
end
Dbv[1]

8Γ—9 BandedMatrix{Float64} with bandwidths (0, 1):
 -1.0   1.0    β‹…     β‹…     β‹…     β‹…     β‹…     β‹…    β‹… 
   β‹…   -1.0   1.0    β‹…     β‹…     β‹…     β‹…     β‹…    β‹… 
   β‹…     β‹…   -1.0   1.0    β‹…     β‹…     β‹…     β‹…    β‹… 
   β‹…     β‹…     β‹…   -1.0   1.0    β‹…     β‹…     β‹…    β‹… 
   β‹…     β‹…     β‹…     β‹…   -1.0   1.0    β‹…     β‹…    β‹… 
   β‹…     β‹…     β‹…     β‹…     β‹…   -1.0   1.0    β‹…    β‹… 
   β‹…     β‹…     β‹…     β‹…     β‹…     β‹…   -1.0   1.0   β‹… 
   β‹…     β‹…     β‹…     β‹…     β‹…     β‹…     β‹…   -1.0  1.0
D = BlockDiagonal(Dbv)

MethodError: no method matching BlockDiagonal(::Vector{Any})
Closest candidates are:
  BlockDiagonal(::BlockDiagonal) at ~/.julia/packages/BlockDiagonals/Y9q2S/src/blockdiagonal.jl:20
  BlockDiagonal(::Vector{V}) where {T, V<:AbstractMatrix{T}} at ~/.julia/packages/BlockDiagonals/Y9q2S/src/blockdiagonal.jl:16

Stacktrace:
 [1] top-level scope
   @ In[3]:1
 [2] eval
   @ ./boot.jl:373 [inlined]
 [3] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)
   @ Base ./loading.jl:1196
D = BlockDiagonal([Dbv[1], Dbv[2], Dbv[3]])

19Γ—22 BlockDiagonal{Float64, BandedMatrix{Float64, Matrix{Float64}, Base.OneTo{Int64}}}:
 -1.0   1.0   0.0   0.0   0.0   0.0  …   0.0   0.0   0.0   0.0   0.0  0.0
  0.0  -1.0   1.0   0.0   0.0   0.0      0.0   0.0   0.0   0.0   0.0  0.0
  0.0   0.0  -1.0   1.0   0.0   0.0      0.0   0.0   0.0   0.0   0.0  0.0
  0.0   0.0   0.0  -1.0   1.0   0.0      0.0   0.0   0.0   0.0   0.0  0.0
  0.0   0.0   0.0   0.0  -1.0   1.0      0.0   0.0   0.0   0.0   0.0  0.0
  0.0   0.0   0.0   0.0   0.0  -1.0  …   0.0   0.0   0.0   0.0   0.0  0.0
  0.0   0.0   0.0   0.0   0.0   0.0      0.0   0.0   0.0   0.0   0.0  0.0
  0.0   0.0   0.0   0.0   0.0   0.0      0.0   0.0   0.0   0.0   0.0  0.0
  0.0   0.0   0.0   0.0   0.0   0.0      0.0   0.0   0.0   0.0   0.0  0.0
  0.0   0.0   0.0   0.0   0.0   0.0      0.0   0.0   0.0   0.0   0.0  0.0
  0.0   0.0   0.0   0.0   0.0   0.0  …   0.0   0.0   0.0   0.0   0.0  0.0
  0.0   0.0   0.0   0.0   0.0   0.0      0.0   0.0   0.0   0.0   0.0  0.0
  0.0   0.0   0.0   0.0   0.0   0.0      0.0   0.0   0.0   0.0   0.0  0.0
  0.0   0.0   0.0   0.0   0.0   0.0      0.0   0.0   0.0   0.0   0.0  0.0
  0.0   0.0   0.0   0.0   0.0   0.0     -1.0   1.0   0.0   0.0   0.0  0.0
  0.0   0.0   0.0   0.0   0.0   0.0  …   0.0  -1.0   1.0   0.0   0.0  0.0
  0.0   0.0   0.0   0.0   0.0   0.0      0.0   0.0  -1.0   1.0   0.0  0.0
  0.0   0.0   0.0   0.0   0.0   0.0      0.0   0.0   0.0  -1.0   1.0  0.0
  0.0   0.0   0.0   0.0   0.0   0.0      0.0   0.0   0.0   0.0  -1.0  1.0

The type information isn’t lost, you explicitly removed it. When you do [], you create a Vector{Any}.

EDIT: actually instead of β€œremoved” I should have said you explicitly asked Julia not to create it in the first place.

Dbv = [] creates a Vector{Any}. You want Dbv = BandedMatrix[]

1 Like

Got it, thanks!