[petsc-users] Ordering with MatCreateNest and MPI

[Alessandro Francavilla]

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.

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