Add support for kron

[devon-research]

Do people think this would be appropriate / useful to include?

And do people have ideas for how to make it efficient?

Here is a rudimentary and inefficient but seemingly functional implementation I used to get something working:

function Base.kron(a::BandedMatrix{T}, b::BandedMatrix{T}) where T
    R = BandedMatrix{T}(undef, (size(a,1) * size(b,1), size(a,2) * size(b,2)), (bandwidth(a,1) * size(b,1) + bandwidth(b,1), bandwidth(a,2) * size(b,2) + bandwidth(b,2)))
    m = 0
    @inbounds for j = 1:size(a,2), l = 1:size(b,2), i = 1:size(a,1)
        aij = a[i,j]
        for k = 1:size(b,1)
            R[m += 1] = aij * b[k,l]
        end
    end
    R
end

[dlfivefifty]

Yep! Note it currently exists in LazyBandedMatrices.jl, but would be nice here.

Ideally it would support an arbitrary number of matrices

[dlfivefifty]

Note quick and dirty for more matrices:

kron(A::BandedMatrix, B::BandedMatrix, C::AbstractMatrix...) = kron(kron(A,B), C...)

[dlfivefifty]

Also, no need to restrict to same type, do:

function Base.kron(a::BandedMatrix{T}, b::BandedMatrix{V}) where {T,V}
   R = BandedMatrix{Base.promote_op(*,T,V)}(undef, (size(a,1) * size(b,1), size(a,2) * size(b,2)), (bandwidth(a,1) * size(b,1) + bandwidth(b,1), bandwidth(a,2) * size(b,2) + bandwidth(b,2)))
   ...
end