[Enhancement] Enable computation of main diagonal of the inverse from factorizations

[MothNik]

Requiring the main diagonal of the inverse is not uncommon when solving pentadiagonal systems, e.g., during Cross-Validation for spline smoothing.
Similar to #28 this should be easily doable from the factorization A=LU where we have:

Since the main diagonal of L is known, this can also be written as A=LDU where

Now, special formulas can be applied for computing selected main diagonal elements of the inverse inv(A) as taken from the master thesis

Fasllija, E., Parallel computation of the diagonal of the inverse of a sparse matrix, 2017, pp. 8-9

All of this applies in our case and it would be a minimum effort to add a function that computes this from a factorization.
Given that one can solely rely on U, the fact that L's lower triangle is dense should not pose a problem.