#include <petscvec.h>

int main(int argc,char **argv)
{
  PetscErrorCode    ierr;
  Vec               x,xseq,perm_vecseq,perm_vecmpi;
  const PetscReal   *xdata,*perm_vecmpi_data;
  PetscInt          i,n,N,*iperm;
  PetscReal         *rperm;
  PetscRandom       rnd;
  VecScatter        vscat;

  ierr = PetscInitialize(&argc,&argv,NULL,NULL);if (ierr) return ierr;

  /* Build a random vector of length 10 */
  ierr = VecCreate(PETSC_COMM_WORLD,&x);CHKERRQ(ierr);
  ierr = VecSetSizes(x,PETSC_DECIDE,10);CHKERRQ(ierr);
  ierr = VecSetFromOptions(x);CHKERRQ(ierr);

  ierr = PetscRandomCreate(PETSC_COMM_WORLD,&rnd);CHKERRQ(ierr);
  ierr = VecSetRandom(x,rnd);CHKERRQ(ierr);
  ierr = PetscRandomDestroy(&rnd);CHKERRQ(ierr);

  /* perm_vecmpi has the same layout as x. It will store permutation of x (in real) */
  ierr = VecDuplicate(x,&perm_vecmpi);CHKERRQ(ierr);

  /* Gather x to xseq on rank 0. xseq on other ranks are actually empty */
  ierr = VecScatterCreateToZero(x,&vscat,&xseq);CHKERRQ(ierr);
  ierr = VecScatterBegin(vscat,x,xseq,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
  ierr = VecScatterEnd(vscat,x,xseq,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);

  /* Do PetscSortRealWithPermutation on xseq */

  /* If want descending order, you need ierr = VecScale(xseq,-1.0);CHKERRQ(ierr); */

  ierr = VecGetLocalSize(xseq,&N);CHKERRQ(ierr); /* N=0 on ranks other than rank 0 */
  ierr = PetscMalloc2(N,&iperm,N,&rperm);CHKERRQ(ierr);
  for (i=0; i<N; i++) iperm[i] = i;
  ierr = VecGetArrayRead(xseq,&xdata);CHKERRQ(ierr);
  ierr = PetscSortRealWithPermutation(N,xdata,iperm);CHKERRQ(ierr);
  ierr = VecRestoreArrayRead(xseq,&xdata);CHKERRQ(ierr);

  /* Convert iperm[] to rperm[], i.e., 1 to 1.0, 2 to 2.0, etc */
  for (i=0; i<N; i++) rperm[i] = (PetscScalar)iperm[i];

  /* Wrap rperm[] in a sequential vector perm_vecseq */
  ierr = VecCreateSeqWithArray(PETSC_COMM_SELF,1,N,rperm,&perm_vecseq);CHKERRQ(ierr);

  /* Do reverse vecscatter to scatter values in perm_vecseq to perm_vecmpi */
  ierr = VecScatterBegin(vscat,perm_vecseq,perm_vecmpi,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
  ierr = VecScatterEnd(vscat,perm_vecseq,perm_vecmpi,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);

  /* Free perm_vecseq, iperm, rperm */
  ierr = VecDestroy(&perm_vecseq);CHKERRQ(ierr);
  ierr = PetscFree2(iperm,rperm);CHKERRQ(ierr);

  /* Read values from perm_vecmpi can convert 1.0 to 1, 2.0 to 2, etc */
  ierr = VecGetArrayRead(perm_vecmpi,&perm_vecmpi_data);CHKERRQ(ierr);
  ierr = VecGetLocalSize(x,&n);CHKERRQ(ierr);
  ierr = PetscMalloc1(n,&iperm);CHKERRQ(ierr); /* Allocate iperm again, in different size */
  for (i=0; i<n; i++) iperm[i] = (PetscInt)perm_vecmpi_data[i];
  ierr = VecRestoreArrayRead(perm_vecmpi,&perm_vecmpi_data);CHKERRQ(ierr);

  /* iperm[] now has the permutation in PetscInt. In a global view, if j = iperm[i], then i-th entry of a sorted x[] would be x[j] */

  /* Print and verify */
  ierr = VecView(x,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
  ierr = VecView(perm_vecmpi,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);

  ierr = VecScatterDestroy(&vscat);CHKERRQ(ierr);
  ierr = VecDestroy(&xseq);CHKERRQ(ierr);
  ierr = VecDestroy(&perm_vecmpi);CHKERRQ(ierr);
  ierr = VecDestroy(&x);CHKERRQ(ierr);
  ierr = PetscFree(iperm);CHKERRQ(ierr);
  ierr = PetscFinalize();
  return ierr;
}

/*TEST

   test:
      nsize: 3

TEST*/