Skip to content

CSR matrix and docs for PETSc extension #949

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
ChrisRackauckas merged 4 commits into from
Apr 10, 2026
Merged

CSR matrix and docs for PETSc extension #949

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
401299f
first push
Iddingsite Apr 10, 2026
2fb123b
small optim and runic
Iddingsite Apr 10, 2026
0557a96
add docs and more tests
Iddingsite Apr 10, 2026
0133e9c
last addition to docs
Iddingsite Apr 10, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
15 changes: 9 additions & 6 deletions Project.toml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -27,14 +27,14 @@ Setfield = "efcf1570-3423-57d1-acb7-fd33fddbac46"
StaticArraysCore = "1e83bf80-4336-4d27-bf5d-d5a4f845583c"

[weakdeps]
AlgebraicMultigrid = "2169fc97-5a83-5252-b627-83903c6c433c"
ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
AMDGPU = "21141c5a-9bdb-4563-92ae-f87d6854732e"
AlgebraicMultigrid = "2169fc97-5a83-5252-b627-83903c6c433c"
BandedMatrices = "aae01518-5342-5314-be14-df237901396f"
BlockDiagonals = "0a1fb500-61f7-11e9-3c65-f5ef3456f9f0"
CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"
CUDSS = "45b445bb-4962-46a0-9369-b4df9d0f772e"
CUSOLVERRF = "a8cc9031-bad2-4722-94f5-40deabb4245c"
ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
CliqueTrees = "60701a23-6482-424a-84db-faee86b9b1f8"
Elemental = "902c3f28-d1ec-5e7e-8399-a24c3845ee38"
EnzymeCore = "f151be2c-9106-41f4-ab19-57ee4f262869"
Expand All @@ -54,21 +54,22 @@ ParU_jll = "9e0b026c-e8ce-559c-a2c4-6a3d5c955bc9"
Pardiso = "46dd5b70-b6fb-5a00-ae2d-e8fea33afaf2"
RecursiveFactorization = "f2c3362d-daeb-58d1-803e-2bc74f2840b4"
SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
SparseMatricesCSR = "a0a7dd2c-ebf4-11e9-1f05-cf50bc540ca1"
Sparspak = "e56a9233-b9d6-4f03-8d0f-1825330902ac"
blis_jll = "6136c539-28a5-5bf0-87cc-b183200dce32"

[extensions]
LinearSolveAlgebraicMultigridExt = "AlgebraicMultigrid"
LinearSolveAMDGPUExt = "AMDGPU"
LinearSolveAlgebraicMultigridExt = "AlgebraicMultigrid"
LinearSolveBLISExt = ["blis_jll", "LAPACK_jll"]
LinearSolveBandedMatricesExt = "BandedMatrices"
LinearSolveBlockDiagonalsExt = "BlockDiagonals"
LinearSolveCUDAExt = "CUDA"
LinearSolveCUDSSExt = "CUDSS"
LinearSolveCUSOLVERRFExt = ["CUSOLVERRF", "SparseArrays"]
LinearSolveChainRulesCoreExt = "ChainRulesCore"
LinearSolveCliqueTreesExt = ["CliqueTrees", "SparseArrays"]
LinearSolveElementalExt = "Elemental"
LinearSolveChainRulesCoreExt = "ChainRulesCore"
LinearSolveEnzymeExt = ["EnzymeCore", "SparseArrays"]
LinearSolveFastAlmostBandedMatricesExt = "FastAlmostBandedMatrices"
LinearSolveFastLapackInterfaceExt = "FastLapackInterface"
Expand All @@ -80,6 +81,7 @@ LinearSolveKernelAbstractionsExt = "KernelAbstractions"
LinearSolveKrylovKitExt = "KrylovKit"
LinearSolveMetalExt = "Metal"
LinearSolveMooncakeExt = "Mooncake"
LinearSolvePETScCSRExt = ["PETSc", "SparseArrays", "SparseMatricesCSR"]
LinearSolvePETScExt = ["PETSc", "SparseArrays"]
LinearSolveParUExt = ["ParU_jll", "SparseArrays"]
LinearSolvePardisoExt = ["Pardiso", "SparseArrays"]
Expand All @@ -88,8 +90,8 @@ LinearSolveSparseArraysExt = "SparseArrays"
LinearSolveSparspakExt = ["SparseArrays", "Sparspak"]

[compat]
AlgebraicMultigrid = "1"
AMDGPU = "1.2, 2"
AlgebraicMultigrid = "1"
AllocCheck = "0.2"
Aqua = "0.8"
ArrayInterface = "7.19"
Expand Down Expand Up @@ -130,7 +132,7 @@ Metal = "1.4"
Mooncake = "0.5.15"
MultiFloats = "2.3, 3"
OpenBLAS_jll = "0.3"
PETSc = "0.4.6"
PETSc = "0.4.8"
ParU_jll = "1"
Pardiso = "1"
Pkg = "1.10"
Expand All @@ -146,6 +148,7 @@ SciMLLogging = "1.7"
SciMLOperators = "1.15"
Setfield = "1.1.1"
SparseArrays = "1.10"
SparseMatricesCSR = "0.6.12"
Sparspak = "0.3.9"
StableRNGs = "1.0"
StaticArrays = "1.9"
Expand Down
192 changes: 192 additions & 0 deletions docs/src/solvers/solvers.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -375,6 +375,198 @@ GinkgoJL_CG
GinkgoJL_GMRES
```

### PETSc.jl

!!! note

Using this solver requires loading PETSc.jl and MPI.jl, and initialising MPI:
```julia
using PETSc, MPI, SparseMatricesCSR # SparseMatricesCSR is optional but recommended for best performance
MPI.Init()
```

!!! warning "Serial only"

The current implementation supports only **single-process** solves (`MPI.COMM_SELF`).
Passing a multi-rank communicator will raise an error. MPI-parallel support is planned
for a future release.

[PETSc](https://petsc.org) (Portable, Extensible Toolkit for Scientific Computation) is a
library for the parallel numerical solution of scientific applications. Its KSP
component provides a comprehensive suite of Krylov iterative solvers paired with a large
selection of preconditioners.

`PETScAlgorithm` wraps PETSc's KSP solvers via [PETSc.jl](https://github.com/JuliaParallel/PETSc.jl)
and exposes the full preconditioner interface. It works with **dense matrices**,
**`SparseMatrixCSC`**, and **`SparseMatrixCSR`** (from
[SparseMatricesCSR.jl](https://github.com/gridap/SparseMatricesCSR.jl)).

**When to choose PETSc over the pure-Julia Krylov solvers:**
- You want to test a wide variety of Krylov methods and preconditioners without needing to add multiple Julia packages.
- You want direct access to PETSc's Options Database to fine-tune solver behavior at runtime
without recompiling.

#### Solver type

The first positional argument selects the KSP algorithm. Any string accepted by
[`KSPSetType`](https://petsc.org/release/manualpages/KSP/KSPSetType/) can be passed as a `Symbol`.
The most commonly used options are:

| Symbol | Method | Notes |
| :--- | :--- | :--- |
| `:gmres` (default) | GMRES | General non-symmetric systems |
| `:fgmres` | Flexible GMRES | Allows variable preconditioner |
| `:lgmres` | LGMRES | Augmented GMRES, better for restarting |
| `:cg` | Conjugate Gradient | SPD systems only |
| `:fcg` | Flexible CG | CG with variable preconditioner |
| `:minres` | MINRES | Symmetric indefinite systems |
| `:symmlq` | SYMMLQ | Symmetric indefinite systems |
| `:bcgs` | BiCGStab | Non-symmetric, more stable than BiCG |
| `:fbcgs` | Flexible BiCGStab | BiCGStab with variable preconditioner |
| `:bcgsl` | BiCGStab(ℓ) | Stabilised BiCGStab variant |
| `:bicg` | BiConjugate Gradient | Non-symmetric |
| `:cgs` | CGS | Non-symmetric, faster but less stable |
| `:tfqmr` | TFQMR | Transpose-free QMR |
| `:tcqmr` | TCQMR | Transpose-free QMR variant |
| `:cr` | Conjugate Residuals | Symmetric systems |
| `:gcr` | GCR | Generalized CR, flexible preconditioner |
| `:chebyshev` | Chebyshev iteration | Requires eigenvalue bounds; good for smoothing |
| `:richardson` | Richardson iteration | Stationary; mainly used as smoother |
| `:lsqr` | LSQR | Least-squares problems |
| `:cgls` | CGLS | Least-squares problems |
| `:preonly` | Preconditioner only | Use with `:lu` for a direct solve |
| `:none` | No solver | Identity; useful for debugging |

#### Preconditioners

Preconditioners are selected via the `pc_type` keyword. Any string accepted by
[`PCSetType`](https://petsc.org/release/manualpages/PC/PCSetType/) can be passed as a `Symbol`.
The most commonly used options are:

| Symbol | Preconditioner | Notes |
| :--- | :--- | :--- |
| `:none` (default) | No preconditioner | Useful for well-conditioned problems |
| `:jacobi` | Diagonal (Jacobi) scaling | Cheap; good for diagonally dominant systems |
| `:pbjacobi` | Point Block Jacobi | Fixed-size dense blocks along the diagonal |
| `:sor` | SOR / Gauss-Seidel | Successive over-relaxation |
| `:eisenstat` | Eisenstat SSOR | Symmetric SOR; cheaper than a full SSOR sweep |
| `:ilu` | Incomplete LU | General sparse systems |
| `:icc` | Incomplete Cholesky | SPD systems; symmetric analogue of ILU |
| `:lu` | Exact LU (direct) | Use with `:preonly` for a direct solve |
| `:cholesky` | Exact Cholesky (direct) | SPD systems; use with `:preonly` |
| `:bjacobi` | Block Jacobi | Applies an independent ILU/LU solve per block |
| `:asm` | Additive Schwarz | Overlapping domain decomposition |
| `:gasm` | Generalized Additive Schwarz | Multi-level ASM variant |
| `:gamg` | Algebraic Multigrid (GAMG) | No hierarchy needed; good for PDEs |
| `:hypre` | Hypre BoomerAMG | Excellent AMG for large ill-conditioned systems |
| `:kaczmarz` | Kaczmarz | Row-projection smoother |

A separate matrix for building the preconditioner can be supplied via `prec_matrix`:

```julia
PETScAlgorithm(:gmres; prec_matrix = P)
```

#### Matrix format recommendations

PETSc operates internally on 0-based CSR arrays. The recommended matrix format is
**`SparseMatrixCSR{0}`** (from SparseMatricesCSR.jl), which matches PETSc's native layout
exactly:

- **`SparseMatrixCSR{0}`** — *fastest*: zero-copy path on construction; direct `copyto!`
on value-only updates.
- **`SparseMatrixCSR{1}`** — slightly slower than `{0}` on construction (index shift on
cold start), same fast value-update path.
- **`SparseMatrixCSC`** — supported, but requires a CSC→CSR permutation and scatter on
every value update.
- **Dense `Matrix`** — supported via `MatSeqDense`; works out of the box.

```julia
using SparseMatricesCSR, SparseArrays

A_csc = spdiagm(-1 => -ones(n-1), 0 => 2ones(n), 1 => -ones(n-1))

# Recommended: one-liner to build SparseMatrixCSR{0} from a CSC matrix.
# Note: this mutates A_csc's internal storage (colptr/rowvals are shifted in-place).
# Use a copy if you need to keep A_csc usable afterwards.
A = SparseMatrixCSR{0}(transpose(sparse(transpose(A_csc))))
```

#### Basic usage

```julia
using LinearSolve, PETSc, MPI, SparseArrays, LinearAlgebra
MPI.Init()

n = 200
A = sprand(n, n, 0.05); A = A + A' + 20I
b = rand(n)

# Simple one-shot solve with ILU preconditioner
sol = solve(LinearProblem(A, b), PETScAlgorithm(:gmres; pc_type = :ilu))
@show norm(A * sol.u - b) / norm(b)
```

#### Repeated solves (same sparsity pattern, values change)

When the sparsity pattern is fixed across calls,
the KSP is reused and only the matrix values are updated.

```julia
using LinearSolve, PETSc, MPI, SparseArrays, SparseMatricesCSR, LinearAlgebra
import SciMLBase
MPI.Init()

n = 200
A_csc = sprand(n, n, 0.05); A_csc = A_csc + A_csc' + 20I
b = rand(n)

# Convert to SparseMatrixCSR{0} once — getrowptr/getcolval require a CSR matrix
A = SparseMatrixCSR{0}(transpose(sparse(transpose(A_csc))))

cache = SciMLBase.init(LinearProblem(A, b), PETScAlgorithm(:gmres; pc_type = :ilu))
solve!(cache)

# Extract the fixed sparsity structure once (0-based row pointers and column indices)
rowptr0 = copy(SparseMatricesCSR.getrowptr(A))
colval0 = copy(SparseMatricesCSR.getcolval(A))

# Iterate: only nzval changes, sparsity pattern is fixed
for t in 1:10
new_vals = A.nzval .* (1 + 0.1 * t) # e.g. time-varying coefficients
A_new = SparseMatrixCSR{0}(n, n, rowptr0, colval0, new_vals)
SciMLBase.reinit!(cache; A = A_new, b = rand(n))
solve!(cache)
end
```

#### Extra PETSc options

Any PETSc Options Database key can be forwarded via `ksp_options`:

```julia
PETScAlgorithm(:gmres;
pc_type = :ilu,
ksp_options = (ksp_monitor = "", ksp_rtol = 1e-12, pc_factor_levels = 2))
```

#### Memory management

PETSc objects live in C-managed memory outside Julia's GC. Call
`cleanup_petsc_cache!` explicitly when finished to release resources promptly:

```julia
PETScExt = Base.get_extension(LinearSolve, :LinearSolvePETScExt)
PETScExt.cleanup_petsc_cache!(sol) # after solve(...)
PETScExt.cleanup_petsc_cache!(cache) # after init/solve! cycle
```

A GC finalizer is registered as a safety net, but explicit cleanup is strongly preferred.

```@docs
PETScAlgorithm
```

### LinearSolvePyAMG.jl

!!! note
Expand Down
Loading
Loading