Actual source code: dapreallocate.c
petsc-3.6.1 2015-08-06
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, "-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: if (debug) {
130: PetscPrintf(comm, "Dof SF for Preallocation:\n");
131: PetscSFView(sfDof, NULL);
132: }
133: /* Create section for dof adjacency (dof ==> # adj dof) */
134: /* FEM: Two points p and q are adjacent if q \in closure(star(p)), preallocCenterDim = dim */
135: /* FVM: Two points p and q are adjacent if q \in star(cone(p)), preallocCenterDim = dim-1 */
136: /* FVM++: Two points p and q are adjacent if q \in star(closure(p)), preallocCenterDim = 0 */
137: DMDAGetPreallocationCenterDimension(dm, ¢erDim);
138: if (centerDim == dim) {
139: useClosure = PETSC_FALSE;
140: } else if (centerDim == 0) {
141: useClosure = PETSC_TRUE;
142: } else SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "Do not support preallocation with center points of dimension %d", centerDim);
144: PetscSectionGetChart(section, &pStart, &pEnd);
145: PetscSectionGetStorageSize(section, &numDof);
146: PetscSectionCreate(comm, &leafSectionAdj);
147: PetscSectionSetChart(leafSectionAdj, 0, numDof);
148: PetscSectionCreate(comm, &rootSectionAdj);
149: PetscSectionSetChart(rootSectionAdj, 0, numDof);
150: /* Fill in the ghost dofs on the interface */
151: PetscSFGetGraph(sf, NULL, &nleaves, &leaves, &remotes);
152: maxConeSize = 6;
153: maxSupportSize = 6;
155: maxClosureSize = 2*PetscMax(PetscPowInt(maxConeSize,depth+1),PetscPowInt(maxSupportSize,depth+1));
156: maxAdjSize = PetscPowInt(maxConeSize,depth+1) * PetscPowInt(maxSupportSize,depth+1) + 1;
158: PetscMalloc2(maxClosureSize,&tmpClosure,maxAdjSize,&tmpAdj);
160: /*
161: ** The bootstrapping process involves six rounds with similar structure of visiting neighbors of each point.
162: 1. Visit unowned points on interface, count adjacencies placing in leafSectionAdj
163: Reduce those counts to rootSectionAdj (now redundantly counting some interface points)
164: 2. Visit owned points on interface, count adjacencies placing in rootSectionAdj
165: Create sfAdj connecting rootSectionAdj and leafSectionAdj
166: 3. Visit unowned points on interface, write adjacencies to adj
167: Gather adj to rootAdj (note that there is redundancy in rootAdj when multiple procs find the same adjacencies)
168: 4. Visit owned points on interface, write adjacencies to rootAdj
169: Remove redundancy in rootAdj
170: ** The last two traversals use transitive closure
171: 5. Visit all owned points in the subdomain, count dofs for each point (sectionAdj)
172: Allocate memory addressed by sectionAdj (cols)
173: 6. Visit all owned points in the subdomain, insert dof adjacencies into cols
174: ** Knowing all the column adjacencies, check ownership and sum into dnz and onz
175: */
177: for (l = 0; l < nleaves; ++l) {
178: PetscInt dof, off, d, q;
179: PetscInt p = leaves[l], numAdj = maxAdjSize;
181: if ((p < pStart) || (p >= pEnd)) continue;
182: PetscSectionGetDof(section, p, &dof);
183: PetscSectionGetOffset(section, p, &off);
184: DMDAGetAdjacency_Internal(dm, p, useClosure, tmpClosure, &numAdj, tmpAdj);
185: for (q = 0; q < numAdj; ++q) {
186: const PetscInt padj = tmpAdj[q];
187: PetscInt ndof, ncdof;
189: if ((padj < pStart) || (padj >= pEnd)) continue;
190: PetscSectionGetDof(section, padj, &ndof);
191: PetscSectionGetConstraintDof(section, padj, &ncdof);
192: for (d = off; d < off+dof; ++d) {
193: PetscSectionAddDof(leafSectionAdj, d, ndof-ncdof);
194: }
195: }
196: }
197: PetscSectionSetUp(leafSectionAdj);
198: if (debug) {
199: PetscPrintf(comm, "Adjacency Section for Preallocation on Leaves:\n");
200: PetscSectionView(leafSectionAdj, PETSC_VIEWER_STDOUT_WORLD);
201: }
202: /* Get maximum remote adjacency sizes for owned dofs on interface (roots) */
203: if (size > 1) {
204: PetscSFReduceBegin(sfDof, MPIU_INT, leafSectionAdj->atlasDof, rootSectionAdj->atlasDof, MPI_SUM);
205: PetscSFReduceEnd(sfDof, MPIU_INT, leafSectionAdj->atlasDof, rootSectionAdj->atlasDof, MPI_SUM);
206: }
207: if (debug) {
208: PetscPrintf(comm, "Adjancency Section for Preallocation on Roots:\n");
209: PetscSectionView(rootSectionAdj, PETSC_VIEWER_STDOUT_WORLD);
210: }
211: /* Add in local adjacency sizes for owned dofs on interface (roots) */
212: for (p = pStart; p < pEnd; ++p) {
213: PetscInt numAdj = maxAdjSize, adof, dof, off, d, q;
215: PetscSectionGetDof(section, p, &dof);
216: PetscSectionGetOffset(section, p, &off);
217: if (!dof) continue;
218: PetscSectionGetDof(rootSectionAdj, off, &adof);
219: if (adof <= 0) continue;
220: DMDAGetAdjacency_Internal(dm, p, useClosure, tmpClosure, &numAdj, tmpAdj);
221: for (q = 0; q < numAdj; ++q) {
222: const PetscInt padj = tmpAdj[q];
223: PetscInt ndof, ncdof;
225: if ((padj < pStart) || (padj >= pEnd)) continue;
226: PetscSectionGetDof(section, padj, &ndof);
227: PetscSectionGetConstraintDof(section, padj, &ncdof);
228: for (d = off; d < off+dof; ++d) {
229: PetscSectionAddDof(rootSectionAdj, d, ndof-ncdof);
230: }
231: }
232: }
233: PetscSectionSetUp(rootSectionAdj);
234: if (debug) {
235: PetscPrintf(comm, "Adjancency Section for Preallocation on Roots after local additions:\n");
236: PetscSectionView(rootSectionAdj, PETSC_VIEWER_STDOUT_WORLD);
237: }
238: /* Create adj SF based on dof SF */
239: PetscSFCreateRemoteOffsets(sfDof, rootSectionAdj, leafSectionAdj, &remoteOffsets);
240: PetscSFCreateSectionSF(sfDof, rootSectionAdj, remoteOffsets, leafSectionAdj, &sfAdj);
241: if (debug) {
242: PetscPrintf(comm, "Adjacency SF for Preallocation:\n");
243: PetscSFView(sfAdj, NULL);
244: }
245: PetscSFDestroy(&sfDof);
246: /* Create leaf adjacency */
247: PetscSectionSetUp(leafSectionAdj);
248: PetscSectionGetStorageSize(leafSectionAdj, &adjSize);
249: PetscMalloc1(adjSize, &adj);
250: PetscMemzero(adj, adjSize * sizeof(PetscInt));
251: for (l = 0; l < nleaves; ++l) {
252: PetscInt dof, off, d, q;
253: PetscInt p = leaves[l], numAdj = maxAdjSize;
255: if ((p < pStart) || (p >= pEnd)) continue;
256: PetscSectionGetDof(section, p, &dof);
257: PetscSectionGetOffset(section, p, &off);
258: DMDAGetAdjacency_Internal(dm, p, useClosure, tmpClosure, &numAdj, tmpAdj);
259: for (d = off; d < off+dof; ++d) {
260: PetscInt aoff, i = 0;
262: PetscSectionGetOffset(leafSectionAdj, d, &aoff);
263: for (q = 0; q < numAdj; ++q) {
264: const PetscInt padj = tmpAdj[q];
265: PetscInt ndof, ncdof, ngoff, nd;
267: if ((padj < pStart) || (padj >= pEnd)) continue;
268: PetscSectionGetDof(section, padj, &ndof);
269: PetscSectionGetConstraintDof(section, padj, &ncdof);
270: PetscSectionGetOffset(sectionGlobal, padj, &ngoff);
271: for (nd = 0; nd < ndof-ncdof; ++nd) {
272: adj[aoff+i] = (ngoff < 0 ? -(ngoff+1) : ngoff) + nd;
273: ++i;
274: }
275: }
276: }
277: }
278: /* Debugging */
279: if (debug) {
280: IS tmp;
281: PetscPrintf(comm, "Leaf adjacency indices\n");
282: ISCreateGeneral(comm, adjSize, adj, PETSC_USE_POINTER, &tmp);
283: ISView(tmp, NULL);
284: ISDestroy(&tmp);
285: }
286: /* Gather adjacenct indices to root */
287: PetscSectionGetStorageSize(rootSectionAdj, &adjSize);
288: PetscMalloc1(adjSize, &rootAdj);
289: for (r = 0; r < adjSize; ++r) rootAdj[r] = -1;
290: if (size > 1) {
291: PetscSFGatherBegin(sfAdj, MPIU_INT, adj, rootAdj);
292: PetscSFGatherEnd(sfAdj, MPIU_INT, adj, rootAdj);
293: }
294: PetscSFDestroy(&sfAdj);
295: PetscFree(adj);
296: /* Debugging */
297: if (debug) {
298: IS tmp;
299: PetscPrintf(comm, "Root adjacency indices after gather\n");
300: ISCreateGeneral(comm, adjSize, rootAdj, PETSC_USE_POINTER, &tmp);
301: ISView(tmp, NULL);
302: ISDestroy(&tmp);
303: }
304: /* Add in local adjacency indices for owned dofs on interface (roots) */
305: for (p = pStart; p < pEnd; ++p) {
306: PetscInt numAdj = maxAdjSize, adof, dof, off, d, q;
308: PetscSectionGetDof(section, p, &dof);
309: PetscSectionGetOffset(section, p, &off);
310: if (!dof) continue;
311: PetscSectionGetDof(rootSectionAdj, off, &adof);
312: if (adof <= 0) continue;
313: DMDAGetAdjacency_Internal(dm, p, useClosure, tmpClosure, &numAdj, tmpAdj);
314: for (d = off; d < off+dof; ++d) {
315: PetscInt adof, aoff, i;
317: PetscSectionGetDof(rootSectionAdj, d, &adof);
318: PetscSectionGetOffset(rootSectionAdj, d, &aoff);
319: i = adof-1;
320: for (q = 0; q < numAdj; ++q) {
321: const PetscInt padj = tmpAdj[q];
322: PetscInt ndof, ncdof, ngoff, nd;
324: if ((padj < pStart) || (padj >= pEnd)) continue;
325: PetscSectionGetDof(section, padj, &ndof);
326: PetscSectionGetConstraintDof(section, padj, &ncdof);
327: PetscSectionGetOffset(sectionGlobal, padj, &ngoff);
328: for (nd = 0; nd < ndof-ncdof; ++nd) {
329: rootAdj[aoff+i] = ngoff < 0 ? -(ngoff+1)+nd : ngoff+nd;
330: --i;
331: }
332: }
333: }
334: }
335: /* Debugging */
336: if (debug) {
337: IS tmp;
338: PetscPrintf(comm, "Root adjacency indices\n");
339: ISCreateGeneral(comm, adjSize, rootAdj, PETSC_USE_POINTER, &tmp);
340: ISView(tmp, NULL);
341: ISDestroy(&tmp);
342: }
343: /* Compress indices */
344: PetscSectionSetUp(rootSectionAdj);
345: for (p = pStart; p < pEnd; ++p) {
346: PetscInt dof, cdof, off, d;
347: PetscInt adof, aoff;
349: PetscSectionGetDof(section, p, &dof);
350: PetscSectionGetConstraintDof(section, p, &cdof);
351: PetscSectionGetOffset(section, p, &off);
352: if (!dof) continue;
353: PetscSectionGetDof(rootSectionAdj, off, &adof);
354: if (adof <= 0) continue;
355: for (d = off; d < off+dof-cdof; ++d) {
356: PetscSectionGetDof(rootSectionAdj, d, &adof);
357: PetscSectionGetOffset(rootSectionAdj, d, &aoff);
358: PetscSortRemoveDupsInt(&adof, &rootAdj[aoff]);
359: PetscSectionSetDof(rootSectionAdj, d, adof);
360: }
361: }
362: /* Debugging */
363: if (debug) {
364: IS tmp;
365: PetscPrintf(comm, "Adjancency Section for Preallocation on Roots after compression:\n");
366: PetscSectionView(rootSectionAdj, PETSC_VIEWER_STDOUT_WORLD);
367: PetscPrintf(comm, "Root adjacency indices after compression\n");
368: ISCreateGeneral(comm, adjSize, rootAdj, PETSC_USE_POINTER, &tmp);
369: ISView(tmp, NULL);
370: ISDestroy(&tmp);
371: }
372: /* Build adjacency section: Maps global indices to sets of adjacent global indices */
373: PetscSectionGetOffsetRange(sectionGlobal, &globalOffStart, &globalOffEnd);
374: PetscSectionCreate(comm, §ionAdj);
375: PetscSectionSetChart(sectionAdj, globalOffStart, globalOffEnd);
376: for (p = pStart; p < pEnd; ++p) {
377: PetscInt numAdj = maxAdjSize, dof, cdof, off, goff, d, q;
378: PetscBool found = PETSC_TRUE;
380: PetscSectionGetDof(section, p, &dof);
381: PetscSectionGetConstraintDof(section, p, &cdof);
382: PetscSectionGetOffset(section, p, &off);
383: PetscSectionGetOffset(sectionGlobal, p, &goff);
384: for (d = 0; d < dof-cdof; ++d) {
385: PetscInt ldof, rdof;
387: PetscSectionGetDof(leafSectionAdj, off+d, &ldof);
388: PetscSectionGetDof(rootSectionAdj, off+d, &rdof);
389: if (ldof > 0) {
390: /* We do not own this point */
391: } else if (rdof > 0) {
392: PetscSectionSetDof(sectionAdj, goff+d, rdof);
393: } else {
394: found = PETSC_FALSE;
395: }
396: }
397: if (found) continue;
398: PetscSectionGetDof(section, p, &dof);
399: PetscSectionGetOffset(sectionGlobal, p, &goff);
400: DMDAGetAdjacency_Internal(dm, p, useClosure, tmpClosure, &numAdj, tmpAdj);
401: for (q = 0; q < numAdj; ++q) {
402: const PetscInt padj = tmpAdj[q];
403: PetscInt ndof, ncdof, noff;
405: if ((padj < pStart) || (padj >= pEnd)) continue;
406: PetscSectionGetDof(section, padj, &ndof);
407: PetscSectionGetConstraintDof(section, padj, &ncdof);
408: PetscSectionGetOffset(section, padj, &noff);
409: for (d = goff; d < goff+dof-cdof; ++d) {
410: PetscSectionAddDof(sectionAdj, d, ndof-ncdof);
411: }
412: }
413: }
414: PetscSectionSetUp(sectionAdj);
415: if (debug) {
416: PetscPrintf(comm, "Adjacency Section for Preallocation:\n");
417: PetscSectionView(sectionAdj, PETSC_VIEWER_STDOUT_WORLD);
418: }
419: /* Get adjacent indices */
420: PetscSectionGetStorageSize(sectionAdj, &numCols);
421: PetscMalloc1(numCols, &cols);
422: for (p = pStart; p < pEnd; ++p) {
423: PetscInt numAdj = maxAdjSize, dof, cdof, off, goff, d, q;
424: PetscBool found = PETSC_TRUE;
426: PetscSectionGetDof(section, p, &dof);
427: PetscSectionGetConstraintDof(section, p, &cdof);
428: PetscSectionGetOffset(section, p, &off);
429: PetscSectionGetOffset(sectionGlobal, p, &goff);
430: for (d = 0; d < dof-cdof; ++d) {
431: PetscInt ldof, rdof;
433: PetscSectionGetDof(leafSectionAdj, off+d, &ldof);
434: PetscSectionGetDof(rootSectionAdj, off+d, &rdof);
435: if (ldof > 0) {
436: /* We do not own this point */
437: } else if (rdof > 0) {
438: PetscInt aoff, roff;
440: PetscSectionGetOffset(sectionAdj, goff+d, &aoff);
441: PetscSectionGetOffset(rootSectionAdj, off+d, &roff);
442: PetscMemcpy(&cols[aoff], &rootAdj[roff], rdof * sizeof(PetscInt));
443: } else {
444: found = PETSC_FALSE;
445: }
446: }
447: if (found) continue;
448: DMDAGetAdjacency_Internal(dm, p, useClosure, tmpClosure, &numAdj, tmpAdj);
449: for (d = goff; d < goff+dof-cdof; ++d) {
450: PetscInt adof, aoff, i = 0;
452: PetscSectionGetDof(sectionAdj, d, &adof);
453: PetscSectionGetOffset(sectionAdj, d, &aoff);
454: for (q = 0; q < numAdj; ++q) {
455: const PetscInt padj = tmpAdj[q];
456: PetscInt ndof, ncdof, ngoff, nd;
457: const PetscInt *ncind;
459: /* Adjacent points may not be in the section chart */
460: if ((padj < pStart) || (padj >= pEnd)) continue;
461: PetscSectionGetDof(section, padj, &ndof);
462: PetscSectionGetConstraintDof(section, padj, &ncdof);
463: PetscSectionGetConstraintIndices(section, padj, &ncind);
464: PetscSectionGetOffset(sectionGlobal, padj, &ngoff);
465: for (nd = 0; nd < ndof-ncdof; ++nd, ++i) {
466: cols[aoff+i] = ngoff < 0 ? -(ngoff+1)+nd : ngoff+nd;
467: }
468: }
469: 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);
470: }
471: }
472: PetscSectionDestroy(&leafSectionAdj);
473: PetscSectionDestroy(&rootSectionAdj);
474: PetscFree(rootAdj);
475: PetscFree2(tmpClosure, tmpAdj);
476: /* Debugging */
477: if (debug) {
478: IS tmp;
479: PetscPrintf(comm, "Column indices\n");
480: ISCreateGeneral(comm, numCols, cols, PETSC_USE_POINTER, &tmp);
481: ISView(tmp, NULL);
482: ISDestroy(&tmp);
483: }
484: /* Create allocation vectors from adjacency graph */
485: MatGetLocalSize(A, &locRows, NULL);
486: PetscLayoutCreate(PetscObjectComm((PetscObject)A), &rLayout);
487: PetscLayoutSetLocalSize(rLayout, locRows);
488: PetscLayoutSetBlockSize(rLayout, 1);
489: PetscLayoutSetUp(rLayout);
490: PetscLayoutGetRange(rLayout, &rStart, &rEnd);
491: PetscLayoutDestroy(&rLayout);
492: /* Only loop over blocks of rows */
493: 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);
494: for (r = rStart/bs; r < rEnd/bs; ++r) {
495: const PetscInt row = r*bs;
496: PetscInt numCols, cStart, c;
498: PetscSectionGetDof(sectionAdj, row, &numCols);
499: PetscSectionGetOffset(sectionAdj, row, &cStart);
500: for (c = cStart; c < cStart+numCols; ++c) {
501: if ((cols[c] >= rStart*bs) && (cols[c] < rEnd*bs)) {
502: ++dnz[r-rStart];
503: if (cols[c] >= row) ++dnzu[r-rStart];
504: } else {
505: ++onz[r-rStart];
506: if (cols[c] >= row) ++onzu[r-rStart];
507: }
508: }
509: }
510: if (bs > 1) {
511: for (r = 0; r < locRows/bs; ++r) {
512: dnz[r] /= bs;
513: onz[r] /= bs;
514: dnzu[r] /= bs;
515: onzu[r] /= bs;
516: }
517: }
518: /* Set matrix pattern */
519: MatXAIJSetPreallocation(A, bs, dnz, onz, dnzu, onzu);
520: MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);
521: /* Check for symmetric storage */
522: MatGetType(A, &mtype);
523: PetscStrcmp(mtype, MATSBAIJ, &isSymBlock);
524: PetscStrcmp(mtype, MATSEQSBAIJ, &isSymSeqBlock);
525: PetscStrcmp(mtype, MATMPISBAIJ, &isSymMPIBlock);
526: if (isSymBlock || isSymSeqBlock || isSymMPIBlock) {MatSetOption(A, MAT_IGNORE_LOWER_TRIANGULAR, PETSC_TRUE);}
527: /* Fill matrix with zeros */
528: if (fillMatrix) {
529: PetscScalar *values;
530: PetscInt maxRowLen = 0;
532: for (r = rStart; r < rEnd; ++r) {
533: PetscInt len;
535: PetscSectionGetDof(sectionAdj, r, &len);
536: maxRowLen = PetscMax(maxRowLen, len);
537: }
538: PetscCalloc1(maxRowLen, &values);
539: for (r = rStart; r < rEnd; ++r) {
540: PetscInt numCols, cStart;
542: PetscSectionGetDof(sectionAdj, r, &numCols);
543: PetscSectionGetOffset(sectionAdj, r, &cStart);
544: MatSetValues(A, 1, &r, numCols, &cols[cStart], values, INSERT_VALUES);
545: }
546: PetscFree(values);
547: MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
548: MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
549: }
550: PetscSectionDestroy(§ionAdj);
551: PetscFree(cols);
552: return(0);
553: }