Actual source code: ex34.c
3: static char help[] = "Demonstrates PetscOpenMPMerge() 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: PetscOpenMPMalloc(PETSC_COMM_LOCAL_WORLD,sizeof(int),&subctx);
53: /* runs MyMult() function on each process, both masters and workers */
54: PetscOpenMPRunCtx(PETSC_COMM_LOCAL_WORLD,(PetscErrorCode (*)(MPI_Comm,void*,void *))MyMult,subctx);
55: PetscOpenMPFree(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(PETSC_NULL,"-nodesize",&nodesize,PETSC_NULL);
77: if (size % nodesize) SETERRQ(PETSC_ERR_SUP,"MPI_COMM_WORLD size must be divisible by nodesize");
79: /* Read matrix */
80: PetscOptionsGetString(PETSC_NULL,"-f",file,PETSC_MAX_PATH_LEN-1,PETSC_NULL);
81: PetscViewerBinaryOpen(PETSC_COMM_WORLD,file,FILE_MODE_READ,&fd);
82: MatLoad(fd,MATMPIAIJ,&ctx.A);
83: PetscViewerDestroy(fd);
85: /* Create work vectors for matrix-vector product */
86: MatGetVecs(ctx.A,&ctx.x,&ctx.y);
87: VecGetOwnershipRanges(ctx.x,&ns);
88: if (!(rank % nodesize)) { /* I am master process; I will own all vector elements on all my worker processes*/
89: for (i=0; i<nodesize; i++) n += ns[rank+i+1] - ns[rank+i];
90: }
91: VecCreateMPI(MPI_COMM_WORLD,n,PETSC_DETERMINE,&ctx.xr);
92: VecDuplicate(ctx.xr,&ctx.yr);
93: /* create scatter from ctx.xr to ctx.x vector */
94: VecGetOwnershipRange(ctx.x,&rstart,PETSC_NULL);
95: ISCreateStride(PETSC_COMM_WORLD,ns[rank],rstart,1,&is);
96: VecScatterCreate(ctx.xr,is,ctx.x,is,&ctx.sct);
97: ISDestroy(is);
99: /*
100: The master nodes call the function MySubsolver() while the worker nodes wait for requests to call functions
101: These requests are triggered by the calls from the masters on PetscOpenMPRunCtx()
102: */
103: PetscOpenMPMerge(nodesize,(PetscErrorCode (*)(void*))MySubsolver,&ctx);
105: PetscOpenMPFinalize();
106: MatDestroy(ctx.A);
107: VecDestroy(ctx.x);
108: VecDestroy(ctx.y);
109: VecDestroy(ctx.xr);
110: VecDestroy(ctx.yr);
111: VecScatterDestroy(ctx.sct);
113: PetscFinalize();
114: return 0;
115: }