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*/