SciML / ExponentialUtilities.jl

Fast and differentiable implementations of matrix exponentials, Krylov exponential matrix-vector multiplications ("expmv"), KIOPS, ExpoKit functions, and more. All your exponential needs in SciML form.
https://docs.sciml.ai/ExponentialUtilities/stable/
Other
93 stars 29 forks source link

Higher allocations using `MMatrix` and `exponential!` #164

Open alblaz opened 7 months ago

alblaz commented 7 months ago

The snippet below shows that for small matrices using MMatrix as argument of exponential! is faster than using a regular Matrix. However the MMatrix allocates more. Is this the expected behaviour?

using StaticArrays, ExponentialUtilities
import ExponentialUtilities.alloc_mem

method=ExpMethodHigham2005();

A = [3. 2; 3 4.];
cacheA = alloc_mem(A,method);
println("Matrix");
@time exponential!(A, method,cacheA);

B  = @MMatrix [1 2; 3 4.];
cacheB = alloc_mem(B,method);
println("MMatrix");
@time exponential!(B, method,cacheB);

Matrix 0.000037 seconds (1 allocation: 80 bytes) MMatrix 0.000006 seconds (6 allocations: 416 bytes)

ChrisRackauckas commented 7 months ago
function ldiv_for_generated!(C, A, B) # C=A\B. Called from generated code
    F = lu!(A) # This allocation is unavoidable, due to the interface of LinearAlgebra
    ldiv!(F, B) # Result stored in B
    if (pointer_from_objref(C) != pointer_from_objref(B)) # Aliasing allowed
        copyto!(C, B)
    end
    return C
end

The array version has an allocation due to BLAS. This can be removed by using LinearSolve.jl. But note that there are sub-allocations in that not tracked by GC because it's in the BLAS workspace.

The MMatrix version is using GenericSchur.jl. You can see the extra allocations here:

https://github.com/RalphAS/GenericSchur.jl/blob/master/src/balance.jl#L45-L49

It would be nice to improve the alloc_mem interface to allocate those and then extend the GenericSchur.jl balance! function so that they can be passed in. If anyone has the time to do that then this is a nice straightforward issue to improve performance a bit.