1 #include <petsc/private/matimpl.h> /*I "petscmat.h" I*/
2 #include <../src/mat/impls/aij/seq/aij.h>
3 #include <../src/mat/impls/aij/mpi/mpiaij.h>
4 #include <petscsf.h>
5
6 typedef struct {
7 PetscBool symmetric;
8 } MC_Greedy;
9
MatColoringDestroy_Greedy(MatColoring mc)10 static PetscErrorCode MatColoringDestroy_Greedy(MatColoring mc)
11 {
12 PetscFunctionBegin;
13 PetscCall(PetscFree(mc->data));
14 PetscFunctionReturn(PETSC_SUCCESS);
15 }
16
GreedyColoringLocalDistanceOne_Private(MatColoring mc,PetscReal * wts,PetscInt * lperm,ISColoringValue * colors)17 static PetscErrorCode GreedyColoringLocalDistanceOne_Private(MatColoring mc, PetscReal *wts, PetscInt *lperm, ISColoringValue *colors)
18 {
19 PetscInt i, j, k, s, e, n, no, nd, nd_global, n_global, idx, ncols, maxcolors, masksize, ccol, *mask;
20 Mat m = mc->mat;
21 Mat_MPIAIJ *aij = (Mat_MPIAIJ *)m->data;
22 Mat md = NULL, mo = NULL;
23 const PetscInt *md_i, *mo_i, *md_j, *mo_j;
24 PetscBool isMPIAIJ, isSEQAIJ;
25 PetscInt pcol;
26 const PetscInt *cidx;
27 PetscInt *lcolors, *ocolors;
28 PetscReal *owts = NULL;
29 PetscSF sf;
30 PetscLayout layout;
31
32 PetscFunctionBegin;
33 PetscCall(MatGetSize(m, &n_global, NULL));
34 PetscCall(MatGetOwnershipRange(m, &s, &e));
35 n = e - s;
36 masksize = 20;
37 nd_global = 0;
38 /* get the matrix communication structures */
39 PetscCall(PetscObjectTypeCompare((PetscObject)m, MATMPIAIJ, &isMPIAIJ));
40 PetscCall(PetscObjectBaseTypeCompare((PetscObject)m, MATSEQAIJ, &isSEQAIJ));
41 PetscCheck(isMPIAIJ || isSEQAIJ, PetscObjectComm((PetscObject)mc), PETSC_ERR_ARG_WRONG, "Matrix must be AIJ for greedy coloring");
42 if (isMPIAIJ) {
43 /* get the CSR data for on and off-diagonal portions of m */
44 Mat_SeqAIJ *dseq;
45 Mat_SeqAIJ *oseq;
46 md = aij->A;
47 dseq = (Mat_SeqAIJ *)md->data;
48 mo = aij->B;
49 oseq = (Mat_SeqAIJ *)mo->data;
50 md_i = dseq->i;
51 md_j = dseq->j;
52 mo_i = oseq->i;
53 mo_j = oseq->j;
54 } else {
55 /* get the CSR data for m */
56 Mat_SeqAIJ *dseq;
57 /* no off-processor nodes */
58 md = m;
59 dseq = (Mat_SeqAIJ *)md->data;
60 mo = NULL;
61 no = 0;
62 md_i = dseq->i;
63 md_j = dseq->j;
64 mo_i = NULL;
65 mo_j = NULL;
66 }
67
68 PetscCall(MatColoringGetMaxColors(mc, &maxcolors));
69 if (mo) {
70 PetscCall(VecGetSize(aij->lvec, &no));
71 PetscCall(PetscMalloc2(no, &ocolors, no, &owts));
72 for (i = 0; i < no; i++) ocolors[i] = maxcolors;
73 }
74
75 PetscCall(PetscMalloc1(masksize, &mask));
76 PetscCall(PetscMalloc1(n, &lcolors));
77 for (i = 0; i < n; i++) lcolors[i] = maxcolors;
78 for (i = 0; i < masksize; i++) mask[i] = -1;
79 if (mo) {
80 /* transfer neighbor weights */
81 PetscCall(PetscSFCreate(PetscObjectComm((PetscObject)m), &sf));
82 PetscCall(MatGetLayouts(m, &layout, NULL));
83 PetscCall(PetscSFSetGraphLayout(sf, layout, no, NULL, PETSC_COPY_VALUES, aij->garray));
84 PetscCall(PetscSFBcastBegin(sf, MPIU_REAL, wts, owts, MPI_REPLACE));
85 PetscCall(PetscSFBcastEnd(sf, MPIU_REAL, wts, owts, MPI_REPLACE));
86 }
87 while (nd_global < n_global) {
88 nd = n;
89 /* assign lowest possible color to each local vertex */
90 PetscCall(PetscLogEventBegin(MATCOLORING_Local, mc, 0, 0, 0));
91 for (i = 0; i < n; i++) {
92 idx = lperm[i];
93 if (lcolors[idx] == maxcolors) {
94 ncols = md_i[idx + 1] - md_i[idx];
95 cidx = &(md_j[md_i[idx]]);
96 for (j = 0; j < ncols; j++) {
97 if (lcolors[cidx[j]] != maxcolors) {
98 ccol = lcolors[cidx[j]];
99 if (ccol >= masksize) {
100 PetscInt *newmask;
101 PetscCall(PetscMalloc1(masksize * 2, &newmask));
102 for (k = 0; k < 2 * masksize; k++) newmask[k] = -1;
103 for (k = 0; k < masksize; k++) newmask[k] = mask[k];
104 PetscCall(PetscFree(mask));
105 mask = newmask;
106 masksize *= 2;
107 }
108 mask[ccol] = idx;
109 }
110 }
111 if (mo) {
112 ncols = mo_i[idx + 1] - mo_i[idx];
113 cidx = &(mo_j[mo_i[idx]]);
114 for (j = 0; j < ncols; j++) {
115 if (ocolors[cidx[j]] != maxcolors) {
116 ccol = ocolors[cidx[j]];
117 if (ccol >= masksize) {
118 PetscInt *newmask;
119 PetscCall(PetscMalloc1(masksize * 2, &newmask));
120 for (k = 0; k < 2 * masksize; k++) newmask[k] = -1;
121 for (k = 0; k < masksize; k++) newmask[k] = mask[k];
122 PetscCall(PetscFree(mask));
123 mask = newmask;
124 masksize *= 2;
125 }
126 mask[ccol] = idx;
127 }
128 }
129 }
130 for (j = 0; j < masksize; j++) {
131 if (mask[j] != idx) break;
132 }
133 pcol = j;
134 if (pcol > maxcolors) pcol = maxcolors;
135 lcolors[idx] = pcol;
136 }
137 }
138 PetscCall(PetscLogEventEnd(MATCOLORING_Local, mc, 0, 0, 0));
139 if (mo) {
140 /* transfer neighbor colors */
141 PetscCall(PetscLogEventBegin(MATCOLORING_Comm, mc, 0, 0, 0));
142 PetscCall(PetscSFBcastBegin(sf, MPIU_INT, lcolors, ocolors, MPI_REPLACE));
143 PetscCall(PetscSFBcastEnd(sf, MPIU_INT, lcolors, ocolors, MPI_REPLACE));
144 /* check for conflicts -- this is merely checking if any adjacent off-processor rows have the same color and marking the ones that are lower weight locally for changing */
145 for (i = 0; i < n; i++) {
146 ncols = mo_i[i + 1] - mo_i[i];
147 cidx = &(mo_j[mo_i[i]]);
148 for (j = 0; j < ncols; j++) {
149 /* in the case of conflicts, the highest weight one stays and the others go */
150 if ((ocolors[cidx[j]] == lcolors[i]) && (owts[cidx[j]] > wts[i]) && lcolors[i] < maxcolors) {
151 lcolors[i] = maxcolors;
152 nd--;
153 }
154 }
155 }
156 nd_global = 0;
157 }
158 PetscCallMPI(MPIU_Allreduce(&nd, &nd_global, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)mc)));
159 }
160 for (i = 0; i < n; i++) colors[i] = (ISColoringValue)lcolors[i];
161 PetscCall(PetscFree(mask));
162 PetscCall(PetscFree(lcolors));
163 if (mo) {
164 PetscCall(PetscFree2(ocolors, owts));
165 PetscCall(PetscSFDestroy(&sf));
166 }
167 PetscFunctionReturn(PETSC_SUCCESS);
168 }
169
GreedyColoringLocalDistanceTwo_Private(MatColoring mc,PetscReal * wts,PetscInt * lperm,ISColoringValue * colors)170 static PetscErrorCode GreedyColoringLocalDistanceTwo_Private(MatColoring mc, PetscReal *wts, PetscInt *lperm, ISColoringValue *colors)
171 {
172 MC_Greedy *gr = (MC_Greedy *)mc->data;
173 PetscInt i, j, k, l, s, e, n, nd, nd_global, n_global, idx, ncols, maxcolors, mcol, mcol_global, nd1cols, *mask, masksize, *d1cols, *bad, *badnext, nbad, badsize, ccol, no, cbad;
174 Mat m = mc->mat, mt;
175 Mat_MPIAIJ *aij = (Mat_MPIAIJ *)m->data;
176 Mat md = NULL, mo = NULL;
177 const PetscInt *md_i, *mo_i, *md_j, *mo_j;
178 const PetscInt *rmd_i, *rmo_i, *rmd_j, *rmo_j;
179 PetscBool isMPIAIJ, isSEQAIJ;
180 PetscInt pcol, *dcolors, *ocolors;
181 ISColoringValue *badidx;
182 const PetscInt *cidx;
183 PetscReal *owts, *colorweights;
184 PetscInt *oconf, *conf;
185 PetscSF sf;
186 PetscLayout layout;
187
188 PetscFunctionBegin;
189 PetscCall(MatGetSize(m, &n_global, NULL));
190 PetscCall(MatGetOwnershipRange(m, &s, &e));
191 n = e - s;
192 nd_global = 0;
193 /* get the matrix communication structures */
194 PetscCall(PetscObjectBaseTypeCompare((PetscObject)m, MATMPIAIJ, &isMPIAIJ));
195 PetscCall(PetscObjectBaseTypeCompare((PetscObject)m, MATSEQAIJ, &isSEQAIJ));
196 PetscCheck(isMPIAIJ || isSEQAIJ, PetscObjectComm((PetscObject)mc), PETSC_ERR_ARG_WRONG, "Matrix must be AIJ for greedy coloring");
197 if (isMPIAIJ) {
198 Mat_SeqAIJ *dseq;
199 Mat_SeqAIJ *oseq;
200 md = aij->A;
201 dseq = (Mat_SeqAIJ *)md->data;
202 mo = aij->B;
203 oseq = (Mat_SeqAIJ *)mo->data;
204 md_i = dseq->i;
205 md_j = dseq->j;
206 mo_i = oseq->i;
207 mo_j = oseq->j;
208 rmd_i = dseq->i;
209 rmd_j = dseq->j;
210 rmo_i = oseq->i;
211 rmo_j = oseq->j;
212 } else {
213 Mat_SeqAIJ *dseq;
214 /* no off-processor nodes */
215 md = m;
216 dseq = (Mat_SeqAIJ *)md->data;
217 md_i = dseq->i;
218 md_j = dseq->j;
219 mo_i = NULL;
220 mo_j = NULL;
221 rmd_i = dseq->i;
222 rmd_j = dseq->j;
223 rmo_i = NULL;
224 rmo_j = NULL;
225 }
226 if (!gr->symmetric) {
227 Mat_SeqAIJ *dseq = NULL;
228
229 PetscCheck(isSEQAIJ, PetscObjectComm((PetscObject)mc), PETSC_ERR_SUP, "Nonsymmetric greedy coloring only works in serial");
230 PetscCall(MatTranspose(m, MAT_INITIAL_MATRIX, &mt));
231 dseq = (Mat_SeqAIJ *)mt->data;
232 rmd_i = dseq->i;
233 rmd_j = dseq->j;
234 rmo_i = NULL;
235 rmo_j = NULL;
236 }
237 /* create the vectors and communication structures if necessary */
238 no = 0;
239 if (mo) {
240 PetscCall(VecGetLocalSize(aij->lvec, &no));
241 PetscCall(PetscSFCreate(PetscObjectComm((PetscObject)m), &sf));
242 PetscCall(MatGetLayouts(m, &layout, NULL));
243 PetscCall(PetscSFSetGraphLayout(sf, layout, no, NULL, PETSC_COPY_VALUES, aij->garray));
244 }
245 PetscCall(MatColoringGetMaxColors(mc, &maxcolors));
246 masksize = n;
247 nbad = 0;
248 badsize = n;
249 PetscCall(PetscMalloc1(masksize, &mask));
250 PetscCall(PetscMalloc4(n, &d1cols, n, &dcolors, n, &conf, n, &bad));
251 PetscCall(PetscMalloc2(badsize, &badidx, badsize, &badnext));
252 for (i = 0; i < masksize; i++) mask[i] = -1;
253 for (i = 0; i < n; i++) {
254 dcolors[i] = maxcolors;
255 bad[i] = -1;
256 }
257 for (i = 0; i < badsize; i++) badnext[i] = -1;
258 if (mo) {
259 PetscCall(PetscMalloc3(no, &owts, no, &oconf, no, &ocolors));
260 PetscCall(PetscSFBcastBegin(sf, MPIU_REAL, wts, owts, MPI_REPLACE));
261 PetscCall(PetscSFBcastEnd(sf, MPIU_REAL, wts, owts, MPI_REPLACE));
262 for (i = 0; i < no; i++) ocolors[i] = maxcolors;
263 } else { /* Appease overzealous -Wmaybe-initialized */
264 owts = NULL;
265 oconf = NULL;
266 ocolors = NULL;
267 }
268 mcol = 0;
269 while (nd_global < n_global) {
270 nd = n;
271 /* assign lowest possible color to each local vertex */
272 mcol_global = 0;
273 PetscCall(PetscLogEventBegin(MATCOLORING_Local, mc, 0, 0, 0));
274 for (i = 0; i < n; i++) {
275 idx = lperm[i];
276 if (dcolors[idx] == maxcolors) {
277 /* entries in bad */
278 cbad = bad[idx];
279 while (cbad >= 0) {
280 ccol = badidx[cbad];
281 if (ccol >= masksize) {
282 PetscInt *newmask;
283 PetscCall(PetscMalloc1(masksize * 2, &newmask));
284 for (k = 0; k < 2 * masksize; k++) newmask[k] = -1;
285 for (k = 0; k < masksize; k++) newmask[k] = mask[k];
286 PetscCall(PetscFree(mask));
287 mask = newmask;
288 masksize *= 2;
289 }
290 mask[ccol] = idx;
291 cbad = badnext[cbad];
292 }
293 /* diagonal distance-one rows */
294 nd1cols = 0;
295 ncols = rmd_i[idx + 1] - rmd_i[idx];
296 cidx = &(rmd_j[rmd_i[idx]]);
297 for (j = 0; j < ncols; j++) {
298 d1cols[nd1cols] = cidx[j];
299 nd1cols++;
300 ccol = dcolors[cidx[j]];
301 if (ccol != maxcolors) {
302 if (ccol >= masksize) {
303 PetscInt *newmask;
304 PetscCall(PetscMalloc1(masksize * 2, &newmask));
305 for (k = 0; k < 2 * masksize; k++) newmask[k] = -1;
306 for (k = 0; k < masksize; k++) newmask[k] = mask[k];
307 PetscCall(PetscFree(mask));
308 mask = newmask;
309 masksize *= 2;
310 }
311 mask[ccol] = idx;
312 }
313 }
314 /* off-diagonal distance-one rows */
315 if (mo) {
316 ncols = rmo_i[idx + 1] - rmo_i[idx];
317 cidx = &(rmo_j[rmo_i[idx]]);
318 for (j = 0; j < ncols; j++) {
319 ccol = ocolors[cidx[j]];
320 if (ccol != maxcolors) {
321 if (ccol >= masksize) {
322 PetscInt *newmask;
323 PetscCall(PetscMalloc1(masksize * 2, &newmask));
324 for (k = 0; k < 2 * masksize; k++) newmask[k] = -1;
325 for (k = 0; k < masksize; k++) newmask[k] = mask[k];
326 PetscCall(PetscFree(mask));
327 mask = newmask;
328 masksize *= 2;
329 }
330 mask[ccol] = idx;
331 }
332 }
333 }
334 /* diagonal distance-two rows */
335 for (j = 0; j < nd1cols; j++) {
336 ncols = md_i[d1cols[j] + 1] - md_i[d1cols[j]];
337 cidx = &(md_j[md_i[d1cols[j]]]);
338 for (l = 0; l < ncols; l++) {
339 ccol = dcolors[cidx[l]];
340 if (ccol != maxcolors) {
341 if (ccol >= masksize) {
342 PetscInt *newmask;
343 PetscCall(PetscMalloc1(masksize * 2, &newmask));
344 for (k = 0; k < 2 * masksize; k++) newmask[k] = -1;
345 for (k = 0; k < masksize; k++) newmask[k] = mask[k];
346 PetscCall(PetscFree(mask));
347 mask = newmask;
348 masksize *= 2;
349 }
350 mask[ccol] = idx;
351 }
352 }
353 }
354 /* off-diagonal distance-two rows */
355 if (mo) {
356 for (j = 0; j < nd1cols; j++) {
357 ncols = mo_i[d1cols[j] + 1] - mo_i[d1cols[j]];
358 cidx = &(mo_j[mo_i[d1cols[j]]]);
359 for (l = 0; l < ncols; l++) {
360 ccol = ocolors[cidx[l]];
361 if (ccol != maxcolors) {
362 if (ccol >= masksize) {
363 PetscInt *newmask;
364 PetscCall(PetscMalloc1(masksize * 2, &newmask));
365 for (k = 0; k < 2 * masksize; k++) newmask[k] = -1;
366 for (k = 0; k < masksize; k++) newmask[k] = mask[k];
367 PetscCall(PetscFree(mask));
368 mask = newmask;
369 masksize *= 2;
370 }
371 mask[ccol] = idx;
372 }
373 }
374 }
375 }
376 /* assign this one the lowest color possible by seeing if there's a gap in the sequence of sorted neighbor colors */
377 for (j = 0; j < masksize; j++) {
378 if (mask[j] != idx) break;
379 }
380 pcol = j;
381 if (pcol > maxcolors) pcol = maxcolors;
382 dcolors[idx] = pcol;
383 if (pcol > mcol) mcol = pcol;
384 }
385 }
386 PetscCall(PetscLogEventEnd(MATCOLORING_Local, mc, 0, 0, 0));
387 if (mo) {
388 /* transfer neighbor colors */
389 PetscCall(PetscSFBcastBegin(sf, MPIU_INT, dcolors, ocolors, MPI_REPLACE));
390 PetscCall(PetscSFBcastEnd(sf, MPIU_INT, dcolors, ocolors, MPI_REPLACE));
391 /* find the maximum color assigned locally and allocate a mask */
392 PetscCallMPI(MPIU_Allreduce(&mcol, &mcol_global, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)mc)));
393 PetscCall(PetscMalloc1(mcol_global + 1, &colorweights));
394 /* check for conflicts */
395 for (i = 0; i < n; i++) conf[i] = PETSC_FALSE;
396 for (i = 0; i < no; i++) oconf[i] = PETSC_FALSE;
397 for (i = 0; i < n; i++) {
398 ncols = mo_i[i + 1] - mo_i[i];
399 cidx = &(mo_j[mo_i[i]]);
400 if (ncols > 0) {
401 /* fill in the mask */
402 for (j = 0; j < mcol_global + 1; j++) colorweights[j] = 0;
403 colorweights[dcolors[i]] = wts[i];
404 /* fill in the off-diagonal part of the mask */
405 for (j = 0; j < ncols; j++) {
406 ccol = ocolors[cidx[j]];
407 if (ccol < maxcolors) {
408 if (colorweights[ccol] < owts[cidx[j]]) colorweights[ccol] = owts[cidx[j]];
409 }
410 }
411 /* fill in the on-diagonal part of the mask */
412 ncols = md_i[i + 1] - md_i[i];
413 cidx = &(md_j[md_i[i]]);
414 for (j = 0; j < ncols; j++) {
415 ccol = dcolors[cidx[j]];
416 if (ccol < maxcolors) {
417 if (colorweights[ccol] < wts[cidx[j]]) colorweights[ccol] = wts[cidx[j]];
418 }
419 }
420 /* go back through and set up on and off-diagonal conflict vectors */
421 ncols = md_i[i + 1] - md_i[i];
422 cidx = &(md_j[md_i[i]]);
423 for (j = 0; j < ncols; j++) {
424 ccol = dcolors[cidx[j]];
425 if (ccol < maxcolors) {
426 if (colorweights[ccol] > wts[cidx[j]]) conf[cidx[j]] = PETSC_TRUE;
427 }
428 }
429 ncols = mo_i[i + 1] - mo_i[i];
430 cidx = &(mo_j[mo_i[i]]);
431 for (j = 0; j < ncols; j++) {
432 ccol = ocolors[cidx[j]];
433 if (ccol < maxcolors) {
434 if (colorweights[ccol] > owts[cidx[j]]) oconf[cidx[j]] = PETSC_TRUE;
435 }
436 }
437 }
438 }
439 nd_global = 0;
440 PetscCall(PetscFree(colorweights));
441 PetscCall(PetscLogEventBegin(MATCOLORING_Comm, mc, 0, 0, 0));
442 PetscCall(PetscSFReduceBegin(sf, MPIU_INT, oconf, conf, MPI_SUM));
443 PetscCall(PetscSFReduceEnd(sf, MPIU_INT, oconf, conf, MPI_SUM));
444 PetscCall(PetscLogEventEnd(MATCOLORING_Comm, mc, 0, 0, 0));
445 /* go through and unset local colors that have conflicts */
446 for (i = 0; i < n; i++) {
447 if (conf[i] > 0) {
448 /* push this color onto the bad stack */
449 PetscCall(ISColoringValueCast(dcolors[i], &badidx[nbad]));
450 badnext[nbad] = bad[i];
451 bad[i] = nbad;
452 nbad++;
453 if (nbad >= badsize) {
454 PetscInt *newbadnext;
455 ISColoringValue *newbadidx;
456 PetscCall(PetscMalloc2(badsize * 2, &newbadidx, badsize * 2, &newbadnext));
457 for (k = 0; k < 2 * badsize; k++) newbadnext[k] = -1;
458 for (k = 0; k < badsize; k++) {
459 newbadidx[k] = badidx[k];
460 newbadnext[k] = badnext[k];
461 }
462 PetscCall(PetscFree2(badidx, badnext));
463 badidx = newbadidx;
464 badnext = newbadnext;
465 badsize *= 2;
466 }
467 dcolors[i] = maxcolors;
468 nd--;
469 }
470 }
471 }
472 PetscCallMPI(MPIU_Allreduce(&nd, &nd_global, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)mc)));
473 }
474 if (mo) {
475 PetscCall(PetscSFDestroy(&sf));
476 PetscCall(PetscFree3(owts, oconf, ocolors));
477 }
478 for (i = 0; i < n; i++) PetscCall(ISColoringValueCast(dcolors[i], colors + i));
479 PetscCall(PetscFree(mask));
480 PetscCall(PetscFree4(d1cols, dcolors, conf, bad));
481 PetscCall(PetscFree2(badidx, badnext));
482 if (!gr->symmetric) PetscCall(MatDestroy(&mt));
483 PetscFunctionReturn(PETSC_SUCCESS);
484 }
485
MatColoringApply_Greedy(MatColoring mc,ISColoring * iscoloring)486 static PetscErrorCode MatColoringApply_Greedy(MatColoring mc, ISColoring *iscoloring)
487 {
488 PetscInt finalcolor, finalcolor_global;
489 ISColoringValue *colors;
490 PetscInt ncolstotal, ncols;
491 PetscReal *wts;
492 PetscInt i, *lperm;
493
494 PetscFunctionBegin;
495 PetscCall(MatGetSize(mc->mat, NULL, &ncolstotal));
496 PetscCall(MatGetLocalSize(mc->mat, NULL, &ncols));
497 if (!mc->user_weights) {
498 PetscCall(MatColoringCreateWeights(mc, &wts, &lperm));
499 } else {
500 wts = mc->user_weights;
501 lperm = mc->user_lperm;
502 }
503 PetscCheck(mc->dist == 1 || mc->dist == 2, PetscObjectComm((PetscObject)mc), PETSC_ERR_ARG_OUTOFRANGE, "Only distance 1 and distance 2 supported by MatColoringGreedy");
504 PetscCall(PetscMalloc1(ncols, &colors));
505 if (mc->dist == 1) {
506 PetscCall(GreedyColoringLocalDistanceOne_Private(mc, wts, lperm, colors));
507 } else {
508 PetscCall(GreedyColoringLocalDistanceTwo_Private(mc, wts, lperm, colors));
509 }
510 finalcolor = 0;
511 for (i = 0; i < ncols; i++) {
512 if (colors[i] > finalcolor) finalcolor = colors[i];
513 }
514 finalcolor_global = 0;
515 PetscCallMPI(MPIU_Allreduce(&finalcolor, &finalcolor_global, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)mc)));
516 PetscCall(PetscLogEventBegin(MATCOLORING_ISCreate, mc, 0, 0, 0));
517 PetscCall(ISColoringCreate(PetscObjectComm((PetscObject)mc), finalcolor_global + 1, ncols, colors, PETSC_OWN_POINTER, iscoloring));
518 PetscCall(PetscLogEventEnd(MATCOLORING_ISCreate, mc, 0, 0, 0));
519 if (!mc->user_weights) {
520 PetscCall(PetscFree(wts));
521 PetscCall(PetscFree(lperm));
522 }
523 PetscFunctionReturn(PETSC_SUCCESS);
524 }
525
MatColoringSetFromOptions_Greedy(MatColoring mc,PetscOptionItems PetscOptionsObject)526 static PetscErrorCode MatColoringSetFromOptions_Greedy(MatColoring mc, PetscOptionItems PetscOptionsObject)
527 {
528 MC_Greedy *gr = (MC_Greedy *)mc->data;
529
530 PetscFunctionBegin;
531 PetscOptionsHeadBegin(PetscOptionsObject, "Greedy options");
532 PetscCall(PetscOptionsBool("-mat_coloring_greedy_symmetric", "Flag for assuming a symmetric matrix", "", gr->symmetric, &gr->symmetric, NULL));
533 PetscOptionsHeadEnd();
534 PetscFunctionReturn(PETSC_SUCCESS);
535 }
536
537 /*MC
538 MATCOLORINGGREEDY - Greedy-with-conflict correction based matrix coloring for distance 1 and 2 {cite}`bozdaug2005parallel`
539
540 Level: beginner
541
542 Notes:
543 These algorithms proceed in two phases -- local coloring and conflict resolution. The local coloring
544 tentatively colors all vertices at the distance required given what's known of the global coloring. Then,
545 the updated colors are transferred to different processors at distance one. In the distance one case, each
546 vertex with nonlocal neighbors is then checked to see if it conforms, with the vertex being
547 marked for recoloring if its lower weight than its same colored neighbor. In the distance two case,
548 each boundary vertex's immediate star is checked for validity of the coloring. Lower-weight conflict
549 vertices are marked, and then the conflicts are gathered back on owning processors. In both cases
550 this is done until each column has received a valid color.
551
552 Supports both distance one and distance two colorings.
553
554 .seealso: `MatColoringType`, `MatColoringCreate()`, `MatColoring`, `MatColoringSetType()`, `MatColoringType`
555 M*/
MatColoringCreate_Greedy(MatColoring mc)556 PETSC_EXTERN PetscErrorCode MatColoringCreate_Greedy(MatColoring mc)
557 {
558 MC_Greedy *gr;
559
560 PetscFunctionBegin;
561 PetscCall(PetscNew(&gr));
562 mc->data = gr;
563 mc->ops->apply = MatColoringApply_Greedy;
564 mc->ops->view = NULL;
565 mc->ops->destroy = MatColoringDestroy_Greedy;
566 mc->ops->setfromoptions = MatColoringSetFromOptions_Greedy;
567
568 gr->symmetric = PETSC_TRUE;
569 PetscFunctionReturn(PETSC_SUCCESS);
570 }
571