SciML / LinearSolve.jl

LinearSolve.jl: High-Performance Unified Interface for Linear Solvers in Julia. Easily switch between factorization and Krylov methods, add preconditioners, and all in one interface.
https://docs.sciml.ai/LinearSolve/stable/
Other
245 stars 52 forks source link

Phantom allocations: need to define effects? #245

Closed ChrisRackauckas closed 1 year ago

ChrisRackauckas commented 1 year ago

After https://github.com/SciML/LinearSolve.jl/pull/244 we see:

using LinearSolve, LinearAlgebra, Random, BenchmarkTools
N = 100

A1 = Diagonal(rand(N))
A2 = copy(A1)
b1 = rand(N)
b2 = copy(b1)
x = similar(b1)

println("\nNaive:")
@btime x = $A2\$b2
x = A2\b2
# 98.419 ns (1 allocation: 896 bytes)

println("\npre-allocated:")
@btime $x .= $A2.diag.\$b2
# 23.370 ns (0 allocations: 0 bytes)

println("\nLinearSolve1:")
sol = @btime begin 
    prob = LinearProblem($A2, $b2)
    cache1 = init(prob,alias_A = true,alias_b = true)
    solve(cache1).u
end
# 181.523 ns (3 allocations: 1.05 KiB)

prob = LinearProblem(A1, b1)
cache1 = init(prob)
sol = solve(cache1)
println("\nLinearSolve2:")
sol = @btime begin 
    cache2 = LinearSolve.set_b($cache1, $b2)
    cache2 = LinearSolve.set_A(cache2, $A2)
    solve(cache2);
end
# 127.021 ns (2 allocations: 176 bytes)

We need to figure out how to get rid of those phantom allocations and why the compiler does not remove them in simple solves like this.

oscardssmith commented 1 year ago

Fixed on Julia 1.10 (and I'm pretty sure 1.9 also)

julia> @btime $x .= $A2.diag.\$b2;
  47.766 ns (0 allocations: 0 bytes)

julia> sol = @btime begin 
           prob = LinearProblem($A2, $b2)
           cache1 = init(prob,alias_A = true,alias_b = true)
           solve(cache1).u
       end;
  117.728 ns (1 allocation: 896 bytes)
julia> sol = @btime begin 
           cache2 = LinearSolve.set_b($cache1, $b2)
           cache2 = LinearSolve.set_A(cache2, $A2)
           solve(cache2);
       end;
  47.767 ns (0 allocations: 0 bytes)