Actual source code: ex3.c

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

      4: /* T
      5:   Concepts: DMNetwork
      6: */
      7: #include <petscdmnetwork.h>

      9: typedef struct{
     10:   PetscInt id;
     11: } Comp0;

     13: typedef struct{
     14:   PetscScalar val;
     15: } Comp1;

     17: int main(int argc,char ** argv)
     18: {
     20:   PetscMPIInt    size,rank;
     21:   DM             dmnetwork;
     22:   PetscInt       i,j,net,Nsubnet,nsubnet,ne,nv,nvar,v,ncomp,compkey0,compkey1,compkey,goffset,row;
     23:   PetscInt       numVertices[10],numEdges[10],*edgelist[10],asvtx,bsvtx;
     24:   const PetscInt *vtx,*edges;
     25:   PetscBool      sharedv,ghost,distribute=PETSC_TRUE,test=PETSC_FALSE;
     26:   Vec            X;
     27:   Comp0          comp0;
     28:   Comp1          comp1;
     29:   PetscScalar    val;

     31:   PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
     32:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
     33:   MPI_Comm_size(PETSC_COMM_WORLD,&size);

     35:   /* Create a network of subnetworks */
     36:   nsubnet = 1;
     37:   if (size == 1) nsubnet = 2;

     39:   /* Create a dmnetwork and register components */
     40:   DMNetworkCreate(PETSC_COMM_WORLD,&dmnetwork);
     41:   DMNetworkRegisterComponent(dmnetwork,"comp0",sizeof(Comp0),&compkey0);
     42:   DMNetworkRegisterComponent(dmnetwork,"comp1",sizeof(Comp1),&compkey1);

     44:   /* Set componnet values - intentionally take rank-dependent value for test */
     45:   comp0.id  = rank;
     46:   comp1.val = 10.0*rank;

     48:   /* Set number of subnetworks, numbers of vertices and edges over each subnetwork */
     49:   DMNetworkSetNumSubNetworks(dmnetwork,nsubnet,PETSC_DECIDE);
     50:   DMNetworkGetNumSubNetworks(dmnetwork,NULL,&Nsubnet);

     52:   /* Input subnetworks; when size>1, process[i] creates subnetwork[i] */
     53:   for (i=0; i<Nsubnet; i++) {numVertices[i] = 0; numEdges[i] = 0;}
     54:   for (i=0; i<Nsubnet; i++) {
     55:     if (i == 0 && (size == 1 || (rank == i && size >1))) {
     56:       numVertices[i] = 4; numEdges[i] = 3;
     57:       PetscMalloc1(2*numEdges[i],&edgelist[i]);
     58:       edgelist[i][0] = 0; edgelist[i][1] = 2;
     59:       edgelist[i][2] = 2; edgelist[i][3] = 1;
     60:       edgelist[i][4] = 1; edgelist[i][5] = 3;

     62:     } else if (i == 1 && (size == 1 || (rank == i && size >1))) {
     63:       numVertices[i] = 4; numEdges[i] = 3;
     64:       PetscMalloc1(2*numEdges[i],&edgelist[i]);
     65:       edgelist[i][0] = 0; edgelist[i][1] = 3;
     66:       edgelist[i][2] = 3; edgelist[i][3] = 2;
     67:       edgelist[i][4] = 2; edgelist[i][5] = 1;

     69:     } else if (i>1 && (size == 1 || (rank == i && size >1))){
     70:       numVertices[i] = 4; numEdges[i] = 3;
     71:       PetscMalloc1(2*numEdges[i],&edgelist[i]);
     72:       for (j=0; j< numEdges[i]; j++) {
     73:         edgelist[i][2*j] = j; edgelist[i][2*j+1] = j+1;
     74:       }
     75:     }
     76:   }

     78:   /* Add subnetworks */
     79:   for (i=0; i<Nsubnet; i++) {
     80:     PetscInt netNum = -1;
     81:     DMNetworkAddSubnetwork(dmnetwork,NULL,numVertices[i],numEdges[i],edgelist[i],&netNum);
     82:   }

     84:   /* Add shared vertices -- all processes hold this info at current implementation */
     85:   asvtx = bsvtx = 0;
     86:   for (j=1; j<Nsubnet; j++) {
     87:     /* vertex subnet[0].0 shares with vertex subnet[j].0 */
     88:     DMNetworkAddSharedVertices(dmnetwork,0,j,1,&asvtx,&bsvtx);
     89:   }

     91:   /* Setup the network layout */
     92:   DMNetworkLayoutSetUp(dmnetwork);

     94:   /* Get Subnetwork(); Add nvar=1 to subnet[0] and nvar=2 to other subnets */
     95:   for (net=0; net<Nsubnet; net++) {
     96:     DMNetworkGetSubnetwork(dmnetwork,net,&nv,&ne,&vtx,&edges);
     97:     for (v=0; v<nv; v++) {
     98:       DMNetworkIsSharedVertex(dmnetwork,vtx[v],&sharedv);
     99:       if (sharedv) {
    100:         #if 0
    101:         /* current version requirs all processess add componenets and nvar at the shared vertices! */
    102:         if (net == 0) {
    103:           //printf("[%d] net %d, v[%d]=%d is a shared vertex, add comp0[0]\n",rank,net,v,vtx[v]);
    104:           DMNetworkAddComponent(dmnetwork,vtx[v],compkey0,&comp0,1);
    105:         } else if (net == 1) {
    106:           //printf("[%d] net %d, v[%d]=%d is a shared vertex, add comp1\n",rank,net,v,vtx[v]);
    107:           DMNetworkAddComponent(dmnetwork,vtx[v],compkey1,&comp1,2);
    108:         }
    109:         #endif
    110:         continue;
    111:       }

    113:       if (!net) {
    114:         DMNetworkAddComponent(dmnetwork,vtx[v],compkey0,&comp0,1);
    115:       } else {
    116:         DMNetworkAddComponent(dmnetwork,vtx[v],compkey1,&comp1,2);
    117:       }
    118:     }
    119:   }

    121:   /* Add components and nvar to shared vertex -- only owner of the vertex does this! */
    122:   DMNetworkGetSharedVertices(dmnetwork,&nv,&vtx);
    123:   for (v=0; v<nv; v++) {
    124:     DMNetworkIsGhostVertex(dmnetwork,vtx[v],&ghost);
    125:     if (ghost) continue;
    126:     DMNetworkAddComponent(dmnetwork,vtx[v],compkey0,&comp0,1);
    127:     DMNetworkAddComponent(dmnetwork,vtx[v],compkey1,&comp1,2);
    128:   }

    130:   /* Enable runtime option of graph partition type -- must be called before DMSetUp() */
    131:   if (size > 1) {
    132:     DM               plexdm;
    133:     PetscPartitioner part;
    134:     DMNetworkGetPlex(dmnetwork,&plexdm);
    135:     DMPlexGetPartitioner(plexdm, &part);
    136:     PetscPartitionerSetType(part,PETSCPARTITIONERSIMPLE);
    137:     PetscOptionsSetValue(NULL,"-dm_plex_csr_via_mat","true"); /* for parmetis */
    138:   }

    140:   /* Setup dmnetwork */
    141:   DMSetUp(dmnetwork);

    143:   /* Redistribute the network layout; use '-distribute false' to skip */
    144:   PetscOptionsGetBool(NULL,NULL,"-distribute",&distribute,NULL);
    145:   if (distribute) {
    146:     DMNetworkDistribute(&dmnetwork,0);
    147:     DMView(dmnetwork,PETSC_VIEWER_STDOUT_WORLD);
    148:   }

    150:   /* Create a global vector */
    151:   DMCreateGlobalVector(dmnetwork,&X);
    152:   VecSet(X,0.0);

    154:   /* Set X values at shared vertex */
    155:   DMNetworkGetSharedVertices(dmnetwork,&nv,&vtx);
    156:   for (v=0; v<nv; v++) {
    157:     DMNetworkIsGhostVertex(dmnetwork,vtx[v],&ghost);
    158:     if (ghost) continue;

    160:     /* only one process holds a non-ghost vertex */
    161:     DMNetworkGetComponent(dmnetwork,vtx[v],ALL_COMPONENTS,NULL,NULL,&nvar);
    162:     DMNetworkGetNumComponents(dmnetwork,vtx[v],&ncomp);
    163:     /* PetscPrintf(PETSC_COMM_SELF,"[%d] shared v %D: nvar %D, ncomp %D\n",rank,vtx[v],nvar,ncomp); */
    164:     for (j=0; j<ncomp; j++) {
    165:       DMNetworkGetComponent(dmnetwork,vtx[v],j,&compkey,NULL,&nvar);
    166:       DMNetworkGetGlobalVecOffset(dmnetwork,vtx[v],j,&goffset);
    167:       for (i=0; i<nvar; i++) {
    168:         row = goffset + i;
    169:         val = compkey + 1.0;
    170:         VecSetValues(X,1,&row,&val,INSERT_VALUES);
    171:       }
    172:     }
    173:   }
    174:   VecAssemblyBegin(X);
    175:   VecAssemblyEnd(X);
    176:   VecView(X,PETSC_VIEWER_STDOUT_WORLD);

    178:   /* Test DMNetworkGetSubnetwork() */
    179:   PetscOptionsGetBool(NULL,NULL,"-test_getsubnet",&test,NULL);
    180:   if (test) {
    181:     net = 0;
    182:     PetscOptionsGetInt(NULL,NULL,"-subnet",&net,NULL);
    183:     DMNetworkGetSubnetwork(dmnetwork,net,&nv,&ne,&vtx,&edges);
    184:     PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d] subnet %D: nv %D, ne %D\n",rank,net,nv,ne);
    185:     PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT);
    186:     MPI_Barrier(PETSC_COMM_WORLD);

    188:     for (i=0; i<nv; i++) {
    189:       DMNetworkIsGhostVertex(dmnetwork,vtx[i],&ghost);
    190:       DMNetworkIsSharedVertex(dmnetwork,vtx[i],&sharedv);

    192:       DMNetworkGetNumComponents(dmnetwork,vtx[i],&ncomp);
    193:       if (sharedv || ghost) {
    194:         PetscPrintf(PETSC_COMM_SELF,"  [%d] v %D is shared %d, is ghost %d, ncomp %D\n",rank,vtx[i],sharedv,ghost,ncomp);
    195:       }

    197:       for (j=0; j<ncomp; j++) {
    198:         void* component;
    199:         DMNetworkGetComponent(dmnetwork,vtx[i],j,&compkey,(void**)&component,NULL);
    200:         if (compkey == 0) {
    201:           Comp0  *mycomp0;
    202:           mycomp0 = (Comp0*)component;
    203:           PetscPrintf(PETSC_COMM_SELF,"  [%d] v %D compkey %D, mycomp0->id %D\n",rank,vtx[i],compkey,mycomp0->id);
    204:         } else if (compkey == 1) {
    205:           Comp1  *mycomp1;
    206:           mycomp1 = (Comp1*)component;
    207:           PetscPrintf(PETSC_COMM_SELF,"  [%d] v %D compkey %D, mycomp1->val %g\n",rank,vtx[i],compkey,mycomp1->val);
    208:         }
    209:       }
    210:     }
    211:   }

    213:   /* Free work space */
    214:   VecDestroy(&X);
    215:   for (i=0; i<Nsubnet; i++) {
    216:     if (size == 1 || rank == i) {PetscFree(edgelist[i]);}
    217:   }

    219:   DMDestroy(&dmnetwork);
    220:   PetscFinalize();
    221:   return ierr;
    222: }

    224: /*TEST

    226:    build:
    227:       requires: !single double define(PETSC_HAVE_ATTRIBUTEALIGNED)

    229:    test:
    230:       args:

    232:    test:
    233:       suffix: 2
    234:       nsize: 2
    235:       args: -options_left no

    237:    test:
    238:       suffix: 3
    239:       nsize: 4
    240:       args: -options_left no

    242: TEST*/