Actual source code: hem.c
2: #include <petsc/private/matimpl.h>
3: #include <../src/mat/impls/aij/seq/aij.h>
4: #include <../src/mat/impls/aij/mpi/mpiaij.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: /* NPDestroy
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;
39: n = n->next;
40: PetscCall(PetscFree(lstn));
41: }
42: if (ail->pool_list.array) PetscCall(PetscFree(ail->pool_list.array));
43: PetscCall(PetscFree(ail->array));
44: if (ail->mat) PetscCall(MatDestroy(&ail->mat));
45: /* delete this (+agg+pool array) */
46: PetscCall(PetscFree(ail));
47: PetscFunctionReturn(PETSC_SUCCESS);
48: }
50: /* PetscCDSetChuckSize
51: */
52: PetscErrorCode PetscCDSetChuckSize(PetscCoarsenData *ail, PetscInt a_sz)
53: {
54: PetscFunctionBegin;
55: ail->chk_sz = a_sz;
56: PetscFunctionReturn(PETSC_SUCCESS);
57: }
59: /* PetscCDGetNewNode
60: */
61: PetscErrorCode PetscCDGetNewNode(PetscCoarsenData *ail, PetscCDIntNd **a_out, PetscInt a_id)
62: {
63: PetscFunctionBegin;
64: *a_out = NULL; /* squelch -Wmaybe-uninitialized */
65: if (ail->extra_nodes) {
66: PetscCDIntNd *node = ail->extra_nodes;
67: ail->extra_nodes = node->next;
68: node->gid = a_id;
69: node->next = NULL;
70: *a_out = node;
71: } else {
72: if (!ail->pool_list.array) {
73: if (!ail->chk_sz) ail->chk_sz = 10; /* use a chuck size of ail->size? */
74: PetscCall(PetscMalloc1(ail->chk_sz, &ail->pool_list.array));
75: ail->new_node = ail->pool_list.array;
76: ail->new_left = ail->chk_sz;
77: ail->new_node->next = NULL;
78: } else if (!ail->new_left) {
79: PetscCDArrNd *node;
80: PetscCall(PetscMalloc(ail->chk_sz * sizeof(PetscCDIntNd) + sizeof(PetscCDArrNd), &node));
81: node->array = (PetscCDIntNd *)(node + 1);
82: node->next = ail->pool_list.next;
83: ail->pool_list.next = node;
84: ail->new_left = ail->chk_sz;
85: ail->new_node = node->array;
86: }
87: ail->new_node->gid = a_id;
88: ail->new_node->next = NULL;
89: *a_out = ail->new_node++;
90: ail->new_left--;
91: }
92: PetscFunctionReturn(PETSC_SUCCESS);
93: }
95: /* PetscCDIntNdSetID
96: */
97: PetscErrorCode PetscCDIntNdSetID(PetscCDIntNd *a_this, PetscInt a_id)
98: {
99: PetscFunctionBegin;
100: a_this->gid = a_id;
101: PetscFunctionReturn(PETSC_SUCCESS);
102: }
104: /* PetscCDIntNdGetID
105: */
106: PetscErrorCode PetscCDIntNdGetID(const PetscCDIntNd *a_this, PetscInt *a_gid)
107: {
108: PetscFunctionBegin;
109: *a_gid = a_this->gid;
110: PetscFunctionReturn(PETSC_SUCCESS);
111: }
113: /* PetscCDGetHeadPos
114: */
115: PetscErrorCode PetscCDGetHeadPos(const PetscCoarsenData *ail, PetscInt a_idx, PetscCDIntNd **pos)
116: {
117: PetscFunctionBegin;
118: PetscCheck(a_idx < ail->size, PETSC_COMM_SELF, PETSC_ERR_PLIB, "a_idx >= ail->size: a_idx=%" PetscInt_FMT ".", a_idx);
119: *pos = ail->array[a_idx];
120: PetscFunctionReturn(PETSC_SUCCESS);
121: }
123: /* PetscCDGetNextPos
124: */
125: PetscErrorCode PetscCDGetNextPos(const PetscCoarsenData *ail, PetscInt l_idx, PetscCDIntNd **pos)
126: {
127: PetscFunctionBegin;
128: PetscCheck((*pos), PETSC_COMM_SELF, PETSC_ERR_PLIB, "NULL input position.");
129: *pos = (*pos)->next;
130: PetscFunctionReturn(PETSC_SUCCESS);
131: }
133: /* PetscCDAppendID
134: */
135: PetscErrorCode PetscCDAppendID(PetscCoarsenData *ail, PetscInt a_idx, PetscInt a_id)
136: {
137: PetscCDIntNd *n, *n2;
139: PetscFunctionBegin;
140: PetscCall(PetscCDGetNewNode(ail, &n, a_id));
141: PetscCheck(a_idx < ail->size, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Index %" PetscInt_FMT " out of range.", a_idx);
142: if (!(n2 = ail->array[a_idx])) ail->array[a_idx] = n;
143: else {
144: do {
145: if (!n2->next) {
146: n2->next = n;
147: PetscCheck(!n->next, PETSC_COMM_SELF, PETSC_ERR_PLIB, "n should not have a next");
148: break;
149: }
150: n2 = n2->next;
151: } while (n2);
152: PetscCheck(n2, PETSC_COMM_SELF, PETSC_ERR_PLIB, "n2 should be non-null");
153: }
154: PetscFunctionReturn(PETSC_SUCCESS);
155: }
157: /* PetscCDAppendNode
158: */
159: PetscErrorCode PetscCDAppendNode(PetscCoarsenData *ail, PetscInt a_idx, PetscCDIntNd *a_n)
160: {
161: PetscCDIntNd *n2;
163: PetscFunctionBegin;
164: PetscCheck(a_idx < ail->size, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Index %" PetscInt_FMT " out of range.", a_idx);
165: if (!(n2 = ail->array[a_idx])) ail->array[a_idx] = a_n;
166: else {
167: do {
168: if (!n2->next) {
169: n2->next = a_n;
170: a_n->next = NULL;
171: break;
172: }
173: n2 = n2->next;
174: } while (n2);
175: PetscCheck(n2, PETSC_COMM_SELF, PETSC_ERR_PLIB, "n2 should be non-null");
176: }
177: PetscFunctionReturn(PETSC_SUCCESS);
178: }
180: /* PetscCDRemoveNextNode: a_last->next, this exposes single linked list structure to API
181: */
182: PetscErrorCode PetscCDRemoveNextNode(PetscCoarsenData *ail, PetscInt a_idx, PetscCDIntNd *a_last)
183: {
184: PetscCDIntNd *del;
186: PetscFunctionBegin;
187: PetscCheck(a_idx < ail->size, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Index %" PetscInt_FMT " out of range.", a_idx);
188: PetscCheck(a_last->next, PETSC_COMM_SELF, PETSC_ERR_PLIB, "a_last should have a next");
189: del = a_last->next;
190: a_last->next = del->next;
191: /* del->next = NULL; -- this still used in a iterator so keep it intact -- need to fix this with a double linked list */
192: /* could reuse n2 but PetscCDAppendNode sometimes uses it */
193: PetscFunctionReturn(PETSC_SUCCESS);
194: }
196: /* PetscCDPrint
197: */
198: PetscErrorCode PetscCDPrint(const PetscCoarsenData *ail, MPI_Comm comm)
199: {
200: PetscCDIntNd *n;
201: PetscInt ii, kk;
202: PetscMPIInt rank;
204: PetscFunctionBegin;
205: PetscCallMPI(MPI_Comm_rank(comm, &rank));
206: for (ii = 0; ii < ail->size; ii++) {
207: kk = 0;
208: n = ail->array[ii];
209: if (n) PetscCall(PetscPrintf(comm, "[%d]%s list %" PetscInt_FMT ":\n", rank, PETSC_FUNCTION_NAME, ii));
210: while (n) {
211: PetscCall(PetscPrintf(comm, "\t[%d] %" PetscInt_FMT ") id %" PetscInt_FMT "\n", rank, ++kk, n->gid));
212: n = n->next;
213: }
214: }
215: PetscFunctionReturn(PETSC_SUCCESS);
216: }
218: /* PetscCDAppendRemove
219: */
220: PetscErrorCode PetscCDAppendRemove(PetscCoarsenData *ail, PetscInt a_destidx, PetscInt a_srcidx)
221: {
222: PetscCDIntNd *n;
224: PetscFunctionBegin;
225: PetscCheck(a_srcidx < ail->size, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Index %" PetscInt_FMT " out of range.", a_srcidx);
226: PetscCheck(a_destidx < ail->size, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Index %" PetscInt_FMT " out of range.", a_destidx);
227: PetscCheck(a_destidx != a_srcidx, PETSC_COMM_SELF, PETSC_ERR_PLIB, "a_destidx==a_srcidx %" PetscInt_FMT ".", a_destidx);
228: n = ail->array[a_destidx];
229: if (!n) ail->array[a_destidx] = ail->array[a_srcidx];
230: else {
231: do {
232: if (!n->next) {
233: n->next = ail->array[a_srcidx];
234: break;
235: }
236: n = n->next;
237: } while (1);
238: }
239: ail->array[a_srcidx] = NULL;
240: PetscFunctionReturn(PETSC_SUCCESS);
241: }
243: /* PetscCDRemoveAll
244: */
245: PetscErrorCode PetscCDRemoveAll(PetscCoarsenData *ail, PetscInt a_idx)
246: {
247: PetscCDIntNd *rem, *n1;
249: PetscFunctionBegin;
250: PetscCheck(a_idx < ail->size, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Index %" PetscInt_FMT " out of range.", a_idx);
251: rem = ail->array[a_idx];
252: ail->array[a_idx] = NULL;
253: if (!(n1 = ail->extra_nodes)) ail->extra_nodes = rem;
254: else {
255: while (n1->next) n1 = n1->next;
256: n1->next = rem;
257: }
258: PetscFunctionReturn(PETSC_SUCCESS);
259: }
261: /* PetscCDSizeAt
262: */
263: PetscErrorCode PetscCDSizeAt(const PetscCoarsenData *ail, PetscInt a_idx, PetscInt *a_sz)
264: {
265: PetscCDIntNd *n1;
266: PetscInt sz = 0;
268: PetscFunctionBegin;
269: PetscCheck(a_idx < ail->size, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Index %" PetscInt_FMT " out of range.", a_idx);
270: n1 = ail->array[a_idx];
271: while (n1) {
272: n1 = n1->next;
273: sz++;
274: }
275: *a_sz = sz;
276: PetscFunctionReturn(PETSC_SUCCESS);
277: }
279: /* PetscCDEmptyAt
280: */
281: PetscErrorCode PetscCDEmptyAt(const PetscCoarsenData *ail, PetscInt a_idx, PetscBool *a_e)
282: {
283: PetscFunctionBegin;
284: PetscCheck(a_idx < ail->size, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Index %" PetscInt_FMT " out of range.", a_idx);
285: *a_e = (PetscBool)(ail->array[a_idx] == NULL);
286: PetscFunctionReturn(PETSC_SUCCESS);
287: }
289: /* PetscCDGetMIS - used for C-F methods
290: */
291: PetscErrorCode PetscCDGetMIS(PetscCoarsenData *ail, IS *a_mis)
292: {
293: PetscCDIntNd *n;
294: PetscInt ii, kk;
295: PetscInt *permute;
297: PetscFunctionBegin;
298: for (ii = kk = 0; ii < ail->size; ii++) {
299: n = ail->array[ii];
300: if (n) kk++;
301: }
302: PetscCall(PetscMalloc1(kk, &permute));
303: for (ii = kk = 0; ii < ail->size; ii++) {
304: n = ail->array[ii];
305: if (n) permute[kk++] = ii;
306: }
307: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, kk, permute, PETSC_OWN_POINTER, a_mis));
308: PetscFunctionReturn(PETSC_SUCCESS);
309: }
311: /* PetscCDGetMat
312: */
313: PetscErrorCode PetscCDGetMat(PetscCoarsenData *ail, Mat *a_mat)
314: {
315: PetscFunctionBegin;
316: *a_mat = ail->mat;
317: ail->mat = NULL; // give it up
318: PetscFunctionReturn(PETSC_SUCCESS);
319: }
321: /* PetscCDSetMat
322: */
323: PetscErrorCode PetscCDSetMat(PetscCoarsenData *ail, Mat a_mat)
324: {
325: PetscFunctionBegin;
326: if (ail->mat) {
327: PetscCall(MatDestroy(&ail->mat)); //should not happen
328: }
329: ail->mat = a_mat;
330: PetscFunctionReturn(PETSC_SUCCESS);
331: }
333: /* PetscCDGetASMBlocks - get IS of aggregates for ASM smoothers
334: */
335: PetscErrorCode PetscCDGetASMBlocks(const PetscCoarsenData *ail, const PetscInt a_bs, Mat mat, PetscInt *a_sz, IS **a_local_is)
336: {
337: PetscCDIntNd *n;
338: PetscInt lsz, ii, kk, *idxs, jj, s, e, gid;
339: IS *is_loc, is_bcs;
341: PetscFunctionBegin;
342: for (ii = kk = 0; ii < ail->size; ii++) {
343: if (ail->array[ii]) kk++;
344: }
345: /* count BCs */
346: PetscCall(MatGetOwnershipRange(mat, &s, &e));
347: for (gid = s, lsz = 0; gid < e; gid++) {
348: PetscCall(MatGetRow(mat, gid, &jj, NULL, NULL));
349: if (jj < 2) lsz++;
350: PetscCall(MatRestoreRow(mat, gid, &jj, NULL, NULL));
351: }
352: if (lsz) {
353: PetscCall(PetscMalloc1(a_bs * lsz, &idxs));
354: for (gid = s, lsz = 0; gid < e; gid++) {
355: PetscCall(MatGetRow(mat, gid, &jj, NULL, NULL));
356: if (jj < 2) {
357: for (jj = 0; jj < a_bs; lsz++, jj++) idxs[lsz] = a_bs * gid + jj;
358: }
359: PetscCall(MatRestoreRow(mat, gid, &jj, NULL, NULL));
360: }
361: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, lsz, idxs, PETSC_OWN_POINTER, &is_bcs));
362: *a_sz = kk + 1; /* out */
363: } else {
364: is_bcs = NULL;
365: *a_sz = kk; /* out */
366: }
367: PetscCall(PetscMalloc1(*a_sz, &is_loc));
369: for (ii = kk = 0; ii < ail->size; ii++) {
370: for (lsz = 0, n = ail->array[ii]; n; lsz++, n = n->next) /* void */
371: ;
372: if (lsz) {
373: PetscCall(PetscMalloc1(a_bs * lsz, &idxs));
374: for (lsz = 0, n = ail->array[ii]; n; n = n->next) {
375: PetscCall(PetscCDIntNdGetID(n, &gid));
376: for (jj = 0; jj < a_bs; lsz++, jj++) idxs[lsz] = a_bs * gid + jj;
377: }
378: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, lsz, idxs, PETSC_OWN_POINTER, &is_loc[kk++]));
379: }
380: }
381: if (is_bcs) is_loc[kk++] = is_bcs;
382: PetscCheck(*a_sz == kk, PETSC_COMM_SELF, PETSC_ERR_PLIB, "*a_sz %" PetscInt_FMT " != kk %" PetscInt_FMT, *a_sz, kk);
383: *a_local_is = is_loc; /* out */
385: PetscFunctionReturn(PETSC_SUCCESS);
386: }
388: /* edge for priority queue */
389: typedef struct edge_tag {
390: PetscReal weight;
391: PetscInt lid0, gid1, cpid1;
392: } Edge;
394: static int gamg_hem_compare(const void *a, const void *b)
395: {
396: PetscReal va = ((Edge *)a)->weight, vb = ((Edge *)b)->weight;
397: return (va < vb) ? 1 : (va == vb) ? 0 : -1; /* 0 for equal */
398: }
400: /*
401: MatCoarsenApply_HEM_private - parallel heavy edge matching
403: Input Parameter:
404: . perm - permutation
405: . a_Gmat - global matrix of the graph
407: Output Parameter:
408: . a_locals_llist - array of list of local nodes rooted at local node
409: */
410: static PetscErrorCode MatCoarsenApply_HEM_private(IS perm, Mat a_Gmat, PetscCoarsenData **a_locals_llist)
411: {
412: PetscBool isMPI;
413: MPI_Comm comm;
414: PetscInt sub_it, kk, n, ix, *idx, *ii, iter, Iend, my0;
415: PetscMPIInt rank, size;
416: const PetscInt nloc = a_Gmat->rmap->n, n_iter = 4; /* need to figure out how to stop this */
417: PetscInt *lid_cprowID, *lid_gid;
418: PetscBool *lid_matched;
419: Mat_SeqAIJ *matA, *matB = NULL;
420: Mat_MPIAIJ *mpimat = NULL;
421: PetscScalar one = 1.;
422: PetscCoarsenData *agg_llists = NULL, *deleted_list = NULL;
423: Mat cMat, tMat, P;
424: MatScalar *ap;
425: PetscMPIInt tag1, tag2;
427: PetscFunctionBegin;
428: PetscCall(PetscObjectGetComm((PetscObject)a_Gmat, &comm));
429: PetscCallMPI(MPI_Comm_rank(comm, &rank));
430: PetscCallMPI(MPI_Comm_size(comm, &size));
431: PetscCall(MatGetOwnershipRange(a_Gmat, &my0, &Iend));
432: PetscCall(PetscCommGetNewTag(comm, &tag1));
433: PetscCall(PetscCommGetNewTag(comm, &tag2));
435: PetscCall(PetscMalloc1(nloc, &lid_gid)); /* explicit array needed */
436: PetscCall(PetscMalloc1(nloc, &lid_cprowID));
437: PetscCall(PetscMalloc1(nloc, &lid_matched));
439: PetscCall(PetscCDCreate(nloc, &agg_llists));
440: /* PetscCall(PetscCDSetChuckSize(agg_llists, nloc+1)); */
441: *a_locals_llist = agg_llists;
442: PetscCall(PetscCDCreate(size, &deleted_list));
443: PetscCall(PetscCDSetChuckSize(deleted_list, 100));
444: /* setup 'lid_gid' for scatters and add self to all lists */
445: for (kk = 0; kk < nloc; kk++) {
446: lid_gid[kk] = kk + my0;
447: PetscCall(PetscCDAppendID(agg_llists, kk, my0 + kk));
448: }
450: /* make a copy of the graph, this gets destroyed in iterates */
451: PetscCall(MatDuplicate(a_Gmat, MAT_COPY_VALUES, &cMat));
452: PetscCall(PetscObjectTypeCompare((PetscObject)a_Gmat, MATMPIAIJ, &isMPI));
453: iter = 0;
454: while (iter++ < n_iter) {
455: PetscScalar *cpcol_gid, *cpcol_max_ew, *cpcol_max_pe, *lid_max_ew;
456: PetscBool *cpcol_matched;
457: PetscMPIInt *cpcol_pe, proc;
458: Vec locMaxEdge, locMaxPE, ghostMaxEdge, ghostMaxPE;
459: PetscInt nEdges, n_nz_row, jj;
460: Edge *Edges;
461: PetscInt gid;
462: const PetscInt *perm_ix, n_sub_its = 120;
464: /* get submatrices of cMat */
465: if (isMPI) {
466: mpimat = (Mat_MPIAIJ *)cMat->data;
467: matA = (Mat_SeqAIJ *)mpimat->A->data;
468: matB = (Mat_SeqAIJ *)mpimat->B->data;
469: if (!matB->compressedrow.use) {
470: /* force construction of compressed row data structure since code below requires it */
471: PetscCall(MatCheckCompressedRow(mpimat->B, matB->nonzerorowcnt, &matB->compressedrow, matB->i, mpimat->B->rmap->n, -1.0));
472: }
473: } else {
474: matA = (Mat_SeqAIJ *)cMat->data;
475: }
477: /* set max edge on nodes */
478: PetscCall(MatCreateVecs(cMat, &locMaxEdge, NULL));
479: PetscCall(MatCreateVecs(cMat, &locMaxPE, NULL));
481: /* get 'cpcol_pe' & 'cpcol_gid' & init. 'cpcol_matched' using 'mpimat->lvec' */
482: if (mpimat) {
483: Vec vec;
484: PetscScalar vval;
486: PetscCall(MatCreateVecs(cMat, &vec, NULL));
487: /* cpcol_pe */
488: vval = (PetscScalar)(rank);
489: for (kk = 0, gid = my0; kk < nloc; kk++, gid++) { PetscCall(VecSetValues(vec, 1, &gid, &vval, INSERT_VALUES)); /* set with GID */ }
490: PetscCall(VecAssemblyBegin(vec));
491: PetscCall(VecAssemblyEnd(vec));
492: PetscCall(VecScatterBegin(mpimat->Mvctx, vec, mpimat->lvec, INSERT_VALUES, SCATTER_FORWARD));
493: PetscCall(VecScatterEnd(mpimat->Mvctx, vec, mpimat->lvec, INSERT_VALUES, SCATTER_FORWARD));
494: PetscCall(VecGetArray(mpimat->lvec, &cpcol_gid)); /* get proc ID in 'cpcol_gid' */
495: PetscCall(VecGetLocalSize(mpimat->lvec, &n));
496: PetscCall(PetscMalloc1(n, &cpcol_pe));
497: for (kk = 0; kk < n; kk++) cpcol_pe[kk] = (PetscMPIInt)PetscRealPart(cpcol_gid[kk]);
498: PetscCall(VecRestoreArray(mpimat->lvec, &cpcol_gid));
500: /* cpcol_gid */
501: for (kk = 0, gid = my0; kk < nloc; kk++, gid++) {
502: vval = (PetscScalar)(gid);
503: PetscCall(VecSetValues(vec, 1, &gid, &vval, INSERT_VALUES)); /* set with GID */
504: }
505: PetscCall(VecAssemblyBegin(vec));
506: PetscCall(VecAssemblyEnd(vec));
507: PetscCall(VecScatterBegin(mpimat->Mvctx, vec, mpimat->lvec, INSERT_VALUES, SCATTER_FORWARD));
508: PetscCall(VecScatterEnd(mpimat->Mvctx, vec, mpimat->lvec, INSERT_VALUES, SCATTER_FORWARD));
509: PetscCall(VecDestroy(&vec));
510: PetscCall(VecGetArray(mpimat->lvec, &cpcol_gid)); /* get proc ID in 'cpcol_gid' */
512: /* cpcol_matched */
513: PetscCall(VecGetLocalSize(mpimat->lvec, &n));
514: PetscCall(PetscMalloc1(n, &cpcol_matched));
515: for (kk = 0; kk < n; kk++) cpcol_matched[kk] = PETSC_FALSE;
516: }
518: /* need an inverse map - locals */
519: for (kk = 0; kk < nloc; kk++) lid_cprowID[kk] = -1;
520: /* set index into compressed row 'lid_cprowID' */
521: if (matB) {
522: for (ix = 0; ix < matB->compressedrow.nrows; ix++) lid_cprowID[matB->compressedrow.rindex[ix]] = ix;
523: }
525: /* compute 'locMaxEdge' & 'locMaxPE', and create list of edges, count edges' */
526: for (nEdges = 0, kk = 0, gid = my0; kk < nloc; kk++, gid++) {
527: PetscReal max_e = 0., tt;
528: PetscScalar vval;
529: PetscInt lid = kk;
530: PetscMPIInt max_pe = rank, pe;
532: ii = matA->i;
533: n = ii[lid + 1] - ii[lid];
534: idx = matA->j + ii[lid];
535: ap = matA->a + ii[lid];
536: for (jj = 0; jj < n; jj++) {
537: PetscInt lidj = idx[jj];
538: if (lidj != lid && PetscRealPart(ap[jj]) > max_e) max_e = PetscRealPart(ap[jj]);
539: if (lidj > lid) nEdges++;
540: }
541: if ((ix = lid_cprowID[lid]) != -1) { /* if I have any ghost neighbors */
542: ii = matB->compressedrow.i;
543: n = ii[ix + 1] - ii[ix];
544: ap = matB->a + ii[ix];
545: idx = matB->j + ii[ix];
546: for (jj = 0; jj < n; jj++) {
547: if ((tt = PetscRealPart(ap[jj])) > max_e) max_e = tt;
548: nEdges++;
549: if ((pe = cpcol_pe[idx[jj]]) > max_pe) max_pe = pe;
550: }
551: }
552: vval = max_e;
553: PetscCall(VecSetValues(locMaxEdge, 1, &gid, &vval, INSERT_VALUES));
555: vval = (PetscScalar)max_pe;
556: PetscCall(VecSetValues(locMaxPE, 1, &gid, &vval, INSERT_VALUES));
557: }
558: PetscCall(VecAssemblyBegin(locMaxEdge));
559: PetscCall(VecAssemblyEnd(locMaxEdge));
560: PetscCall(VecAssemblyBegin(locMaxPE));
561: PetscCall(VecAssemblyEnd(locMaxPE));
563: /* get 'cpcol_max_ew' & 'cpcol_max_pe' */
564: if (mpimat) {
565: PetscCall(VecDuplicate(mpimat->lvec, &ghostMaxEdge));
566: PetscCall(VecScatterBegin(mpimat->Mvctx, locMaxEdge, ghostMaxEdge, INSERT_VALUES, SCATTER_FORWARD));
567: PetscCall(VecScatterEnd(mpimat->Mvctx, locMaxEdge, ghostMaxEdge, INSERT_VALUES, SCATTER_FORWARD));
568: PetscCall(VecGetArray(ghostMaxEdge, &cpcol_max_ew));
570: PetscCall(VecDuplicate(mpimat->lvec, &ghostMaxPE));
571: PetscCall(VecScatterBegin(mpimat->Mvctx, locMaxPE, ghostMaxPE, INSERT_VALUES, SCATTER_FORWARD));
572: PetscCall(VecScatterEnd(mpimat->Mvctx, locMaxPE, ghostMaxPE, INSERT_VALUES, SCATTER_FORWARD));
573: PetscCall(VecGetArray(ghostMaxPE, &cpcol_max_pe));
574: }
576: /* setup sorted list of edges */
577: PetscCall(PetscMalloc1(nEdges, &Edges));
578: PetscCall(ISGetIndices(perm, &perm_ix));
579: for (nEdges = n_nz_row = kk = 0; kk < nloc; kk++) {
580: PetscInt nn, lid = perm_ix[kk];
581: ii = matA->i;
582: nn = n = ii[lid + 1] - ii[lid];
583: idx = matA->j + ii[lid];
584: ap = matA->a + ii[lid];
585: for (jj = 0; jj < n; jj++) {
586: PetscInt lidj = idx[jj];
587: if (lidj > lid) {
588: Edges[nEdges].lid0 = lid;
589: Edges[nEdges].gid1 = lidj + my0;
590: Edges[nEdges].cpid1 = -1;
591: Edges[nEdges].weight = PetscRealPart(ap[jj]);
592: nEdges++;
593: }
594: }
595: if ((ix = lid_cprowID[lid]) != -1) { /* if I have any ghost neighbors */
596: ii = matB->compressedrow.i;
597: n = ii[ix + 1] - ii[ix];
598: ap = matB->a + ii[ix];
599: idx = matB->j + ii[ix];
600: nn += n;
601: for (jj = 0; jj < n; jj++) {
602: Edges[nEdges].lid0 = lid;
603: Edges[nEdges].gid1 = (PetscInt)PetscRealPart(cpcol_gid[idx[jj]]);
604: Edges[nEdges].cpid1 = idx[jj];
605: Edges[nEdges].weight = PetscRealPart(ap[jj]);
606: nEdges++;
607: }
608: }
609: if (nn > 1) n_nz_row++;
610: else if (iter == 1) {
611: /* should select this because it is technically in the MIS but lets not */
612: PetscCall(PetscCDRemoveAll(agg_llists, lid));
613: }
614: }
615: PetscCall(ISRestoreIndices(perm, &perm_ix));
617: qsort(Edges, nEdges, sizeof(Edge), gamg_hem_compare);
619: /* projection matrix */
620: PetscCall(MatCreateAIJ(comm, nloc, nloc, PETSC_DETERMINE, PETSC_DETERMINE, 1, NULL, 1, NULL, &P));
622: /* clear matched flags */
623: for (kk = 0; kk < nloc; kk++) lid_matched[kk] = PETSC_FALSE;
624: /* process - communicate - process */
625: for (sub_it = 0; sub_it < n_sub_its; sub_it++) {
626: PetscInt nactive_edges;
628: PetscCall(VecGetArray(locMaxEdge, &lid_max_ew));
629: for (kk = nactive_edges = 0; kk < nEdges; kk++) {
630: /* HEM */
631: const Edge *e = &Edges[kk];
632: const PetscInt lid0 = e->lid0, gid1 = e->gid1, cpid1 = e->cpid1, gid0 = lid0 + my0, lid1 = gid1 - my0;
633: PetscBool isOK = PETSC_TRUE;
635: /* skip if either (local) vertex is done already */
636: if (lid_matched[lid0] || (gid1 >= my0 && gid1 < Iend && lid_matched[gid1 - my0])) continue;
638: /* skip if ghost vertex is done */
639: if (cpid1 != -1 && cpcol_matched[cpid1]) continue;
641: nactive_edges++;
642: /* skip if I have a bigger edge someplace (lid_max_ew gets updated) */
643: if (PetscRealPart(lid_max_ew[lid0]) > e->weight + PETSC_SMALL) continue;
645: if (cpid1 == -1) {
646: if (PetscRealPart(lid_max_ew[lid1]) > e->weight + PETSC_SMALL) continue;
647: } else {
648: /* see if edge might get matched on other proc */
649: PetscReal g_max_e = PetscRealPart(cpcol_max_ew[cpid1]);
650: if (g_max_e > e->weight + PETSC_SMALL) continue;
651: else if (e->weight > g_max_e - PETSC_SMALL && (PetscMPIInt)PetscRealPart(cpcol_max_pe[cpid1]) > rank) {
652: /* check for max_e == to this edge and larger processor that will deal with this */
653: continue;
654: }
655: }
657: /* check ghost for v0 */
658: if (isOK) {
659: PetscReal max_e, ew;
660: if ((ix = lid_cprowID[lid0]) != -1) { /* if I have any ghost neighbors */
661: ii = matB->compressedrow.i;
662: n = ii[ix + 1] - ii[ix];
663: ap = matB->a + ii[ix];
664: idx = matB->j + ii[ix];
665: for (jj = 0; jj < n && isOK; jj++) {
666: PetscInt lidj = idx[jj];
667: if (cpcol_matched[lidj]) continue;
668: ew = PetscRealPart(ap[jj]);
669: max_e = PetscRealPart(cpcol_max_ew[lidj]);
670: /* check for max_e == to this edge and larger processor that will deal with this */
671: if (ew > max_e - PETSC_SMALL && ew > PetscRealPart(lid_max_ew[lid0]) - PETSC_SMALL && (PetscMPIInt)PetscRealPart(cpcol_max_pe[lidj]) > rank) isOK = PETSC_FALSE;
672: }
673: }
675: /* for v1 */
676: if (cpid1 == -1 && isOK) {
677: if ((ix = lid_cprowID[lid1]) != -1) { /* if I have any ghost neighbors */
678: ii = matB->compressedrow.i;
679: n = ii[ix + 1] - ii[ix];
680: ap = matB->a + ii[ix];
681: idx = matB->j + ii[ix];
682: for (jj = 0; jj < n && isOK; jj++) {
683: PetscInt lidj = idx[jj];
684: if (cpcol_matched[lidj]) continue;
685: ew = PetscRealPart(ap[jj]);
686: max_e = PetscRealPart(cpcol_max_ew[lidj]);
687: /* check for max_e == to this edge and larger processor that will deal with this */
688: if (ew > max_e - PETSC_SMALL && ew > PetscRealPart(lid_max_ew[lid1]) - PETSC_SMALL && (PetscMPIInt)PetscRealPart(cpcol_max_pe[lidj]) > rank) isOK = PETSC_FALSE;
689: }
690: }
691: }
692: }
694: /* do it */
695: if (isOK) {
696: if (cpid1 == -1) {
697: lid_matched[lid1] = PETSC_TRUE; /* keep track of what we've done this round */
698: PetscCall(PetscCDAppendRemove(agg_llists, lid0, lid1));
699: } else if (sub_it != n_sub_its - 1) {
700: /* add gid1 to list of ghost deleted by me -- I need their children */
701: proc = cpcol_pe[cpid1];
703: cpcol_matched[cpid1] = PETSC_TRUE; /* messing with VecGetArray array -- needed??? */
705: PetscCall(PetscCDAppendID(deleted_list, proc, cpid1)); /* cache to send messages */
706: PetscCall(PetscCDAppendID(deleted_list, proc, lid0));
707: } else continue;
709: lid_matched[lid0] = PETSC_TRUE; /* keep track of what we've done this round */
710: /* set projection */
711: PetscCall(MatSetValues(P, 1, &gid0, 1, &gid0, &one, INSERT_VALUES));
712: PetscCall(MatSetValues(P, 1, &gid1, 1, &gid0, &one, INSERT_VALUES));
713: } /* matched */
714: } /* edge loop */
716: /* deal with deleted ghost on first pass */
717: if (size > 1 && sub_it != n_sub_its - 1) {
718: #define REQ_BF_SIZE 100
719: PetscCDIntNd *pos;
720: PetscBool ise = PETSC_FALSE;
721: PetscInt nSend1, **sbuffs1, nSend2;
722: MPI_Request *sreqs2[REQ_BF_SIZE], *rreqs2[REQ_BF_SIZE];
723: MPI_Status status;
725: /* send request */
726: for (proc = 0, nSend1 = 0; proc < size; proc++) {
727: PetscCall(PetscCDEmptyAt(deleted_list, proc, &ise));
728: if (!ise) nSend1++;
729: }
730: PetscCall(PetscMalloc1(nSend1, &sbuffs1));
731: for (proc = 0, nSend1 = 0; proc < size; proc++) {
732: /* count ghosts */
733: PetscCall(PetscCDSizeAt(deleted_list, proc, &n));
734: if (n > 0) {
735: #define CHUNCK_SIZE 100
736: PetscInt *sbuff, *pt;
737: MPI_Request *request;
738: n /= 2;
739: PetscCall(PetscMalloc1(2 + 2 * n + n * CHUNCK_SIZE + 2, &sbuff));
740: /* save requests */
741: sbuffs1[nSend1] = sbuff;
742: request = (MPI_Request *)sbuff;
743: sbuff = pt = (PetscInt *)(request + 1);
744: *pt++ = n;
745: *pt++ = rank;
747: PetscCall(PetscCDGetHeadPos(deleted_list, proc, &pos));
748: while (pos) {
749: PetscInt lid0, cpid, gid;
750: PetscCall(PetscCDIntNdGetID(pos, &cpid));
751: gid = (PetscInt)PetscRealPart(cpcol_gid[cpid]);
752: PetscCall(PetscCDGetNextPos(deleted_list, proc, &pos));
753: PetscCall(PetscCDIntNdGetID(pos, &lid0));
754: PetscCall(PetscCDGetNextPos(deleted_list, proc, &pos));
755: *pt++ = gid;
756: *pt++ = lid0;
757: }
758: /* send request tag1 [n, proc, n*[gid1,lid0] ] */
759: PetscCallMPI(MPI_Isend(sbuff, 2 * n + 2, MPIU_INT, proc, tag1, comm, request));
760: /* post receive */
761: request = (MPI_Request *)pt;
762: rreqs2[nSend1] = request; /* cache recv request */
763: pt = (PetscInt *)(request + 1);
764: PetscCallMPI(MPI_Irecv(pt, n * CHUNCK_SIZE, MPIU_INT, proc, tag2, comm, request));
765: /* clear list */
766: PetscCall(PetscCDRemoveAll(deleted_list, proc));
767: nSend1++;
768: }
769: }
770: /* receive requests, send response, clear lists */
771: kk = nactive_edges;
772: PetscCall(MPIU_Allreduce(&kk, &nactive_edges, 1, MPIU_INT, MPI_SUM, comm)); /* not correct synchronization and global */
773: nSend2 = 0;
774: while (1) {
775: #define BF_SZ 10000
776: PetscMPIInt flag, count;
777: PetscInt rbuff[BF_SZ], *pt, *pt2, *pt3, count2, *sbuff, count3;
778: MPI_Request *request;
780: PetscCallMPI(MPI_Iprobe(MPI_ANY_SOURCE, tag1, comm, &flag, &status));
781: if (!flag) break;
782: PetscCallMPI(MPI_Get_count(&status, MPIU_INT, &count));
783: PetscCheck(count <= BF_SZ, PETSC_COMM_SELF, PETSC_ERR_SUP, "buffer too small for receive: %d", count);
784: proc = status.MPI_SOURCE;
785: /* receive request tag1 [n, proc, n*[gid1,lid0] ] */
786: PetscCallMPI(MPI_Recv(rbuff, count, MPIU_INT, proc, tag1, comm, &status));
787: /* count sends */
788: pt = rbuff;
789: count3 = count2 = 0;
790: n = *pt++;
791: kk = *pt++;
792: while (n--) {
793: PetscInt gid1 = *pt++, lid1 = gid1 - my0;
794: kk = *pt++;
795: PetscCheck(!lid_matched[lid1], PETSC_COMM_SELF, PETSC_ERR_PLIB, "Received deleted gid %" PetscInt_FMT ", deleted by (lid) %" PetscInt_FMT " from proc %" PetscInt_FMT, sub_it, gid1, kk);
796: lid_matched[lid1] = PETSC_TRUE; /* keep track of what we've done this round */
797: PetscCall(PetscCDSizeAt(agg_llists, lid1, &kk));
798: count2 += kk + 2;
799: count3++; /* number of verts requested (n) */
800: }
801: PetscCheck(count2 <= count3 * CHUNCK_SIZE, PETSC_COMM_SELF, PETSC_ERR_SUP, "Irecv will be too small: %" PetscInt_FMT, count2);
802: /* send tag2 *[lid0, n, n*[gid] ] */
803: PetscCall(PetscMalloc(count2 * sizeof(PetscInt) + sizeof(MPI_Request), &sbuff));
804: request = (MPI_Request *)sbuff;
805: sreqs2[nSend2++] = request; /* cache request */
806: PetscCheck(nSend2 != REQ_BF_SIZE, PETSC_COMM_SELF, PETSC_ERR_SUP, "buffer too small for requests: %" PetscInt_FMT, nSend2);
807: pt2 = sbuff = (PetscInt *)(request + 1);
808: pt = rbuff;
809: n = *pt++;
810: kk = *pt++;
811: while (n--) {
812: /* read [n, proc, n*[gid1,lid0] */
813: PetscInt gid1 = *pt++, lid1 = gid1 - my0, lid0 = *pt++;
814: /* write [lid0, n, n*[gid] ] */
815: *pt2++ = lid0;
816: pt3 = pt2++; /* save pointer for later */
817: PetscCall(PetscCDGetHeadPos(agg_llists, lid1, &pos));
818: while (pos) {
819: PetscInt gid;
820: PetscCall(PetscCDIntNdGetID(pos, &gid));
821: PetscCall(PetscCDGetNextPos(agg_llists, lid1, &pos));
822: *pt2++ = gid;
823: }
824: *pt3 = (pt2 - pt3) - 1;
825: /* clear list */
826: PetscCall(PetscCDRemoveAll(agg_llists, lid1));
827: }
828: /* send requested data tag2 *[lid0, n, n*[gid1] ] */
829: PetscCallMPI(MPI_Isend(sbuff, count2, MPIU_INT, proc, tag2, comm, request));
830: }
832: /* receive tag2 *[lid0, n, n*[gid] ] */
833: for (kk = 0; kk < nSend1; kk++) {
834: PetscMPIInt count;
835: MPI_Request *request;
836: PetscInt *pt, *pt2;
838: request = rreqs2[kk]; /* no need to free -- buffer is in 'sbuffs1' */
839: PetscCallMPI(MPI_Wait(request, &status));
840: PetscCallMPI(MPI_Get_count(&status, MPIU_INT, &count));
841: pt = pt2 = (PetscInt *)(request + 1);
842: while (pt - pt2 < count) {
843: PetscInt lid0 = *pt++, n = *pt++;
844: while (n--) {
845: PetscInt gid1 = *pt++;
846: PetscCall(PetscCDAppendID(agg_llists, lid0, gid1));
847: }
848: }
849: }
851: /* wait for tag1 isends */
852: while (nSend1--) {
853: MPI_Request *request;
854: request = (MPI_Request *)sbuffs1[nSend1];
855: PetscCallMPI(MPI_Wait(request, &status));
856: PetscCall(PetscFree(request));
857: }
858: PetscCall(PetscFree(sbuffs1));
860: /* wait for tag2 isends */
861: while (nSend2--) {
862: MPI_Request *request = sreqs2[nSend2];
863: PetscCallMPI(MPI_Wait(request, &status));
864: PetscCall(PetscFree(request));
865: }
867: PetscCall(VecRestoreArray(ghostMaxEdge, &cpcol_max_ew));
868: PetscCall(VecRestoreArray(ghostMaxPE, &cpcol_max_pe));
870: /* get 'cpcol_matched' - use locMaxPE, ghostMaxEdge, cpcol_max_ew */
871: for (kk = 0, gid = my0; kk < nloc; kk++, gid++) {
872: PetscScalar vval = lid_matched[kk] ? 1.0 : 0.0;
873: PetscCall(VecSetValues(locMaxPE, 1, &gid, &vval, INSERT_VALUES)); /* set with GID */
874: }
875: PetscCall(VecAssemblyBegin(locMaxPE));
876: PetscCall(VecAssemblyEnd(locMaxPE));
877: PetscCall(VecScatterBegin(mpimat->Mvctx, locMaxPE, ghostMaxEdge, INSERT_VALUES, SCATTER_FORWARD));
878: PetscCall(VecScatterEnd(mpimat->Mvctx, locMaxPE, ghostMaxEdge, INSERT_VALUES, SCATTER_FORWARD));
879: PetscCall(VecGetArray(ghostMaxEdge, &cpcol_max_ew));
880: PetscCall(VecGetLocalSize(mpimat->lvec, &n));
881: for (kk = 0; kk < n; kk++) cpcol_matched[kk] = (PetscBool)(PetscRealPart(cpcol_max_ew[kk]) != 0.0);
882: PetscCall(VecRestoreArray(ghostMaxEdge, &cpcol_max_ew));
883: } /* size > 1 */
885: /* compute 'locMaxEdge' */
886: PetscCall(VecRestoreArray(locMaxEdge, &lid_max_ew));
887: for (kk = 0, gid = my0; kk < nloc; kk++, gid++) {
888: PetscReal max_e = 0., tt;
889: PetscScalar vval;
890: PetscInt lid = kk;
892: if (lid_matched[lid]) vval = 0.;
893: else {
894: ii = matA->i;
895: n = ii[lid + 1] - ii[lid];
896: idx = matA->j + ii[lid];
897: ap = matA->a + ii[lid];
898: for (jj = 0; jj < n; jj++) {
899: PetscInt lidj = idx[jj];
900: if (lid_matched[lidj]) continue; /* this is new - can change local max */
901: if (lidj != lid && PetscRealPart(ap[jj]) > max_e) max_e = PetscRealPart(ap[jj]);
902: }
903: if (lid_cprowID && (ix = lid_cprowID[lid]) != -1) { /* if I have any ghost neighbors */
904: ii = matB->compressedrow.i;
905: n = ii[ix + 1] - ii[ix];
906: ap = matB->a + ii[ix];
907: idx = matB->j + ii[ix];
908: for (jj = 0; jj < n; jj++) {
909: PetscInt lidj = idx[jj];
910: if (cpcol_matched[lidj]) continue;
911: if ((tt = PetscRealPart(ap[jj])) > max_e) max_e = tt;
912: }
913: }
914: }
915: vval = (PetscScalar)max_e;
916: PetscCall(VecSetValues(locMaxEdge, 1, &gid, &vval, INSERT_VALUES)); /* set with GID */
917: }
918: PetscCall(VecAssemblyBegin(locMaxEdge));
919: PetscCall(VecAssemblyEnd(locMaxEdge));
921: if (size > 1 && sub_it != n_sub_its - 1) {
922: /* compute 'cpcol_max_ew' */
923: PetscCall(VecScatterBegin(mpimat->Mvctx, locMaxEdge, ghostMaxEdge, INSERT_VALUES, SCATTER_FORWARD));
924: PetscCall(VecScatterEnd(mpimat->Mvctx, locMaxEdge, ghostMaxEdge, INSERT_VALUES, SCATTER_FORWARD));
925: PetscCall(VecGetArray(ghostMaxEdge, &cpcol_max_ew));
926: PetscCall(VecGetArray(locMaxEdge, &lid_max_ew));
928: /* compute 'cpcol_max_pe' */
929: for (kk = 0, gid = my0; kk < nloc; kk++, gid++) {
930: PetscInt lid = kk;
931: PetscReal ew, v1_max_e, v0_max_e = PetscRealPart(lid_max_ew[lid]);
932: PetscScalar vval;
933: PetscMPIInt max_pe = rank, pe;
935: if (lid_matched[lid]) vval = (PetscScalar)rank;
936: else if ((ix = lid_cprowID[lid]) != -1) { /* if I have any ghost neighbors */ ii = matB->compressedrow.i;
937: n = ii[ix + 1] - ii[ix];
938: ap = matB->a + ii[ix];
939: idx = matB->j + ii[ix];
940: for (jj = 0; jj < n; jj++) {
941: PetscInt lidj = idx[jj];
942: if (cpcol_matched[lidj]) continue;
943: ew = PetscRealPart(ap[jj]);
944: v1_max_e = PetscRealPart(cpcol_max_ew[lidj]);
945: /* get max pe that has a max_e == to this edge w */
946: if ((pe = cpcol_pe[idx[jj]]) > max_pe && ew > v1_max_e - PETSC_SMALL && ew > v0_max_e - PETSC_SMALL) max_pe = pe;
947: }
948: vval = (PetscScalar)max_pe;
949: }
950: PetscCall(VecSetValues(locMaxPE, 1, &gid, &vval, INSERT_VALUES));
951: }
952: PetscCall(VecAssemblyBegin(locMaxPE));
953: PetscCall(VecAssemblyEnd(locMaxPE));
955: PetscCall(VecScatterBegin(mpimat->Mvctx, locMaxPE, ghostMaxPE, INSERT_VALUES, SCATTER_FORWARD));
956: PetscCall(VecScatterEnd(mpimat->Mvctx, locMaxPE, ghostMaxPE, INSERT_VALUES, SCATTER_FORWARD));
957: PetscCall(VecGetArray(ghostMaxPE, &cpcol_max_pe));
958: PetscCall(VecRestoreArray(locMaxEdge, &lid_max_ew));
959: } /* deal with deleted ghost */
960: PetscCall(PetscInfo(a_Gmat, "\t %" PetscInt_FMT ".%" PetscInt_FMT ": %" PetscInt_FMT " active edges.\n", iter, sub_it, nactive_edges));
961: if (!nactive_edges) break;
962: } /* sub_it loop */
964: /* clean up iteration */
965: PetscCall(PetscFree(Edges));
966: if (mpimat) {
967: PetscCall(VecRestoreArray(ghostMaxEdge, &cpcol_max_ew));
968: PetscCall(VecDestroy(&ghostMaxEdge));
969: PetscCall(VecRestoreArray(ghostMaxPE, &cpcol_max_pe));
970: PetscCall(VecDestroy(&ghostMaxPE));
971: PetscCall(PetscFree(cpcol_pe));
972: PetscCall(PetscFree(cpcol_matched));
973: }
975: PetscCall(VecDestroy(&locMaxEdge));
976: PetscCall(VecDestroy(&locMaxPE));
978: if (mpimat) PetscCall(VecRestoreArray(mpimat->lvec, &cpcol_gid));
980: /* create next G if needed */
981: if (iter == n_iter) { /* hard wired test - need to look at full surrounded nodes or something */
982: PetscCall(MatDestroy(&P));
983: PetscCall(MatDestroy(&cMat));
984: break;
985: } else {
986: Vec diag;
987: /* add identity for unmatched vertices so they stay alive */
988: for (kk = 0, gid = my0; kk < nloc; kk++, gid++) {
989: if (!lid_matched[kk]) {
990: gid = kk + my0;
991: PetscCall(MatGetRow(cMat, gid, &n, NULL, NULL));
992: if (n > 1) PetscCall(MatSetValues(P, 1, &gid, 1, &gid, &one, INSERT_VALUES));
993: PetscCall(MatRestoreRow(cMat, gid, &n, NULL, NULL));
994: }
995: }
996: PetscCall(MatAssemblyBegin(P, MAT_FINAL_ASSEMBLY));
997: PetscCall(MatAssemblyEnd(P, MAT_FINAL_ASSEMBLY));
999: /* project to make new graph with colapsed edges */
1000: PetscCall(MatPtAP(cMat, P, MAT_INITIAL_MATRIX, 1.0, &tMat));
1001: PetscCall(MatDestroy(&P));
1002: PetscCall(MatDestroy(&cMat));
1003: cMat = tMat;
1004: PetscCall(MatCreateVecs(cMat, &diag, NULL));
1005: PetscCall(MatGetDiagonal(cMat, diag)); /* effectively PCJACOBI */
1006: PetscCall(VecReciprocal(diag));
1007: PetscCall(VecSqrtAbs(diag));
1008: PetscCall(MatDiagonalScale(cMat, diag, diag));
1009: PetscCall(VecDestroy(&diag));
1010: }
1011: } /* coarsen iterator */
1013: /* make fake matrix */
1014: if (size > 1) {
1015: Mat mat;
1016: PetscCDIntNd *pos;
1017: PetscInt gid, NN, MM, jj = 0, mxsz = 0;
1019: for (kk = 0; kk < nloc; kk++) {
1020: PetscCall(PetscCDSizeAt(agg_llists, kk, &jj));
1021: if (jj > mxsz) mxsz = jj;
1022: }
1023: PetscCall(MatGetSize(a_Gmat, &MM, &NN));
1024: if (mxsz > MM - nloc) mxsz = MM - nloc;
1026: PetscCall(MatCreateAIJ(comm, nloc, nloc, PETSC_DETERMINE, PETSC_DETERMINE, 0, NULL, mxsz, NULL, &mat));
1028: for (kk = 0, gid = my0; kk < nloc; kk++, gid++) {
1029: /* for (pos=PetscCDGetHeadPos(agg_llists,kk) ; pos ; pos=PetscCDGetNextPos(agg_llists,kk,pos)) { */
1030: PetscCall(PetscCDGetHeadPos(agg_llists, kk, &pos));
1031: while (pos) {
1032: PetscInt gid1;
1033: PetscCall(PetscCDIntNdGetID(pos, &gid1));
1034: PetscCall(PetscCDGetNextPos(agg_llists, kk, &pos));
1036: if (gid1 < my0 || gid1 >= my0 + nloc) PetscCall(MatSetValues(mat, 1, &gid, 1, &gid1, &one, ADD_VALUES));
1037: }
1038: }
1039: PetscCall(MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY));
1040: PetscCall(MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY));
1041: PetscCall(PetscCDSetMat(agg_llists, mat));
1042: }
1044: PetscCall(PetscFree(lid_cprowID));
1045: PetscCall(PetscFree(lid_gid));
1046: PetscCall(PetscFree(lid_matched));
1047: PetscCall(PetscCDDestroy(deleted_list));
1048: PetscFunctionReturn(PETSC_SUCCESS);
1049: }
1051: /*
1052: HEM coarsen, simple greedy.
1053: */
1054: static PetscErrorCode MatCoarsenApply_HEM(MatCoarsen coarse)
1055: {
1056: Mat mat = coarse->graph;
1058: PetscFunctionBegin;
1059: if (!coarse->perm) {
1060: IS perm;
1061: PetscInt n, m;
1063: PetscCall(MatGetLocalSize(mat, &m, &n));
1064: PetscCall(ISCreateStride(PetscObjectComm((PetscObject)mat), m, 0, 1, &perm));
1065: PetscCall(MatCoarsenApply_HEM_private(perm, mat, &coarse->agg_lists));
1066: PetscCall(ISDestroy(&perm));
1067: } else {
1068: PetscCall(MatCoarsenApply_HEM_private(coarse->perm, mat, &coarse->agg_lists));
1069: }
1070: PetscFunctionReturn(PETSC_SUCCESS);
1071: }
1073: static PetscErrorCode MatCoarsenView_HEM(MatCoarsen coarse, PetscViewer viewer)
1074: {
1075: PetscMPIInt rank;
1076: PetscBool iascii;
1078: PetscFunctionBegin;
1079: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)coarse), &rank));
1080: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
1081: if (iascii) {
1082: PetscCall(PetscViewerASCIIPushSynchronized(viewer));
1083: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, " [%d] HEM aggregator\n", rank));
1084: PetscCall(PetscViewerFlush(viewer));
1085: PetscCall(PetscViewerASCIIPopSynchronized(viewer));
1086: }
1087: PetscFunctionReturn(PETSC_SUCCESS);
1088: }
1090: /*MC
1091: MATCOARSENHEM - A coarsener that uses HEM a simple greedy coarsener
1093: Level: beginner
1095: .seealso: `MatCoarsen`, `MatCoarsenSetType()`, `MatCoarsenGetData()`, `MatCoarsenType`, `MatCoarsenCreate()`
1096: M*/
1098: PETSC_EXTERN PetscErrorCode MatCoarsenCreate_HEM(MatCoarsen coarse)
1099: {
1100: PetscFunctionBegin;
1101: coarse->ops->apply = MatCoarsenApply_HEM;
1102: coarse->ops->view = MatCoarsenView_HEM;
1103: PetscFunctionReturn(PETSC_SUCCESS);
1104: }