Passing a PETScPreconditioner to a PETScKrylovSolver segfaults

[Patrick Farrell]

Consider the following code:

###############
from dolfin import *

mesh = UnitSquare(32, 32)
V = FunctionSpace(mesh, 'CG', 1)
bc = DirichletBC(V, Constant(0.0), lambda x, on_boundary: on_boundary)
u = TrialFunction(V); v = TestFunction(V);
A, b = assemble_system( inner(grad(u), grad(v))*dx, Constant(1.0)*v*dx, bc)

KrylovSolver = PETScKrylovSolver("cg", "ilu")
KrylovSolver = PETScKrylovSolver("cg", PETScPreconditioner("ilu")) # comment this out and it works

x_petsc = PETScVector(V.dim())
KrylovSolver.solve(A, x_petsc, down_cast(b))
###############

With the fenics PPA (1.0.0-2~ppa1~natty1), this segfaults:

[pef@aislinn:~/src/dolfin/userpc]$ python petscpc.py
[0]PETSC ERROR: ------------------------------------------------------------------------
[0]PETSC ERROR: Caught signal number 11 SEGV: Segmentation Violation, probably memory access out of range
...

If you comment out the second KrylovSolver constructor, it works.

Taking a look at the code, I suspect the problem is something to do with line 123 of dolfin/la/PETScKrylovSolver.cpp, where the preconditioner of the PETScKrylovSolver is made to point to the user-supplied preconditioner, which disappears from Python as soon as the call is over. In fact, if you run this slightly modified code:

 ###############
from dolfin import *

mesh = UnitSquare(32, 32)
V = FunctionSpace(mesh, 'CG', 1)
bc = DirichletBC(V, Constant(0.0), lambda x, on_boundary: on_boundary)
u = TrialFunction(V); v = TestFunction(V);
A, b = assemble_system( inner(grad(u), grad(v))*dx, Constant(1.0)*v*dx, bc)

p = PETScPreconditioner("ilu")
KrylovSolver = PETScKrylovSolver("cg", p)

x_petsc = PETScVector(V.dim())
KrylovSolver.solve(A, x_petsc, down_cast(b))
 ###############

it works, since p is still around in memory. Unfortunately, I don't know anything about C++/boost/shared_ptrs/swig/whatever to say what line 123 actually /should/ be ...

[Johan Hake (johan-hake)]

Thanks, you were spot on!

These errors tend to be common in the Python interface of DOLFIN previously. However, after we added shared_ptr storage for most data types in Py DOLFIN this has changed. But we still have some untested parts of the interface which still tend to blow up.