Actual source code: ex8.c

petsc-3.3-p7 2013-05-11
  2: static char help[] = "Demonstrates using a local ordering to set values into a parallel vector.\n\n";

  4: /*T
  5:    Concepts: vectors^assembling vectors with local ordering;
  6:    Processors: n
  7: T*/

  9: /* 
 10:   Include "petscvec.h" so that we can use vectors.  Note that this file
 11:   automatically includes:
 12:      petscsys.h       - base PETSc routines   petscis.h     - index sets
 13:      petscviewer.h - viewers
 14: */
 15: #include <petscvec.h>

 19: int main(int argc,char **argv)
 20: {
 22:   PetscMPIInt    rank;
 23:   PetscInt       i,N,ng,*gindices,rstart,rend,M;
 24:   PetscScalar    one = 1.0;
 25:   Vec            x;

 27:   PetscInitialize(&argc,&argv,(char *)0,help);
 28:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);

 30:   /*
 31:      Create a parallel vector.
 32:       - In this case, we specify the size of each processor's local
 33:         portion, and PETSc computes the global size.  Alternatively,
 34:         PETSc could determine the vector's distribution if we specify
 35:         just the global size.
 36:   */
 37:   VecCreate(PETSC_COMM_WORLD,&x);
 38:   VecSetSizes(x,rank+1,PETSC_DECIDE);
 39:   VecSetFromOptions(x);
 40:   VecGetSize(x,&N);
 41:   VecSet(x,one);

 43:   /*
 44:      Set the local to global ordering for the vector. Each processor 
 45:      generates a list of the global indices for each local index. Note that
 46:      the local indices are just whatever is convenient for a particular application.
 47:      In this case we treat the vector as lying on a one dimensional grid and 
 48:      have one ghost point on each end of the blocks owned by each processor. 
 49:   */

 51:   VecGetSize(x,&M);
 52:   VecGetOwnershipRange(x,&rstart,&rend);
 53:   ng   = rend - rstart + 2;
 54:   PetscMalloc(ng*sizeof(PetscInt),&gindices);
 55:   gindices[0] = rstart - 1;
 56:   for (i=0; i<ng-1; i++) {
 57:     gindices[i+1] = gindices[i] + 1;
 58:   }
 59:   /* map the first and last point as periodic */
 60:   if (gindices[0]    == -1) gindices[0]    = M - 1;
 61:   if (gindices[ng-1] == M)  gindices[ng-1] = 0;
 62:   {
 63:     ISLocalToGlobalMapping ltog;
 64:     ISLocalToGlobalMappingCreate(PETSC_COMM_SELF,ng,gindices,PETSC_COPY_VALUES,&ltog);
 65:     VecSetLocalToGlobalMapping(x,ltog);
 66:     ISLocalToGlobalMappingDestroy(&ltog);
 67:   }
 68:   PetscFree(gindices);

 70:   /*
 71:      Set the vector elements.
 72:       - In this case set the values using the local ordering
 73:       - Each processor can contribute any vector entries,
 74:         regardless of which processor "owns" them; any nonlocal
 75:         contributions will be transferred to the appropriate processor
 76:         during the assembly process.
 77:       - In this example, the flag ADD_VALUES indicates that all
 78:         contributions will be added together.
 79:   */
 80:   for (i=0; i<ng; i++) {
 81:     VecSetValuesLocal(x,1,&i,&one,ADD_VALUES);
 82:   }

 84:   /* 
 85:      Assemble vector, using the 2-step process:
 86:        VecAssemblyBegin(), VecAssemblyEnd()
 87:      Computations can be done while messages are in transition
 88:      by placing code between these two statements.
 89:   */
 90:   VecAssemblyBegin(x);
 91:   VecAssemblyEnd(x);

 93:   /*
 94:       View the vector; then destroy it.
 95:   */
 96:   VecView(x,PETSC_VIEWER_STDOUT_WORLD);
 97:   VecDestroy(&x);

 99:   PetscFinalize();
100:   return 0;
101: }
102: