Get rid of A*_mul_B!

[dkarrasch]

But this time it works! I'll comment on the challenging parts directly in the code. This is BREAKING, because anybody with custom LinearMap types will have the action implemented in A[*]_mul_B!s, which are no longer called from the code base here. So, after this, people would need to write mul! methods instead, which I believe would be easier to communicate.

[dkarrasch]

Here are some more specific comments for anybody who wishes to review this. I'm still amazed how all this turned out to work without changing the tests (except for removing the stupid "4-arg mul! test, which is inconsistent with LinearAlgebra anyway).

[codecov]

Codecov Report

Merging #97 into master will increase coverage by 1.03%.
The diff coverage is 100.00%.

@@            Coverage Diff             @@
##           master      #97      +/-   ##
==========================================
+ Coverage   98.29%   99.33%   +1.03%     
==========================================
  Files          12       12              
  Lines         822      904      +82     
==========================================
+ Hits          808      898      +90     
+ Misses         14        6       -8     

Impacted Files	Coverage Δ
src/LinearMaps.jl	98.76% <100.00%> (+0.25%)	⬆️
src/blockmap.jl	100.00% <100.00%> (+0.96%)	⬆️
src/composition.jl	100.00% <100.00%> (ø)
src/conversion.jl	100.00% <100.00%> (ø)
src/functionmap.jl	97.87% <100.00%> (+0.11%)	⬆️
src/kronecker.jl	100.00% <100.00%> (ø)
src/left.jl	100.00% <100.00%> (ø)
src/linearcombination.jl	100.00% <100.00%> (ø)
src/scaledmap.jl	90.00% <100.00%> (+10.00%)	⬆️
src/transpose.jl	100.00% <100.00%> (+7.14%)	⬆️
... and 13 more

---

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update 3c2e531...b1ec064. Read the comment docs.

[dkarrasch]

It may look like this is a major overhaul, but actually this only does a few distinct things, so instead of asking for a line-by-line code review, let me just tell you what is done here. Let me know if you disagree with something fundamentally, and I'll take care of adjusting the code. Besides the goal stated in the OP, I've also included the following things:

• allow a new signature mul!(::AbstractVecOrMat,::LinearMap,::AbstractVetor,[alpha,beta]), consistently with LinearAlgebra;
• move dimension checking to the top and ignoring dim-checking subsequently; in all LinearMap subtypes we do check for dimension consistency at construction, so applying the individual maps does not require additional size checking;
• import LinearAlgebra: mul! and remove the specifier in the code; this added significantly to the line length, and sometimes I forgot the specifier and suddenly defined "a new" mul! method that didn't extend LinearAlgebras. Surely, we do want to extend LinearAlgebra's mul! and nothing else;
• import Base: @propagate_inbounds to keep the code a bit shorter.

[dkarrasch]

Ok, I think I'm done with this one. I'll leave it hanging for a little while, in case somebody wants to take a closer look. Passing tests and coverage at 100% - eps() look great, don't they? 😉

[Jutho]

I'll try to give this a good look asap. Thanks for all this amazing work.

[Jutho]

Looking great already. I've left a few comments, not all of them have to do with the current changes, some are also about older changes. I have not been giving this package sufficient attention :-).

The major comment is all the @propagate_inbounds machinery, which I think is a bit abused for a purpose for which it was not intended. Typically, the dimension check does not matter in comparison to the actual multiplication, which is different for indexing, as I wrote in one of the comments.

However, because of the different nested levels of multiplication that a typical LinearMap might have, I do understand that you want to disable checking proper dimensions at all levels over and over again. How about a design where there is a single entry point mul!(y::AbstractVecOrMat, L::LinearMap, x::AbstractVecOrMat) (as well as the one with the scalar coefficients) which does the dimension check (without giving the user a possibility to disable it), and which then calls some internal method _unsafe_mul! which does no check. All the current mul! implementations for the different types of linear maps could become _unsafe_mul! methods with all dimension checks removed.

[Jutho]

When all of this is merged, I will probably go through the complete package code locally at my computer before we tag version 3. I've now skipped the parts I am less familiar with because you implemented it (blockmap, kronecker stuff).

[dkarrasch]

Thanks @Jutho for the thorough review and your comments. Perhaps I went overboard with the whole at-inbounds stuff. I agree that this is not crucial for performance, so I removed all of it.

I thought a lot about the "single entry point" approach, especially in the context of leveraging the MulStyle trait to passing intermediate storage place, but I couldn't wrap my head around it. The idea of my general approach here is to have only mul!s, and do everything else via dispatch. This means in particular, that users who have custom map types write mul! methods as if they had custom matrix types. Period. I find it somewhat ugly to ask them to write _unsafe_mul! (or any other _helper_function_name!) methods and leave mul! itself to us. I guess size checking is not such a big deal that it's worth introducing such a "workaround", and the intermediate storage issue may be something that we may wish to never actually do.

I will probably go through the complete package code locally at my computer before we tag version 3.

Sure, that's a good opportunity now.

[Jutho]

I thought a lot about the "single entry point" approach, especially in the context of leveraging the MulStyle trait to passing intermediate storage place, but I couldn't wrap my head around it. The idea of my general approach here is to have only mul!s, and do everything else via dispatch. This means in particular, that users who have custom map types write mul! methods as if they had custom matrix types. Period. I find it somewhat ugly to ask them to write _unsafe_mul! (or any other _helper_function_name!) methods and leave mul! itself to us. I guess size checking is not such a big deal that it's worth introducing such a "workaround", and the intermediate storage issue may be something that we may wish to never actually do.

Not saying that you have to, but if you would want to get rid of the extra dimension checkings due to the nested structure of mul!s, you could have the following:

• a single entry point mul!(Y::AbstractVecOrMat, L::LinearMap, X::AbstractVecOrMat) that does dimension checking and calls _unsafe_mul!

• an implementation of _unsafe_mul! for all of the specific LinearMap subtypes defined in this package

• a fallback _unsafe_mul!(Y::AbstractVecOrMat, L::LinearMap, X::AbstractVecOrMat) that calls back to mul!, and would thus never be called with the types from within this package, but ensures that new subtypes only need to define mul!

There was a discussion, or at least a post, about disconnecting method tables:
https://discourse.julialang.org/t/tiered-dispatch-algebra/45866

Unfortunately no one responded yet, cause I am interested in this. Not for this package in particular, but for other packages where I am also facing issues with excessive compilation times.

However, I am not sure that the premise made in that post is correct, and I will probably draft a response now.

[dkarrasch]

Aha, that proposal sounds clean and easy enough, I'll give that a try. I was thinking about such a redirection approach, but always tried to mix that with MulStyle and intermediate storage issue... and eventually failed. But without the storage stuff, this seems perfectly doable.

[dkarrasch]

Argh, that doesn't help. Within the mul! methods for nested types, we need to call mul!, in case somebody has an own map type in a linear combination. But these calls again would do a size check, so this approach doesn't eliminate downstream size checking. Ok, no worries, I'll just calm down as this is not a performance issue.

[Jutho]

I am not sure I am following; though I haven't thought this through as much. Within the _unsafe_mul! methods, you just keep calling _unsafe_mul!. As long as the linear maps involved are the ones in this package, this will not perform additional dimension checks. If a user-defined LinearMap type is encountered, the fallback _unsafe_mul! kicks in which calls mul! on that user defined type. So that one will do an additional dimension check probably, but that one is indeed unavoidable.

[dkarrasch]

I forgot the fallback. 🤦 I'll try again.

[dkarrasch]

Alright, this has the acclaimed feature, plus the minimal size checking thanks to @Jutho's suggestion and is basically 100% tested. Soon, we will have a code coverage of 110%, at least! 🤣 I think I'll merge this and then @Jutho can take a fresh look at it, meanwhile I'll work on the custom map documentation and the show methods.