[petsc-users] Ordering with MatCreateNest and MPI

Mon Apr 11 08:16:18 CDT 2016

Dear Jed,

thanks for the reply. Despite my efforts, I'm still unable to use the 
information carried by the index set. I'm able to extract the IS through 
MatNestGetISs but I don't really know how to interprete it and to exploit it 
to pursue my goal.

I don't seem to find a tutorial explaining how to do the following:

"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."

Can you please give me some references for the above operation?

Thanks,
Alessandro

On Fri, 08 Apr 2016 09:59:28 +0100
  Jed Brown <jed at jedbrown.org> wrote:
> 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.

___________________________________
Matteo Alessandro Francavilla, PhD
Antenna and EMC Lab (LACE)
Istituto Superiore Mario Boella (ISMB)
Torino, Italy
phone +39 011 2276 711