Fix Radau methods to support KLU with sparse jacobian prototypes #2893
+115
−12
Conversation
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
- Added is_sparse and nonzeros imports to OrdinaryDiffEqFIRK - Added SparseArrays as dependency - Modified RadauIIA3, RadauIIA5, RadauIIA9, and AdaptiveRadau cache initialization to preserve sparsity pattern from jacobian prototype for complex matrices - Added test cases to reproduce the issue Addresses SciML#2892 where KLUFactorization fails with jacobian prototypes. The fix preserves the sparsity structure but further work may be needed to ensure KLU can properly initialize with the complex sparse matrices. Issue: Complex sparse matrices initialized with zeros prevent KLU from creating proper factorization cache, resulting in 'type Nothing has no field colptr' error.
ChrisRackauckas-Claude
changed the title
| J, W1_temp = build_J_W(alg, u, uprev, p, t, dt, f, jac_config, uEltypeNoUnits, Val(true)) | ||
| # For sparse matrices, preserve sparsity pattern for KLU compatibility | ||
| if is_sparse(J) | ||
| W1 = similar(J, Complex{eltype(W1_temp)}) |
There was a problem hiding this comment.
this line is the same between branches.
Comment on lines
+72
to
+73
| SparseConnectivityTracer = "9f842d2f-2579-4b1d-911e-f412cf18a3f5" | ||
| Sparspak = "e56a9233-b9d6-4f03-8d0f-1825330902ac" |
There was a problem hiding this comment.
these should be in extras.
|
I think the real answer is that we should generate the W via |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
3 participants
Add this suggestion to a batch that can be applied as a single commit.
This suggestion is invalid because no changes were made to the code.
Suggestions cannot be applied while the pull request is closed.
Suggestions cannot be applied while viewing a subset of changes.
Only one suggestion per line can be applied in a batch.
Add this suggestion to a batch that can be applied as a single commit.
Applying suggestions on deleted lines is not supported.
You must change the existing code in this line in order to create a valid suggestion.
Outdated suggestions cannot be applied.
This suggestion has been applied or marked resolved.
Suggestions cannot be applied from pending reviews.
Suggestions cannot be applied on multi-line comments.
Suggestions cannot be applied while the pull request is queued to merge.
Suggestion cannot be applied right now. Please check back later.
Fixes #2892
Summary
This PR fixes an issue where Radau methods (RadauIIA3, RadauIIA5, RadauIIA9, and AdaptiveRadau) would fail when using
KLUFactorizationwith a sparsejac_prototype.Root Cause
The issue was caused by a limitation in LinearSolve's KLU interface: KLU doesn't properly support complex sparse matrices. When initializing a LinearSolve cache with a complex sparse matrix and KLU, the factorization is never created (remains
Nothing), leading to aFieldError: type Nothing has no field colptrwhen trying to solve.Solution
This PR implements a two-part fix:
Proper initialization of complex W matrices: Use
f.jac_prototype .* (1 + 0im)(when available) to create complex W matrices with the correct sparsity pattern from the user-provided prototype.Automatic fallback to UMFPACK: When KLU is specified as the linear solver, automatically use
UMFPACKFactorizationfor the complex W matrices (W2/W3), while still using KLU for the real W1 matrix. UMFPACK fully supports complex sparse matrices.Changes
lib/OrdinaryDiffEqFIRK/src/firk_caches.jl: Updated cache initialization for all four Radau methods to use jac_prototype and fallback to UMFPACK for complex matriceslib/OrdinaryDiffEqFIRK/src/OrdinaryDiffEqFIRK.jl: Added imports forKLUFactorizationandUMFPACKFactorizationTesting
Verified that all Radau methods now work correctly with:
KLUFactorizationwith sparsejac_prototype✓LUFactorization✓SparspakFactorization✓The fix preserves the sparsity pattern correctly and allows users to benefit from KLU's performance for real matrices while seamlessly handling the complex matrix requirements of Radau methods.
🤖 Generated with Claude Code
Co-Authored-By: Claude [email protected]