Actual source code: plexpreallocate.c
petsc-3.5.4 2015-05-23
1: #include <petsc-private/dmpleximpl.h> /*I "petscdmplex.h" I*/
2: #include <petsc-private/isimpl.h>
3: #include <petscsf.h>
7: PetscErrorCode DMPlexPreallocateOperator(DM dm, PetscInt bs, PetscSection section, PetscSection sectionGlobal, PetscInt dnz[], PetscInt onz[], PetscInt dnzu[], PetscInt onzu[], Mat A, PetscBool fillMatrix)
8: {
9: MPI_Comm comm;
10: MatType mtype;
11: PetscSF sf, sfDof, sfAdj;
12: PetscSection leafSectionAdj, rootSectionAdj, sectionAdj;
13: PetscInt nroots, nleaves, l, p;
14: const PetscInt *leaves;
15: const PetscSFNode *remotes;
16: PetscInt dim, pStart, pEnd, numDof, globalOffStart, globalOffEnd, numCols;
17: PetscInt *tmpAdj = NULL, *adj, *rootAdj, *cols, *remoteOffsets;
18: PetscInt adjSize;
19: PetscLayout rLayout;
20: PetscInt locRows, rStart, rEnd, r;
21: PetscMPIInt size;
22: PetscBool doCommLocal, doComm, debug = PETSC_FALSE, isSymBlock, isSymSeqBlock, isSymMPIBlock;
23: PetscErrorCode ierr;
34: PetscLogEventBegin(DMPLEX_Preallocate,dm,0,0,0);
35: PetscObjectGetComm((PetscObject)dm,&comm);
36: PetscOptionsGetBool(NULL, "-dm_view_preallocation", &debug, NULL);
37: MPI_Comm_size(comm, &size);
38: DMPlexGetDimension(dm, &dim);
39: DMGetPointSF(dm, &sf);
40: PetscSFGetGraph(sf, &nroots, NULL, NULL, NULL);
41: doCommLocal = (size > 1) && (nroots >= 0) ? PETSC_TRUE : PETSC_FALSE;
42: MPI_Allreduce(&doCommLocal, &doComm, 1, MPIU_BOOL, MPI_LAND, comm);
43: /* Create dof SF based on point SF */
44: if (debug) {
45: PetscPrintf(comm, "Input Section for Preallocation:\n");
46: PetscSectionView(section, PETSC_VIEWER_STDOUT_WORLD);
47: PetscPrintf(comm, "Input Global Section for Preallocation:\n");
48: PetscSectionView(sectionGlobal, PETSC_VIEWER_STDOUT_WORLD);
49: PetscPrintf(comm, "Input SF for Preallocation:\n");
50: PetscSFView(sf, NULL);
51: }
52: PetscSFCreateRemoteOffsets(sf, section, section, &remoteOffsets);
53: PetscSFCreateSectionSF(sf, section, remoteOffsets, section, &sfDof);
54: if (debug) {
55: PetscPrintf(comm, "Dof SF for Preallocation:\n");
56: PetscSFView(sfDof, NULL);
57: }
58: /* Create section for dof adjacency (dof ==> # adj dof) */
59: PetscSectionGetChart(section, &pStart, &pEnd);
60: PetscSectionGetStorageSize(section, &numDof);
61: PetscSectionCreate(comm, &leafSectionAdj);
62: PetscSectionSetChart(leafSectionAdj, 0, numDof);
63: PetscSectionCreate(comm, &rootSectionAdj);
64: PetscSectionSetChart(rootSectionAdj, 0, numDof);
65: /* Fill in the ghost dofs on the interface */
66: PetscSFGetGraph(sf, NULL, &nleaves, &leaves, &remotes);
68: /*
69: section - maps points to (# dofs, local dofs)
70: sectionGlobal - maps points to (# dofs, global dofs)
71: leafSectionAdj - maps unowned local dofs to # adj dofs
72: rootSectionAdj - maps owned local dofs to # adj dofs
73: adj - adj global dofs indexed by leafSectionAdj
74: rootAdj - adj global dofs indexed by rootSectionAdj
75: sf - describes shared points across procs
76: sfDof - describes shared dofs across procs
77: sfAdj - describes shared adjacent dofs across procs
78: ** The bootstrapping process involves six rounds with similar structure of visiting neighbors of each point.
79: 1. Visit unowned points on interface, count adjacencies placing in leafSectionAdj
80: Reduce those counts to rootSectionAdj (now redundantly counting some interface points)
81: 2. Visit owned points on interface, count adjacencies placing in rootSectionAdj
82: Create sfAdj connecting rootSectionAdj and leafSectionAdj
83: 3. Visit unowned points on interface, write adjacencies to adj
84: Gather adj to rootAdj (note that there is redundancy in rootAdj when multiple procs find the same adjacencies)
85: 4. Visit owned points on interface, write adjacencies to rootAdj
86: Remove redundancy in rootAdj
87: ** The last two traversals use transitive closure
88: 5. Visit all owned points in the subdomain, count dofs for each point (sectionAdj)
89: Allocate memory addressed by sectionAdj (cols)
90: 6. Visit all owned points in the subdomain, insert dof adjacencies into cols
91: ** Knowing all the column adjacencies, check ownership and sum into dnz and onz
92: */
94: for (l = 0; l < nleaves; ++l) {
95: PetscInt dof, off, d, q;
96: PetscInt p = leaves[l], numAdj = PETSC_DETERMINE;
98: if ((p < pStart) || (p >= pEnd)) continue;
99: PetscSectionGetDof(section, p, &dof);
100: PetscSectionGetOffset(section, p, &off);
101: DMPlexGetAdjacency(dm, p, &numAdj, &tmpAdj);
102: for (q = 0; q < numAdj; ++q) {
103: const PetscInt padj = tmpAdj[q];
104: PetscInt ndof, ncdof;
106: if ((padj < pStart) || (padj >= pEnd)) continue;
107: PetscSectionGetDof(section, padj, &ndof);
108: PetscSectionGetConstraintDof(section, padj, &ncdof);
109: for (d = off; d < off+dof; ++d) {
110: PetscSectionAddDof(leafSectionAdj, d, ndof-ncdof);
111: }
112: }
113: }
114: PetscSectionSetUp(leafSectionAdj);
115: if (debug) {
116: PetscPrintf(comm, "Adjacency Section for Preallocation on Leaves:\n");
117: PetscSectionView(leafSectionAdj, PETSC_VIEWER_STDOUT_WORLD);
118: }
119: /* Get maximum remote adjacency sizes for owned dofs on interface (roots) */
120: if (doComm) {
121: PetscSFReduceBegin(sfDof, MPIU_INT, leafSectionAdj->atlasDof, rootSectionAdj->atlasDof, MPI_SUM);
122: PetscSFReduceEnd(sfDof, MPIU_INT, leafSectionAdj->atlasDof, rootSectionAdj->atlasDof, MPI_SUM);
123: }
124: if (debug) {
125: PetscPrintf(comm, "Adjancency Section for Preallocation on Roots:\n");
126: PetscSectionView(rootSectionAdj, PETSC_VIEWER_STDOUT_WORLD);
127: }
128: /* Add in local adjacency sizes for owned dofs on interface (roots) */
129: for (p = pStart; p < pEnd; ++p) {
130: PetscInt numAdj = PETSC_DETERMINE, adof, dof, off, d, q;
132: PetscSectionGetDof(section, p, &dof);
133: PetscSectionGetOffset(section, p, &off);
134: if (!dof) continue;
135: PetscSectionGetDof(rootSectionAdj, off, &adof);
136: if (adof <= 0) continue;
137: DMPlexGetAdjacency(dm, p, &numAdj, &tmpAdj);
138: for (q = 0; q < numAdj; ++q) {
139: const PetscInt padj = tmpAdj[q];
140: PetscInt ndof, ncdof;
142: if ((padj < pStart) || (padj >= pEnd)) continue;
143: PetscSectionGetDof(section, padj, &ndof);
144: PetscSectionGetConstraintDof(section, padj, &ncdof);
145: for (d = off; d < off+dof; ++d) {
146: PetscSectionAddDof(rootSectionAdj, d, ndof-ncdof);
147: }
148: }
149: }
150: PetscSectionSetUp(rootSectionAdj);
151: if (debug) {
152: PetscPrintf(comm, "Adjancency Section for Preallocation on Roots after local additions:\n");
153: PetscSectionView(rootSectionAdj, PETSC_VIEWER_STDOUT_WORLD);
154: }
155: /* Create adj SF based on dof SF */
156: PetscSFCreateRemoteOffsets(sfDof, rootSectionAdj, leafSectionAdj, &remoteOffsets);
157: PetscSFCreateSectionSF(sfDof, rootSectionAdj, remoteOffsets, leafSectionAdj, &sfAdj);
158: if (debug) {
159: PetscPrintf(comm, "Adjacency SF for Preallocation:\n");
160: PetscSFView(sfAdj, NULL);
161: }
162: PetscSFDestroy(&sfDof);
163: /* Create leaf adjacency */
164: PetscSectionSetUp(leafSectionAdj);
165: PetscSectionGetStorageSize(leafSectionAdj, &adjSize);
166: PetscCalloc1(adjSize, &adj);
167: for (l = 0; l < nleaves; ++l) {
168: PetscInt dof, off, d, q;
169: PetscInt p = leaves[l], numAdj = PETSC_DETERMINE;
171: if ((p < pStart) || (p >= pEnd)) continue;
172: PetscSectionGetDof(section, p, &dof);
173: PetscSectionGetOffset(section, p, &off);
174: DMPlexGetAdjacency(dm, p, &numAdj, &tmpAdj);
175: for (d = off; d < off+dof; ++d) {
176: PetscInt aoff, i = 0;
178: PetscSectionGetOffset(leafSectionAdj, d, &aoff);
179: for (q = 0; q < numAdj; ++q) {
180: const PetscInt padj = tmpAdj[q];
181: PetscInt ndof, ncdof, ngoff, nd;
183: if ((padj < pStart) || (padj >= pEnd)) continue;
184: PetscSectionGetDof(section, padj, &ndof);
185: PetscSectionGetConstraintDof(section, padj, &ncdof);
186: PetscSectionGetOffset(sectionGlobal, padj, &ngoff);
187: for (nd = 0; nd < ndof-ncdof; ++nd) {
188: adj[aoff+i] = (ngoff < 0 ? -(ngoff+1) : ngoff) + nd;
189: ++i;
190: }
191: }
192: }
193: }
194: /* Debugging */
195: if (debug) {
196: IS tmp;
197: PetscPrintf(comm, "Leaf adjacency indices\n");
198: ISCreateGeneral(comm, adjSize, adj, PETSC_USE_POINTER, &tmp);
199: ISView(tmp, NULL);
200: ISDestroy(&tmp);
201: }
202: /* Gather adjacenct indices to root */
203: PetscSectionGetStorageSize(rootSectionAdj, &adjSize);
204: PetscMalloc1(adjSize, &rootAdj);
205: for (r = 0; r < adjSize; ++r) rootAdj[r] = -1;
206: if (doComm) {
207: PetscSFGatherBegin(sfAdj, MPIU_INT, adj, rootAdj);
208: PetscSFGatherEnd(sfAdj, MPIU_INT, adj, rootAdj);
209: }
210: PetscSFDestroy(&sfAdj);
211: PetscFree(adj);
212: /* Debugging */
213: if (debug) {
214: IS tmp;
215: PetscPrintf(comm, "Root adjacency indices after gather\n");
216: ISCreateGeneral(comm, adjSize, rootAdj, PETSC_USE_POINTER, &tmp);
217: ISView(tmp, NULL);
218: ISDestroy(&tmp);
219: }
220: /* Add in local adjacency indices for owned dofs on interface (roots) */
221: for (p = pStart; p < pEnd; ++p) {
222: PetscInt numAdj = PETSC_DETERMINE, adof, dof, off, d, q;
224: PetscSectionGetDof(section, p, &dof);
225: PetscSectionGetOffset(section, p, &off);
226: if (!dof) continue;
227: PetscSectionGetDof(rootSectionAdj, off, &adof);
228: if (adof <= 0) continue;
229: DMPlexGetAdjacency(dm, p, &numAdj, &tmpAdj);
230: for (d = off; d < off+dof; ++d) {
231: PetscInt adof, aoff, i;
233: PetscSectionGetDof(rootSectionAdj, d, &adof);
234: PetscSectionGetOffset(rootSectionAdj, d, &aoff);
235: i = adof-1;
236: for (q = 0; q < numAdj; ++q) {
237: const PetscInt padj = tmpAdj[q];
238: PetscInt ndof, ncdof, ngoff, nd;
240: if ((padj < pStart) || (padj >= pEnd)) continue;
241: PetscSectionGetDof(section, padj, &ndof);
242: PetscSectionGetConstraintDof(section, padj, &ncdof);
243: PetscSectionGetOffset(sectionGlobal, padj, &ngoff);
244: for (nd = 0; nd < ndof-ncdof; ++nd) {
245: rootAdj[aoff+i] = ngoff < 0 ? -(ngoff+1)+nd : ngoff+nd;
246: --i;
247: }
248: }
249: }
250: }
251: /* Debugging */
252: if (debug) {
253: IS tmp;
254: PetscPrintf(comm, "Root adjacency indices\n");
255: ISCreateGeneral(comm, adjSize, rootAdj, PETSC_USE_POINTER, &tmp);
256: ISView(tmp, NULL);
257: ISDestroy(&tmp);
258: }
259: /* Compress indices */
260: PetscSectionSetUp(rootSectionAdj);
261: for (p = pStart; p < pEnd; ++p) {
262: PetscInt dof, cdof, off, d;
263: PetscInt adof, aoff;
265: PetscSectionGetDof(section, p, &dof);
266: PetscSectionGetConstraintDof(section, p, &cdof);
267: PetscSectionGetOffset(section, p, &off);
268: if (!dof) continue;
269: PetscSectionGetDof(rootSectionAdj, off, &adof);
270: if (adof <= 0) continue;
271: for (d = off; d < off+dof-cdof; ++d) {
272: PetscSectionGetDof(rootSectionAdj, d, &adof);
273: PetscSectionGetOffset(rootSectionAdj, d, &aoff);
274: PetscSortRemoveDupsInt(&adof, &rootAdj[aoff]);
275: PetscSectionSetDof(rootSectionAdj, d, adof);
276: }
277: }
278: /* Debugging */
279: if (debug) {
280: IS tmp;
281: PetscPrintf(comm, "Adjancency Section for Preallocation on Roots after compression:\n");
282: PetscSectionView(rootSectionAdj, PETSC_VIEWER_STDOUT_WORLD);
283: PetscPrintf(comm, "Root adjacency indices after compression\n");
284: ISCreateGeneral(comm, adjSize, rootAdj, PETSC_USE_POINTER, &tmp);
285: ISView(tmp, NULL);
286: ISDestroy(&tmp);
287: }
288: /* Build adjacency section: Maps global indices to sets of adjacent global indices */
289: PetscSectionGetOffsetRange(sectionGlobal, &globalOffStart, &globalOffEnd);
290: PetscSectionCreate(comm, §ionAdj);
291: PetscSectionSetChart(sectionAdj, globalOffStart, globalOffEnd);
292: for (p = pStart; p < pEnd; ++p) {
293: PetscInt numAdj = PETSC_DETERMINE, dof, cdof, off, goff, d, q;
294: PetscBool found = PETSC_TRUE;
296: PetscSectionGetDof(section, p, &dof);
297: PetscSectionGetConstraintDof(section, p, &cdof);
298: PetscSectionGetOffset(section, p, &off);
299: PetscSectionGetOffset(sectionGlobal, p, &goff);
300: for (d = 0; d < dof-cdof; ++d) {
301: PetscInt ldof, rdof;
303: PetscSectionGetDof(leafSectionAdj, off+d, &ldof);
304: PetscSectionGetDof(rootSectionAdj, off+d, &rdof);
305: if (ldof > 0) {
306: /* We do not own this point */
307: } else if (rdof > 0) {
308: PetscSectionSetDof(sectionAdj, goff+d, rdof);
309: } else {
310: found = PETSC_FALSE;
311: }
312: }
313: if (found) continue;
314: PetscSectionGetDof(section, p, &dof);
315: PetscSectionGetOffset(sectionGlobal, p, &goff);
316: DMPlexGetAdjacency(dm, p, &numAdj, &tmpAdj);
317: for (q = 0; q < numAdj; ++q) {
318: const PetscInt padj = tmpAdj[q];
319: PetscInt ndof, ncdof, noff;
321: if ((padj < pStart) || (padj >= pEnd)) continue;
322: PetscSectionGetDof(section, padj, &ndof);
323: PetscSectionGetConstraintDof(section, padj, &ncdof);
324: PetscSectionGetOffset(section, padj, &noff);
325: for (d = goff; d < goff+dof-cdof; ++d) {
326: PetscSectionAddDof(sectionAdj, d, ndof-ncdof);
327: }
328: }
329: }
330: PetscSectionSetUp(sectionAdj);
331: if (debug) {
332: PetscPrintf(comm, "Adjacency Section for Preallocation:\n");
333: PetscSectionView(sectionAdj, PETSC_VIEWER_STDOUT_WORLD);
334: }
335: /* Get adjacent indices */
336: PetscSectionGetStorageSize(sectionAdj, &numCols);
337: PetscMalloc1(numCols, &cols);
338: for (p = pStart; p < pEnd; ++p) {
339: PetscInt numAdj = PETSC_DETERMINE, dof, cdof, off, goff, d, q;
340: PetscBool found = PETSC_TRUE;
342: PetscSectionGetDof(section, p, &dof);
343: PetscSectionGetConstraintDof(section, p, &cdof);
344: PetscSectionGetOffset(section, p, &off);
345: PetscSectionGetOffset(sectionGlobal, p, &goff);
346: for (d = 0; d < dof-cdof; ++d) {
347: PetscInt ldof, rdof;
349: PetscSectionGetDof(leafSectionAdj, off+d, &ldof);
350: PetscSectionGetDof(rootSectionAdj, off+d, &rdof);
351: if (ldof > 0) {
352: /* We do not own this point */
353: } else if (rdof > 0) {
354: PetscInt aoff, roff;
356: PetscSectionGetOffset(sectionAdj, goff+d, &aoff);
357: PetscSectionGetOffset(rootSectionAdj, off+d, &roff);
358: PetscMemcpy(&cols[aoff], &rootAdj[roff], rdof * sizeof(PetscInt));
359: } else {
360: found = PETSC_FALSE;
361: }
362: }
363: if (found) continue;
364: DMPlexGetAdjacency(dm, p, &numAdj, &tmpAdj);
365: for (d = goff; d < goff+dof-cdof; ++d) {
366: PetscInt adof, aoff, i = 0;
368: PetscSectionGetDof(sectionAdj, d, &adof);
369: PetscSectionGetOffset(sectionAdj, d, &aoff);
370: for (q = 0; q < numAdj; ++q) {
371: const PetscInt padj = tmpAdj[q];
372: PetscInt ndof, ncdof, ngoff, nd;
373: const PetscInt *ncind;
375: /* Adjacent points may not be in the section chart */
376: if ((padj < pStart) || (padj >= pEnd)) continue;
377: PetscSectionGetDof(section, padj, &ndof);
378: PetscSectionGetConstraintDof(section, padj, &ncdof);
379: PetscSectionGetConstraintIndices(section, padj, &ncind);
380: PetscSectionGetOffset(sectionGlobal, padj, &ngoff);
381: for (nd = 0; nd < ndof-ncdof; ++nd, ++i) {
382: cols[aoff+i] = ngoff < 0 ? -(ngoff+1)+nd : ngoff+nd;
383: }
384: }
385: if (i != adof) SETERRQ4(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid number of entries %D != %D for dof %D (point %D)", i, adof, d, p);
386: }
387: }
388: PetscSectionDestroy(&leafSectionAdj);
389: PetscSectionDestroy(&rootSectionAdj);
390: PetscFree(rootAdj);
391: PetscFree(tmpAdj);
392: /* Debugging */
393: if (debug) {
394: IS tmp;
395: PetscPrintf(comm, "Column indices\n");
396: ISCreateGeneral(comm, numCols, cols, PETSC_USE_POINTER, &tmp);
397: ISView(tmp, NULL);
398: ISDestroy(&tmp);
399: }
400: /* Create allocation vectors from adjacency graph */
401: MatGetLocalSize(A, &locRows, NULL);
402: PetscLayoutCreate(PetscObjectComm((PetscObject)A), &rLayout);
403: PetscLayoutSetLocalSize(rLayout, locRows);
404: PetscLayoutSetBlockSize(rLayout, 1);
405: PetscLayoutSetUp(rLayout);
406: PetscLayoutGetRange(rLayout, &rStart, &rEnd);
407: PetscLayoutDestroy(&rLayout);
408: /* Only loop over blocks of rows */
409: if (rStart%bs || rEnd%bs) SETERRQ3(PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Invalid layout [%d, %d) for matrix, must be divisible by block size %d", rStart, rEnd, bs);
410: for (r = rStart/bs; r < rEnd/bs; ++r) {
411: const PetscInt row = r*bs;
412: PetscInt numCols, cStart, c;
414: PetscSectionGetDof(sectionAdj, row, &numCols);
415: PetscSectionGetOffset(sectionAdj, row, &cStart);
416: for (c = cStart; c < cStart+numCols; ++c) {
417: if ((cols[c] >= rStart*bs) && (cols[c] < rEnd*bs)) {
418: ++dnz[r-rStart];
419: if (cols[c] >= row) ++dnzu[r-rStart];
420: } else {
421: ++onz[r-rStart];
422: if (cols[c] >= row) ++onzu[r-rStart];
423: }
424: }
425: }
426: if (bs > 1) {
427: for (r = 0; r < locRows/bs; ++r) {
428: dnz[r] /= bs;
429: onz[r] /= bs;
430: dnzu[r] /= bs;
431: onzu[r] /= bs;
432: }
433: }
434: /* Set matrix pattern */
435: MatXAIJSetPreallocation(A, bs, dnz, onz, dnzu, onzu);
436: MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);
437: /* Check for symmetric storage */
438: MatGetType(A, &mtype);
439: PetscStrcmp(mtype, MATSBAIJ, &isSymBlock);
440: PetscStrcmp(mtype, MATSEQSBAIJ, &isSymSeqBlock);
441: PetscStrcmp(mtype, MATMPISBAIJ, &isSymMPIBlock);
442: if (isSymBlock || isSymSeqBlock || isSymMPIBlock) {MatSetOption(A, MAT_IGNORE_LOWER_TRIANGULAR, PETSC_TRUE);}
443: /* Fill matrix with zeros */
444: if (fillMatrix) {
445: PetscScalar *values;
446: PetscInt maxRowLen = 0;
448: for (r = rStart; r < rEnd; ++r) {
449: PetscInt len;
451: PetscSectionGetDof(sectionAdj, r, &len);
452: maxRowLen = PetscMax(maxRowLen, len);
453: }
454: PetscCalloc1(maxRowLen, &values);
455: for (r = rStart; r < rEnd; ++r) {
456: PetscInt numCols, cStart;
458: PetscSectionGetDof(sectionAdj, r, &numCols);
459: PetscSectionGetOffset(sectionAdj, r, &cStart);
460: MatSetValues(A, 1, &r, numCols, &cols[cStart], values, INSERT_VALUES);
461: }
462: PetscFree(values);
463: MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
464: MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
465: }
466: PetscSectionDestroy(§ionAdj);
467: PetscFree(cols);
468: PetscLogEventEnd(DMPLEX_Preallocate,dm,0,0,0);
469: return(0);
470: }
472: #if 0
475: PetscErrorCode DMPlexPreallocateOperator_2(DM dm, PetscInt bs, PetscSection section, PetscSection sectionGlobal, PetscInt dnz[], PetscInt onz[], PetscInt dnzu[], PetscInt onzu[], Mat A, PetscBool fillMatrix)
476: {
477: PetscInt *tmpClosure,*tmpAdj,*visits;
478: PetscInt c,cStart,cEnd,pStart,pEnd;
479: PetscErrorCode ierr;
482: DMPlexGetDimension(dm, &dim);
483: DMPlexGetDepth(dm, &depth);
484: DMPlexGetMaxSizes(dm, &maxConeSize, &maxSupportSize);
486: maxClosureSize = 2*PetscMax(PetscPowInt(mesh->maxConeSize,depth+1),PetscPowInt(mesh->maxSupportSize,depth+1));
488: PetscSectionGetChart(section, &pStart, &pEnd);
489: npoints = pEnd - pStart;
491: PetscMalloc3(maxClosureSize,&tmpClosure,npoints,&lvisits,npoints,&visits);
492: PetscMemzero(lvisits,(pEnd-pStart)*sizeof(PetscInt));
493: PetscMemzero(visits,(pEnd-pStart)*sizeof(PetscInt));
494: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
495: for (c=cStart; c<cEnd; c++) {
496: PetscInt *support = tmpClosure;
497: DMPlexGetTransitiveClosure(dm, c, PETSC_FALSE, &supportSize, (PetscInt**)&support);
498: for (p=0; p<supportSize; p++) lvisits[support[p]]++;
499: }
500: PetscSFReduceBegin(sf,MPIU_INT,lvisits,visits,MPI_SUM);
501: PetscSFReduceEnd (sf,MPIU_INT,lvisits,visits,MPI_SUM);
502: PetscSFBcastBegin(sf,MPIU_INT,visits,lvisits);
503: PetscSFBcastEnd (sf,MPIU_INT,visits,lvisits);
505: PetscSFGetRanks();
508: PetscMalloc2(maxClosureSize*maxClosureSize,&cellmat,npoints,&owner);
509: for (c=cStart; c<cEnd; c++) {
510: PetscMemzero(cellmat,maxClosureSize*maxClosureSize*sizeof(PetscInt));
511: /*
512: Depth-first walk of transitive closure.
513: At each leaf frame f of transitive closure that we see, add 1/visits[f] to each pair (p,q) not marked as done in cellmat.
514: This contribution is added to dnz if owning ranks of p and q match, to onz otherwise.
515: */
516: }
518: PetscSFReduceBegin(sf,MPIU_INT,ldnz,dnz,MPI_SUM);
519: PetscSFReduceEnd (sf,MPIU_INT,lonz,onz,MPI_SUM);
520: return(0);
521: }
522: #endif