fix perf in exp(::Matrix) (for smallish matrices)

[KristofferC]

Fixes JuliaLang/LinearAlgebra.jl#565

Before:

julia> @btime exp([1. 0; 2 0])
  18.912 μs (73 allocations: 4.80 KiB)

After:

julia> @btime exp([1. 0; 2 0])
  3.529 μs (25 allocations: 2.30 KiB)

Edit: The text below is no longer relevant:

Should the broadcasting be able to do better here? cc @mbauman Here is a MWE

f(A, A2, A4, A6, CC, Inn)      =  A * (A6 * (CC[14]*A6  + CC[12]*A4  + CC[10]*A2)  + CC[8]*A6  + CC[6]*A4  + CC[4]*A2  + CC[2]*Inn)
f_dots(A, A2, A4, A6, CC, Inn) =  A * (A6 * (CC[14]*A6 .+ CC[12]*A4 .+ CC[10]*A2) .+ CC[8]*A6 .+ CC[6]*A4 .+ CC[4]*A2 .+ CC[2]*Inn)
CC = rand(100); A = rand(2,2); A2 = A*A; A4 = A2*A2; A6 = A2 * A4; Inn = rand(2,2);

@btime f($A, $A2, $A4, $A6, $CC, $Inn)
  6.949 μs (28 allocations: 1.69 KiB)

@btime f_dots($A, $A2, $A4, $A6, $CC, $Inn)
  565.278 ns (11 allocations: 1.20 KiB)

Update since it is being discussed: using full fusion:

f_ddots(A, A2, A4, A6, CC, Inn) =  A * (A6 * @.((CC[14]*A6  + CC[12]*A4  + CC[10]*A2)  + CC[8]*A6  + CC[6]*A4  + CC[4]*A2  + CC[2]*Inn))

@btime f_ddots($A, $A2, $A4, $A6, $CC, $Inn)
  1.787 μs (7 allocations: 464 bytes)

[stevengj]

CC[14]*A6 etcetera should be CC[14].*A6 in order to fuse. Probably cleaner to just do @.(CC[14]*A6 + CC[12]*A4 + …)

[stevengj]

(I'm not sure why CC is an array rather than a bunch of scalars. Not that any of these micro-optimizations matter for large arrays.)

[KristofferC]

That makes the code three times slower so I rather not do that.

[KristofferC]

Like, I put a MWE on the top so we can talk about actual performance instead of theoretical performance :P

[mbauman]

That looks like #28126, particularly since when I tried it I saw a 2x gain. :-\

julia> f_dots2(A, A2, A4, A6, CC, Inn) =  A * (A6 * (CC[14].*A6 .+ CC[12].*A4 .+ CC[10].*A2) .+ CC[8]*A6 .+ CC[6].*A4 .+ CC[4].*A2 .+ CC[2].*Inn)
f_dots2 (generic function with 1 method)

julia> @btime f_dots2($A, $A2, $A4, $A6, $CC, $Inn)
  254.430 ns (5 allocations: 560 bytes)
2×2 Array{Float64,2}:
 0.516647  0.311652
 0.186972  0.116884

[KristofferC]

Ok, I also get that, but not when using the macro... Did I mess up?

[KristofferC]

Anyway, just pushed it with explicit dots.

[mbauman]

Ah, the difference appears to be in the n-ary operations. The parser doesn't use special n-ary parsing for .+ but the @. macro does.

[stevengj]

If we really wanted to we could split out the + CC[2]*Inn operation into a separate loop that just updates the diagonals in-place.

In general, it looks like there are lots of opportunities for avoiding temporary array allocations in this routine if we care, both by re-using arrays and by using things like rmul!.

It will only matter for tiny arrays like the one in #29113, however, since for larger arrays the Θ(n³) operations will dominate.

[Keno]

Do we have this as a nanosoldier benchmark and didn't notice? If not, add it?

[KristofferC]

Benchmark label is on there so anyone feel free to add it.