Actual source code: ex26.c
1: /*$Id: ex26.c,v 1.13 2001/09/11 16:32:10 bsmith Exp $*/
2: /*
4: Test program follows. Writing it I realised that
5: 1/ I built the pipeline object around an MPI-to-MPI vector scatter.
6: 2/ That necessitates the 'little trick' below.
7: 3/ This trick will always be necessary, since the code inside the
8: pipe is a critical section.
9: 4/ Hence I really should have used an MPI-to-Seq scatter.
10: 5/ This shouldn't be too hard to fix in the implementation you
11: guys are making,right? :-) <-- smiley just in case.
13: If this is not clear, I'll try to elaborate.
15: */
16: /* Example of pipeline code:
17: accumulation of the sum $s_p=\sum_{q\leq p} (q+1)^2$.
18: E.g., processor 3 computes 1^2+2^2+3^+4^2 = 30.
19: Every processor computes its term, then passes it on to the next.
20: */
21: #include petscvec.h
25: int main(int Argc,char **Args)
26: {
27: Vec src_v,tar_v,loc_v;
28: IS src_idx,tar_idx;
29: VecPipeline pipe;
30: MPI_Comm comm;
31: int size,rank,src_loc,tar_loc,ierr,zero_loc=0;
32: PetscScalar zero=0,my_value,*vec_values,*loc_ar;
34: PetscInitialize(&Argc,&Args,PETSC_NULL,PETSC_NULL);
36: comm = MPI_COMM_WORLD;
37: MPI_Comm_size(comm,&size);
38: MPI_Comm_rank(comm,&rank);
39:
40: /* Create the necessary vectors; one element per processor */
41: VecCreate(comm,&tar_v);
42: VecSetSizes(tar_v,1,PETSC_DECIDE);
43: VecSetFromOptions(tar_v);
44: VecSet(&zero,tar_v);
45: VecCreate(comm,&src_v);
46: VecSetSizes(src_v,1,PETSC_DECIDE);
47: VecSetFromOptions(src_v);
48: VecCreateSeq(MPI_COMM_SELF,1,&loc_v);
49: /* -- little trick: we need a distributed and a local vector
50: that share each other's data; see below for application */
51: VecGetArray(loc_v,&loc_ar);
52: VecPlaceArray(src_v,loc_ar);
54: /* Create the pipeline data: we write into our own location,
55: and read one location from the left */
56: tar_loc = rank;
57: if (tar_loc>0) src_loc = tar_loc-1; else src_loc = tar_loc;
58: ISCreateGeneral(MPI_COMM_SELF,1,&tar_loc,&tar_idx);
59: ISCreateGeneral(MPI_COMM_SELF,1,&src_loc,&src_idx);
60: VecPipelineCreate(comm,src_v,src_idx,tar_v,tar_idx,&pipe);
61: VecPipelineSetType(pipe,PIPELINE_SEQUENTIAL,PETSC_NULL);
62: VecPipelineSetup(pipe);
64: /* The actual pipe:
65: receive accumulated value from previous processor,
66: add the square of your own value, and send on. */
67: VecPipelineBegin(src_v,tar_v,INSERT_VALUES,SCATTER_FORWARD,PIPELINE_UP,pipe);
68: VecGetArray(tar_v,&vec_values);
69: my_value = vec_values[0] + (PetscReal)((rank+1)*(rank+1));
71: VecRestoreArray(tar_v,&vec_values);CHKERRQ(ierr)
72: /* -- little trick: we have to be able to call VecAssembly,
73: but since this code executed sequentially (critical section!),
74: we have a local vector with data aliased to the distributed one */
75: VecSetValues(loc_v,1,&zero_loc,&my_value,INSERT_VALUES);
76: VecAssemblyBegin(loc_v);
77: VecAssemblyEnd(loc_v);
78: VecPipelineEnd(src_v,tar_v,INSERT_VALUES,SCATTER_FORWARD,PIPELINE_UP,pipe);
80: PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d] value=%d\n",rank,(int)PetscRealPart(my_value));
81: PetscSynchronizedFlush(PETSC_COMM_WORLD);
83: /* Clean up */
84: VecPipelineDestroy(pipe);
85: VecDestroy(src_v);
86: VecDestroy(tar_v);
87: VecDestroy(loc_v);
88: ISDestroy(src_idx);
89: ISDestroy(tar_idx);
91: PetscFinalize();
93: return 0;
94: }