scipy sparse matrices as an input for petsc4py

Question

I can't seem to find a way how to efficiently load scipy sparse matrices, e.g. csr_matrix, into a petsc4py matrix, e.g. PETSc.Mat().createAIJ. I found this thread, but I'm not able to apply it.

I would also appreciate a pointer where this stuff is actually documented. The examples in the demo directory only explain a part, and I can't see any docstrings.

Phyllisphylloclade

What exactly have you tried, and how has it not worked? I don't have petsc4py on my system, but the instructions in your link seem to me pretty straightforward to follow. – Buckskin

Answer 1

Your link says that to create a sparse matrix in PETSc, you should use a command like this:

PETSc.Mat().createAIJ(size=(nrows,ncols), csr=(ai,aj,aa))

According to this, the ai, aj and aa are, in PETSc-speak:

> i - row indices
> j - column indices
> a - matrix values

These are equivalent, respectively, to the .indptr, .indices and .data attributes of a scypy.sparse.csr_matrix, see the docs for details.

So, if your link is right, the following should work:

>>> from petsc4py import PETSc
>>> import scipy.sparse
>>> csr_mat = scipy.sparse.rand(1000, 1000, density=0.001, format='csr')
>>> petsc_mat = PETSc.Mat().createAIJ(size=csr_mat.shape,
...                                   csr=(csr_mat.indptr, csr_mat.indices,
...                                        csr_mat.data))

Unfortunately, I cannot test it myself.

Buckskin

This completely answers my question. You also lifted my confusion about documentation because the Mat().createAIJ() function obviously calls the PETSc function (easily readable code) which is indeed documented. I just didn't realize that this function call is the direct call of the PETSc function. – Phyllisphylloclade

@RafaelReiter You can actually look at the source code of petsc4py here, the function createAIJ is here. It seems to be a little bit more complicated than that, a Mat_AllocAIJ_CSR which I haven't been able to find gets called, so I am not really sure what goes on under the hood. But if this works for you, then it works for you... – Buckskin

You're right, it's not that simple. I thought I'd seen MatCreateSeqAIJWithArrays a few hours ago in the code, but I didn't. Do you have any hint how to come up with the right way to use all petsc4py functionality beside the demos? – Phyllisphylloclade

I found the call hierarchy: Mat_AllocAIJ_CSR is in src/PETSc/petscmat.pxi, and it calls MatAnyAIJSetPreallocationCSR in src/include/custom.h, which calls various PETSc functions, e.g. MatMPIAIJSetPreallocationCSR which is actually documented in the PETSc documentation. Phew... – Phyllisphylloclade

Some one know how to get an parallel solver for a sparse sistem with petsc4py ? – Duckweed