Bunchkaufman- and LU-decomposition based generalized eigenvalues and eigenvectors

[aravindh-krishnamoorthy]

Introduction

This PR implements Bunchkaufman (BK)- and LU-decomposition based generalized eigenvalues and eigenvectors of $\boldsymbol{A}$ and $\boldsymbol{B}$ when $\boldsymbol{A}$ is arbitrary and $\boldsymbol{B}$ is (Hermitian) symmetric in the following new methods in stdlib/LinearAlgebra/src/symmetriceigen.jl:

1. Functions eigvals and eigen
2. In-place versions: eigvals! and eigen!

For $N \geq 32,$ these implementations outperform the equivalent native (eigen.jl) functions in both execution time and memory usage. Test and benchmark results will be updated here regularly.

Comments are suggestions are welcome.

• Note: This work is an overspill from PR Rework symmetric generalized eigen/eigvals #49673. The comments there may also be relevant for this PR.

Open points

• Outputs when $B$ is singular. (Output is analogous to the Cholesky solution.)

[aravindh-krishnamoorthy]

Benchmark Results

Benchmarking file and instructions

• The benchmarking file is here (rename to bkmark.jl: bkmark.txt
• Execute after checking out the PR branch:

julia> using Revise
julia> using LinearAlgebra
julia> Revise.track(LinearAlgebra)
julia> using BenchmarkTools
julia> N=32; include("/home/ark/temp/bkmark.jl")
julia> N=1000; include("/home/ark/temp/bkmark.jl")

Results (06.12.2023)

Please see the comment diff for a comparison with the previous benchmark results. Please note that the absolute values of these figures may depend on external factors, while relative values and memory figures may provide intuition into the algorithm complexity.

N=32

Function	Native f = (x->x)	BK f = bunchkaufman	LU f = lu
eigvals(A, f(B)) (Real)	165.657 μs (28 allocations: 68.41 KiB)	138.446 μs (51 allocations: 78.81 KiB)	134.494 μs (41 allocations: 77.53 KiB)
eigen(A, f(B)) (Real)	268.748 μs (33 allocations: 92.02 KiB)	227.372 μs (67 allocations: 119.25 KiB)	228.709 μs (57 allocations: 117.97 KiB)
eigvals!(A, f!(B)) (Real)	185.929 μs (22 allocations: 52.25 KiB)	131.394 μs (39 allocations: 54.22 KiB)	126.212 μs (29 allocations: 52.95 KiB)
eigen!(A , f!(B)) (Real)	279.662 μs (27 allocations: 75.86 KiB)	224.781 μs (55 allocations: 94.66 KiB)	211.973 μs (45 allocations: 93.39 KiB)
eigvals(A, f(B)) (Complex)	601.306 μs (26 allocations: 134.33 KiB)	417.711 μs (49 allocations: 153.64 KiB)	411.837 μs (39 allocations: 151.31 KiB)
eigen(A , f(B)) (Complex)	983.353 μs (30 allocations: 150.17 KiB)	875.138 μs (55 allocations: 169.78 KiB)	794.085 μs (44 allocations: 167.44 KiB)
eigvals!(A, f!(B)) (Complex)	677.062 μs (26 allocations: 134.33 KiB)	490.176 μs (37 allocations: 105.05 KiB)	481.066 μs (27 allocations: 102.73 KiB)
eigen!(A , f!(B)) (Complex)	1.039 ms (30 allocations: 150.17 KiB)	851.008 μs (43 allocations: 121.19 KiB)	800.396 μs (32 allocations: 118.86 KiB)

N=1000

Function	Native f = (x->x)	BK f = bunchkaufman	LU f = lu
eigvals(A, f(B)) (Real)	2.008 s (36 allocations: 16.83 MiB)	492.375 ms (64 allocations: 23.74 MiB)	495.089 ms (50 allocations: 38.48 MiB)
eigen(A, f(B)) (Real)	2.592 s (40 allocations: 39.70 MiB)	655.926 ms (80 allocations: 62.62 MiB)	636.457 ms (66 allocations: 77.36 MiB)
eigvals!(A, f!(B)) (Real)	2.004 s (30 allocations: 1.57 MiB)	443.697 ms (51 allocations: 860.56 KiB)	448.144 ms (37 allocations: 15.58 MiB)
eigen!(A , f!(B)) (Real)	2.687 s (34 allocations: 24.45 MiB)	699.842 ms (67 allocations: 39.73 MiB)	706.874 ms (53 allocations: 54.47 MiB)
eigvals(A, f(B)) (Complex)	6.734 s (32 allocations: 33.62 MiB)	1.546 s (61 allocations: 47.40 MiB)	1.465 s (47 allocations: 76.87 MiB)
eigen(A , f(B)) (Complex)	7.470 s (36 allocations: 48.87 MiB)	1.967 s (67 allocations: 64.12 MiB)	2.075 s (52 allocations: 93.60 MiB)
eigvals!(A, f!(B)) (Complex)	6.261 s (32 allocations: 33.62 MiB)	1.432 s (48 allocations: 1.61 MiB)	1.415 s (34 allocations: 31.09 MiB)
eigen!(A , f!(B)) (Complex)	7.949 s (36 allocations: 48.87 MiB)	2.061 s (54 allocations: 18.33 MiB)	2.039 s (39 allocations: 47.81 MiB)

Test Results (06.12.2023)

Testing file and instructions

• Tests are located in stdlib/LinearAlgebra/test/symmetriceigen.jl
• Execute after checking out the PR branch:

julia> using Revise
julia> using LinearAlgebra
julia> Revise.track(LinearAlgebra)
julia> include("symmetriceigen.jl")

Results

WARNING: replacing module TestSymmetricEigen.
Test Summary:      | Pass  Total  Time
chol-eigen-eigvals |    4      4  0.9s
Test Summary: | Pass  Total  Time
issue #49533  |   12     12  0.4s
Test Summary:       | Pass  Total  Time
bk-lu-eigen-eigvals |   24     24  1.7s
Main.TestSymmetricEigen

[aravindh-krishnamoorthy]

Hello @dkarrasch, @jishnub, and @schneiderfelipe, I look forward to your comments and suggestions on this new PR.

[stevengj]

Is there a reference for this algorithm?

[aravindh-krishnamoorthy]

Is there a reference for this algorithm?

Hello @stevengj. Unfortunately, I don't have a reference handy. But, the solution follows the same method as that of the Cholesky decomposition one in eigen/Cholesky and that in LAPACK (NIST)'s chegv(), and is explained below.

Bunch-Kaufman Decomposition

Consider an arbitrary matrix $\boldsymbol{A}$ and an invertible(1) Hermitian symmetric matrix $\boldsymbol{B}$ with Bunch-Kaufman (BK) decomposition given by $$\boldsymbol{B} = \boldsymbol{P}^{\mathrm{H}}\boldsymbol{U}\boldsymbol{D}\boldsymbol{U}^{\mathrm{H}}\boldsymbol{P},$$ where $\boldsymbol{P}$ is a permutation matrix, $\boldsymbol{U}$ is an invertible upper (lower) triangular matrix, and $\boldsymbol{D}$ is a Hermitian symmetric $2\times 2$-block diagonal matrix (a special case of tridiagonal), as given in LinearAlgebra.bunchkaufman.

Generalised eigenvalues

We note that the roots of the generalised eigenvalue problem $$\det(\boldsymbol{A} - \lambda\boldsymbol{B})$$ and those of $$\det\big(\boldsymbol{D}^{-1}\boldsymbol{U}^{-1}(\boldsymbol{P}\boldsymbol{A}\boldsymbol{P}^{\mathrm{H}})(\boldsymbol{U}^{\mathrm{H}})^{-1} - \lambda\boldsymbol{I}\big)$$ are equivalent. Hence, the generalised eigenvalues of $\boldsymbol{A}$ and $\boldsymbol{B}$ are the eigenvalues of $\boldsymbol{A}' = \boldsymbol{D}^{-1}\boldsymbol{U}^{-1}(\boldsymbol{P}\boldsymbol{A}\boldsymbol{P}^{\mathrm{H}})(\boldsymbol{U}^{\mathrm{H}})^{-1}.$ Let $\boldsymbol{\Lambda}$ be a matrix with the eigenvalues of $\boldsymbol{A}'$ on the main diagonal.

Generalised eigenvectors

Now, let $\boldsymbol{V}'$ contain the eigenvectors of $\boldsymbol{A}'$ in its columns. Then, the generalised eigenvectors, based on the above, are contained in the columns of $$\boldsymbol{V} = \boldsymbol{P}^{\mathrm{H}}(\boldsymbol{U}^{\mathrm{H}})^{-1}\boldsymbol{V}'.$$

LU decomposition

Let $\boldsymbol{B}$ be a square matrix (generalization of Hermitian and symmetric matrices) with LU decomposition given by
$$\boldsymbol{B} = \boldsymbol{P}\boldsymbol{L}\boldsymbol{U},$$ where $\boldsymbol{P}$ is a permutation matrix, and $\boldsymbol{L}$ and $\boldsymbol{U}$ are invertible lower and upper triangular matrices, respectively.

Generalised eigenvalues and eigenvectors

Analogous to the BK-decomposition case, the generalised eigenvalues are the eigenvalues of eigenvalues of $$\boldsymbol{A}' = \boldsymbol{L}^{-1}(\boldsymbol{P}\boldsymbol{A})\boldsymbol{U}^{-1}.$$ Now, let $\boldsymbol{V}'$ contain the eigenvectors of $\boldsymbol{A}'$ in its columns. Then, the generalised eigenvectors, analogous to the BK-decomposition case, are contained in the columns of $$\boldsymbol{V} = \boldsymbol{U}^{-1}\boldsymbol{V}'.$$

Validation

For both decompositions, the eigenvalues and eigenvectors satisfy the relation $$\boldsymbol{A}\boldsymbol{V} = \boldsymbol{B}\boldsymbol{V}\boldsymbol{\Lambda},$$ which is utilized for testing the functions.

Comments and suggestions are welcome.

(1) For non-invertible $\boldsymbol{B},$ the usual convention is to mark the corresponding eigenvalue with infinity. Julia (complex double) returns NaN+Nan*im. This algorithm returns numbers of the order $10^{16}$ in place of infinity. A discussion on this is still open, please see checkbox 3 in the PR text.

[dkarrasch]

Thanks for tackling this. I see this is addressing one of my major concerns, that the algorithm is not really in-place of the left-hand matrix! Cool. Here are a couple of first comments.

[aravindh-krishnamoorthy]

Thanks for tackling this. I see this is addressing one of my major concerns, that the algorithm is not really in-place of the left-hand matrix! Cool. Here are a couple of first comments.

Hello @dkarrasch. Thank you for the first round of comments. I've responded to the comments individually and incorporated the changes via a single commit aravindh-krishnamoorthy@8b8f556.

Below, I'd like to address the following comment.

Why changing the sorting function in unrelated tests?

In the tests in test/symmetriceigen.jl, besides the test $\boldsymbol{A}\boldsymbol{V} = \boldsymbol{B}\boldsymbol{V}\boldsymbol{\Lambda},$ the outputs of dense.jl's generalised eigenvalues and those from BK are compared. Gaussian random matrices are used as test matrices (1). Hence, the test matrices are almost surely full-rank and have eigenvalue separation (magnitude) $> \epsilon.$ Furthermore, most generalized eigenvalues are in pairs $(a \pm \mathrm{i}b).$

However, due to numerical issues, the eigenvalue pairs between the functions differ as $(a + \mathrm{i}b)$ and $(a + \epsilon - \mathrm{i}b),$ which reverses the ordering causing the comparisons to fail. Arranging them by x->(imag(x),real(x)) solves this issue.

(1) However, even though using Gaussian random matrices as test matrices is advantageous (greater coverage as the matrices vary every test), it is still possible that we end up with rank deficient or closely spaced eigenvalues. In such cases, the resulting test failures will be difficult to reproduce. If you do not prefer to have this uncertainty, I can locally check with Gaussian random matrices but have a fixed matrix in the official test.

Kindly let me know your opinion.

[dkarrasch]

I think this is a great opportunity to clean-up/extend a few things. I'd suggest to use more generic functions to achieve the same thing, but for specific cases use LAPACK functions. For instance:

• regarding the row and column permutation, we already have

  julia/stdlib/LinearAlgebra/src/lu.jl

  Lines 373 to 418 in 8dd58a4

  	_apply_ipiv_rows!(A::LU, B::AbstractVecOrMat) = _ipiv_rows!(A, 1 : length(A.ipiv), B)
  	_apply_inverse_ipiv_rows!(A::LU, B::AbstractVecOrMat) = _ipiv_rows!(A, length(A.ipiv) : -1 : 1, B)
  	function _ipiv_rows!(A::LU, order::OrdinalRange, B::AbstractVecOrMat)
  	for i = order
  	if i != A.ipiv[i]
  	_swap_rows!(B, i, A.ipiv[i])
  	end
  	end
  	B
  	end
  	function _swap_rows!(B::AbstractVector, i::Integer, j::Integer)
  	B[i], B[j] = B[j], B[i]
  	B
  	end
  	function _swap_rows!(B::AbstractMatrix, i::Integer, j::Integer)
  	for col = 1 : size(B, 2)
  	B[i,col], B[j,col] = B[j,col], B[i,col]
  	end
  	B
  	end
  	_apply_ipiv_cols!(A::LU, B::AbstractVecOrMat) = _ipiv_cols!(A, 1 : length(A.ipiv), B)
  	_apply_inverse_ipiv_cols!(A::LU, B::AbstractVecOrMat) = _ipiv_cols!(A, length(A.ipiv) : -1 : 1, B)
  	function _ipiv_cols!(A::LU, order::OrdinalRange, B::AbstractVecOrMat)
  	for i = order
  	if i != A.ipiv[i]
  	_swap_cols!(B, i, A.ipiv[i])
  	end
  	end
  	B
  	end
  	function _swap_cols!(B::AbstractVector, i::Integer, j::Integer)
  	_swap_rows!(B, i, j)
  	end
  	function _swap_cols!(B::AbstractMatrix, i::Integer, j::Integer)
  	for row = 1 : size(B, 1)
  	B[row,i], B[row,j] = B[row,j], B[row,i]
  	end
  	B
  	end

  My suspicion is that it works the same for both LU and BunchKaufman (to be checked). In that case, one could simply extend the method signature. The question is whether LAPACK's version is any faster, in which case we might redirect the case of BlasFloat-StridedMatrixes to your added functions. If not, then we may not need to add the LAPACK functions.
• regarding the definition of M, du etc.: This looks like code duplication to getproperty of BunchKaufman. Then why not simply define D = B.D and M = B.uplo == 'U' ? B.U : B.L? If the issue is later that we need un-aliased off-diagonals, then I think we should change that in getproperty and make it un-aliased in all cases, and not sometimes yes and sometimes not. The extra effort is relatively small.
• The next lines have well-established generic function calls:
  • "# In-place A .= inv(M)*A" is ldiv!(M, A)
  • "# In-place A .= A*inv(M^H)" is rdiv!(A, M')
  • "# In-place A .= inv(Tridiagonal(dl, d = diag(LUD), du))*A" could be ldiv!(lu!(D), A).
  • "vecs .= inv(M^H)*vecs" is again just ldiv!(M', vecs)
  • "# In-place vecs .= vecs[invperm(B.p), :]" might be done with one of the apply_... functions referenced above.

[dkarrasch]

A few more suggestions for simplification and genericity. It would be good to use only generic functions (i.e., not use any LAPACK functions directly), and then not restrict on BlasFloat unnecessarily.

[aravindh-krishnamoorthy]

Hello @dkarrasch. Thank you very much for your support with this PR and your comments. Sorry for the delay in my response.

Please note that the intermediate version was inadvertently pushed to GitHub. It was primarily for my benchmarking use. I've removed the function and reverted the signatures.

I think this is a great opportunity to clean-up/extend a few things. I'd suggest to use more generic functions to achieve the same thing, but for specific cases use LAPACK functions.

Thank you. Indeed, my idea is to first understand the most optimal performance via BLAS/LAPACK and then ensure that the final Julia code can get there. I've incorporated most of the suggestions below. Unfortunately, in a few cases, the LAPACK routines were superior either in terms of complexity or memory usage. Any suggestions to rewrite them in Julia will be welcome. Details are as follows.

For instance:

• regarding the row and column permutation, we already have

  julia/stdlib/LinearAlgebra/src/lu.jl

  Lines 373 to 418 in 8dd58a4

  	_apply_ipiv_rows!(A::LU, B::AbstractVecOrMat) = _ipiv_rows!(A, 1 : length(A.ipiv), B)
  	_apply_inverse_ipiv_rows!(A::LU, B::AbstractVecOrMat) = _ipiv_rows!(A, length(A.ipiv) : -1 : 1, B)
  	function _ipiv_rows!(A::LU, order::OrdinalRange, B::AbstractVecOrMat)
  	for i = order
  	if i != A.ipiv[i]
  	_swap_rows!(B, i, A.ipiv[i])
  	end
  	end
  	B
  	end
  	function _swap_rows!(B::AbstractVector, i::Integer, j::Integer)
  	B[i], B[j] = B[j], B[i]
  	B
  	end
  	function _swap_rows!(B::AbstractMatrix, i::Integer, j::Integer)
  	for col = 1 : size(B, 2)
  	B[i,col], B[j,col] = B[j,col], B[i,col]
  	end
  	B
  	end
  	_apply_ipiv_cols!(A::LU, B::AbstractVecOrMat) = _ipiv_cols!(A, 1 : length(A.ipiv), B)
  	_apply_inverse_ipiv_cols!(A::LU, B::AbstractVecOrMat) = _ipiv_cols!(A, length(A.ipiv) : -1 : 1, B)
  	function _ipiv_cols!(A::LU, order::OrdinalRange, B::AbstractVecOrMat)
  	for i = order
  	if i != A.ipiv[i]
  	_swap_cols!(B, i, A.ipiv[i])
  	end
  	end
  	B
  	end
  	function _swap_cols!(B::AbstractVector, i::Integer, j::Integer)
  	_swap_rows!(B, i, j)
  	end
  	function _swap_cols!(B::AbstractMatrix, i::Integer, j::Integer)
  	for row = 1 : size(B, 1)
  	B[row,i], B[row,j] = B[row,j], B[row,i]
  	end
  	B
  	end

  
  My suspicion is that it works the same for both LU and BunchKaufman (to be checked). In that case, one could simply extend the method signature. The question is whether LAPACK's version is any faster, in which case we might redirect the case of BlasFloat-StridedMatrixes to your added functions. If not, then we may not need to add the LAPACK functions.

Unfortunately, the routines are for the LAPACK-style swap permutations and are best suited for the IPIV-style permutation vectors. Using these methods for permuting $A$ will entail additional complexity and memory.

For permuting $A,$ I feel that the LAPACK routines are the superior alternative. If necessary, these LAPACK routines can also be wrapped in a permute! :: StridedMatrix -> Vector -> StridedMatrix method to make them look more Julia-like.

In our case, indeed, obtaining B.L/U from B.LD in BK necessitates an IPIV permutation. However, presently, BK is only supported via LAPACK and, in this case, this job is best done by syconvf/_rook!, which additionally also extract the super/sub diagonals.

Hence, these routines may not be beneficial for our case.

• regarding the definition of M, du etc.: This looks like code duplication to getproperty of BunchKaufman. Then why not simply define D = B.D and M = B.uplo == 'U' ? B.U : B.L? If the issue is later that we need un-aliased off-diagonals, then I think we should change that in getproperty and make it un-aliased in all cases, and not sometimes yes and sometimes not. The extra effort is relatively small.

Unfortunately, the getproperty code of BK is inefficient. It calls syconvf/_rook! over-and-over again for each invocation, and additionally, utilizes memory, as seen below.

julia> B = complex.(randn(1000,1000),randn(1000,1000)) ; B = (B+B')/2 ;
julia> K = bunchkaufman(B) ;
julia> @time K.D ;
  0.002521 seconds (7 allocations: 15.320 MiB)
julia> @time K.uplo == 'U' ? K.U : K.L ;
  0.010232 seconds (5 allocations: 15.274 MiB)

On the other hand, as seen from the benchmark in the first comment above, the present code with duplicated getproperty is almost completely in place and requires less memory than even the generic counterpart.

• The next lines have well-established generic function calls:

  • "# In-place A .= inv(M)*A" is ldiv!(M, A)
  • "# In-place A .= A*inv(M^H)" is rdiv!(A, M')
  • "# In-place A .= inv(Tridiagonal(dl, d = diag(LUD), du))*A" could be ldiv!(lu!(D), A).
  • "vecs .= inv(M^H)*vecs" is again just ldiv!(M', vecs)

The above four suggestions are incorporated!

• "# In-place vecs .= vecs[invperm(B.p), :]" might be done with one of the apply_... functions referenced above.

Unfortunately, as explained earlier, this may not be possible.

[dkarrasch]

One quick suggestion, but actually, I think we should implement these in julia. After all, this is just a helper function. I'll post a suggestion shortly.

[dkarrasch]

Hi! I'd like to propose to add the following methods for column and row permutations in combinatorics.jl in Base:

permutecols!(a::AbstractMatrix, p::AbstractVector{<:Integer}) =
    _permute!(a, p, Base.swapcols!)
permuterows!(a::AbstractMatrix, p::AbstractVector{<:Integer}) =
    _permute!(a, p, Base.swaprows!)
@inline function _permute!(a::AbstractMatrix, p::AbstractVector{<:Integer}, swapfun!::F) where {F}
    require_one_based_indexing(a, p)
    p .= .-p
    for i in 1:length(p)
        p[i] > 0 && continue
        j = i
        in = p[j] = -p[j]
        while p[in] < 0
            swapfun!(a, in, j)
            j = in
            in = p[in] = -p[in]
        end
    end
    a
end
invpermutecols!(a::AbstractMatrix, p::AbstractVector{<:Integer}) =
    _invpermute!(a, p, Base.swapcols!)
invpermuterows!(a::AbstractMatrix, p::AbstractVector{<:Integer}) =
    _invpermute!(a, p, Base.swaprows!)
@inline function _invpermute!(a::AbstractMatrix, p::AbstractVector{<:Integer}, swapfun!::F) where {F}
    require_one_based_indexing(a, p)
    p .= .-p
    for i in 1:length(p)
        p[i] > 0 && continue
        j = p[i] = -p[i]
        while j != i
           swapfun!(a, j, i)
           j = p[j] = -p[j]
        end
     end
    a
end

They are straightforward translations of the LAPACK functions, and the row permutation (IIRC) is even faster than LAPACK. They are written such that p doesn't get destroyed in passing, so can be reused, in contrast to the currently used permutecols!!. We may decide to remove it or perhaps to redirect it to permutecols!.

[dkarrasch]

A few more comments. Moreover, there's a generic BunchKaufman decomposition in the pipeline (#51487), so we should keep this as generic as possible w.r.t. the eltype. After all, this is a generic linear algebra algorithm, not a LAPACK specific one.

[dkarrasch]

What a beautifully generic piece of code!!! If we keep the LAPACK permutation functions, then they need to be tested. Since they are new and still unused internally, perhaps we could leave them out here?

[aravindh-krishnamoorthy]

Thank you, @dkarrasch. Now, the code indeed is beautiful and generic. The main changes are done. I will review and benchmark everything again (esp with the new permutation routines) and incorporate other cosmetic changes soon.

[aravindh-krishnamoorthy]

@dkarrasch Perhaps it makes sense to wait for #51763 to be merged so that getproperties! in bunchkaufman.jl can be cleaned up in this PR a bit? The other points are resolved now.

[aravindh-krishnamoorthy]

Hello @dkarrasch. I see that #51763 is merged. Now, I'll resume work on this one incorporating those changes. Hope to finalise this PR soon.

[dkarrasch]

Last "hygiene" comment, and then this is finally ready to go. Please add an entry in NEWS.md that announces these features.

[aravindh-krishnamoorthy]

@dkarrasch This now looks good to go from my side. I looked at the benchmarks and they seem Ok now. For both cases (N=32, 1000), the new functions significantly outperform the generic function (not optimised for symmetric matrices). Furthermore, I feel it makes sense to retain both Bunchkaufman and LU functions since each is optimal in a certain scenario. Also, I've completed a few last-minute checks and updated NEWS.md.

Nevertheless, I do welcome any last-minute comments from you and will be happy to incorporate them -- ark

Edit: Sadly, there's always a residual bug :) Added commit 485756d.

[dkarrasch]

To be "fair", one should perhaps include the factorization of B in the measurement. Did you do that? In any case, we're not changing any defaults here, so users can benchmark their use cases and make a call. Having the factorization available may be helpful or even necessary in some cases, depending on the context.

[aravindh-krishnamoorthy]

To be "fair", one should perhaps include the factorization of B in the measurement. Did you do that? In any case, we're not changing any defaults here, so users can benchmark their use cases and make a call. Having the factorization available may be helpful or even necessary in some cases, depending on the context.

Indeed, the factorisation overheads due to bunchkaufman and lu are included in the runtime and memory benchmark figures given above. The new functions do well also when the calls the factorisation routines are considered!