use LinearAlgebra's new generic in-place 5-arg mul! (requires Julia 1.3)?

[chriscoey]

see JuliaLang/julia#29634 - this may be relevant for the mul! definition at https://github.com/Jutho/LinearMaps.jl/blob/a024a4dd26cd4752e43142feda5cf3da8b5dc601/src/LinearMaps.jl#L23 and maybe some of the definitions at https://github.com/Jutho/LinearMaps.jl/blob/a024a4dd26cd4752e43142feda5cf3da8b5dc601/src/wrappedmap.jl#L27

ultimately I'm interested in using LinearMaps.jl with IterativeSolvers.jl in a way that avoids unnecessary allocations and dispatches to the fastest methods possible, for any linear maps that I construct.

[Jutho]

Yes indeed, we will want to overwrite the fallback method for 5-arg mul! of linear maps for the WrappedMap{T,<:AbstractMatrix} case. Currently only UniformScalingMap has a specialized 5-arg mul! implementation.

[dkarrasch]

We may actually want to consider making this „type“ invariant under some operations:

const MatrixMap{T} = WrappedMap{T,<:AbstractMatrix}
transpose(A::MatrixMap) = LinearMap(transpose(A.lmap))

because transpose is as lazy in Base as is our TransposeMap, and same with adjoint... and potentially +, - and scalar multiplication, so relatively cheap matrix operations. Cf. JuliaLang/julia#21, which had a different motivation originally, though.

[andreasnoack]

Notice that the tests currently fail on the release branch of Julia 1.3 because of the direct use of generic_mul!

[dkarrasch]

Thanks, @andreasnoack, I know. We have JuliaLang/julia#62 to restrict our current version to v1.2. Maybe you could provide some advice on what you think is reasonable there? Has there emerged a "standard" way of handling the v1.0 LTS vs. v1.3 generic_mul! issue in linear algebra related packages?

[andreasnoack]

I've made this solution https://github.com/JuliaLinearAlgebra/GenericLinearAlgebra.jl/blob/master/src/juliaBLAS.jl#L87 in GenericLinearAlgebra and something like that might work here as well. It's quite simple and the only downside seems to be that it's bad for your coverage stats.

[tkf]

Is this the only place you need to use generic_mul!?

https://github.com/Jutho/LinearMaps.jl/blob/1555f5658cee7df2ad84aecf65acaa07fd4a5156/src/uniformscalingmap.jl#L27-L36

Then why not use broadcasting? As you've already found out, it's faster (JuliaLang/LinearAlgebra.jl#661) even after a workaround JuliaLang/julia#33108.

[dkarrasch]

This is the pre-1.3 version, which works fine. For the 1.3 version, I do use broadcasting in JuliaLang/julia#62. BTW, is broadcasting also faster for larger vectors/matrices? If that was the case, why not use it in mul! along the lines of
https://github.com/Jutho/LinearMaps.jl/blob/1555f5658cee7df2ad84aecf65acaa07fd4a5156/src/uniformscalingmap.jl#L40-L70
where the 0-1-logic is probably handled by the MulAddMul-thing.

[tkf]

is broadcasting also faster for larger vectors/matrices?

It looks like so. I think it makes sense to use broadcasting. I tried to do a minimal possible change in 5-arg mul! and there is no other reasons why I didn't do it. The diagonal matrices are actually handled via broadcasting. I guess that's probably because mul! by a number predates broadcasting or something?

Just FYI, if you want to avoid writing all the if branches, look at how mul!(out::AbstractVector, A::Diagonal, in::AbstractVector) etc. are implemented. I found that it was a neat trick :)
https://github.com/JuliaLang/julia/blob/b5f4e877743f484de43dc4a8d30f74499ec96016/stdlib/LinearAlgebra/src/diagonal.jl#L279-L290

[tkf]

...although it is probably unnecessary if you just want to define one function