Actual source code: ex34.c

petsc-3.4.2 2013-07-02

  3: static char help[] = "Demonstrates PetscHMPIMerge() usage\n\n";

  5: #include <petscmat.h>
  6: #include <petscksp.h>

  8: typedef struct {
  9:   MPI_Comm   comm;
 10:   Mat        A;
 11:   Vec        x,y;      /* contains the vector values spread across all the processes */
 12:   Vec        xr,yr;    /* contains the vector values on the master processes, all other processes have zero length */
 13:   VecScatter sct;
 14: } MyMultCtx;

 18: /*
 19:     This is called on ALL processess, master and worker
 20: */
 21: PetscErrorCode MyMult(MPI_Comm comm,MyMultCtx *ctx,void *dummy)
 22: {

 26:   PetscSynchronizedPrintf(ctx->comm,"Doing multiply\n");
 27:   PetscSynchronizedFlush(ctx->comm);
 28:   /* moves data that lives only on master processes to all processes */
 29:   VecScatterBegin(ctx->sct,ctx->xr,ctx->x,INSERT_VALUES,SCATTER_FORWARD);
 30:   VecScatterEnd(ctx->sct,ctx->xr,ctx->x,INSERT_VALUES,SCATTER_FORWARD);
 31:   MatMult(ctx->A,ctx->x,ctx->y);
 32:   /* moves data that lives on all processes to master processes */
 33:   VecScatterBegin(ctx->sct,ctx->y,ctx->yr,INSERT_VALUES,SCATTER_REVERSE);
 34:   VecScatterEnd(ctx->sct,ctx->y,ctx->yr,INSERT_VALUES,SCATTER_REVERSE);
 35:   return(0);
 36: }

 40: /*
 41:     This is called only on the master processes
 42: */
 43: PetscErrorCode MySubsolver(MyMultCtx *ctx)
 44: {
 46:   void           *subctx;

 49:   PetscSynchronizedPrintf(PETSC_COMM_WORLD,"MySubsolver\n");
 50:   PetscSynchronizedFlush(PETSC_COMM_WORLD);
 51:   /* allocates memory on each process, both masters and workers */
 52:   PetscHMPIMalloc(PETSC_COMM_LOCAL_WORLD,sizeof(int),&subctx);
 53:   /* runs MyMult() function on each process, both masters and workers */
 54:   PetscHMPIRunCtx(PETSC_COMM_LOCAL_WORLD,(PetscErrorCode (*)(MPI_Comm,void*,void*))MyMult,subctx);
 55:   PetscHMPIFree(PETSC_COMM_LOCAL_WORLD,subctx);
 56:   return(0);
 57: }

 61: int main(int argc,char **args)
 62: {
 64:   char           file[PETSC_MAX_PATH_LEN];
 65:   PetscViewer    fd;
 66:   PetscMPIInt    rank,size,nodesize = 1;
 67:   MyMultCtx      ctx;
 68:   const PetscInt *ns; /* length of vector ctx.x on all process */
 69:   PetscInt       i,rstart,n = 0;   /* length of vector ctx.xr on this process */
 70:   IS             is;

 72:   PetscInitialize(&argc,&args,(char*)0,help);
 73:   ctx.comm = PETSC_COMM_WORLD;
 74:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
 75:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
 76:   PetscOptionsGetInt(NULL,"-nodesize",&nodesize,NULL);
 77:   if (size % nodesize) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"MPI_COMM_WORLD size must be divisible by nodesize");

 79:   /* Read matrix */
 80:   PetscOptionsGetString(NULL,"-f",file,PETSC_MAX_PATH_LEN,NULL);
 81:   PetscViewerBinaryOpen(PETSC_COMM_WORLD,file,FILE_MODE_READ,&fd);
 82:   MatCreate(PETSC_COMM_WORLD,&ctx.A);
 83:   MatSetType(ctx.A,MATMPIAIJ);
 84:   MatLoad(ctx.A,fd);
 85:   PetscViewerDestroy(&fd);

 87:   /* Create work vectors for matrix-vector product */
 88:   MatGetVecs(ctx.A,&ctx.x,&ctx.y);
 89:   VecGetOwnershipRanges(ctx.x,&ns);
 90:   if (!(rank % nodesize)) { /* I am master process; I will own all vector elements on all my worker processes*/
 91:     for (i=0; i<nodesize; i++) n += ns[rank+i+1] - ns[rank+i];
 92:   }
 93:   VecCreateMPI(MPI_COMM_WORLD,n,PETSC_DETERMINE,&ctx.xr);
 94:   VecDuplicate(ctx.xr,&ctx.yr);
 95:   /* create scatter from ctx.xr to ctx.x vector */
 96:   VecGetOwnershipRange(ctx.x,&rstart,NULL);
 97:   ISCreateStride(PETSC_COMM_WORLD,ns[rank],rstart,1,&is);
 98:   VecScatterCreate(ctx.xr,is,ctx.x,is,&ctx.sct);
 99:   ISDestroy(&is);

101:   /*
102:      The master nodes call the function MySubsolver() while the worker nodes wait for requests to call functions
103:      These requests are triggered by the calls from the masters on PetscHMPIRunCtx()
104:   */
105:   PetscHMPIMerge(nodesize,(PetscErrorCode (*)(void*))MySubsolver,&ctx);

107:   PetscHMPIFinalize();
108:   MatDestroy(&ctx.A);
109:   VecDestroy(&ctx.x);
110:   VecDestroy(&ctx.y);
111:   VecDestroy(&ctx.xr);
112:   VecDestroy(&ctx.yr);
113:   VecScatterDestroy(&ctx.sct);

115:   PetscFinalize();
116:   return 0;
117: }