Add LU factorization and ldiv!

[ctessum]

In response to SciML/LinearSolve.jl#575

[mipals]

Hey Christopher,

I do not know if you have realized but this package is not registered. The reason for this is that I in the end did not need block diagonals for my application and I stopped the development.

However, if the current features together with your changes suits your needs I would be happy to include them as well as registering the package.

Cheers,
Mikkel

[ChrisRackauckas]

If this is useful, can we move to JuliaLinearAlgebra and reigster+maintain it? It's worth testing first @ctessum

[ctessum]

Thanks, I realized it wasn't registered after I tried to actually use it 😅.

My case has square blocks that are all the same size, which makes fast indexing much easier (which is the reason why this package was created in the first place AFAIK), so I just ended up creating a special implementation for that case and putting it here: https://github.com/EarthSciML/EarthSciMLBase.jl/blob/main/src/blockdiagonal.jl

@ChrisRackauckas there is a registered package here but the indexing is too slow for it to be useful, and it doesn't seem like anyone is responding to issues/PRs. But I'd be happy to work my code for the special case into that one (or this one), just let me know what would be useful.

[mipals]

That’s correct. I made this package because I believe that the current implementation in BlockDiagonals.jl is flawed as it does not include the sizes of the blocks in the BlockDiagonal struct. This flaw makes e.g. indexing a O(n) operation and mul/div serial as you always have to loop through the blocks from the start in order to get the index of the blocks (in the case of equally sized blocks this can be avoided, but in general it is not the case). I did try to get this fixed but the issue i created seemed to go stale.

If you want I would be happy to transfer this package. I would also be up for doing a little maintaining/getting it up to date if needed.

[ChrisRackauckas]

Yeah let's get this revived and in JuliaLinearAlgebra. It seems like we already found a use case for it.

[mipals]

Cool. As said I would be happy to help in the revitalization of the package.

I don't know the procedure of getting a package into an organization. Is just that a fork is created in JuliaLinearAlgebra?

[ctessum]

Sounds good! To make things work for my use case I ended up making some changes that I temporarily put here: https://github.com/EarthSciML/EarthSciMLBase.jl/blob/main/src/blockdiagonal.jl

That version:

• Only handles the special case where all of the blocks are the same size. This allows zero-time-order indexing rather than the (i think) linear time indexing here. My version assumes all the blocks are square, but it would also work if they were same-size rectangles.
• Allows the user to specify a strategy for applying operations to blocks. The currently available options are to use broadcasting or threads.

I'd be happy to contribute some of those changes here (or in the new fork) if they are deemed useful.