Specialized `default_blasmul!` for `Diagonal`, `Tridiagonal`, and `SymTridiagonal`?

[jagot]

Since

julia> Base.IndexStyle(Diagonal)
IndexCartesian()

julia> Base.IndexStyle(Tridiagonal)
IndexCartesian()

julia> Base.IndexStyle(SymTridiagonal)
IndexCartesian()

default_blasmul! will fallback to a dense algorithm: https://github.com/JuliaMatrices/ArrayLayouts.jl/blob/a5767f77a867758956c58078442bc61ca8c345c1/src/muladd.jl#L202-L221
which is suboptimal.

Would we like something like the following in ArrayLayouts.jl?:

for MT in [:Diagonal, :SymTridiagonal, :Tridiagonal]
    @eval begin
        ArrayLayouts.default_blasmul!(α, A::$MT, B::AbstractVector, β, C::AbstractVector) =
            mul!(C, A, B, α, β)
    end
end

[dlfivefifty]

That's not true: _default_blasmul! already takes advantage of sparsity via colsupport.