Actual source code: ex8.c
petsc-3.3-p6 2013-02-11
1:
2: static char help[] = "Demonstrates generating a slice from a DMDA Vector.\n\n";
4: #include <petscdmda.h>
5: #include <petscao.h>
9: /*
10: Given a DMDA 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(DM 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: DMDAGetAO(da,&ao);
32: DMDAGetInfo(da,0,&M,&N,&P,0,0,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: DMCreateGlobalVector(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,PETSC_OWN_POINTER,&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: return 0;
76: }
81: int main(int argc,char **argv)
82: {
83: PetscMPIInt rank;
84: PetscInt m = PETSC_DECIDE,n = PETSC_DECIDE,p = PETSC_DECIDE,M = 3,N = 5,P=3,s=1;
85: PetscInt *lx = PETSC_NULL,*ly = PETSC_NULL,*lz = PETSC_NULL;
87: PetscBool flg = PETSC_FALSE;
88: DM da;
89: Vec local,global,vslice;
90: PetscScalar value;
91: DMDABoundaryType wrap = DMDA_XYPERIODIC;
92: DMDAStencilType stencil_type = DMDA_STENCIL_BOX;
93: VecScatter scatter;
95: PetscInitialize(&argc,&argv,(char*)0,help);
96: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
98: /* Read options */
99: PetscOptionsGetInt(PETSC_NULL,"-M",&M,PETSC_NULL);
100: PetscOptionsGetInt(PETSC_NULL,"-N",&N,PETSC_NULL);
101: PetscOptionsGetInt(PETSC_NULL,"-P",&P,PETSC_NULL);
102: PetscOptionsGetInt(PETSC_NULL,"-m",&m,PETSC_NULL);
103: PetscOptionsGetInt(PETSC_NULL,"-n",&n,PETSC_NULL);
104: PetscOptionsGetInt(PETSC_NULL,"-p",&p,PETSC_NULL);
105: PetscOptionsGetInt(PETSC_NULL,"-s",&s,PETSC_NULL);
106: PetscOptionsGetBool(PETSC_NULL,"-star",&flg,PETSC_NULL);
107: if (flg) stencil_type = DMDA_STENCIL_STAR;
109: /* Create distributed array and get vectors */
110: DMDACreate3d(PETSC_COMM_WORLD,wrap,stencil_type,M,N,P,m,n,p,1,s,
111: lx,ly,lz,&da);
112: DMView(da,PETSC_VIEWER_DRAW_WORLD);
113: DMCreateGlobalVector(da,&global);
114: DMCreateLocalVector(da,&local);
116: GenerateSliceScatter(da,&scatter,&vslice);
118: /* Put the value rank+1 into all locations of vslice and transfer back to global vector */
119: value = 1.0 + rank;
120: VecSet(vslice,value);
121: VecScatterBegin(scatter,vslice,global,INSERT_VALUES,SCATTER_REVERSE);
122: VecScatterEnd(scatter,vslice,global,INSERT_VALUES,SCATTER_REVERSE);
124: VecView(global,PETSC_VIEWER_DRAW_WORLD);
126: VecDestroy(&local);
127: VecDestroy(&global);
128: DMDestroy(&da);
129: PetscFinalize();
130: return 0;
131: }