[petsc-users] Ordering of preallocation and OwnershipRange

[Matthew Knepley]

On Sat, Aug 13, 2011 at 2:57 PM, Josh Hykes <jmhykes at ncsu.edu> wrote:

> Hello,
>
> I'm just starting to experiment with PETSc (v3.1), and I like the Python
> bindings provided by petsc4py (v1.1.2). So far things seem fairly
> straightforward, but I'm stumped on a small issue.
>
> While creating a parallel AIJ matrix, I'd like to preallocate it using
> arrays d_nnz and o_nnz. As I understand it, these arrays correspond to the
> processor's local rows.
>
> Currently I specify the global matrix size, and let PETSc decide on the
> decomposition of the rows. I'd like to ask PETSc what rows each processor
> has with the getOwnershipRange() function, and then do the preallocation.
> However, according to the error message
>
> > [1] MatAnyAIJSetPreallocation() line 393 in
> petsc4py-1.1.2/src/include/custom.h
> > [1] Operation done in wrong order
> > [1] matrix is already preallocated
>
> I'm not allowed to do it in this order.
>
> Thus, my question is: is it possible to let PETSc figure out the row
> decomposition while still using d_nnz and o_nnz for the preallocation? I
> figure that I could resolve the problem by doing my own decomposition, but
> it'd be nice if I could let those details up to PETSc.
>

You are correct. We require that preallocation is done at the same time as
decomposition. There
are tricky dependencies in matrix creation. However, an easy workaround is
to create a Vec at
the same time with the same global size, since it is guaranteed to have the
same layout. I will look
into simplifying this if it is possible.

  Thanks,

     Matt


> I'm including an example using petsc4py of what I'd like to do, run with 2
> MPI processes.
>
> I apologize if this is a dumb question. Thank you for your help.
>
> -Josh
>
> # -----------------------------------------------
> from petsc4py import PETSc as petsc
>
> M, N = 4, 6
>
> global_d_nnz = [2, 1, 1, 2]
> global_o_nnz = [1, 3, 2, 1]
>
>  A = petsc.Mat()
> A.create(petsc.COMM_WORLD)
> A.setSizes([M, N])
> A.setType('aij')
>
> i_start, i_end = A.getOwnershipRange()
>
> A.setPreallocationNNZ([global_d_nnz[i_start:i_end],
>                        global_o_nnz[i_start:i_end]]) # error occurs here
>
>


-- 
What most experimenters take for granted before they begin their experiments
is infinitely more interesting than any results to which their experiments
lead.
-- Norbert Wiener
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20110813/50716e6c/attachment.htm>