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: }