Speed of multiplication of Diagonal

[vincentcp]

The Mul of two Diagonals might create a full matrix somewhere when applied to a vector.

julia> using LinearAlgebra, LazyArrays
julia> d = Diagonal(rand(1000));
julia> b = rand(1000);
julia> d*d*b;@time d*d*b;
  0.000004 seconds (7 allocations: 16.047 KiB)
julia> M = Mul(d,d);
julia> M*b;@time M*b;
  0.691181 seconds (18 allocations: 7.638 MiB, 0.84% gc time)
julia> LazyArrays.materialize(Mul(M.factors[1], M.factors[2])) isa Array
true

[vincentcp]

Could

LazyArrays.jl/src/linalg/blasmul.jl

Line 145 in 4559ba2

function default_blasmul!(α, A::AbstractMatrix, B::AbstractMatrix, β, C::AbstractMatrix)

and

LazyArrays.jl/src/linalg/blasmul.jl

Line 162 in 4559ba2

function default_blasmul!(α, A::AbstractMatrix, B::AbstractVector, β, C::AbstractVector)

use mul! instead?

Or should materialize! be implemented for Diagonal in order to get an efficient materialize?

[dlfivefifty]

No. The whole point is to support the 5-argument version (αAx + β*y) which mul! does not.

[vincentcp]

Ok, I understand. Which functions need to be specialised for Diagonal/DiagonalLayout to get the LazyArrays machinery working?

[dlfivefifty]

It loosely follows the broadcast interface, so you need to

1. Overload similar
2. Overload _copyto!

Banded matrices give a nice example: https://github.com/JuliaMatrices/BandedMatrices.jl/blob/master/src/generic/matmul.jl

Note that many cases first lower to the 5-argument version α*A*x + β*y, which is implemented via MulAdd. For these you overload materialize! and overwrite the y argument with the value.

[dlfivefifty]

This is closer to being fixed, though now the issue is that its creating a (slow) SparseMatrixCSC:

julia> apply(*,d,d) |> typeof
SparseMatrixCSC{Float64,Int64}

And we still neeed a special materialize! for diagonal*diagonal.

[dlfivefifty]

This has probably been fixed and is now in ArrayLayouts.jl