Dense mat keep alive

CHANGELOG.md

@@ -1,5 +1,12 @@
 # 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.
+- 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
 
 ### Behaviour changes

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

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
 
 """

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)

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)

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)

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