Surprising denseness?

[jagot]

I just converted my code to use BlockBandedMatrices (the latest release) instead of my homegrown BlockMaps and noticed a substantial performance degradation. I found out that although my operators where block-banded, anytime I used them, they had become dense (due to multiplications by scalars, etc). In summary, these operations, for T::BlockBandedMatrix, yield dense matrices:

• T/2
• 2T (also 2 .* T)
• T *= 2
• T /= 2
• T + I
• copy(T)
• similar(T)

All of the above work as I would expect for e.g. Tridiagonal matrices.

T ./= 2 does not become dense, neither does lmul!(2, T), but sometimes I would like to have a scaled copy of a matrix, without modifying the original one. Is this all by design?

[dlfivefifty]

This should be an easy fix: all these operations are fast for BandedBlockBanded via a fast broadcast mechanism that reduces to the underlying data.

[dlfivefifty]

Are you on master:

julia> A = BandedBlockBandedMatrix{Float64}(undef, (1:5,1:5), (1,1), (1,1));

julia> typeof(2A)
BandedBlockBandedMatrix{Float64}

julia> A = BlockBandedMatrix{Float64}(undef, (1:5,1:5), (1,1));

julia> typeof(2A)
BlockBandedMatrix{Float64}

[jagot]

Alas, I am not, due to the stack-overflow in matrix–vector multiplications. Can confirm that all the operations above work on master though.

[dlfivefifty]

Unfortunately I forgot to push a change at home to start tagging everything. I'll try to get everything updated (including your PR) over the weekend. Hope you can wait!

[jagot]

I have since discovered that the stack overflow only occurs if either the matrix or vector has eltype Int (I have only checked Float64 and Int64, and combinations thereof).

[dlfivefifty]

Ah OK. If you give a MWE I can have a look. (The new LazyArrays.jl is a bit better at only dispatching to BLAS when all types are correct, though this needs more testing.)