xref: /petsc/src/dm/partitioner/impls/matpart/partmatpart.c (revision abe9303e6325b68e0d8957680d58b261e7a295d5)
1*abe9303eSLisandro Dalcin #include <petscmat.h>
2*abe9303eSLisandro Dalcin #include <petsc/private/partitionerimpl.h>   /*I      "petscpartitioner.h"   I*/
3*abe9303eSLisandro Dalcin 
4*abe9303eSLisandro Dalcin typedef struct {
5*abe9303eSLisandro Dalcin   MatPartitioning mp;
6*abe9303eSLisandro Dalcin } PetscPartitioner_MatPartitioning;
7*abe9303eSLisandro Dalcin 
8*abe9303eSLisandro Dalcin static PetscErrorCode PetscPartitionerMatPartitioningGetMatPartitioning_MatPartitioning(PetscPartitioner part, MatPartitioning *mp)
9*abe9303eSLisandro Dalcin {
10*abe9303eSLisandro Dalcin   PetscPartitioner_MatPartitioning  *p = (PetscPartitioner_MatPartitioning *) part->data;
11*abe9303eSLisandro Dalcin 
12*abe9303eSLisandro Dalcin   PetscFunctionBegin;
13*abe9303eSLisandro Dalcin   *mp = p->mp;
14*abe9303eSLisandro Dalcin   PetscFunctionReturn(0);
15*abe9303eSLisandro Dalcin }
16*abe9303eSLisandro Dalcin 
17*abe9303eSLisandro Dalcin /*@C
18*abe9303eSLisandro Dalcin   PetscPartitionerMatPartitioningGetMatPartitioning - Get a MatPartitioning instance wrapped by this PetscPartitioner.
19*abe9303eSLisandro Dalcin 
20*abe9303eSLisandro Dalcin   Not Collective
21*abe9303eSLisandro Dalcin 
22*abe9303eSLisandro Dalcin   Input Parameters:
23*abe9303eSLisandro Dalcin . part     - The PetscPartitioner
24*abe9303eSLisandro Dalcin 
25*abe9303eSLisandro Dalcin   Output Parameters:
26*abe9303eSLisandro Dalcin . mp       - The MatPartitioning
27*abe9303eSLisandro Dalcin 
28*abe9303eSLisandro Dalcin   Level: developer
29*abe9303eSLisandro Dalcin 
30*abe9303eSLisandro Dalcin .seealso DMPlexDistribute(), PetscPartitionerCreate()
31*abe9303eSLisandro Dalcin @*/
32*abe9303eSLisandro Dalcin PetscErrorCode PetscPartitionerMatPartitioningGetMatPartitioning(PetscPartitioner part, MatPartitioning *mp)
33*abe9303eSLisandro Dalcin {
34*abe9303eSLisandro Dalcin   PetscErrorCode          ierr;
35*abe9303eSLisandro Dalcin 
36*abe9303eSLisandro Dalcin   PetscFunctionBegin;
37*abe9303eSLisandro Dalcin   PetscValidHeaderSpecific(part, PETSCPARTITIONER_CLASSID, 1);
38*abe9303eSLisandro Dalcin   PetscValidPointer(mp, 2);
39*abe9303eSLisandro Dalcin   ierr = PetscUseMethod(part,"PetscPartitionerMatPartitioningGetMatPartitioning_C",(PetscPartitioner,MatPartitioning*),(part,mp));CHKERRQ(ierr);
40*abe9303eSLisandro Dalcin   PetscFunctionReturn(0);
41*abe9303eSLisandro Dalcin }
42*abe9303eSLisandro Dalcin 
43*abe9303eSLisandro Dalcin static PetscErrorCode PetscPartitionerDestroy_MatPartitioning(PetscPartitioner part)
44*abe9303eSLisandro Dalcin {
45*abe9303eSLisandro Dalcin   PetscPartitioner_MatPartitioning  *p = (PetscPartitioner_MatPartitioning *) part->data;
46*abe9303eSLisandro Dalcin   PetscErrorCode                    ierr;
47*abe9303eSLisandro Dalcin 
48*abe9303eSLisandro Dalcin   PetscFunctionBegin;
49*abe9303eSLisandro Dalcin   ierr = MatPartitioningDestroy(&p->mp);CHKERRQ(ierr);
50*abe9303eSLisandro Dalcin   ierr = PetscFree(part->data);CHKERRQ(ierr);
51*abe9303eSLisandro Dalcin   PetscFunctionReturn(0);
52*abe9303eSLisandro Dalcin }
53*abe9303eSLisandro Dalcin 
54*abe9303eSLisandro Dalcin static PetscErrorCode PetscPartitionerView_MatPartitioning_ASCII(PetscPartitioner part, PetscViewer viewer)
55*abe9303eSLisandro Dalcin {
56*abe9303eSLisandro Dalcin   PetscPartitioner_MatPartitioning  *p = (PetscPartitioner_MatPartitioning *) part->data;
57*abe9303eSLisandro Dalcin   PetscViewerFormat                 format;
58*abe9303eSLisandro Dalcin   PetscErrorCode                    ierr;
59*abe9303eSLisandro Dalcin 
60*abe9303eSLisandro Dalcin   PetscFunctionBegin;
61*abe9303eSLisandro Dalcin   ierr = PetscViewerGetFormat(viewer, &format);CHKERRQ(ierr);
62*abe9303eSLisandro Dalcin   ierr = PetscViewerASCIIPrintf(viewer, "MatPartitioning Graph Partitioner:\n");CHKERRQ(ierr);
63*abe9303eSLisandro Dalcin   ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr);
64*abe9303eSLisandro Dalcin   if (p->mp) {ierr = MatPartitioningView(p->mp, viewer);CHKERRQ(ierr);}
65*abe9303eSLisandro Dalcin   ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr);
66*abe9303eSLisandro Dalcin   PetscFunctionReturn(0);
67*abe9303eSLisandro Dalcin }
68*abe9303eSLisandro Dalcin 
69*abe9303eSLisandro Dalcin static PetscErrorCode PetscPartitionerView_MatPartitioning(PetscPartitioner part, PetscViewer viewer)
70*abe9303eSLisandro Dalcin {
71*abe9303eSLisandro Dalcin   PetscBool      iascii;
72*abe9303eSLisandro Dalcin   PetscErrorCode ierr;
73*abe9303eSLisandro Dalcin 
74*abe9303eSLisandro Dalcin   PetscFunctionBegin;
75*abe9303eSLisandro Dalcin   PetscValidHeaderSpecific(part, PETSCPARTITIONER_CLASSID, 1);
76*abe9303eSLisandro Dalcin   PetscValidHeaderSpecific(viewer, PETSC_VIEWER_CLASSID, 2);
77*abe9303eSLisandro Dalcin   ierr = PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);CHKERRQ(ierr);
78*abe9303eSLisandro Dalcin   if (iascii) {ierr = PetscPartitionerView_MatPartitioning_ASCII(part, viewer);CHKERRQ(ierr);}
79*abe9303eSLisandro Dalcin   PetscFunctionReturn(0);
80*abe9303eSLisandro Dalcin }
81*abe9303eSLisandro Dalcin 
82*abe9303eSLisandro Dalcin static PetscErrorCode PetscPartitionerSetFromOptions_MatPartitioning(PetscOptionItems *PetscOptionsObject, PetscPartitioner part)
83*abe9303eSLisandro Dalcin {
84*abe9303eSLisandro Dalcin   PetscPartitioner_MatPartitioning  *p = (PetscPartitioner_MatPartitioning *) part->data;
85*abe9303eSLisandro Dalcin   PetscErrorCode                    ierr;
86*abe9303eSLisandro Dalcin 
87*abe9303eSLisandro Dalcin   PetscFunctionBegin;
88*abe9303eSLisandro Dalcin   ierr = PetscObjectSetOptionsPrefix((PetscObject)p->mp,((PetscObject)part)->prefix);CHKERRQ(ierr);
89*abe9303eSLisandro Dalcin   ierr = MatPartitioningSetFromOptions(p->mp);CHKERRQ(ierr);
90*abe9303eSLisandro Dalcin   PetscFunctionReturn(0);
91*abe9303eSLisandro Dalcin }
92*abe9303eSLisandro Dalcin 
93*abe9303eSLisandro Dalcin static PetscErrorCode PetscPartitionerPartition_MatPartitioning(PetscPartitioner part, PetscInt nparts, PetscInt numVertices, PetscInt start[], PetscInt adjacency[], PetscSection vertSection, PetscSection targetSection, PetscSection partSection, IS *is)
94*abe9303eSLisandro Dalcin {
95*abe9303eSLisandro Dalcin   PetscPartitioner_MatPartitioning  *p = (PetscPartitioner_MatPartitioning *) part->data;
96*abe9303eSLisandro Dalcin   Mat                               matadj;
97*abe9303eSLisandro Dalcin   IS                                is1, is2, is3;
98*abe9303eSLisandro Dalcin   PetscReal                         *tpwgts = NULL;
99*abe9303eSLisandro Dalcin   PetscInt                          numVerticesGlobal, numEdges;
100*abe9303eSLisandro Dalcin   PetscInt                          *i, *j, *vwgt = NULL;
101*abe9303eSLisandro Dalcin   MPI_Comm                          comm;
102*abe9303eSLisandro Dalcin   PetscErrorCode                    ierr;
103*abe9303eSLisandro Dalcin 
104*abe9303eSLisandro Dalcin   PetscFunctionBegin;
105*abe9303eSLisandro Dalcin   ierr = PetscObjectGetComm((PetscObject)part, &comm);CHKERRQ(ierr);
106*abe9303eSLisandro Dalcin 
107*abe9303eSLisandro Dalcin   /* TODO: MatCreateMPIAdj should maybe take global number of ROWS */
108*abe9303eSLisandro Dalcin   /* TODO: And vertex distribution in PetscPartitionerPartition_ParMetis should be done using PetscSplitOwnership */
109*abe9303eSLisandro Dalcin   numVerticesGlobal = PETSC_DECIDE;
110*abe9303eSLisandro Dalcin   ierr = PetscSplitOwnership(comm, &numVertices, &numVerticesGlobal);CHKERRQ(ierr);
111*abe9303eSLisandro Dalcin 
112*abe9303eSLisandro Dalcin   /* copy arrays to avoid memory errors because MatMPIAdjSetPreallocation copies just pointers */
113*abe9303eSLisandro Dalcin   numEdges = start[numVertices];
114*abe9303eSLisandro Dalcin   ierr = PetscMalloc1(numVertices+1, &i);CHKERRQ(ierr);
115*abe9303eSLisandro Dalcin   ierr = PetscMalloc1(numEdges, &j);CHKERRQ(ierr);
116*abe9303eSLisandro Dalcin   ierr = PetscArraycpy(i, start, numVertices+1);CHKERRQ(ierr);
117*abe9303eSLisandro Dalcin   ierr = PetscArraycpy(j, adjacency, numEdges);CHKERRQ(ierr);
118*abe9303eSLisandro Dalcin 
119*abe9303eSLisandro Dalcin   /* construct the adjacency matrix */
120*abe9303eSLisandro Dalcin   ierr = MatCreateMPIAdj(comm, numVertices, numVerticesGlobal, i, j, NULL, &matadj);CHKERRQ(ierr);
121*abe9303eSLisandro Dalcin   ierr = MatPartitioningSetAdjacency(p->mp, matadj);CHKERRQ(ierr);
122*abe9303eSLisandro Dalcin   ierr = MatPartitioningSetNParts(p->mp, nparts);CHKERRQ(ierr);
123*abe9303eSLisandro Dalcin 
124*abe9303eSLisandro Dalcin   /* calculate partition weights */
125*abe9303eSLisandro Dalcin   if (targetSection) {
126*abe9303eSLisandro Dalcin     PetscReal sumt;
127*abe9303eSLisandro Dalcin     PetscInt  p;
128*abe9303eSLisandro Dalcin 
129*abe9303eSLisandro Dalcin     sumt = 0.0;
130*abe9303eSLisandro Dalcin     ierr = PetscMalloc1(nparts,&tpwgts);CHKERRQ(ierr);
131*abe9303eSLisandro Dalcin     for (p = 0; p < nparts; ++p) {
132*abe9303eSLisandro Dalcin       PetscInt tpd;
133*abe9303eSLisandro Dalcin 
134*abe9303eSLisandro Dalcin       ierr = PetscSectionGetDof(targetSection,p,&tpd);CHKERRQ(ierr);
135*abe9303eSLisandro Dalcin       sumt += tpd;
136*abe9303eSLisandro Dalcin       tpwgts[p] = tpd;
137*abe9303eSLisandro Dalcin     }
138*abe9303eSLisandro Dalcin     if (sumt) { /* METIS/ParMETIS do not like exactly zero weight */
139*abe9303eSLisandro Dalcin       for (p = 0, sumt = 0.0; p < nparts; ++p) {
140*abe9303eSLisandro Dalcin         tpwgts[p] = PetscMax(tpwgts[p],PETSC_SMALL);
141*abe9303eSLisandro Dalcin         sumt += tpwgts[p];
142*abe9303eSLisandro Dalcin       }
143*abe9303eSLisandro Dalcin       for (p = 0; p < nparts; ++p) tpwgts[p] /= sumt;
144*abe9303eSLisandro Dalcin       for (p = 0, sumt = 0.0; p < nparts-1; ++p) sumt += tpwgts[p];
145*abe9303eSLisandro Dalcin       tpwgts[nparts - 1] = 1. - sumt;
146*abe9303eSLisandro Dalcin     } else {
147*abe9303eSLisandro Dalcin       ierr = PetscFree(tpwgts);CHKERRQ(ierr);
148*abe9303eSLisandro Dalcin     }
149*abe9303eSLisandro Dalcin   }
150*abe9303eSLisandro Dalcin   ierr = MatPartitioningSetPartitionWeights(p->mp, tpwgts);CHKERRQ(ierr);
151*abe9303eSLisandro Dalcin 
152*abe9303eSLisandro Dalcin   /* calculate vertex weights */
153*abe9303eSLisandro Dalcin   if (vertSection) {
154*abe9303eSLisandro Dalcin     PetscInt v;
155*abe9303eSLisandro Dalcin 
156*abe9303eSLisandro Dalcin     ierr = PetscMalloc1(numVertices,&vwgt);CHKERRQ(ierr);
157*abe9303eSLisandro Dalcin     for (v = 0; v < numVertices; ++v) {
158*abe9303eSLisandro Dalcin       ierr = PetscSectionGetDof(vertSection, v, &vwgt[v]);CHKERRQ(ierr);
159*abe9303eSLisandro Dalcin     }
160*abe9303eSLisandro Dalcin   }
161*abe9303eSLisandro Dalcin   ierr = MatPartitioningSetVertexWeights(p->mp, vwgt);CHKERRQ(ierr);
162*abe9303eSLisandro Dalcin 
163*abe9303eSLisandro Dalcin   /* apply the partitioning */
164*abe9303eSLisandro Dalcin   ierr = MatPartitioningApply(p->mp, &is1);CHKERRQ(ierr);
165*abe9303eSLisandro Dalcin 
166*abe9303eSLisandro Dalcin   /* construct the PetscSection */
167*abe9303eSLisandro Dalcin   {
168*abe9303eSLisandro Dalcin     PetscInt v;
169*abe9303eSLisandro Dalcin     const PetscInt *assignment_arr;
170*abe9303eSLisandro Dalcin 
171*abe9303eSLisandro Dalcin     ierr = ISGetIndices(is1, &assignment_arr);CHKERRQ(ierr);
172*abe9303eSLisandro Dalcin     for (v = 0; v < numVertices; ++v) {ierr = PetscSectionAddDof(partSection, assignment_arr[v], 1);CHKERRQ(ierr);}
173*abe9303eSLisandro Dalcin     ierr = ISRestoreIndices(is1, &assignment_arr);CHKERRQ(ierr);
174*abe9303eSLisandro Dalcin   }
175*abe9303eSLisandro Dalcin 
176*abe9303eSLisandro Dalcin   /* convert assignment IS to global numbering IS */
177*abe9303eSLisandro Dalcin   ierr = ISPartitioningToNumbering(is1, &is2);CHKERRQ(ierr);
178*abe9303eSLisandro Dalcin   ierr = ISDestroy(&is1);CHKERRQ(ierr);
179*abe9303eSLisandro Dalcin 
180*abe9303eSLisandro Dalcin   /* renumber IS into local numbering */
181*abe9303eSLisandro Dalcin   ierr = ISOnComm(is2, PETSC_COMM_SELF, PETSC_USE_POINTER, &is1);CHKERRQ(ierr);
182*abe9303eSLisandro Dalcin   ierr = ISRenumber(is1, NULL, NULL, &is3);CHKERRQ(ierr);
183*abe9303eSLisandro Dalcin   ierr = ISDestroy(&is1);CHKERRQ(ierr);
184*abe9303eSLisandro Dalcin   ierr = ISDestroy(&is2);CHKERRQ(ierr);
185*abe9303eSLisandro Dalcin 
186*abe9303eSLisandro Dalcin   /* invert IS */
187*abe9303eSLisandro Dalcin   ierr = ISSetPermutation(is3);CHKERRQ(ierr);
188*abe9303eSLisandro Dalcin   ierr = ISInvertPermutation(is3, numVertices, &is1);CHKERRQ(ierr);
189*abe9303eSLisandro Dalcin   ierr = ISDestroy(&is3);CHKERRQ(ierr);
190*abe9303eSLisandro Dalcin 
191*abe9303eSLisandro Dalcin   ierr = MatDestroy(&matadj);CHKERRQ(ierr);
192*abe9303eSLisandro Dalcin   *is = is1;
193*abe9303eSLisandro Dalcin   PetscFunctionReturn(0);
194*abe9303eSLisandro Dalcin }
195*abe9303eSLisandro Dalcin 
196*abe9303eSLisandro Dalcin static PetscErrorCode PetscPartitionerInitialize_MatPartitioning(PetscPartitioner part)
197*abe9303eSLisandro Dalcin {
198*abe9303eSLisandro Dalcin   PetscErrorCode                    ierr;
199*abe9303eSLisandro Dalcin 
200*abe9303eSLisandro Dalcin   PetscFunctionBegin;
201*abe9303eSLisandro Dalcin   part->ops->view           = PetscPartitionerView_MatPartitioning;
202*abe9303eSLisandro Dalcin   part->ops->setfromoptions = PetscPartitionerSetFromOptions_MatPartitioning;
203*abe9303eSLisandro Dalcin   part->ops->destroy        = PetscPartitionerDestroy_MatPartitioning;
204*abe9303eSLisandro Dalcin   part->ops->partition      = PetscPartitionerPartition_MatPartitioning;
205*abe9303eSLisandro Dalcin   ierr = PetscObjectComposeFunction((PetscObject)part,"PetscPartitionerMatPartitioningGetMatPartitioning_C",PetscPartitionerMatPartitioningGetMatPartitioning_MatPartitioning);CHKERRQ(ierr);
206*abe9303eSLisandro Dalcin   PetscFunctionReturn(0);
207*abe9303eSLisandro Dalcin }
208*abe9303eSLisandro Dalcin 
209*abe9303eSLisandro Dalcin /*MC
210*abe9303eSLisandro Dalcin   PETSCPARTITIONERMATPARTITIONING = "matpartitioning" - A PetscPartitioner object
211*abe9303eSLisandro Dalcin 
212*abe9303eSLisandro Dalcin   Level: developer
213*abe9303eSLisandro Dalcin 
214*abe9303eSLisandro Dalcin .seealso: PetscPartitionerType, PetscPartitionerCreate(), PetscPartitionerSetType()
215*abe9303eSLisandro Dalcin M*/
216*abe9303eSLisandro Dalcin 
217*abe9303eSLisandro Dalcin PETSC_EXTERN PetscErrorCode PetscPartitionerCreate_MatPartitioning(PetscPartitioner part)
218*abe9303eSLisandro Dalcin {
219*abe9303eSLisandro Dalcin   PetscPartitioner_MatPartitioning  *p;
220*abe9303eSLisandro Dalcin   PetscErrorCode                    ierr;
221*abe9303eSLisandro Dalcin 
222*abe9303eSLisandro Dalcin   PetscFunctionBegin;
223*abe9303eSLisandro Dalcin   PetscValidHeaderSpecific(part, PETSCPARTITIONER_CLASSID, 1);
224*abe9303eSLisandro Dalcin   ierr       = PetscNewLog(part, &p);CHKERRQ(ierr);
225*abe9303eSLisandro Dalcin   part->data = p;
226*abe9303eSLisandro Dalcin   ierr = PetscPartitionerInitialize_MatPartitioning(part);CHKERRQ(ierr);
227*abe9303eSLisandro Dalcin   ierr = MatPartitioningCreate(PetscObjectComm((PetscObject)part), &p->mp);CHKERRQ(ierr);
228*abe9303eSLisandro Dalcin   PetscFunctionReturn(0);
229*abe9303eSLisandro Dalcin }
230*abe9303eSLisandro Dalcin 
231