[SOLVED] Petsc - combine distributed vector into local vector

Petsc - combine distributed vector into local vector

I am using Petsc, and I would like to combine a distributed Vec so that every process has a complete copy of the Vec. I have a minimal example which starts with an array of data, constructs an MPI Vec from it, and then tries to use a VecScatter to combine the vectors from multiple processes. When I do this, the local vector only receives the values that were stored on the 0th process, it does not receive the information from the other process. How do I combine the distributed vector to produce a complete, local vector?

#include <petscvec.h>

double primes[] = {2,3,5,7,11,13,17};
int nprimes = 7;

int main(int argc,char **argv)
{
    PetscInitialize(&argc,&argv, NULL,NULL);

    MPI_Comm       comm=MPI_COMM_WORLD;
    Vec xpar,xseq;
    PetscInt low,high;
    IS index_set_global, index_set_local;
    const PetscInt *indices;
    VecScatter vc;
    PetscErrorCode ierr;

    //Set up parallel vector
    ierr = VecCreateMPI(comm, PETSC_DETERMINE, nprimes, &xpar); CHKERRQ(ierr);
    ierr = VecGetOwnershipRange(xpar, &low, &high); CHKERRQ(ierr);
    ierr = ISCreateStride(comm, high - low, low, 1, &index_set_global); CHKERRQ(ierr);
    ierr = ISGetIndices(index_set_global, &indices); CHKERRQ(ierr);
    ierr = ISView(index_set_global, PETSC_VIEWER_STDOUT_WORLD); CHKERRQ(ierr);
    ierr = VecSetValues(xpar, high - low, indices, primes + low, INSERT_VALUES);CHKERRQ(ierr);
    ierr = VecAssemblyBegin(xpar); CHKERRQ(ierr);
    ierr = VecAssemblyEnd(xpar); CHKERRQ(ierr);
    ierr = VecView(xpar, PETSC_VIEWER_STDOUT_WORLD); CHKERRQ(ierr);

    //Scatter parallel vector so all processes have full vector
    ierr = VecCreateSeq(PETSC_COMM_SELF, nprimes, &xseq); CHKERRQ(ierr);
    //ierr = VecCreateMPI(comm, high - low, nprimes, &xseq); CHKERRQ(ierr);
    ierr = ISCreateStride(comm, high - low, 0, 1, &index_set_local); CHKERRQ(ierr);
    ierr = VecScatterCreate(xpar, index_set_local, xseq, index_set_global, &vc); CHKERRQ(ierr);
    ierr = VecScatterBegin(vc, xpar, xseq, ADD_VALUES, SCATTER_FORWARD); CHKERRQ(ierr);
    ierr = VecScatterEnd(vc, xpar, xseq, ADD_VALUES, SCATTER_FORWARD); CHKERRQ(ierr);
    ierr = PetscPrintf(PETSC_COMM_SELF, "\nPrinting out scattered vector\n"); CHKERRQ(ierr);
    ierr = VecView(xseq, PETSC_VIEWER_STDOUT_WORLD); CHKERRQ(ierr);

    PetscFinalize();
}

OUTPUT:

mpiexec -n 2 ./test

IS Object: 2 MPI processes
type: stride
[0] Index set is permutation
[0] Number of indices in (stride) set 4
[0] 0 0
[0] 1 1
[0] 2 2
[0] 3 3
[1] Number of indices in (stride) set 3
[1] 0 4
[1] 1 5
[1] 2 6
Vec Object: 2 MPI processes
type: mpi
Process [0]
2.
3.
5.
7.
Process [1]
11.
13.
17.

Printing out scattered vector

Printing out scattered vector
Vec Object: 1 MPI processes
type: seq
2.
3.
5.
7.
0.
0.
0.

Solution

VecScatterCreateToAll() is exactly what you need:

Creates a vector and a scatter context that copies all vector values to each processor

It is used in ksp/.../ex49.c. Lasty, it is implemented in vecmpitoseq.c.

The naming convention is likely inspired by MPI functions such as MPI_Allgather(), which distribute the gathered data to all processes while MPI_Gather() only gathers the data on the specified root process.