speed up multiplication with BlockArrays

[dkarrasch]

I had this little dirty piece of code, and thought I'll leave it here for when there is some common ground between BlockArrays.jl and BlockDiagonals.jl.

using LinearAlgebra, BlockArrays, BlockDiagonals
Base.@propagate_inbounds function LinearAlgebra.mul!(y::BlockMatrix, A::BlockDiagonal, x::BlockMatrix,
					α::Number=true, β::Number=false)
	yblocks = y.blocks
	xblocks = x.blocks
	blockmaps = BlockDiagonals.blocks(A)
	ym, yn = size(yblocks)
	xm, xn = size(xblocks)
	Am = length(blockmaps)
	# second condition could probably be relaxed
	(ym == Am == xm  && yn == xn == 1) || error("blocks mismatching")
	for (yi, Ai, xi) in zip(yblocks, blockmaps, xblocks)
		# some kind of elegant size checking
		(size(yi, 1) == size(Ai, 1) && size(xi, 1) == size(Ai, 2) &&
			size(yi, 2) == size(xi, 2)) || throw(DimensionMismatch())
	end
	for (yi, Ai, xi) in zip(yblocks, blockmaps, xblocks)
		# this currently only works if none of the blocks is a BlockDiagonal matrix
		mul!(yi, Ai, xi, α, β)
	end
	return y
end

Then

nblocks = 20
As = [rand(50, 50) for _ in 1:nblocks]
B = BlockDiagonal(As)
Xb = BlockArray(rand(50nblocks, 20), fill(50, nblocks), [20])
Yb = BlockArray(zeros(50nblocks, 20), fill(50, nblocks), [20])
mul!(Yb, B, Xb)
using BenchmarkTools
@btime mul!($Yb, $B, $Xb) # 80.115 μs (0 allocations: 0 bytes)
@btime $B * $(Matrix(Xb)) # 114.894 μs (162 allocations: 321.41 KiB)

is probably nice to have.

[nickrobinson251]

Thanks! I'm traveling right now but will take a proper look when I can. Perhaps related to #29 and of course #19

[nickrobinson251]

yeah, since you were using this seems to motivate #29 (which would be a prerequiste for including this, to avoid the BlockArrays.jl dependency)

[nickrobinson251]

related #37