Actual source code: partmatpart.c

  1: #include <petscmat.h>
  2: #include <petsc/private/partitionerimpl.h>

  4: typedef struct {
  5:   MatPartitioning mp;
  6: } PetscPartitioner_MatPartitioning;

  8: static PetscErrorCode PetscPartitionerMatPartitioningGetMatPartitioning_MatPartitioning(PetscPartitioner part, MatPartitioning *mp)
  9: {
 10:   PetscPartitioner_MatPartitioning  *p = (PetscPartitioner_MatPartitioning *) part->data;

 12:   *mp = p->mp;
 13:   return 0;
 14: }

 16: /*@C
 17:   PetscPartitionerMatPartitioningGetMatPartitioning - Get a MatPartitioning instance wrapped by this PetscPartitioner.

 19:   Not Collective

 21:   Input Parameters:
 22: . part     - The PetscPartitioner

 24:   Output Parameters:
 25: . mp       - The MatPartitioning

 27:   Level: developer

 29: .seealso DMPlexDistribute(), PetscPartitionerCreate()
 30: @*/
 31: PetscErrorCode PetscPartitionerMatPartitioningGetMatPartitioning(PetscPartitioner part, MatPartitioning *mp)
 32: {
 35:   PetscUseMethod(part,"PetscPartitionerMatPartitioningGetMatPartitioning_C",(PetscPartitioner,MatPartitioning*),(part,mp));
 36:   return 0;
 37: }

 39: static PetscErrorCode PetscPartitionerDestroy_MatPartitioning(PetscPartitioner part)
 40: {
 41:   PetscPartitioner_MatPartitioning  *p = (PetscPartitioner_MatPartitioning *) part->data;

 43:   MatPartitioningDestroy(&p->mp);
 44:   PetscFree(part->data);
 45:   return 0;
 46: }

 48: static PetscErrorCode PetscPartitionerView_MatPartitioning_ASCII(PetscPartitioner part, PetscViewer viewer)
 49: {
 50:   PetscPartitioner_MatPartitioning  *p = (PetscPartitioner_MatPartitioning *) part->data;
 51:   PetscViewerFormat                 format;

 53:   PetscViewerGetFormat(viewer, &format);
 54:   PetscViewerASCIIPrintf(viewer, "MatPartitioning Graph Partitioner:\n");
 55:   PetscViewerASCIIPushTab(viewer);
 56:   if (p->mp) MatPartitioningView(p->mp, viewer);
 57:   PetscViewerASCIIPopTab(viewer);
 58:   return 0;
 59: }

 61: static PetscErrorCode PetscPartitionerView_MatPartitioning(PetscPartitioner part, PetscViewer viewer)
 62: {
 63:   PetscBool      iascii;

 67:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
 68:   if (iascii) PetscPartitionerView_MatPartitioning_ASCII(part, viewer);
 69:   return 0;
 70: }

 72: static PetscErrorCode PetscPartitionerSetFromOptions_MatPartitioning(PetscOptionItems *PetscOptionsObject, PetscPartitioner part)
 73: {
 74:   PetscPartitioner_MatPartitioning  *p = (PetscPartitioner_MatPartitioning *) part->data;

 76:   PetscObjectSetOptionsPrefix((PetscObject)p->mp,((PetscObject)part)->prefix);
 77:   MatPartitioningSetFromOptions(p->mp);
 78:   return 0;
 79: }

 81: static PetscErrorCode PetscPartitionerPartition_MatPartitioning(PetscPartitioner part, PetscInt nparts, PetscInt numVertices, PetscInt start[], PetscInt adjacency[], PetscSection vertSection, PetscSection targetSection, PetscSection partSection, IS *is)
 82: {
 83:   PetscPartitioner_MatPartitioning  *p = (PetscPartitioner_MatPartitioning *) part->data;
 84:   Mat                               matadj;
 85:   IS                                is1, is2, is3;
 86:   PetscReal                         *tpwgts = NULL;
 87:   PetscInt                          numVerticesGlobal, numEdges;
 88:   PetscInt                          *i, *j, *vwgt = NULL;
 89:   MPI_Comm                          comm;

 91:   PetscObjectGetComm((PetscObject)part, &comm);

 93:   /* TODO: MatCreateMPIAdj should maybe take global number of ROWS */
 94:   /* TODO: And vertex distribution in PetscPartitionerPartition_ParMetis should be done using PetscSplitOwnership */
 95:   numVerticesGlobal = PETSC_DECIDE;
 96:   PetscSplitOwnership(comm, &numVertices, &numVerticesGlobal);

 98:   /* copy arrays to avoid memory errors because MatMPIAdjSetPreallocation copies just pointers */
 99:   numEdges = start[numVertices];
100:   PetscMalloc1(numVertices+1, &i);
101:   PetscMalloc1(numEdges, &j);
102:   PetscArraycpy(i, start, numVertices+1);
103:   PetscArraycpy(j, adjacency, numEdges);

105:   /* construct the adjacency matrix */
106:   MatCreateMPIAdj(comm, numVertices, numVerticesGlobal, i, j, NULL, &matadj);
107:   MatPartitioningSetAdjacency(p->mp, matadj);
108:   MatPartitioningSetNParts(p->mp, nparts);

110:   /* calculate partition weights */
111:   if (targetSection) {
112:     PetscReal sumt;
113:     PetscInt  p;

115:     sumt = 0.0;
116:     PetscMalloc1(nparts,&tpwgts);
117:     for (p = 0; p < nparts; ++p) {
118:       PetscInt tpd;

120:       PetscSectionGetDof(targetSection,p,&tpd);
121:       sumt += tpd;
122:       tpwgts[p] = tpd;
123:     }
124:     if (sumt) { /* METIS/ParMETIS do not like exactly zero weight */
125:       for (p = 0, sumt = 0.0; p < nparts; ++p) {
126:         tpwgts[p] = PetscMax(tpwgts[p],PETSC_SMALL);
127:         sumt += tpwgts[p];
128:       }
129:       for (p = 0; p < nparts; ++p) tpwgts[p] /= sumt;
130:       for (p = 0, sumt = 0.0; p < nparts-1; ++p) sumt += tpwgts[p];
131:       tpwgts[nparts - 1] = 1. - sumt;
132:     } else {
133:       PetscFree(tpwgts);
134:     }
135:   }
136:   MatPartitioningSetPartitionWeights(p->mp, tpwgts);

138:   /* calculate vertex weights */
139:   if (vertSection) {
140:     PetscInt v;

142:     PetscMalloc1(numVertices,&vwgt);
143:     for (v = 0; v < numVertices; ++v) {
144:       PetscSectionGetDof(vertSection, v, &vwgt[v]);
145:     }
146:   }
147:   MatPartitioningSetVertexWeights(p->mp, vwgt);

149:   /* apply the partitioning */
150:   MatPartitioningApply(p->mp, &is1);

152:   /* construct the PetscSection */
153:   {
154:     PetscInt v;
155:     const PetscInt *assignment_arr;

157:     ISGetIndices(is1, &assignment_arr);
158:     for (v = 0; v < numVertices; ++v) PetscSectionAddDof(partSection, assignment_arr[v], 1);
159:     ISRestoreIndices(is1, &assignment_arr);
160:   }

162:   /* convert assignment IS to global numbering IS */
163:   ISPartitioningToNumbering(is1, &is2);
164:   ISDestroy(&is1);

166:   /* renumber IS into local numbering */
167:   ISOnComm(is2, PETSC_COMM_SELF, PETSC_USE_POINTER, &is1);
168:   ISRenumber(is1, NULL, NULL, &is3);
169:   ISDestroy(&is1);
170:   ISDestroy(&is2);

172:   /* invert IS */
173:   ISSetPermutation(is3);
174:   ISInvertPermutation(is3, numVertices, &is1);
175:   ISDestroy(&is3);

177:   MatDestroy(&matadj);
178:   *is = is1;
179:   return 0;
180: }

182: static PetscErrorCode PetscPartitionerInitialize_MatPartitioning(PetscPartitioner part)
183: {
184:   part->ops->view           = PetscPartitionerView_MatPartitioning;
185:   part->ops->setfromoptions = PetscPartitionerSetFromOptions_MatPartitioning;
186:   part->ops->destroy        = PetscPartitionerDestroy_MatPartitioning;
187:   part->ops->partition      = PetscPartitionerPartition_MatPartitioning;
188:   PetscObjectComposeFunction((PetscObject)part,"PetscPartitionerMatPartitioningGetMatPartitioning_C",PetscPartitionerMatPartitioningGetMatPartitioning_MatPartitioning);
189:   return 0;
190: }

192: /*MC
193:   PETSCPARTITIONERMATPARTITIONING = "matpartitioning" - A PetscPartitioner object

195:   Level: developer

197: .seealso: PetscPartitionerType, PetscPartitionerCreate(), PetscPartitionerSetType()
198: M*/

200: PETSC_EXTERN PetscErrorCode PetscPartitionerCreate_MatPartitioning(PetscPartitioner part)
201: {
202:   PetscPartitioner_MatPartitioning  *p;

205:   PetscNewLog(part, &p);
206:   part->data = p;
207:   PetscPartitionerInitialize_MatPartitioning(part);
208:   MatPartitioningCreate(PetscObjectComm((PetscObject)part), &p->mp);
209:   return 0;
210: }