Skip to content

Factor narrow-band sparse Jacobians with LU and order ILU by RCM - #325

Merged
singhharsh1708 merged 10 commits into
mainfrom
sparse-solver-defaults
Oct 24, 2026
Merged

singhharsh1708 merged 10 commits into
mainfrom
sparse-solver-defaults

Conversation

@singhharsh1708

@singhharsh1708 singhharsh1708 commented Oct 24, 2026 •

Copy link
Copy Markdown
Member

Pick a better default linear solver for serial solves with a sparse jac_prototype when the user names nothing that chooses one. A Jacobian whose bandwidth is at most 4 is factored with LU in the natural ordering, which keeps all fill inside the band. TSRosW and a ksponly solve wrap that LU in a Richardson iteration, whose residual norm catches a NaN right-hand side as GMRES's does on main; Newton solves use preonly. A wider pattern keeps GMRES with ILU(0) but orders it by reverse Cuthill-McKee. Any -ksp_ or -pc_ option except viewers, monitors and converged reasons, or -snes_mf, keeps PETSc's default instead. TSIRK, TSAlpha2, DM and distributed solves are unchanged. solver_options.md and adjoint.md describe the new defaults. Stacked on #323.

On the 1-D Brusselator page problem (500 points, bandwidth 2), from abstol 1e-5 to 1e-9 TSRosW is 1.18 to 1.77 times faster, BDF 1.55 to 2.30 and TSARKIMEX 1.42 to 2.20. On the docs' 2-D Brusselator (32 by 32 periodic grid, n = 2048, alpha = 0.002 N^2, t in [0, 5], the same tolerances) TSRosW is 1.38 to 1.47 times faster, BDF 1.22 to 1.38 and TSARKIMEX 1.12 to 1.21. Errors agree to three digits and all 150 retcodes are Success (Julia 1.12, PETSc 3.22, interleaved runs, figures include #323).

The bandwidth cut comes from a sweep of strip problems: LU in the natural order beat GMRES with ILU(0) at every tolerance up to bandwidth 6 (1.16 to 5.8 times), was mixed at 8 and lost at 12 and 16; the Richardson-wrapped LU that TSRosW uses was 1.09 to 1.84 times faster than GMRES with ILU(0) at bandwidth 2 up to n = 20480. Ordering ILU(0) by RCM was 1.02 to 1.43 times faster on every pattern wider than 4, including the 2-D Brusselator at 64 by 64. MUMPS, UMFPACK, KLU and SuperLU_dist were all slower than either default.

A new testset checks that -ksp_converged_reason, -ksp_monitor and -ksp_view leave a banded solve's steps and states unchanged, and that a NaN right-hand side stops an adaptive ksponly BDF and a TSRosW solve where a dense Jacobian does. Two tests that compared a serial sparse solve with a DMDA or one-rank communicator solve using == now pass -ksp_type gmres to both sides, which keeps the solver both used before.

AI Disclosure: Claude assisted with the code. Testing and review by @singhharsh1708.

@singhharsh1708
singhharsh1708 merged commit fc5cb83 into main Oct 24, 2026
29 of 32 checks passed
singhharsh1708 referenced this pull request Oct 24, 2026
* Read the option names PETSc 3.25 renamed

* Bump to 0.12.3
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

1 participant