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 <petscdm.h>
5
6 /* linked list methods
7 *
8 * PetscCDCreate
9 */
PetscCDCreate(PetscInt a_size,PetscCoarsenData ** a_out)10 PetscErrorCode PetscCDCreate(PetscInt a_size, PetscCoarsenData **a_out)
11 {
12 PetscCoarsenData *ail;
13
14 PetscFunctionBegin;
15 /* allocate pool, partially */
16 PetscCall(PetscNew(&ail));
17 *a_out = ail;
18 ail->pool_list.next = NULL;
19 ail->pool_list.array = NULL;
20 ail->chk_sz = 0;
21 /* allocate array */
22 ail->size = a_size;
23 PetscCall(PetscCalloc1(a_size, &ail->array));
24 ail->extra_nodes = NULL;
25 ail->mat = NULL;
26 PetscFunctionReturn(PETSC_SUCCESS);
27 }
28
29 /* PetscCDDestroy
30 */
PetscCDDestroy(PetscCoarsenData * ail)31 PetscErrorCode PetscCDDestroy(PetscCoarsenData *ail)
32 {
33 PetscCDArrNd *n = &ail->pool_list;
34
35 PetscFunctionBegin;
36 n = n->next;
37 while (n) {
38 PetscCDArrNd *lstn = n;
39
40 n = n->next;
41 PetscCall(PetscFree(lstn));
42 }
43 if (ail->pool_list.array) PetscCall(PetscFree(ail->pool_list.array));
44 PetscCall(PetscFree(ail->array));
45 if (ail->mat) PetscCall(MatDestroy(&ail->mat));
46 /* delete this (+agg+pool array) */
47 PetscCall(PetscFree(ail));
48 PetscFunctionReturn(PETSC_SUCCESS);
49 }
50
51 /* PetscCDSetChunkSize
52 */
PetscCDSetChunkSize(PetscCoarsenData * ail,PetscInt a_sz)53 PetscErrorCode PetscCDSetChunkSize(PetscCoarsenData *ail, PetscInt a_sz)
54 {
55 PetscFunctionBegin;
56 ail->chk_sz = a_sz;
57 PetscFunctionReturn(PETSC_SUCCESS);
58 }
59
60 /* PetscCDGetNewNode
61 */
PetscCDGetNewNode(PetscCoarsenData * ail,PetscCDIntNd ** a_out,PetscInt a_id)62 static PetscErrorCode PetscCDGetNewNode(PetscCoarsenData *ail, PetscCDIntNd **a_out, PetscInt a_id)
63 {
64 PetscFunctionBegin;
65 *a_out = NULL; /* squelch -Wmaybe-uninitialized */
66 if (ail->extra_nodes) {
67 PetscCDIntNd *node = ail->extra_nodes;
68
69 ail->extra_nodes = node->next;
70 node->gid = a_id;
71 node->next = NULL;
72 *a_out = node;
73 } else {
74 if (!ail->pool_list.array) {
75 if (!ail->chk_sz) ail->chk_sz = 10; /* use a chuck size of ail->size? */
76 PetscCall(PetscMalloc1(ail->chk_sz, &ail->pool_list.array));
77 ail->new_node = ail->pool_list.array;
78 ail->new_left = ail->chk_sz;
79 ail->new_node->next = NULL;
80 } else if (!ail->new_left) {
81 PetscCDArrNd *node;
82
83 PetscCall(PetscMalloc(ail->chk_sz * sizeof(PetscCDIntNd) + sizeof(PetscCDArrNd), &node));
84 node->array = (PetscCDIntNd *)(node + 1);
85 node->next = ail->pool_list.next;
86 ail->pool_list.next = node;
87 ail->new_left = ail->chk_sz;
88 ail->new_node = node->array;
89 }
90 ail->new_node->gid = a_id;
91 ail->new_node->next = NULL;
92 *a_out = ail->new_node++;
93 ail->new_left--;
94 }
95 PetscFunctionReturn(PETSC_SUCCESS);
96 }
97
98 /* PetscCDIntNdSetID
99 */
PetscCDIntNdSetID(PetscCDIntNd * a_this,PetscInt a_id)100 PetscErrorCode PetscCDIntNdSetID(PetscCDIntNd *a_this, PetscInt a_id)
101 {
102 PetscFunctionBegin;
103 a_this->gid = a_id;
104 PetscFunctionReturn(PETSC_SUCCESS);
105 }
106
107 /* PetscCDIntNdGetID
108 */
PetscCDIntNdGetID(const PetscCDIntNd * a_this,PetscInt * a_gid)109 PetscErrorCode PetscCDIntNdGetID(const PetscCDIntNd *a_this, PetscInt *a_gid)
110 {
111 PetscFunctionBegin;
112 *a_gid = a_this->gid;
113 PetscFunctionReturn(PETSC_SUCCESS);
114 }
115
116 /* PetscCDGetHeadPos
117 */
PetscCDGetHeadPos(const PetscCoarsenData * ail,PetscInt a_idx,PetscCDIntNd ** pos)118 PetscErrorCode PetscCDGetHeadPos(const PetscCoarsenData *ail, PetscInt a_idx, PetscCDIntNd **pos)
119 {
120 PetscFunctionBegin;
121 PetscCheck(a_idx < ail->size, PETSC_COMM_SELF, PETSC_ERR_PLIB, "a_idx >= ail->size: a_idx=%" PetscInt_FMT ".", a_idx);
122 *pos = ail->array[a_idx];
123 PetscFunctionReturn(PETSC_SUCCESS);
124 }
125
126 /* PetscCDGetNextPos
127 */
PetscCDGetNextPos(const PetscCoarsenData * ail,PetscInt l_idx,PetscCDIntNd ** pos)128 PetscErrorCode PetscCDGetNextPos(const PetscCoarsenData *ail, PetscInt l_idx, PetscCDIntNd **pos)
129 {
130 PetscFunctionBegin;
131 PetscCheck(*pos, PETSC_COMM_SELF, PETSC_ERR_PLIB, "NULL input position.");
132 *pos = (*pos)->next;
133 PetscFunctionReturn(PETSC_SUCCESS);
134 }
135
136 /* PetscCDAppendID
137 */
PetscCDAppendID(PetscCoarsenData * ail,PetscInt a_idx,PetscInt a_id)138 PetscErrorCode PetscCDAppendID(PetscCoarsenData *ail, PetscInt a_idx, PetscInt a_id)
139 {
140 PetscCDIntNd *n, *n2;
141
142 PetscFunctionBegin;
143 PetscCall(PetscCDGetNewNode(ail, &n, a_id));
144 PetscCheck(a_idx < ail->size, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Index %" PetscInt_FMT " out of range.", a_idx);
145 if (!(n2 = ail->array[a_idx])) ail->array[a_idx] = n;
146 else {
147 do {
148 if (!n2->next) {
149 n2->next = n;
150 PetscCheck(!n->next, PETSC_COMM_SELF, PETSC_ERR_PLIB, "n should not have a next");
151 break;
152 }
153 n2 = n2->next;
154 } while (n2);
155 PetscCheck(n2, PETSC_COMM_SELF, PETSC_ERR_PLIB, "n2 should be non-null");
156 }
157 PetscFunctionReturn(PETSC_SUCCESS);
158 }
159
160 /* PetscCDAppendNode
161 */
PetscCDAppendNode(PetscCoarsenData * ail,PetscInt a_idx,PetscCDIntNd * a_n)162 PetscErrorCode PetscCDAppendNode(PetscCoarsenData *ail, PetscInt a_idx, PetscCDIntNd *a_n)
163 {
164 PetscCDIntNd *n2;
165
166 PetscFunctionBegin;
167 PetscCheck(a_idx < ail->size, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Index %" PetscInt_FMT " out of range.", a_idx);
168 if (!(n2 = ail->array[a_idx])) ail->array[a_idx] = a_n;
169 else {
170 do {
171 if (!n2->next) {
172 n2->next = a_n;
173 a_n->next = NULL;
174 break;
175 }
176 n2 = n2->next;
177 } while (n2);
178 PetscCheck(n2, PETSC_COMM_SELF, PETSC_ERR_PLIB, "n2 should be non-null");
179 }
180 PetscFunctionReturn(PETSC_SUCCESS);
181 }
182
183 /* PetscCDRemoveNextNode: a_last->next, this exposes single linked list structure to API (not used)
184 */
PetscCDRemoveNextNode(PetscCoarsenData * ail,PetscInt a_idx,PetscCDIntNd * a_last)185 PetscErrorCode PetscCDRemoveNextNode(PetscCoarsenData *ail, PetscInt a_idx, PetscCDIntNd *a_last)
186 {
187 PetscCDIntNd *del;
188
189 PetscFunctionBegin;
190 PetscCheck(a_idx < ail->size, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Index %" PetscInt_FMT " out of range.", a_idx);
191 PetscCheck(a_last->next, PETSC_COMM_SELF, PETSC_ERR_PLIB, "a_last should have a next");
192 del = a_last->next;
193 a_last->next = del->next;
194 /* del->next = NULL; -- this still used in a iterator so keep it intact -- need to fix this with a double linked list */
195 /* could reuse n2 but PetscCDAppendNode sometimes uses it */
196 PetscFunctionReturn(PETSC_SUCCESS);
197 }
198
199 /* PetscCDPrint
200 */
PetscCDPrint(const PetscCoarsenData * ail,PetscInt Istart,MPI_Comm comm)201 PetscErrorCode PetscCDPrint(const PetscCoarsenData *ail, PetscInt Istart, MPI_Comm comm)
202 {
203 PetscCDIntNd *n, *n2;
204 PetscInt ii;
205
206 PetscFunctionBegin;
207 for (ii = 0; ii < ail->size; ii++) {
208 n2 = n = ail->array[ii];
209 if (n) PetscCall(PetscSynchronizedPrintf(comm, "list %" PetscInt_FMT ":", ii + Istart));
210 while (n) {
211 PetscCall(PetscSynchronizedPrintf(comm, " %" PetscInt_FMT, n->gid));
212 n = n->next;
213 }
214 if (n2) PetscCall(PetscSynchronizedPrintf(comm, "\n"));
215 }
216 PetscCall(PetscSynchronizedFlush(comm, PETSC_STDOUT));
217 PetscFunctionReturn(PETSC_SUCCESS);
218 }
219
220 /* PetscCDMoveAppend - take list in a_srcidx and appends to destidx
221 */
PetscCDMoveAppend(PetscCoarsenData * ail,PetscInt a_destidx,PetscInt a_srcidx)222 PetscErrorCode PetscCDMoveAppend(PetscCoarsenData *ail, PetscInt a_destidx, PetscInt a_srcidx)
223 {
224 PetscCDIntNd *n;
225
226 PetscFunctionBegin;
227 PetscCheck(a_srcidx < ail->size, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Index %" PetscInt_FMT " out of range.", a_srcidx);
228 PetscCheck(a_destidx < ail->size, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Index %" PetscInt_FMT " out of range.", a_destidx);
229 PetscCheck(a_destidx != a_srcidx, PETSC_COMM_SELF, PETSC_ERR_PLIB, "a_destidx==a_srcidx %" PetscInt_FMT ".", a_destidx);
230 n = ail->array[a_destidx];
231 if (!n) ail->array[a_destidx] = ail->array[a_srcidx];
232 else {
233 do {
234 if (!n->next) {
235 n->next = ail->array[a_srcidx]; // append
236 break;
237 }
238 n = n->next;
239 } while (1);
240 }
241 ail->array[a_srcidx] = NULL; // empty
242 PetscFunctionReturn(PETSC_SUCCESS);
243 }
244
245 /* PetscCDRemoveAllAt - empty one list and move data to cache
246 */
PetscCDRemoveAllAt(PetscCoarsenData * ail,PetscInt a_idx)247 PetscErrorCode PetscCDRemoveAllAt(PetscCoarsenData *ail, PetscInt a_idx)
248 {
249 PetscCDIntNd *rem, *n1;
250
251 PetscFunctionBegin;
252 PetscCheck(a_idx < ail->size, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Index %" PetscInt_FMT " out of range.", a_idx);
253 rem = ail->array[a_idx];
254 ail->array[a_idx] = NULL;
255 if (!(n1 = ail->extra_nodes)) ail->extra_nodes = rem;
256 else {
257 while (n1->next) n1 = n1->next;
258 n1->next = rem;
259 }
260 PetscFunctionReturn(PETSC_SUCCESS);
261 }
262
263 /* PetscCDCountAt
264 */
PetscCDCountAt(const PetscCoarsenData * ail,PetscInt a_idx,PetscInt * a_sz)265 PetscErrorCode PetscCDCountAt(const PetscCoarsenData *ail, PetscInt a_idx, PetscInt *a_sz)
266 {
267 PetscCDIntNd *n1;
268 PetscInt sz = 0;
269
270 PetscFunctionBegin;
271 PetscCheck(a_idx < ail->size, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Index %" PetscInt_FMT " out of range.", a_idx);
272 n1 = ail->array[a_idx];
273 while (n1) {
274 n1 = n1->next;
275 sz++;
276 }
277 *a_sz = sz;
278 PetscFunctionReturn(PETSC_SUCCESS);
279 }
280
281 /* PetscCDSize
282 */
PetscCDCount(const PetscCoarsenData * ail,PetscInt * a_sz)283 PetscErrorCode PetscCDCount(const PetscCoarsenData *ail, PetscInt *a_sz)
284 {
285 PetscInt sz = 0;
286
287 PetscFunctionBegin;
288 for (PetscInt ii = 0; ii < ail->size; ii++) {
289 PetscCDIntNd *n1 = ail->array[ii];
290
291 while (n1) {
292 n1 = n1->next;
293 sz++;
294 }
295 }
296 *a_sz = sz;
297 PetscFunctionReturn(PETSC_SUCCESS);
298 }
299
300 /* PetscCDIsEmptyAt - Is the list empty? (not used)
301 */
PetscCDIsEmptyAt(const PetscCoarsenData * ail,PetscInt a_idx,PetscBool * a_e)302 PetscErrorCode PetscCDIsEmptyAt(const PetscCoarsenData *ail, PetscInt a_idx, PetscBool *a_e)
303 {
304 PetscFunctionBegin;
305 PetscCheck(a_idx < ail->size, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Index %" PetscInt_FMT " out of range.", a_idx);
306 *a_e = (PetscBool)(ail->array[a_idx] == NULL);
307 PetscFunctionReturn(PETSC_SUCCESS);
308 }
309
310 /* PetscCDGetNonemptyIS - used for C-F methods
311 */
PetscCDGetNonemptyIS(PetscCoarsenData * ail,IS * a_mis)312 PetscErrorCode PetscCDGetNonemptyIS(PetscCoarsenData *ail, IS *a_mis)
313 {
314 PetscCDIntNd *n;
315 PetscInt ii, kk;
316 PetscInt *permute;
317
318 PetscFunctionBegin;
319 for (ii = kk = 0; ii < ail->size; ii++) {
320 n = ail->array[ii];
321 if (n) kk++;
322 }
323 PetscCall(PetscMalloc1(kk, &permute));
324 for (ii = kk = 0; ii < ail->size; ii++) {
325 n = ail->array[ii];
326 if (n) permute[kk++] = ii;
327 }
328 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, kk, permute, PETSC_OWN_POINTER, a_mis));
329 PetscFunctionReturn(PETSC_SUCCESS);
330 }
331
332 /* PetscCDGetMat
333 */
PetscCDGetMat(PetscCoarsenData * ail,Mat * a_mat)334 PetscErrorCode PetscCDGetMat(PetscCoarsenData *ail, Mat *a_mat)
335 {
336 PetscFunctionBegin;
337 *a_mat = ail->mat;
338 PetscFunctionReturn(PETSC_SUCCESS);
339 }
340
341 /* PetscCDSetMat
342 */
PetscCDSetMat(PetscCoarsenData * ail,Mat a_mat)343 PetscErrorCode PetscCDSetMat(PetscCoarsenData *ail, Mat a_mat)
344 {
345 PetscFunctionBegin;
346 if (ail->mat) PetscCall(MatDestroy(&ail->mat)); //should not happen
347 ail->mat = a_mat;
348 PetscFunctionReturn(PETSC_SUCCESS);
349 }
350
351 /* PetscCDClearMat
352 */
PetscCDClearMat(PetscCoarsenData * ail)353 PetscErrorCode PetscCDClearMat(PetscCoarsenData *ail)
354 {
355 PetscFunctionBegin;
356 ail->mat = NULL;
357 PetscFunctionReturn(PETSC_SUCCESS);
358 }
359
360 /* PetscCDGetASMBlocks - get IS of aggregates for ASM smoothers
361 */
PetscCDGetASMBlocks(const PetscCoarsenData * ail,const PetscInt a_bs,PetscInt * a_sz,IS ** a_local_is)362 PetscErrorCode PetscCDGetASMBlocks(const PetscCoarsenData *ail, const PetscInt a_bs, PetscInt *a_sz, IS **a_local_is)
363 {
364 PetscCDIntNd *n;
365 PetscInt lsz, ii, kk, *idxs, jj, gid;
366 IS *is_loc = NULL;
367
368 PetscFunctionBegin;
369 for (ii = kk = 0; ii < ail->size; ii++) {
370 if (ail->array[ii]) kk++;
371 }
372 *a_sz = kk;
373 PetscCall(PetscMalloc1(kk, &is_loc));
374 for (ii = kk = 0; ii < ail->size; ii++) {
375 for (lsz = 0, n = ail->array[ii]; n; lsz++, n = n->next) /* void */
376 ;
377 if (lsz) {
378 PetscCall(PetscMalloc1(a_bs * lsz, &idxs));
379 for (lsz = 0, n = ail->array[ii]; n; n = n->next) {
380 PetscCall(PetscCDIntNdGetID(n, &gid));
381 for (jj = 0; jj < a_bs; lsz++, jj++) idxs[lsz] = a_bs * gid + jj;
382 }
383 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, lsz, idxs, PETSC_OWN_POINTER, &is_loc[kk++]));
384 }
385 }
386 PetscCheck(*a_sz == kk, PETSC_COMM_SELF, PETSC_ERR_PLIB, "*a_sz %" PetscInt_FMT " != kk %" PetscInt_FMT, *a_sz, kk);
387 *a_local_is = is_loc; /* out */
388 PetscFunctionReturn(PETSC_SUCCESS);
389 }
390
391 /* edge for priority queue */
392 typedef struct edge_tag {
393 PetscReal weight;
394 PetscInt lid0, gid1, ghost1_idx;
395 } Edge;
396
397 #define MY_MEPS (PETSC_MACHINE_EPSILON * 100)
gamg_hem_compare(const void * a,const void * b)398 static int gamg_hem_compare(const void *a, const void *b)
399 {
400 PetscReal va = ((Edge *)a)->weight, vb = ((Edge *)b)->weight;
401 return (va <= vb - MY_MEPS) ? 1 : (va > vb + MY_MEPS) ? -1 : 0; /* 0 for equal */
402 }
403
404 /*
405 MatCoarsenApply_HEM_private - parallel heavy edge matching
406
407 Input Parameter:
408 . a_Gmat - global matrix of the graph
409 . n_iter - number of matching iterations
410 . threshold - threshold for filtering graphs
411
412 Output Parameter:
413 . a_locals_llist - array of list of local nodes rooted at local node
414 */
MatCoarsenApply_HEM_private(Mat a_Gmat,const PetscInt n_iter,const PetscReal threshold,PetscCoarsenData ** a_locals_llist)415 static PetscErrorCode MatCoarsenApply_HEM_private(Mat a_Gmat, const PetscInt n_iter, const PetscReal threshold, PetscCoarsenData **a_locals_llist)
416 {
417 #define REQ_BF_SIZE 100
418 PetscBool isMPI;
419 MPI_Comm comm;
420 PetscInt ix, *ii, *aj, Istart, bc_agg = -1, *rbuff = NULL, rbuff_sz = 0;
421 PetscMPIInt rank, size, comm_procs[REQ_BF_SIZE], ncomm_procs, *lid_max_pe;
422 const PetscInt nloc = a_Gmat->rmap->n, request_size = PetscCeilInt((int)sizeof(MPI_Request), (int)sizeof(PetscInt));
423 PetscInt *lid_cprowID;
424 PetscBool *lid_matched;
425 Mat_SeqAIJ *matA, *matB = NULL;
426 Mat_MPIAIJ *mpimat = NULL;
427 PetscScalar one = 1.;
428 PetscCoarsenData *agg_llists = NULL, *ghost_deleted_list = NULL, *bc_list = NULL;
429 Mat cMat, tMat, P;
430 MatScalar *ap;
431 IS info_is;
432
433 PetscFunctionBegin;
434 PetscCall(PetscObjectGetComm((PetscObject)a_Gmat, &comm));
435 PetscCallMPI(MPI_Comm_rank(comm, &rank));
436 PetscCallMPI(MPI_Comm_size(comm, &size));
437 PetscCall(MatGetOwnershipRange(a_Gmat, &Istart, NULL));
438 PetscCall(ISCreate(comm, &info_is));
439 PetscCall(PetscInfo(info_is, "Start %" PetscInt_FMT " iterations of HEM.\n", n_iter));
440
441 PetscCall(PetscMalloc3(nloc, &lid_matched, nloc, &lid_cprowID, nloc, &lid_max_pe));
442 PetscCall(PetscCDCreate(nloc, &agg_llists));
443 PetscCall(PetscCDSetChunkSize(agg_llists, nloc + 1));
444 *a_locals_llist = agg_llists;
445 /* add self to all lists */
446 for (PetscInt kk = 0; kk < nloc; kk++) PetscCall(PetscCDAppendID(agg_llists, kk, Istart + kk));
447 /* make a copy of the graph, this gets destroyed in iterates */
448 PetscCall(MatDuplicate(a_Gmat, MAT_COPY_VALUES, &cMat));
449 PetscCall(MatConvert(cMat, MATAIJ, MAT_INPLACE_MATRIX, &cMat));
450 isMPI = (PetscBool)(size > 1);
451 if (isMPI) {
452 /* list of deleted ghosts, should compress this */
453 PetscCall(PetscCDCreate(size, &ghost_deleted_list));
454 PetscCall(PetscCDSetChunkSize(ghost_deleted_list, 100));
455 }
456 for (PetscInt iter = 0; iter < n_iter; iter++) {
457 const PetscScalar *lghost_max_ew, *lid_max_ew;
458 PetscBool *lghost_matched;
459 PetscMPIInt *lghost_pe, *lghost_max_pe;
460 Vec locMaxEdge, ghostMaxEdge, ghostMaxPE, locMaxPE;
461 PetscInt *lghost_gid, nEdges, nEdges0, num_ghosts = 0;
462 Edge *Edges;
463 const PetscInt n_sub_its = 1000; // in case of a bug, stop at some point
464
465 /* get submatrices of cMat */
466 for (PetscInt kk = 0; kk < nloc; kk++) lid_cprowID[kk] = -1;
467 if (isMPI) {
468 mpimat = (Mat_MPIAIJ *)cMat->data;
469 matA = (Mat_SeqAIJ *)mpimat->A->data;
470 matB = (Mat_SeqAIJ *)mpimat->B->data;
471 if (!matB->compressedrow.use) {
472 /* force construction of compressed row data structure since code below requires it */
473 PetscCall(MatCheckCompressedRow(mpimat->B, matB->nonzerorowcnt, &matB->compressedrow, matB->i, mpimat->B->rmap->n, -1.0));
474 }
475 /* set index into compressed row 'lid_cprowID' */
476 for (ix = 0; ix < matB->compressedrow.nrows; ix++) {
477 PetscInt *ridx = matB->compressedrow.rindex, lid = ridx[ix];
478 if (ridx[ix] >= 0) lid_cprowID[lid] = ix;
479 }
480 } else {
481 matA = (Mat_SeqAIJ *)cMat->data;
482 }
483 /* set matched flags: true for empty list */
484 for (PetscInt kk = 0; kk < nloc; kk++) {
485 PetscCall(PetscCDCountAt(agg_llists, kk, &ix));
486 if (ix > 0) lid_matched[kk] = PETSC_FALSE;
487 else lid_matched[kk] = PETSC_TRUE; // call deleted gids as matched
488 }
489 /* max edge and pe vecs */
490 PetscCall(MatCreateVecs(cMat, &locMaxEdge, NULL));
491 PetscCall(MatCreateVecs(cMat, &locMaxPE, NULL));
492 /* get 'lghost_pe' & 'lghost_gid' & init. 'lghost_matched' using 'mpimat->lvec' */
493 if (isMPI) {
494 Vec vec;
495 PetscScalar vval;
496 const PetscScalar *buf;
497
498 PetscCall(MatCreateVecs(cMat, &vec, NULL));
499 PetscCall(VecGetLocalSize(mpimat->lvec, &num_ghosts));
500 /* lghost_matched */
501 for (PetscInt kk = 0, gid = Istart; kk < nloc; kk++, gid++) {
502 PetscScalar vval = lid_matched[kk] ? 1.0 : 0.0;
503
504 PetscCall(VecSetValues(vec, 1, &gid, &vval, INSERT_VALUES));
505 }
506 PetscCall(VecAssemblyBegin(vec));
507 PetscCall(VecAssemblyEnd(vec));
508 PetscCall(VecScatterBegin(mpimat->Mvctx, vec, mpimat->lvec, INSERT_VALUES, SCATTER_FORWARD));
509 PetscCall(VecScatterEnd(mpimat->Mvctx, vec, mpimat->lvec, INSERT_VALUES, SCATTER_FORWARD));
510 PetscCall(VecGetArrayRead(mpimat->lvec, &buf)); /* get proc ID in 'buf' */
511 PetscCall(PetscMalloc4(num_ghosts, &lghost_matched, num_ghosts, &lghost_pe, num_ghosts, &lghost_gid, num_ghosts, &lghost_max_pe));
512
513 for (PetscInt kk = 0; kk < num_ghosts; kk++) lghost_matched[kk] = (PetscBool)(PetscRealPart(buf[kk]) != 0); // the proc of the ghost for now
514 PetscCall(VecRestoreArrayRead(mpimat->lvec, &buf));
515 /* lghost_pe */
516 vval = (PetscScalar)rank;
517 for (PetscInt kk = 0, gid = Istart; kk < nloc; kk++, gid++) PetscCall(VecSetValues(vec, 1, &gid, &vval, INSERT_VALUES)); /* set with GID */
518 PetscCall(VecAssemblyBegin(vec));
519 PetscCall(VecAssemblyEnd(vec));
520 PetscCall(VecScatterBegin(mpimat->Mvctx, vec, mpimat->lvec, INSERT_VALUES, SCATTER_FORWARD));
521 PetscCall(VecScatterEnd(mpimat->Mvctx, vec, mpimat->lvec, INSERT_VALUES, SCATTER_FORWARD));
522 PetscCall(VecGetArrayRead(mpimat->lvec, &buf)); /* get proc ID in 'buf' */
523 for (PetscInt kk = 0; kk < num_ghosts; kk++) lghost_pe[kk] = (PetscMPIInt)PetscRealPart(buf[kk]); // the proc of the ghost for now
524 PetscCall(VecRestoreArrayRead(mpimat->lvec, &buf));
525 /* lghost_gid */
526 for (PetscInt kk = 0, gid = Istart; kk < nloc; kk++, gid++) {
527 vval = (PetscScalar)gid;
528
529 PetscCall(VecSetValues(vec, 1, &gid, &vval, INSERT_VALUES)); /* set with GID */
530 }
531 PetscCall(VecAssemblyBegin(vec));
532 PetscCall(VecAssemblyEnd(vec));
533 PetscCall(VecScatterBegin(mpimat->Mvctx, vec, mpimat->lvec, INSERT_VALUES, SCATTER_FORWARD));
534 PetscCall(VecScatterEnd(mpimat->Mvctx, vec, mpimat->lvec, INSERT_VALUES, SCATTER_FORWARD));
535 PetscCall(VecDestroy(&vec));
536 PetscCall(VecGetArrayRead(mpimat->lvec, &buf)); /* get proc ID in 'lghost_gid' */
537 for (PetscInt kk = 0; kk < num_ghosts; kk++) lghost_gid[kk] = (PetscInt)PetscRealPart(buf[kk]);
538 PetscCall(VecRestoreArrayRead(mpimat->lvec, &buf));
539 }
540 // get 'comm_procs' (could hoist)
541 for (PetscInt kk = 0; kk < REQ_BF_SIZE; kk++) comm_procs[kk] = -1;
542 for (ix = 0, ncomm_procs = 0; ix < num_ghosts; ix++) {
543 PetscMPIInt proc = lghost_pe[ix], idx = -1;
544
545 for (PetscMPIInt k = 0; k < ncomm_procs && idx == -1; k++)
546 if (comm_procs[k] == proc) idx = k;
547 if (idx == -1) comm_procs[ncomm_procs++] = proc;
548 PetscCheck(ncomm_procs != REQ_BF_SIZE, PETSC_COMM_SELF, PETSC_ERR_SUP, "Receive request array too small: %d", ncomm_procs);
549 }
550 /* count edges, compute initial 'locMaxEdge', 'locMaxPE' */
551 nEdges0 = 0;
552 for (PetscInt kk = 0, gid = Istart; kk < nloc; kk++, gid++) {
553 PetscReal max_e = 0., tt;
554 PetscScalar vval;
555 PetscInt lid = kk, max_pe = rank, pe, n;
556
557 ii = matA->i;
558 n = ii[lid + 1] - ii[lid];
559 aj = PetscSafePointerPlusOffset(matA->j, ii[lid]);
560 ap = PetscSafePointerPlusOffset(matA->a, ii[lid]);
561 for (PetscInt jj = 0; jj < n; jj++) {
562 PetscInt lidj = aj[jj];
563
564 if ((tt = PetscRealPart(ap[jj])) > threshold && lidj != lid) {
565 if (tt > max_e) max_e = tt;
566 if (lidj > lid) nEdges0++;
567 }
568 }
569 if ((ix = lid_cprowID[lid]) != -1) { /* if I have any ghost neighbors */
570 ii = matB->compressedrow.i;
571 n = ii[ix + 1] - ii[ix];
572 ap = matB->a + ii[ix];
573 aj = matB->j + ii[ix];
574 for (PetscInt jj = 0; jj < n; jj++) {
575 if ((tt = PetscRealPart(ap[jj])) > threshold) {
576 if (tt > max_e) max_e = tt;
577 nEdges0++;
578 if ((pe = lghost_pe[aj[jj]]) > max_pe) max_pe = pe;
579 }
580 }
581 }
582 vval = max_e;
583 PetscCall(VecSetValues(locMaxEdge, 1, &gid, &vval, INSERT_VALUES));
584 vval = (PetscScalar)max_pe;
585 PetscCall(VecSetValues(locMaxPE, 1, &gid, &vval, INSERT_VALUES));
586 if (iter == 0 && max_e <= MY_MEPS) { // add BCs to fake aggregate
587 lid_matched[lid] = PETSC_TRUE;
588 if (bc_agg == -1) {
589 bc_agg = lid;
590 PetscCall(PetscCDCreate(1, &bc_list));
591 }
592 PetscCall(PetscCDRemoveAllAt(agg_llists, lid));
593 PetscCall(PetscCDAppendID(bc_list, 0, Istart + lid));
594 }
595 }
596 PetscCall(VecAssemblyBegin(locMaxEdge));
597 PetscCall(VecAssemblyEnd(locMaxEdge));
598 PetscCall(VecAssemblyBegin(locMaxPE));
599 PetscCall(VecAssemblyEnd(locMaxPE));
600 /* make 'ghostMaxEdge_max_ew', 'lghost_max_pe' */
601 if (isMPI) {
602 const PetscScalar *buf;
603
604 PetscCall(VecDuplicate(mpimat->lvec, &ghostMaxEdge));
605 PetscCall(VecScatterBegin(mpimat->Mvctx, locMaxEdge, ghostMaxEdge, INSERT_VALUES, SCATTER_FORWARD));
606 PetscCall(VecScatterEnd(mpimat->Mvctx, locMaxEdge, ghostMaxEdge, INSERT_VALUES, SCATTER_FORWARD));
607
608 PetscCall(VecDuplicate(mpimat->lvec, &ghostMaxPE));
609 PetscCall(VecScatterBegin(mpimat->Mvctx, locMaxPE, ghostMaxPE, INSERT_VALUES, SCATTER_FORWARD));
610 PetscCall(VecScatterEnd(mpimat->Mvctx, locMaxPE, ghostMaxPE, INSERT_VALUES, SCATTER_FORWARD));
611 PetscCall(VecGetArrayRead(ghostMaxPE, &buf));
612 for (PetscInt kk = 0; kk < num_ghosts; kk++) lghost_max_pe[kk] = (PetscMPIInt)PetscRealPart(buf[kk]); // the MAX proc of the ghost now
613 PetscCall(VecRestoreArrayRead(ghostMaxPE, &buf));
614 }
615 { // make lid_max_pe
616 const PetscScalar *buf;
617
618 PetscCall(VecGetArrayRead(locMaxPE, &buf));
619 for (PetscInt kk = 0; kk < nloc; kk++) lid_max_pe[kk] = (PetscMPIInt)PetscRealPart(buf[kk]); // the MAX proc of the ghost now
620 PetscCall(VecRestoreArrayRead(locMaxPE, &buf));
621 }
622 /* setup sorted list of edges, and make 'Edges' */
623 PetscCall(PetscMalloc1(nEdges0, &Edges));
624 nEdges = 0;
625 for (PetscInt kk = 0, n; kk < nloc; kk++) {
626 const PetscInt lid = kk;
627 PetscReal tt;
628
629 ii = matA->i;
630 n = ii[lid + 1] - ii[lid];
631 aj = PetscSafePointerPlusOffset(matA->j, ii[lid]);
632 ap = PetscSafePointerPlusOffset(matA->a, ii[lid]);
633 for (PetscInt jj = 0; jj < n; jj++) {
634 PetscInt lidj = aj[jj];
635
636 if ((tt = PetscRealPart(ap[jj])) > threshold && lidj != lid) {
637 if (lidj > lid) {
638 Edges[nEdges].lid0 = lid;
639 Edges[nEdges].gid1 = lidj + Istart;
640 Edges[nEdges].ghost1_idx = -1;
641 Edges[nEdges].weight = tt;
642 nEdges++;
643 }
644 }
645 }
646 if ((ix = lid_cprowID[lid]) != -1) { /* if I have any ghost neighbor */
647 ii = matB->compressedrow.i;
648 n = ii[ix + 1] - ii[ix];
649 ap = matB->a + ii[ix];
650 aj = matB->j + ii[ix];
651 for (PetscInt jj = 0; jj < n; jj++) {
652 if ((tt = PetscRealPart(ap[jj])) > threshold) {
653 Edges[nEdges].lid0 = lid;
654 Edges[nEdges].gid1 = lghost_gid[aj[jj]];
655 Edges[nEdges].ghost1_idx = aj[jj];
656 Edges[nEdges].weight = tt;
657 nEdges++;
658 }
659 }
660 }
661 }
662 PetscCheck(nEdges == nEdges0, PETSC_COMM_SELF, PETSC_ERR_SUP, "nEdges != nEdges0: %" PetscInt_FMT " %" PetscInt_FMT, nEdges0, nEdges);
663 if (Edges) qsort(Edges, nEdges, sizeof(Edge), gamg_hem_compare);
664
665 PetscCall(PetscInfo(info_is, "[%d] HEM iteration %" PetscInt_FMT " with %" PetscInt_FMT " edges\n", rank, iter, nEdges));
666
667 /* projection matrix */
668 PetscCall(MatCreate(comm, &P));
669 PetscCall(MatSetType(P, MATAIJ));
670 PetscCall(MatSetSizes(P, nloc, nloc, PETSC_DETERMINE, PETSC_DETERMINE));
671 PetscCall(MatMPIAIJSetPreallocation(P, 1, NULL, 1, NULL));
672 PetscCall(MatSeqAIJSetPreallocation(P, 1, NULL));
673 PetscCall(MatSetUp(P));
674 /* process - communicate - process */
675 for (PetscInt sub_it = 0, old_num_edge = 0; /* sub_it < n_sub_its */; /* sub_it++ */) {
676 PetscInt nactive_edges = 0, n_act_n[3], gn_act_n[3];
677 PetscMPIInt tag1, tag2;
678
679 PetscCall(VecGetArrayRead(locMaxEdge, &lid_max_ew));
680 if (isMPI) {
681 PetscCall(VecGetArrayRead(ghostMaxEdge, &lghost_max_ew));
682 PetscCall(PetscCommGetNewTag(comm, &tag1));
683 PetscCall(PetscCommGetNewTag(comm, &tag2));
684 }
685 for (PetscInt kk = 0; kk < nEdges; kk++) {
686 const Edge *e = &Edges[kk];
687 const PetscInt lid0 = e->lid0, gid1 = e->gid1, ghost1_idx = e->ghost1_idx, gid0 = lid0 + Istart, lid1 = gid1 - Istart;
688 PetscBool isOK = PETSC_TRUE, print = PETSC_FALSE;
689
690 if (print)
691 PetscCall(PetscSynchronizedPrintf(comm, "\t[%d] edge (%" PetscInt_FMT " %" PetscInt_FMT "), %s %s %s\n", rank, gid0, gid1, lid_matched[lid0] ? "true" : "false", (ghost1_idx != -1 && lghost_matched[ghost1_idx]) ? "true" : "false", (ghost1_idx == -1 && lid_matched[lid1]) ? "true" : "false"));
692 /* skip if either vertex is matched already */
693 if (lid_matched[lid0] || (ghost1_idx != -1 && lghost_matched[ghost1_idx]) || (ghost1_idx == -1 && lid_matched[lid1])) continue;
694
695 nactive_edges++;
696 PetscCheck(PetscRealPart(lid_max_ew[lid0]) >= e->weight - MY_MEPS, PETSC_COMM_SELF, PETSC_ERR_SUP, "edge weight %e > max %e", (double)e->weight, (double)PetscRealPart(lid_max_ew[lid0]));
697 if (print) PetscCall(PetscSynchronizedPrintf(comm, "\t[%d] active edge (%" PetscInt_FMT " %" PetscInt_FMT "), diff0 = %10.4e\n", rank, gid0, gid1, (double)(PetscRealPart(lid_max_ew[lid0]) - e->weight)));
698 // smaller edge, lid_max_ew get updated - e0
699 if (PetscRealPart(lid_max_ew[lid0]) > e->weight + MY_MEPS) {
700 if (print)
701 PetscCall(PetscSynchronizedPrintf(comm, "\t\t[%d] 1) e0 SKIPPING small edge %20.14e edge (%" PetscInt_FMT " %" PetscInt_FMT "), diff = %10.4e to proc %d. max = %20.14e, w = %20.14e\n", rank, (double)e->weight, gid0, gid1, (double)(PetscRealPart(lid_max_ew[lid0]) - e->weight), ghost1_idx != -1 ? lghost_pe[ghost1_idx] : rank, (double)PetscRealPart(lid_max_ew[lid0]),
702 (double)e->weight));
703 continue; // we are basically filter edges here
704 }
705 // e1 - local
706 if (ghost1_idx == -1) {
707 if (PetscRealPart(lid_max_ew[lid1]) > e->weight + MY_MEPS) {
708 if (print)
709 PetscCall(PetscSynchronizedPrintf(comm, "\t\t%c[%d] 2) e1 SKIPPING small local edge %20.14e edge (%" PetscInt_FMT " %" PetscInt_FMT "), diff = %10.4e\n", ghost1_idx != -1 ? '\t' : ' ', rank, (double)e->weight, gid0, gid1, (double)(PetscRealPart(lid_max_ew[lid1]) - e->weight)));
710 continue; // we are basically filter edges here
711 }
712 } else { // e1 - ghost
713 /* see if edge might get matched on other proc */
714 PetscReal g_max_e1 = PetscRealPart(lghost_max_ew[ghost1_idx]);
715
716 if (print)
717 PetscCall(PetscSynchronizedPrintf(comm, "\t\t\t[%d] CHECK GHOST e1, edge (%" PetscInt_FMT " %" PetscInt_FMT "), E0 MAX EDGE WEIGHT = %10.4e, EDGE WEIGHT = %10.4e, diff1 = %10.4e, ghost proc %d with max pe %d on e0 and %d on e1\n", rank, gid0, gid1, (double)PetscRealPart(lid_max_ew[lid0]),
718 (double)e->weight, (double)(PetscRealPart(lghost_max_ew[ghost1_idx]) - e->weight), lghost_pe[ghost1_idx], lid_max_pe[lid0], lghost_max_pe[ghost1_idx]));
719 if (g_max_e1 > e->weight + MY_MEPS) {
720 /* PetscCall(PetscSynchronizedPrintf(comm,"\t\t\t\t[%d] 3) ghost e1 SKIPPING small edge (%d %d), diff = %10.4e from proc %d with max pe %d. max = %20.14e, w = %20.14e\n", rank, gid0, gid1, g_max_e1 - e->weight, lghost_pe[ghost1_idx], lghost_max_pe[ghost1_idx], g_max_e1, e->weight )); */
721 continue;
722 } else if (g_max_e1 >= e->weight - MY_MEPS && lghost_pe[ghost1_idx] > rank) { // is 'lghost_max_pe[ghost1_idx] > rank' needed?
723 /* check for max_ea == to this edge and larger processor that will deal with this */
724 if (print)
725 PetscCall(PetscSynchronizedPrintf(comm, "\t\t\t[%d] ghost e1 SKIPPING EQUAL (%" PetscInt_FMT " %" PetscInt_FMT "), diff = %10.4e from larger proc %d with max pe %d. max = %20.14e, w = %20.14e\n", rank, gid0, gid1, (double)(PetscRealPart(lid_max_ew[lid0]) - e->weight), lghost_pe[ghost1_idx], lghost_max_pe[ghost1_idx], (double)g_max_e1,
726 (double)e->weight));
727 continue;
728 } else {
729 /* PetscCall(PetscSynchronizedPrintf(comm,"\t[%d] Edge (%d %d) passes gid0 tests, diff = %10.4e from proc %d with max pe %d. max = %20.14e, w = %20.14e\n", rank, gid0, gid1, g_max_e1 - e->weight, lghost_pe[ghost1_idx], lghost_max_pe[ghost1_idx], g_max_e1, e->weight )); */
730 }
731 }
732 /* check ghost for v0 */
733 if (isOK) {
734 PetscReal max_e, ew;
735
736 if ((ix = lid_cprowID[lid0]) != -1) { /* if I have any ghost neighbors */
737 PetscInt n;
738
739 ii = matB->compressedrow.i;
740 n = ii[ix + 1] - ii[ix];
741 ap = matB->a + ii[ix];
742 aj = matB->j + ii[ix];
743 for (PetscInt jj = 0; jj < n && isOK; jj++) {
744 PetscInt lidj = aj[jj];
745
746 if (lghost_matched[lidj]) continue;
747 ew = PetscRealPart(ap[jj]);
748 if (ew <= threshold) continue;
749 max_e = PetscRealPart(lghost_max_ew[lidj]);
750
751 /* check for max_e == to this edge and larger processor that will deal with this */
752 if (ew >= PetscRealPart(lid_max_ew[lid0]) - MY_MEPS && lghost_max_pe[lidj] > rank) isOK = PETSC_FALSE;
753 PetscCheck(ew <= max_e + MY_MEPS, PETSC_COMM_SELF, PETSC_ERR_SUP, "edge weight %e > max %e. ncols = %" PetscInt_FMT ", gid0 = %" PetscInt_FMT ", gid1 = %" PetscInt_FMT, (double)PetscRealPart(ew), (double)PetscRealPart(max_e), n, lid0 + Istart, lghost_gid[lidj]);
754 if (print)
755 PetscCall(PetscSynchronizedPrintf(comm, "\t\t\t\t[%d] e0: looked at ghost adj (%" PetscInt_FMT " %" PetscInt_FMT "), diff = %10.4e, ghost on proc %d (max %d). isOK = %d, %d %d %d; ew = %e, lid0 max ew = %e, diff = %e, eps = %e\n", rank, gid0, lghost_gid[lidj], (double)(max_e - ew), lghost_pe[lidj], lghost_max_pe[lidj], isOK, (double)ew >= (double)(max_e - MY_MEPS), ew >= PetscRealPart(lid_max_ew[lid0]) - MY_MEPS, lghost_pe[lidj] > rank, (double)ew, (double)PetscRealPart(lid_max_ew[lid0]), (double)(ew - PetscRealPart(lid_max_ew[lid0])), (double)MY_MEPS));
756 }
757 if (!isOK && print) PetscCall(PetscSynchronizedPrintf(comm, "\t\t[%d] skip edge (%" PetscInt_FMT " %" PetscInt_FMT ") from ghost inspection\n", rank, gid0, gid1));
758 }
759 /* check local v1 */
760 if (ghost1_idx == -1) {
761 if ((ix = lid_cprowID[lid1]) != -1) { /* if I have any ghost neighbors */
762 PetscInt n;
763
764 ii = matB->compressedrow.i;
765 n = ii[ix + 1] - ii[ix];
766 ap = matB->a + ii[ix];
767 aj = matB->j + ii[ix];
768 for (PetscInt jj = 0; jj < n && isOK; jj++) {
769 PetscInt lidj = aj[jj];
770
771 if (lghost_matched[lidj]) continue;
772 ew = PetscRealPart(ap[jj]);
773 if (ew <= threshold) continue;
774 max_e = PetscRealPart(lghost_max_ew[lidj]);
775 /* check for max_e == to this edge and larger processor that will deal with this */
776 if (ew >= PetscRealPart(lid_max_ew[lid1]) - MY_MEPS && lghost_max_pe[lidj] > rank) isOK = PETSC_FALSE;
777 PetscCheck(ew <= max_e + MY_MEPS, PETSC_COMM_SELF, PETSC_ERR_SUP, "edge weight %e > max %e", (double)PetscRealPart(ew), (double)PetscRealPart(max_e));
778 if (print)
779 PetscCall(PetscSynchronizedPrintf(comm, "\t\t\t\t\t[%d] e1: looked at ghost adj (%" PetscInt_FMT " %" PetscInt_FMT "), diff = %10.4e, ghost on proc %d (max %d)\n", rank, gid0, lghost_gid[lidj], (double)(max_e - ew), lghost_pe[lidj], lghost_max_pe[lidj]));
780 }
781 }
782 if (!isOK && print) PetscCall(PetscSynchronizedPrintf(comm, "\t\t[%d] skip edge (%" PetscInt_FMT " %" PetscInt_FMT ") from ghost inspection\n", rank, gid0, gid1));
783 }
784 }
785 PetscReal e1_max_w = (ghost1_idx == -1 ? PetscRealPart(lid_max_ew[lid0]) : PetscRealPart(lghost_max_ew[ghost1_idx]));
786 if (print)
787 PetscCall(PetscSynchronizedPrintf(comm, "\t[%d] MATCHING (%" PetscInt_FMT " %" PetscInt_FMT ") e1 max weight = %e, e1 weight diff %e, %s. isOK = %d\n", rank, gid0, gid1, (double)e1_max_w, (double)(e1_max_w - e->weight), ghost1_idx == -1 ? "local" : "ghost", isOK));
788 /* do it */
789 if (isOK) {
790 if (ghost1_idx == -1) {
791 PetscCheck(!lid_matched[lid1], PETSC_COMM_SELF, PETSC_ERR_SUP, "local %" PetscInt_FMT " is matched", gid1);
792 lid_matched[lid1] = PETSC_TRUE; /* keep track of what we've done this round */
793 PetscCall(PetscCDMoveAppend(agg_llists, lid0, lid1)); // takes lid1's list and appends to lid0's
794 } else {
795 /* add gid1 to list of ghost deleted by me -- I need their children */
796 PetscMPIInt proc = lghost_pe[ghost1_idx];
797 PetscCheck(!lghost_matched[ghost1_idx], PETSC_COMM_SELF, PETSC_ERR_SUP, "ghost %" PetscInt_FMT " is matched", lghost_gid[ghost1_idx]);
798 lghost_matched[ghost1_idx] = PETSC_TRUE;
799 PetscCall(PetscCDAppendID(ghost_deleted_list, proc, ghost1_idx)); /* cache to send messages */
800 PetscCall(PetscCDAppendID(ghost_deleted_list, proc, lid0));
801 }
802 lid_matched[lid0] = PETSC_TRUE; /* keep track of what we've done this round */
803 /* set projection */
804 PetscCall(MatSetValues(P, 1, &gid0, 1, &gid0, &one, INSERT_VALUES));
805 PetscCall(MatSetValues(P, 1, &gid1, 1, &gid0, &one, INSERT_VALUES));
806 //PetscCall(PetscPrintf(comm,"\t %" PetscInt_FMT ".%" PetscInt_FMT ") match active EDGE %" PetscInt_FMT " : (%" PetscInt_FMT " %" PetscInt_FMT ")\n",iter,sub_it, nactive_edges, gid0, gid1));
807 } /* matched */
808 } /* edge loop */
809 PetscCall(PetscSynchronizedFlush(comm, PETSC_STDOUT));
810 if (isMPI) PetscCall(VecRestoreArrayRead(ghostMaxEdge, &lghost_max_ew));
811 PetscCall(VecRestoreArrayRead(locMaxEdge, &lid_max_ew));
812 // count active for test, latter, update deleted ghosts
813 n_act_n[0] = nactive_edges;
814 if (ghost_deleted_list) PetscCall(PetscCDCount(ghost_deleted_list, &n_act_n[2]));
815 else n_act_n[2] = 0;
816 PetscCall(PetscCDCount(agg_llists, &n_act_n[1]));
817 PetscCallMPI(MPIU_Allreduce(n_act_n, gn_act_n, 3, MPIU_INT, MPI_SUM, comm));
818 PetscCall(PetscInfo(info_is, "[%d] %" PetscInt_FMT ".%" PetscInt_FMT ") nactive edges=%" PetscInt_FMT ", ncomm_procs=%d, nEdges=%" PetscInt_FMT ", %" PetscInt_FMT " deleted ghosts, N=%" PetscInt_FMT "\n", rank, iter, sub_it, gn_act_n[0], ncomm_procs, nEdges, gn_act_n[2], gn_act_n[1]));
819 /* deal with deleted ghost */
820 if (isMPI) {
821 PetscCDIntNd *pos;
822 PetscInt *sbuffs1[REQ_BF_SIZE], ndel;
823 PetscInt *sbuffs2[REQ_BF_SIZE];
824 MPI_Status status;
825
826 /* send deleted ghosts */
827 for (PetscInt proc_idx = 0; proc_idx < ncomm_procs; proc_idx++) {
828 const PetscMPIInt proc = comm_procs[proc_idx];
829 PetscInt *sbuff, *pt, scount;
830 MPI_Request *request;
831
832 /* count ghosts */
833 PetscCall(PetscCDCountAt(ghost_deleted_list, proc, &ndel));
834 ndel /= 2; // two entries for each proc
835 scount = 2 + 2 * ndel;
836 PetscCall(PetscMalloc1(scount + request_size, &sbuff));
837 /* save requests */
838 sbuffs1[proc_idx] = sbuff;
839 request = (MPI_Request *)sbuff;
840 sbuff = pt = sbuff + request_size;
841 /* write [ndel, proc, n*[gid1,gid0] */
842 *pt++ = ndel; // number of deleted to send
843 *pt++ = rank; // proc (not used)
844 PetscCall(PetscCDGetHeadPos(ghost_deleted_list, proc, &pos));
845 while (pos) {
846 PetscInt lid0, ghost_idx, gid1;
847
848 PetscCall(PetscCDIntNdGetID(pos, &ghost_idx));
849 gid1 = lghost_gid[ghost_idx];
850 PetscCall(PetscCDGetNextPos(ghost_deleted_list, proc, &pos));
851 PetscCall(PetscCDIntNdGetID(pos, &lid0));
852 PetscCall(PetscCDGetNextPos(ghost_deleted_list, proc, &pos));
853 *pt++ = gid1;
854 *pt++ = lid0 + Istart; // gid0
855 }
856 PetscCheck(pt - sbuff == scount, PETSC_COMM_SELF, PETSC_ERR_SUP, "sbuff-pt != scount: %zu", pt - sbuff);
857 /* MPIU_Isend: tag1 [ndel, proc, n*[gid1,gid0] ] */
858 PetscCallMPI(MPIU_Isend(sbuff, scount, MPIU_INT, proc, tag1, comm, request));
859 PetscCall(PetscCDRemoveAllAt(ghost_deleted_list, proc)); // done with this list
860 }
861 /* receive deleted, send back partial aggregates, clear lists */
862 for (PetscInt proc_idx = 0; proc_idx < ncomm_procs; proc_idx++) {
863 PetscCallMPI(MPI_Probe(comm_procs[proc_idx] /* MPI_ANY_SOURCE */, tag1, comm, &status));
864 {
865 PetscInt *pt, *pt2, *pt3, *sbuff, tmp;
866 MPI_Request *request;
867 PetscMPIInt rcount, scount;
868 const PetscMPIInt proc = status.MPI_SOURCE;
869
870 PetscCallMPI(MPI_Get_count(&status, MPIU_INT, &rcount));
871 if (rcount > rbuff_sz) {
872 if (rbuff) PetscCall(PetscFree(rbuff));
873 PetscCall(PetscMalloc1(rcount, &rbuff));
874 rbuff_sz = rcount;
875 }
876 /* MPI_Recv: tag1 [ndel, proc, ndel*[gid1,gid0] ] */
877 PetscCallMPI(MPI_Recv(rbuff, rcount, MPIU_INT, proc, tag1, comm, &status));
878 /* read and count sends *[lid0, n, n*[gid] ] */
879 pt = rbuff;
880 scount = 0;
881 ndel = *pt++; // number of deleted to recv
882 tmp = *pt++; // proc (not used)
883 while (ndel--) {
884 PetscInt gid1 = *pt++, lid1 = gid1 - Istart;
885 PetscInt gh_gid0 = *pt++; // gid on other proc (not used here to count)
886
887 PetscCheck(lid1 >= 0 && lid1 < nloc, PETSC_COMM_SELF, PETSC_ERR_SUP, "received ghost deleted %" PetscInt_FMT, gid1);
888 PetscCheck(!lid_matched[lid1], PETSC_COMM_SELF, PETSC_ERR_PLIB, "%" PetscInt_FMT ") received matched local gid %" PetscInt_FMT ",%" PetscInt_FMT ", with ghost (lid) %" PetscInt_FMT " from proc %d", sub_it, gid1, gh_gid0, tmp, proc);
889 lid_matched[lid1] = PETSC_TRUE; /* keep track of what we've done this round */
890 PetscCall(PetscCDCountAt(agg_llists, lid1, &tmp)); // n
891 scount += tmp + 2; // lid0, n, n*[gid]
892 }
893 PetscCheck((pt - rbuff) == (ptrdiff_t)rcount, PETSC_COMM_SELF, PETSC_ERR_SUP, "receive buffer size != num read: %zu; rcount: %d", pt - rbuff, rcount);
894 /* send tag2: *[gid0, n, n*[gid] ] */
895 PetscCall(PetscMalloc1(scount + request_size, &sbuff));
896 sbuffs2[proc_idx] = sbuff; /* cache request */
897 request = (MPI_Request *)sbuff;
898 pt2 = sbuff = sbuff + request_size;
899 // read again: n, proc, n*[gid1,gid0]
900 pt = rbuff;
901 ndel = *pt++;
902 tmp = *pt++; // proc (not used)
903 while (ndel--) {
904 PetscInt gid1 = *pt++, lid1 = gid1 - Istart, gh_gid0 = *pt++;
905
906 /* write [gid0, aggSz, aggSz[gid] ] */
907 *pt2++ = gh_gid0;
908 pt3 = pt2++; /* save pointer for later */
909 PetscCall(PetscCDGetHeadPos(agg_llists, lid1, &pos));
910 while (pos) {
911 PetscInt gid;
912
913 PetscCall(PetscCDIntNdGetID(pos, &gid));
914 PetscCall(PetscCDGetNextPos(agg_llists, lid1, &pos));
915 *pt2++ = gid;
916 }
917 PetscCall(PetscIntCast(pt2 - pt3 - 1, pt3));
918 /* clear list */
919 PetscCall(PetscCDRemoveAllAt(agg_llists, lid1));
920 }
921 PetscCheck((pt2 - sbuff) == (ptrdiff_t)scount, PETSC_COMM_SELF, PETSC_ERR_SUP, "buffer size != num write: %zu %d", pt2 - sbuff, scount);
922 /* MPIU_Isend: requested data tag2 *[lid0, n, n*[gid1] ] */
923 PetscCallMPI(MPIU_Isend(sbuff, scount, MPIU_INT, proc, tag2, comm, request));
924 }
925 } // proc_idx
926 /* receive tag2 *[gid0, n, n*[gid] ] */
927 for (PetscMPIInt proc_idx = 0; proc_idx < ncomm_procs; proc_idx++) {
928 PetscMPIInt proc;
929 PetscInt *pt;
930 int rcount;
931
932 PetscCallMPI(MPI_Probe(comm_procs[proc_idx] /* MPI_ANY_SOURCE */, tag2, comm, &status));
933 PetscCallMPI(MPI_Get_count(&status, MPIU_INT, &rcount));
934 if (rcount > rbuff_sz) {
935 if (rbuff) PetscCall(PetscFree(rbuff));
936 PetscCall(PetscMalloc1(rcount, &rbuff));
937 rbuff_sz = rcount;
938 }
939 proc = status.MPI_SOURCE;
940 /* MPI_Recv: tag1 [n, proc, n*[gid1,lid0] ] */
941 PetscCallMPI(MPI_Recv(rbuff, rcount, MPIU_INT, proc, tag2, comm, &status));
942 pt = rbuff;
943 while (pt - rbuff < rcount) {
944 PetscInt gid0 = *pt++, n = *pt++;
945
946 while (n--) {
947 PetscInt gid1 = *pt++;
948
949 PetscCall(PetscCDAppendID(agg_llists, gid0 - Istart, gid1));
950 }
951 }
952 PetscCheck((pt - rbuff) == (ptrdiff_t)rcount, PETSC_COMM_SELF, PETSC_ERR_SUP, "recv buffer size != num read: %zu %d", pt - rbuff, rcount);
953 }
954 /* wait for tag1 isends */
955 for (PetscMPIInt proc_idx = 0; proc_idx < ncomm_procs; proc_idx++) {
956 MPI_Request *request = (MPI_Request *)sbuffs1[proc_idx];
957
958 PetscCallMPI(MPI_Wait(request, &status));
959 PetscCall(PetscFree(sbuffs1[proc_idx]));
960 }
961 /* wait for tag2 isends */
962 for (PetscMPIInt proc_idx = 0; proc_idx < ncomm_procs; proc_idx++) {
963 MPI_Request *request = (MPI_Request *)sbuffs2[proc_idx];
964
965 PetscCallMPI(MPI_Wait(request, &status));
966 PetscCall(PetscFree(sbuffs2[proc_idx]));
967 }
968 } /* MPI */
969 /* set 'lghost_matched' - use locMaxEdge, ghostMaxEdge (recomputed next) */
970 if (isMPI) {
971 const PetscScalar *sbuff;
972
973 for (PetscInt kk = 0, gid = Istart; kk < nloc; kk++, gid++) {
974 PetscScalar vval = lid_matched[kk] ? 1.0 : 0.0;
975
976 PetscCall(VecSetValues(locMaxEdge, 1, &gid, &vval, INSERT_VALUES)); /* set with GID */
977 }
978 PetscCall(VecAssemblyBegin(locMaxEdge));
979 PetscCall(VecAssemblyEnd(locMaxEdge));
980 PetscCall(VecScatterBegin(mpimat->Mvctx, locMaxEdge, ghostMaxEdge, INSERT_VALUES, SCATTER_FORWARD));
981 PetscCall(VecScatterEnd(mpimat->Mvctx, locMaxEdge, ghostMaxEdge, INSERT_VALUES, SCATTER_FORWARD));
982 PetscCall(VecGetArrayRead(ghostMaxEdge, &sbuff));
983 for (PetscInt kk = 0; kk < num_ghosts; kk++) lghost_matched[kk] = (PetscBool)(PetscRealPart(sbuff[kk]) != 0.0);
984 PetscCall(VecRestoreArrayRead(ghostMaxEdge, &sbuff));
985 }
986 /* compute 'locMaxEdge' inside sub iteration b/c max weight can drop as neighbors are matched */
987 for (PetscInt kk = 0, gid = Istart; kk < nloc; kk++, gid++) {
988 PetscReal max_e = 0., tt;
989 PetscScalar vval;
990 const PetscInt lid = kk;
991 PetscMPIInt max_pe = rank, pe, n;
992
993 ii = matA->i;
994 PetscCall(PetscMPIIntCast(ii[lid + 1] - ii[lid], &n));
995 aj = PetscSafePointerPlusOffset(matA->j, ii[lid]);
996 ap = PetscSafePointerPlusOffset(matA->a, ii[lid]);
997 for (PetscMPIInt jj = 0; jj < n; jj++) {
998 PetscInt lidj = aj[jj];
999
1000 if (lid_matched[lidj]) continue; /* this is new - can change local max */
1001 if (lidj != lid && PetscRealPart(ap[jj]) > max_e) max_e = PetscRealPart(ap[jj]);
1002 }
1003 if (lid_cprowID && (ix = lid_cprowID[lid]) != -1) { /* if I have any ghost neighbors */
1004 ii = matB->compressedrow.i;
1005 PetscCall(PetscMPIIntCast(ii[ix + 1] - ii[ix], &n));
1006 ap = matB->a + ii[ix];
1007 aj = matB->j + ii[ix];
1008 for (PetscMPIInt jj = 0; jj < n; jj++) {
1009 PetscInt lidj = aj[jj];
1010
1011 if (lghost_matched[lidj]) continue;
1012 if ((tt = PetscRealPart(ap[jj])) > max_e) max_e = tt;
1013 }
1014 }
1015 vval = max_e;
1016 PetscCall(VecSetValues(locMaxEdge, 1, &gid, &vval, INSERT_VALUES)); /* set with GID */
1017 // max PE with max edge
1018 if (lid_cprowID && (ix = lid_cprowID[lid]) != -1) { /* if I have any ghost neighbors */
1019 ii = matB->compressedrow.i;
1020 PetscCall(PetscMPIIntCast(ii[ix + 1] - ii[ix], &n));
1021 ap = matB->a + ii[ix];
1022 aj = matB->j + ii[ix];
1023 for (PetscInt jj = 0; jj < n; jj++) {
1024 PetscInt lidj = aj[jj];
1025
1026 if (lghost_matched[lidj]) continue;
1027 if ((pe = lghost_pe[aj[jj]]) > max_pe && PetscRealPart(ap[jj]) >= max_e - MY_MEPS) max_pe = pe;
1028 }
1029 }
1030 vval = max_pe;
1031 PetscCall(VecSetValues(locMaxPE, 1, &gid, &vval, INSERT_VALUES));
1032 }
1033 PetscCall(VecAssemblyBegin(locMaxEdge));
1034 PetscCall(VecAssemblyEnd(locMaxEdge));
1035 PetscCall(VecAssemblyBegin(locMaxPE));
1036 PetscCall(VecAssemblyEnd(locMaxPE));
1037 /* compute 'lghost_max_ew' and 'lghost_max_pe' to get ready for next iteration*/
1038 if (isMPI) {
1039 const PetscScalar *buf;
1040
1041 PetscCall(VecScatterBegin(mpimat->Mvctx, locMaxEdge, ghostMaxEdge, INSERT_VALUES, SCATTER_FORWARD));
1042 PetscCall(VecScatterEnd(mpimat->Mvctx, locMaxEdge, ghostMaxEdge, INSERT_VALUES, SCATTER_FORWARD));
1043 PetscCall(VecScatterBegin(mpimat->Mvctx, locMaxPE, ghostMaxPE, INSERT_VALUES, SCATTER_FORWARD));
1044 PetscCall(VecScatterEnd(mpimat->Mvctx, locMaxPE, ghostMaxPE, INSERT_VALUES, SCATTER_FORWARD));
1045 PetscCall(VecGetArrayRead(ghostMaxPE, &buf));
1046 for (PetscInt kk = 0; kk < num_ghosts; kk++) lghost_max_pe[kk] = (PetscMPIInt)PetscRealPart(buf[kk]); // the MAX proc of the ghost now
1047 PetscCall(VecRestoreArrayRead(ghostMaxPE, &buf));
1048 }
1049 // if no active edges, stop
1050 if (gn_act_n[0] < 1) break;
1051 // inc and check (self stopping iteration
1052 PetscCheck(old_num_edge != gn_act_n[0], PETSC_COMM_SELF, PETSC_ERR_SUP, "HEM stalled step %" PetscInt_FMT "/%" PetscInt_FMT, sub_it + 1, n_sub_its);
1053 sub_it++;
1054 PetscCheck(sub_it < n_sub_its, PETSC_COMM_SELF, PETSC_ERR_SUP, "failed to finish HEM step %" PetscInt_FMT "/%" PetscInt_FMT, sub_it + 1, n_sub_its);
1055 old_num_edge = gn_act_n[0];
1056 } /* sub_it loop */
1057 /* clean up iteration */
1058 PetscCall(PetscFree(Edges));
1059 if (isMPI) { // can be hoisted
1060 PetscCall(VecRestoreArrayRead(ghostMaxEdge, &lghost_max_ew));
1061 PetscCall(VecDestroy(&ghostMaxEdge));
1062 PetscCall(VecDestroy(&ghostMaxPE));
1063 PetscCall(PetscFree4(lghost_matched, lghost_pe, lghost_gid, lghost_max_pe));
1064 }
1065 PetscCall(VecDestroy(&locMaxEdge));
1066 PetscCall(VecDestroy(&locMaxPE));
1067 /* create next graph */
1068 {
1069 Vec diag;
1070
1071 /* add identity for unmatched vertices so they stay alive */
1072 for (PetscInt kk = 0, gid1, gid = Istart; kk < nloc; kk++, gid++) {
1073 if (!lid_matched[kk]) {
1074 const PetscInt lid = kk;
1075 PetscCDIntNd *pos;
1076
1077 PetscCall(PetscCDGetHeadPos(agg_llists, lid, &pos));
1078 PetscCheck(pos, PETSC_COMM_SELF, PETSC_ERR_PLIB, "empty list in singleton: %" PetscInt_FMT, gid);
1079 PetscCall(PetscCDIntNdGetID(pos, &gid1));
1080 PetscCheck(gid1 == gid, PETSC_COMM_SELF, PETSC_ERR_PLIB, "first in list (%" PetscInt_FMT ") in singleton not %" PetscInt_FMT, gid1, gid);
1081 PetscCall(MatSetValues(P, 1, &gid, 1, &gid, &one, INSERT_VALUES));
1082 }
1083 }
1084 PetscCall(MatAssemblyBegin(P, MAT_FINAL_ASSEMBLY));
1085 PetscCall(MatAssemblyEnd(P, MAT_FINAL_ASSEMBLY));
1086
1087 /* project to make new graph with collapsed edges */
1088 PetscCall(MatPtAP(cMat, P, MAT_INITIAL_MATRIX, 1.0, &tMat));
1089 PetscCall(MatDestroy(&P));
1090 PetscCall(MatDestroy(&cMat));
1091 cMat = tMat;
1092 PetscCall(MatCreateVecs(cMat, &diag, NULL));
1093 PetscCall(MatGetDiagonal(cMat, diag));
1094 PetscCall(VecReciprocal(diag));
1095 PetscCall(VecSqrtAbs(diag));
1096 PetscCall(MatDiagonalScale(cMat, diag, diag));
1097 PetscCall(VecDestroy(&diag));
1098 }
1099 } /* coarsen iterator */
1100
1101 /* make fake matrix with Mat->B only for smoothed agg QR. Need this if we make an aux graph (ie, PtAP) with k > 1 */
1102 if (size > 1) {
1103 Mat mat;
1104 PetscCDIntNd *pos;
1105 PetscInt NN, MM, jj = 0, mxsz = 0;
1106
1107 for (PetscInt kk = 0; kk < nloc; kk++) {
1108 PetscCall(PetscCDCountAt(agg_llists, kk, &jj));
1109 if (jj > mxsz) mxsz = jj;
1110 }
1111 PetscCall(MatGetSize(a_Gmat, &MM, &NN));
1112 if (mxsz > MM - nloc) mxsz = MM - nloc;
1113 /* matrix of ghost adj for square graph */
1114 PetscCall(MatCreateAIJ(comm, nloc, nloc, PETSC_DETERMINE, PETSC_DETERMINE, 0, NULL, mxsz, NULL, &mat));
1115 for (PetscInt lid = 0, gid = Istart; lid < nloc; lid++, gid++) {
1116 PetscCall(PetscCDGetHeadPos(agg_llists, lid, &pos));
1117 while (pos) {
1118 PetscInt gid1;
1119
1120 PetscCall(PetscCDIntNdGetID(pos, &gid1));
1121 PetscCall(PetscCDGetNextPos(agg_llists, lid, &pos));
1122 if (gid1 < Istart || gid1 >= Istart + nloc) PetscCall(MatSetValues(mat, 1, &gid, 1, &gid1, &one, ADD_VALUES));
1123 }
1124 }
1125 PetscCall(MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY));
1126 PetscCall(MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY));
1127 PetscCall(PetscCDSetMat(agg_llists, mat));
1128 PetscCall(PetscCDDestroy(ghost_deleted_list));
1129 if (rbuff_sz) PetscCall(PetscFree(rbuff)); // always true
1130 }
1131 // move BCs into some node
1132 if (bc_list) {
1133 PetscCDIntNd *pos;
1134
1135 PetscCall(PetscCDGetHeadPos(bc_list, 0, &pos));
1136 while (pos) {
1137 PetscInt gid1;
1138
1139 PetscCall(PetscCDIntNdGetID(pos, &gid1));
1140 PetscCall(PetscCDGetNextPos(bc_list, 0, &pos));
1141 PetscCall(PetscCDAppendID(agg_llists, bc_agg, gid1));
1142 }
1143 PetscCall(PetscCDRemoveAllAt(bc_list, 0));
1144 PetscCall(PetscCDDestroy(bc_list));
1145 }
1146 {
1147 // check sizes -- all vertices must get in graph
1148 PetscInt sz, globalsz, MM;
1149
1150 PetscCall(MatGetSize(a_Gmat, &MM, NULL));
1151 PetscCall(PetscCDCount(agg_llists, &sz));
1152 PetscCallMPI(MPIU_Allreduce(&sz, &globalsz, 1, MPIU_INT, MPI_SUM, comm));
1153 PetscCheck(MM == globalsz, comm, PETSC_ERR_SUP, "lost %" PetscInt_FMT " equations ?", MM - globalsz);
1154 }
1155 // cleanup
1156 PetscCall(MatDestroy(&cMat));
1157 PetscCall(PetscFree3(lid_matched, lid_cprowID, lid_max_pe));
1158 PetscCall(ISDestroy(&info_is));
1159 PetscFunctionReturn(PETSC_SUCCESS);
1160 }
1161
1162 /*
1163 HEM coarsen, simple greedy.
1164 */
MatCoarsenApply_HEM(MatCoarsen coarse)1165 static PetscErrorCode MatCoarsenApply_HEM(MatCoarsen coarse)
1166 {
1167 Mat mat = coarse->graph;
1168
1169 PetscFunctionBegin;
1170 PetscCall(MatCoarsenApply_HEM_private(mat, coarse->max_it, coarse->threshold, &coarse->agg_lists));
1171 PetscFunctionReturn(PETSC_SUCCESS);
1172 }
1173
MatCoarsenView_HEM(MatCoarsen coarse,PetscViewer viewer)1174 static PetscErrorCode MatCoarsenView_HEM(MatCoarsen coarse, PetscViewer viewer)
1175 {
1176 PetscMPIInt rank;
1177 PetscBool isascii;
1178
1179 PetscFunctionBegin;
1180 PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)coarse), &rank));
1181 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
1182 if (isascii) {
1183 PetscCDIntNd *pos, *pos2;
1184 PetscViewerFormat format;
1185
1186 PetscCall(PetscViewerASCIIPrintf(viewer, "%" PetscInt_FMT " matching steps with threshold = %g\n", coarse->max_it, (double)coarse->threshold));
1187 PetscCall(PetscViewerGetFormat(viewer, &format));
1188 if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
1189 if (coarse->agg_lists) {
1190 PetscCall(PetscViewerASCIIPushSynchronized(viewer));
1191 for (PetscInt kk = 0; kk < coarse->agg_lists->size; kk++) {
1192 PetscCall(PetscCDGetHeadPos(coarse->agg_lists, kk, &pos));
1193 if ((pos2 = pos)) PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "selected local %" PetscInt_FMT ": ", kk));
1194 while (pos) {
1195 PetscInt gid1;
1196
1197 PetscCall(PetscCDIntNdGetID(pos, &gid1));
1198 PetscCall(PetscCDGetNextPos(coarse->agg_lists, kk, &pos));
1199 PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, " %" PetscInt_FMT " ", gid1));
1200 }
1201 if (pos2) PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "\n"));
1202 }
1203 PetscCall(PetscViewerFlush(viewer));
1204 PetscCall(PetscViewerASCIIPopSynchronized(viewer));
1205 } else {
1206 PetscCall(PetscViewerASCIIPrintf(viewer, " HEM aggregator lists are not available\n"));
1207 }
1208 }
1209 }
1210 PetscFunctionReturn(PETSC_SUCCESS);
1211 }
1212
1213 /*MC
1214 MATCOARSENHEM - A coarsener that uses HEM a simple greedy coarsener
1215
1216 Level: beginner
1217
1218 .seealso: `MatCoarsen`, `MatCoarsenMISKSetDistance()`, `MatCoarsenApply()`, `MatCoarsenSetType()`, `MatCoarsenType`, `MatCoarsenCreate()`, `MATCOARSENMISK`, `MATCOARSENMIS`
1219 M*/
1220
MatCoarsenCreate_HEM(MatCoarsen coarse)1221 PETSC_EXTERN PetscErrorCode MatCoarsenCreate_HEM(MatCoarsen coarse)
1222 {
1223 PetscFunctionBegin;
1224 coarse->ops->apply = MatCoarsenApply_HEM;
1225 coarse->ops->view = MatCoarsenView_HEM;
1226 coarse->max_it = 4;
1227 PetscFunctionReturn(PETSC_SUCCESS);
1228 }
1229