From 919de4c4ff76d29e88d7c1e8a959105bcf85674e Mon Sep 17 00:00:00 2001 From: filoferra <66474175+filoferra@users.noreply.github.com> Date: Mon, 19 Oct 2026 13:53:02 +0200 Subject: [PATCH 1/2] Keep dense matrix storage alive while PETSc uses it A dense matrix made from a Julia array kept that array alive only as long as its Julia handle. A KSP holds the matrix past that point, so once the handle was collected the next solve read freed memory: a wrong answer, or a PETSc zero-pivot error. PetscMat(petsclib, A::Matrix) now keeps the array with the PETSc object, as CSR matrices and vectors on Julia arrays already did. PetscMat(petsclib, m, n; type = :dense) lets PETSc allocate the storage, so no Julia array is involved. --- CHANGELOG.md | 6 ++++++ src/mat.jl | 16 +++++++--------- test/core/lifetimes.jl | 27 +++++++++++++++++++++++++++ 3 files changed, 40 insertions(+), 9 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index f5da51a82..ab5b55563 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,5 +1,11 @@ # Changelog +## Unreleased + +### Fixed + +- A dense matrix made with `PetscMat(petsclib, A::Matrix)` or `PetscMat(petsclib, m, n; type = :dense)` could lose its storage while a KSP still used it: once the Julia handle was collected, the array went with it, and the next solve read freed memory (a wrong answer or a PETSc zero-pivot error). The array now lives as long as the PETSc matrix, as it already did for CSR matrices and vectors on Julia arrays. `type = :dense` lets PETSc allocate the storage. + ## v0.5.3 ### Behaviour changes diff --git a/src/mat.jl b/src/mat.jl index 6dc4ed3b4..fdc1f4453 100644 --- a/src/mat.jl +++ b/src/mat.jl @@ -233,14 +233,13 @@ function LibPETSc.PetscMat( ) PetscInt = petsclib.PetscInt - # PETSc stores the data pointer directly without copying, so we must keep the - # backing array alive for the entire lifetime of the PETSc Mat. The finalizer - # closure captures `data`, making it reachable (and therefore uncollectable) for - # as long as `mat` is alive. + # PETSc uses `A`'s memory without copying it, and a KSP can hold the Mat + # after this handle is gone, so the array lives with the PETSc object data = vec(A) mat = LibPETSc.MatCreateSeqDense(petsclib, comm, PetscInt(size(A, 1)), PetscInt(size(A, 2)), data) + keep_alive!(mat, data) - finalizer(m -> (destroy!(m); data), mat) + finalizer(destroy!, mat) return mat end @@ -268,17 +267,16 @@ function LibPETSc.PetscMat( comm = MPI.COMM_SELF check_initialized(petsclib) PetscInt = petsclib.PetscInt - PetscScalar = petsclib.PetscScalar if type === :dense || type === :seqdense - data = zeros(PetscScalar, Int(num_rows) * Int(num_cols)) + # Null data: PETSc allocates zeroed storage and frees it with the Mat mat = LibPETSc.MatCreateSeqDense( petsclib, comm, PetscInt(num_rows), PetscInt(num_cols), - data, + C_NULL, ) - finalizer(m -> (destroy!(m); data), mat) + finalizer(destroy!, mat) return mat elseif type === :seqaij || type === :aij return LibPETSc.PetscMat(petsclib, num_rows, num_cols, 0) diff --git a/test/core/lifetimes.jl b/test/core/lifetimes.jl index 32258436f..fd9c45be8 100644 --- a/test/core/lifetimes.jl +++ b/test/core/lifetimes.jl @@ -54,6 +54,33 @@ end PETSc.destroy!(ksp) end + @testset "dense PetscMat on a Julia array, outliving its handle inside a KSP" begin + # 2I on 3×3, stored in a Julia array the PETSc Mat uses as is + function dense_solver() + A = PETSc.PetscMat(petsclib, PetscScalar[2 0 0; 0 2 0; 0 0 2]) + ksp = PETSc.KSP(A) + PETSc.destroy!(A) # the KSP still holds the matrix + return ksp + end + function empty_dense_solver() + A = PETSc.PetscMat(petsclib, 3, 3; type = :dense) + for i in 1:3 + A[i, i] = 2 + end + PETSc.assemble!(A) + ksp = PETSc.KSP(A) + PETSc.destroy!(A) + return ksp + end + for make_solver in (dense_solver, empty_dense_solver) + ksp = make_solver() + scrub!() + x = ksp \ PetscScalar[2, 4, 6] + @test x ≈ PetscScalar[1, 2, 3] + PETSc.destroy!(ksp) + end + end + @testset "temporaries survive a collection during the operation" begin S = PETSc.PetscMat(petsclib, PetscScalar[2 0 0; 0 2 0; 0 0 2]) ksp = PETSc.KSP(S) From a1d7e73d79261ba56fbd178e870676413fd5d032 Mon Sep 17 00:00:00 2001 From: filoferra <66474175+filoferra@users.noreply.github.com> Date: Mon, 19 Oct 2026 14:03:55 +0200 Subject: [PATCH 2/2] Pop DM constructor options when PETSc throws The DMDA, DMStag and DMPlex constructors and fe_create_default pushed their options onto PETSc's options stack around a PETSc call and popped them afterwards, but not when the call threw. The options then stayed on the stack: every object created afterwards read them, and they could no longer be destroyed. The options are now popped whether or not the call throws. Assisted-by: Claude Opus 5.5 --- CHANGELOG.md | 1 + src/dmda.jl | 6 +++--- src/dmplex.jl | 31 +++++++++++++++---------------- src/dmstag.jl | 12 ++++++------ test/core/options.jl | 29 +++++++++++++++++++++++++++++ 5 files changed, 54 insertions(+), 25 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index ab5b55563..66194f9b7 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -5,6 +5,7 @@ ### Fixed - A dense matrix made with `PetscMat(petsclib, A::Matrix)` or `PetscMat(petsclib, m, n; type = :dense)` could lose its storage while a KSP still used it: once the Julia handle was collected, the array went with it, and the next solve read freed memory (a wrong answer or a PETSc zero-pivot error). The array now lives as long as the PETSc matrix, as it already did for CSR matrices and vectors on Julia arrays. `type = :dense` lets PETSc allocate the storage. +- A `DMDA`, `DMStag` or `DMPlex` constructor, or `fe_create_default`, that threw inside PETSc (for example on `da_refine = "abc"`) left its options pushed on PETSc's options stack. Every object created afterwards read them, and the options could no longer be destroyed. The options are now popped whether or not the call throws. ## v0.5.3 diff --git a/src/dmda.jl b/src/dmda.jl index dc65683e4..af5a7da96 100644 --- a/src/dmda.jl +++ b/src/dmda.jl @@ -153,9 +153,9 @@ function DMDA( if setfromoptions # set options (if any) opts = PetscOptions(petsclib; options...) - push!(opts) - LibPETSc.DMSetFromOptions(PetscLib, da) - pop!(opts) + with_options(opts) do + LibPETSc.DMSetFromOptions(PetscLib, da) + end end if dmsetup diff --git a/src/dmplex.jl b/src/dmplex.jl index 261f88595..e241d884b 100644 --- a/src/dmplex.jl +++ b/src/dmplex.jl @@ -298,9 +298,9 @@ function DMPlex( if setfromoptions opts = PetscOptions(petsclib; options...) - push!(opts) - LibPETSc.DMSetFromOptions(petsclib, dm) - pop!(opts) + with_options(opts) do + LibPETSc.DMSetFromOptions(petsclib, dm) + end end if dmsetup @@ -387,9 +387,9 @@ function DMPlex( if setfromoptions opts = PetscOptions(petsclib; options...) - push!(opts) - LibPETSc.DMSetFromOptions(petsclib, dm) - pop!(opts) + with_options(opts) do + LibPETSc.DMSetFromOptions(petsclib, dm) + end end if dmsetup @@ -537,16 +537,15 @@ function fe_create_default( # picks up the correct space even when no command-line option is given. opt_key = isempty(prefix) ? "petscspace_degree" : "$(prefix)petscspace_degree" opts = PetscOptions(petsclib; Symbol(opt_key) => degree) - push!(opts) - fe = LibPETSc.PetscFECreateDefault( - petsclib, comm, - PetscInt(dim), PetscInt(Nc), - LibPETSc.PetscBool(simplex), - String(prefix), - PetscInt(qorder), - ) - pop!(opts) - return fe + return with_options(opts) do + LibPETSc.PetscFECreateDefault( + petsclib, comm, + PetscInt(dim), PetscInt(Nc), + LibPETSc.PetscBool(simplex), + String(prefix), + PetscInt(qorder), + ) + end end """ diff --git a/src/dmstag.jl b/src/dmstag.jl index 09fcbada0..abdc5763b 100644 --- a/src/dmstag.jl +++ b/src/dmstag.jl @@ -155,9 +155,9 @@ function DMStag( if setfromoptions # set options (if any) opts = PetscOptions(petsclib; options...) - push!(opts) - LibPETSc.DMSetFromOptions(PetscLib, da) - pop!(opts) + with_options(opts) do + LibPETSc.DMSetFromOptions(PetscLib, da) + end end if dmsetup @@ -209,9 +209,9 @@ function DMStag( if setfromoptions opts = PetscOptions(petsclib; options...) - push!(opts) - LibPETSc.DMSetFromOptions(PetscLib, dmnew) - pop!(opts) + with_options(opts) do + LibPETSc.DMSetFromOptions(PetscLib, dmnew) + end end dmsetup && setup!(dmnew) diff --git a/test/core/options.jl b/test/core/options.jl index 24aafbb3e..6bb1de463 100644 --- a/test/core/options.jl +++ b/test/core/options.jl @@ -123,3 +123,32 @@ end @test_throws ArgumentError PETSc.parse_options(["-a=b=c"]) @test PETSc.parse_options(["-pc_type", "mg"]) === (pc_type = "mg",) end + +@testset "a DM constructor that throws leaves the options stack as it was" begin + petsclib = PETSc.petsclibs[1] + PETSc.initialize(petsclib) + comm = PETSc.MPI.COMM_SELF + PetscLib = typeof(petsclib) + nb = PETSc.DM_BOUNDARY_NONE + # With nothing pushed, PETSc reads the global database + global_opts = PETSc.LibPETSc.PetscOptions{PetscLib}(C_NULL, petsclib.age; own = false) + visible(key) = PETSc.LibPETSc.PetscOptionsHasName(petsclib, global_opts, "", key) + + base = PETSc.DMStag(petsclib, comm, (nb,), (10,), (1, 1), 1) + # Each call fails inside PETSc, on an option value that isn't an integer + cases = [ + "-da_refine" => () -> PETSc.DMDA(petsclib, comm, (nb,), (10,), 1, 1; da_refine = "abc"), + "-stag_grid_x" => () -> PETSc.DMStag(petsclib, comm, (nb,), (10,), (1, 1), 1; stag_grid_x = "abc"), + "-stag_grid_x" => () -> PETSc.DMStag(base, (1, 1); stag_grid_x = "abc"), + "-dm_refine" => () -> PETSc.DMPlex(petsclib, comm; dm_refine = "abc"), + "-dm_refine" => () -> PETSc.DMPlex(petsclib, comm, 2, false, [2, 2]; dm_refine = "abc"), + "-petscspace_degree" => () -> PETSc.fe_create_default(petsclib, comm, 2, 1, false; degree = -3), + ] + for (key, construct) in cases + @test_throws PETSc.LibPETSc.PetscError construct() + @test isempty(PETSc.options_stack(PetscLib)) + @test !visible(key) + end + PETSc.destroy!(base) + PETSc.finalize(petsclib) +end