Reusing nonzero patterns

Question

I am trying to solve the same PDE repeatedly, with slightly different parameters. This results in the problem $Ax = B$ where $A = C + dD_i$ which I solve for plenty of different parameters $d$ and their corresponding matrix $D_i$.

The nonzero pattern of these $D_i$ is the same, however, and I wanted to use this to my advantage. The nonzero pattern of C and the D's are different. I tried:

# C and all D's have been assembled at this point and are matrices

# initialize the nonzero pattern of A as the combination of C and D
D1 = Ds[0]
A = C.copy()
A.axpy(1.0, D1, False) # `False` here refers to same_nonzero_pattern

for d, D in zip(ds, Ds):
    A.zero()
    A.axpy(1.0, C, True)
    A.axpy(d, D, True)
    solve(A, solution.vector(), B)

At this point I get the error

*** -------------------------------------------------------------------------
*** Error:   Unable to successfully call PETSc function 'KSPSolve'.
*** Reason:  PETSc error code is: 72.
*** Where:   This error was encountered inside <path>/dolfin-1.4.0/dolfin/la/PETScLUSolver.cpp.
*** Process: unknown
*** 
*** DOLFIN version: 1.4.0
*** Git changeset:  
*** -------------------------------------------------------------------------

EDIT
Here is a small example. Did I misunderstand something, because I would expect A and C to be the same at the end.

from dolfin import *

# setup
mesh = UnitSquareMesh(1, 1)
V = FunctionSpace(mesh, "Lagrange", 1)

class LeftBoundary(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], 0.0)
lbnd = LeftBoundary()

class RightBoundary(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], 1.0)
rbnd = RightBoundary()

boundaries = FacetFunction('size_t', mesh)
boundaries.set_all(0)
lbnd.mark(boundaries, 1)
rbnd.mark(boundaries, 2)
ds = Measure('ds')[boundaries]
u = TrialFunction(V)
v = TestFunction(V)

# the interesting part
c = u * v * ds(1)
d = u * v * ds(2)

# dolfin API
C = assemble(c)
D = assemble(d)
A = assemble(c + d)
A.zero()
A.axpy(1.0, C, True) # wrong
print A.array(), '\n'

# petsc4py API
Cp = as_backend_type(C).mat()
Dp = as_backend_type(D).mat()
Ap = as_backend_type(A).mat()
Ap.zeroEntries()
# Ap.axpy(1.0, Cp, structure=0) # correct
Ap.axpy(1.0, Cp, structure=1) # also correct but faster
# Ap.axpy(1.0, Cp, structure=2) # wrong, same as `A` above
print PETScMatrix(Ap).array(), '\n'

print C.array()

From PETSc documentation ( http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Mat/MatAXPY.html ) the structure can be

either SAME_NONZERO_PATTERN, DIFFERENT_NONZERO_PATTERN or
SUBSET_NONZERO_PATTERN (nonzeros of X is a subset of Y's)

I believe that this corresponds to the structure= keyword argument of the petsc4py wrapping, but I haven't found which corresponds to which, neither why one doesn't work.

Comments

Error code 72 is floating point exception and may be unrelated to the matrix axpy operations. You can check that matrices are correct before solving (perform matrix-vector products or look at the numpy arrays for small problem sizes).

---

I provided a small working example that shows the problem. Either I misunderstand the method axpy or it sets the wrong nonzero elements (hence leaving zeros where the solver does not expect them, which leads to the error). Does that make sense?

Another way to create a matrix with a given nonzero pattern seems to be by using a TensorLayout object (see also the PS at the end of this question http://fenicsproject.org/qa/2951/how-to-populate-a-dolfin-petscmatrix-using-petsc4py?show=2951) but I did not find an example of how to use it?

A last option would be to use petsc4py to create a matrix with a sparsity pattern and use that one to initialize a dolfin matrix. I assume this must be possible somehow.

---

I updated the working example to show the differences between dolfin and petsc4py. It seems that I can use a workaround using petsc4py, but some parts still remain a mystery to me.
If there is a better way I would be happy to hear.

---

A and C does not have the same sparsity. The nonzero entries of C is a proper subset of the nonzero entries of A, which is why structure=1 can be used here, but not structure=2.

You can also write

Ap.axpy(1.0, Cp, structure = Ap.Structure.SUBSET_NONZERO_PATTERN)

---

Thanks, that seems more clear! So if I understand things now, same_nonzero_pattern means exactly the same; but since dolfin only offers exactly the same or different, I have to use petsc4py with the subset_nonzero_pattern keyword.

Then, from the discussion above, I assume that structure=0 corresponds to different_nonzero_pattern (correct but slow), structure=1 to subset_nonzero_pattern and structure=2 to same_nonzero_pattern?