avoid allocation for length-1 linear combinations

[dkarrasch]

No description provided.

[codecov]

Codecov Report

Merging #34 into master will increase coverage by 15.2%.
The diff coverage is 68.42%.

@@            Coverage Diff             @@
##           master      #34      +/-   ##
==========================================
+ Coverage   60.94%   76.14%   +15.2%     
==========================================
  Files           7        7              
  Lines         338      348      +10     
==========================================
+ Hits          206      265      +59     
+ Misses        132       83      -49

Impacted Files	Coverage Δ
src/linearcombination.jl	82.05% <68.42%> (+19.55%)	⬆️
src/LinearMaps.jl	67.6% <0%> (+0.46%)	⬆️
src/composition.jl	76.53% <0%> (+6.12%)	⬆️
src/functionmap.jl	75.67% <0%> (+16.21%)	⬆️
src/transpose.jl	70.37% <0%> (+29.62%)	⬆️
src/identitymap.jl	75% <0%> (+45.58%)	⬆️
src/wrappedmap.jl	94.11% <0%> (+52.94%)	⬆️

---

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 b425454...df4fcb3. Read the comment docs.

[Jutho]

Looks good; if Julia would expose the true BLAS functionality of gemv, which can do y = beta*y+ alpha* A * x , we wouldn't need a temporary at all.

[dkarrasch]

I think we need the temporary for FunctionMaps anyway, don't we? We can't compute the image vector "componentwise", so we need to store the A*x somewhere. But now that you say it, I think there is a 5-args mul! somewhere in Base, at least there used to be. Certainly, we have it for LinearMaps. I'll rewrite it and see if it works with WrappedMaps. I can't seem to find a 5-args method for mul! in LinearAlgebra, and I remember an issue in IterativeSolvers.jl.

EDIT: There is a discussion here and a PR, so once this is finished and merged, we may need to adopt our own API to that. Then, we may be able to avoid allocations for LinearMaps built from matrices.

[Jutho]

Indeed for FunctionMaps we do, unless we require the user to implement a 5 argument mul!. I was a bit too fast here :)

[Jutho]

This looks good. Is this ready to be merged?

[chriscoey]

Discussion on the 5-arg mul! seems to have stalled, which makes me sad

[dkarrasch]

Yes, I think this can be merged. We were seriously hit by code coverage changes, so we dropped down to 60%. I started a major revision of the tests, which has significantly improved coverage. I wanted to take another look at it, but I can do it in a separate PR (I have some other minor things to suggest).

[dkarrasch]

Regarding the 5-arg mul!, we will need to be very careful. Our current mul! for LinearCombinations allocates exactly once (z = similar(y)) when at least two LinearMaps are combined (given that the LinearMap is built from matrices and/or mutating FunctionMaps). If we replace our current "3-arg mul" A_mul_B!-axpy!-combination by a 5-arg mul! in each step, we may get rid of any extra allocation in case everything is built from matrices (and the 5-arg mul! is provided, of course), but we would allocate once per linear map in the linear combination, see our own current 5-arg mul! which allocates A*x at some point. So, either we find a nice workaround, or live with our current version, which would be not too bad.