Introduce lu! for UmfpackLU to make use of its symbolic factorization

[goggle]

This PR introduces a lu! method for UmfpackLU to make use of a previously computed symbolic factorization (closes #33323).

Here are some benchmarks:

using LinearAlgebra
using SparseArrays

nSamples = 10000
n = 40
A = sparse(sprandn(n, n, 0.1) + Matrix(10.0*I, n, n))
matrices = [sparse(SparseMatrixCSC(n, n, copy(A.colptr), copy(A.rowval), randn(nnz(A))) + Matrix(10.0*I, n, n)) for i=1:nSamples]

function naive_lu(matrices)
    for A in matrices 
        F = lu(A)
    end
end

function reuse_lu(matrices)
    F = lu(matrices[1])
    for i=2:length(matrices)
        lu!(F, matrices[i])
    end
end

julia> using BenchmarkTools

julia> @btime naive_lu(matrices)
  1.584 s (623885 allocations: 862.68 MiB)

julia> @btime reuse_lu(matrices)
  1.182 s (344645 allocations: 707.05 MiB)

[dkarrasch]

Very nice! This is outside of what I typically work with, so I only have some confusion to share with you.

[goggle]

Here are some more explanation:
The UMFPACK QuickStart Guide explains that there are 5 primary routines for the factorization of a sparse matrix A in the context of A * x = b. Important for us are the routines umfpack_di_symbolic and umfpack_di_numeric. The routine umfpack_di_symbolic performs a symbolic factorization of A. The routine umfpack_di_numeric uses the information from the previously performed symbolic factorization and calculates the numeric factorization of A.

Prior to this PR, the calls to these two routines were coupled: Every time we want to perform a LU factorization of A, we had to first call umfpack_di_symbolic and then umfpack_di_numeric. This PR adds the ability to skip the call for umfpack_di_symbolic if we want to perform a LU factorization of a matrix B with the same sparsity pattern as the matrix A. This is done by modifying the Julia method umfpack_numeric!(U::UmfpackLU), such that it does not return anymore when a numeric factorization is already stored in U, but instead recomputes the numeric factorization without calling umfpack_di_symbolic first. It makes use of the already stored symbolic factorization stored in U and UMFPACK itself handles errors (like different sparsity patterns) by returning a status code.

[dkarrasch]

I see, so the sparsity pattern check is internal in UMFPACK.

[goggle]

@dkarrasch Exactly, yes.

[dkarrasch]

LGTM, except for the minor style comment, but that style is probably not applied elsewhere in this file, so we could as well leave it.

[dkarrasch]

I'm not sure, but is there a place in the Julia docs to show the docstring of the sparse lu!. I couldn't find a page specifically on "sparse linear algebra". It's just when you ask the help, you get all those "sparse" docstrings, but they don't seem to appear in the julia manual.

[ViralBShah]

What's the status on this? Does this need another review before we merge?

[ViralBShah]

Does lu! mutate one of the input arguments? I wonder if a better interface might be to use a named argument which would make it clearer:

lu(A; reuse=F)

[dkarrasch]

Does lu! mutate one of the input arguments?

Yes, it does mutate the first argument F::UmfpackLU IIUC.

[ViralBShah]

In that case is this just a case of rebase and merge?

[dkarrasch]

In that case is this just a case of rebase and merge?

I think so. I wasn't confident enough to simply merge at the time, and was hoping someone else would have a look. But the feature was asked for, so we should let people use and test it in the wild.

[j-fu]

Very nice to see this. This is certainly important for implicit Euler timestepping and path following for PDEs. Will try this out as soon as it becomes available and I can afford the time (unfortunately not before mid February).

[ViralBShah]

Perhaps @ChrisRackauckas and @YingboMa may also want it in DiffEq.

[ChrisRackauckas]

yes I do.

[ViralBShah]

Will be out in 1.5. would be good to test it on master.

[YingboMa]

This is very nice! Thanks for your work!

julia> @btime naive_lu(matrices)
  1.584 s (623885 allocations: 862.68 MiB)

julia> @btime reuse_lu(matrices)
  1.182 s (344645 allocations: 707.05 MiB)

I am surprised that the improvement isn’t very dramatic. Did Julia cache something before?

[ViralBShah]

Symbolic factorization is usually only 10-20% of the factorization time. We already had the capability of reusing the numeric factorization.

We are not caching anything automatically.

[j-fu]

May be the picture changes a bit when setting the pivot tolerance to zero (if you can afford this for your problem, certainly this cannot be the default). The hypothesis is that there would be no reordering steps from partial pivoting, making the numerical factorization relatively cheaper. Not sure though.

[YingboMa]

Thanks for the insight, I will keep that in mind when I implement the feature in JuliaDiffEq.