From 401299fd55e54a3d44ed8504b7ddea20f8b4386d Mon Sep 17 00:00:00 2001 From: Iddingsite Date: Fri, 10 Apr 2026 11:46:35 +0200 Subject: [PATCH 1/4] first push --- ext/LinearSolvePETScCSRExt.jl | 124 ++++++++++++++++++++++++++++ ext/LinearSolvePETScExt.jl | 78 ++++++++++-------- test/petsctests.jl | 149 ++++++++++++++++++++++++++++++++++ 3 files changed, 317 insertions(+), 34 deletions(-) create mode 100644 ext/LinearSolvePETScCSRExt.jl diff --git a/ext/LinearSolvePETScCSRExt.jl b/ext/LinearSolvePETScCSRExt.jl new file mode 100644 index 000000000..c2fda39eb --- /dev/null +++ b/ext/LinearSolvePETScCSRExt.jl @@ -0,0 +1,124 @@ +module LinearSolvePETScCSRExt + +using SparseMatricesCSR: SparseMatrixCSR, getrowptr, getcolval +using PETSc +import PETSc: MPI, LibPETSc +using LinearSolve: PETScAlgorithm, LinearCache, LinearSolve + +const PETScExt = Base.get_extension(LinearSolve, :LinearSolvePETScExt) + +# ── Sparse trait ───────────────────────────────────────────────────────────── + +PETScExt.is_sparse_petsc(::SparseMatrixCSR) = true + +# ── Canonical-type conversion ──────────────────────────────────────────────── +# +# PETSc's MatCreateSeqAIJWithArrays requires 0-based Int and the correct scalar +# type. Following GridapPETSc.jl we convert A to SparseMatrixCSR{0,PetscScalar, +# PetscInt} (the "canonical" form) before handing arrays to PETSc. +# +# Cost breakdown: +# Bi=0, types match → zero allocation; the original CSR is returned as-is. +# Bi=0, index type differs → type-cast rowptr/colval only (no shift). +# Bi≠0 (e.g. default Bi=1) → shift every index by -Bi and type-cast if needed. +# +# In all cases nzval is copied so that the user's Julia matrix is never mutated +# by the Case-3 in-place value update. + +function _to_petsc_canonical( + A::SparseMatrixCSR{Bi, Tv, Ti}, + ::Type{PetscInt}, ::Type{PetscScalar}) where {Bi, Tv, Ti, PetscInt, PetscScalar} + m, n = size(A) + src_rp = getrowptr(A) + src_cv = getcolval(A) + src_nz = A.nzval + + # ── index arrays ───────────────────────────────────────────────────────── + if Bi == 0 + if Ti == PetscInt + # Perfect-match path: borrow the original arrays — zero allocation. + Tv == PetscScalar && return A + rp = src_rp + cv = src_cv + else + rp = PetscInt.(src_rp) + cv = PetscInt.(src_cv) + end + else + rp = PetscInt[v - Bi for v in src_rp] + cv = PetscInt[v - Bi for v in src_cv] + end + + # ── value array (always copied to avoid mutating the user's matrix) ────── + nz = Tv == PetscScalar ? copy(src_nz) : PetscScalar.(src_nz) + return SparseMatrixCSR{0}(m, n, rp, cv, nz) +end + +# ── Matrix construction ─────────────────────────────────────────────────────── +# +# Following GridapPETSc.jl: +# 1. Convert A to the canonical SparseMatrixCSR{0,PetscScalar,PetscInt}. +# 2. Pass its arrays directly to MatCreateSeqAIJWithArrays (zero-copy for +# indices when types already match). +# 3. Store the canonical CSR object — not just the arrays — as the GC anchor +# so that PETSc's borrowed pointers remain valid for the lifetime of PA. + +function PETScExt.to_petsc_mat(petsclib, A::SparseMatrixCSR{Bi}) where {Bi} + PetscInt = PETSc.inttype(petsclib) + PetscScalar = PETSc.scalartype(petsclib) + + canonical = _to_petsc_canonical(A, PetscInt, PetscScalar) + m, n = size(canonical) + + mat = LibPETSc.MatCreateSeqAIJWithArrays( + petsclib, MPI.COMM_SELF, + PetscInt(m), PetscInt(n), + getrowptr(canonical), getcolval(canonical), canonical.nzval) + + # Store the whole canonical CSR as GC anchor (GridapPETSc pattern). + # PETSc.destroy will pop! this automatically when the mat is destroyed. + PETSc._MATSEQAIJ_WITHARRAYS_STORAGE[mat.ptr] = canonical + return mat +end + +# ── Sparsity pattern storage ───────────────────────────────────────────────── +# +# CSR is already in PETSc's native row-major order — no permutation needed. +# Reuse prev_colptr / prev_rowval to hold rowptr / colval for change detection. + +function PETScExt.store_sparse_pattern!(pcache, A::SparseMatrixCSR) + pcache.sparse_perm = nothing + pcache.sparse_scratch = nothing + pcache.prev_colptr = getrowptr(A) + pcache.prev_rowval = getcolval(A) +end + +# ── Sparsity pattern change detection ──────────────────────────────────────── + +function PETScExt.check_pattern_changed(pcache, A::SparseMatrixCSR) + pcache.prev_colptr === nothing && return true + old_rp, old_cv = pcache.prev_colptr, pcache.prev_rowval + new_rp, new_cv = getrowptr(A), getcolval(A) + (length(old_rp) != length(new_rp) || length(old_cv) != length(new_cv)) && return true + return old_rp != new_rp || old_cv != new_cv +end + +# ── In-place value update (Case 3) ─────────────────────────────────────────── +# +# CSR nzval is already in PETSc's row-major order, so a direct copyto! suffices +# — no scatter/permutation needed. MatSeqAIJGetArray returns a pointer to the +# buffer that was passed to MatCreateSeqAIJWithArrays (the canonical nzval). + +function PETScExt.update_sparse_values!( + petsclib, PA, pcache, A::SparseMatrixCSR; assemble::Bool = true) + vals = LibPETSc.MatSeqAIJGetArray(petsclib, PA) + try + copyto!(vals, A.nzval) + finally + LibPETSc.MatSeqAIJRestoreArray(petsclib, PA, vals) + end + assemble && PETSc.assemble!(PA) + return nothing +end + +end diff --git a/ext/LinearSolvePETScExt.jl b/ext/LinearSolvePETScExt.jl index b8257a1f8..30604c296 100644 --- a/ext/LinearSolvePETScExt.jl +++ b/ext/LinearSolvePETScExt.jl @@ -159,15 +159,19 @@ end # ── Matrix conversion ───────────────────────────────────────────────────────── +# Trait function: is A a sparse format that PETSc can handle as an AIJ matrix? +# Overloaded by LinearSolvePETScCSRExt to also return true for SparseMatrixCSR. +is_sparse_petsc(A) = A isa SparseMatrixCSC + # Convert a Julia matrix to a PETSc sequential matrix. # # Sparse (SparseMatrixCSC) → MatSeqAIJWithArrays : -# Internally converts CSR format, so quite inefficient to construct from CSC. -# +# Internally converts to CSR format. # Dense (AbstractMatrix) → MatSeqDense: # PETSc allocates its own buffer and copies the values in on construction. # This protects the Julia array even when a factorising preconditioner is # used without a separate `prec_matrix`. +# SparseMatrixCSR → handled by LinearSolvePETScCSRExt via to_petsc_mat overload. function to_petsc_mat(petsclib, A) if A isa SparseMatrixCSC return PETSc.MatSeqAIJWithArrays(petsclib, MPI.COMM_SELF, A) @@ -262,6 +266,33 @@ function update_mat_values!(petsclib, PA, A::AbstractMatrix) return PETSc.assemble!(PA) end +# ── Extensible sparse dispatch ─────────────────────────────────────────────── +# +# These functions are overloaded by LinearSolvePETScCSRExt to add CSR support. +# They abstract over the sparse format so build_ksp! and solve! don't need +# hardcoded `isa SparseMatrixCSC` checks. + +function store_sparse_pattern!(pcache, A::SparseMatrixCSC) + nnz_val = length(A.nzval) + m = size(A, 1) + if pcache.sparse_perm === nothing || length(pcache.sparse_perm) != nnz_val + pcache.sparse_perm = Vector{Int}(undef, nnz_val) + pcache.sparse_scratch = Vector{Int}(undef, 2 * (m + 1)) + end + _csc_to_csr_perm!(pcache.sparse_perm, pcache.sparse_scratch, A) + pcache.prev_colptr = A.colptr + pcache.prev_rowval = A.rowval +end + +function check_pattern_changed(pcache, A::SparseMatrixCSC) + pcache.prev_colptr === nothing && return true + return sparsity_pattern_changed(pcache.prev_colptr, pcache.prev_rowval, A) +end + +function update_sparse_values!(petsclib, PA, pcache, A::SparseMatrixCSC; assemble::Bool = true) + return update_mat_values!(petsclib, PA, A, pcache.sparse_perm; assemble = assemble) +end + # ── Vector I/O ──────────────────────────────────────────────────────────────── # # VecGetArray / VecGetArrayRead return Julia Arrays wrapping PETSc's internal @@ -350,9 +381,8 @@ end # optionally attach a null space. Called on Case 1 (first solve) and Case 2 # (structure changed). # -# For sparse matrices, the CSC→CSR permutation is computed here once and stored -# in pcache for reuse on all subsequent Case 3 value-only updates. The buffers -# are only reallocated when nnz changes (which implies a Case 2 rebuild anyway). +# For sparse matrices, the sparsity pattern and any permutation data are stored +# in pcache for reuse on all subsequent Case 3 value-only updates. # # vec_n is reset to 0 so ensure_vecs! unconditionally recreates petsc_x and # petsc_b after every rebuild — required when the problem size changes. @@ -360,18 +390,8 @@ function build_ksp!(pcache, petsclib, cache, alg) pcache.vec_n = 0 pcache.petsc_A = to_petsc_mat(petsclib, cache.A) - if cache.A isa SparseMatrixCSC - nnz = length(cache.A.nzval) - m = size(cache.A, 1) - # Reallocate only when nnz changes (implies a structural change). - if pcache.sparse_perm === nothing || length(pcache.sparse_perm) != nnz - pcache.sparse_perm = Vector{Int}(undef, nnz) - pcache.sparse_scratch = Vector{Int}(undef, 2 * (m + 1)) - end - _csc_to_csr_perm!(pcache.sparse_perm, pcache.sparse_scratch, cache.A) - # Store the current sparsity pattern - pcache.prev_colptr = cache.A.colptr - pcache.prev_rowval = cache.A.rowval + if is_sparse_petsc(cache.A) + store_sparse_pattern!(pcache, cache.A) else pcache.sparse_perm = nothing pcache.sparse_scratch = nothing @@ -440,17 +460,8 @@ function SciMLBase.solve!(cache::LinearCache, alg::PETScAlgorithm; kwargs...) # ── Decide: rebuild, update in-place, or reuse ──────────────────────────── rebuild_ksp = pcache.ksp === nothing # Case 1 if !rebuild_ksp && cache.isfresh - if cache.A isa SparseMatrixCSC - if pcache.prev_colptr === nothing || pcache.prev_rowval === nothing - rebuild_ksp = true - else - # Compare stored pointers/arrays directly against the new matrix - rebuild_ksp = sparsity_pattern_changed( - pcache.prev_colptr, - pcache.prev_rowval, - cache.A - ) - end + if is_sparse_petsc(cache.A) + rebuild_ksp = check_pattern_changed(pcache, cache.A) else # For dense matrices, check if the size has changed rebuild_ksp = pcache.prev_size !== nothing && size(cache.A) != pcache.prev_size @@ -461,18 +472,17 @@ function SciMLBase.solve!(cache::LinearCache, alg::PETScAlgorithm; kwargs...) build_ksp!(pcache, petsclib, cache, alg) elseif cache.isfresh # Case 3: same structure — update values without touching the KSP object. - if cache.A isa SparseMatrixCSC - update_mat_values!( - petsclib, pcache.petsc_A, cache.A, pcache.sparse_perm; + if is_sparse_petsc(cache.A) + update_sparse_values!( + petsclib, pcache.petsc_A, pcache, cache.A; assemble = cache.precsisfresh ) else update_mat_values!(petsclib, pcache.petsc_A, cache.A) end if alg.prec_matrix !== nothing - if alg.prec_matrix isa SparseMatrixCSC - # Preconditioner shares the pattern with A — reuse the permutation. - update_mat_values!(petsclib, pcache.petsc_P, alg.prec_matrix, pcache.sparse_perm) + if is_sparse_petsc(alg.prec_matrix) + update_sparse_values!(petsclib, pcache.petsc_P, pcache, alg.prec_matrix) else update_mat_values!(petsclib, pcache.petsc_P, alg.prec_matrix) end diff --git a/test/petsctests.jl b/test/petsctests.jl index c1e82f2f5..8f8e8552a 100644 --- a/test/petsctests.jl +++ b/test/petsctests.jl @@ -3,6 +3,7 @@ using LinearSolve using MPI using PETSc using SparseArrays +using SparseMatricesCSR using Random using Test @@ -502,3 +503,151 @@ end PETScExt.cleanup_petsc_cache!(sol) end end + +# ══════════════════════════════════════════════════════════════════════════════ +# CSR MATRIX TESTS (SparseMatrixCSR via SparseMatricesCSR.jl) +# ══════════════════════════════════════════════════════════════════════════════ + +@testset "Serial: CSR Basic Solvers" begin + n = 100 + A_csc = sprand(n, n, 0.05) + 10I + A_csc = A_csc'A_csc + A_csr = SparseMatrixCSR(A_csc) + b = rand(n) + + @testset "GMRES (CSR)" begin + sol = solve( + LinearProblem(A_csr, b), PETScAlgorithm(:gmres); abstol = 1.0e-10, reltol = 1.0e-10 + ) + @test norm(A_csc * sol.u - b) / norm(b) < 1.0e-6 + PETScExt.cleanup_petsc_cache!(sol) + end + + @testset "CG (CSR)" begin + sol = solve( + LinearProblem(A_csr, b), PETScAlgorithm(:cg); abstol = 1.0e-10, reltol = 1.0e-10 + ) + @test norm(A_csc * sol.u - b) / norm(b) < 1.0e-6 + PETScExt.cleanup_petsc_cache!(sol) + end + + @testset "GMRES + ILU (CSR)" begin + sol = solve( + LinearProblem(A_csr, b), PETScAlgorithm(:gmres; pc_type = :ilu); + abstol = 1.0e-10, reltol = 1.0e-10 + ) + @test norm(A_csc * sol.u - b) / norm(b) < 1.0e-6 + PETScExt.cleanup_petsc_cache!(sol) + end +end + +@testset "Serial: CSR Cache Interface" begin + n = 100 + A_csc = sprand(n, n, 0.05) + 10I; A_csc = A_csc'A_csc + A_csr = SparseMatrixCSR(A_csc) + b = rand(n) + + cache = SciMLBase.init( + LinearProblem(A_csr, b), PETScAlgorithm(:gmres); abstol = 1.0e-10, reltol = 1.0e-10 + ) + sol1 = solve!(cache) + @test norm(A_csc * sol1.u - b) / norm(b) < 1.0e-6 + + b2 = rand(n) + cache.b = b2 + sol2 = solve!(cache) + @test norm(A_csc * sol2.u - b2) / norm(b2) < 1.0e-6 + PETScExt.cleanup_petsc_cache!(cache) +end + +@testset "Serial: CSR Matrix Rebuild Logic" begin + n = 50 + Random.seed!(123) + + base = sprand(n, n, 0.1); base = base + base' + 10I + pattern_I, pattern_J, _ = findnz(base) + + function make_csr_spd_with_pattern(scale) + vals = abs.(randn(length(pattern_I))) .* scale + A = sparse(pattern_I, pattern_J, vals, n, n) + A = A + A' + d = vec(sum(abs.(A), dims = 2)) .+ scale + for i in 1:n + A[i, i] = d[i] + end + return SparseMatrixCSR(A) + end + + @testset "Case 1: first solve builds KSP (CSR)" begin + A = make_csr_spd_with_pattern(10.0) + A_dense = Matrix(A) + b = rand(n) + cache = SciMLBase.init( + LinearProblem(A, b), PETScAlgorithm(:cg; pc_type = :jacobi); abstol = 1.0e-10 + ) + @test cache.cacheval.ksp === nothing + sol = solve!(cache) + @test cache.cacheval.ksp !== nothing + @test norm(A_dense * sol.u - b) / norm(b) < 1.0e-6 + PETScExt.cleanup_petsc_cache!(cache) + end + + @testset "Case 3: same sparsity — KSP reused, values updated (CSR)" begin + A1 = make_csr_spd_with_pattern(10.0) + b1 = rand(n) + cache = SciMLBase.init( + LinearProblem(A1, b1), PETScAlgorithm(:cg; pc_type = :jacobi); abstol = 1.0e-10 + ) + solve!(cache) + ksp_before = cache.cacheval.ksp + + A2 = make_csr_spd_with_pattern(20.0) + A2_dense = Matrix(A2) + b2 = rand(n) + SciMLBase.reinit!(cache; A = A2, b = b2) + sol2 = solve!(cache) + + @test cache.cacheval.ksp === ksp_before # KSP reused + @test norm(A2_dense * sol2.u - b2) / norm(b2) < 1.0e-6 + PETScExt.cleanup_petsc_cache!(cache) + end + + @testset "Case 2: sparsity changed — KSP rebuilt (CSR)" begin + A1 = make_csr_spd_with_pattern(10.0) + b1 = rand(n) + cache = SciMLBase.init( + LinearProblem(A1, b1), PETScAlgorithm(:cg; pc_type = :jacobi); abstol = 1.0e-10 + ) + solve!(cache) + ksp_before = cache.cacheval.ksp + + A3_csc = sprand(n, n, 0.2) + 15I; A3_csc = A3_csc + A3_csc' + A3 = SparseMatrixCSR(A3_csc) + A3_dense = Matrix(A3) + b3 = rand(n) + SciMLBase.reinit!(cache; A = A3, b = b3) + sol3 = solve!(cache) + + @test cache.cacheval.ksp !== ksp_before # KSP rebuilt + @test norm(A3_dense * sol3.u - b3) / norm(b3) < 1.0e-6 + PETScExt.cleanup_petsc_cache!(cache) + end + + @testset "Case 3: correctness over multiple value-only updates (CSR)" begin + A = make_csr_spd_with_pattern(10.0) + b = rand(n) + cache = SciMLBase.init( + LinearProblem(A, b), PETScAlgorithm(:cg; pc_type = :jacobi); abstol = 1.0e-10 + ) + + for i in 1:5 + A_new = make_csr_spd_with_pattern(10.0 + i) + A_new_dense = Matrix(A_new) + b_new = rand(n) + SciMLBase.reinit!(cache; A = A_new, b = b_new) + sol = solve!(cache) + @test norm(A_new_dense * sol.u - b_new) / norm(b_new) < 1.0e-6 + end + PETScExt.cleanup_petsc_cache!(cache) + end +end From 2fb123be241aa043c27a221b10c1b803791d128b Mon Sep 17 00:00:00 2001 From: Iddingsite Date: Fri, 10 Apr 2026 14:22:10 +0200 Subject: [PATCH 2/4] small optim and runic --- Project.toml | 15 ++++--- ext/LinearSolvePETScCSRExt.jl | 54 +++++++++++++---------- ext/LinearSolvePETScExt.jl | 82 +++++++++++++++++------------------ test/petsc/Project.toml | 1 + 4 files changed, 81 insertions(+), 71 deletions(-) diff --git a/Project.toml b/Project.toml index 5b6b3eaee..6277beb2d 100644 --- a/Project.toml +++ b/Project.toml @@ -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" @@ -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" @@ -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"] @@ -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" @@ -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" @@ -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" diff --git a/ext/LinearSolvePETScCSRExt.jl b/ext/LinearSolvePETScCSRExt.jl index c2fda39eb..875fe47b9 100644 --- a/ext/LinearSolvePETScCSRExt.jl +++ b/ext/LinearSolvePETScCSRExt.jl @@ -13,25 +13,27 @@ PETScExt.is_sparse_petsc(::SparseMatrixCSR) = true # ── Canonical-type conversion ──────────────────────────────────────────────── # -# PETSc's MatCreateSeqAIJWithArrays requires 0-based Int and the correct scalar -# type. Following GridapPETSc.jl we convert A to SparseMatrixCSR{0,PetscScalar, -# PetscInt} (the "canonical" form) before handing arrays to PETSc. +# PETSc's MatCreateSeqAIJWithArrays requires 0-based indices and the correct +# scalar type. Convert A to SparseMatrixCSR{0,PetscScalar,PetscInt} (the +# "canonical" form) before handing arrays to PETSc. # # Cost breakdown: # Bi=0, types match → zero allocation; the original CSR is returned as-is. # Bi=0, index type differs → type-cast rowptr/colval only (no shift). # Bi≠0 (e.g. default Bi=1) → shift every index by -Bi and type-cast if needed. # -# In all cases nzval is copied so that the user's Julia matrix is never mutated -# by the Case-3 in-place value update. +# nzval is borrowed (zero-copy) on the perfect-match path (Bi=0, types match); +# otherwise it is copied so that the user's Julia matrix is never mutated by +# the Case-3 in-place value update. function _to_petsc_canonical( A::SparseMatrixCSR{Bi, Tv, Ti}, - ::Type{PetscInt}, ::Type{PetscScalar}) where {Bi, Tv, Ti, PetscInt, PetscScalar} - m, n = size(A) - src_rp = getrowptr(A) - src_cv = getcolval(A) - src_nz = A.nzval + ::Type{PetscInt}, ::Type{PetscScalar} + ) where {Bi, Tv, Ti, PetscInt, PetscScalar} + m, n = size(A) + src_rp = getrowptr(A) + src_cv = getcolval(A) + src_nz = A.nzval # ── index arrays ───────────────────────────────────────────────────────── if Bi == 0 @@ -56,24 +58,24 @@ end # ── Matrix construction ─────────────────────────────────────────────────────── # -# Following GridapPETSc.jl: -# 1. Convert A to the canonical SparseMatrixCSR{0,PetscScalar,PetscInt}. -# 2. Pass its arrays directly to MatCreateSeqAIJWithArrays (zero-copy for -# indices when types already match). -# 3. Store the canonical CSR object — not just the arrays — as the GC anchor -# so that PETSc's borrowed pointers remain valid for the lifetime of PA. +# 1. Convert A to the canonical SparseMatrixCSR{0,PetscScalar,PetscInt}. +# 2. Pass its arrays directly to MatCreateSeqAIJWithArrays (zero-copy for +# indices when types already match). +# 3. Store the canonical CSR object — not just the arrays — as the GC anchor +# so that PETSc's borrowed pointers remain valid for the lifetime of PA. function PETScExt.to_petsc_mat(petsclib, A::SparseMatrixCSR{Bi}) where {Bi} - PetscInt = PETSc.inttype(petsclib) + PetscInt = PETSc.inttype(petsclib) PetscScalar = PETSc.scalartype(petsclib) canonical = _to_petsc_canonical(A, PetscInt, PetscScalar) - m, n = size(canonical) + m, n = size(canonical) mat = LibPETSc.MatCreateSeqAIJWithArrays( petsclib, MPI.COMM_SELF, PetscInt(m), PetscInt(n), - getrowptr(canonical), getcolval(canonical), canonical.nzval) + getrowptr(canonical), getcolval(canonical), canonical.nzval + ) # Store the whole canonical CSR as GC anchor (GridapPETSc pattern). # PETSc.destroy will pop! this automatically when the mat is destroyed. @@ -87,10 +89,10 @@ end # Reuse prev_colptr / prev_rowval to hold rowptr / colval for change detection. function PETScExt.store_sparse_pattern!(pcache, A::SparseMatrixCSR) - pcache.sparse_perm = nothing + pcache.sparse_perm = nothing pcache.sparse_scratch = nothing - pcache.prev_colptr = getrowptr(A) - pcache.prev_rowval = getcolval(A) + pcache.prev_colptr = getrowptr(A) + return pcache.prev_rowval = getcolval(A) end # ── Sparsity pattern change detection ──────────────────────────────────────── @@ -99,6 +101,8 @@ function PETScExt.check_pattern_changed(pcache, A::SparseMatrixCSR) pcache.prev_colptr === nothing && return true old_rp, old_cv = pcache.prev_colptr, pcache.prev_rowval new_rp, new_cv = getrowptr(A), getcolval(A) + # Cheap identity check first — same arrays means pattern is unchanged. + old_rp === new_rp && old_cv === new_cv && return false (length(old_rp) != length(new_rp) || length(old_cv) != length(new_cv)) && return true return old_rp != new_rp || old_cv != new_cv end @@ -110,10 +114,12 @@ end # buffer that was passed to MatCreateSeqAIJWithArrays (the canonical nzval). function PETScExt.update_sparse_values!( - petsclib, PA, pcache, A::SparseMatrixCSR; assemble::Bool = true) + petsclib, PA, pcache, A::SparseMatrixCSR; assemble::Bool = true + ) vals = LibPETSc.MatSeqAIJGetArray(petsclib, PA) try - copyto!(vals, A.nzval) + # Skip the copy when PETSc already owns A.nzval (zero-copy fast path). + pointer(vals) != pointer(A.nzval) && copyto!(vals, A.nzval) finally LibPETSc.MatSeqAIJRestoreArray(petsclib, PA, vals) end diff --git a/ext/LinearSolvePETScExt.jl b/ext/LinearSolvePETScExt.jl index 30604c296..285cbedba 100644 --- a/ext/LinearSolvePETScExt.jl +++ b/ext/LinearSolvePETScExt.jl @@ -281,7 +281,7 @@ function store_sparse_pattern!(pcache, A::SparseMatrixCSC) end _csc_to_csr_perm!(pcache.sparse_perm, pcache.sparse_scratch, A) pcache.prev_colptr = A.colptr - pcache.prev_rowval = A.rowval + return pcache.prev_rowval = A.rowval end function check_pattern_changed(pcache, A::SparseMatrixCSC) @@ -293,49 +293,47 @@ function update_sparse_values!(petsclib, PA, pcache, A::SparseMatrixCSC; assembl return update_mat_values!(petsclib, PA, A, pcache.sparse_perm; assemble = assemble) end -# ── Vector I/O ──────────────────────────────────────────────────────────────── +# ── Vector I/O (VecPlaceArray pattern) ─────────────────────────────────────── # -# VecGetArray / VecGetArrayRead return Julia Arrays wrapping PETSc's internal -# memory — no extra allocation. The matching Restore call is mandatory and -# must always be reached (hence the try/finally pattern). +# Use VecPlaceArray / VecResetArray to temporarily redirect PETSc's internal +# Vec buffer to point at Julia arrays — zero copy in both directions. # -# Note: VecAssemblyBegin/End is NOT needed after direct pointer writes. -# Assembly is only required after VecSetValues, which uses an off-process -# communication stash. Direct pointer access bypasses the stash entirely. +# VecPlaceArray(petsc_b, cache.b) # redirect pointer — no copy +# VecPlaceArray(petsc_x, cache.u) # redirect pointer — initial guess in-place +# KSPSolve(...) # writes solution directly into cache.u +# VecResetArray(petsc_x/petsc_b) # restore PETSc's own buffer +# +# This avoids 3 × O(n) memcpy per solve (write b, write x, read x back). +# VecPlaceArray requires eltype == PetscScalar, which is guaranteed because +# petsclib is chosen by eltype(cache.A) and u/b share that type. +# VecResetArray is always called in a finally block. + +# Allocate a PETSc sequential Vec of length n. +# Contents are uninitialised; VecPlaceArray redirects its buffer before each use. +function create_seq_vec(petsclib, n::Int) + return LibPETSc.VecCreateSeq(petsclib, MPI.COMM_SELF, petsclib.PetscInt(n)) +end -function write_vec!(pv, src::AbstractVector, petsclib) +# Allocate and fill a Vec from a Julia vector (used for nullspace basis vectors). +function create_seq_vec(petsclib, src::AbstractVector) + pv = create_seq_vec(petsclib, length(src)) ptr = LibPETSc.VecGetArray(petsclib, pv) - return try - copyto!(ptr, src) # Use vectorized copy + try + copyto!(ptr, src) finally LibPETSc.VecRestoreArray(petsclib, pv, ptr) end -end - -function read_vec!(dst::AbstractVector, pv, petsclib) - ptr = LibPETSc.VecGetArrayRead(petsclib, pv) - return try - copyto!(dst, ptr) - finally - LibPETSc.VecRestoreArrayRead(petsclib, pv, ptr) - end -end - -# Allocate a PETSc sequential Vec of length n and fill it from src. -function create_seq_vec(petsclib, src::AbstractVector) - pv = LibPETSc.VecCreateSeq(petsclib, MPI.COMM_SELF, petsclib.PetscInt(length(src))) - write_vec!(pv, src, petsclib) return pv end # Recreate petsc_x and petsc_b only when the problem size changes. -function ensure_vecs!(pcache, petsclib, u, b) +function ensure_vecs!(pcache, petsclib, b) n = length(b) if pcache.vec_n != n || pcache.petsc_x === nothing || pcache.petsc_b === nothing pcache.petsc_x !== nothing && PETSc.destroy(pcache.petsc_x) pcache.petsc_b !== nothing && PETSc.destroy(pcache.petsc_b) - pcache.petsc_x = create_seq_vec(petsclib, u) - pcache.petsc_b = create_seq_vec(petsclib, b) + pcache.petsc_x = create_seq_vec(petsclib, n) + pcache.petsc_b = create_seq_vec(petsclib, n) pcache.vec_n = n end return pcache.petsc_x, pcache.petsc_b @@ -497,20 +495,22 @@ function SciMLBase.solve!(cache::LinearCache, alg::PETScAlgorithm; kwargs...) end cache.isfresh = false - # ── Vectors ─────────────────────────────────────────────────────────────── - petsc_x, petsc_b = ensure_vecs!(pcache, petsclib, cache.u, cache.b) - write_vec!(petsc_x, cache.u, petsclib) - write_vec!(petsc_b, cache.b, petsclib) - - # ── Solve ───────────────────────────────────────────────────────────────── - if alg.transposed - LibPETSc.KSPSolveTranspose(petsclib, pcache.ksp, petsc_b, petsc_x) - else - LibPETSc.KSPSolve(petsclib, pcache.ksp, petsc_b, petsc_x) + # ── Vectors + Solve (VecPlaceArray — zero copy) ─────────────────────────── + petsc_x, petsc_b = ensure_vecs!(pcache, petsclib, cache.b) + LibPETSc.VecPlaceArray(petsclib, petsc_b, cache.b) + LibPETSc.VecPlaceArray(petsclib, petsc_x, cache.u) + try + if alg.transposed + LibPETSc.KSPSolveTranspose(petsclib, pcache.ksp, petsc_b, petsc_x) + else + LibPETSc.KSPSolve(petsclib, pcache.ksp, petsc_b, petsc_x) + end + # Solution is written directly into cache.u — no read-back needed. + finally + LibPETSc.VecResetArray(petsclib, petsc_x) + LibPETSc.VecResetArray(petsclib, petsc_b) end - read_vec!(cache.u, petsc_x, petsclib) - # ── Convergence metadata ────────────────────────────────────────────────── iters = Int(LibPETSc.KSPGetIterationNumber(petsclib, pcache.ksp)) reason = Int(LibPETSc.KSPGetConvergedReason(petsclib, pcache.ksp)) diff --git a/test/petsc/Project.toml b/test/petsc/Project.toml index da266566d..18973d03e 100644 --- a/test/petsc/Project.toml +++ b/test/petsc/Project.toml @@ -4,4 +4,5 @@ LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae" MPI = "da04e1cc-30fd-572f-bb4f-1f8673147195" PETSc = "ace2c81b-2b5f-4b1e-a30d-d662738edfe0" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" +SparseMatricesCSR = "a0a7dd2c-ebf4-11e9-1f05-cf50bc540ca1" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" From 0557a960681df720bd0f6e4a57edbc58195db78f Mon Sep 17 00:00:00 2001 From: Iddingsite Date: Fri, 10 Apr 2026 15:16:36 +0200 Subject: [PATCH 3/4] add docs and more tests --- docs/src/solvers/solvers.md | 181 ++++++++++++++++++++++++++++++++ test/petsctests.jl | 203 ++++++++++++++++++++++++++++++++++++ 2 files changed, 384 insertions(+) diff --git a/docs/src/solvers/solvers.md b/docs/src/solvers/solvers.md index 91a63f26f..47619d8f9 100644 --- a/docs/src/solvers/solvers.md +++ b/docs/src/solvers/solvers.md @@ -375,6 +375,187 @@ 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. + +`PETScAlgorithm` wraps PETSc's KSP (Krylov subspace) solvers and exposes the full PETSc +preconditioner interface using [PETSc.jl](https://github.com/JuliaParallel/PETSc.jl). It works with **dense matrices**, **`SparseMatrixCSC`**, and +**`SparseMatrixCSR`** (from [SparseMatricesCSR.jl](https://github.com/gridap/SparseMatricesCSR.jl)). + +#### 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 diff --git a/test/petsctests.jl b/test/petsctests.jl index 8f8e8552a..d0d958cba 100644 --- a/test/petsctests.jl +++ b/test/petsctests.jl @@ -4,6 +4,7 @@ using MPI using PETSc using SparseArrays using SparseMatricesCSR +import SparseMatricesCSR: getrowptr, getcolval using Random using Test @@ -651,3 +652,205 @@ end PETScExt.cleanup_petsc_cache!(cache) end end + +# ══════════════════════════════════════════════════════════════════════════════ +# CSR{0} MATRIX TESTS (SparseMatrixCSR{0} — 0-based indices, zero-copy path) +# SparseMatrixCSR{0} matches PETSc's native CSR layout exactly: +# • Cold start — no index shift, arrays borrowed directly. +# • Value update — direct copyto! into PETSc's buffer (no permutation). +# ══════════════════════════════════════════════════════════════════════════════ + +# Helper: build a SparseMatrixCSR{0} from a CSC matrix. +# The 1→0 shift is paid once here; all subsequent PETSc calls are zero-copy. +function csr0_from_csc(A_csc::SparseMatrixCSC) + csr1 = SparseMatrixCSR(A_csc) + rowptr = getrowptr(csr1) .- one(eltype(getrowptr(csr1))) + colval = getcolval(csr1) .- one(eltype(getcolval(csr1))) + return SparseMatrixCSR{0}(size(A_csc, 1), size(A_csc, 2), rowptr, colval, copy(csr1.nzval)) +end + +@testset "Serial: CSR{0} Basic Solvers" begin + n = 100 + A_csc = sprand(n, n, 0.05) + 10I; A_csc = A_csc'A_csc + A_csr0 = csr0_from_csc(A_csc) + b = rand(n) + + @testset "GMRES (CSR{0})" begin + sol = solve( + LinearProblem(A_csr0, b), PETScAlgorithm(:gmres); abstol = 1.0e-10, reltol = 1.0e-10 + ) + @test norm(A_csc * sol.u - b) / norm(b) < 1.0e-6 + PETScExt.cleanup_petsc_cache!(sol) + end + + @testset "CG (CSR{0})" begin + sol = solve( + LinearProblem(A_csr0, b), PETScAlgorithm(:cg); abstol = 1.0e-10, reltol = 1.0e-10 + ) + @test norm(A_csc * sol.u - b) / norm(b) < 1.0e-6 + PETScExt.cleanup_petsc_cache!(sol) + end + + @testset "GMRES + ILU (CSR{0})" begin + sol = solve( + LinearProblem(A_csr0, b), PETScAlgorithm(:gmres; pc_type = :ilu); + abstol = 1.0e-10, reltol = 1.0e-10 + ) + @test norm(A_csc * sol.u - b) / norm(b) < 1.0e-6 + PETScExt.cleanup_petsc_cache!(sol) + end +end + +@testset "Serial: CSR{0} Cache Interface" begin + n = 100 + A_csc = sprand(n, n, 0.05) + 10I; A_csc = A_csc'A_csc + A_csr0 = csr0_from_csc(A_csc) + b = rand(n) + + cache = SciMLBase.init( + LinearProblem(A_csr0, b), PETScAlgorithm(:gmres); abstol = 1.0e-10, reltol = 1.0e-10 + ) + sol1 = solve!(cache) + @test norm(A_csc * sol1.u - b) / norm(b) < 1.0e-6 + + b2 = rand(n) + cache.b = b2 + sol2 = solve!(cache) + @test norm(A_csc * sol2.u - b2) / norm(b2) < 1.0e-6 + PETScExt.cleanup_petsc_cache!(cache) +end + +@testset "Serial: CSR{0} Matrix Rebuild Logic" begin + n = 50 + Random.seed!(456) + + base = sprand(n, n, 0.1); base = base + base' + 10I + pattern_I, pattern_J, _ = findnz(base) + + # Shared 0-based index arrays — same sparsity pattern for all Case-3 solves. + base_csr1 = SparseMatrixCSR(base) + rowptr0 = getrowptr(base_csr1) .- 1 + colval0 = getcolval(base_csr1) .- 1 + + function make_csr0_spd_with_pattern(scale) + vals = abs.(randn(length(pattern_I))) .* scale + A = sparse(pattern_I, pattern_J, vals, n, n); A = A + A' + d = vec(sum(abs.(A), dims = 2)) .+ scale + for i in 1:n + A[i, i] = d[i] + end + csr1 = SparseMatrixCSR(A) + # Reuse the shared rowptr0/colval0 only when pattern matches; + # otherwise build fresh 0-based arrays (required for Case-2 test). + rp = getrowptr(csr1) .- 1 + cv = getcolval(csr1) .- 1 + return SparseMatrixCSR{0}(n, n, rp, cv, copy(csr1.nzval)) + end + + @testset "Case 1: first solve builds KSP (CSR{0})" begin + A = make_csr0_spd_with_pattern(10.0) + b = rand(n) + cache = SciMLBase.init( + LinearProblem(A, b), PETScAlgorithm(:cg; pc_type = :jacobi); abstol = 1.0e-10 + ) + @test cache.cacheval.ksp === nothing + sol = solve!(cache) + @test cache.cacheval.ksp !== nothing + @test norm(Matrix(A) * sol.u - b) / norm(b) < 1.0e-6 + PETScExt.cleanup_petsc_cache!(cache) + end + + @testset "Case 3: same sparsity — KSP reused, values updated (CSR{0})" begin + # make_csr0_spd_with_pattern builds valid SPD matrices with the same sparsity + # pattern (same rowptr/colval *values*), so check_pattern_changed returns false + # via element comparison and the KSP is reused. + A1 = make_csr0_spd_with_pattern(10.0) + b1 = rand(n) + cache = SciMLBase.init( + LinearProblem(A1, b1), PETScAlgorithm(:cg; pc_type = :jacobi); abstol = 1.0e-10 + ) + solve!(cache) + ksp_before = cache.cacheval.ksp + + A2 = make_csr0_spd_with_pattern(20.0) + b2 = rand(n) + SciMLBase.reinit!(cache; A = A2, b = b2) + sol2 = solve!(cache) + + @test cache.cacheval.ksp === ksp_before # KSP reused + @test norm(Matrix(A2) * sol2.u - b2) / norm(b2) < 1.0e-6 + PETScExt.cleanup_petsc_cache!(cache) + end + + @testset "Case 2: sparsity changed — KSP rebuilt (CSR{0})" begin + A1 = make_csr0_spd_with_pattern(10.0) + b1 = rand(n) + cache = SciMLBase.init( + LinearProblem(A1, b1), PETScAlgorithm(:cg; pc_type = :jacobi); abstol = 1.0e-10 + ) + solve!(cache) + ksp_before = cache.cacheval.ksp + + A3_csc = sprand(n, n, 0.2) + 15I; A3_csc = A3_csc + A3_csc' + A3 = csr0_from_csc(A3_csc) + b3 = rand(n) + SciMLBase.reinit!(cache; A = A3, b = b3) + sol3 = solve!(cache) + + @test cache.cacheval.ksp !== ksp_before # KSP rebuilt + @test norm(A3_csc * sol3.u - b3) / norm(b3) < 1.0e-6 + PETScExt.cleanup_petsc_cache!(cache) + end + + @testset "Case 3: correctness over multiple value-only updates (CSR{0})" begin + A = make_csr0_spd_with_pattern(10.0) + b = rand(n) + cache = SciMLBase.init( + LinearProblem(A, b), PETScAlgorithm(:cg; pc_type = :jacobi); abstol = 1.0e-10 + ) + solve!(cache) + + for i in 1:5 + A_new = make_csr0_spd_with_pattern(10.0 + i) + b_new = rand(n) + SciMLBase.reinit!(cache; A = A_new, b = b_new) + sol = solve!(cache) + @test norm(Matrix(A_new) * sol.u - b_new) / norm(b_new) < 1.0e-6 + end + PETScExt.cleanup_petsc_cache!(cache) + end +end + +@testset "Serial: CSR{0} zero-copy fast path" begin + # When Bi=0, Ti==PetscInt, Tv==PetscScalar the canonical conversion returns + # A itself — no allocation. Verify the backing arrays are not copied by + # checking that pointer identity is preserved after construction and after + # a Case-3 value-only update. + n = 50 + A_csc = sprand(n, n, 0.1) + 10I; A_csc = A_csc'A_csc + A0 = csr0_from_csc(A_csc) + b = rand(n) + + ptr_rp = pointer(getrowptr(A0)) + ptr_cv = pointer(getcolval(A0)) + ptr_nz = pointer(A0.nzval) + + cache = SciMLBase.init( + LinearProblem(A0, b), PETScAlgorithm(:gmres); abstol = 1.0e-10 + ) + sol = solve!(cache) + @test norm(A_csc * sol.u - b) / norm(b) < 1.0e-6 + + # Pointers must be unchanged — no hidden copy of the index or value arrays. + @test pointer(getrowptr(A0)) === ptr_rp + @test pointer(getcolval(A0)) === ptr_cv + @test pointer(A0.nzval) === ptr_nz + + # Case-3 update: scale A0's nzval — reuses the same index arrays and preserves SPD. + A0_new = SparseMatrixCSR{0}(n, n, getrowptr(A0), getcolval(A0), A0.nzval .* 1.5) + SciMLBase.reinit!(cache; A = A0_new, b = rand(n)) + sol2 = solve!(cache) + @test norm(Matrix(A0_new) * sol2.u - cache.b) / norm(cache.b) < 1.0e-6 + + PETScExt.cleanup_petsc_cache!(cache) +end From 0133e9cf46eb67e0db652957efab8acca5e8d13b Mon Sep 17 00:00:00 2001 From: Iddingsite Date: Fri, 10 Apr 2026 15:27:42 +0200 Subject: [PATCH 4/4] last addition to docs --- docs/src/solvers/solvers.md | 17 ++++++++++++++--- 1 file changed, 14 insertions(+), 3 deletions(-) diff --git a/docs/src/solvers/solvers.md b/docs/src/solvers/solvers.md index 47619d8f9..ee9d54aa1 100644 --- a/docs/src/solvers/solvers.md +++ b/docs/src/solvers/solvers.md @@ -391,9 +391,20 @@ GinkgoJL_GMRES Passing a multi-rank communicator will raise an error. MPI-parallel support is planned for a future release. -`PETScAlgorithm` wraps PETSc's KSP (Krylov subspace) solvers and exposes the full PETSc -preconditioner interface using [PETSc.jl](https://github.com/JuliaParallel/PETSc.jl). It works with **dense matrices**, **`SparseMatrixCSC`**, and -**`SparseMatrixCSR`** (from [SparseMatricesCSR.jl](https://github.com/gridap/SparseMatricesCSR.jl)). +[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