This is a read only copy of the old FEniCS QA forum. Please visit the new QA forum to ask questions

How to populate a dolfin.PETScMatrix using petsc4py?

+2 votes

I would like to explicitly populate a PETScMatrix from given entries (not via automatic assembly in dolfin). Looking through some older question, the consensus seems to be that it's best to use petsc4py for this. I understand that I can get the underlying petsc4py.PETSc.Mat object as follows:

import dolfin as df
A = df.PETScMatrix()
B = A.mat()  # B is a petsc4py.PETSc.Mat object
print(type(B))

Output:

<type 'petsc4py.PETSc.Mat'>

Now, I can fill the matrix B with values using the regular petsc4py functionality, for example:

B.create()
B.setSizes((3, 3))
B.setType('aij')
B.setUp()
B.setValue(0, 2, 42.0)
B.assemble()
print B.getValues(range(3), range(3))

Output:

[[  0.,   0.,  42.],
 [  0.,   0.,   0.],
 [  0.,   0.,   0.]]

However, these changes don't seem to be propagated into A:

print A.array()

Output:

array([], shape=(0, 0), dtype=float64)

So it seems I'm missing the step to "push" the changes made in B back into A. How can I do this?

Many thanks!

P.S.: In the question at [1] Garth Wells suggested an alternative option to construct a TensorLayout object to specify a sparsity pattern for A. However, there doesn't seem to be a Python interface for TensorLayout. Is this possible to do from Python (and if so, is there example code somewhere), or is it best to use petsc4py directly?

[1] http://fenicsproject.org/qa/2008/how-to-initialize-a-petscmatrix-with-given-sizes?show=2074#a2074

asked Mar 11, 2014 by Maximilian Albert FEniCS User (1,750 points)

To (partially) answer my own question, I just found this bug report which explains why I couldn't just initialise the PETScMatrix with the petsc4py.PETSc.Mat object. So this particular issue should be solved with the next release of dolfin. It would still be interesting to know how to directly modify the underlying petsc4py object after the matrix has been created, though (also to have a workaround until the next release).

1 Answer

+1 vote

I would recommend using the development version of DOLFIN to use petsc4py. The way the PETSc Mat object is wrapped by PETScMatrix has been cleaned up in the DOLFIN dev version and will very likely fix the problem that you're seeing.

answered Mar 12, 2014 by Garth N. Wells FEniCS Expert (35,930 points)

Thanks for your answer, Garth! After discovering the bug report mentioned in my other comment I tried out the development version of dolfin and what works indeed is to initialise a PETScMatrix with a petsc4py.PETSc.Mat in the constructor. However, the code example from my question still doesn't initialise the values of A by setting them in B. Is this expected to work automatically, or am I missing any steps? If the former, should I file a bug report about it?

Go ahead a file a bug report.

...