Hi Barry, I did what you said but now I get the error in ISSetPermutation that index set is not a permutation. ı think it is because since I selected indices on each process as you said. I did the following, did I understand you wrong:
PetscMPIInt rank,size; MPI_Comm_rank(PETSC_COMM_WORLD, &rank); MPI_Comm_size(PETSC_COMM_WORLD, &size); PetscInt *idxx,ss; ss = siz/size; PetscMalloc1(ss,&idxx); if (rank != size-1) { j =0; for (i=rank*ss; i<(rank+1)*ss; i++) { idxx[j] = idx[i]; j++; } } else { j =0; for (i=rank*ss; i<siz; i++) { idxx[j] = idx[i]; j++; } } Smith, Barry F. <bsm...@mcs.anl.gov>, 31 Mar 2019 Paz, 01:47 tarihinde şunu yazdı: > > > > On Mar 27, 2019, at 7:33 AM, Eda Oktay via petsc-users < > petsc-users@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); > >