New implementation of _exp! based on generated code

[jarlebring]

This is a new implementation of non-allocating matrix exponential _exp!. This implementation is better in terms of efficiency since it fixes #63, but also makes more use of dot fusion.

julia> using BenchmarkTools, LinearAlgebra, Random
julia> Random.seed!(0);
julia> n=1024*2;
julia> A=randn(n,n);
julia> A=0.85*A/opnorm(A,1);
julia> cache=(similar(A),similar(A),similar(A),similar(A),similar(A));
julia> Eold=@btime ExponentialUtilities._baseexp!(copy($A),caches=$cache);
  1.614 s (7 allocations: 32.03 MiB)
julia> Enew=@btime ExponentialUtilities._exp!(copy($A),caches=$cache);
  1.415 s (53 allocations: 32.04 MiB)
julia> norm(Eold-Enew)
1.696649934880783e-15

The implementation is based on generated code. It is generated using the package https://github.com/matrixfunctions/GraphMatFun.jl and the files exp_XX.jl in exp_generated are generated by the script https://github.com/jarlebring/ExponentialUtilities.jl/blob/master/src/exp_generated/generate_exp_code.jl

[YingboMa]

Haha, this is probably good for precompiles. @ChrisRackauckas will like this :-p

[YingboMa]

I assume this is an internal function that shouldn't be called besides testing?

[jarlebring]

I see _exp2! as equivalent to _exp! or a replacement thereof. When we are done comparing _exp! and _exp2!, we could consider renaming _exp! <- _exp2!. The function _exp! is currently not exported, although it is documented on the main page.

[jarlebring]

Any specific thoughts about using generated code?

The functions exp_6!,... exp_13! are very similar, and one could reduce them to one by using the identity exp_7!(A)=exp_6!(A/2)^2 etc. It would probably require more computation (a scaling of the matrix) but the added computation might be acceptable to get less code.

[ChrisRackauckas]

It's fine, used the generated code.

[ChrisRackauckas]

Test failure looks real

[jarlebring]

I believe E=_exp!(A) writes the solution into A, whereas E=_exp2!(A) has destroyed the contents of A but it does not contain the solution. E points elsewhere. When I run tests locally (julia 1.7-dev), this seems to be the issue. Update: Fixed.

The CI seems to complain about something related to findfirst...

[jarlebring]

Thoughts about replacing _exp with _exp2? Any important users/use-cases that should be checked?

[ChrisRackauckas]

Check whether the generated kernels for CUDA work (probably not without changing the LAPACK calls?).

The main use cases are the exponential integrators https://github.com/SciML/ExponentialUtilities.jl/blob/5460b95594a9e23f418710e568db3a92d7bcfd26/src/krylov_phiv.jl .

It might be interesting to run the exprk benchmark in https://benchmarks.sciml.ai/html/MOLPDE/allen_cahn_fdm_wpd.html and see how that's effected.

[jarlebring]

Currently, CI is complaining about mul!. This a version issue. The code does identity addition with JuliaLang/julia#40731, which is not in the CI version of julia (1.6.2). If we want support for 1.6.2, I could change it by basically including the diagonal add in 40731 as a function uniformscaling_add!.

[jarlebring]

Check whether the generated kernels for CUDA work (probably not without changing the LAPACK calls?).

If we disable balancing, there is only one LAPACK call: the linear solve LAPACK.gesv!.

https://github.com/jarlebring/ExponentialUtilities.jl/blob/src/exp_generated/exp_5.jl#L77

What should the be in a CUDA setting?

[jarlebring]

Okay. lu!+ldiv! seems to work in general. It doesn't work for MMatrix since lu! is not defined for that type but it could be easily be fixed with a dispatch on ldiv_for_generated!. For the moment, I will not add that since I'm unsure how important MMatrix is here.

[ChrisRackauckas]

rename the other to baseexp!?

And do a test on the exponential integrators and this is good to go.

[jarlebring]

rename the other to baseexp!?

You mean _exp_old! -> _baseexp!? Does something speak against removing the old implementation? Less code to maintain. I don't see much advantage of the old code, except maybe that it is not generated.

[jarlebring]

If we keep both, and see the exponential integrators as the main application, I think we need the feature to specify which of the two non-allocating exp implementations should be used, eg here

ExponentialUtilities.jl/src/krylov_phiv.jl

Lines 85 to 90 in 5460b95

	F = eigen!(SymTridiagonal(cache))
	expHe = F.vectors * (exp.(lmul!(t,F.values)) .* @view(F.vectors[1, :]))
	else
	lmul!(t, cache); expH = cache
	_exp!(expH)
	expHe = @view(expH[:, 1])

I think such features should only be added if they add value.

[jarlebring]

When we are anyway messing with this code: _exp! is not exported but it is described on the main documentation page. How about renaming it to exp_noalloc! and exporting it?

[ChrisRackauckas]

I think we mainly want to keep it around to potentially upstream it.

[jarlebring]

Alright. Fine with me. Changed to _baseexp!.

And do a test on the exponential integrators and this is good to go.

But aren't they already tested like in the test Phi?

[ChrisRackauckas]

Oh I just noticed, the README should get an update too. Follow up PR?

[jarlebring]

Do you mean documentation of the kwargs in the README.md?

[ChrisRackauckas]

that no longer describes the method of _exp! and _baseexp! is not documented, right?

[ChrisRackauckas]

or is the new one based off of the 2008 paper?

[jarlebring]

They are both implementations of the Higham 2008 paper.

[ChrisRackauckas]

ahh okay, so all that's needed is a baseexp section then.