Actual source code: ex4.c

  1: static char help[] = "This example tests subnetwork coupling without component. \n\n";

  3: #include <petscdmnetwork.h>

  5: int main(int argc,char ** argv)
  6: {
  8:   PetscMPIInt    size,rank;
  9:   DM             dmnetwork;
 10:   PetscInt       i,j,net,Nsubnet,ne,nv,nvar,v,goffset,row;
 11:   PetscInt       *numVertices,*numEdges,**edgelist,asvtx[2],bsvtx[2];
 12:   const PetscInt *vtx,*edges;
 13:   PetscBool      ghost,distribute=PETSC_TRUE;
 14:   Vec            X;
 15:   PetscScalar    val;

 17:   PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
 18:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
 19:   MPI_Comm_size(PETSC_COMM_WORLD,&size);

 21:   /* Create a network of subnetworks */
 22:   if (size == 1) Nsubnet = 2;
 23:   else Nsubnet = (PetscInt)size;
 24:   PetscCalloc3(Nsubnet,&numVertices,Nsubnet,&numEdges,Nsubnet,&edgelist);

 26:   /* when size>1, process[i] creates subnetwork[i] */
 27:   for (i=0; i<Nsubnet; i++) {
 28:     if (i == 0 && (size == 1 || (rank == i && size >1))) {
 29:       numVertices[i] = 4; numEdges[i] = 3;
 30:       PetscMalloc1(2*numEdges[i],&edgelist[i]);
 31:       edgelist[i][0] = 0; edgelist[i][1] = 1;
 32:       edgelist[i][2] = 1; edgelist[i][3] = 2;
 33:       edgelist[i][4] = 2; edgelist[i][5] = 3;

 35:     } else if (i == 1 && (size == 1 || (rank == i && size >1))) {
 36:       numVertices[i] = 4; numEdges[i] = 3;
 37:       PetscMalloc1(2*numEdges[i],&edgelist[i]);
 38:       edgelist[i][0] = 0; edgelist[i][1] = 1;
 39:       edgelist[i][2] = 1; edgelist[i][3] = 2;
 40:       edgelist[i][4] = 2; edgelist[i][5] = 3;

 42:     } else if (i>1 && (size == 1 || (rank == i && size >1))){
 43:       numVertices[i] = 4; numEdges[i] = 3;
 44:       PetscMalloc1(2*numEdges[i],&edgelist[i]);
 45:       for (j=0; j< numEdges[i]; j++) {
 46:         edgelist[i][2*j] = j; edgelist[i][2*j+1] = j+1;
 47:       }
 48:     }
 49:   }

 51:   /* Create a dmnetwork */
 52:   DMNetworkCreate(PETSC_COMM_WORLD,&dmnetwork);

 54:   /* Set number of subnetworks, numbers of vertices and edges over each subnetwork */
 55:   DMNetworkSetNumSubNetworks(dmnetwork,PETSC_DECIDE,Nsubnet);

 57:   for (i=0; i<Nsubnet; i++) {
 58:     PetscInt netNum = -1;
 59:     DMNetworkAddSubnetwork(dmnetwork,NULL,numVertices[i],numEdges[i],edgelist[i],&netNum);
 60:   }

 62:   /* Add shared vertices -- all processes hold this info at current implementation
 63:        net[0].0 -> net[j].0, j=0,...,Nsubnet-1
 64:        net[0].1 -> net[j].1, j=0,...,Nsubnet-1 */
 65:   asvtx[0] = bsvtx[0] = 0;
 66:   asvtx[1] = bsvtx[1] = 1;
 67:   for (j=Nsubnet-1; j>=1; j--) {
 68:     DMNetworkAddSharedVertices(dmnetwork,0,j,2,asvtx,bsvtx);
 69:   }

 71:   /* Setup the network layout */
 72:   DMNetworkLayoutSetUp(dmnetwork);

 74:   /* Get Subnetwork(); Add nvar=1 to subnet[0] and nvar=2 to other subnets */
 75:   for (net=0; net<Nsubnet; net++) {
 76:     DMNetworkGetSubnetwork(dmnetwork,net,&nv,&ne,&vtx,&edges);
 77:     for (v=0; v<nv; v++) {
 78:       if (!net) {
 79:         /* Set nvar = 2 for subnet0 */
 80:         DMNetworkAddComponent(dmnetwork,vtx[v],-1,NULL,2);
 81:       } else {
 82:         /* Set nvar = 1 for other subnets */
 83:         DMNetworkAddComponent(dmnetwork,vtx[v],-1,NULL,1);
 84:       }
 85:     }
 86:   }

 88:   /* Enable runtime option of graph partition type -- must be called before DMSetUp() */
 89:   if (size > 1) {
 90:     DM               plexdm;
 91:     PetscPartitioner part;
 92:     DMNetworkGetPlex(dmnetwork,&plexdm);
 93:     DMPlexGetPartitioner(plexdm, &part);
 94:     PetscPartitionerSetType(part,PETSCPARTITIONERSIMPLE);
 95:     PetscOptionsSetValue(NULL,"-dm_plex_csr_via_mat","true"); /* for parmetis */
 96:   }

 98:   /* Setup dmnetwork */
 99:   DMSetUp(dmnetwork);

101:   /* Redistribute the network layout; use '-distribute false' to skip */
102:   PetscOptionsGetBool(NULL,NULL,"-distribute",&distribute,NULL);
103:   if (distribute) {
104:     DMNetworkDistribute(&dmnetwork,0);
105:     DMView(dmnetwork,PETSC_VIEWER_STDOUT_WORLD);
106:   }

108:   /* Create a global vector */
109:   DMCreateGlobalVector(dmnetwork,&X);
110:   VecSet(X,0.0);

112:   /* Set X values at shared vertex */
113:   DMNetworkGetSharedVertices(dmnetwork,&nv,&vtx);
114:   for (v=0; v<nv; v++) {
115:     DMNetworkIsGhostVertex(dmnetwork,vtx[v],&ghost);
116:     if (ghost) continue;

118:     /* only one process holds a non-ghost vertex */
119:     DMNetworkGetComponent(dmnetwork,vtx[v],ALL_COMPONENTS,NULL,NULL,&nvar);
120:     DMNetworkGetGlobalVecOffset(dmnetwork,vtx[v],ALL_COMPONENTS,&goffset);
121:     for (i=0; i<nvar; i++) {
122:       row = goffset + i;
123:       val = (PetscScalar)rank + 1.0;
124:       VecSetValues(X,1,&row,&val,INSERT_VALUES);
125:     }
126:   }
127:   VecAssemblyBegin(X);
128:   VecAssemblyEnd(X);
129:   VecView(X,PETSC_VIEWER_STDOUT_WORLD);

131:   /* Free work space */
132:   VecDestroy(&X);
133:   for (i=0; i<Nsubnet; i++) {
134:     if (size == 1 || rank == i) {PetscFree(edgelist[i]);}
135:   }
136:   PetscFree3(numVertices,numEdges,edgelist);
137:   DMDestroy(&dmnetwork);
138:   PetscFinalize();
139:   return ierr;
140: }

142: /*TEST

144:    build:
145:       requires: !single double define(PETSC_HAVE_ATTRIBUTEALIGNED)

147:    test:
148:       args:

150:    test:
151:       suffix: 2
152:       nsize: 2
153:       args: -options_left no

155:    test:
156:       suffix: 3
157:       nsize: 4
158:       args: -options_left no

160: TEST*/