[petsc-users] Ordering with MatCreateNest and MPI

[Jed Brown]

Alessandro Francavilla <matteo.francavilla at polito.it> writes:

> Thanks for the reply; probably I've not been clear enough. I don't need to
> separate out components, I just need to understand the underlying orderings.
>
> For instance, suppose that I want to multiply the matrix with a null vector
> with a 1 in the i-th location (the result is basically the i-th column of
> the matrix):
>
> CALL
> MatCreateNest(PETSC_COMM_WORLD,2,PETSC_NULL_OBJECT,2,PETSC_NULL_OBJECT,matAr
> ray,A,IERR)
> CALL MatCreateVecs(A,x,b,IERR)
> CALL VecSet(x,0.+0.*PETSC_i,IERR)
> CALL VecSetValues(x,1,i-1,(1.+0.*PETSC_i),INSERT_VALUES,IERR)
> CALL VecAssemblyBegin(x,IERR)
> CALL VecAssemblyEnd(x,IERR)
> CALL MatMult(A,x,b,IERR)
>
> The above code works perfectly if I execute it single process (mpirun -np
> 1); if I run it on >1 MPI processes, the result is not what I would expect:
> rather than returning the i-th column, it returns (e.g., if I call VecView
> on b) some permutation of a different column (say a permutation of column
> j). 
>
> Also, if I substitute the call to MatCreateNest() with
> A = matArray(1)

But this is a different matrix.  Why would you expect that the MatNest
has the same indices as its first Mat?

> Then everything works as expected, for arbitrary numbers of processes. By
> working as expected I mean that calling VecView on b prints the i-th column
> of A with its natural ordering. 
>
> So, my question is:
> 1) how do I generate the input vector x such that the MatMult operation
> returns the i-th column of A?

It is returning the ith column of A.  A is not just matArray(1).

> 2) how do I reorder the result of MatMult such that I get the i-th column
> with its "natural" ordering?
>
> Maybe I'm even misunderstanding the output of VecView, which prints the
> elements process by process, and thus the order corresponds to the natural
> ordering only if the elements are distributed consecutively among processes?
> Which would explain why I'm seeing a permutation of the column. However, I'm
> seeing the permutation of a DIFFERENT column of the matrix, which means that
> the input vector does not have a 1 in the i-th position when I fill it with
> VecSetValues(x,1,i-1,(1.+0.*PETSC_i),INSERT_VALUES,IERR). 

PETSc Vecs are distributed with contiguous blocks assigned to each
process.  Mat's column and row indices are induced by their action on
Vecs.  MatNest does not move matrix data around, so its blocks produce
the portions of Vecs that are owned on the current process.
Consequently, MatNest induces an ordering in which one field gets the
first chunk of a contiguous space on each process, with the intervening
entries for the other field(s).  If you have two matrices, B and C, each
distributed over the same two processes, then the MatNest

  [B 0; 0 C]

will keep both B and C distributed.  I.e., it does not move all of B to
rank 0 and all of C to rank 1, or some funny business like that.  If you
don't want to think about this ordering, extract the subvector the way I
described in the previous email.
-------------- next part --------------
A non-text attachment was scrubbed...
Name: signature.asc
Type: application/pgp-signature
Size: 818 bytes
Desc: not available
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20160407/80c6f2c9/attachment.pgp>