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

How to turn a numpy array to PETScMatrix

0 votes

Hi,
I define a numpy array as

C = np.zeros((numver,numver),dtype=np.float_)

Now, I want to turn it to a PETScMatrix. How could I do it?

Or, how to define a PETScMatrix directly? I tried the following code:

C = Matrix(domain.mpi_comm(),(numver,numver))

but it doesn't work.

asked Apr 19, 2017 by xuanyuan9288 FEniCS Novice (290 points)

1 Answer

+3 votes

Try this:

from petsc4py import PETSc
import numpy as np

A = PETSc.Mat().create()
A.setSizes([10, 10])
A.setType("aij")
A.setUp()

# First arg is list of row indices, second list of column indices
A.setValues([1,2,3], [0,5,9], np.ones((3,3)))
A.assemble()

B = A.convert("dense")
B.getDenseArray()

See:

Note that there is a bug in the automatically generated documentation which affects GenericMatrix.set() and get(). Only one function with three parameters is actually exposed by SWIG and the second and last are actually the indices of the rows and columns to modify.

Inserting new data

If we try to to set() or add() rows/cols of an already initialised sparse PETSc matrix, PETSc will complain with an "out of range" error (code #63) if the insertions don't agree with the sparsity pattern. From the mailing list:

"Preallocation routines now automatically set MAT_NEW_NONZERO_ALLOCATION_ERR,
if you intentionally preallocate less than necessary then use
MatSetOption(mat,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE) to disable the
error generation"

Here are some constants: PETScBool, MatOption, PETSc error codes. And here is how to insert new data:

PETSc_MAT_NEW_NONZERO_ALLOCATION_ERR = 19
PETSc_PETSc_FALSE = 0

# mock example:
A = assemble(a)
#A.set(block, row_indices, col_indices)   # error
A.mat().setOption(PETSc_MAT_NEW_NONZERO_ALLOCATION_ERR, PETSc_FALSE)
A.set(block, row_indices, col_indices)   # all good
A.apply('insert')
answered Apr 20, 2017 by mdbenito FEniCS User (4,530 points)

Thank you very much.

...