Actual source code: ex8.c
1: /*$Id: ex8.c,v 1.28 2001/08/07 03:04:42 balay Exp $*/
2:
3: static char help[] = "Demonstrates generating a slice from a DA Vector.\n\n";
5: #include petscda.h
6: #include petscsys.h
7: #include petscao.h
11: /*
12: Given a DA generates a VecScatter context that will deliver a slice
13: of the global vector to each processor. In this example, each processor
14: receives the values i=*, j=*, k=rank, i.e. one z plane.
16: Note: This code is written assuming only one degree of freedom per node.
17: For multiple degrees of freedom per node use ISCreateBlock()
18: instead of ISCreateGeneral().
19: */
20: int GenerateSliceScatter(DA da,VecScatter *scatter,Vec *vslice)
21: {
22: AO ao;
23: int M,N,P,nslice,rank,*sliceindices,count,ierr,i,j;
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(int),&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);
73: VecScatterCreate(vglobal,isfrom,*vslice,isto,scatter);
75: ISDestroy(isfrom);
76: ISDestroy(isto);
78: PetscFree(sliceindices);
79: return 0;
80: }
85: int main(int argc,char **argv)
86: {
87: int rank,M = 3,N = 5,P=3,s=1;
88: int m = PETSC_DECIDE,n = PETSC_DECIDE,p = PETSC_DECIDE,ierr;
89: int *lx = PETSC_NULL,*ly = PETSC_NULL,*lz = PETSC_NULL;
90: PetscTruth flg;
91: DA da;
92: Vec local,global,vslice;
93: PetscScalar value;
94: DAPeriodicType wrap = DA_XYPERIODIC;
95: DAStencilType stencil_type = DA_STENCIL_BOX;
96: VecScatter scatter;
98: PetscInitialize(&argc,&argv,(char*)0,help);
99: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
101: /* Read options */
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,"-m",&m,PETSC_NULL);
106: PetscOptionsGetInt(PETSC_NULL,"-n",&n,PETSC_NULL);
107: PetscOptionsGetInt(PETSC_NULL,"-p",&p,PETSC_NULL);
108: PetscOptionsGetInt(PETSC_NULL,"-s",&s,PETSC_NULL);
109: PetscOptionsHasName(PETSC_NULL,"-star",&flg);
110: if (flg) stencil_type = DA_STENCIL_STAR;
112: /* Create distributed array and get vectors */
113: DACreate3d(PETSC_COMM_WORLD,wrap,stencil_type,M,N,P,m,n,p,1,s,
114: lx,ly,lz,&da);
115: DAView(da,PETSC_VIEWER_DRAW_WORLD);
116: DACreateGlobalVector(da,&global);
117: DACreateLocalVector(da,&local);
119: GenerateSliceScatter(da,&scatter,&vslice);
121: /* Put the value rank+1 into all locations of vslice and transfer back to global vector */
122: value = 1.0 + rank;
123: VecSet(&value,vslice);
124: VecScatterBegin(vslice,global,INSERT_VALUES,SCATTER_REVERSE,scatter);
125: VecScatterEnd(vslice,global,INSERT_VALUES,SCATTER_REVERSE,scatter);
127: VecView(global,PETSC_VIEWER_DRAW_WORLD);
129: VecDestroy(local);
130: VecDestroy(global);
131: DADestroy(da);
132: PetscFinalize();
133: return 0;
134: }