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