Cholesky and Inverse Cholesky in one pass

[Marco-Congedo]

Hello,

I wrote a function for computing the Cholesky factor and its inverse in one pass. For small matrices this appears considerably faster than computing the Cholesky factor and inverting it with LinearAlgebra. Would that be interesting for this package?

For the moment being i am hosting it in the PosDefManifold.jl package
as the choInv function

Here is a quick benchmark:

using PosDefManifold, BenchmarkTools
P=randP(10)

using LinearAlgebra

function t(P)
       C=cholesky(P)
       return C.L, inv(C.L)
end

julia> @benchmark(t(P))
BenchmarkTools.Trial:
  memory estimate:  3.69 KiB
  allocs estimate:  12
  --------------
  minimum time:     40.800 μs (0.00% GC)
  median time:      47.201 μs (0.00% GC)
  mean time:        48.117 μs (0.00% GC)
  maximum time:     291.300 μs (0.00% GC)
  --------------
  samples:          10000
  evals/sample:     1

using function choInv

julia> @benchmark(choInv(P))
BenchmarkTools.Trial:
  memory estimate:  6.09 KiB
  allocs estimate:  23
  --------------
  minimum time:     1.530 μs (0.00% GC)
  median time:      1.690 μs (0.00% GC)
  mean time:        2.883 μs (14.91% GC)
  maximum time:     563.430 μs (98.69% GC)
  --------------
  samples:          10000
  evals/sample:     10

[dlfivefifty]

Please use triple tics (`) when you quote code, I've fixed your comment.

Yes this would be reasonable to add here, but it should be refactored as a Factorization object. Also the name choInv is inconsistent with Julia's naming conventions.

[Marco-Congedo]

What name would be consistent with Julia's convention?

[dlfivefifty]

In Julia functions should be all lower case and normally without capitalization. A reasonable name may be choleskyinv

[Marco-Congedo]

I have been investigating further the phenomenon. After some search and studies, i am now convinced that Julia's cholesky + inv is slower for small matrices because LAPACK's Cholesky method uses BLAS multi-threaded matrix-matrix operations and this is not efficient unless the matrices are large. The problem was found also on LU decomposition (see this post). So, if i force BLAS to use only one thread, Julia's cholesky + inv methods is actually a little bit better than my choleskyinv method. The figures here below show the benchmarks with and without multi-threading.
So, sorry for that, but i am afraid there is no need of the choleskyinv function, provided the user knows that for small matrices it is better to force BLAS to use only one thread!
Should i fork again and delete the add-on of the choleskyinv function?