[petsc-users] Ordering with MatCreateNest and MPI

[Alessandro Francavilla]

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

Thanks,
Alessandro


-----Original Message-----
From: Jed Brown [mailto:jed at jedbrown.org] 
Sent: giovedì 7 aprile 2016 15:26
To: FRANCAVILLA MATTEO ALESSANDRO; petsc-users at mcs.anl.gov
Subject: Re: [petsc-users] Ordering with MatCreateNest and MPI

FRANCAVILLA MATTEO ALESSANDRO <d019598 at polito.it> writes:

> Hi,
>
> I'm trying to use a matrix composed of different sub-matrix blocks; I 
> create the matrix with MatCreateNest(). I'm finding it hard to 
> understand the ordering (in MPI applications) of the result of MatMult
operations.
> Apparently the result of the multiplication is correct, except that 
> both input and output vectors go through some sort of permutation. How 
> should I use the MatMult routine, i.e. how to pass the input vector 
> with the correct ordering and how to get the correct ordering of the
result?
>
> Below a simple Fortran example to demonstrate what I'm trying to do 
> (running the example with 1 process works fine, with  2 or more 
> processes I don't understand how to get the correct ordering). When 
> using standard matrices (i.e., not MatNest) I've never had to care about
orderings.

MatNest is semantically identical to other matrix formats.  In fact, we
strongly recommend that you write your code so that your assembly and all
subsequent operations do not depend in any way on the matrix being MATNEST
of some other format.  See src/snes/examples/tutorials/ex28.c
for an example of this.

But anyway, you can use VecGetSubVector to separate out the components
(MatNestGetISs to get the index sets if you don't already have them).