9

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.

ali_m
  • 71,714
  • 23
  • 223
  • 298
zonksoft
  • 2,400
  • 1
  • 24
  • 36
  • 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. – Jaime Mar 15 '13 at 22:23

1 Answers1

13

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.

Jaime
  • 65,696
  • 17
  • 124
  • 159
  • 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. – zonksoft Mar 15 '13 at 23:15
  • @RafaelReiter You can actually look at the source code of petsc4py [here](https://code.google.com/p/petsc4py/source/browse/), the function `createAIJ` is [here](https://code.google.com/p/petsc4py/source/browse/src/PETSc/Mat.pyx#242). 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... – Jaime Mar 15 '13 at 23:24
  • 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? – zonksoft Mar 15 '13 at 23:32
  • 1
    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... – zonksoft Mar 15 '13 at 23:45
  • Some one know how to get an parallel solver for a sparse sistem with petsc4py ? – ljofre Jan 19 '15 at 18:43