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