[petsc-users] Local and global size of IS

[Smith, Barry F.]

> On Mar 27, 2019, at 7:33 AM, Eda Oktay via petsc-users <petsc-users at mcs.anl.gov> wrote:
> 
> Hello,
> 
> I am trying to permute a matrix A(of size 2n*2n) according to a sorted eigenvector vr (of size 2n) in parallel using 2 processors (processor number can change).

    I don't understand your rational for wanting to permute the matrix rows/columns based on a sorted eigenvector? 

> However, I get an error in MatPermute line stating that arguments are out of range and a new nonzero caused a malloc even if I used MatSetOption. 
> 
> I discovered that this may be because of the unequality of local sizes of is (and newIS) and local size of A. 
> 
> Since I allocate index set idx according to size of the vector vr, global size of is becomes 2n and the local size is also 2n (I think it should be n since both A and vr has local sizes n because of the number of processors). If I change the size of idx, then because of VecGetArrayRead, I think is is created wrongly.
> 
> So, how can I make both global and local sizes of is,newIS and A?
> 
> Below, you can see the part of my program.
> 
> Thanks,
> 
> Eda
>  
>  ierr = VecGetSize(vr,&siz);CHKERRQ(ierr);                            
>   ierr = PetscMalloc1(siz,&idx);CHKERRQ(ierr);
>   for (i=0; i<siz;i++) idx[i] = i;

   VecGetArray only gives you access to the local portion of the vector, not the entire vector so you cannot do the next two lines of code
>   ierr = VecGetArrayRead(vr,&avr);CHKERRQ(ierr);

   The sort routine is sequential; it knows nothing about parallelism so each process is just sorting its part. Hence this code won't work as expected

>   ierr = PetscSortRealWithPermutation(siz,avr,idx);CHKERRQ(ierr); 

   If you need to sort the parallel vector and get the entire permutation then you need to do the following 

1)   communicate the entire vector to each process, you can use VecScatterCreateToAll() for this
2)   call VecGetArray on the new vector (that contains all entries on each process)
3)   Call PetscSortRealWithPermutation() on the entire vector on each process
4)   select out a piece of the resulting indices idx on each process; for example with two processes I think rank = 0 would get the first half of the idx and rank = 1 would get the second half.
>                  
>   ierr = ISCreateGeneral(PETSC_COMM_SELF,siz,idx,PETSC_COPY_VALUES,&is);CHKERRQ(ierr);  
>   ierr = ISSetPermutation(is);CHKERRQ(ierr);

  You don't need to duplicate the is, just pass it in twice.

>   ierr = ISDuplicate(is,&newIS);CHKERRQ(ierr);

   You should definitely not need the next line. The A matrix is untouched by the subroutine call so you don't need to change its allocation options

> MatSetOption(A,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);CHKERRQ(ierr);  
>   ierr = MatPermute(A,is,newIS,&PL);CHKERRQ(ierr);