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: 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: }
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: if (ail->mat) 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: /* lghost_gid */
526: for (PetscInt kk = 0, gid = Istart; kk < nloc; kk++, gid++) {
527: vval = (PetscScalar)gid;
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;
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;
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];
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;
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));
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;
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;
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];
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);
665: PetscCall(PetscInfo(info_is, "[%d] HEM iteration %" PetscInt_FMT " with %" PetscInt_FMT " edges\n", rank, iter, nEdges));
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;
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;
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;
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]);
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;
736: if ((ix = lid_cprowID[lid0]) != -1) { /* if I have any ghost neighbors */
737: PetscInt n;
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];
746: if (lghost_matched[lidj]) continue;
747: ew = PetscRealPart(ap[jj]);
748: if (ew <= threshold) continue;
749: max_e = PetscRealPart(lghost_max_ew[lidj]);
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;
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];
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;
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;
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;
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;
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)
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++;
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;
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;
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++;
946: while (n--) {
947: PetscInt gid1 = *pt++;
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];
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];
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;
973: for (PetscInt kk = 0, gid = Istart; kk < nloc; kk++, gid++) {
974: PetscScalar vval = lid_matched[kk] ? 1.0 : 0.0;
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;
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];
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];
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];
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;
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;
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;
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));
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 */
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;
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;
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;
1135: PetscCall(PetscCDGetHeadPos(bc_list, 0, &pos));
1136: while (pos) {
1137: PetscInt gid1;
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;
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: }
1162: /*
1163: HEM coarsen, simple greedy.
1164: */
1165: static PetscErrorCode MatCoarsenApply_HEM(MatCoarsen coarse)
1166: {
1167: Mat mat = coarse->graph;
1169: PetscFunctionBegin;
1170: PetscCall(MatCoarsenApply_HEM_private(mat, coarse->max_it, coarse->threshold, &coarse->agg_lists));
1171: PetscFunctionReturn(PETSC_SUCCESS);
1172: }
1174: static PetscErrorCode MatCoarsenView_HEM(MatCoarsen coarse, PetscViewer viewer)
1175: {
1176: PetscMPIInt rank;
1177: PetscBool isascii;
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;
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;
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: }
1213: /*MC
1214: MATCOARSENHEM - A coarsener that uses HEM a simple greedy coarsener
1216: Level: beginner
1218: .seealso: `MatCoarsen`, `MatCoarsenMISKSetDistance()`, `MatCoarsenApply()`, `MatCoarsenSetType()`, `MatCoarsenType`, `MatCoarsenCreate()`, `MATCOARSENMISK`, `MATCOARSENMIS`
1219: M*/
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: }