Actual source code: hem.c
1: #include <petsc/private/matimpl.h>
2: #include <../src/mat/impls/aij/seq/aij.h>
3: #include <../src/mat/impls/aij/mpi/mpiaij.h>
4: #include <petscdm.h>
6: /* linked list methods
7: *
8: * PetscCDCreate
9: */
10: PetscErrorCode PetscCDCreate(PetscInt a_size, PetscCoarsenData **a_out)
11: {
12: PetscCoarsenData *ail;
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: }
29: /* PetscCDDestroy
30: */
31: PetscErrorCode PetscCDDestroy(PetscCoarsenData *ail)
32: {
33: PetscCDArrNd *n = &ail->pool_list;
35: PetscFunctionBegin;
36: n = n->next;
37: while (n) {
38: PetscCDArrNd *lstn = n;
40: n = n->next;
41: PetscCall(PetscFree(lstn));
42: }
43: PetscCall(PetscFree(ail->pool_list.array));
44: PetscCall(PetscFree(ail->array));
45: PetscCall(MatDestroy(&ail->mat));
46: /* delete this (+agg+pool array) */
47: PetscCall(PetscFree(ail));
48: PetscFunctionReturn(PETSC_SUCCESS);
49: }
51: /* PetscCDSetChunkSize
52: */
53: PetscErrorCode PetscCDSetChunkSize(PetscCoarsenData *ail, PetscInt a_sz)
54: {
55: PetscFunctionBegin;
56: ail->chk_sz = a_sz;
57: PetscFunctionReturn(PETSC_SUCCESS);
58: }
60: /* PetscCDGetNewNode
61: */
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;
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;
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: }
98: /* PetscCDIntNdSetID
99: */
100: PetscErrorCode PetscCDIntNdSetID(PetscCDIntNd *a_this, PetscInt a_id)
101: {
102: PetscFunctionBegin;
103: a_this->gid = a_id;
104: PetscFunctionReturn(PETSC_SUCCESS);
105: }
107: /* PetscCDIntNdGetID
108: */
109: PetscErrorCode PetscCDIntNdGetID(const PetscCDIntNd *a_this, PetscInt *a_gid)
110: {
111: PetscFunctionBegin;
112: *a_gid = a_this->gid;
113: PetscFunctionReturn(PETSC_SUCCESS);
114: }
116: /* PetscCDGetHeadPos
117: */
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: }
126: /* PetscCDGetNextPos
127: */
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: }
136: /* PetscCDAppendID
137: */
138: PetscErrorCode PetscCDAppendID(PetscCoarsenData *ail, PetscInt a_idx, PetscInt a_id)
139: {
140: PetscCDIntNd *n, *n2;
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: }
160: /* PetscCDAppendNode
161: */
162: PetscErrorCode PetscCDAppendNode(PetscCoarsenData *ail, PetscInt a_idx, PetscCDIntNd *a_n)
163: {
164: PetscCDIntNd *n2;
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: }
183: /* PetscCDRemoveNextNode: a_last->next, this exposes single linked list structure to API (not used)
184: */
185: PetscErrorCode PetscCDRemoveNextNode(PetscCoarsenData *ail, PetscInt a_idx, PetscCDIntNd *a_last)
186: {
187: PetscCDIntNd *del;
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: }
199: /* PetscCDPrint
200: */
201: PetscErrorCode PetscCDPrint(const PetscCoarsenData *ail, PetscInt Istart, MPI_Comm comm)
202: {
203: PetscCDIntNd *n, *n2;
204: PetscInt ii;
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: }
220: /* PetscCDMoveAppend - take list in a_srcidx and appends to destidx
221: */
222: PetscErrorCode PetscCDMoveAppend(PetscCoarsenData *ail, PetscInt a_destidx, PetscInt a_srcidx)
223: {
224: PetscCDIntNd *n;
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: }
245: /* PetscCDRemoveAllAt - empty one list and move data to cache
246: */
247: PetscErrorCode PetscCDRemoveAllAt(PetscCoarsenData *ail, PetscInt a_idx)
248: {
249: PetscCDIntNd *rem, *n1;
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: }
263: /* PetscCDCountAt
264: */
265: PetscErrorCode PetscCDCountAt(const PetscCoarsenData *ail, PetscInt a_idx, PetscInt *a_sz)
266: {
267: PetscCDIntNd *n1;
268: PetscInt sz = 0;
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: }
281: /* PetscCDSize
282: */
283: PetscErrorCode PetscCDCount(const PetscCoarsenData *ail, PetscInt *a_sz)
284: {
285: PetscInt sz = 0;
287: PetscFunctionBegin;
288: for (PetscInt ii = 0; ii < ail->size; ii++) {
289: PetscCDIntNd *n1 = ail->array[ii];
291: while (n1) {
292: n1 = n1->next;
293: sz++;
294: }
295: }
296: *a_sz = sz;
297: PetscFunctionReturn(PETSC_SUCCESS);
298: }
300: /* PetscCDIsEmptyAt - Is the list empty? (not used)
301: */
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: }
310: /* PetscCDGetNonemptyIS - used for C-F methods
311: */
312: PetscErrorCode PetscCDGetNonemptyIS(PetscCoarsenData *ail, IS *a_mis)
313: {
314: PetscCDIntNd *n;
315: PetscInt ii, kk;
316: PetscInt *permute;
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: }
332: /* PetscCDGetMat
333: */
334: PetscErrorCode PetscCDGetMat(PetscCoarsenData *ail, Mat *a_mat)
335: {
336: PetscFunctionBegin;
337: *a_mat = ail->mat;
338: PetscFunctionReturn(PETSC_SUCCESS);
339: }
341: /* PetscCDSetMat
342: */
343: PetscErrorCode PetscCDSetMat(PetscCoarsenData *ail, Mat a_mat)
344: {
345: PetscFunctionBegin;
346: PetscCall(MatDestroy(&ail->mat)); //should not happen
347: ail->mat = a_mat;
348: PetscFunctionReturn(PETSC_SUCCESS);
349: }
351: /* PetscCDClearMat
352: */
353: PetscErrorCode PetscCDClearMat(PetscCoarsenData *ail)
354: {
355: PetscFunctionBegin;
356: ail->mat = NULL;
357: PetscFunctionReturn(PETSC_SUCCESS);
358: }
360: /* PetscCDGetASMBlocks - get IS of aggregates for ASM smoothers
361: */
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;
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: }
391: /* edge for priority queue */
392: typedef struct edge_tag {
393: PetscReal weight;
394: PetscInt lid0, gid1, ghost1_idx;
395: } Edge;
397: #define MY_MEPS (PETSC_MACHINE_EPSILON * 100)
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: }
404: /*
405: MatCoarsenApply_HEM_private - parallel heavy edge matching
407: Input Parameter:
408: . a_Gmat - global matrix of the graph
409: . n_iter - number of matching iterations
410: . threshold - threshold for filtering graphs
412: Output Parameter:
413: . a_locals_llist - array of list of local nodes rooted at local node
414: */
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;
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));
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
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;
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;
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));
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: PetscCall(VecDestroy(&vec));
526: /* lghost_gid: the global ID of each ghost is given by garray */
527: PetscCall(PetscArraycpy(lghost_gid, mpimat->garray, num_ghosts));
528: }
529: // get 'comm_procs' (could hoist)
530: for (PetscInt kk = 0; kk < REQ_BF_SIZE; kk++) comm_procs[kk] = -1;
531: for (ix = 0, ncomm_procs = 0; ix < num_ghosts; ix++) {
532: PetscMPIInt proc = lghost_pe[ix], idx = -1;
534: for (PetscMPIInt k = 0; k < ncomm_procs && idx == -1; k++)
535: if (comm_procs[k] == proc) idx = k;
536: if (idx == -1) comm_procs[ncomm_procs++] = proc;
537: PetscCheck(ncomm_procs != REQ_BF_SIZE, PETSC_COMM_SELF, PETSC_ERR_SUP, "Receive request array too small: %d", ncomm_procs);
538: }
539: /* count edges, compute initial 'locMaxEdge', 'locMaxPE' */
540: nEdges0 = 0;
541: for (PetscInt kk = 0, gid = Istart; kk < nloc; kk++, gid++) {
542: PetscReal max_e = 0., tt;
543: PetscScalar vval;
544: PetscInt lid = kk, max_pe = rank, pe, n;
546: ii = matA->i;
547: n = ii[lid + 1] - ii[lid];
548: aj = PetscSafePointerPlusOffset(matA->j, ii[lid]);
549: ap = PetscSafePointerPlusOffset(matA->a, ii[lid]);
550: for (PetscInt jj = 0; jj < n; jj++) {
551: PetscInt lidj = aj[jj];
553: if ((tt = PetscRealPart(ap[jj])) > threshold && lidj != lid) {
554: if (tt > max_e) max_e = tt;
555: if (lidj > lid) nEdges0++;
556: }
557: }
558: if ((ix = lid_cprowID[lid]) != -1) { /* if I have any ghost neighbors */
559: ii = matB->compressedrow.i;
560: n = ii[ix + 1] - ii[ix];
561: ap = matB->a + ii[ix];
562: aj = matB->j + ii[ix];
563: for (PetscInt jj = 0; jj < n; jj++) {
564: if ((tt = PetscRealPart(ap[jj])) > threshold) {
565: if (tt > max_e) max_e = tt;
566: nEdges0++;
567: if ((pe = lghost_pe[aj[jj]]) > max_pe) max_pe = pe;
568: }
569: }
570: }
571: vval = max_e;
572: PetscCall(VecSetValues(locMaxEdge, 1, &gid, &vval, INSERT_VALUES));
573: vval = (PetscScalar)max_pe;
574: PetscCall(VecSetValues(locMaxPE, 1, &gid, &vval, INSERT_VALUES));
575: if (iter == 0 && max_e <= MY_MEPS) { // add BCs to fake aggregate
576: lid_matched[lid] = PETSC_TRUE;
577: if (bc_agg == -1) {
578: bc_agg = lid;
579: PetscCall(PetscCDCreate(1, &bc_list));
580: }
581: PetscCall(PetscCDRemoveAllAt(agg_llists, lid));
582: PetscCall(PetscCDAppendID(bc_list, 0, Istart + lid));
583: }
584: }
585: PetscCall(VecAssemblyBegin(locMaxEdge));
586: PetscCall(VecAssemblyEnd(locMaxEdge));
587: PetscCall(VecAssemblyBegin(locMaxPE));
588: PetscCall(VecAssemblyEnd(locMaxPE));
589: /* make 'ghostMaxEdge_max_ew', 'lghost_max_pe' */
590: if (isMPI) {
591: const PetscScalar *buf;
593: PetscCall(VecDuplicate(mpimat->lvec, &ghostMaxEdge));
594: PetscCall(VecScatterBegin(mpimat->Mvctx, locMaxEdge, ghostMaxEdge, INSERT_VALUES, SCATTER_FORWARD));
595: PetscCall(VecScatterEnd(mpimat->Mvctx, locMaxEdge, ghostMaxEdge, INSERT_VALUES, SCATTER_FORWARD));
597: PetscCall(VecDuplicate(mpimat->lvec, &ghostMaxPE));
598: PetscCall(VecScatterBegin(mpimat->Mvctx, locMaxPE, ghostMaxPE, INSERT_VALUES, SCATTER_FORWARD));
599: PetscCall(VecScatterEnd(mpimat->Mvctx, locMaxPE, ghostMaxPE, INSERT_VALUES, SCATTER_FORWARD));
600: PetscCall(VecGetArrayRead(ghostMaxPE, &buf));
601: for (PetscInt kk = 0; kk < num_ghosts; kk++) lghost_max_pe[kk] = (PetscMPIInt)PetscRealPart(buf[kk]); // the MAX proc of the ghost now
602: PetscCall(VecRestoreArrayRead(ghostMaxPE, &buf));
603: }
604: { // make lid_max_pe
605: const PetscScalar *buf;
607: PetscCall(VecGetArrayRead(locMaxPE, &buf));
608: for (PetscInt kk = 0; kk < nloc; kk++) lid_max_pe[kk] = (PetscMPIInt)PetscRealPart(buf[kk]); // the MAX proc of the ghost now
609: PetscCall(VecRestoreArrayRead(locMaxPE, &buf));
610: }
611: /* setup sorted list of edges, and make 'Edges' */
612: PetscCall(PetscMalloc1(nEdges0, &Edges));
613: nEdges = 0;
614: for (PetscInt kk = 0, n; kk < nloc; kk++) {
615: const PetscInt lid = kk;
616: PetscReal tt;
618: ii = matA->i;
619: n = ii[lid + 1] - ii[lid];
620: aj = PetscSafePointerPlusOffset(matA->j, ii[lid]);
621: ap = PetscSafePointerPlusOffset(matA->a, ii[lid]);
622: for (PetscInt jj = 0; jj < n; jj++) {
623: PetscInt lidj = aj[jj];
625: if ((tt = PetscRealPart(ap[jj])) > threshold && lidj != lid) {
626: if (lidj > lid) {
627: Edges[nEdges].lid0 = lid;
628: Edges[nEdges].gid1 = lidj + Istart;
629: Edges[nEdges].ghost1_idx = -1;
630: Edges[nEdges].weight = tt;
631: nEdges++;
632: }
633: }
634: }
635: if ((ix = lid_cprowID[lid]) != -1) { /* if I have any ghost neighbor */
636: ii = matB->compressedrow.i;
637: n = ii[ix + 1] - ii[ix];
638: ap = matB->a + ii[ix];
639: aj = matB->j + ii[ix];
640: for (PetscInt jj = 0; jj < n; jj++) {
641: if ((tt = PetscRealPart(ap[jj])) > threshold) {
642: Edges[nEdges].lid0 = lid;
643: Edges[nEdges].gid1 = lghost_gid[aj[jj]];
644: Edges[nEdges].ghost1_idx = aj[jj];
645: Edges[nEdges].weight = tt;
646: nEdges++;
647: }
648: }
649: }
650: }
651: PetscCheck(nEdges == nEdges0, PETSC_COMM_SELF, PETSC_ERR_SUP, "nEdges != nEdges0: %" PetscInt_FMT " %" PetscInt_FMT, nEdges0, nEdges);
652: if (Edges) qsort(Edges, nEdges, sizeof(Edge), gamg_hem_compare);
654: PetscCall(PetscInfo(info_is, "[%d] HEM iteration %" PetscInt_FMT " with %" PetscInt_FMT " edges\n", rank, iter, nEdges));
656: /* projection matrix */
657: PetscCall(MatCreate(comm, &P));
658: PetscCall(MatSetType(P, MATAIJ));
659: PetscCall(MatSetSizes(P, nloc, nloc, PETSC_DETERMINE, PETSC_DETERMINE));
660: PetscCall(MatMPIAIJSetPreallocation(P, 1, NULL, 1, NULL));
661: PetscCall(MatSeqAIJSetPreallocation(P, 1, NULL));
662: PetscCall(MatSetUp(P));
663: /* process - communicate - process */
664: for (PetscInt sub_it = 0, old_num_edge = 0; /* sub_it < n_sub_its */; /* sub_it++ */) {
665: PetscInt nactive_edges = 0, gn_act_n[3];
666: PetscMPIInt tag1, tag2;
668: PetscCall(VecGetArrayRead(locMaxEdge, &lid_max_ew));
669: if (isMPI) {
670: PetscCall(VecGetArrayRead(ghostMaxEdge, &lghost_max_ew));
671: PetscCall(PetscCommGetNewTag(comm, &tag1));
672: PetscCall(PetscCommGetNewTag(comm, &tag2));
673: }
674: for (PetscInt kk = 0; kk < nEdges; kk++) {
675: const Edge *e = &Edges[kk];
676: const PetscInt lid0 = e->lid0, gid1 = e->gid1, ghost1_idx = e->ghost1_idx, gid0 = lid0 + Istart, lid1 = gid1 - Istart;
677: PetscBool isOK = PETSC_TRUE, print = PETSC_FALSE;
679: if (print)
680: 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"));
681: /* skip if either vertex is matched already */
682: if (lid_matched[lid0] || (ghost1_idx != -1 && lghost_matched[ghost1_idx]) || (ghost1_idx == -1 && lid_matched[lid1])) continue;
684: nactive_edges++;
685: 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]));
686: 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)));
687: // smaller edge, lid_max_ew get updated - e0
688: if (PetscRealPart(lid_max_ew[lid0]) > e->weight + MY_MEPS) {
689: if (print)
690: 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]),
691: (double)e->weight));
692: continue; // we are basically filter edges here
693: }
694: // e1 - local
695: if (ghost1_idx == -1) {
696: if (PetscRealPart(lid_max_ew[lid1]) > e->weight + MY_MEPS) {
697: if (print)
698: 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)));
699: continue; // we are basically filter edges here
700: }
701: } else { // e1 - ghost
702: /* see if edge might get matched on other proc */
703: PetscReal g_max_e1 = PetscRealPart(lghost_max_ew[ghost1_idx]);
705: if (print)
706: 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]),
707: (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]));
708: if (g_max_e1 > e->weight + MY_MEPS) {
709: /* 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 )); */
710: continue;
711: } else if (g_max_e1 >= e->weight - MY_MEPS && lghost_pe[ghost1_idx] > rank) { // is 'lghost_max_pe[ghost1_idx] > rank' needed?
712: /* check for max_ea == to this edge and larger processor that will deal with this */
713: if (print)
714: 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,
715: (double)e->weight));
716: continue;
717: } else {
718: /* 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 )); */
719: }
720: }
721: /* check ghost for v0 */
722: if (isOK) {
723: PetscReal max_e, ew;
725: if ((ix = lid_cprowID[lid0]) != -1) { /* if I have any ghost neighbors */
726: PetscInt n;
728: ii = matB->compressedrow.i;
729: n = ii[ix + 1] - ii[ix];
730: ap = matB->a + ii[ix];
731: aj = matB->j + ii[ix];
732: for (PetscInt jj = 0; jj < n && isOK; jj++) {
733: PetscInt lidj = aj[jj];
735: if (lghost_matched[lidj]) continue;
736: ew = PetscRealPart(ap[jj]);
737: if (ew <= threshold) continue;
738: max_e = PetscRealPart(lghost_max_ew[lidj]);
740: /* check for max_e == to this edge and larger processor that will deal with this */
741: if (ew >= PetscRealPart(lid_max_ew[lid0]) - MY_MEPS && lghost_max_pe[lidj] > rank) isOK = PETSC_FALSE;
742: 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]);
743: if (print)
744: 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));
745: }
746: if (!isOK && print) PetscCall(PetscSynchronizedPrintf(comm, "\t\t[%d] skip edge (%" PetscInt_FMT " %" PetscInt_FMT ") from ghost inspection\n", rank, gid0, gid1));
747: }
748: /* check local v1 */
749: if (ghost1_idx == -1) {
750: if ((ix = lid_cprowID[lid1]) != -1) { /* if I have any ghost neighbors */
751: PetscInt n;
753: ii = matB->compressedrow.i;
754: n = ii[ix + 1] - ii[ix];
755: ap = matB->a + ii[ix];
756: aj = matB->j + ii[ix];
757: for (PetscInt jj = 0; jj < n && isOK; jj++) {
758: PetscInt lidj = aj[jj];
760: if (lghost_matched[lidj]) continue;
761: ew = PetscRealPart(ap[jj]);
762: if (ew <= threshold) continue;
763: max_e = PetscRealPart(lghost_max_ew[lidj]);
764: /* check for max_e == to this edge and larger processor that will deal with this */
765: if (ew >= PetscRealPart(lid_max_ew[lid1]) - MY_MEPS && lghost_max_pe[lidj] > rank) isOK = PETSC_FALSE;
766: PetscCheck(ew <= max_e + MY_MEPS, PETSC_COMM_SELF, PETSC_ERR_SUP, "edge weight %e > max %e", (double)PetscRealPart(ew), (double)PetscRealPart(max_e));
767: if (print)
768: 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]));
769: }
770: }
771: if (!isOK && print) PetscCall(PetscSynchronizedPrintf(comm, "\t\t[%d] skip edge (%" PetscInt_FMT " %" PetscInt_FMT ") from ghost inspection\n", rank, gid0, gid1));
772: }
773: }
774: PetscReal e1_max_w = (ghost1_idx == -1 ? PetscRealPart(lid_max_ew[lid0]) : PetscRealPart(lghost_max_ew[ghost1_idx]));
775: if (print)
776: 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));
777: /* do it */
778: if (isOK) {
779: if (ghost1_idx == -1) {
780: PetscCheck(!lid_matched[lid1], PETSC_COMM_SELF, PETSC_ERR_SUP, "local %" PetscInt_FMT " is matched", gid1);
781: lid_matched[lid1] = PETSC_TRUE; /* keep track of what we've done this round */
782: PetscCall(PetscCDMoveAppend(agg_llists, lid0, lid1)); // takes lid1's list and appends to lid0's
783: } else {
784: /* add gid1 to list of ghost deleted by me -- I need their children */
785: PetscMPIInt proc = lghost_pe[ghost1_idx];
786: PetscCheck(!lghost_matched[ghost1_idx], PETSC_COMM_SELF, PETSC_ERR_SUP, "ghost %" PetscInt_FMT " is matched", lghost_gid[ghost1_idx]);
787: lghost_matched[ghost1_idx] = PETSC_TRUE;
788: PetscCall(PetscCDAppendID(ghost_deleted_list, proc, ghost1_idx)); /* cache to send messages */
789: PetscCall(PetscCDAppendID(ghost_deleted_list, proc, lid0));
790: }
791: lid_matched[lid0] = PETSC_TRUE; /* keep track of what we've done this round */
792: /* set projection */
793: PetscCall(MatSetValues(P, 1, &gid0, 1, &gid0, &one, INSERT_VALUES));
794: PetscCall(MatSetValues(P, 1, &gid1, 1, &gid0, &one, INSERT_VALUES));
795: //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));
796: } /* matched */
797: } /* edge loop */
798: PetscCall(PetscSynchronizedFlush(comm, PETSC_STDOUT));
799: if (isMPI) PetscCall(VecRestoreArrayRead(ghostMaxEdge, &lghost_max_ew));
800: PetscCall(VecRestoreArrayRead(locMaxEdge, &lid_max_ew));
801: // count active for test, latter, update deleted ghosts
802: gn_act_n[0] = nactive_edges;
803: if (ghost_deleted_list) PetscCall(PetscCDCount(ghost_deleted_list, &gn_act_n[2]));
804: else gn_act_n[2] = 0;
805: PetscCall(PetscCDCount(agg_llists, &gn_act_n[1]));
806: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, gn_act_n, 3, MPIU_INT, MPI_SUM, comm));
807: 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]));
808: /* deal with deleted ghost */
809: if (isMPI) {
810: PetscCDIntNd *pos;
811: PetscInt *sbuffs1[REQ_BF_SIZE], ndel;
812: PetscInt *sbuffs2[REQ_BF_SIZE];
813: MPI_Status status;
815: /* send deleted ghosts */
816: for (PetscInt proc_idx = 0; proc_idx < ncomm_procs; proc_idx++) {
817: const PetscMPIInt proc = comm_procs[proc_idx];
818: PetscInt *sbuff, *pt, scount;
819: MPI_Request *request;
821: /* count ghosts */
822: PetscCall(PetscCDCountAt(ghost_deleted_list, proc, &ndel));
823: ndel /= 2; // two entries for each proc
824: scount = 2 + 2 * ndel;
825: PetscCall(PetscMalloc1(scount + request_size, &sbuff));
826: /* save requests */
827: sbuffs1[proc_idx] = sbuff;
828: request = (MPI_Request *)sbuff;
829: sbuff = pt = sbuff + request_size;
830: /* write [ndel, proc, n*[gid1,gid0] */
831: *pt++ = ndel; // number of deleted to send
832: *pt++ = rank; // proc (not used)
833: PetscCall(PetscCDGetHeadPos(ghost_deleted_list, proc, &pos));
834: while (pos) {
835: PetscInt lid0, ghost_idx, gid1;
837: PetscCall(PetscCDIntNdGetID(pos, &ghost_idx));
838: gid1 = lghost_gid[ghost_idx];
839: PetscCall(PetscCDGetNextPos(ghost_deleted_list, proc, &pos));
840: PetscCall(PetscCDIntNdGetID(pos, &lid0));
841: PetscCall(PetscCDGetNextPos(ghost_deleted_list, proc, &pos));
842: *pt++ = gid1;
843: *pt++ = lid0 + Istart; // gid0
844: }
845: PetscCheck(pt - sbuff == scount, PETSC_COMM_SELF, PETSC_ERR_SUP, "sbuff-pt != scount: %zu", pt - sbuff);
846: /* MPIU_Isend: tag1 [ndel, proc, n*[gid1,gid0] ] */
847: PetscCallMPI(MPIU_Isend(sbuff, scount, MPIU_INT, proc, tag1, comm, request));
848: PetscCall(PetscCDRemoveAllAt(ghost_deleted_list, proc)); // done with this list
849: }
850: /* receive deleted, send back partial aggregates, clear lists */
851: for (PetscInt proc_idx = 0; proc_idx < ncomm_procs; proc_idx++) {
852: PetscCallMPI(MPI_Probe(comm_procs[proc_idx] /* MPI_ANY_SOURCE */, tag1, comm, &status));
853: {
854: PetscInt *pt, *pt2, *pt3, *sbuff, tmp;
855: MPI_Request *request;
856: PetscMPIInt rcount, scount;
857: const PetscMPIInt proc = status.MPI_SOURCE;
859: PetscCallMPI(MPI_Get_count(&status, MPIU_INT, &rcount));
860: if (rcount > rbuff_sz) {
861: PetscCall(PetscFree(rbuff));
862: PetscCall(PetscMalloc1(rcount, &rbuff));
863: rbuff_sz = rcount;
864: }
865: /* MPI_Recv: tag1 [ndel, proc, ndel*[gid1,gid0] ] */
866: PetscCallMPI(MPI_Recv(rbuff, rcount, MPIU_INT, proc, tag1, comm, &status));
867: /* read and count sends *[lid0, n, n*[gid] ] */
868: pt = rbuff;
869: scount = 0;
870: ndel = *pt++; // number of deleted to recv
871: tmp = *pt++; // proc (not used)
872: while (ndel--) {
873: PetscInt gid1 = *pt++, lid1 = gid1 - Istart;
874: PetscInt gh_gid0 = *pt++; // gid on other proc (not used here to count)
876: PetscCheck(lid1 >= 0 && lid1 < nloc, PETSC_COMM_SELF, PETSC_ERR_SUP, "received ghost deleted %" PetscInt_FMT, gid1);
877: 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);
878: lid_matched[lid1] = PETSC_TRUE; /* keep track of what we've done this round */
879: PetscCall(PetscCDCountAt(agg_llists, lid1, &tmp)); // n
880: scount += tmp + 2; // lid0, n, n*[gid]
881: }
882: PetscCheck((pt - rbuff) == (ptrdiff_t)rcount, PETSC_COMM_SELF, PETSC_ERR_SUP, "receive buffer size != num read: %zu; rcount: %d", pt - rbuff, rcount);
883: /* send tag2: *[gid0, n, n*[gid] ] */
884: PetscCall(PetscMalloc1(scount + request_size, &sbuff));
885: sbuffs2[proc_idx] = sbuff; /* cache request */
886: request = (MPI_Request *)sbuff;
887: pt2 = sbuff = sbuff + request_size;
888: // read again: n, proc, n*[gid1,gid0]
889: pt = rbuff;
890: ndel = *pt++;
891: tmp = *pt++; // proc (not used)
892: while (ndel--) {
893: PetscInt gid1 = *pt++, lid1 = gid1 - Istart, gh_gid0 = *pt++;
895: /* write [gid0, aggSz, aggSz[gid] ] */
896: *pt2++ = gh_gid0;
897: pt3 = pt2++; /* save pointer for later */
898: PetscCall(PetscCDGetHeadPos(agg_llists, lid1, &pos));
899: while (pos) {
900: PetscInt gid;
902: PetscCall(PetscCDIntNdGetID(pos, &gid));
903: PetscCall(PetscCDGetNextPos(agg_llists, lid1, &pos));
904: *pt2++ = gid;
905: }
906: PetscCall(PetscIntCast(pt2 - pt3 - 1, pt3));
907: /* clear list */
908: PetscCall(PetscCDRemoveAllAt(agg_llists, lid1));
909: }
910: PetscCheck((pt2 - sbuff) == (ptrdiff_t)scount, PETSC_COMM_SELF, PETSC_ERR_SUP, "buffer size != num write: %zu %d", pt2 - sbuff, scount);
911: /* MPIU_Isend: requested data tag2 *[lid0, n, n*[gid1] ] */
912: PetscCallMPI(MPIU_Isend(sbuff, scount, MPIU_INT, proc, tag2, comm, request));
913: }
914: } // proc_idx
915: /* receive tag2 *[gid0, n, n*[gid] ] */
916: for (PetscMPIInt proc_idx = 0; proc_idx < ncomm_procs; proc_idx++) {
917: PetscMPIInt proc;
918: PetscInt *pt;
919: int rcount;
921: PetscCallMPI(MPI_Probe(comm_procs[proc_idx] /* MPI_ANY_SOURCE */, tag2, comm, &status));
922: PetscCallMPI(MPI_Get_count(&status, MPIU_INT, &rcount));
923: if (rcount > rbuff_sz) {
924: PetscCall(PetscFree(rbuff));
925: PetscCall(PetscMalloc1(rcount, &rbuff));
926: rbuff_sz = rcount;
927: }
928: proc = status.MPI_SOURCE;
929: /* MPI_Recv: tag1 [n, proc, n*[gid1,lid0] ] */
930: PetscCallMPI(MPI_Recv(rbuff, rcount, MPIU_INT, proc, tag2, comm, &status));
931: pt = rbuff;
932: while (pt - rbuff < rcount) {
933: PetscInt gid0 = *pt++, n = *pt++;
935: while (n--) {
936: PetscInt gid1 = *pt++;
938: PetscCall(PetscCDAppendID(agg_llists, gid0 - Istart, gid1));
939: }
940: }
941: PetscCheck((pt - rbuff) == (ptrdiff_t)rcount, PETSC_COMM_SELF, PETSC_ERR_SUP, "recv buffer size != num read: %zu %d", pt - rbuff, rcount);
942: }
943: /* wait for tag1 isends */
944: for (PetscMPIInt proc_idx = 0; proc_idx < ncomm_procs; proc_idx++) {
945: MPI_Request *request = (MPI_Request *)sbuffs1[proc_idx];
947: PetscCallMPI(MPI_Wait(request, &status));
948: PetscCall(PetscFree(sbuffs1[proc_idx]));
949: }
950: /* wait for tag2 isends */
951: for (PetscMPIInt proc_idx = 0; proc_idx < ncomm_procs; proc_idx++) {
952: MPI_Request *request = (MPI_Request *)sbuffs2[proc_idx];
954: PetscCallMPI(MPI_Wait(request, &status));
955: PetscCall(PetscFree(sbuffs2[proc_idx]));
956: }
957: } /* MPI */
958: /* set 'lghost_matched' - use locMaxEdge, ghostMaxEdge (recomputed next) */
959: if (isMPI) {
960: const PetscScalar *sbuff;
962: for (PetscInt kk = 0, gid = Istart; kk < nloc; kk++, gid++) {
963: PetscScalar vval = lid_matched[kk] ? 1.0 : 0.0;
965: PetscCall(VecSetValues(locMaxEdge, 1, &gid, &vval, INSERT_VALUES)); /* set with GID */
966: }
967: PetscCall(VecAssemblyBegin(locMaxEdge));
968: PetscCall(VecAssemblyEnd(locMaxEdge));
969: PetscCall(VecScatterBegin(mpimat->Mvctx, locMaxEdge, ghostMaxEdge, INSERT_VALUES, SCATTER_FORWARD));
970: PetscCall(VecScatterEnd(mpimat->Mvctx, locMaxEdge, ghostMaxEdge, INSERT_VALUES, SCATTER_FORWARD));
971: PetscCall(VecGetArrayRead(ghostMaxEdge, &sbuff));
972: for (PetscInt kk = 0; kk < num_ghosts; kk++) lghost_matched[kk] = (PetscBool)(PetscRealPart(sbuff[kk]) != 0.0);
973: PetscCall(VecRestoreArrayRead(ghostMaxEdge, &sbuff));
974: }
975: /* compute 'locMaxEdge' inside sub iteration b/c max weight can drop as neighbors are matched */
976: for (PetscInt kk = 0, gid = Istart; kk < nloc; kk++, gid++) {
977: PetscReal max_e = 0., tt;
978: PetscScalar vval;
979: const PetscInt lid = kk;
980: PetscMPIInt max_pe = rank, pe, n;
982: ii = matA->i;
983: PetscCall(PetscMPIIntCast(ii[lid + 1] - ii[lid], &n));
984: aj = PetscSafePointerPlusOffset(matA->j, ii[lid]);
985: ap = PetscSafePointerPlusOffset(matA->a, ii[lid]);
986: for (PetscMPIInt jj = 0; jj < n; jj++) {
987: PetscInt lidj = aj[jj];
989: if (lid_matched[lidj]) continue; /* this is new - can change local max */
990: if (lidj != lid && PetscRealPart(ap[jj]) > max_e) max_e = PetscRealPart(ap[jj]);
991: }
992: if (lid_cprowID && (ix = lid_cprowID[lid]) != -1) { /* if I have any ghost neighbors */
993: ii = matB->compressedrow.i;
994: PetscCall(PetscMPIIntCast(ii[ix + 1] - ii[ix], &n));
995: ap = matB->a + ii[ix];
996: aj = matB->j + ii[ix];
997: for (PetscMPIInt jj = 0; jj < n; jj++) {
998: PetscInt lidj = aj[jj];
1000: if (lghost_matched[lidj]) continue;
1001: if ((tt = PetscRealPart(ap[jj])) > max_e) max_e = tt;
1002: }
1003: }
1004: vval = max_e;
1005: PetscCall(VecSetValues(locMaxEdge, 1, &gid, &vval, INSERT_VALUES)); /* set with GID */
1006: // max PE with max edge
1007: if (lid_cprowID && (ix = lid_cprowID[lid]) != -1) { /* if I have any ghost neighbors */
1008: ii = matB->compressedrow.i;
1009: PetscCall(PetscMPIIntCast(ii[ix + 1] - ii[ix], &n));
1010: ap = matB->a + ii[ix];
1011: aj = matB->j + ii[ix];
1012: for (PetscInt jj = 0; jj < n; jj++) {
1013: PetscInt lidj = aj[jj];
1015: if (lghost_matched[lidj]) continue;
1016: if ((pe = lghost_pe[aj[jj]]) > max_pe && PetscRealPart(ap[jj]) >= max_e - MY_MEPS) max_pe = pe;
1017: }
1018: }
1019: vval = max_pe;
1020: PetscCall(VecSetValues(locMaxPE, 1, &gid, &vval, INSERT_VALUES));
1021: }
1022: PetscCall(VecAssemblyBegin(locMaxEdge));
1023: PetscCall(VecAssemblyEnd(locMaxEdge));
1024: PetscCall(VecAssemblyBegin(locMaxPE));
1025: PetscCall(VecAssemblyEnd(locMaxPE));
1026: /* compute 'lghost_max_ew' and 'lghost_max_pe' to get ready for next iteration*/
1027: if (isMPI) {
1028: const PetscScalar *buf;
1030: PetscCall(VecScatterBegin(mpimat->Mvctx, locMaxEdge, ghostMaxEdge, INSERT_VALUES, SCATTER_FORWARD));
1031: PetscCall(VecScatterEnd(mpimat->Mvctx, locMaxEdge, ghostMaxEdge, INSERT_VALUES, SCATTER_FORWARD));
1032: PetscCall(VecScatterBegin(mpimat->Mvctx, locMaxPE, ghostMaxPE, INSERT_VALUES, SCATTER_FORWARD));
1033: PetscCall(VecScatterEnd(mpimat->Mvctx, locMaxPE, ghostMaxPE, INSERT_VALUES, SCATTER_FORWARD));
1034: PetscCall(VecGetArrayRead(ghostMaxPE, &buf));
1035: for (PetscInt kk = 0; kk < num_ghosts; kk++) lghost_max_pe[kk] = (PetscMPIInt)PetscRealPart(buf[kk]); // the MAX proc of the ghost now
1036: PetscCall(VecRestoreArrayRead(ghostMaxPE, &buf));
1037: }
1038: // if no active edges, stop
1039: if (gn_act_n[0] < 1) break;
1040: // inc and check (self stopping iteration
1041: 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);
1042: sub_it++;
1043: 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);
1044: old_num_edge = gn_act_n[0];
1045: } /* sub_it loop */
1046: /* clean up iteration */
1047: PetscCall(PetscFree(Edges));
1048: if (isMPI) { // can be hoisted
1049: PetscCall(VecRestoreArrayRead(ghostMaxEdge, &lghost_max_ew));
1050: PetscCall(VecDestroy(&ghostMaxEdge));
1051: PetscCall(VecDestroy(&ghostMaxPE));
1052: PetscCall(PetscFree4(lghost_matched, lghost_pe, lghost_gid, lghost_max_pe));
1053: }
1054: PetscCall(VecDestroy(&locMaxEdge));
1055: PetscCall(VecDestroy(&locMaxPE));
1056: /* create next graph */
1057: {
1058: Vec diag;
1060: /* add identity for unmatched vertices so they stay alive */
1061: for (PetscInt kk = 0, gid1, gid = Istart; kk < nloc; kk++, gid++) {
1062: if (!lid_matched[kk]) {
1063: const PetscInt lid = kk;
1064: PetscCDIntNd *pos;
1066: PetscCall(PetscCDGetHeadPos(agg_llists, lid, &pos));
1067: PetscCheck(pos, PETSC_COMM_SELF, PETSC_ERR_PLIB, "empty list in singleton: %" PetscInt_FMT, gid);
1068: PetscCall(PetscCDIntNdGetID(pos, &gid1));
1069: PetscCheck(gid1 == gid, PETSC_COMM_SELF, PETSC_ERR_PLIB, "first in list (%" PetscInt_FMT ") in singleton not %" PetscInt_FMT, gid1, gid);
1070: PetscCall(MatSetValues(P, 1, &gid, 1, &gid, &one, INSERT_VALUES));
1071: }
1072: }
1073: PetscCall(MatAssemblyBegin(P, MAT_FINAL_ASSEMBLY));
1074: PetscCall(MatAssemblyEnd(P, MAT_FINAL_ASSEMBLY));
1076: /* project to make new graph with collapsed edges */
1077: PetscCall(MatPtAP(cMat, P, MAT_INITIAL_MATRIX, 1.0, &tMat));
1078: PetscCall(MatDestroy(&P));
1079: PetscCall(MatDestroy(&cMat));
1080: cMat = tMat;
1081: PetscCall(MatCreateVecs(cMat, &diag, NULL));
1082: PetscCall(MatGetDiagonal(cMat, diag));
1083: PetscCall(VecReciprocal(diag));
1084: PetscCall(VecSqrtAbs(diag));
1085: PetscCall(MatDiagonalScale(cMat, diag, diag));
1086: PetscCall(VecDestroy(&diag));
1087: }
1088: } /* coarsen iterator */
1090: /* make fake matrix with Mat->B only for smoothed agg QR. Need this if we make an aux graph (ie, PtAP) with k > 1 */
1091: if (size > 1) {
1092: Mat mat;
1093: PetscCDIntNd *pos;
1094: PetscInt NN, MM, jj = 0, mxsz = 0;
1096: for (PetscInt kk = 0; kk < nloc; kk++) {
1097: PetscCall(PetscCDCountAt(agg_llists, kk, &jj));
1098: if (jj > mxsz) mxsz = jj;
1099: }
1100: PetscCall(MatGetSize(a_Gmat, &MM, &NN));
1101: if (mxsz > MM - nloc) mxsz = MM - nloc;
1102: /* matrix of ghost adj for square graph */
1103: PetscCall(MatCreateAIJ(comm, nloc, nloc, PETSC_DETERMINE, PETSC_DETERMINE, 0, NULL, mxsz, NULL, &mat));
1104: for (PetscInt lid = 0, gid = Istart; lid < nloc; lid++, gid++) {
1105: PetscCall(PetscCDGetHeadPos(agg_llists, lid, &pos));
1106: while (pos) {
1107: PetscInt gid1;
1109: PetscCall(PetscCDIntNdGetID(pos, &gid1));
1110: PetscCall(PetscCDGetNextPos(agg_llists, lid, &pos));
1111: if (gid1 < Istart || gid1 >= Istart + nloc) PetscCall(MatSetValues(mat, 1, &gid, 1, &gid1, &one, ADD_VALUES));
1112: }
1113: }
1114: PetscCall(MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY));
1115: PetscCall(MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY));
1116: PetscCall(PetscCDSetMat(agg_llists, mat));
1117: PetscCall(PetscCDDestroy(ghost_deleted_list));
1118: if (rbuff_sz) PetscCall(PetscFree(rbuff)); // always true
1119: }
1120: // move BCs into some node
1121: if (bc_list) {
1122: PetscCDIntNd *pos;
1124: PetscCall(PetscCDGetHeadPos(bc_list, 0, &pos));
1125: while (pos) {
1126: PetscInt gid1;
1128: PetscCall(PetscCDIntNdGetID(pos, &gid1));
1129: PetscCall(PetscCDGetNextPos(bc_list, 0, &pos));
1130: PetscCall(PetscCDAppendID(agg_llists, bc_agg, gid1));
1131: }
1132: PetscCall(PetscCDRemoveAllAt(bc_list, 0));
1133: PetscCall(PetscCDDestroy(bc_list));
1134: }
1135: {
1136: // check sizes -- all vertices must get in graph
1137: PetscInt globalsz, MM;
1139: PetscCall(MatGetSize(a_Gmat, &MM, NULL));
1140: PetscCall(PetscCDCount(agg_llists, &globalsz));
1141: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &globalsz, 1, MPIU_INT, MPI_SUM, comm));
1142: PetscCheck(MM == globalsz, comm, PETSC_ERR_SUP, "lost %" PetscInt_FMT " equations ?", MM - globalsz);
1143: }
1144: // cleanup
1145: PetscCall(MatDestroy(&cMat));
1146: PetscCall(PetscFree3(lid_matched, lid_cprowID, lid_max_pe));
1147: PetscCall(ISDestroy(&info_is));
1148: PetscFunctionReturn(PETSC_SUCCESS);
1149: }
1151: /*
1152: HEM coarsen, simple greedy.
1153: */
1154: static PetscErrorCode MatCoarsenApply_HEM(MatCoarsen coarse)
1155: {
1156: Mat mat = coarse->graph;
1158: PetscFunctionBegin;
1159: PetscCall(MatCoarsenApply_HEM_private(mat, coarse->max_it, coarse->threshold, &coarse->agg_lists));
1160: PetscFunctionReturn(PETSC_SUCCESS);
1161: }
1163: static PetscErrorCode MatCoarsenView_HEM(MatCoarsen coarse, PetscViewer viewer)
1164: {
1165: PetscMPIInt rank;
1166: PetscBool isascii;
1168: PetscFunctionBegin;
1169: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)coarse), &rank));
1170: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
1171: if (isascii) {
1172: PetscCDIntNd *pos, *pos2;
1173: PetscViewerFormat format;
1175: PetscCall(PetscViewerASCIIPrintf(viewer, "%" PetscInt_FMT " matching steps with threshold = %g\n", coarse->max_it, (double)coarse->threshold));
1176: PetscCall(PetscViewerGetFormat(viewer, &format));
1177: if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
1178: if (coarse->agg_lists) {
1179: PetscCall(PetscViewerASCIIPushSynchronized(viewer));
1180: for (PetscInt kk = 0; kk < coarse->agg_lists->size; kk++) {
1181: PetscCall(PetscCDGetHeadPos(coarse->agg_lists, kk, &pos));
1182: if ((pos2 = pos)) PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "selected local %" PetscInt_FMT ": ", kk));
1183: while (pos) {
1184: PetscInt gid1;
1186: PetscCall(PetscCDIntNdGetID(pos, &gid1));
1187: PetscCall(PetscCDGetNextPos(coarse->agg_lists, kk, &pos));
1188: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, " %" PetscInt_FMT " ", gid1));
1189: }
1190: if (pos2) PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "\n"));
1191: }
1192: PetscCall(PetscViewerFlush(viewer));
1193: PetscCall(PetscViewerASCIIPopSynchronized(viewer));
1194: } else {
1195: PetscCall(PetscViewerASCIIPrintf(viewer, " HEM aggregator lists are not available\n"));
1196: }
1197: }
1198: }
1199: PetscFunctionReturn(PETSC_SUCCESS);
1200: }
1202: /*MC
1203: MATCOARSENHEM - A coarsener that uses HEM a simple greedy coarsener
1205: Level: beginner
1207: .seealso: `MatCoarsen`, `MatCoarsenMISKSetDistance()`, `MatCoarsenApply()`, `MatCoarsenSetType()`, `MatCoarsenType`, `MatCoarsenCreate()`, `MATCOARSENMISK`, `MATCOARSENMIS`
1208: M*/
1210: PETSC_EXTERN PetscErrorCode MatCoarsenCreate_HEM(MatCoarsen coarse)
1211: {
1212: PetscFunctionBegin;
1213: coarse->ops->apply = MatCoarsenApply_HEM;
1214: coarse->ops->view = MatCoarsenView_HEM;
1215: coarse->max_it = 4;
1216: PetscFunctionReturn(PETSC_SUCCESS);
1217: }