Add mutating versions of functions that reuse external storage for sigma points and weights

[bvdmitri]

Hey! Thanks for this amazing package. Would it be possible to add mutating version of all functions that would support and reuse external storage for sigma points and weights instead of allocating a new one each time, e.g:

n = 100

sigmap  = Vector{Float64}(undef, n)
weights = Vector{Float64}(undef, n)

alpha = 1.0

gausslaguerre!(sigmap, weights, n, alpha)

Usually it is a good practice in Julia to have such functions, e.g. we have map!, filter! etc.

The benefit of this would be improved performance in cases where you need continuously recompute sigma points and weights for different alpha parameter for example (for gausslaguerre). Non-mutating versions would be simply:

# Float64 here is just an example, maybe support type for storage as well?
gausslaguerre(n, alpha) = gausslaguerre!(Vector{Float64}(undef, n), Vector{Float64}(undef, n), n, alpha)

[dlfivefifty]

We would be happy to accept a PR for this.

[bvdmitri]

👍 I'll try to find time to implement this feature and return to you if I need help.

[bvdmitri]

@dlfivefifty I started to implement this feature for gausslaguerre and noticed that some of your tests are incorrect. E.g. this one:

FastGaussQuadrature.jl/test/test_gausslaguerre.jl

Line 82 in 7b256f3

@test abs(x[37] - 98.388267163326702) < tol

It should be

@test abs(x_gw[37] - 98.388267163326702) < tol

Notice different x and x_gw. Do you want me to open another issue or it can be a part of my PR for in-place mutating functions. I'm not sure though how long it would take me to go through your codebase and add mutating versions for all of your functions.

P.S. fixed tests do pass

[bvdmitri]

Updated Gauss-Laguerre looks some-what like this:

bvdmitri@3146d9d

I also added more place for external storage re-usage and as a result we save some allocations:

function test()
     dummy = 0.0
     for n in 1:1000
         x, w = gausslaguerre(n)
         dummy += sum(w)
     end
     return dummy
 end

# master
julia> @btime test()
  762.632 ms (54755 allocations: 12.85 MiB)
999.9999999999991
# after update
julia> @btime test()
  744.996 ms (53507 allocations: 12.60 MiB)

It is also interesting to see how much do we safe with in-place versions:

function test_reuse(x, w)
    dummy = 0.0
    for n in 1:1000
        x, w = FastGaussQuadrature.gausslaguerre!(x, w, n)
        dummy += sum(view(w, 1:n))
    end
    return dummy
end

julia> @btime test_reuse(s[1], s[2]) setup=(s=(Vector{Float64}(undef, 1000), Vector{Float64}(undef, 1000)))
  743.897 ms (51507 allocations: 4.69 MiB) # third of the original allocations
999.9999999999991

I'm not sure where do this remaining allocations come from. Probably there are some other places that I overlooked. It looks that saved extra allocations do not contribute much in the overall performance, but it may affect performance in larger programs with a lot of GC pressure I suppose.

[dlfivefifty]

The code is mostly a Matlab port so there’s a lot of room food improvement