SciML / LinearSolve.jl

LinearSolve.jl: High-Performance Unified Interface for Linear Solvers in Julia. Easily switch between factorization and Krylov methods, add preconditioners, and all in one interface.
https://docs.sciml.ai/LinearSolve/stable/
Other
241 stars 52 forks source link

remove `KLU` and `UMFPACK` from the default #456

Closed oscardssmith closed 7 months ago

oscardssmith commented 7 months ago

UMFPACK and KLU both throw errors with singular matrices that are not especially easy to silence ldiv! for sparse UMFPACK and klu do not support passing check=false

Checklist

Additional context

Add any other context about the problem here.

codecov[bot] commented 7 months ago

Codecov Report

Attention: 1 lines in your changes are missing coverage. Please review.

Comparison is base (49183f7) 66.06% compared to head (7df519c) 24.88%.

Files Patch % Lines
src/factorization.jl 0.00% 1 Missing :warning:
Additional details and impacted files ```diff @@ Coverage Diff @@ ## main #456 +/- ## =========================================== - Coverage 66.06% 24.88% -41.19% =========================================== Files 27 27 Lines 2122 2110 -12 =========================================== - Hits 1402 525 -877 - Misses 720 1585 +865 ```

:umbrella: View full report in Codecov by Sentry.
:loudspeaker: Have feedback on the report? Share it here.

ChrisRackauckas commented 7 months ago

Let's be serious...

rayegun commented 7 months ago

Why exactly do these errors need to be suppressed and what is the preferred behavior here? At least KLU I can modify, but I'm somewhat confused.

oscardssmith commented 7 months ago

the sciml style is to return error codes rather than throw errors because that way, higher level interfaces can check the results better. for example, this PR was created because I was getting a SingularException from the nonlinear solve that initializes a DAE solve which would have initialized properly if the linear solve had returned with a failing error code rather than erroring.

rayegun commented 7 months ago

Surely LinSolve should just wrap the solve calls in try catch then? KLU would be sort of non-conformant to the LinearAlgebra factorization "interface" if it returned error codes instead of exceptional errors.

ChrisRackauckas commented 7 months ago

No, try/catch is really slow and is generally not recommended in any numerical code because it introduces dynamic behavior. In general, mathematical functions shouldn't throw but instead report errors. That's why LAPACK's interface has return codes instead of throws, and SuiteSparse does this as well, Sundials, it's why NaNMath.jl exists, etc. all are setup to not throw because that interferes with the recovery behavior of solvers. Julia's interface throws by default, but it's an important enough feature that check=false is standard in LinearAlgebra.jl. It's just missing from the SuiteSparse bindings. They are the odd one out here.

rayegun commented 7 months ago

Alright that's easy enough to fix, I'll try to do that tonight if this flight isn't any further delayed

ChrisRackauckas commented 7 months ago

Thanks!

rayegun commented 7 months ago

UMFPACK already has it :P, I'll check on KLU and then make a PR

rayegun commented 7 months ago

Wait. ldiv! doesn't have a check keyword for any Julia stdlib solvers. Just the factorization call? Is the factorization supposed to store the option, or is this just to disable checking on the factorization call.

ChrisRackauckas commented 7 months ago

None of the operations should error throw. If ldiv! has the ability to throw an error then it should have a way to turn that off, or we can directly do a ccall

rayegun commented 7 months ago

Any errors? Even OOM?

oscardssmith commented 7 months ago

OOM is acceptable (since an OOM isn't really recoverable anyway)

ChrisRackauckas commented 7 months ago

Yeah things that prevent segfaults, sure that needs to error. But it's the errors that are non-essential. Think for example NaNMath.jl. Julia throws an error if it enters the complex plane but that's not necessary, it could just be a NaN as a message that the solution doesn't exist and then a runtime check for NaNs allows optimization routines to continue. This is why JuMP for example has to translate all math functions to NaNMath functions internally.

rayegun commented 7 months ago

Yeah I see at least a few possible throws in LinearAlgebra/cholesky.jl. ldiv! does not respect check = false, but I believe most of these checks guard against malformed decompositions anyway.

The codes that at least Cholesky would throw at ldiv! time are:

  1. CHOLMOD_INVALID: I don't think these should be reachable if the factorization is already valid. Just datatype checking. If you get this either promotion didn't occur correctly (a bug), or it's being called incorrectly.
  2. CHOLMOD_OOM

I'm going to operate under the assumption that these are both acceptable (I imagine similar is true for UMFPACK and KLU). And it's only factorization input issues (singularity, issues with the diagonal, etc) which should be disabled for possible recovery.

rayegun commented 6 months ago

I've added support for check = false to KLU's klu function. However I can't see a path forward to adding this keyword to ldiv!. When I load using LinearAlgebra 0 methods of ldiv! have this keyword. And the contract for ldiv! is to return the result. There is no way for to return an error code unless you want to use the low level calls in this library (which is fine, they have reasonable wrappers, and I can add more check keywords to them if you desire).

These errors must be caught after the factorization and before the solve using issuccess. I'm not 100% sure but I believe ldiv! will not throw numerical errors which should always be found in the factorization call.

ChrisRackauckas commented 6 months ago

Okay, if ldiv! isn't throwing errors then I think we're good?