[petsc-users] Ordering with MatCreateNest and MPI

[Jed Brown]

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

> In fact I don't expect anything special about indices, I just can't
> understand how to use them.
>
> If I set a matrix [A B; C D], and each of the blocks is distributed over
> some arbitrary number of processes, I want to know how to compute the
> product
> y = [A B; C D]*x
> and be able to see the result. I'm not concerned about how the elements are
> distributed, I just want to be able to execute the matmult operation and
> view the result.

MatNestGetISs give you the ISs to extract or set the portions of x and y
corresponding to the first and second block.  The index sets are not
globally contiguous and can't be if you want sane parallelism.

>> It is returning the ith column of A.  A is not just matArray(1).
>
> No, it is not. I know very well that A is not matArray(1), and this
> operation is returning one column of A, but definitely not i-th
> column. 

Evidently the ith column is not defined the same way as your internal
mental model.  Use the index sets.

> I'm definitely sure it is my mistake: 1) building the vector x the way
> I am setting it does not actually do what I'm expecting it to do; 2)
> even if I fix (1) I still need to permute the result in order to
> obtain y as defined by y=[A B;C D]*x .
>
> This is what's going on, which summarizes my 2 mistakes: (1) the result
> needs be permuted, and (2) the input vector x is not correct
>
> This is what the matrix looks like:
>
> ! A =            [   1   2 | 101 102 ;
> !          	          3   4 | 102 103 ;
> !	___________________
> !                 201 202 | 401 402 ;
> !	   203 204 | 403 404 ]
>
>
> If I set the x vector with CALL
> VecSetValues(x,1,0,(1.+0.*PETSC_i),INSERT_VALUES,IERR)
>
> The result with 1 MPI process is what I expect (the first column of A):
>
> Vec Object: 1 MPI processes
>   type: seq
> 1
> 3
> 201
> 203
>
> With 2 MPI processes the results are distributed: the result is correct but
> I need to reorder it to get [1 3 201 203]'  but don't know how:
>
> Vec Object: 2 MPI processes
>   type: mpi
> Process [0]
> 1
> 201
> Process [1]
> 3
> 203
>
>
> If I set the x vector with CALL
> VecSetValues(x,1,1,(1.+0.*PETSC_i),INSERT_VALUES,IERR)
>
> 1 MPI process correctly yields the 2nd column of A:
>
> Vec Object: 1 MPI processes
>   type: seq
> 2
> 4
> 202
> 204
>
> 2 MPI processes DO NOT YIELD the 2nd column, but the 3rd column (apart from
> a permutation similarly to the previous case):
>
> Vec Object: 2 MPI processes
>   type: mpi
> Process [0]
> 101
> 301
> Process [1]
> 103
> 303
>
>
> -----Original Message-----
> From: Jed Brown [mailto:jed at jedbrown.org] 
> Sent: giovedì 7 aprile 2016 20:04
> To: Alessandro Francavilla
> Cc: petsc-users at mcs.anl.gov
> Subject: RE: [petsc-users] Ordering with MatCreateNest and MPI
>
> 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/20160408/c2bbec3d/attachment.pgp>