Actual source code: ex8.c
1:
2: static char help[] = "Demonstrates generating a slice from a DA Vector.\n\n";
4: #include petscda.h
5: #include petscao.h
9: /*
10: Given a DA generates a VecScatter context that will deliver a slice
11: of the global vector to each processor. In this example, each processor
12: receives the values i=*, j=*, k=rank, i.e. one z plane.
14: Note: This code is written assuming only one degree of freedom per node.
15: For multiple degrees of freedom per node use ISCreateBlock()
16: instead of ISCreateGeneral().
17: */
18: PetscErrorCode GenerateSliceScatter(DA da,VecScatter *scatter,Vec *vslice)
19: {
20: AO ao;
21: PetscInt M,N,P,nslice,*sliceindices,count,i,j;
22: PetscMPIInt rank;
24: MPI_Comm comm;
25: Vec vglobal;
26: IS isfrom,isto;
28: PetscObjectGetComm((PetscObject)da,&comm);
29: MPI_Comm_rank(comm,&rank);
31: DAGetAO(da,&ao);
32: DAGetInfo(da,0,&M,&N,&P,0,0,0,0,0,0,0);
34: /*
35: nslice is number of degrees of freedom in this processors slice
36: if there are more processors then z plans the extra processors get 0
37: elements in their slice.
38: */
39: if (rank < P) {nslice = M*N;} else nslice = 0;
41: /*
42: Generate the local vector to hold this processors slice
43: */
44: VecCreateSeq(PETSC_COMM_SELF,nslice,vslice);
45: DACreateGlobalVector(da,&vglobal);
47: /*
48: Generate the indices for the slice in the "natural" global ordering
49: Note: this is just an example, one could select any subset of nodes
50: on each processor. Just list them in the global natural ordering.
52: */
53: PetscMalloc((nslice+1)*sizeof(PetscInt),&sliceindices);
54: count = 0;
55: if (rank < P) {
56: for (j=0; j<N; j++) {
57: for (i=0; i<M; i++) {
58: sliceindices[count++] = rank*M*N + j*M + i;
59: }
60: }
61: }
62: /*
63: Convert the indices to the "PETSc" global ordering
64: */
65: AOApplicationToPetsc(ao,nslice,sliceindices);
66:
67: /* Create the "from" and "to" index set */
68: /* This is to scatter from the global vector */
69: ISCreateGeneral(PETSC_COMM_SELF,nslice,sliceindices,&isfrom);
70: /* This is to gather into the local vector */
71: ISCreateStride(PETSC_COMM_SELF,nslice,0,1,&isto);
72: VecScatterCreate(vglobal,isfrom,*vslice,isto,scatter);
73: ISDestroy(isfrom);
74: ISDestroy(isto);
75: PetscFree(sliceindices);
76: return 0;
77: }
82: int main(int argc,char **argv)
83: {
84: PetscMPIInt rank;
85: PetscInt m = PETSC_DECIDE,n = PETSC_DECIDE,p = PETSC_DECIDE,M = 3,N = 5,P=3,s=1;
86: PetscInt *lx = PETSC_NULL,*ly = PETSC_NULL,*lz = PETSC_NULL;
88: PetscTruth flg = PETSC_FALSE;
89: DA da;
90: Vec local,global,vslice;
91: PetscScalar value;
92: DAPeriodicType wrap = DA_XYPERIODIC;
93: DAStencilType stencil_type = DA_STENCIL_BOX;
94: VecScatter scatter;
96: PetscInitialize(&argc,&argv,(char*)0,help);
97: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
99: /* Read options */
100: PetscOptionsGetInt(PETSC_NULL,"-M",&M,PETSC_NULL);
101: PetscOptionsGetInt(PETSC_NULL,"-N",&N,PETSC_NULL);
102: PetscOptionsGetInt(PETSC_NULL,"-P",&P,PETSC_NULL);
103: PetscOptionsGetInt(PETSC_NULL,"-m",&m,PETSC_NULL);
104: PetscOptionsGetInt(PETSC_NULL,"-n",&n,PETSC_NULL);
105: PetscOptionsGetInt(PETSC_NULL,"-p",&p,PETSC_NULL);
106: PetscOptionsGetInt(PETSC_NULL,"-s",&s,PETSC_NULL);
107: PetscOptionsGetTruth(PETSC_NULL,"-star",&flg,PETSC_NULL);
108: if (flg) stencil_type = DA_STENCIL_STAR;
110: /* Create distributed array and get vectors */
111: DACreate3d(PETSC_COMM_WORLD,wrap,stencil_type,M,N,P,m,n,p,1,s,
112: lx,ly,lz,&da);
113: DAView(da,PETSC_VIEWER_DRAW_WORLD);
114: DACreateGlobalVector(da,&global);
115: DACreateLocalVector(da,&local);
117: GenerateSliceScatter(da,&scatter,&vslice);
119: /* Put the value rank+1 into all locations of vslice and transfer back to global vector */
120: value = 1.0 + rank;
121: VecSet(vslice,value);
122: VecScatterBegin(scatter,vslice,global,INSERT_VALUES,SCATTER_REVERSE);
123: VecScatterEnd(scatter,vslice,global,INSERT_VALUES,SCATTER_REVERSE);
125: VecView(global,PETSC_VIEWER_DRAW_WORLD);
127: VecDestroy(local);
128: VecDestroy(global);
129: DADestroy(da);
130: PetscFinalize();
131: return 0;
132: }