search

JuliaSparse / SparseArrays.jl

SparseArrays.jl is a Julia stdlib
https://sparsearrays.juliasparse.org/
Other
90 stars 51 forks source link

Pure-Julia implementation of LU for sparse matrices #568

Closed johnomotani closed 1 week ago

johnomotani commented 1 week ago

It would be nice (?) to have a pure-Julia implementation of LU factorisation and solve for sparse matrices.

I made a start on this (because I want to parallelise ldiv!()) here https://github.com/johnomotani/SharedMemSparseLU.jl/tree/serial-version There's a serial version of ldiv!(), but it relies on SparseArrays.jl/UMFPACK for the factorisation. ldiv!() turned out to be pretty simple in the end - the whole package is only 262 lines - although I didn't worry about trying to use BLAS or anything else fancy. I realise that the factorisation is the hard part (!), but even just having ldiv!() gives options, e.g. adding threads, etc., if you want to reuse the factorisation many times and care about scaling.

I have no plans to work on factorisation myself, but if anyone wants to adopt that code, then go ahead! I'm also very much not an expert on numerical methods, so if anyone wants to look at the code and see if I've done anything stupid, input would be very welcome :slightly_smiling_face:

j-fu commented 1 week ago

Just to make you aware: there is Sparspak.jl which has a pure Julia implementation path for number types besides Float64 and Float32. However its structures had been implemented without parallelization in mind. Also there is Pardiso.jl which is calls the multithreaded Pardiso solver (but restricted to Float32/Float64). A Julia implementation from scratch would have to compete with these sophisticated packages, and also UMFPACK, so the "hard part" is indeed a quite nontrivial endeavour.

johnomotani commented 1 week ago

Thanks @j-fu! I think Sparspak.jl was what I was looking for to start with (i.e. pure-Julia LU). Guess I should've asked for help earlier!