Actual source code: dapreallocate.c
petsc-3.7.7 2017-09-25
1: #include <petsc/private/dmdaimpl.h> /*I "petscdmda.h" I*/
2: #include <petsc/private/isimpl.h>
3: #include <petscsf.h>
7: static PetscErrorCode DMDAGetAdjacency_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: DMDAGetTransitiveClosure(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: DMDAGetTransitiveClosure(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: DMDARestoreTransitiveClosure(dm, star[s], (PetscBool)!useClosure, &closureSize, (PetscInt**) &closure);
27: }
28: *adjSize = numAdj;
29: return(0);
30: }
34: /*@
35: DMDASetPreallocationCenterDimension - Determine the topology used to determine adjacency
37: Input Parameters:
38: + dm - The DM object
39: - preallocCenterDim - The dimension of points which connect adjacent entries
41: Level: developer
43: Notes:
44: $ FEM: Two points p and q are adjacent if q \in closure(star(p)), preallocCenterDim = dim
45: $ FVM: Two points p and q are adjacent if q \in star(cone(p)), preallocCenterDim = dim-1
46: $ FVM++: Two points p and q are adjacent if q \in star(closure(p)), preallocCenterDim = 0
48: .seealso: DMCreateMatrix(), DMDAPreallocateOperator()
49: @*/
50: PetscErrorCode DMDASetPreallocationCenterDimension(DM dm, PetscInt preallocCenterDim)
51: {
52: DM_DA *mesh = (DM_DA*) dm->data;
56: mesh->preallocCenterDim = preallocCenterDim;
57: return(0);
58: }
62: /*@
63: DMDAGetPreallocationCenterDimension - Return the topology used to determine adjacency
65: Input Parameter:
66: . dm - The DM object
68: Output Parameter:
69: . preallocCenterDim - The dimension of points which connect adjacent entries
71: Level: developer
73: Notes:
74: $ FEM: Two points p and q are adjacent if q \in closure(star(p)), preallocCenterDim = dim
75: $ FVM: Two points p and q are adjacent if q \in star(cone(p)), preallocCenterDim = dim-1
76: $ FVM++: Two points p and q are adjacent if q \in star(closure(p)), preallocCenterDim = 0
78: .seealso: DMCreateMatrix(), DMDAPreallocateOperator(), DMDASetPreallocationCenterDimension()
79: @*/
80: PetscErrorCode DMDAGetPreallocationCenterDimension(DM dm, PetscInt *preallocCenterDim)
81: {
82: DM_DA *mesh = (DM_DA*) dm->data;
87: *preallocCenterDim = mesh->preallocCenterDim;
88: return(0);
89: }
93: PetscErrorCode DMDAPreallocateOperator(DM dm, PetscInt bs, PetscSection section, PetscSection sectionGlobal, PetscInt dnz[], PetscInt onz[], PetscInt dnzu[], PetscInt onzu[], Mat A, PetscBool fillMatrix)
94: {
95: MPI_Comm comm;
96: MatType mtype;
97: PetscSF sf, sfDof, sfAdj;
98: PetscSection leafSectionAdj, rootSectionAdj, sectionAdj;
99: PetscInt nleaves, l, p;
100: const PetscInt *leaves;
101: const PetscSFNode *remotes;
102: PetscInt dim, pStart, pEnd, numDof, globalOffStart, globalOffEnd, numCols;
103: PetscInt *tmpClosure, *tmpAdj, *adj, *rootAdj, *cols, *remoteOffsets;
104: PetscInt depth, centerDim, maxConeSize, maxSupportSize, maxClosureSize, maxAdjSize, adjSize;
105: PetscLayout rLayout;
106: PetscInt locRows, rStart, rEnd, r;
107: PetscMPIInt size;
108: PetscBool useClosure, debug = PETSC_FALSE, isSymBlock, isSymSeqBlock, isSymMPIBlock;
109: PetscErrorCode ierr;
112: PetscObjectGetComm((PetscObject)dm,&comm);
113: PetscOptionsGetBool(NULL,NULL, "-dm_view_preallocation", &debug, NULL);
114: MPI_Comm_size(comm, &size);
115: DMDAGetInfo(dm, &dim,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL);
116: depth = dim;
117: DMGetPointSF(dm, &sf);
118: /* Create dof SF based on point SF */
119: if (debug) {
120: PetscPrintf(comm, "Input Section for Preallocation:\n");
121: PetscSectionView(section, PETSC_VIEWER_STDOUT_WORLD);
122: PetscPrintf(comm, "Input Global Section for Preallocation:\n");
123: PetscSectionView(sectionGlobal, PETSC_VIEWER_STDOUT_WORLD);
124: PetscPrintf(comm, "Input SF for Preallocation:\n");
125: PetscSFView(sf, NULL);
126: }
127: PetscSFCreateRemoteOffsets(sf, section, section, &remoteOffsets);
128: PetscSFCreateSectionSF(sf, section, remoteOffsets, section, &sfDof);
129: PetscFree(remoteOffsets);
130: if (debug) {
131: PetscPrintf(comm, "Dof SF for Preallocation:\n");
132: PetscSFView(sfDof, NULL);
133: }
134: /* Create section for dof adjacency (dof ==> # adj dof) */
135: /* FEM: Two points p and q are adjacent if q \in closure(star(p)), preallocCenterDim = dim */
136: /* FVM: Two points p and q are adjacent if q \in star(cone(p)), preallocCenterDim = dim-1 */
137: /* FVM++: Two points p and q are adjacent if q \in star(closure(p)), preallocCenterDim = 0 */
138: DMDAGetPreallocationCenterDimension(dm, ¢erDim);
139: if (centerDim == dim) {
140: useClosure = PETSC_FALSE;
141: } else if (centerDim == 0) {
142: useClosure = PETSC_TRUE;
143: } else SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "Do not support preallocation with center points of dimension %d", centerDim);
145: PetscSectionGetChart(section, &pStart, &pEnd);
146: PetscSectionGetStorageSize(section, &numDof);
147: PetscSectionCreate(comm, &leafSectionAdj);
148: PetscSectionSetChart(leafSectionAdj, 0, numDof);
149: PetscSectionCreate(comm, &rootSectionAdj);
150: PetscSectionSetChart(rootSectionAdj, 0, numDof);
151: /* Fill in the ghost dofs on the interface */
152: PetscSFGetGraph(sf, NULL, &nleaves, &leaves, &remotes);
153: maxConeSize = 6;
154: maxSupportSize = 6;
156: maxClosureSize = 2*PetscMax(PetscPowInt(maxConeSize,depth+1),PetscPowInt(maxSupportSize,depth+1));
157: maxAdjSize = PetscPowInt(maxConeSize,depth+1) * PetscPowInt(maxSupportSize,depth+1) + 1;
159: PetscMalloc2(maxClosureSize,&tmpClosure,maxAdjSize,&tmpAdj);
161: /*
162: ** The bootstrapping process involves six rounds with similar structure of visiting neighbors of each point.
163: 1. Visit unowned points on interface, count adjacencies placing in leafSectionAdj
164: Reduce those counts to rootSectionAdj (now redundantly counting some interface points)
165: 2. Visit owned points on interface, count adjacencies placing in rootSectionAdj
166: Create sfAdj connecting rootSectionAdj and leafSectionAdj
167: 3. Visit unowned points on interface, write adjacencies to adj
168: Gather adj to rootAdj (note that there is redundancy in rootAdj when multiple procs find the same adjacencies)
169: 4. Visit owned points on interface, write adjacencies to rootAdj
170: Remove redundancy in rootAdj
171: ** The last two traversals use transitive closure
172: 5. Visit all owned points in the subdomain, count dofs for each point (sectionAdj)
173: Allocate memory addressed by sectionAdj (cols)
174: 6. Visit all owned points in the subdomain, insert dof adjacencies into cols
175: ** Knowing all the column adjacencies, check ownership and sum into dnz and onz
176: */
178: for (l = 0; l < nleaves; ++l) {
179: PetscInt dof, off, d, q;
180: PetscInt p = leaves[l], numAdj = maxAdjSize;
182: if ((p < pStart) || (p >= pEnd)) continue;
183: PetscSectionGetDof(section, p, &dof);
184: PetscSectionGetOffset(section, p, &off);
185: DMDAGetAdjacency_Internal(dm, p, useClosure, tmpClosure, &numAdj, tmpAdj);
186: for (q = 0; q < numAdj; ++q) {
187: const PetscInt padj = tmpAdj[q];
188: PetscInt ndof, ncdof;
190: if ((padj < pStart) || (padj >= pEnd)) continue;
191: PetscSectionGetDof(section, padj, &ndof);
192: PetscSectionGetConstraintDof(section, padj, &ncdof);
193: for (d = off; d < off+dof; ++d) {
194: PetscSectionAddDof(leafSectionAdj, d, ndof-ncdof);
195: }
196: }
197: }
198: PetscSectionSetUp(leafSectionAdj);
199: if (debug) {
200: PetscPrintf(comm, "Adjacency Section for Preallocation on Leaves:\n");
201: PetscSectionView(leafSectionAdj, PETSC_VIEWER_STDOUT_WORLD);
202: }
203: /* Get maximum remote adjacency sizes for owned dofs on interface (roots) */
204: if (size > 1) {
205: PetscSFReduceBegin(sfDof, MPIU_INT, leafSectionAdj->atlasDof, rootSectionAdj->atlasDof, MPI_SUM);
206: PetscSFReduceEnd(sfDof, MPIU_INT, leafSectionAdj->atlasDof, rootSectionAdj->atlasDof, MPI_SUM);
207: }
208: if (debug) {
209: PetscPrintf(comm, "Adjancency Section for Preallocation on Roots:\n");
210: PetscSectionView(rootSectionAdj, PETSC_VIEWER_STDOUT_WORLD);
211: }
212: /* Add in local adjacency sizes for owned dofs on interface (roots) */
213: for (p = pStart; p < pEnd; ++p) {
214: PetscInt numAdj = maxAdjSize, adof, dof, off, d, q;
216: PetscSectionGetDof(section, p, &dof);
217: PetscSectionGetOffset(section, p, &off);
218: if (!dof) continue;
219: PetscSectionGetDof(rootSectionAdj, off, &adof);
220: if (adof <= 0) continue;
221: DMDAGetAdjacency_Internal(dm, p, useClosure, tmpClosure, &numAdj, tmpAdj);
222: for (q = 0; q < numAdj; ++q) {
223: const PetscInt padj = tmpAdj[q];
224: PetscInt ndof, ncdof;
226: if ((padj < pStart) || (padj >= pEnd)) continue;
227: PetscSectionGetDof(section, padj, &ndof);
228: PetscSectionGetConstraintDof(section, padj, &ncdof);
229: for (d = off; d < off+dof; ++d) {
230: PetscSectionAddDof(rootSectionAdj, d, ndof-ncdof);
231: }
232: }
233: }
234: PetscSectionSetUp(rootSectionAdj);
235: if (debug) {
236: PetscPrintf(comm, "Adjancency Section for Preallocation on Roots after local additions:\n");
237: PetscSectionView(rootSectionAdj, PETSC_VIEWER_STDOUT_WORLD);
238: }
239: /* Create adj SF based on dof SF */
240: PetscSFCreateRemoteOffsets(sfDof, rootSectionAdj, leafSectionAdj, &remoteOffsets);
241: PetscSFCreateSectionSF(sfDof, rootSectionAdj, remoteOffsets, leafSectionAdj, &sfAdj);
242: PetscFree(remoteOffsets);
243: if (debug) {
244: PetscPrintf(comm, "Adjacency SF for Preallocation:\n");
245: PetscSFView(sfAdj, NULL);
246: }
247: PetscSFDestroy(&sfDof);
248: /* Create leaf adjacency */
249: PetscSectionSetUp(leafSectionAdj);
250: PetscSectionGetStorageSize(leafSectionAdj, &adjSize);
251: PetscMalloc1(adjSize, &adj);
252: PetscMemzero(adj, adjSize * sizeof(PetscInt));
253: for (l = 0; l < nleaves; ++l) {
254: PetscInt dof, off, d, q;
255: PetscInt p = leaves[l], numAdj = maxAdjSize;
257: if ((p < pStart) || (p >= pEnd)) continue;
258: PetscSectionGetDof(section, p, &dof);
259: PetscSectionGetOffset(section, p, &off);
260: DMDAGetAdjacency_Internal(dm, p, useClosure, tmpClosure, &numAdj, tmpAdj);
261: for (d = off; d < off+dof; ++d) {
262: PetscInt aoff, i = 0;
264: PetscSectionGetOffset(leafSectionAdj, d, &aoff);
265: for (q = 0; q < numAdj; ++q) {
266: const PetscInt padj = tmpAdj[q];
267: PetscInt ndof, ncdof, ngoff, nd;
269: if ((padj < pStart) || (padj >= pEnd)) continue;
270: PetscSectionGetDof(section, padj, &ndof);
271: PetscSectionGetConstraintDof(section, padj, &ncdof);
272: PetscSectionGetOffset(sectionGlobal, padj, &ngoff);
273: for (nd = 0; nd < ndof-ncdof; ++nd) {
274: adj[aoff+i] = (ngoff < 0 ? -(ngoff+1) : ngoff) + nd;
275: ++i;
276: }
277: }
278: }
279: }
280: /* Debugging */
281: if (debug) {
282: IS tmp;
283: PetscPrintf(comm, "Leaf adjacency indices\n");
284: ISCreateGeneral(comm, adjSize, adj, PETSC_USE_POINTER, &tmp);
285: ISView(tmp, NULL);
286: ISDestroy(&tmp);
287: }
288: /* Gather adjacenct indices to root */
289: PetscSectionGetStorageSize(rootSectionAdj, &adjSize);
290: PetscMalloc1(adjSize, &rootAdj);
291: for (r = 0; r < adjSize; ++r) rootAdj[r] = -1;
292: if (size > 1) {
293: PetscSFGatherBegin(sfAdj, MPIU_INT, adj, rootAdj);
294: PetscSFGatherEnd(sfAdj, MPIU_INT, adj, rootAdj);
295: }
296: PetscSFDestroy(&sfAdj);
297: PetscFree(adj);
298: /* Debugging */
299: if (debug) {
300: IS tmp;
301: PetscPrintf(comm, "Root adjacency indices after gather\n");
302: ISCreateGeneral(comm, adjSize, rootAdj, PETSC_USE_POINTER, &tmp);
303: ISView(tmp, NULL);
304: ISDestroy(&tmp);
305: }
306: /* Add in local adjacency indices for owned dofs on interface (roots) */
307: for (p = pStart; p < pEnd; ++p) {
308: PetscInt numAdj = maxAdjSize, adof, dof, off, d, q;
310: PetscSectionGetDof(section, p, &dof);
311: PetscSectionGetOffset(section, p, &off);
312: if (!dof) continue;
313: PetscSectionGetDof(rootSectionAdj, off, &adof);
314: if (adof <= 0) continue;
315: DMDAGetAdjacency_Internal(dm, p, useClosure, tmpClosure, &numAdj, tmpAdj);
316: for (d = off; d < off+dof; ++d) {
317: PetscInt adof, aoff, i;
319: PetscSectionGetDof(rootSectionAdj, d, &adof);
320: PetscSectionGetOffset(rootSectionAdj, d, &aoff);
321: i = adof-1;
322: for (q = 0; q < numAdj; ++q) {
323: const PetscInt padj = tmpAdj[q];
324: PetscInt ndof, ncdof, ngoff, nd;
326: if ((padj < pStart) || (padj >= pEnd)) continue;
327: PetscSectionGetDof(section, padj, &ndof);
328: PetscSectionGetConstraintDof(section, padj, &ncdof);
329: PetscSectionGetOffset(sectionGlobal, padj, &ngoff);
330: for (nd = 0; nd < ndof-ncdof; ++nd) {
331: rootAdj[aoff+i] = ngoff < 0 ? -(ngoff+1)+nd : ngoff+nd;
332: --i;
333: }
334: }
335: }
336: }
337: /* Debugging */
338: if (debug) {
339: IS tmp;
340: PetscPrintf(comm, "Root adjacency indices\n");
341: ISCreateGeneral(comm, adjSize, rootAdj, PETSC_USE_POINTER, &tmp);
342: ISView(tmp, NULL);
343: ISDestroy(&tmp);
344: }
345: /* Compress indices */
346: PetscSectionSetUp(rootSectionAdj);
347: for (p = pStart; p < pEnd; ++p) {
348: PetscInt dof, cdof, off, d;
349: PetscInt adof, aoff;
351: PetscSectionGetDof(section, p, &dof);
352: PetscSectionGetConstraintDof(section, p, &cdof);
353: PetscSectionGetOffset(section, p, &off);
354: if (!dof) continue;
355: PetscSectionGetDof(rootSectionAdj, off, &adof);
356: if (adof <= 0) continue;
357: for (d = off; d < off+dof-cdof; ++d) {
358: PetscSectionGetDof(rootSectionAdj, d, &adof);
359: PetscSectionGetOffset(rootSectionAdj, d, &aoff);
360: PetscSortRemoveDupsInt(&adof, &rootAdj[aoff]);
361: PetscSectionSetDof(rootSectionAdj, d, adof);
362: }
363: }
364: /* Debugging */
365: if (debug) {
366: IS tmp;
367: PetscPrintf(comm, "Adjancency Section for Preallocation on Roots after compression:\n");
368: PetscSectionView(rootSectionAdj, PETSC_VIEWER_STDOUT_WORLD);
369: PetscPrintf(comm, "Root adjacency indices after compression\n");
370: ISCreateGeneral(comm, adjSize, rootAdj, PETSC_USE_POINTER, &tmp);
371: ISView(tmp, NULL);
372: ISDestroy(&tmp);
373: }
374: /* Build adjacency section: Maps global indices to sets of adjacent global indices */
375: PetscSectionGetOffsetRange(sectionGlobal, &globalOffStart, &globalOffEnd);
376: PetscSectionCreate(comm, §ionAdj);
377: PetscSectionSetChart(sectionAdj, globalOffStart, globalOffEnd);
378: for (p = pStart; p < pEnd; ++p) {
379: PetscInt numAdj = maxAdjSize, dof, cdof, off, goff, d, q;
380: PetscBool found = PETSC_TRUE;
382: PetscSectionGetDof(section, p, &dof);
383: PetscSectionGetConstraintDof(section, p, &cdof);
384: PetscSectionGetOffset(section, p, &off);
385: PetscSectionGetOffset(sectionGlobal, p, &goff);
386: for (d = 0; d < dof-cdof; ++d) {
387: PetscInt ldof, rdof;
389: PetscSectionGetDof(leafSectionAdj, off+d, &ldof);
390: PetscSectionGetDof(rootSectionAdj, off+d, &rdof);
391: if (ldof > 0) {
392: /* We do not own this point */
393: } else if (rdof > 0) {
394: PetscSectionSetDof(sectionAdj, goff+d, rdof);
395: } else {
396: found = PETSC_FALSE;
397: }
398: }
399: if (found) continue;
400: PetscSectionGetDof(section, p, &dof);
401: PetscSectionGetOffset(sectionGlobal, p, &goff);
402: DMDAGetAdjacency_Internal(dm, p, useClosure, tmpClosure, &numAdj, tmpAdj);
403: for (q = 0; q < numAdj; ++q) {
404: const PetscInt padj = tmpAdj[q];
405: PetscInt ndof, ncdof, noff;
407: if ((padj < pStart) || (padj >= pEnd)) continue;
408: PetscSectionGetDof(section, padj, &ndof);
409: PetscSectionGetConstraintDof(section, padj, &ncdof);
410: PetscSectionGetOffset(section, padj, &noff);
411: for (d = goff; d < goff+dof-cdof; ++d) {
412: PetscSectionAddDof(sectionAdj, d, ndof-ncdof);
413: }
414: }
415: }
416: PetscSectionSetUp(sectionAdj);
417: if (debug) {
418: PetscPrintf(comm, "Adjacency Section for Preallocation:\n");
419: PetscSectionView(sectionAdj, PETSC_VIEWER_STDOUT_WORLD);
420: }
421: /* Get adjacent indices */
422: PetscSectionGetStorageSize(sectionAdj, &numCols);
423: PetscMalloc1(numCols, &cols);
424: for (p = pStart; p < pEnd; ++p) {
425: PetscInt numAdj = maxAdjSize, dof, cdof, off, goff, d, q;
426: PetscBool found = PETSC_TRUE;
428: PetscSectionGetDof(section, p, &dof);
429: PetscSectionGetConstraintDof(section, p, &cdof);
430: PetscSectionGetOffset(section, p, &off);
431: PetscSectionGetOffset(sectionGlobal, p, &goff);
432: for (d = 0; d < dof-cdof; ++d) {
433: PetscInt ldof, rdof;
435: PetscSectionGetDof(leafSectionAdj, off+d, &ldof);
436: PetscSectionGetDof(rootSectionAdj, off+d, &rdof);
437: if (ldof > 0) {
438: /* We do not own this point */
439: } else if (rdof > 0) {
440: PetscInt aoff, roff;
442: PetscSectionGetOffset(sectionAdj, goff+d, &aoff);
443: PetscSectionGetOffset(rootSectionAdj, off+d, &roff);
444: PetscMemcpy(&cols[aoff], &rootAdj[roff], rdof * sizeof(PetscInt));
445: } else {
446: found = PETSC_FALSE;
447: }
448: }
449: if (found) continue;
450: DMDAGetAdjacency_Internal(dm, p, useClosure, tmpClosure, &numAdj, tmpAdj);
451: for (d = goff; d < goff+dof-cdof; ++d) {
452: PetscInt adof, aoff, i = 0;
454: PetscSectionGetDof(sectionAdj, d, &adof);
455: PetscSectionGetOffset(sectionAdj, d, &aoff);
456: for (q = 0; q < numAdj; ++q) {
457: const PetscInt padj = tmpAdj[q];
458: PetscInt ndof, ncdof, ngoff, nd;
459: const PetscInt *ncind;
461: /* Adjacent points may not be in the section chart */
462: if ((padj < pStart) || (padj >= pEnd)) continue;
463: PetscSectionGetDof(section, padj, &ndof);
464: PetscSectionGetConstraintDof(section, padj, &ncdof);
465: PetscSectionGetConstraintIndices(section, padj, &ncind);
466: PetscSectionGetOffset(sectionGlobal, padj, &ngoff);
467: for (nd = 0; nd < ndof-ncdof; ++nd, ++i) {
468: cols[aoff+i] = ngoff < 0 ? -(ngoff+1)+nd : ngoff+nd;
469: }
470: }
471: 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);
472: }
473: }
474: PetscSectionDestroy(&leafSectionAdj);
475: PetscSectionDestroy(&rootSectionAdj);
476: PetscFree(rootAdj);
477: PetscFree2(tmpClosure, tmpAdj);
478: /* Debugging */
479: if (debug) {
480: IS tmp;
481: PetscPrintf(comm, "Column indices\n");
482: ISCreateGeneral(comm, numCols, cols, PETSC_USE_POINTER, &tmp);
483: ISView(tmp, NULL);
484: ISDestroy(&tmp);
485: }
486: /* Create allocation vectors from adjacency graph */
487: MatGetLocalSize(A, &locRows, NULL);
488: PetscLayoutCreate(PetscObjectComm((PetscObject)A), &rLayout);
489: PetscLayoutSetLocalSize(rLayout, locRows);
490: PetscLayoutSetBlockSize(rLayout, 1);
491: PetscLayoutSetUp(rLayout);
492: PetscLayoutGetRange(rLayout, &rStart, &rEnd);
493: PetscLayoutDestroy(&rLayout);
494: /* Only loop over blocks of rows */
495: 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);
496: for (r = rStart/bs; r < rEnd/bs; ++r) {
497: const PetscInt row = r*bs;
498: PetscInt numCols, cStart, c;
500: PetscSectionGetDof(sectionAdj, row, &numCols);
501: PetscSectionGetOffset(sectionAdj, row, &cStart);
502: for (c = cStart; c < cStart+numCols; ++c) {
503: if ((cols[c] >= rStart*bs) && (cols[c] < rEnd*bs)) {
504: ++dnz[r-rStart];
505: if (cols[c] >= row) ++dnzu[r-rStart];
506: } else {
507: ++onz[r-rStart];
508: if (cols[c] >= row) ++onzu[r-rStart];
509: }
510: }
511: }
512: if (bs > 1) {
513: for (r = 0; r < locRows/bs; ++r) {
514: dnz[r] /= bs;
515: onz[r] /= bs;
516: dnzu[r] /= bs;
517: onzu[r] /= bs;
518: }
519: }
520: /* Set matrix pattern */
521: MatXAIJSetPreallocation(A, bs, dnz, onz, dnzu, onzu);
522: MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);
523: /* Check for symmetric storage */
524: MatGetType(A, &mtype);
525: PetscStrcmp(mtype, MATSBAIJ, &isSymBlock);
526: PetscStrcmp(mtype, MATSEQSBAIJ, &isSymSeqBlock);
527: PetscStrcmp(mtype, MATMPISBAIJ, &isSymMPIBlock);
528: if (isSymBlock || isSymSeqBlock || isSymMPIBlock) {MatSetOption(A, MAT_IGNORE_LOWER_TRIANGULAR, PETSC_TRUE);}
529: /* Fill matrix with zeros */
530: if (fillMatrix) {
531: PetscScalar *values;
532: PetscInt maxRowLen = 0;
534: for (r = rStart; r < rEnd; ++r) {
535: PetscInt len;
537: PetscSectionGetDof(sectionAdj, r, &len);
538: maxRowLen = PetscMax(maxRowLen, len);
539: }
540: PetscCalloc1(maxRowLen, &values);
541: for (r = rStart; r < rEnd; ++r) {
542: PetscInt numCols, cStart;
544: PetscSectionGetDof(sectionAdj, r, &numCols);
545: PetscSectionGetOffset(sectionAdj, r, &cStart);
546: MatSetValues(A, 1, &r, numCols, &cols[cStart], values, INSERT_VALUES);
547: }
548: PetscFree(values);
549: MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
550: MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
551: }
552: PetscSectionDestroy(§ionAdj);
553: PetscFree(cols);
554: return(0);
555: }