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